Covariance pattern models
Ching Wen, Su
2016/11/7
require(pacman)## Loading required package: pacmanpacman::p_load(tidyverse, lme4, nlme, magrittr)
#theme set
theme_set(theme_bw())Q2
#data
dta2<-read.csv("claudication.csv", h=T)
str(dta2) 'data.frame': 38 obs. of 8 variables:
$ Treatment: Factor w/ 2 levels "ACT","PBO": 1 1 1 1 1 1 1 1 1 1 ...
$ PID : Factor w/ 38 levels "P101","P102",..: 1 5 9 12 17 22 23 28 29 32 ...
$ Gender : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 2 2 2 2 ...
$ D0 : int 190 98 155 245 182 140 196 162 195 167 ...
$ D1 : int 212 137 145 228 205 138 185 176 232 187 ...
$ D2 : int 213 185 196 280 218 187 185 192 199 228 ...
$ D3 : int 195 215 189 274 194 195 227 230 185 192 ...
$ D4 : int 248 225 176 260 193 205 180 215 200 210 ...head(dta2) Treatment PID Gender D0 D1 D2 D3 D4
1 ACT P101 M 190 212 213 195 248
2 ACT P105 M 98 137 185 215 225
3 ACT P109 M 155 145 196 189 176
4 ACT P112 M 245 228 280 274 260
5 ACT P117 M 182 205 218 194 193
6 ACT P122 M 140 138 187 195 205dta2L <- gather(dta2,key=Month, value=Distance, D0:D4)
head(dta2L) Treatment PID Gender Month Distance
1 ACT P101 M D0 190
2 ACT P105 M D0 98
3 ACT P109 M D0 155
4 ACT P112 M D0 245
5 ACT P117 M D0 182
6 ACT P122 M D0 140str(dta2L) 'data.frame': 190 obs. of 5 variables:
$ Treatment: Factor w/ 2 levels "ACT","PBO": 1 1 1 1 1 1 1 1 1 1 ...
$ PID : Factor w/ 38 levels "P101","P102",..: 1 5 9 12 17 22 23 28 29 32 ...
$ Gender : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 2 2 2 2 ...
$ Month : chr "D0" "D0" "D0" "D0" ...
$ Distance : int 190 98 155 245 182 140 196 162 195 167 ...dta2L$Month <- as.factor(dta2L$Month)
dta2L$Month_s<-dta2L$Month
dta2L<-separate(dta2L, col=Month_s, into=c("Prefix", "Time"), sep = 1)
dta2L$Time<-as.numeric(dta2L$Time)
#plot
ggplot(data = dta2L, aes(x = Month, y = Distance, group=Treatment, color=Treatment)) +
facet_grid(.~Gender)+
geom_point(alpha = .3)+
stat_summary( fun.y = mean, geom = "point", size = 2) +
stat_summary( fun.y = mean, geom = "line") +
stat_summary( fun.data = mean_se, geom = "errorbar",
linetype = "solid", width = .5) +
labs(x = "Time(in months)", y = "Distance(in meters)") +
theme(legend.position = c(.15,.85))
#mean & sd
with(dta2,aggregate(cbind(D0,D1,D2,D3,D4)~Gender+Treatment,FUN=mean)) Gender Treatment D0 D1 D2 D3 D4
1 F ACT 180.3750 193.3750 207.8750 196.7500 208.5000
2 M ACT 163.1667 179.5000 195.5833 204.0833 207.6667
3 F PBO 165.5556 170.1111 175.7778 168.1111 171.2222
4 M PBO 178.3333 165.5556 173.3333 193.6667 194.2222with(dta2,aggregate(cbind(D0,D1,D2,D3,D4)~Gender+Treatment,FUN=sd)) Gender Treatment D0 D1 D2 D3 D4
1 F ACT 42.18306 33.10994 58.21006 41.65762 57.10892
2 M ACT 42.34669 34.51350 40.14623 28.40441 28.04001
3 F PBO 41.52743 39.13261 47.25404 33.43443 45.34252
4 M PBO 45.00556 45.53051 43.68352 55.37824 47.35182#cor
cor(subset(dta2, Treatment == "ACT" & Gender == "F")[,-(1:3)]) D0 D1 D2 D3 D4
D0 1.0000000 0.8021966 0.8724733 0.8483852 0.8757233
D1 0.8021966 1.0000000 0.8645096 0.7073820 0.9054398
D2 0.8724733 0.8645096 1.0000000 0.8937510 0.9159571
D3 0.8483852 0.7073820 0.8937510 1.0000000 0.9251102
D4 0.8757233 0.9054398 0.9159571 0.9251102 1.0000000cor(subset(dta2, Treatment == "ACT" & Gender == "M")[,-(1:3)]) D0 D1 D2 D3 D4
D0 1.0000000 0.8624819 0.8052568 0.5928303 0.4068980
D1 0.8624819 1.0000000 0.6464277 0.3498338 0.4496799
D2 0.8052568 0.6464277 1.0000000 0.6704930 0.6323603
D3 0.5928303 0.3498338 0.6704930 1.0000000 0.5941441
D4 0.4068980 0.4496799 0.6323603 0.5941441 1.0000000cor(subset(dta2, Treatment == "PBO" & Gender == "F")[,-(1:3)]) D0 D1 D2 D3 D4
D0 1.0000000 0.9269134 0.6057892 0.8657560 0.7776916
D1 0.9269134 1.0000000 0.6826171 0.7707920 0.7579993
D2 0.6057892 0.6826171 1.0000000 0.3713199 0.4846551
D3 0.8657560 0.7707920 0.3713199 1.0000000 0.7931043
D4 0.7776916 0.7579993 0.4846551 0.7931043 1.0000000cor(subset(dta2, Treatment == "PBO" & Gender == "M")[,-(1:3)]) D0 D1 D2 D3 D4
D0 1.0000000 0.9271229 0.8530641 0.8833606 0.9202039
D1 0.9271229 1.0000000 0.7036016 0.7846666 0.9031326
D2 0.8530641 0.7036016 1.0000000 0.8292800 0.8073104
D3 0.8833606 0.7846666 0.8292800 1.0000000 0.9197535
D4 0.9202039 0.9031326 0.8073104 0.9197535 1.0000000#model
summary(m0<-gls(Distance ~ Time + Gender * Treatment + Time:Treatment,
correlation = corSymm(form = ~1|PID),
weights = varIdent(form = ~1|Month), data=dta2L)) Generalized least squares fit by REML
Model: Distance ~ Time + Gender * Treatment + Time:Treatment
Data: dta2L
AIC BIC logLik
1781.019 1848.533 -869.5095
Correlation Structure: General
Formula: ~1 | PID
Parameter estimate(s):
Correlation:
1 2 3 4
2 0.876
3 0.774 0.721
4 0.777 0.640 0.672
5 0.750 0.730 0.704 0.846
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | Month
Parameter estimates:
D0 D1 D2 D3 D4
1.0000000 0.8838348 1.0627783 0.9500921 1.0337160
Coefficients:
Value Std.Error t-value p-value
(Intercept) 177.00300 12.719156 13.916254 0.0000
Time 10.28326 1.679140 6.124126 0.0000
GenderM 0.01961 15.640605 0.001254 0.9990
TreatmentPBO -11.14603 17.575891 -0.634166 0.5268
GenderM:TreatmentPBO 0.20352 22.484786 0.009051 0.9928
Time:TreatmentPBO -7.78950 2.439733 -3.192766 0.0017
Correlation:
(Intr) Time GendrM TrtPBO GM:TPB
Time -0.304
GenderM -0.738 0.000
TreatmentPBO -0.724 0.220 0.534
GenderM:TreatmentPBO 0.513 0.000 -0.696 -0.702
Time:TreatmentPBO 0.210 -0.688 0.000 -0.320 0.000
Standardized residuals:
Min Q1 Med Q3 Max
-1.82526759 -0.58373590 -0.06888869 0.53349861 3.70599149
Residual standard error: 43.29371
Degrees of freedom: 190 total; 184 residualQ3
#data
dta3<-read.table("http://titan.ccunix.ccu.edu.tw/~psycfs/lmm/Data/adas.txt", header = TRUE, na.strings = ".")
head(dta3) Treatment PID adas02 adas04 adas06 adas08 adas10 adas12
1 L 1 22 30 NA 33 28 30
2 L 5 34 35 46 37 31 35
3 L 8 40 41 41 46 52 48
4 L 12 24 NA 21 28 30 27
5 L 13 29 26 29 26 NA 36
6 L 15 31 36 41 46 52 57str(dta3) 'data.frame': 80 obs. of 8 variables:
$ Treatment: Factor w/ 3 levels "H","L","P": 2 2 2 2 2 2 2 2 2 2 ...
$ PID : int 1 5 8 12 13 15 19 21 24 28 ...
$ adas02 : int 22 34 40 24 29 31 22 43 18 25 ...
$ adas04 : int 30 35 41 NA 26 36 27 49 28 24 ...
$ adas06 : int NA 46 41 21 29 41 28 42 29 27 ...
$ adas08 : int 33 37 46 28 26 46 24 48 NA 18 ...
$ adas10 : int 28 31 52 30 NA 52 27 48 25 21 ...
$ adas12 : int 30 35 48 27 36 57 28 46 28 22 ...dta3$Treatment <- factor(dta3$Treatment, levels = c("P", "L", "H"))
dta3$PID<-as.factor(dta3$PID)
dta3L<-gather(dta3,key=Month, value=adas, adas02:adas12, na.rm = T)
head(dta3L) Treatment PID Month adas
1 L 1 adas02 22
2 L 5 adas02 34
3 L 8 adas02 40
4 L 12 adas02 24
5 L 13 adas02 29
6 L 15 adas02 31str(dta3L) 'data.frame': 454 obs. of 4 variables:
$ Treatment: Factor w/ 3 levels "P","L","H": 2 2 2 2 2 2 2 2 2 2 ...
$ PID : Factor w/ 80 levels "1","2","3","4",..: 1 5 8 12 13 15 19 21 24 28 ...
$ Month : chr "adas02" "adas02" "adas02" "adas02" ...
$ adas : int 22 34 40 24 29 31 22 43 18 25 ...#plot
ggplot(dta3L, aes(x=Month, y=adas, group=Treatment, color=Treatment))+
geom_point(alpha=.3)+
stat_summary( fun.y = mean, geom = "point", size = 2) +
stat_summary( fun.y = mean, geom = "line") +
stat_summary( fun.data = mean_se, geom = "errorbar",
linetype = "solid", width = .5) +
labs(x = "Month", y = "ADAS") 
#mean & sd
with(dta3,aggregate(cbind(adas02,adas04,adas06,adas08,adas10,adas12)~Treatment,FUN=mean)) Treatment adas02 adas04 adas06 adas08 adas10 adas12
1 P 29.61905 33.42857 34.85714 36.38095 37.57143 39.04762
2 L 32.82353 35.11765 37.64706 36.41176 38.35294 39.17647
3 H 30.22222 32.44444 33.55556 33.33333 33.72222 34.27778with(dta3,aggregate(cbind(adas02,adas04,adas06,adas08,adas10,adas12)~Treatment,FUN=sd)) Treatment adas02 adas04 adas06 adas08 adas10 adas12
1 P 6.778467 6.415383 6.428730 6.344101 7.032577 8.800433
2 L 7.510287 8.637725 9.089425 11.090868 11.191186 10.765891
3 H 8.558144 7.326200 8.555853 9.055385 7.442371 7.102655#model
summary(m0<-lme(adas ~ Treatment,random = ~1|Month, data = dta3L )) Linear mixed-effects model fit by REML
Data: dta3L
AIC BIC logLik
3279.055 3299.612 -1634.527
Random effects:
Formula: ~1 | Month
(Intercept) Residual
StdDev: 2.258585 8.835074
Fixed effects: adas ~ Treatment
Value Std.Error DF t-value p-value
(Intercept) 36.92824 1.1455550 446 32.23611 0.0000
TreatmentL -1.67993 0.9965295 446 -1.68578 0.0925
TreatmentH -2.61476 1.0137870 446 -2.57920 0.0102
Correlation:
(Intr) TrtmnL
TreatmentL -0.405
TreatmentH -0.398 0.457
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.0470165 -0.7501237 -0.1456095 0.6286689 3.2113600
Number of Observations: 454
Number of Groups: 6