如何在 R 中将主题拟合为随机的随机效应模型？

How to fit a random effects model with Subject as random in R?

提问人：Dan Goldstein 提问时间：9/5/2009 最后编辑：Dan Goldstein 更新时间：8/3/2011 访问量：15750

问：

给定以下形式的数据

myDat = structure(list(Score = c(1.84, 2.24, 3.8, 2.3, 3.8, 4.55, 1.13, 
2.49, 3.74, 2.84, 3.3, 4.82, 1.74, 2.89, 3.39, 2.08, 3.99, 4.07, 
1.93, 2.39, 3.63, 2.55, 3.09, 4.76), Subject = c(1L, 1L, 1L, 
2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 6L, 7L, 
7L, 7L, 8L, 8L, 8L), Condition = c(0L, 0L, 0L, 1L, 1L, 1L, 0L, 
0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 
1L), Time = c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 
1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L)), .Names = c("Score", 
"Subject", "Condition", "Time"), class = "data.frame", row.names = c(NA, 
-24L))

我想将分数建模为主题、条件和时间的函数。每个（人类）受试者的分数被测量了三次，由变量时间表示，所以我重复了测量。

如何在 R 中构建一个随机效应模型，将主体效应拟合为随机？

ADDENDUM：有人问我是如何生成这些数据的。你猜对了，数据是假的，因为一天很长。分数是时间加上随机噪声，处于条件 1 会给分数加一分。作为典型的 Psych 设置，这很有启发性。你有一项任务，人们的分数会随着练习（时间）和提高分数的药物（条件==1）而变得更好。

以下是一些更现实的数据，用于本次讨论。现在，模拟参与者有一个随机的“技能”级别，该级别会添加到他们的分数中。此外，因子现在是字符串。

myDat = structure(list(Score = c(1.62, 2.18, 2.3, 3.46, 3.85, 4.7, 1.41, 
2.21, 3.32, 2.73, 3.34, 3.27, 2.14, 2.73, 2.74, 3.39, 3.59, 4.01, 
1.81, 1.83, 3.22, 3.64, 3.51, 4.26), Subject = structure(c(1L, 
1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 
6L, 7L, 7L, 7L, 8L, 8L, 8L), .Label = c("A", "B", "C", "D", "E", 
"F", "G", "H"), class = "factor"), Condition = structure(c(1L, 
1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 
2L, 1L, 1L, 1L, 2L, 2L, 2L), .Label = c("No", "Yes"), class = "factor"), 
    Time = structure(c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), .Label = c("1PM", 
    "2PM", "3PM"), class = "factor")), .Names = c("Score", "Subject", 
"Condition", "Time"), class = "data.frame", row.names = c(NA, 
-24L))

看它：

library(ggplot2)
qplot(Time, Score, data = myDat, geom = "line", group = Subject, colour = factor(Condition))
R 模型 统计建模

评论

1赞 Jonathan Chang 9/5/2009
您可以使用以下函数更简洁地构建数据框： myDat <- data.frame（Score = c（1.84， 2.24， 3.8， 2.3， 3.8， 4.55， 1.13， 2.49， 3.74， 2.84， 3.3， 4.82， 1.74， 2.89， 3.39， 2.08， 3.99， 4.07， 1.93， 2.39， 3.63， 2.55， 3.09， 4.76）， 主题 = c（1L， 1L， 1L， 1L， 1L， 1L， 1L， 2L， 2L， 2L， 2L， 2L， 2L， 2L， 3L， 3L， 3L， 3L， 3L， 3L， 4L， 4L， 4L， 4L， 4L， 4L）， 条件 = c（0L， 0L， 0L， 1L， 1L， 1L， 0L， 0L， 0L， 1L， 1L， 0L， 0L， 0L， 0L， 1L， 1L， 0L， 0L， 0L， 1L， 1L， 1L）， 时间 = c（1L， 2L， 3L， 1L， 2L， 3L， 1L， 2L， 3L， 1L， 2L， 3L， 1L， 3L， 1L， 1L， 1L， 1L， 1L， 1L， 1L， 1L， 1L， 1L， 1L， 1L， 1L， 1L， 1L， 1L， 1L， 1L， 1L， 1L， 1L， 1L， 1L， 1L， 1L， 1L、2L、3L、1L、2L、3L））data.frame
1赞 Eduardo Leoni 9/5/2009
@Chang。结构语法是您在 data.frame 中使用时得到的语法。dput
0赞 Jonathan Chang 9/5/2009
@Leoni。每天都能学到新东西！

