Using R to bootstrap the variance difference

Question

I asked this question on Stack Exchange, but I think it might be too specialized. Hopefully someone in the mixed model group can help me out.

I want to be able to bootstrap the variance differences between two data sets obtained at different times while taking out the error in a random effect.

I have 2 sets of experimental data, where the data was measured at 2 time points (initial and final). I also have a set of simulation data. I want to compare the variance of the simulated date with the variance difference between the experimental data (final - initial). The idea is to get confidence intervals from the bootstrap to compare the experimental data with the simulation.

I am having trouble making the statistic for the bootstrap function in the boot package for R. So far I have.

varcomp <- function ( formula, data, indices ) {
    d <- data[indices,] #sample for boot
    fit <- lmer(formula, data=d) #linear model
    res.var = (attr (VarCorr(fit), "sc")^2) # variance estimation
    return(res.var)
    }

But this function only returns the variance of a single data set. I want to be able to input 2 sets of data and have it return the difference between the two data sets' variance.

When I try something like:

varcomp <- function ( formula, data1, data2, indices ) {
d1 <- data1[indices,] #sample for boot
d2 <- data2[indices,] #sample for boot
fit1 <- lmer(formula, data=d1) #linear model
fit2 <- lmer(formula, data=d2) #linear model
a = (attr (VarCorr(fit1), "sc")^2) #output variance estimation
b = (attr (VarCorr(fit2), "sc")^2) #output variance estimation
drv = a - b #difference between the variance estimations
return(drv)
}

I would then put it into boot such as:

ip1.boot <- boot ( data = ip1, statistic=varcomp, R=100, formula=CNPC~(1|Cell.line:DNA.extract)+Cell.line)

I can't do it this way because the boot function only allows for one data set to be inputted.

Does anyone know how to create the correct statistic function for this?

An example of the data can also be downloaded here (2 csv files zipped 1.22KB.)

My data looks something like the following:

Initial

       Cell.line    Time DNA.extract   Gene      CNPC
1          9 initial           1 atubP1 1778.4589
2          9 initial           1 atubP1 2108.0552
3          9 initial           1 atubP1 2118.6725
4          9 initial           2 atubP1 2018.6593
5          9 initial           2 atubP1 1935.9008
6          9 initial           2 atubP1 1749.9158
7          9 initial           3 atubP1 1524.7475
8          9 initial           3 atubP1 1532.9781
9          9 initial           3 atubP1 1693.3098
10        17 initial           1 atubP1 1076.4720
11        17 initial           1 atubP1 1101.3315
12        17 initial           1 atubP1 1185.3606
13        17 initial           2 atubP1 1131.1118
14        17 initial           2 atubP1  892.7087
15        17 initial           2 atubP1 1028.5465
16        17 initial           3 atubP1  887.9972
17        17 initial           3 atubP1  732.9646
18        17 initial           3 atubP1  680.6724

Final

   Cell.line  Time DNA.extract   Gene      CNPC
1          9 final           1 atubP1 1262.2378
2          9 final           1 atubP1 1261.9858
3          9 final           1 atubP1 1390.6873
4          9 final           2 atubP1 1539.7180
5          9 final           2 atubP1 1510.5405
6          9 final           2 atubP1 1443.1767
7          9 final           3 atubP1 1456.2050
8          9 final           3 atubP1 1578.6396
9          9 final           3 atubP1 1656.1822
10        17 final           1 atubP1 1462.5179
11        17 final           1 atubP1 1580.9956
12        17 final           1 atubP1 1255.9020
13        17 final           2 atubP1  886.7579
14        17 final           2 atubP1  581.8116
15        17 final           2 atubP1  722.0526
16        17 final           3 atubP1 4168.7895
17        17 final           3 atubP1 3266.2105
18        17 final           3 atubP1 4219.5645

Answer 1

I believe that the key is that you can pass in a list, so varcomp would become:

varcomp <- function (data, indices, formula) {
d1 <- data$d1[indices,] #sample for boot
d2 <- data$d2[indices,] #sample for boot
fit1 <- lmer(formula, data=d1) #linear model
fit2 <- lmer(formula, data=d2) #linear model
a = (attr (VarCorr(fit1), "sc")^2) #output variance estimation
b = (attr (VarCorr(fit2), "sc")^2) #output variance estimation
drv = a - b #difference between the variance estimations
return(drv)
}

And you call boot with:

ip1.boot <- boot ( data = list (d1=dataframe1, d2=dataframe2), statistic=varcomp, R=100, formula=CNPC~(1|Cell.line:DNA.extract)+Cell.line)

That is, you'd stuff both sets of data (I assume data.frames) into a list and pull them out as needed. This does assume that both sets of data have the same number of observations, of course. At least this should put you on the right path, I hope...