6

线性混合效应模型传统上按以下方式制定。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 拟合模型中提取该矩阵。

4

2 回答 2

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

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

于 2013-11-23T03:13:07.780 回答
2

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

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

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

于 2013-10-31T16:32:42.637 回答