答：

7赞 Ian Fellows 9/5/2009 #1

使用 NLME 库...

回答您陈述的问题，您可以使用以下代码创建随机 intecept 混合效应模型：

> library(nlme)
> m1 <- lme(Score ~ Condition + Time + Condition*Time,
+ data = myDat, random = ~ 1 | Subject)
> summary(m1)
Linear mixed-effects model fit by REML
 Data: myDat 
       AIC      BIC    logLik
  31.69207 37.66646 -9.846036

Random effects:
 Formula: ~1 | Subject
         (Intercept)  Residual
StdDev: 5.214638e-06 0.3151035

Fixed effects: Score ~ Condition + Time + Condition * Time 
                   Value Std.Error DF  t-value p-value
(Intercept)    0.6208333 0.2406643 14 2.579666  0.0218
Condition      0.7841667 0.3403507  6 2.303996  0.0608
Time           0.9900000 0.1114059 14 8.886423  0.0000
Condition:Time 0.0637500 0.1575517 14 0.404629  0.6919
 Correlation: 
               (Intr) Condtn Time  
Condition      -0.707              
Time           -0.926  0.655       
Condition:Time  0.655 -0.926 -0.707

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.5748794 -0.6704147  0.2069426  0.7467785  1.5153752 

Number of Observations: 24
Number of Groups: 8

截距方差基本为0，表示没有主体效应，因此该模型不能很好地捕捉时间关系。随机截距模型很少是重复测量设计所需的模型类型。随机截距模型假设所有时间点之间的相关性相等。即时间 1 和时间 2 之间的相关性与时间 1 和时间 3 之间的相关性相同。在正常情况下（也许不是那些生成虚假数据的人），我们预计后者会少于前者。自动回归结构通常是更好的方法。

> m2<-gls(Score ~ Condition + Time + Condition*Time,
+ data = myDat, correlation = corAR1(form = ~ Time | Subject))
> summary(m2)
Generalized least squares fit by REML
  Model: Score ~ Condition + Time + Condition * Time 
  Data: myDat 
       AIC      BIC    logLik
  25.45446 31.42886 -6.727232

Correlation Structure: AR(1)
 Formula: ~Time | Subject 
 Parameter estimate(s):
       Phi 
-0.5957973 

Coefficients:
                   Value Std.Error   t-value p-value
(Intercept)    0.6045402 0.1762743  3.429543  0.0027
Condition      0.8058448 0.2492895  3.232566  0.0042
Time           0.9900000 0.0845312 11.711652  0.0000
Condition:Time 0.0637500 0.1195452  0.533271  0.5997

 Correlation: 
               (Intr) Condtn Time  
Condition      -0.707              
Time           -0.959  0.678       
Condition:Time  0.678 -0.959 -0.707

Standardized residuals:
       Min         Q1        Med         Q3        Max 
-1.6850557 -0.6730898  0.2373639  0.8269703  1.5858942 

Residual standard error: 0.2976964 
Degrees of freedom: 24 total; 20 residual

您的数据显示时间点之间的相关性为 -.596，这似乎很奇怪。通常，时间点之间至少应该有正相关关系。这些数据是如何产生的？

补遗：

有了您的新数据，我们知道数据生成过程等同于随机截距模型（尽管对于纵向研究来说，这不是最现实的。可视化显示，时间的影响似乎是相当线性的，因此我们应该放心地将其视为数值变量。

