Стандартная ошибка случайных эффектов в R (lme4) против Stata (xtmixed)

Пожалуйста, рассмотрите эти данные:

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

Мы подходим к простой модели компонентов дисперсии. В R имеем:

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

Затем мы производим гусеничный участок:

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

Теперь мы соответствуем той же модели в Stata. Сначала напишите в формат Stata из R:

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

В стате

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

Вывод согласуется с выводом R (не показан), и мы пытаемся создать тот же график гусеницы:

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 использует стандартную ошибку, отличную от R. На самом деле Stata использует 2.13, тогда как R использует 1.32.

Из того, что я могу сказать, 1,32 в 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

хотя я не могу сказать, что я действительно понимаю, что это делает. Может кто-нибудь объяснить?

И я понятия не имею, откуда взялся 2.13 от Stata, за исключением того, что если я изменю метод оценки на максимальную вероятность:

xtmixed g || id:, ml variance

.... тогда он, кажется, использует 1,32 в качестве стандартной ошибки и дает те же результаты, что и R ....

.... но тогда оценка для случайной дисперсии эффекта больше не согласуется с R (35,04 против 31,97).

Так что, похоже, что-то связано с ML против REML: если я запускаю REML в обеих системах, выходные данные модели согласуются, но стандартные ошибки, используемые на графиках гусеницы, не согласуются, тогда как если я запускаю REML в R и ML в Stata , гусеничные участки согласны, но модельные оценки - нет.

Кто-нибудь может объяснить, что происходит?

Согласно [XT]инструкции для Stata 11:

Стандартные ошибки для BLUP рассчитываются на основе итерационной методики Бейтса и Пинейро (1998, с. 3.3) для оценки самих BLUP. Если оценка выполняется REML, эти стандартные ошибки учитывают неопределенность в оценке , в то время как для ML стандартные ошибки трактуют как известные. Таким образом, стандартные ошибки BLUP на основе REML обычно будут больше.$\beta$$\beta$

Поскольку стандартные ошибки Stata ML соответствуют стандартным ошибкам из R в вашем примере, кажется, что R не учитывает неопределенность в оценке . Должно ли это быть, я не знаю.$\beta$

Исходя из вашего вопроса, вы пробовали REML в Stata и R, и ML в Stata с REML в R. Если вы попробуете ML в обоих, вы должны получить одинаковые результаты в обоих.