从nlme中提取随机效应设计矩阵

Question

传统上，线性混合效应模型的公式如下。Ri = Xi × β + Zi × bi + εi，其中β代表估计的固定效应，Z代表随机效应。因此，X是经典设计矩阵。使用R语言，我想在使用nlme包的lme拟合模型后能够提取这两个矩阵。例如，数据集“Rails”（也可以在nlme包中找到）包含对6条随机选定的铁路轨道进行三次超声波旅行时间测量的结果。我可以使用以下代码拟合一个简单模型，其中包括截距固定效应和每条轨道的随机效应。

library(nlme)
lmemodel<-lme(travel ~ 1, random = ~ 1 | Rail, data=Rail)

X设计矩阵只是一个由18个1组成的矩阵（6个轨道 * 3个测量值），可以通过以下方式轻松提取：

 model.matrix(lmemodel, data=Rail)
   (Intercept)
1            1
2            1
3            1
4            1
5            1
6            1
7            1
8            1
9            1
10           1
11           1
12           1
13           1
14           1
15           1
16           1
17           1
18           1
attr(,"assign")
[1] 0

我希望做的是提取随机效应设计矩阵Z。我意识到，如果我使用lme4包拟合相同的模型，可以通过以下方式完成：

library(lme4)
lmermodel<-lmer(travel ~ 1 + (1|Rail),data=Rail) 
t(lmermodel@Zt)  ##takes the transpose of lmermodel@Zt
lmermodel@X  ## extracts the X matrix

然而，我不知道如何从已安装的lme模型中提取此矩阵。

- iantist

---

最好使用 getME(lmermodel,"Z") 或 getME(lmermodel,"X") 而不是直接访问对象中的插槽，以防将来内部表示发生更改。 - Ben Bolker

Answer 1

model.matrix(formula(lmemodel$modelStruct$reStr)[[1]],data=lmemodel$data)

数字1在这个例子中比较特殊，因为只有一个随机效应。当你有多个随机效应时，可以进行一些更自动化的编程来堆叠不同的Z_i。

- user3023981

Answer 2

据我所见，在lme对象中没有存储Z矩阵的任何地方。最好的方法是在modelStruct$reStruct组件中查找（尝试使用names(modelfit); str(modelfit); sapply(modelfit,class)等方法进行探索），但据我所知，它并不存在于那里。实际上，对lme.default内部进行一些探索表明，Z矩阵可能从未被显式构建；在内部，lme似乎使用分组结构来工作。当然，您可以这样做。

Z <- model.matrix(~Rail-1,data=Rail)

但这可能不是你想要的……

- Ben Bolker