> library(nlme)
> m1 <- lme(Score ~ Condition + as.numeric(Time) + Condition*as.numeric(Time),
+ data = myDat, random = ~ 1 | Subject)
> summary(m1)
Linear mixed-effects model fit by REML
 Data: myDat 
       AIC      BIC    logLik
  38.15055 44.12494 -13.07527

Random effects:
 Formula: ~1 | Subject
        (Intercept)  Residual
StdDev:   0.2457355 0.3173421

Fixed effects: Score ~ Condition + as.numeric(Time) + Condition * as.numeric(Time) 
                                  Value Std.Error DF   t-value p-value
(Intercept)                    1.142500 0.2717382 14  4.204415  0.0009
ConditionYes                   1.748333 0.3842958  6  4.549447  0.0039
as.numeric(Time)               0.575000 0.1121974 14  5.124898  0.0002
ConditionYes:as.numeric(Time) -0.197500 0.1586710 14 -1.244714  0.2337
 Correlation: 
                              (Intr) CndtnY as.(T)
ConditionYes                  -0.707              
as.numeric(Time)              -0.826  0.584       
ConditionYes:as.numeric(Time)  0.584 -0.826 -0.707

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.44560367 -0.65018585  0.01864079  0.52930925  1.40824838 

Number of Observations: 24
Number of Groups: 8

我们看到一个显着的条件效应，表明“是”条件往往具有更高的分数（约1.7），以及一个显着的时间效应，表明两组都随着时间的推移而上升。支持该图，我们发现两组之间没有时间差异（交互作用）。即斜坡是相同的。

评论

0赞 Dan Goldstein 9/5/2009
@Ian - 请参阅我添加到问题中的注释
0赞 Amyunimus 7/14/2014
这里的自由度从何而来？我知道 lme4 更复杂，但这应该很简单，不是吗？
0赞 Javier Hernando 2/8/2023
@Ian研究员，如果条件由 3 组组成，它会如何工作？Condition1:as.numeric(Time) & Condition2:as.numeric(Time) & Condition3:as.numeric(Time)
6赞 hadley 9/5/2009 #2

这不是您问题的答案，但您可能会发现这种数据可视化信息丰富。

library(ggplot2)
qplot(Time, Score, data = myDat, geom = "line", 
  group = Subject, colour = factor(Condition))

Data visulation

3赞 John 8/3/2011 #3

（使用 LME4 库） 这既适合你的主题效果是随机的，也适合你的随机效果所分组的变量。在这个模型中，随机效应是因受试者而异的截距。

m <- lmer( Score ~ Condition + Time + (1|Subject), data=myDat )

要查看随机效果，您可以使用

ranef(m)

正如 Ian Fellows 所提到的，您的数据也可能具有随机的条件和时间分量。您可以使用其他模型进行测试。在下面的一个中，条件、时间和截距允许因受试者而随机变化。它还评估了它们的相关性。

mi <- lmer( Score ~ Condition + Time + (Condition + Time|Subject), data=myDat )

并尝试

summary(mi)
ranef(mi)

您还可以在不与截距、条件和时间之间的交互以及许多其他模型的相关性的情况下对此进行测试，以查看哪个最适合您的数据和/或您的理论。你的问题有点含糊不清，但这几个命令应该能让你入门。

请注意，主体是您的分组因素，因此它是您将其其他效果随机归入的因素。它也不是你明确适合作为预测因子的东西。

评论

0赞 John 8/4/2011
我将保持答案不变，因为 ranef（） 通常与 mer 对象一起使用。不幸的是，在某些版本中它没有。在这种情况下，请尝试myModel@ranef。
0赞 Ben Bolker 10/8/2011
你能给我发一个例子（我是谷歌可搜索的）什么时候不能与对象一起使用吗......？？ranefmer
0赞 John 10/9/2011
当他们更新 LME4 并忘记让它工作时..:)...在我写这篇文章的时候，这种情况正在发生。

上一个：如何删除绘图轴和 R 中的绘图内容之间的绘图边距？

下一个：开发 Sweave 文档的建议