Prueba t emparejada como un caso especial de modelado lineal de efectos mixtos

20

Sabemos que una prueba t emparejada es solo un caso especial de ANOVA de medidas repetidas unidireccionales (o dentro del sujeto), así como el modelo lineal de efectos mixtos, que se puede demostrar con la función lme () del paquete nlme en R Como se muestra abajo.

#response data from 10 subjects under two conditions
x1<-rnorm(10)
x2<-1+rnorm(10)

# Now create a dataframe for lme
myDat <- data.frame(c(x1,x2), c(rep("x1", 10), rep("x2", 10)), rep(paste("S", seq(1,10), sep=""), 2))
names(myDat) <- c("y", "x", "subj")

Cuando ejecuto la siguiente prueba t emparejada:

t.test(x1, x2, paired = TRUE)

Obtuve este resultado (obtendrá un resultado diferente debido al generador aleatorio):

t = -2.3056, df = 9, p-value = 0.04657

Con el enfoque ANOVA podemos obtener el mismo resultado:

summary(aov(y ~ x + Error(subj/x), myDat))

# the F-value below is just the square of the t-value from paired t-test:
          Df  F value Pr(>F)
x          1  5.3158  0.04657

Ahora puedo obtener el mismo resultado en lme con el siguiente modelo, suponiendo una matriz de correlación simétrica definida positiva para las dos condiciones:

summary(fm1 <- lme(y ~ x, random=list(subj=pdSymm(form=~x-1)), data=myDat))

# the 2nd row in the following agrees with the paired t-test
# (Intercept) -0.2488202 0.3142115  9 -0.7918878  0.4488
# xx2          1.3325786 0.5779727  9  2.3056084  0.0466

U otro modelo, suponiendo una simetría compuesta para la matriz de correlación de las dos condiciones:

summary(fm2 <- lme(y ~ x, random=list(subj=pdCompSymm(form=~x-1)), data=myDat))

# the 2nd row in the following agrees with the paired t-test
# (Intercept) -0.2488202 0.4023431  9 -0.618428  0.5516
# xx2          1.3325786 0.5779727  9  2.305608  0.0466

Con la prueba t pareada y el ANOVA de medidas repetidas unidireccionales, puedo escribir el modelo de media celular tradicional como

Yij = μ + αi + βj + εij, i = 1, 2; j = 1, ..., 10

donde i indexa la condición, j indexa el sujeto, Y ij es la variable de respuesta, μ es constante para el efecto fijo para la media general, α i es el efecto fijo para la condición, β j es el efecto aleatorio para el sujeto después de N (0, σ p 2 ) (σ p 2 es la varianza de la población), y ε ij es residual después de N (0, σ 2 ) (σ 2 es la varianza dentro del sujeto).

Pensé que el modelo de media de celda anterior no sería apropiado para los modelos lme, pero el problema es que no puedo encontrar un modelo razonable para los dos enfoques lme () con el supuesto de estructura de correlación. La razón es que el modelo lme parece tener más parámetros para los componentes aleatorios que los que ofrece el modelo de media celular anterior. Al menos el modelo lme proporciona exactamente el mismo valor F, grados de libertad y valor p también, que gls no puede. Más específicamente, gls da DF incorrectos debido al hecho de que no tiene en cuenta el hecho de que cada sujeto tiene dos observaciones, lo que lleva a DF muy inflados. Lo más probable es que el modelo lme esté sobre-parametrizado al especificar los efectos aleatorios, pero no sé cuál es el modelo y cuáles son los parámetros. Entonces el problema aún no está resuelto para mí.

