Calculate the mean and 95% confidence interval of a sample of numbers under a transform

sample.mean.CI.transform <- function(x, transform.func = identity, anti.transform.func = identity)
{
    ## Purpose: 
    ##   Calculate the mean and 95% CI of a sample of numbers under a transform. 
    ## 
    ## Author: Feiming Chen
    ## 
    ## Arguments: 
    ##   - x : a vector of positive numerical values
    ##   - transform.func: a function for transforming the original values to a new scale.
    ##   - anti.transform.func: a function for back-transforming the values to the original scale.
    ## 
    ## Return: 
    ##   - m : mean based on a particular transform function. 
    ## 
    ## ________________________________________________

    y <- transform.func(x)              # transform to a new scale
    a <- t.test(y)                      # estimate the mean and 95% CI based on t-distribution
    c(anti.transform.func(a$estimate), CI = anti.transform.func(a$conf.int))
}
if (F) {                                # Unit Test
    x <- exp(1:10)
    sample.mean.CI.transform(x)             # no transform
    ## estimate the mean under the log-normal distributional assumption
    sample.mean.CI.transform(x, transform.func = log, anti.transform.func = exp) 
}

Use simulation to estimate the coverage probability of the quantile rank-based confidence interval

The original code comes from https://stats.stackexchange.com/questions/99829/how-to-obtain-a-confidence-interval-for-a-percentile/284970#284970

test.coverage.via.simulation <- function(lu, n, p, alpha = 0.1)
{
    ## Purpose: 
    ##   Use simulation to estimate the coverage probability of the quantile rank-based CI.
    ##   
    ## Reference (code origin):
    ##   https://stats.stackexchange.com/questions/99829/how-to-obtain-a-confidence-interval-for-a-percentile/284970#284970
    ## 
    ## Arguments:
    ##   - lu: a vector of two numbers for the lower and upper limit of the CI in ranks. 
    ##   - n: sample size
    ##   - p: quantile probability cutpoint (between 0 and 1)
    ##   - alpha: type I error level (default to 0.1 so a 90% CI is calculated)
    ## 
    ## Return: 
    ##   Average coverage probability from 10,000 simulated samples.
    ## ________________________________________________

    ## Generate many random samples of size "n" from a known distribution
    ## and compute actual CI's from those samples using the given rank-based CI. 
    set.seed(1)
    n.sim <- 1e4
    index <- function(x, i) ifelse(i == Inf, Inf, ifelse(i == -Inf, -Inf, x[i]))
    sim <- replicate(n.sim, index(sort(rnorm(n)), lu))

    ## Compute the actual proportion of those intervals that cover the theoretical quantile.
    F.q <- qnorm(p)
    covers <- sim[1, ] <= F.q & F.q <= sim[2, ]
    mean.coverage <- mean(covers)
    message("Mean Coverage Over 10,000 Simulated Samples = ", signif(mean.coverage, 4))
}
if (F) {                                # Unit Test
    lu <- c(85, 97)
    n <- 100
    p <- 0.9
    alpha <- 0.05
    test.coverage.via.simulation(lu, n, p, alpha)
    ## Mean Coverage Over 10,000 Simulated Samples = 0.9528
}

Rank-based Confidence Interval for a Quantile

The original code comes from https://stats.stackexchange.com/questions/99829/how-to-obtain-a-confidence-interval-for-a-percentile/284970#284970

quantile.CI.ranks <- function(n, p, alpha = 0.1) {
    ## Purpose:
    ##   Calculate a two-sided near-symmetric distribution-free
    ##   confidence interval with confidence level of (1 - alpha) for
    ##   a quantile, by searching over a small range of upper and
    ##   lower order statistics for the closest coverage to (1 -
    ##   alpha) (but not less than it, if possible).
    ##
    ## Reference (code origin):
    ##   https://stats.stackexchange.com/questions/99829/how-to-obtain-a-confidence-interval-for-a-percentile/284970#284970
    ## 
    ## Arguments:
    ##   - n: sample size
    ##   - p: quantile probability cutpoint (between 0 and 1)
    ##   - alpha: type I error level (default to 0.1 so a 90% CI is calculated)
    ##
    ## Return:
    ##   - CI: two indices into the order statistics for the two-sided CI of the quantile
    ##   - coverage: theoretical coverage probability of the CI, which should be close to (1 - alpha).

    ## a small candidate list of order statistics for the lower/upper limits of the CI of a quantile
    l <- qbinom(alpha/2, n, p) + (-2:2) + 1 # lower limit candidates
    u <- qbinom(1 - alpha/2, n, p) + (-2:2)  # upper limit candidates

    ## out-of-bound order statistics correspond to no limit
    l[l < 0] <- -Inf                    # no lower limit
    u[u > n] <- Inf                     # no upper limit

    ## for each pair of candidate lower/upper limit, calculate the coverage probability
    coverage <- outer(l, u, function(l, u) pbinom(u - 1, n, p) - pbinom(l - 1, n, p))

    ## if no coverage is above (1 - alpha), choose the max coverage; 
    ## otherwise, look for the smallest coverage that is above (1 - alpha). 
    if (max(coverage) < 1 - alpha) {
        i <- which(coverage == max(coverage)) 
    } else {                            
        i <- which(coverage == min(coverage[coverage >= 1 - alpha]))
    }

    j <- i[1]                           # in case there are multiple candidates

    ## identify the order statistics and its coverage.
    L <- rep(l, 5)[j]
    U <- rep(u, each = 5)[j]
    return(list(CI = c(L, U), coverage = coverage[j]))
}
if (F) {                                # Unit Test
    alpha <- 0.05
    n <- 100
    p <- 0.9
    quantile.CI(n, p, alpha)
    ## $CI
    ## [1] 85 97

    ## $coverage
    ## [1] 0.95227
}

