La prueba t pareada como caso especial de la modelización lineal de efectos mixtos

Question

Sabemos que una pareja t -es sólo un caso especial de ANOVA unidireccional de medidas repetidas (o dentro de un sujeto), así como de modelo lineal de efectos mixtos, que puede demostrarse con la función lme() del paquete nlme en R, como se muestra a continuación.

#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 pareada:

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

Obtuve este resultado (usted 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, asumiendo una matriz de correlación simétrica positiva-definida 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, asumiendo 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 emparejada y el ANOVA unidireccional de medidas repetidas, puedo escribir el modelo tradicional de media celular 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 de la media global, _i es el efecto fijo de la condición, _j es el efecto aleatorio para el sujeto siguiente 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 intra-sujeto).

Pensaba que el modelo de media de celda anterior no sería apropiado para los modelos lme, pero el problema es que no consigo un modelo razonable para las dos aproximaciones 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, lo que gls no puede. Más concretamente, gls proporciona DF incorrectos debido a 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é sobreparametrizado al especificar los efectos aleatorios, pero no sé cuál es el modelo ni cuáles son los parámetros. Así que la cuestión sigue sin resolverse para mí.

Answer 1

La equivalencia de los modelos puede observarse calculando la correlación entre dos observaciones de un mismo individuo, como se indica a continuación:

Como en su notación, dejemos que $Y_{ij} = \mu + \alpha_i + \beta_j + \epsilon_{ij}$ donde $\beta_j \sim N(0, \sigma_p^2)$ y $\epsilon_{ij} \sim N(0, \sigma^2)$ . Entonces $Cov(y_{ik}, y_{jk}) = Cov(\mu + \alpha_i + \beta_k + \epsilon_{ik}, \mu + \alpha_j + \beta_k + \epsilon_{jk}) = Cov(\beta_k, \beta_k) = \sigma_p^2$ porque todos los demás términos son independientes o fijos, y $Var(y_{ik}) = Var(y_{jk}) = \sigma_p^2 + \sigma^2$ , por lo que la correlación es $\sigma_p^2/(\sigma_p^2 + \sigma^2)$ .

Sin embargo, los modelos no son totalmente equivalentes, ya que el modelo de efectos aleatorios obliga a que la correlación sea positiva. El modelo CS y el modelo t-test/anova no lo hacen.

EDIT: También hay otras dos diferencias. En primer lugar, los modelos CS y de efectos aleatorios asumen la 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 verosimilitud, mientras que el anova se ajusta utilizando cuadrados medios; cuando todo está equilibrado coincidirán, pero no necesariamente en situaciones más complejas. Por último, yo desconfiaría de utilizar los valores F/df/p de los distintos ajustes como medidas del grado de concordancia de los modelos; para más detalles, véase el famoso artículo de Doug Bates sobre df. (FIN EDIT)

El problema con su R código es que no está especificando la estructura de correlación correctamente. Debe utilizar gls con el corCompSymm estructura de correlación.

Generar 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")

Así es como se ajustarían los modelos de efectos aleatorios y 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 varianza residual del modelo CS es:

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

Y están a la altura de 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 no suelen ajustarse con efectos aleatorios, sino simplemente especificando la estructura deseada; una excepción común es el modelo AR(1) + efecto aleatorio, que tiene un efecto aleatorio y una correlación AR(1) entre observaciones sobre el mismo efecto aleatorio.

EDIT2: Cuando ajusto 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 **

(El intercepto es diferente aquí porque con la codificación por defecto, no es la media de todos los sujetos sino la media del primer sujeto).

También es interesante señalar que los nuevos lme4 da los mismos resultados, pero ni siquiera intenta 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

Answer 2

También puede utilizar la función mixed en el paquete afex para devolver valores p con la aproximación df de Kenward-Roger, que devuelve valores p idénticos a los de 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")