如何在nlme和lme4中指定不同的随机效果? [英] How to specify different random effects in nlme vs. lme4?

问题描述

我想使用nlme::lme在模型中指定不同的随机效果(底部数据).随机效应是:1)intercept和position在subject上变化； 2)intercept与comparison不同.使用lme4::lmer很简单:

I want to specify different random effects in a model using nlme::lme (data at the bottom). The random effects are: 1) intercept and position varies over subject; 2) intercept varies over comparison. This is straightforward using 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与上述相同?我在下面的尝试中嵌套了效果，这不是我想要的.

However, I want to stick to lme as I also want to model the autocorrelation structure (position is a time variable). How can I do the same as above using lme? My try below nests the effect, which is not what I want.

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也有一些类似的问题此处和在这里，但我要么不明白答案，要么建议使用lmer，在这里不算作;)

Note: There are some similar questions on SO and CV here, here, and here but I either did not understand the answer or the suggestion was to use lmer which not count here ;)

示例中使用的数据

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完全相同的模型，但是我可以接近.

I've been meaning to try to figure this out for a while. Without a lot more work I don't think I can get exactly the same model as in lme4, but I can get close.

## source("SO36643713.dat")
library(nlme)
library(lme4)

这是您想要的模型，对于subject具有完整的随机斜率项(相关的斜率和截距)，对于comparison具有随机截距:

This is the model you wanted, with a full random-slopes term for subject (correlated slope and intercept) and a random-intercept for comparison:

m1 <- lmer(rating ~ 1 + position + 
               (1 + position | subject) + 
               (1 | comparison), data=d)

我可以弄清楚如何在lme中进行复制:独立的截距和斜率. (我不特别喜欢这些模型，但是它们作为人们简化过于复杂的随机效应模型的一种相当普遍的用法.)

This is the one I can figure out how to replicate in lme: independent intercepts and slopes. (I don't like these models particularly, but they are in fairly common use as a way for people to simplify too-complex random-effects models.)

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 

对于这个特定的数据集，无论如何，随机斜率估计都具有零方差.

For this particular data set, the random slopes are estimated to have zero variance anyway.

现在，我们将其设置为lme.关键的洞察力是pdBlocked()矩阵中的所有术语必须嵌套在同一分组变量中.例如，Pinheiro和Bates的pp.163ff上的交叉随机效应示例具有块，块内的行和块内的列作为随机效果.由于没有将comparison和subject都嵌套在其中的分组因子，因此我将组成一个dummy因子"，该因子将整个数据集包含在一个块中:

Now let's set it up for lme. The key (???) insight is that all the terms inside a pdBlocked() matrix must be nested inside the same grouping variable. For example the crossed-random-effect example on pp. 163ff of Pinheiro and Bates has blocks, rows within blocks, and columns within blocks as the random effects. Since there is no grouping factor within which comparison and subject are both nested, I'm just going to make up a dummy "factor" that includes the whole data set in a single block:

d$dummy <- factor(1)

现在我们可以拟合模型了.

Now we can fit the model.

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

我们在随机效应方差-协方差矩阵中具有三个块:一个用于subject，一个用于position -by-c3>交互，另一个用于comparison.没有定义一个全新的pdMat类，我想不出一种简单的方法来允许每个斜率(position:subjectXX)与对应的截距(subjectXX)相关. (您可能会认为可以使用pdBlocked结构进行设置，但是我看不出有任何方法可以将方差估计值约束为在pdBlocked对象中的多个块之间相同.)

We have three blocks in the random-effects variance-covariance matrix: one for subject, one for the position-by-subject interaction, and one for comparison. Short of defining a brand-new pdMat class, I couldn't figure out an easy way to allow each slope (position:subjectXX) to be correlated with its corresponding intercept (subjectXX). (You might think you could set this up with a pdBlocked structure, but I don't see any way to constrain the variance estimates to be the same across multiple blocks within a pdBlocked object.)

尽管报告的结果不同，但结果几乎相同.

The results are pretty much identical, although they're reported differently.

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和lme4中指定不同的随机效果?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！