Standardowy błąd efektów losowych w R (Lme4) vs Stata (Xtmixed)

Proszę wziąć pod uwagę te dane:

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

Dopasowujemy prosty model komponentów wariancyjnych. W R mamy:

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

Następnie produkujemy działkę gąsienicową:

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

Teraz pasujemy do tego samego modelu w Stata. Najpierw napisz do formatu Stata z R:

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

W Stacie

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

Wyjście zgadza się z wyjściem R (nie pokazano), a my próbujemy wyprodukować ten sam wykres gąsienicy:

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 używa innego standardowego błędu niż R. W rzeczywistości Stata używa 2.13, podczas gdy R używa 1.32.

Z tego, co mogę powiedzieć, pochodzi 1,32 w R.

> 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

chociaż nie mogę powiedzieć, że naprawdę rozumiem, co to robi. Czy ktoś może wyjaśnić?

I nie mam pojęcia, skąd pochodzi 2.13 ze Staty, z wyjątkiem tego, że jeśli zmienię metodę szacowania na maksymalne prawdopodobieństwo:

xtmixed g || id:, ml variance

.... to wydaje się używać 1,32 jako błędu standardowego i dawać takie same wyniki jak R ....

.... ale wtedy oszacowanie wariancji efektu losowego nie zgadza się już z R (35,04 vs 31,97).

Wygląda więc na to, że ma to coś wspólnego z ML vs REML: jeśli uruchomię REML w obu systemach, dane wyjściowe modelu są zgodne, ale standardowe błędy użyte na wykresach gąsienicowych się nie zgadzają, natomiast jeśli uruchomię REML w R i ML w Stata , wykresy gąsienicowe są zgodne, ale szacunki modelu nie.

Czy ktoś może wyjaśnić, co się dzieje?

Zgodnie z [XT]instrukcją dla Stata 11:

Standardowe błędy dla BLUP są obliczane na podstawie iteracyjnej techniki Batesa i Pinheiro (1998, pkt. 3.3) w celu oszacowania samych BLUP. Jeżeli REML dokonuje oszacowania, te standardowe błędy uwzględniają niepewność oszacowania$\beta$, natomiast w przypadku ML traktowane są standardowe błędy $\beta$jak wiadomo W związku z tym standardowe błędy BLUP opartych na REML będą zwykle większe.

Ponieważ standardowe błędy Stata ML pasują do standardowych błędów z R w twoim przykładzie, wydaje się, że R nie uwzględnia niepewności przy szacowaniu $\beta$. Czy to powinno, nie wiem.

Z twojego pytania wypróbowałeś REML zarówno w Stata, jak i R, a ML w Stata z REML w R. Jeśli spróbujesz ML w obu, powinieneś uzyskać te same wyniki w obu.