如何在 nlme 和 lme4 中指定不同的随机效应？
How to specify different random effects in nlme vs. lme4?
我想使用 nlme::lme（底部的数据）在模型中指定不同的随机效应。随机效应是：1) intercept 和 position 在 subject 内变化； 2) intercept 在 comparison 内变化。使用 lme4::lmer 很简单：
lmer(rating ~ 1 + position +
(1 + position | subject) +
(1 | comparison), data=d)
> ...
Random effects:
Groups Name Std.Dev. Corr
comparison (Intercept) 0.31877
subject (Intercept) 0.63289
position 0.06254 -1.00
Residual 0.91458
...
但是，我想坚持使用 lme，因为我还想对自相关结构建模（position 是一个时间变量）。 如何使用lme实现与上面相同的效果？我在下面尝试嵌套效果，这不是我想要的。
lme(rating ~ 1 + position,
random = list( ~ 1 + position | subject,
~ 1 | comparison), data=d)
> ...
Random effects:
Formula: ~1 + position | subject
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.53817955 (Intr)
position 0.04847635 -1
Formula: ~1 | comparison %in% subject # NESTED :(
(Intercept) Residual
StdDev: 0.9707665 0.0002465237
...
注意: SO 和 CV 上有一些类似的问题 here, here, and here 但我要么不明白答案要么建议使用 lmer这不算在这里 ;)
示例中使用的数据
d <- structure(list(rating = c(2, 3, 4, 3, 2, 4, 4, 3, 2, 1, 3, 2,
2, 2, 4, 2, 4, 3, 2, 2, 3, 5, 3, 4, 4, 4, 3, 2, 3, 5, 4, 5, 2,
3, 4, 2, 4, 4, 1, 2, 4, 5, 4, 2, 3, 4, 3, 2, 2, 2, 4, 5, 4, 4,
5, 2, 3, 4, 3, 2), subject = structure(c(1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 6L), .Label = c("1", "2", "3", "4", "5",
"6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16",
"17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27",
"28", "29", "30", "31", "32", "33", "34", "35", "36", "37", "38",
"39", "40", "41", "42", "43", "44", "45", "46", "47", "48", "49",
"50", "51", "52", "53", "54", "55", "56", "57", "58", "59", "60",
"61", "62", "63"), class = "factor"), position = c(1, 2, 3, 4,
5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4,
5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4,
5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10), comparison = structure(c(1L,
7L, 9L, 8L, 3L, 4L, 10L, 2L, 5L, 6L, 2L, 6L, 4L, 5L, 8L, 10L,
7L, 3L, 1L, 9L, 3L, 9L, 10L, 1L, 5L, 7L, 6L, 8L, 2L, 4L, 4L,
2L, 8L, 6L, 7L, 5L, 1L, 10L, 9L, 3L, 5L, 10L, 6L, 3L, 2L, 9L,
4L, 1L, 8L, 7L, 6L, 5L, 2L, 10L, 4L, 3L, 8L, 9L, 7L, 1L), contrasts = structure(c(1,
0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1, 0,
0, 1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0,
0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0,
0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 0, -1, 0,
0, 0, 0, 0, 0, 0, 0, 1, -1), .Dim = c(10L, 9L), .Dimnames = list(
c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10"), NULL)), .Label = c("1",
"2", "3", "4", "5", "6", "7", "8", "9", "10"), class = "factor")), .Names = c("rating",
"subject", "position", "comparison"), row.names = c(1L, 2L, 3L,
4L, 5L, 6L, 7L, 8L, 9L, 10L, 111L, 112L, 113L, 114L, 115L, 116L,
117L, 118L, 119L, 120L, 221L, 222L, 223L, 224L, 225L, 226L, 227L,
228L, 229L, 230L, 331L, 332L, 333L, 334L, 335L, 336L, 337L, 338L,
339L, 340L, 441L, 442L, 443L, 444L, 445L, 446L, 447L, 448L, 449L,
450L, 551L, 552L, 553L, 554L, 555L, 556L, 557L, 558L, 559L, 560L
), class = "data.frame")
一段时间以来，我一直想弄清楚这个问题。如果没有更多的工作，我认为我无法获得与 lme4 中完全相同的模型，但我可以接近。
## source("SO36643713.dat")
library(nlme)
library(lme4)
这是您想要的模型，具有 subject 的完整随机斜率项（相关斜率和截距）和 comparison 的随机截距：
m1 <- lmer(rating ~ 1 + position +
(1 + position | subject) +
(1 | comparison), data=d)
这是我可以弄清楚如何在 lme 中复制的那个：独立的截距和斜率。 （我不是特别喜欢这些模型，但它们作为人们简化过于复杂的随机效应模型的一种方式相当普遍。）
m2 <- lmer(rating ~ 1 + position +
(1 + position || subject) +
(1 | comparison), data=d)
结果：
VarCorr(m2)
## Groups Name Std.Dev.
## comparison (Intercept) 0.28115
## subject position 0.00000
## subject.1 (Intercept) 0.28015
## Residual 0.93905
对于这个特定的数据集，估计随机斜率无论如何都具有零方差。
现在让我们为 lme 设置它。关键 (???) 见解是 pdBlocked() 矩阵中的所有项必须 嵌套在同一分组变量 中。例如，Pinheiro 和 Bates 第 163ff 页上的交叉随机效应示例将块、块内的行和块内的列作为随机效应。由于没有 comparison 和 subject 都嵌套在其中的分组因子，我将组成一个 dummy "factor" 将整个数据集包含在一个单块：
d$dummy <- factor(1)
现在我们可以拟合模型了。
m3 <- lme(rating~1+position,
random=list(dummy =
pdBlocked(list(pdIdent(~subject-1),
pdIdent(~position:subject),
pdIdent(~comparison-1)))),
data=d)
我们在随机效应方差-协方差矩阵中有三个块：一个用于 subject，一个用于 position-by-subject 交互，一个用于 comparison。除了定义全新的 pdMat class，我想不出一个简单的方法让每个斜率 (position:subjectXX) 与其对应的截距 (subjectXX). （您可能认为可以使用 pdBlocked 结构来设置它，但我看不到任何方法可以将方差估计值限制为在 pdBlocked 对象内的多个块之间相同。）
结果几乎相同，尽管报告不同。
vv <- VarCorr(m3)
vv2 <- vv[c("subject1","position:subject1","comparison1","Residual"),]
storage.mode(vv2) <- "numeric"
print(vv2,digits=4)
Variance StdDev
subject1 7.849e-02 2.802e-01
position:subject1 4.681e-11 6.842e-06
comparison1 7.905e-02 2.812e-01
Residual 8.818e-01 9.390e-01
我想使用 nlme::lme（底部的数据）在模型中指定不同的随机效应。随机效应是：1) intercept 和 position 在 subject 内变化； 2) intercept 在 comparison 内变化。使用 lme4::lmer 很简单：
lmer(rating ~ 1 + position +
(1 + position | subject) +
(1 | comparison), data=d)
> ...
Random effects:
Groups Name Std.Dev. Corr
comparison (Intercept) 0.31877
subject (Intercept) 0.63289
position 0.06254 -1.00
Residual 0.91458
...
但是，我想坚持使用 lme，因为我还想对自相关结构建模（position 是一个时间变量）。 如何使用lme实现与上面相同的效果？我在下面尝试嵌套效果，这不是我想要的。
lme(rating ~ 1 + position,
random = list( ~ 1 + position | subject,
~ 1 | comparison), data=d)
> ...
Random effects:
Formula: ~1 + position | subject
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.53817955 (Intr)
position 0.04847635 -1
Formula: ~1 | comparison %in% subject # NESTED :(
(Intercept) Residual
StdDev: 0.9707665 0.0002465237
...
注意: SO 和 CV 上有一些类似的问题 here, here, and here 但我要么不明白答案要么建议使用 lmer这不算在这里 ;)
示例中使用的数据
d <- structure(list(rating = c(2, 3, 4, 3, 2, 4, 4, 3, 2, 1, 3, 2,
2, 2, 4, 2, 4, 3, 2, 2, 3, 5, 3, 4, 4, 4, 3, 2, 3, 5, 4, 5, 2,
3, 4, 2, 4, 4, 1, 2, 4, 5, 4, 2, 3, 4, 3, 2, 2, 2, 4, 5, 4, 4,
5, 2, 3, 4, 3, 2), subject = structure(c(1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 6L), .Label = c("1", "2", "3", "4", "5",
"6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16",
"17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27",
"28", "29", "30", "31", "32", "33", "34", "35", "36", "37", "38",
"39", "40", "41", "42", "43", "44", "45", "46", "47", "48", "49",
"50", "51", "52", "53", "54", "55", "56", "57", "58", "59", "60",
"61", "62", "63"), class = "factor"), position = c(1, 2, 3, 4,
5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4,
5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4,
5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10), comparison = structure(c(1L,
7L, 9L, 8L, 3L, 4L, 10L, 2L, 5L, 6L, 2L, 6L, 4L, 5L, 8L, 10L,
7L, 3L, 1L, 9L, 3L, 9L, 10L, 1L, 5L, 7L, 6L, 8L, 2L, 4L, 4L,
2L, 8L, 6L, 7L, 5L, 1L, 10L, 9L, 3L, 5L, 10L, 6L, 3L, 2L, 9L,
4L, 1L, 8L, 7L, 6L, 5L, 2L, 10L, 4L, 3L, 8L, 9L, 7L, 1L), contrasts = structure(c(1,
0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1, 0,
0, 1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0,
0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0,
0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 0, -1, 0,
0, 0, 0, 0, 0, 0, 0, 1, -1), .Dim = c(10L, 9L), .Dimnames = list(
c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10"), NULL)), .Label = c("1",
"2", "3", "4", "5", "6", "7", "8", "9", "10"), class = "factor")), .Names = c("rating",
"subject", "position", "comparison"), row.names = c(1L, 2L, 3L,
4L, 5L, 6L, 7L, 8L, 9L, 10L, 111L, 112L, 113L, 114L, 115L, 116L,
117L, 118L, 119L, 120L, 221L, 222L, 223L, 224L, 225L, 226L, 227L,
228L, 229L, 230L, 331L, 332L, 333L, 334L, 335L, 336L, 337L, 338L,
339L, 340L, 441L, 442L, 443L, 444L, 445L, 446L, 447L, 448L, 449L,
450L, 551L, 552L, 553L, 554L, 555L, 556L, 557L, 558L, 559L, 560L
), class = "data.frame")
一段时间以来，我一直想弄清楚这个问题。如果没有更多的工作，我认为我无法获得与 lme4 中完全相同的模型，但我可以接近。
## source("SO36643713.dat")
library(nlme)
library(lme4)
这是您想要的模型，具有 subject 的完整随机斜率项（相关斜率和截距）和 comparison 的随机截距：
m1 <- lmer(rating ~ 1 + position +
(1 + position | subject) +
(1 | comparison), data=d)
这是我可以弄清楚如何在 lme 中复制的那个：独立的截距和斜率。 （我不是特别喜欢这些模型，但它们作为人们简化过于复杂的随机效应模型的一种方式相当普遍。）
m2 <- lmer(rating ~ 1 + position +
(1 + position || subject) +
(1 | comparison), data=d)
结果：
VarCorr(m2)
## Groups Name Std.Dev.
## comparison (Intercept) 0.28115
## subject position 0.00000
## subject.1 (Intercept) 0.28015
## Residual 0.93905
对于这个特定的数据集，估计随机斜率无论如何都具有零方差。
现在让我们为 lme 设置它。关键 (???) 见解是 pdBlocked() 矩阵中的所有项必须 嵌套在同一分组变量 中。例如，Pinheiro 和 Bates 第 163ff 页上的交叉随机效应示例将块、块内的行和块内的列作为随机效应。由于没有 comparison 和 subject 都嵌套在其中的分组因子，我将组成一个 dummy "factor" 将整个数据集包含在一个单块：
d$dummy <- factor(1)
现在我们可以拟合模型了。
m3 <- lme(rating~1+position,
random=list(dummy =
pdBlocked(list(pdIdent(~subject-1),
pdIdent(~position:subject),
pdIdent(~comparison-1)))),
data=d)
我们在随机效应方差-协方差矩阵中有三个块：一个用于 subject，一个用于 position-by-subject 交互，一个用于 comparison。除了定义全新的 pdMat class，我想不出一个简单的方法让每个斜率 (position:subjectXX) 与其对应的截距 (subjectXX). （您可能认为可以使用 pdBlocked 结构来设置它，但我看不到任何方法可以将方差估计值限制为在 pdBlocked 对象内的多个块之间相同。）
结果几乎相同，尽管报告不同。
vv <- VarCorr(m3)
vv2 <- vv[c("subject1","position:subject1","comparison1","Residual"),]
storage.mode(vv2) <- "numeric"
print(vv2,digits=4)
Variance StdDev
subject1 7.849e-02 2.802e-01
position:subject1 4.681e-11 6.842e-06
comparison1 7.905e-02 2.812e-01
Residual 8.818e-01 9.390e-01