如何在 nlme 与 lme4 中指定不同的随机效果？

问题描述 投票：0回答：2

我想使用 

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 和 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")
r lme4 mixed-models nlme
2个回答
10
投票

我一直想尝试解决这个问题有一段时间了。如果没有更多的工作，我认为我无法获得与 

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
“因子”，将整个数据集包含在单个块中：

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
类，我无法找到一种简单的方法来让每个斜率 (
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
0
投票

作为对 Ben Bolker 答案的补充，为了允许主体位置的相关截距和斜率，我们应该使用 

pdLogChol()
或 
pdSymm()
来构造一般正定对称方差分量矩阵，而不是 
pdIdent()
来强制恒定对角线（不同随机效应预测变量之间的方差相同）和零非对角线元素（不同随机效应预测变量之间没有相关性）。我们还可以简单地删除 
pdIdent()
中的所有 
pdBlocked(list())
，因为如果仅提供公式，则 
pdLogChol()
是默认方法（请参阅 
?reStruct
）。

m3 <- lme(rating ~ 1 + position, random = list(dummy = pdBlocked(list(
    pdLogChol(~ 0 + subject + position), pdIdent(~ 0 + comparison)))),
    data = d)
m3 <- lme(rating ~ 1 + position, random = list(dummy = pdBlocked(list(
    ~ 0 + subject + position, ~ 0 + comparison))),
    data = d)

尽管 

nlme
非常强大，包括随机效应、误差相关性和方差结构，但估计这些元素的倍数可能具有挑战性。指定随机效应已经允许同一组中的观察结果以特定模式关联。因此，在存在随机效应的情况下，残差的相关结构通常是不必要的。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.