Bluepole
fuente
2
No estoy muy seguro de lo que estás preguntando. El modelo que anotó es precisamente el modelo para el modelo de efectos aleatorios; La estructura de correlación es inducida por el efecto aleatorio.
Aaron - Restablece a Monica
@ Aaron: el efecto aleatorio βj en el modelo de media celular se supone que sigue a N (0, σp2). Mi confusión es, ¿cómo se asocia este término (con un solo parámetro σp2) con la estructura de correlación especificada por simetría compuesta o una matriz simétrica simple en el modelo lme?
bluepole
Cuando calcula la correlación entre las dos observaciones sobre el mismo tema, la correlación es sigma_p ^ 2 / (sigma_p ^ 2 + sigma ^ 2) porque comparten la misma beta_j. Ver Pinheiro / Bates p.8. Además, el modelo de efectos aleatorios tal como lo escribió es equivalente a la simetría compuesta; Otras estructuras de correlación son más complejas.
Aaron - Restablece a Mónica
@ Aaron: ¡Gracias! Ya leí el libro de Pinheiro / Bates sobre esto, y todavía no podía entender los detalles sobre los efectos aleatorios. Las páginas más relevantes parecen ser el ejemplo en P.160-161. Además, la salida de efectos aleatorios de lme () con suposición de simetría compuesta no parece estar de acuerdo con la correlación de σp2 / (σp2 + σ2) en el modelo de media celular. Todavía desconcertado sobre la estructura del modelo.
bluepole
Bueno, casi equivalente a la simetría compuesta; en CS la correlación puede ser negativa pero no con efectos aleatorios. Quizás ahí es donde surge tu diferencia. Consulte stats.stackexchange.com/a/14185/3601 para más detalles.
Aaron - Restablece a Monica

Respuestas:

16

La equivalencia de los modelos se puede observar calculando la correlación entre dos observaciones del mismo individuo, de la siguiente manera:

Yyoj=μ+αyo+βj+ϵyojβjnorte(0 0,σpag2)ϵyojnorte(0 0,σ2)Cov(yyok,yjk)=Cov(μ+αyo+βk+ϵyok,μ+αj+βk+ϵjk)=Cov(βk,βk)=σpag2Vunr(yyok)=Vunr(yjk)=σpag2+σ2σpag2/ /(σpag2+σ2)

Sin embargo, tenga en cuenta que los modelos no son del todo equivalentes, ya que el modelo de efectos aleatorios obliga a la correlación a ser positiva. El modelo CS y el modelo t-test / anova no.

EDITAR: También hay otras dos diferencias. Primero, los modelos CS y de efectos aleatorios suponen normalidad para el efecto aleatorio, pero el modelo t-test / anova no. En segundo lugar, los modelos CS y de efectos aleatorios se ajustan utilizando la máxima probabilidad, mientras que el anova se ajusta utilizando cuadrados medios; cuando todo esté equilibrado estarán de acuerdo, pero no necesariamente en situaciones más complejas. Finalmente, desconfiaría de usar los valores de F / df / p de los diversos ajustes como medidas de cuánto coinciden los modelos; vea la famosa regla de Doug Bates en df para más detalles. (EDICIÓN FINAL)

El problema con su Rcódigo es que no está especificando la estructura de correlación correctamente. Debe usar glscon la corCompSymmestructura de correlación.

Genere datos para que haya un efecto sujeto:

set.seed(5)
x <- rnorm(10)
x1<-x+rnorm(10)
x2<-x+1 + rnorm(10)
myDat <- data.frame(c(x1,x2), c(rep("x1", 10), rep("x2", 10)), 
                    rep(paste("S", seq(1,10), sep=""), 2))
names(myDat) <- c("y", "x", "subj")

Entonces, así es como encajaría los efectos aleatorios y los modelos de simetría compuesta.

library(nlme)
fm1 <- lme(y ~ x, random=~1|subj, data=myDat)
fm2 <- gls(y ~ x, correlation=corCompSymm(form=~1|subj), data=myDat)

Los errores estándar del modelo de efectos aleatorios son:

m1.varp <- 0.5453527^2
m1.vare <- 1.084408^2

Y la correlación y la varianza residual del modelo CS es:

m2.rho <- 0.2018595
m2.var <- 1.213816^2

Y son iguales a lo que se espera:

> m1.varp/(m1.varp+m1.vare)
[1] 0.2018594
> sqrt(m1.varp + m1.vare)
[1] 1.213816

