1106_covar_model
#peijun
Refer to the questions here.
Q1 Weight Lifting R/W
data <- read.table("http://140.116.183.121/~sheu/lmm/Data/weightlifting.txt", header = T)
# fill in NA
data[data == "."] <- NA
# long
dataL <- gather(data, key = "Time", value = "Weight", 3:9) %>%
mutate(Program = as.factor(Program), Weight = as.numeric(Weight), Time_f = as.factor(Time))Warning: attributes are not identical across measure variables;
they will be droppeddataL <- dataL[complete.cases(dataL), ] #remove rows with NA
# covariance matrices
data %>%
filter( Program == "R") %>%
select( -(1:2)) %>%
var(na.rm = T) Day0 Day2 Day4 Day6 Day8 Day10 Day12
Day0 11.06667 10.86667 10.42222 10.06667 10.55556 11.82222 11.82222
Day2 10.86667 11.65556 10.48889 10.14444 11.05556 12.13333 12.30000
Day4 10.42222 10.48889 11.37778 10.64444 11.77778 12.46667 12.91111
Day6 10.06667 10.14444 10.64444 10.23333 11.16667 12.15556 12.43333
Day8 10.55556 11.05556 11.77778 11.16667 12.72222 13.66667 14.05556
Day10 11.82222 12.13333 12.46667 12.15556 13.66667 15.51111 15.73333
Day12 11.82222 12.30000 12.91111 12.43333 14.05556 15.73333 16.45556data %>%
filter( Program == "W") %>%
select( -(1:2)) %>%
var(na.rm = T) Day0 Day2 Day4 Day6 Day8 Day10 Day12
Day0 11.192308 10.948718 10.506410 8.346154 9.750000 8.955128 8.923077
Day2 10.948718 11.269231 10.737179 8.224359 9.583333 8.923077 8.653846
Day4 10.506410 10.737179 11.064103 8.711538 10.000000 9.384615 9.064103
Day6 8.346154 8.224359 8.711538 7.756410 8.416667 7.852564 7.878205
Day8 9.750000 9.583333 10.000000 8.416667 10.166667 9.250000 9.250000
Day10 8.955128 8.923077 9.384615 7.852564 9.250000 8.974359 8.634615
Day12 8.923077 8.653846 9.064103 7.878205 9.250000 8.634615 9.230769# means with error bars
theme_set(theme_bw())
ggplot(dataL, aes(Time, Weight,
group = Program, shape = Program, linetype = Program)) +
stat_summary(fun.y = mean, geom = "line") +
stat_summary(fun.y = mean, geom = "point") +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
scale_shape_manual(values = c(1, 2)) +
labs(x = "Time (days since baseline)", y = "Mean Weight",
linetype = "Program", shape = "Program") +
theme(legend.justification = c(1, 1), legend.position = c(1, 1),
legend.background = element_rect(fill = "white", color = "black"))
# random intercept model
m0 <- gls(Weight ~ Program*Time_f, data = dataL)
m1 <- gls(Weight ~ Program*Time_f,
weights = varIdent(form = ~ 1 | Time_f*Program),
data = dataL)
m2 <- gls(Weight ~ Program*Time_f,
correlation = corSymm(form = ~ 1 | ID),
weights = varIdent(form = ~ 1 | Time_f*Program),
data = dataL)
anova(m0, m1, m2) Model df AIC BIC logLik Test L.Ratio p-value
m0 1 15 1260.1128 1311.354 -615.0564
m1 2 28 1278.9418 1374.593 -611.4709 1 vs 2 7.1710 0.8931
m2 3 49 848.4102 1015.799 -375.2051 2 vs 3 472.5317 <.0001# get m2
summary(m2)Generalized least squares fit by REML
Model: Weight ~ Program * Time_f
Data: dataL
AIC BIC logLik
848.4102 1015.799 -375.2051
Correlation Structure: General
Formula: ~1 | ID
Parameter estimate(s):
Correlation:
1 2 3 4 5 6
2 0.968
3 0.926 0.936
4 0.880 0.870 0.960
5 0.868 0.875 0.955 0.940
6 0.814 0.814 0.900 0.896 0.948
7 0.841 0.835 0.919 0.930 0.959 0.952
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | Time_f * Program
Parameter estimates:
Day0*R Day4*R Day6*R Day8*R Day10*R Day12*R Day2*R
1.0000000 1.1320878 1.1213941 1.2588616 1.3599243 1.3142947 1.0721688
Day0*W Day2*W Day4*W Day6*W Day8*W Day10*W Day12*W
0.9038511 0.9686848 1.0354647 0.8599033 1.0081359 0.9319670 0.9351870
Coefficients:
Value Std.Error t-value p-value
(Intercept) 79.92691 0.8227333 97.14803 0.0000
ProgramW 1.10221 1.0493577 1.05036 0.2947
Time_fDay10 0.94820 0.6360308 1.49081 0.1374
Time_fDay12 1.15376 0.5824175 1.98098 0.0488
Time_fDay2 0.88387 0.2272122 3.89007 0.0001
Time_fDay4 1.07717 0.3386506 3.18078 0.0017
Time_fDay6 1.34615 0.4266733 3.15498 0.0018
Time_fDay8 1.17608 0.5011235 2.34689 0.0198
ProgramW:Time_fDay10 0.79505 0.7532226 1.05553 0.2923
ProgramW:Time_fDay12 0.92893 0.6967468 1.33324 0.1838
ProgramW:Time_fDay2 -0.23153 0.2869728 -0.80682 0.4206
ProgramW:Time_fDay4 -0.13160 0.4423320 -0.29751 0.7663
ProgramW:Time_fDay6 0.19642 0.5242986 0.37464 0.7083
ProgramW:Time_fDay8 0.36008 0.6185880 0.58210 0.5611
Correlation:
(Intr) PrgrmW Tm_D10 Tm_D12 Tm_fD2 Tm_fD4 Tm_fD6
ProgramW -0.784
Time_fDay10 0.173 -0.136
Time_fDay12 0.159 -0.124 0.822
Time_fDay2 0.138 -0.108 0.174 0.150
Time_fDay4 0.130 -0.102 0.670 0.675 0.403
Time_fDay6 -0.011 0.009 0.647 0.709 0.128 0.790
Time_fDay8 0.171 -0.134 0.808 0.824 0.255 0.795 0.734
ProgramW:Time_fDay10 -0.146 0.034 -0.844 -0.694 -0.147 -0.566 -0.547
ProgramW:Time_fDay12 -0.133 0.033 -0.687 -0.836 -0.126 -0.564 -0.592
ProgramW:Time_fDay2 -0.109 0.135 -0.138 -0.119 -0.792 -0.319 -0.101
ProgramW:Time_fDay4 -0.100 0.136 -0.513 -0.517 -0.308 -0.766 -0.605
ProgramW:Time_fDay6 0.009 -0.129 -0.527 -0.577 -0.104 -0.643 -0.814
ProgramW:Time_fDay8 -0.138 0.091 -0.654 -0.668 -0.206 -0.644 -0.595
Tm_fD8 PW:T_D10 PW:T_D12 PW:T_D2 PW:T_D4 PW:T_D6
ProgramW
Time_fDay10
Time_fDay12
Time_fDay2
Time_fDay4
Time_fDay6
Time_fDay8
ProgramW:Time_fDay10 -0.682
ProgramW:Time_fDay12 -0.689 0.818
ProgramW:Time_fDay2 -0.202 0.166 0.147
ProgramW:Time_fDay4 -0.609 0.642 0.650 0.418
ProgramW:Time_fDay6 -0.597 0.659 0.717 0.130 0.757
ProgramW:Time_fDay8 -0.810 0.801 0.817 0.262 0.790 0.732
Standardized residuals:
Min Q1 Med Q3 Max
-2.35009552 -0.67840990 -0.01854778 0.61703103 2.63395367
Residual standard error: 3.309166
Degrees of freedom: 239 total; 225 residualQ3 Claudication & ACT/PBO
dta <- read.table("E:/201709/multilevel_analysis/1106/claudication.csv", header = T, sep = ",")
# long
dtaL <- gather(dta, key = "Month", value = "Distance", 4:8)
dtaL <- dtaL %>%
mutate(Month = as.factor(Month), Distance = as.numeric(Distance), Time = as.numeric(substr(dtaL$Month, 2 , 2)))
#plot
ggplot(dtaL, aes(Time, Distance, color = Treatment),
group = Treatment) +
geom_point(alpha = 0.5) +
stat_summary(fun.y = mean, geom = "line", aes(group=Treatment)) +
stat_summary(fun.y = mean, geom = "point", aes(group=Treatment)) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.1) +
facet_wrap(~ Gender) +
labs(x = "Time(in months)", y = "Distance(in meters)") +
theme(legend.position = c(.1, .8))
#calculate
aggregate(dta[, 4:8], list(dta$Gender, dta$Treatment), mean) Group.1 Group.2 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.2222aggregate(dta[, 4:8], list(dta$Gender, dta$Treatment), sd) Group.1 Group.2 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# correlation matrix
dta %>%
filter(Treatment == "ACT" & Gender == "F") %>%
select( -(1:3)) %>%
cor() 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.0000000dta %>%
filter(Treatment == "ACT" & Gender == "M") %>%
select( -(1:3)) %>%
cor() 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.0000000dta %>%
filter(Treatment == "PBO" & Gender == "F") %>%
select( -(1:3)) %>%
cor() 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.0000000dta %>%
filter(Treatment == "PBO" & Gender == "M") %>%
select( -(1:3)) %>%
cor() 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
dtaL <- dtaL %>% mutate(Month = as.numeric(Month))
m0 <- gls(Distance ~ Month + Gender*Treatment + Month:Treatment,
correlation = corSymm(form = ~ 1 | PID),
weights = varIdent(form = ~ 1 |Time),
data = dtaL)
summary(m0)Generalized least squares fit by REML
Model: Distance ~ Month + Gender * Treatment + Month:Treatment
Data: dtaL
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 | Time
Parameter estimates:
0 1 2 3 4
1.0000000 0.8838348 1.0627779 0.9500948 1.0337142
Coefficients:
Value Std.Error t-value p-value
(Intercept) 166.71981 13.326806 12.510110 0.0000
Month 10.28325 1.679138 6.124124 0.0000
GenderM 0.01944 15.640654 0.001243 0.9990
TreatmentPBO -3.35659 18.501998 -0.181418 0.8562
GenderM:TreatmentPBO 0.20371 22.484857 0.009060 0.9928
Month:TreatmentPBO -7.78949 2.439731 -3.192765 0.0017
Correlation:
(Intr) Month GendrM TrtPBO GM:TPB
Month -0.417
GenderM -0.704 0.000
TreatmentPBO -0.720 0.300 0.507
GenderM:TreatmentPBO 0.490 0.000 -0.696 -0.666
Month:TreatmentPBO 0.287 -0.688 0.000 -0.436 0.000
Standardized residuals:
Min Q1 Med Q3 Max
-1.82526479 -0.58373609 -0.06888758 0.53349971 3.70597932
Residual standard error: 43.29372
Degrees of freedom: 190 total; 184 residualQ4 Alzheimer (ADAS)
dta <- read.table("http://140.116.183.121/~sheu/lmm/Data/adas.txt", header = T)
dta[dta == "."] <- NA
# long
dtaL <- gather(dta, key = "Time", value = "Score", 3:8) %>%
mutate(Score = as.numeric(Score))Warning: attributes are not identical across measure variables;
they will be dropped#random intercept model
dtaL <- dtaL %>%
mutate(Month = as.numeric(substr(dtaL$Time, 5 , 6)))
dtaL <- dtaL[complete.cases(dtaL), ]
m0 <- gls(Score ~ Treatment*Month,
correlation = corSymm(form = ~ 1 | PID),
weights = varIdent(form = ~ 1 | Month),
data = dtaL)
summary(m0)Generalized least squares fit by REML
Model: Score ~ Treatment * Month
Data: dtaL
AIC BIC logLik
2758.115 2868.944 -1352.058
Correlation Structure: General
Formula: ~1 | PID
Parameter estimate(s):
Correlation:
1 2 3 4 5
2 0.888
3 0.778 0.824
4 0.803 0.791 0.865
5 0.793 0.760 0.739 0.868
6 0.766 0.683 0.718 0.810 0.916
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | Month
Parameter estimates:
2 4 8 10 12 6
1.0000000 0.9326570 1.0210262 1.0059898 1.0502727 0.9436511
Coefficients:
Value Std.Error t-value p-value
(Intercept) 32.17435 1.8357962 17.526099 0.0000
TreatmentL -0.68093 2.5491781 -0.267119 0.7895
TreatmentP -1.59821 2.4650451 -0.648350 0.5171
Month 0.44234 0.1266709 3.492017 0.0005
TreatmentL:Month 0.22198 0.1754074 1.265513 0.2063
TreatmentP:Month 0.52206 0.1710956 3.051296 0.0024
Correlation:
(Intr) TrtmnL TrtmnP Month TrtL:M
TreatmentL -0.720
TreatmentP -0.745 0.536
Month -0.455 0.328 0.339
TreatmentL:Month 0.328 -0.455 -0.245 -0.722
TreatmentP:Month 0.337 -0.242 -0.455 -0.740 0.535
Standardized residuals:
Min Q1 Med Q3 Max
-2.0259504 -0.7883728 -0.1975711 0.5308529 3.0726640
Residual standard error: 9.092343
Degrees of freedom: 454 total; 448 residualQ5 Grip strength with treat and gender
dta <- read.table("http://140.116.183.121/~sheu/lmm/Data/grip_strength.txt", header = T)
dta <- dta %>%
mutate(ID = as.factor(ID), Treat = as.factor(Treat))
#plot
ggplot(dta, aes(Time, Strength, color = ID), group = ID) +
geom_point(alpha = 0.5) +
geom_line() +
facet_wrap(~ Gender + Treat) +
labs(x = "Time", y = "Grip strength (in mmHg)")Warning: Removed 12 rows containing missing values (geom_point).Warning: Removed 12 rows containing missing values (geom_path).
# model
dta <- dta[complete.cases(dta), ]
m0 <- gls(Strength ~ Time + Baseline + Treat + Gender + Treat:Gender + Treat:Time + Gender:Time + Treat*Gender*Time + Baseline:Treat + Baseline:Gender + Baseline:Time,
correlation = corSymm(form = ~ 1 | ID),
weights = varIdent(form = ~ 1 | Time ), data = dta)
summary(m0)Generalized least squares fit by REML
Model: Strength ~ Time + Baseline + Treat + Gender + Treat:Gender + Treat:Time + Gender:Time + Treat * Gender * Time + Baseline:Treat + Baseline:Gender + Baseline:Time
Data: dta
AIC BIC logLik
1757.138 1814.309 -860.569
Correlation Structure: General
Formula: ~1 | ID
Parameter estimate(s):
Correlation:
1 2
2 0.414
3 0.332 0.797
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | Time
Parameter estimates:
1 2 3
1.000000 1.257335 1.475943
Coefficients:
Value Std.Error t-value p-value
(Intercept) 24.85306 13.699843 1.814113 0.0714
Time 5.59494 6.204543 0.901749 0.3684
Baseline 0.76187 0.130612 5.833103 0.0000
Treat2 1.17375 15.553154 0.075467 0.9399
GenderM 33.70866 25.766909 1.308215 0.1925
Treat2:GenderM 19.34235 22.197900 0.871360 0.3847
Time:Treat2 -2.30594 6.240323 -0.369522 0.7122
Time:GenderM -7.62308 8.577623 -0.888717 0.3754
Baseline:Treat2 -0.03765 0.129945 -0.289726 0.7724
Baseline:GenderM 0.02715 0.134783 0.201450 0.8406
Time:Baseline -0.00968 0.048761 -0.198482 0.8429
Time:Treat2:GenderM -0.44934 9.260537 -0.048522 0.9614
Correlation:
(Intr) Time Baseln Treat2 GendrM Tr2:GM Tm:Tr2
Time -0.603
Baseline -0.804 0.335
Treat2 -0.522 0.300 0.264
GenderM -0.151 -0.019 -0.043 -0.160
Treat2:GenderM 0.048 -0.201 0.212 0.062 -0.603
Time:Treat2 0.338 -0.567 -0.015 -0.542 -0.155 0.380
Time:GenderM -0.018 0.023 0.311 -0.192 -0.458 0.339 0.355
Baseline:Treat2 0.276 0.006 -0.346 -0.683 0.477 -0.637 -0.001
Baseline:GenderM 0.402 0.011 -0.504 0.176 -0.732 0.207 -0.010
Time:Baseline 0.414 -0.681 -0.504 -0.022 0.281 -0.005 0.051
Time:Treat2:GenderM -0.211 0.355 -0.010 0.363 0.267 -0.569 -0.672
Tm:GnM Bsl:T2 Bsl:GM Tm:Bsl
Time
Baseline
Treat2
GenderM
Treat2:GenderM
Time:Treat2
Time:GenderM
Baseline:Treat2 -0.003
Baseline:GenderM 0.002 -0.280
Time:Baseline -0.602 -0.003 -0.008
Time:Treat2:GenderM -0.593 0.003 0.005 0.005
Standardized residuals:
Min Q1 Med Q3 Max
-3.48804793 -0.49642940 0.03161909 0.51041336 3.24261761
Residual standard error: 25.30277
Degrees of freedom: 189 total; 177 residual