在nlme中提取随机效果设计矩阵

Question

传统上，线性混合效应模型是按以下方式制定的。Ri = Xi×？+ Zi×bi +？i在哪里？代表估计的固定效应，Z代表随机效应。因此，X是经典设计矩阵。使用R，在使用nlme包中的lme拟合模型后，我希望能够提取这两个矩阵。例如，也可以在nlme软件包中找到的数据集“ Rails”包含在6条随机选择的铁路轨道上的三个超声行进时间的单独测量值。我可以通过以下方法为每个轨道设置截距固定效果和随机效果的简单模型。

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

X设计矩阵只是一个18x1矩阵（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拟合模型中提取此矩阵。

Answer 1

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

1 有点特定于此示例，因为只有一个随机效应。当你有多个随机效果时，你可以做一些更多的自动编程来将不同的 Z_i 堆叠在一起。

Answer 2

据我所知，Z矩阵没有存储在lme对象中的任何位置。最好的选择是在modelStruct$reStruct组件中（尝试names(modelfit); str(modelfit); sapply(modelfit,class)等探索），但据我所知，它并不存在。事实上，一些深入研究lme.default表明该Z矩阵实际上可能永远不会被显式构建；在内部lme似乎可以使用分组结构。你当然可以

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

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