0
$\begingroup$

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!

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

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[2]
            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[2]
                    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)

}

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.