Bootstrap Confidence Interval for a Quantile

quantile.CI.via.bootstrap <- function(x, p, alpha = 0.1) {
    ## Purpose:
    ##   Calculate a two-sided confidence interval with confidence level of (1 - alpha) for
    ##   a quantile, based on the (computing intensive) bootstrap resampling method. 
    ##
    ## Arguments:
    ##   - x: a vector of values, representing a data sample. 
    ##   - p: probability cutpoint for the quantile (between 0 and 1).
    ##   - alpha: type I error level (default to 0.1 so a 90% CI is calculated)
    ##
    ## Return:
    ##   - CI: the lower and upper limits of the two-sided CI. 

    q <- quantile(x, probs = p)         
    message("Quantile Point Estimate = ", q, " (Probability Cutpoint = ", p, ")\n")

    ## Bootstrap resampling with 2000 replications
    library(boot)
    set.seed(1)
    b <- boot(x, function(x, i) quantile(x[i], probs = p), R = 2000)

    boot.ci(b, conf = 1 - alpha, type = c("norm", "basic", "perc", "bca"))

}
if (F) {                                # Unit Test
    x <- 1:100
    p <- 0.9
    alpha <- 0.05
    quantile.CI.via.bootstrap(x, p, alpha)
    ## Intervals : 
    ## Level      Normal              Basic         
    ## 95%   (84.50, 96.34 )   (85.10, 97.10 )  

    ## Level     Percentile            BCa          
    ## 95%   (83.1, 95.1 )   (83.3, 95.1 )  
}

Make a plot to compare point estimates and their confidence intervals

plot.compare.CI <- function(x, CI, ylab = "Estimates", ...)
{
    ## Purpose: Make a plot to compare point estimates and their confidence intervals
    ##          (presumably derived from different methods)
    ## Arguments:
    ##   x: a vector of point estimates
    ##   CI: a matrix of CI's. Each column is a CI pair.
    ##       1st row is lower CI. 2nd row is upper CI
    ##   ...: pass to the plot function. 
    ## Return: a plot. 
    ## Author: Feiming Chen
    ## ________________________________________________

    require(gplots)

    ## Determine how best to put X-axis Labels --------------------
    lx <- length(x)                     # number of point estimates

    if (lx > 7) {
        xlas <- 2
        o <- par(mar = c(10, 4, 4, 2))
        on.exit(par(o))
    }

    li <- CI[1,]
    ui <- CI[2,]
    plotCI(x, li = li, ui = ui, barcol="blue", xaxt="n", xlab = "", ylab = ylab, sfrac = 0.005, xlim = c(0.7, lx + 0.3), ...)

    s <- 1:lx
    ss <- names(x)
    axis(side=1, at = s, labels= ss, tick = F, las = 0)

    ## put on numerical labels on the CI's. 
    s1 <- s + 0.18                       # offset
    ## median & CI labels
    text(s1, y = li, labels = round(li, 2))
    text(s1, y = ui, labels = round(ui, 2))
    text(s1, y = x,  labels = round(x, 2))
}
if (F) {                                # Unit Test
    x <- 1:3
    names(x) <- c("A", "B", "C")
    CI <- matrix(c(0.8, 1.2, 1.7, 2.3, 2.9, 3.1), nrow = 2)
    plot.compare.CI(x, CI)
}

Preprocess clustered binary data before inference on sample proportion

