机器算法验证 - 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"]]

R中的毛毛虫图

现在我们在 Stata 中拟合相同的模型。首先从 R 写入 Stata 格式：

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

在Stata中

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。

据我所知，R 中的 1.32 来自

> 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

虽然我不能说我真的明白这是在做什么。有人可以解释吗？

而且我不知道来自 Stata 的 2.13 来自哪里，除了，如果我将估计方法更改为最大似然：

xtmixed g || id:, ml variance

....那么它似乎使用 1.32 作为标准误差并产生与 R 相同的结果....

在此处输入图像描述

.... 但随后随机效应方差的估计值不再与 R 一致（35.04 对 31.97）。

所以它似乎与 ML 与 REML 有关：如果我在两个系统中运行 REML，模型输出一致，但在毛毛虫图中使用的标准错误不一致，而如果我在 R 中运行 REML，在 Stata 中运行 ML ，毛虫图同意，但模型估计不同意。

谁能解释发生了什么？