## ----include=FALSE------------------------------------------------------------ ##### _Advanced Data Analysis from an Elementary Point of View_ ##### # R code for the chapter "Bootstrap" # Please do not re-distribute or use without attribution # http://www.stat.cmu.edu/~cshalizi/ADAfaEPoV/ # All demos omitted from this version, for re-use in later assignments, etc. # For the full version, see the URL above ## ----rboot-bootstrap-bootstrap.se--------------------------------------------- # Generate random values of a statistic by repeatedly running a simulator # Inputs: function to calculate the statistic (statistic) # function to run the simulation (simulator) # number of replicates (B) # Output: array of bootstrapped values of the statistic, with B columns # To work more nicely with other functions, a vector is converted to an # array of dimensions 1*B rboot <- function(statistic, simulator, B) { tboots <- replicate(B, statistic(simulator())) if(is.null(dim(tboots))) { tboots <- array(tboots, dim=c(1, B)) } return(tboots) } # Summarize the sampling distribution of a statistic, obtained by repeatedly # running a simultor # Inputs: array of bootstrapped statistics values (tboots) # function that summarizes the distribution (summarizer) # optional additional arguments to summarizer (...) # Output: vector giving a summary of the statistic # Presumes: tboots is an array with one column per simulation # each row of tboots is a separate component of the statistic # applying the summarizer to each row separately makes sense bootstrap <- function(tboots, summarizer, ...) { summaries <- apply(tboots, 1, summarizer, ...) # using apply() like this has interchanged rows and columns # because each chunk processed by apply() results in a new column, but # here those chunks are the rows of tboots # therefore use transpose to restore original orientation return(t(summaries)) } # Calculate a bootstrap standard error from scratch # Inputs: function to calculate the statistic (statistic) # function to run the simulation (simulator) # number of replicates (B) # Output: standard error for each bootstrap.se <- function(statistic, simulator, B) { bootstrap(rboot(statistic, simulator, B), summarizer=sd) } ## ----bootstrap.bias----------------------------------------------------------- # Calculate bootstrap biases # Inputs: function to run the simulation (simulator) # function to calculate the statistic (statistic) # number of replicates (B) # observed value of the statistic (t.hat) # Outputs: difference between mean of replicates and observed value bootstrap.bias <- function(simulator, statistic, B, t.hat) { # What's the expected value of the statistic, according to the bootstrap? expect <- bootstrap(rboot(statistic, simulator, B), summarizer=mean) # Bias is expected value minus truth return(expect-t.hat) } ## ----bootstrap.ci------------------------------------------------------------- # Find equal-tail interval with specified probability # Inputs: vector of values to sort (x) # total tail probability (alpha) # Output: length-two vector, giving interval of probability 1-alpha, with # probability alpha/2 in each tail equitails <- function(x, alpha) { lower <- quantile(x, alpha/2) upper <- quantile(x, 1-alpha/2) return(c(lower, upper)) } # Calculate (basic or pivotal) bootstrap confidence interval # Inputs: function to calculate the statistic (statistic) # function to run the simulation (simulator) # optional array of bootstrapped values (tboots) # if this is not NULL, over-rides the statistic & simulator arguments # number of replicates (B) # observed value of statistic (t.hat) # confidence level (level) # Outputs: two-column array with lower and upper confidence limits bootstrap.ci <- function(statistic=NULL, simulator=NULL, tboots=NULL, B=if(!is.null(tboots)) { ncol(tboots) }, t.hat, level) { # draw the bootstrap values, if not already provided if (is.null(tboots)) { # panic if we're not given an array of simulated values _and_ also lack # the means to calculate it for ourselves stopifnot(!is.null(statistic)) stopifnot(!is.null(simulator)) stopifnot(!is.null(B)) tboots <- rboot(statistic, simulator, B) } # easier to work with error probability than confidence level alpha <- 1-level # Calculate probability intervals for each coordinate intervals <- bootstrap(tboots, summarizer=equitails, alpha=alpha) # Re-center the intervals around the observed values upper <- t.hat + (t.hat - intervals[,1]) lower <- t.hat + (t.hat - intervals[,2]) # calculate CIs, centered on observed value plus bootstrap fluctuations # around it CIs <- cbind(lower=lower, upper=upper) return(CIs) } ## ----bootstrap-p-value-------------------------------------------------------- # Calculate a bootstrap p-value # Inputs: function to calculate a test statistic (test) # function to run the simulation (simulator) # number of replicates (B) # observed value of the test statistic (testhat) # Outputs: p-value for the hypothesis test # Presumes: larger values of the test statistic are stronger evidence against # the null hypothesis boot.pvalue <- function(test,simulator,B,testhat) { # bootstrap B values of the test statistic testboot <- rboot(B=B, statistic=test, simulator=simulator) # What proportion of simulated test statistics are at least as extreme as # the observed? # The +1 in numerator and denominator avoids the embarrassment of claiming # a p-value is exactly 0 on the basis of finite simulations p <- (sum(testboot >= testhat)+1)/(B+1) return(p) } ## ----double-bootstrap--------------------------------------------------------- # Calculate a p-value by two levels of bootstrapping # Useful when an estimated parameter affects the distribution of the test # Inputs: function to calculate a test statistic (test) # function to run the simulation (simulator) # number of replicates for top-level bootstrap (B1) # number of replicates per top-level replicate (B2) # function to estimate parameters (estimator) # estimate of parameter on actual data (thetahat) # observed value of the test statistic (testhat) # optional additional arguments to simulator (...) # Outputs: p-value for the hypothesis test # Presumes: larger values of the test statistic are stronger evidence against # the null hypothesis # simulator() can take thetahat as an argument # estimator() returns a value which simulator() can take as an argument doubleboot.pvalue <- function(test, simulator, B1, B2, estimator, thetahat, testhat, ...) { # For each top-level or outer replicate for (i in 1:B1) { # Run the simulator at the estimated parameter value xboot <- simulator(theta=thetahat, ...) # Re-estimate and re-calculate the test statistic thetaboot <- estimator(xboot) testboot[i] <- test(xboot) # Calculate a bootstrapped p-value for _that_ test pboot[i] <- boot.pvalue(test, simulator, B2, testhat=testboot[i], theta=thetaboot) } # EXERCISE for the reader: replace that for() loop with something vectorized # Get an unadjusted p-value for our observed test statistic p <- (sum(testboot >= testhat)+1)/(B1+1) # How extreme is our un-adjusted p-value? p.adj <- (sum(pboot <= p)+1)/(B1+1) return(p.adj) } ## ----resample----------------------------------------------------------------- # Resample a vector # That is, treat a sample as though it were a whole population, and draw # from it by sampling-with-replacement until we have a simulated data set # as big as the original # Equivalently, do IID draws from the empirical distribution # Inputs: vector to resample (x) # Outputs: vector of resampled values resample <- function(x) { sample(x,size=length(x),replace=TRUE) } # Resample whole rows from a data frame # That is, treat the rows as a population, and sample them with replacement # until we have a new data frame the same size as the original # Equivalently, draw IIDly from the joint empirical distribution over all # variables/columns # Inputs: data frame to resample (data) # Outputs: new data frame resample.data.frame <- function(data) { # Resample the row indices sample.rows <- resample(1:nrow(data)) # Return a new data frame with those rows in that order return(data[sample.rows,]) }