Error estándar de efectos aleatorios en R (lme4) frente a Stata (xtmixed)

Question

Tenga en cuenta estos datos:

dt.m <- structure(list(id = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12), occasion = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("g1", "g2"), class = "factor"),     g = c(12, 8, 22, 10, 10, 6, 8, 4, 14, 6, 2, 22, 12, 7, 24, 14, 8, 4, 5, 6, 14, 5, 5, 16)), .Names = c("id", "occasion", "g"), row.names = c(NA, -24L), class = "data.frame")

Ajustamos un modelo simple de componentes de la varianza. En R tenemos:

require(lme4)
fit.vc <- lmer( g ~ (1|id), data=dt.m )

Entonces producimos una trama de oruga:

rr1 <- ranef(fit.vc, postVar = TRUE)
dotplot(rr1, scales = list(x = list(relation = 'free')))[["id"]]

Ahora ajustamos el mismo modelo en Stata. Primero escriba al formato de Stata desde R:

require(foreign)
write.dta(dt.m, "dt.m.dta")

En Stata

use "dt.m.dta"
xtmixed g || id:, reml variance

La salida concuerda con la salida de R (no se muestra), e intentamos producir el mismo gráfico de oruga:

predict u_plus_e, residuals
predict u, reffects
gen e = u_plus_e – u
predict u_se, reses

egen tag = tag(id)
sort u
gen u_rank = sum(tag)

serrbar u u_se u_rank if tag==1, scale(1.96) yline(0)

Claramente Stata está usando un error estándar diferente al de R. De hecho Stata está usando 2.13 mientras que R está usando 1.32.

Por lo que puedo decir, el 1,32 en R viene de

> sqrt(attr(ranef(fit.vc, postVar = TRUE)[[1]], "postVar")[1, , ])
 [1] 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977 1.319977

aunque no puedo decir que entienda realmente lo que hace esto. ¿Puede alguien explicarlo?

Y no tengo ni idea de dónde viene el 2,13 de Stata, salvo que, si cambio el método de estimación a máxima verosimilitud:

xtmixed g || id:, ml variance

.... entonces parece que utiliza 1,32 como error estándar y produce los mismos resultados que R....

.... pero entonces la estimación de la varianza del efecto aleatorio ya no coincide con R (35,04 frente a 31,97).

Así que parece tener algo que ver con ML vs REML: Si ejecuto REML en ambos sistemas, la salida del modelo está de acuerdo pero los errores estándar utilizados en los gráficos de oruga no están de acuerdo, mientras que si ejecuto REML en R y ML en Stata, los gráficos de oruga están de acuerdo, pero las estimaciones del modelo no.

¿Alguien puede explicar qué está pasando?

Answer 1

Según el [XT] manual para Stata 11:

Los errores estándar de los BLUPs se calculan según la técnica iterativa de Bates y Pinheiro (1998). técnica iterativa de Bates y Pinheiro (1998, sec. 3.3) para estimar los BLUPs. Si la estimación se realiza mediante REML, estos errores estándar tienen en cuenta la incertidumbre en la estimación de $\beta$ mientras que en el caso de ML los errores estándar tratan $\beta$ como se sabe. Como tal, los errores estándar de BLUPs basados en REML suelen ser mayores.

Como los errores estándar de Stata ML coinciden con los errores estándar de R en su ejemplo, parece que R no está teniendo en cuenta la incertidumbre en la estimación $\beta$ . Si debería, no lo sé.

Por su pregunta, usted ha probado REML tanto en Stata como en R, y ML en Stata con REML en R. Si prueba ML en ambos, debería obtener los mismos resultados en ambos.