如何在 R 中拟合随机效应模型，主题为随机? [英] How to fit a random effects model with Subject as random in R?
问题描述
给定以下形式的数据
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))
我想将分数建模为主题、条件和时间的函数.每个(人类)受试者的分数被测量了 3 次，用变量 Time 表示，所以我有重复测量.
I would like to model Score as a function of Subject, Condition and Time. Each (human) Subject's score was measured three times, indicated by the variable Time, so I have repeated measures.
如何在 R 中构建随机效应模型，其中主题效应拟合为随机?
How can I build in R a random effects model with Subject effects fitted as random?
附录:有人问我是如何生成这些数据的.你猜对了，数据是假的，因为一天很长.分数是时间加上随机噪声，处于条件 1 会为分数增加一分.作为典型的 Psych 设置，它具有指导意义.你有一项任务，人们的分数随着练习(时间)和提高分数的药物(条件==1)而提高.
ADDENDUM: It's been asked how I generated these data. You guessed it, the data are fake as the day is long. Score is time plus random noise and being in Condition 1 adds a point to Score. It's instructive as a typical Psych setup. You have a task where people's score gets better with practice (time) and a drug (condition==1) that enhances score.
为了讨论的目的，这里有一些更现实的数据.现在模拟参与者有一个随机的技能"级别，可以添加到他们的分数中.此外，因子现在是字符串.
Here are some more realistic data for the purposes of this discussion. Now simulated participants have a random "skill" level that is added to their scores. Also, the factors are now strings.
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))
推荐答案
使用 nlme 库...
using the nlme library...
回答您提出的问题，您可以使用以下代码创建随机拦截混合效应模型:
Answering your stated question, you can create a random intecept mixed effect model using the following code:
> 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 之间的相关性相同.在正常情况下(可能不是那些生成假数据的情况)，我们希望后者小于前者.自回归结构通常是更好的方法.
The intercept variance is basically 0, indicating no within subject effect, so this model is not capturing the between time relationship well. A random intercept model is rarely the type of model you want for a repeated measures design. A random intercept model assumes that the correlations between all time points are equal. i.e. the correlation between time 1 and time 2 is the same as between time 1 and time 3. Under normal circumstances (perhaps not those generating your fake data) we would expect the later to be less than the former. An auto regressive structure is usually a better way to go.
> 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，这似乎很奇怪.通常，时间点之间至少应该存在正相关.这些数据是如何生成的?
Your data is showing a -.596 between time point correlation, which seems odd. normally there should, at a minimum be a positive correlation between time points. How was this data generated?
附录:
根据您的新数据，我们知道数据生成过程相当于一个随机截距模型(尽管这对于纵向研究来说不是最现实的.可视化显示时间的影响似乎是相当线性的，因此我们将其视为数字变量应该会很舒服.
With your new data we know that the data generating process is equivalent to a random intercept model (though that is not the most realistic for a longitudinal study. The visualization shows that the effect of time seems to be fairly linear, so we should feel comfortable treating it as a numeric variable.
> 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)，以及显着的时间效应，表明两组都随着时间的推移而上升.支持情节，我们发现两组之间没有时间差异效应(交互作用).即斜率相同.
We see a significant Condition effect, indicating that the 'yes' condition tends to have higher scores (by about 1.7), and a significant time effect, indicating that both groups go up over time. Supporting the plot, we find no differential effect of time between the two groups (the interaction). i.e. the slopes are the same.
这篇关于如何在 R 中拟合随机效应模型，主题为随机?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
