Erro padrão de efeitos aleatórios em R (lme4) vs Stata (xtmixed)

Por favor, considere estes dados:

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 um modelo simples de componentes de variação. Em R temos:

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

Em seguida, produzimos um gráfico de lagarta:

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

Agora ajustamos o mesmo modelo no Stata. Primeiro, escreva para o formato Stata a partir de R:

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

Na Stata

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

A saída concorda com a saída R (nem mostrada) e tentamos produzir o mesmo gráfico de lagarta:

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)

O Clearty Stata está usando um erro padrão diferente de R. Na verdade, o Stata está usando 2.13 enquanto R está usando 1.32.

Pelo que sei, o 1,32 em R vem 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

embora eu não possa dizer que realmente entendo o que isso está fazendo. Alguém pode explicar?

E não tenho idéia de onde o 2,13 do Stata é proveniente, exceto que, se eu mudar o método de estimativa para a máxima probabilidade:

xtmixed g || id:, ml variance

.... então parece usar 1,32 como erro padrão e produzir os mesmos resultados que R ....

.... mas a estimativa para a variação de efeito aleatório não concorda mais com R (35,04 vs 31,97).

Portanto, parece ter algo a ver com ML vs REML: se eu executar REML em ambos os sistemas, a saída do modelo concorda, mas os erros padrão usados ​​nas plotagens da lagarta não concordam, enquanto que se eu executar REML em R e ML em Stata , as parcelas da lagarta concordam, mas as estimativas do modelo não.

Alguém pode explicar o que está acontecendo?

De acordo com o [XT]manual para Stata 11:

$\beta$$\beta$

$\beta$

Com a sua pergunta, você tentou REML em Stata e R e ML em Stata com REML em R. Se você tentar ML em ambos, deverá obter os mesmos resultados em ambos.