如何在R中将Subject设置为Random的随机效应模型? [英] 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))
我想将Score建模为Subject，Condition和Time的函数.每个(人类)受试者的得分都被测量了三次，用变量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设置，它具有启发性.您需要完成一项任务，使人们的评分随着实践(时间)的提高而提高，并且使用一种药物(condition == 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))
推荐答案
...
回答您提出的问题，您可以使用以下代码创建随机的incept混合效果模型:
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中将Subject设置为Random的随机效应模型?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
