Random coefficients models
Ching Wen, Su
2016/10/30
require(pacman)## Loading required package: pacmanpacman::p_load(ggplot2, lme4, nlme)Q7
#data
dta<- read.table("http://titan.ccunix.ccu.edu.tw/~psycfs/lmm/Data/lecithin.txt", h=T)
head(dta) Group Visit Recall Patient
1 P 1 20 P01
2 P 1 14 P02
3 P 1 7 P03
4 P 1 6 P04
5 P 1 9 P05
6 P 1 9 P06str(dta) 'data.frame': 235 obs. of 4 variables:
$ Group : Factor w/ 2 levels "L","P": 2 2 2 2 2 2 2 2 2 2 ...
$ Visit : int 1 1 1 1 1 1 1 1 1 1 ...
$ Recall : int 20 14 7 6 9 9 7 18 6 10 ...
$ Patient: Factor w/ 48 levels "0P2","L26","L27",..: 24 25 26 27 28 29 30 31 32 33 ...# plot
ggplot(dta, aes(Visit,Recall, color=Patient) ) +
facet_grid(.~Group)+
geom_line(aes(group = Patient)) +
geom_point(alpha = 0.5) +
theme(legend.position="none")+
labs(x = "Visit", y = "Recall score")
# plot by group
ggplot(dta, aes(Visit,Recall, group=Group, linetype=Group) ) +
stat_summary(fun.y = mean, geom = "line") +
stat_summary(fun.y = mean, geom = "point") +
labs(x = "Visit", y = "Recall score")
#model
options(contrasts = c("contr.sum", "contr.poly"))
summary(m0 <- lmer(Recall ~ Visit + Group + (1 | Patient), data = dta)) Linear mixed model fit by REML ['lmerMod']
Formula: Recall ~ Visit + Group + (1 | Patient)
Data: dta
REML criterion at convergence: 1272.4
Scaled residuals:
Min 1Q Median 3Q Max
-2.7252 -0.5828 -0.1048 0.5776 3.1179
Random effects:
Groups Name Variance Std.Dev.
Patient (Intercept) 15.686 3.961
Residual 8.277 2.877
Number of obs: 235, groups: Patient, 48
Fixed effects:
Estimate Std. Error t value
(Intercept) 8.4357 0.7263 11.615
Visit 0.5064 0.1332 3.802
Group1 1.5087 0.6048 2.495
Correlation of Fixed Effects:
(Intr) Visit
Visit -0.554
Group1 0.062 0.00415.686/(15.686+8.277) [1] 0.6545925summary(m1 <- lmer(Recall ~ Visit + Group+Visit:Group + (1 | Patient), data = dta)) Linear mixed model fit by REML ['lmerMod']
Formula: Recall ~ Visit + Group + Visit:Group + (1 | Patient)
Data: dta
REML criterion at convergence: 1143.3
Scaled residuals:
Min 1Q Median 3Q Max
-2.5563 -0.5452 -0.0379 0.5492 3.2297
Random effects:
Groups Name Variance Std.Dev.
Patient (Intercept) 16.524 4.065
Residual 4.097 2.024
Number of obs: 235, groups: Patient, 48
Fixed effects:
Estimate Std. Error t value
(Intercept) 8.24808 0.66747 12.357
Visit 0.57904 0.09393 6.165
Group1 -2.41626 0.66747 -3.620
Visit:Group1 1.29823 0.09393 13.821
Correlation of Fixed Effects:
(Intr) Visit Group1
Visit -0.425
Group1 0.073 -0.020
Visit:Grop1 -0.020 0.055 -0.42516.524/(16.524+4.97) [1] 0.7687727Q6
# data
dta<-sleepstudy
head(dta) Reaction Days Subject
1 249.5600 0 308
2 258.7047 1 308
3 250.8006 2 308
4 321.4398 3 308
5 356.8519 4 308
6 414.6901 5 308str(dta) 'data.frame': 180 obs. of 3 variables:
$ Reaction: num 250 259 251 321 357 ...
$ Days : num 0 1 2 3 4 5 6 7 8 9 ...
$ Subject : Factor w/ 18 levels "308","309","310",..: 1 1 1 1 1 1 1 1 1 1 ...dta$Days <- factor(dta$Days)
# plot
ggplot(dta, aes(Days, Reaction, color=Subject) ) +
stat_smooth(aes(group = Subject), method = "lm", se = F)+
geom_point(alpha = 0.5) +
theme_bw()+
labs(x = "Days", y = "Reaction")
# model
options(contrasts = c("contr.sum", "contr.poly"))
summary(m0 <- lmer(Reaction ~ Days + (1 | Subject), data = dta)) Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (1 | Subject)
Data: dta
REML criterion at convergence: 1734.1
Scaled residuals:
Min 1Q Median 3Q Max
-3.3473 -0.5293 0.0317 0.5328 4.2570
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 1375.5 37.09
Residual 987.6 31.43
Number of obs: 180, groups: Subject, 18
Fixed effects:
Estimate Std. Error t value
(Intercept) 298.508 9.050 32.98
Days1 -41.856 7.027 -5.96
Days2 -34.012 7.027 -4.84
Days3 -33.146 7.027 -4.72
Days4 -15.516 7.027 -2.21
Days5 -9.858 7.027 -1.40
Days6 10.011 7.027 1.42
Days7 13.670 7.027 1.95
Days8 20.243 7.027 2.88
Days9 38.122 7.027 5.42
Correlation of Fixed Effects:
(Intr) Days1 Days2 Days3 Days4 Days5 Days6 Days7 Days8
Days1 0.000
Days2 0.000 -0.111
Days3 0.000 -0.111 -0.111
Days4 0.000 -0.111 -0.111 -0.111
Days5 0.000 -0.111 -0.111 -0.111 -0.111
Days6 0.000 -0.111 -0.111 -0.111 -0.111 -0.111
Days7 0.000 -0.111 -0.111 -0.111 -0.111 -0.111 -0.111
Days8 0.000 -0.111 -0.111 -0.111 -0.111 -0.111 -0.111 -0.111
Days9 0.000 -0.111 -0.111 -0.111 -0.111 -0.111 -0.111 -0.111 -0.111Q5
# data
dta<-read.table("http://titan.ccunix.ccu.edu.tw/~psycfs/lmm/Data/alcohol_use.txt", h=T)
head(dta) sid coa sex age14 alcuse peer cpeer ccoa
1 1 1 0 0 1.73205 1.26491 0.24691 0.549
2 1 1 0 1 2.00000 1.26491 0.24691 0.549
3 1 1 0 2 2.00000 1.26491 0.24691 0.549
4 2 1 1 0 0.00000 0.89443 -0.12357 0.549
5 2 1 1 1 0.00000 0.89443 -0.12357 0.549
6 2 1 1 2 1.00000 0.89443 -0.12357 0.549str(dta) 'data.frame': 246 obs. of 8 variables:
$ sid : int 1 1 1 2 2 2 3 3 3 4 ...
$ coa : int 1 1 1 1 1 1 1 1 1 1 ...
$ sex : int 0 0 0 1 1 1 1 1 1 1 ...
$ age14 : int 0 1 2 0 1 2 0 1 2 0 ...
$ alcuse: num 1.73 2 2 0 0 ...
$ peer : num 1.265 1.265 1.265 0.894 0.894 ...
$ cpeer : num 0.247 0.247 0.247 -0.124 -0.124 ...
$ ccoa : num 0.549 0.549 0.549 0.549 0.549 0.549 0.549 0.549 0.549 0.549 ...dta$sid<-factor(dta$sid)
# plot
ggplot(dta, aes(age14, alcuse, color=sid) ) +
stat_smooth(aes(group = sid), method = "lm", se = F)+
geom_point(alpha = 0.5) +
theme(legend.position="none")+
labs(x = "ages-14", y = "alcohol use")
ggplot(dta, aes(age14, alcuse, group=sid, color=coa) ) +
facet_grid(.~coa)+
stat_smooth(aes(group = sid), method = "lm", se = F)+
geom_point(alpha = 0.5) +
theme(legend.position="none")+
labs(x = "ages-14", y = "alcohol use")
# model
options(contrasts = c("contr.sum", "contr.poly"))
summary(m0 <- lmer(alcuse ~ coa + peer + age14 + peer:age14 + (1 | sid), data = dta)) Linear mixed model fit by REML ['lmerMod']
Formula: alcuse ~ coa + peer + age14 + peer:age14 + (1 | sid)
Data: dta
REML criterion at convergence: 621.6
Scaled residuals:
Min 1Q Median 3Q Max
-2.30588 -0.66394 -0.05233 0.55906 2.60999
Random effects:
Groups Name Variance Std.Dev.
sid (Intercept) 0.3377 0.5811
Residual 0.4823 0.6945
Number of obs: 246, groups: sid, 82
Fixed effects:
Estimate Std. Error t value
(Intercept) -0.31177 0.17225 -1.810
coa 0.56513 0.15878 3.559
peer 0.69586 0.13219 5.264
age14 0.42469 0.09346 4.544
peer:age14 -0.15138 0.07481 -2.024
Correlation of Fixed Effects:
(Intr) coa peer age14
coa -0.312
peer -0.725 -0.134
age14 -0.543 0.000 0.461
peer:age14 0.442 0.000 -0.566 -0.814summary(m1 <- lmer(alcuse ~ coa + peer + age14 + peer:age14 + (1 + age14 | sid), data = dta)) Linear mixed model fit by REML ['lmerMod']
Formula: alcuse ~ coa + peer + age14 + peer:age14 + (1 + age14 | sid)
Data: dta
REML criterion at convergence: 603.7
Scaled residuals:
Min 1Q Median 3Q Max
-2.59995 -0.40432 -0.07739 0.44372 2.27436
Random effects:
Groups Name Variance Std.Dev. Corr
sid (Intercept) 0.2595 0.5094
age14 0.1469 0.3832 -0.05
Residual 0.3373 0.5808
Number of obs: 246, groups: sid, 82
Fixed effects:
Estimate Std. Error t value
(Intercept) -0.31382 0.14871 -2.110
coa 0.57120 0.14898 3.834
peer 0.69518 0.11322 6.140
age14 0.42469 0.10690 3.973
peer:age14 -0.15138 0.08556 -1.769
Correlation of Fixed Effects:
(Intr) coa peer age14
coa -0.339
peer -0.708 -0.146
age14 -0.408 0.000 0.350
peer:age14 0.332 0.000 -0.429 -0.814Q3
# data
dta<-read.csv("visual_search.csv", h=T)
head(dta) Participant Dx Set.Size RT
1 42 Aphasic 1 1946.00
2 42 Aphasic 5 2371.17
3 42 Aphasic 15 4354.83
4 42 Aphasic 30 6825.60
5 44 Aphasic 1 1202.67
6 44 Aphasic 5 1796.67str(dta) 'data.frame': 132 obs. of 4 variables:
$ Participant: int 42 42 42 42 44 44 44 44 83 83 ...
$ Dx : Factor w/ 2 levels "Aphasic","Control": 1 1 1 1 1 1 1 1 1 1 ...
$ Set.Size : int 1 5 15 30 1 5 15 30 1 5 ...
$ RT : num 1946 2371 4355 6826 1203 ...dta$Participant<-factor(dta$Participant)
# plot by subject
ggplot(dta, aes(Set.Size, RT))+
geom_point(alpha = .3)+
geom_line(aes(group = Participant)) +
stat_smooth(method="lm", col = "gray") +
facet_grid(. ~ Dx)+
labs(x = "Set Size", y = "Reaction Time(ms)") 
# model
summary(m0 <- lmer(RT ~ Set.Size*Dx + ( Set.Size | Participant), data=dta)) Linear mixed model fit by REML ['lmerMod']
Formula: RT ~ Set.Size * Dx + (Set.Size | Participant)
Data: dta
REML criterion at convergence: 2188.8
Scaled residuals:
Min 1Q Median 3Q Max
-3.7578 -0.3130 -0.0750 0.3117 6.1372
Random effects:
Groups Name Variance Std.Dev. Corr
Participant (Intercept) 657549.5 810.89
Set.Size 407.4 20.18 1.00
Residual 772433.6 878.88
Number of obs: 132, groups: Participant, 33
Fixed effects:
Estimate Std. Error t value
(Intercept) 1525.722 183.452 8.317
Set.Size 62.625 7.719 8.113
Dx1 553.027 183.452 3.015
Set.Size:Dx1 10.869 7.719 1.408
Correlation of Fixed Effects:
(Intr) Set.Sz Dx1
Set.Size -0.071
Dx1 0.091 -0.006
Set.Siz:Dx1 -0.006 0.091 -0.071# plot by group
ggplot(data = dta, aes(x = Set.Size, y = RT, group=Dx, shape = Dx, linetype = Dx)) +
stat_summary( fun.y = mean, geom = "point", size = 2) +
stat_summary( fun.data = mean_se, geom = "errorbar",
linetype = "solid", width = .5) +
stat_smooth(method="lm", se=F, color="black") +
labs(x = "Set Size", y = "Reaction Time(ms)", shape = "Group", linetype = "Group") +
theme(legend.position = c(.2,.8))
# pearson residuals plot
plot(m0, resid(.,type = "pearson")~fitted(.)|Dx,
layout=c(2,1), abline = 0,
xlab = "Fitted Values", ylab = "Pearson Residuals", pch = 20, cex = 0.8, aspect = 1.2)
Q1
# data
dta<- read.table("http://titan.ccunix.ccu.edu.tw/~psycfs/lmm/Data/reading_piat.txt", h=T)
head(dta) ID Wave Age_G Age PIAT
1 1 1 6.5 6.0000 18
2 1 2 8.5 8.3333 35
3 1 3 10.5 10.3333 59
4 2 1 6.5 6.0000 18
5 2 2 8.5 8.5000 25
6 2 3 10.5 10.5833 28str(dta) 'data.frame': 267 obs. of 5 variables:
$ ID : int 1 1 1 2 2 2 3 3 3 4 ...
$ Wave : int 1 2 3 1 2 3 1 2 3 1 ...
$ Age_G: num 6.5 8.5 10.5 6.5 8.5 10.5 6.5 8.5 10.5 6.5 ...
$ Age : num 6 8.33 10.33 6 8.5 ...
$ PIAT : int 18 35 59 18 25 28 18 23 32 18 ...dta$age65<- dta$Age_G-6.5
# plot
ggplot(dta, aes(x = age65, y = PIAT, group = ID))+
geom_point(alpha = 0.3)+
geom_line(size = 0.5)+
labs(x = "Age - 6.5", y = "PIAT Score")
# model
summary(m0 <- lmer(PIAT~age65+(1+age65|ID), data = dta)) Linear mixed model fit by REML ['lmerMod']
Formula: PIAT ~ age65 + (1 + age65 | ID)
Data: dta
REML criterion at convergence: 1819.8
Scaled residuals:
Min 1Q Median 3Q Max
-2.6539 -0.5413 -0.1497 0.3854 3.2811
Random effects:
Groups Name Variance Std.Dev. Corr
ID (Intercept) 11.427 3.380
age65 4.486 2.118 0.22
Residual 27.043 5.200
Number of obs: 267, groups: ID, 89
Fixed effects:
Estimate Std. Error t value
(Intercept) 21.1629 0.6177 34.26
age65 5.0309 0.2973 16.92
Correlation of Fixed Effects:
(Intr)
age65 -0.316# residual plot
plot(m0, pch = 20, cex = 0.5, col = "steelblue",
xlab = "Fitted values", ylab = "Pearson residuals")