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

Por favor considere 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")

Nos ajustamos a un modelo de componentes de varianza simple. En R tenemos:

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

Luego producimos un diagrama de oruga:

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

Ahora encajamos el mismo modelo en Stata. Primero escriba en formato Stata desde R:

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

En Stata

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

El resultado coincide con el resultado R (ninguno de los dos se muestra) e intentamos producir el mismo diagrama 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)

Clearty Stata está usando un error estándar diferente a 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 realmente entiendo lo que está haciendo esto. Alguien puede explicar?

Y no tengo idea de dónde proviene el 2.13 de Stata, excepto que, si cambio el método de estimación a la máxima probabilidad:

xtmixed g || id:, ml variance

.... entonces parece usar 1.32 como el error estándar y producir los mismos resultados que R ...

.... pero luego la estimación de la varianza del efecto aleatorio ya no concuerda con R (35.04 vs 31.97).

Por lo tanto, 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 las parcelas de oruga no están de acuerdo, mientras que si ejecuto REML en R y ML en Stata , las parcelas de oruga están de acuerdo, pero las estimaciones del modelo no.

Puede alguien explicar qué está pasando ?

De acuerdo con el [XT]manual de Stata 11:

$\beta$$\beta$

$\beta$

A partir de su pregunta, ha probado REML en Stata y R, y ML en Stata con REML en R. Si prueba ML en ambos, debería obtener los mismos resultados en ambos.