preprocess.clustered.binary.data <- function(x, n)
{
    ## Purpose: Preprocess clustered binary data before inference on sample proportion
    ##          Reference: Rao, J. N. K., & Scott, A. J. (1992). A simple method for
    ##                     the analysis of clustered binary data. Biometrics, 577-585.
    ## Keywords: variance adjustment, ratio estimator, correlated data
    ## Arguments:
    ##   x: a vector for the numerators across clusters.
    ##   n: a vector for the denominators across clusters.
    ## Return: a pre-processed numerator and denominator for overall sample proportion
    ## Author: Feiming Chen
    ## ________________________________________________

    m <- length(x)
    cat("Number of Clusters =", m, "\n")

    n0 <- sum(n)
    cat("Raw Sample Size =", n0, "\n")

    x0 <- sum(x)
    cat("Raw Incidence Count =", x0, "\n")

    cat("\nEstimates Regarding the Overall Sample Proportion:\n")

    p <- x0 / n0
    cat("  Point Estimate =", round(p, 4), "\n")

    v0 <- p * (1 - p) / n0
    cat("  Naive Binomial Variance =", round(v0, 4), "\n")

    r <- x - n * p
    v <- m * sum(r^2) / (m - 1) / n0^2
    cat("  Correct Variance  =", round(v, 4), "\n")

    d <- v / v0
    cat("\nDesign Effect (Variance Inflation Factor due to Clustering) =", round(d, 4), "\n")

    n1 <- n0 / d
    cat("Effective Sample Size (n) =", round(n1), "\n")

    x1 <- x0 / d
    cat("Effective Incidence Count (x) =", round(x1), "\n")

    ## Return the pre-processed numerator (x) and denominator (n) for
    ## the overall sample proportion (p = x / n). 
    list(x = x1, n = n1)

}
if (F) {                                # Unit Test
    x = c(1,1,2,0,5,0,1,4,0,1,0,0,0,0,0,4,3,0,1,1,1,2,0,1,0,2,1,0,1,5,0,0,0,0,0,0,2,0,2,0,0,2,1,2,2,0,0,1,0,1) 
    n = c(2,2,3,0,7,1,2,5,1,2,0,0,4,2,0,7,4,1,1,1,4,2,3,1,0,2,1,0,1,5,3,0,1,3,0,0,2,1,2,0,0,5,3,2,2,1,0,2,1,1) 

    a <- preprocess.clustered.binary.data(x, n)
    ## Number of Clusters = 50 
    ## Raw Sample Size = 93 
    ## Raw Incidence Count = 50 

    ## Estimates Regarding the Overall Sample Proportion:
    ##   Point Estimate = 0.5376 
    ##   Naive Binomial Variance = 0.0027 
    ##   Correct Variance  = 0.004 

    ## Design Effect (Variance Inflation Factor due to Clustering) = 1.4895 
    ## Effective Sample Size (n) = 62 
    ## Effective Incidence Count (x) = 34 

    prop.test(a$x, a$n)
    ## 95 percent confidence interval:
    ##  0.40774 0.66290
    ## sample estimates:
    ##       p 
    ## 0.53763 

    ## Compare to Bootstrap Confidence Intervals
    library(boot)
    r <- boot(data.frame(x=x, n=n), function(dat, idx) { d <- dat[idx,]; sum(d$x)/sum(d$n)}, R = 10000)
    boot.ci(r)
    ## Bootstrap Percentile CI: ( 0.4118,  0.6569 ), which is a bit tighter than the VIF method. 

    ## Compare to Unadjusted (Wrong) CI: 
    prop.test(sum(x), sum(n))
    ## 95 percent confidence interval:
    ##  0.43159 0.64053, which is too narrow and is wrong. 

}

Bootstrap Confidence Interval for Data with Repeated Measures

stat.bootstrap.cluster <- function(id, val, func = mean, boot.size = 10000)
{
    ## Purpose: Calculate bootstrap-based 95% confidence interval for data with repeated measures
    ## Arguments:
    ##   id: uniquely identifies a subject (patient)
    ##   val: a numeric vector to be summarized.  It contains repeated measures per subject.
    ##   func: a function for calculating the summary statistic from a numeric vector.  Default to "mean". 
    ##   boot.size: number of bootstrap samples.  Default to 10000. 
    ## Return: Point estimate, bootstrap percentile 95% CI, histogram for the bootstrap distribution of the target statistic
    ## Author: Feiming Chen
    ## ________________________________________________

    fname <- deparse(substitute(func))
    ans <- func(val)
    cat("Summary Statistic:", fname, "=", ans, "\n")

    unique.ID <- unique(id)

    set.seed(1)
    replicate(n = boot.size, {
        s <- sample(unique.ID, replace = TRUE)  # a bootstrap sample of patient ID's
        ## find all rows with the ID's in the bootstrap sample of ID's 
        v <- c()
        for (j in s) v <- c(v, val[id == j]) # a new bootstrap sample
        func(v)             #  and its statistic
    }) -> est.boot

    hist(est.boot, xlab = fname, main = paste("Bootstrap Distribution: ", fname))
    cat("95% Bootstrap Percentile Confidence Interval:\n")
    quantile(est.boot, c(0.025, 0.975))
}
if (F) {                                # Unit Test
    id <- c(1, 1, 1, 2, 2, 3)
    val <- c(3, 3, 3, 4, 4, 5)
    stat.bootstrap.cluster(id, val)          # compare with t-test based CI: (2.8098 4.5235)
    ## Summary Statistic: mean = 3.6667 
    ## 95% Bootstrap Percentile Confidence Interval:
    ##  2.5% 97.5% 
    ##     3     5 
    stat.bootstrap.cluster(id, val, func = sd)
    ## Summary Statistic: sd = 0.8165 
    ## 95% Bootstrap Percentile Confidence Interval:
    ##   2.5%  97.5% 
    ## 0.0000 1.095
}