Bootstrapped Hypothesis Test Power with Covariate Adjustment

I am simulating data in which I ultimately test a difference in means (i.e difference between drug and placebo). I test this using a standard T-Test, and then also with a bootstrapped T-Test. I am doing this using 3 differently distributed outcomes - Normal outcome, Exponential outcome, and Gamma Outcome. And I am also interested in adjusting for a covariate.

I fit a linear regression to the data and then select the estimates for treatment effect and it's standard error to construct the test statistic. I do this same thing when I adjust for a covariate.

I am finding that the power from the bootstrapped procedure is the same as if I just computed the T-Test (without bootstrapping) when I do not adjust for a covariate and for all outcome types, but when I adjust for a covariate, the bootstrapping procedure has much less power than the T-Test (without bootstrapping). Is this to be expected? Number of simulations and sample size remain the same, the only thing that changes is the addition of a covariate.

Some background: According to the bootstrap procedure for hypothesis testing a difference in means (Efron and Tibshirani), we must "re-center" the dataset from which we bootstrap. So we compute the test statistic, and then re-center that data to create a more appropriate null distribution when we select the bootstrapped samples. I don't know if this is relevant to answering my questions, but I thought that I should give some background

Thanks for any suggestions you may have!

• Can you say more about how you did this? Can you provide pseudocode or real code? – gung Aug 17 '17 at 15:41

This is a modified version of my actual code - I changed it since I am calling some other user defined functions, but it should give you the jist of what I am doing.

n <- c(1200)
dd = 1000
ff = 8000
set.seed(62)
pval <- rep(NA,length(n))

for (i in 1:length(n)) {
pr.more.extreme <- foreach(j=1:ff,.packages = c("MASS","RcppArmadillo","dplyr")) %dopar% {

## 1) Create Data with columns trt1 (values P and D), visit (outcome of interest), cov (covariate)
## The first 2/3 of the dataset are P and the last 1/3 is D
## 2) Analysis
X <- cbind(1,as.integer(Data$trt == "P"),Data$cov)
Y <- Data$visit lmStage1 <- fastLm(X,Y) coef <- lmStage1$coefficients
var <- lmStage1$stderr[2,1]^2 DataTS <- coef/sqrt(var) # Standardizing Dataset Pmean = mean(Data[which(Data[,1]=="P"),2]) Amean = mean(Data[which(Data[,1]=="A"),2]) S1mean = mean(Data[,2]) Data_stand <- Data Data_stand[(1:(2*n[i]/3)),4] = Data[(1:(2*n[i]/3)),2] - Pmean + S1mean Data_stand[((2*n[i]/3+1):n[i]),4] = Data[((2*n[i]/3+1):n[i]),2] - Amean + S1mean Data <- Data_stand # Standardizing Dataset END local.indicators <- as.data.frame(matrix(NA, nrow(Data), 1)) colnames(local.indicators) <- c("trt1.P") local.indicators[,"trt1.P"] <- Data$trt1 == "P"
X1 <- cbind(1, local.indicators$trt1.P) coef.boot <- TS.boot <- var.boot <- rep(NA,dd) ## 3) Bootstrap for (mm in 1:dd) { #if (mm %% 1000 ==0) {print(mm/1000)} tmp <- nrow(Data)/3 samp <- c(sample(1:(2*tmp), (2*tmp), replace=T), sample((2*tmp+1):(3*tmp), tmp, replace=T)) # Analysis Y1 <- Data[samp,2] lm.boot <- fastLm(cbind(X1,Data$cov[samp]), Y1)

## Coefficients, Variances, and Test Statistics
coef.boot[mm] <- lm.boot$coefficients var.boot[mm] <- lm.boot$stderr[2,1]^2

TS.boot[mm] <- coef.boot[mm]/sqrt(var.boot[mm])

}

(sum (abs(TS.boot) >= abs(DataTS)) + 1)/ (dd+1))

}

pval[i] <- sum(unlist(pr.more.extreme)<0.05)/ff
print(pval)

}

• This should be an addition to your question and not given as an answer. – Michael Chernick Aug 17 '17 at 17:28