bootstrapCIs.R

getCI <- function(B, muH0, sdH0, N) {
    getM <- function(orgDV, idx) {
        bsM   <- mean(orgDV[idx])                       # M*
        bsS2M <- (((N-1) / N) * var(orgDV[idx])) / N    # S^2*(M)
        c(bsM, bsS2M)
    }

    DV  <- rnorm(N, muH0, sdH0)            # simulated data: original sample
    M   <- mean(DV)                        # M from original sample
    S2M <- (((N-1)/N) * var(DV)) / N       # S^2(M) from original sample

    # bootstrap
    boots   <- t(replicate(B, getM(DV, sample(seq(along=DV), replace=TRUE))))
    Mstar   <- boots[ , 1]                 # M* for each replicate
    S2Mstar <- boots[ , 2]                 # S^2*(M) for each replicate
    biasM   <- mean(Mstar) - M             # bias of estimator M

    # indices for sorted vector of estimates
    idx   <- trunc((B + 1) * c(0.05/2, 1 - 0.05/2))
    zCrit <- qnorm(c(1 - 0.05/2, 0.05/2))  # z-quantiles from std-normal distribution
    tStar <- (Mstar-M) / sqrt(S2Mstar)     # t*
    tCrit <- sort(tStar)[rev(idx)]         # t-quantiles from empirical t* distribution

    ciBasic <- 2*M - sort(Mstar)[rev(idx)] # basic CI
    ciPerc  <- sort(Mstar)[idx]            # percentile CI
    ciNorm  <- M-biasM - zCrit*sd(Mstar)   # normal CI
    ciT     <- M - tCrit * sqrt(S2M)       # studentized t-CI

    c(basic=ciBasic, percentile=ciPerc, normal=ciNorm, t=ciT)
}

## 1000 bootstraps - this will take a while
B    <- 999                  # number of replicates
muH0 <- 100                  # for generating data: true mean
sdH0 <- 40                   # for generating data: true sd
N    <- 200                  # sample size
DV   <- rnorm(N, muH0, sdH0) # simulated data: original sample
Nrep <- 1000                 # number of bootstraps
CIs  <- t(replicate(Nrep, getCI(B=B, muH0=muH0, sdH0=sdH0, N=N)))

## coverage probabilities
sum((CIs[ , "basic1"]      < muH0) & (CIs[ , "basic2"]      > muH0)) / Nrep
sum((CIs[ , "percentile1"] < muH0) & (CIs[ , "percentile2"] > muH0)) / Nrep
sum((CIs[ , "normal1"]     < muH0) & (CIs[ , "normal2"]     > muH0)) / Nrep
sum((CIs[ , "t1"]          < muH0) & (CIs[ , "t2"]          > muH0)) / Nrep

I just found this to be a very helpful function for my own research. So thank you! Just a quick question: what is the difference between B and Nrep? Is B just Nrep minus 1?