R: interpretation of boot() output

Question

I'm trying to run an ANOVA using bootstrapped data (because my data are not normally distributed) but I don't really know if I did this correctly & how to make sense of my output.

This is what I'm trying to do: I conducted the same experiment with the same subjects online & in the lab (= independent variable "testing situation" with 2 factor levels). In the experiment, I manipulate cognitive load as an independent variable with 4 factor levels (called "no-back", "zero-back", "one-back" and "two-back") and I measure the reaction times (in ms) as a dependent variable.

This means I have a 2 x 4 within-subjects-design with the reaction times as the outcome variable and want to know if there are main or interaction effects.

What I tried to do is the following:

# write regression function
bootReg <- function(formula, # Formula of the regression
                    data, 
                    indices)
{
  d <- data[indices,]
  fit <- lm(formula, data=d)
  return(coef(fit))  
}

# bootstrap the data 
boot.object <- boot(statistic = bootReg, formula = lm(RT ~ Code + Situation + Block, data = dataframe), data = dataframe, R = 2000)

My output looks like this:

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = NBACK_DESCR, statistic = bootReg, R = 2000, formula = lm(NBACK_Median_RT ~ 
    Code + Situation + Block, data = NBACK_DESCR))


Bootstrap Statistics :
        original     bias    std. error
t1*   322.927313 -9.0002985    79.96588
t2*   -12.014833  5.6209447   117.02878
t3*   109.197500  0.8386920   120.86134
t4*   338.548500  1.0563602   123.06327
t5*   212.354750  0.5961423   307.84955
t6*   115.336083  1.0862478    78.74367
t7*   204.884583  0.6035880    94.50454
t8*  -119.986083  2.2980845    72.79074
t9*   -93.026833  3.3750698    79.26258
t10*    0.311750  7.5767305   183.46302
t11*  200.108625 -1.8049229   371.22341
t12*  -53.072917  0.2976762    95.20676
t13*  126.300083  3.3038699   107.50477
t14*   -3.794000  2.6890971    85.11730
t15*   68.130917  0.1380621   109.92370
t16* -144.711750  1.6015020    74.13766
t17*    0.920000  0.8054492    98.44356
t18* -120.711167  0.7836202    78.31914
t19*   10.794083 -0.6042305    98.66546
t20*  519.203600  9.8466741   571.22411
t21*   90.910500 -0.2344282    90.77725
t22*  108.026250  1.1320475    77.27769
t23*   16.168000  0.3672671   126.07834
t24*  284.315333 -2.4115301   287.93144
t25*  198.447917  2.9121272   112.64016
t26*   37.165250  1.5303775    94.42860
t27*  -98.688833  3.0493664    79.98359
t28*   45.922417  2.0774330    74.65226
t29*   -6.227517  3.8654708   166.54048
t30*   50.998118  2.9716901    49.62328
t31*  -23.885188  6.9669819    64.99859
t32*   59.188070 10.5457197    73.22344

Does anyone know what this means and how I can see if I have significant interaction or main effects?

I guess t1* is the original test statistic & the other t*s are the bootstrapped test statistics, but even if that's right, that doesn't really help me with understanding what this output is trying to tell me.

My idea was to count how many ts are > t1, divide that by the number of samples (in this case 2000? Or 31?) to get p-values. Also I thought about doing that for different kinds of models with different combinations of predictors to see which are significant. Does that make sense?! I really don't know. Also I guess I should apply a correction?

It would be really great if anyone could help me with this - I'm an undergrad currently trying to learn R programming and I'm completely lost! Thanks in advance!

user_12258459 · Answer

t1-32 are your coefficients. Since you did not provide the data.frame, I use an example below:

library(boot)
set.seed(111)

dataframe = data.frame(RT=rnorm(100),Code=rbinom(100,1,0.5),
Situation=factor(sample(1:3,100,replace=TRUE)),
Block=sample(letters[1:3],100,replace=TRUE))

dataframe$RT[dataframe$Code==1] = dataframe$RT[dataframe$Code==1] + 1

Before we do the bootstrap, if we run the linear model, we expect the output to be like this:

(Intercept)        Code  Situation2  Situation3      Blockb      Blockc 
-0.02469464  0.78758240 -0.10768677  0.32013080  0.29325885 -0.05753515

You have a vector of 6, 1 intercept 5 coefficients.. since some of them are factors with > 2 levels. Now we bootstrap:

# bootstrap the data 
boot.object <- boot(statistic = bootReg, formula = lm(RT ~ Code + Situation + Block, data = dataframe), data = dataframe, R = 2000)

    Call:
boot(data = dataframe, statistic = bootReg, R = 2000, formula = lm(RT ~ 
    Code + Situation + Block, data = dataframe))


Bootstrap Statistics :
       original       bias    std. error
t1* -0.02469464  0.004240916   0.2588214
t2*  0.78758240  0.005184110   0.2130492
t3* -0.10768677  0.002780318   0.2479104
t4*  0.32013080  0.001647815   0.2734137
t5*  0.29325885 -0.005158781   0.2510207
t6* -0.05753515 -0.012981051   0.2718108

You can see that the printed output corresponds to the exact coefficients when you run lm on the original data. Another way to look at it:

(Intercept)        Code  Situation2  Situation3      Blockb      Blockc 
-0.02469464  0.78758240 -0.10768677  0.32013080  0.29325885 -0.05753515

The bootstrapped values are stored under $t, here you can see there are 6 columns, one for every coefficient and each row is a bootstrap:

head(boot.object$t)
            [,1]      [,2]         [,3]       [,4]        [,5]        [,6]
[1,]  0.37081996 0.3009173  0.307350121  0.2271736 -0.02898838 -0.38958835
[2,] -0.09836689 1.0306144 -0.272608134  0.1617208  0.30521958 -0.09391564
[3,] -0.24588583 0.9835756 -0.416804093  0.1581820  0.28454367  0.24282730
[4,]  0.29403111 0.5777657 -0.283601680 -0.1328344  0.20086620 -0.20614676
[5,] -0.00692040 0.6228231 -0.150136418  0.3648773  0.42969597 -0.07899494
[6,] -0.24859844 0.8226603  0.008036868  0.6543648  0.43781238  0.25347543

Your bootstrapped values should hover around your observed values, here I plot the coefficient "Code" :

hist(boot.object$t[,2],br=50)
abline(v=boot.object$t0[2],col="blue")

From the bootstrap, you can estimate the standard error of your coefficient terms and use it to construct a confidence interval. It's not used to test the hypothesis that your coefficient is not zero.

You mixed up permutation test with bootstrap. What you need is to construct a similar test where you swap the labels. You can check out posts such as this or maybe this by Ben Bolker

R: interpretation of boot() output

1 Answers