Standardfehler von zufälligen Effekten in R (lme4) vs Stata (xtmixed)

Bitte beachten Sie diese Daten:

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

Wir passen ein einfaches Varianzkomponentenmodell an. In R haben wir:

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

Dann produzieren wir eine Raupenparzelle:

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

Jetzt passen wir das gleiche Modell in Stata. Schreiben Sie zuerst in das Stata-Format von R:

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

In Stata

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

Die Ausgabe stimmt mit der R-Ausgabe überein (nicht gezeigt), und wir versuchen, die gleiche Raupenzeichnung zu erzeugen:

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 verwendet einen anderen Standardfehler als R. Tatsächlich verwendet Stata 2.13, während R 1.32 verwendet.

Soweit ich das beurteilen kann, kommt der 1.32 in R von

> 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

obwohl ich nicht sagen kann, dass ich wirklich verstehe, was das tut. Kann das jemand erklären?

Und ich habe keine Ahnung, woher der 2.13 von Stata kommt, außer wenn ich die Schätzmethode auf maximale Wahrscheinlichkeit ändere:

xtmixed g || id:, ml variance

.... dann scheint es 1,32 als Standardfehler zu verwenden und die gleichen Ergebnisse wie R zu liefern ....

.... aber dann stimmt die Schätzung für die zufällige Effektvarianz nicht mehr mit R (35.04 vs 31.97) überein.

Es scheint also etwas mit ML vs REML zu tun zu haben: Wenn ich REML in beiden Systemen ausführe, stimmt die Modellausgabe überein, aber die in den Raupendiagrammen verwendeten Standardfehler stimmen nicht überein, wohingegen, wenn ich REML in R und ML in Stata ausführe stimmen die Raupenplots überein, die Modellschätzungen jedoch nicht.

Kann mir jemand erklären, was los ist?

Nach dem [XT]Handbuch für Stata 11:

$\beta$$\beta$

$\beta$

Aufgrund Ihrer Frage haben Sie REML in Stata und R und ML in Stata mit REML in R ausprobiert. Wenn Sie ML in beiden ausprobieren, sollten Sie in beiden die gleichen Ergebnisse erzielen.