如何在R中将随机效应模型与Subject随机拟合？
Dan*_*ein 17 statistics modeling model 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))
我想将得分作为主题,条件和时间的函数来建模.每个(人)受试者的分数测量三次,由变量Time表示,所以我有重复测量.
如何构建R随机效应模型,主题效果随机拟合？
ADDENDUM:有人问我是如何生成这些数据的.你猜对了,数据是假的,因为这一天很长.得分是时间加上随机噪音,在条件1中加分得分.它作为典型的心理学设置具有指导意义.你有一个任务,通过练习(时间)和提高分数的药物(条件== 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))
使用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),并且具有显着的时间效应,表明两个组随着时间的推移而上升.支持该情节,我们发现两组之间没有时间差异(相互作用).即斜坡是相同的.
这不是您的问题的答案,但您可能会发现您的数据可视化信息.
library(ggplot2)
qplot(Time, Score, data = myDat, geom = "line",
group = Subject, colour = factor(Condition))
