r - 在 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 堆叠在一起。