Erreur standard des effets aléatoires dans R (lme4) vs Stata (xtmixed)

Veuillez considérer ces données:

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

Nous adaptons un modèle de composants de variance simple. Dans R, nous avons:

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

Ensuite, nous produisons un tracé de chenille:

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

Nous adaptons maintenant le même modèle dans Stata. Écrivez d'abord au format Stata à partir de R:

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

In Stata

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

La sortie est en accord avec la sortie R (non illustrée), et nous essayons de produire le même tracé de chenille:

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 utilise une erreur standard différente de R. En fait, Stata utilise 2.13 tandis que R utilise 1.32.

D'après ce que je peux dire, le 1,32 en R vient 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

même si je ne peux pas dire que je comprends vraiment ce que cela fait. Quelqu'un peut-il expliquer?

Et je n'ai aucune idée d'où vient le 2.13 de Stata, sauf que, si je change la méthode d'estimation au maximum de vraisemblance:

xtmixed g || id:, ml variance

.... alors il semble utiliser 1,32 comme erreur standard et produire les mêmes résultats que R ....

.... mais alors l'estimation de la variance à effet aléatoire n'est plus en accord avec R (35,04 vs 31,97).

Il semble donc que cela ait quelque chose à voir avec ML vs REML: si j'exécute REML dans les deux systèmes, la sortie du modèle est d'accord mais les erreurs standard utilisées dans les tracés de chenilles ne sont pas d'accord, tandis que si j'exécute REML dans R et ML dans Stata , les parcelles de chenilles sont d'accord, mais les estimations du modèle ne le sont pas.

Quelqu'un peut-il expliquer ce qui se passe?

Selon le [XT]manuel de Stata 11:

$\beta$$\beta$

$\beta$

D'après votre question, vous avez essayé REML dans Stata et R, et ML dans Stata avec REML dans R. Si vous essayez ML dans les deux, vous devriez obtenir les mêmes résultats dans les deux.