Uso de R para el bootstrap de la diferencia de varianza

Question

Hice esta pregunta en Stack Exchange, pero creo que puede ser demasiado especializada. Espero que alguien del grupo de modelos mixtos pueda ayudarme.

Quiero ser capaz de bootstrap las diferencias de varianza entre dos conjuntos de datos obtenidos en diferentes momentos, mientras que sacar el error en un efecto aleatorio.

Tengo 2 conjuntos de datos experimentales, en los que los datos se midieron en 2 puntos temporales (inicial y final). También tengo un conjunto de datos de simulación. Quiero comparar la varianza del dato simulado con la diferencia de varianza entre los datos experimentales (final - inicial). La idea es obtener intervalos de confianza a partir del bootstrap para comparar los datos experimentales con la simulación.

Estoy teniendo problemas para hacer la estadística para la función bootstrap en el boot para R. Hasta ahora lo he hecho.

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)
    }

Pero esta función sólo devuelve la varianza de un único conjunto de datos. Quiero ser capaz de introducir 2 conjuntos de datos y que devuelva la diferencia entre la varianza de los dos conjuntos de datos.

Cuando intento algo como:

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)
}

Luego lo pondría en arranque como:

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

No puedo hacerlo así porque la función de arranque sólo permite introducir un conjunto de datos.

¿Alguien sabe cómo crear la función estadística correcta para esto?

También puede descargarse aquí un ejemplo de los datos (2 archivos csv comprimidos en 1,22KB.)

Mis datos son parecidos a los siguientes:

Inicial

       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

Creo que la clave está en que se puede pasar una lista, por lo que varcomp se convertiría:

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)
}

Y llamas boot con:

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

Es decir, meterías ambos conjuntos de datos (supongo que data.frames) en una lista y los sacarías cuando fuera necesario. Esto supone que ambos conjuntos de datos tienen el mismo número de observaciones, por supuesto. Al menos esto debería ponerte en el camino correcto, espero...