Otras estructuras de correlación generalmente no se ajustan a los efectos aleatorios, sino simplemente especificando la estructura deseada; Una excepción común es el modelo de efectos aleatorios AR (1) +, que tiene un efecto aleatorio y la correlación AR (1) entre observaciones sobre el mismo efecto aleatorio.

EDIT2: cuando me ajusto a las tres opciones, obtengo exactamente los mismos resultados, excepto que gls no intenta adivinar el df para el término de interés.

> summary(fm1)
...
Fixed effects: y ~ x 
                 Value Std.Error DF   t-value p-value
(Intercept) -0.5611156 0.3838423  9 -1.461839  0.1778
xx2          2.0772757 0.4849618  9  4.283380  0.0020

> summary(fm2)
...
                 Value Std.Error   t-value p-value
(Intercept) -0.5611156 0.3838423 -1.461839  0.1610
xx2          2.0772757 0.4849618  4.283380  0.0004

> m1 <- lm(y~ x + subj, data=myDat)
> summary(m1)
...
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -0.3154     0.8042  -0.392  0.70403   
xx2           2.0773     0.4850   4.283  0.00204 **

(La intersección es diferente aquí porque con la codificación predeterminada, no es la media de todos los sujetos, sino la media del primer sujeto).

También es interesante notar que el lme4paquete más nuevo da los mismos resultados pero ni siquiera trata de calcular un valor p.

> mm1 <- lmer(y ~ x + (1|subj), data=myDat)
> summary(mm1)
...
            Estimate Std. Error t value
(Intercept)  -0.5611     0.3838  -1.462
xx2           2.0773     0.4850   4.283
Aaron - Restablece a Monica
fuente
¡Gracias de nuevo por la ayuda! Conozco esta parte desde la perspectiva del modelo de media celular. Sin embargo, con el siguiente resultado de lme () con simetría compuesta: Efectos aleatorios: Fórmula: ~ x - 1 | Subj Estructura: Compuesto Simetría StdDev xx1 1.1913363 xx2 1.1913363 Corr: -0.036 Residual 0.4466733. Todavía no puedo conciliar estos números con el modelo de media celular. ¿Quizás pueda ayudarme más a resolver estos números?
bluepole
Además, ¿alguna idea sobre la formulación del modelo con otras estructuras de correlación como la matriz simétrica simple?
bluepole
¡Veo! Debería haber leído su respuesta en el otro hilo con más cuidado. Pensé en usar gls () antes, pero no pude entender la especificación de correlación. Es interesante que lme () con estructura de simetría compuesta para el efecto aleatorio siga representando el mismo valor t, pero parece que las variaciones para los efectos aleatorios no son directamente interpretables. ¡Realmente aprecio tu ayuda!
bluepole
Después de pensarlo dos veces, siento que mi confusión original aún no se ha resuelto. Sí, los gls se pueden usar para demostrar la estructura de correlación y los rones cuadrados medios, pero el modelo que se encuentra debajo no es exactamente el mismo que la prueba t pareada (o ANOVA de medidas repetidas unidireccionales en general), y dicha evaluación es más respaldado por los DF incorrectos y el valor p de gls. Por el contrario, mi comando lme con simetría compuesta proporciona los mismos F, DF y valor p. Lo único que me desconcierta es cómo se parametriza el modelo lme como se indica en mi publicación original. ¿Alguna ayuda por ahí?
bluepole
No estoy seguro de cómo ayudarte. ¿Podrías escribir cuáles crees que son los dos modelos diferentes? Algo está mal en cómo estás pensando en uno de ellos.
Aaron - Restablece a Mónica
3

También puede considerar el uso de la función mixeden el paquete afexpara devolver los valores de p con la aproximación df de Kenward-Roger, que devuelve valores de p idénticos como una prueba t emparejada:

library(afex)
mixed(y ~ x + (1|subj), type=3,method="KR",data=myDat) 

O

library(lmerTest)
options(contrasts=c('contr.sum', 'contr.poly'))
anova(lmer(y ~ x + (1|subj),data=myDat),ddf="Kenward-Roger")
Tom Wenseleers
fuente