## 'data.frame': 38 obs. of 8 variables:
## $ Treatment: chr "ACT" "ACT" "ACT" "ACT" ...
## $ PID : chr "P101" "P105" "P109" "P112" ...
## $ Gender : chr "M" "M" "M" "M" ...
## $ 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 ...
## 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 205
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
dta1 <- mutate(dta,
Treatment = factor(Treatment),
Gender = factor(Gender),
PID = factor(PID))
dtaL <- gather(data=dta1, key=Time, value=Distance, D0:D4, factor_key=TRUE)
str(dtaL)
## '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 ...
## $ Time : Factor w/ 5 levels "D0","D1","D2",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Distance : int 190 98 155 245 182 140 196 162 195 167 ...
## Treatment PID Gender Time 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 140
library(ggplot2)
ggplot(dtaL, aes(Time, Distance,
group = Treatment,
shape = Treatment,
linetype = Treatment,
color = Treatment)) +
geom_point() +
stat_summary(fun = mean, geom = "point") +
stat_summary(fun = mean, geom = "line") +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.3) +
scale_shape_manual(values = c(1, 1)) +
facet_wrap( ~ Gender)+
labs(x = "TIME (in months)",
y = "Walking distance (in meters)",
linetype = "Treatment", shape = "Treatment") +
theme(legend.justification = c(7.7, 1),
legend.position = c(1, 1),
legend.background = element_rect(fill = "white",color = "black"))
dta1_nid <- dta1 %>% select(-2)
#mean
aggregate(.~ Treatment + Gender, data = dta1_nid, mean, na.rm=T)
## Treatment Gender D0 D1 D2 D3 D4
## 1 ACT F 180.3750 193.3750 207.8750 196.7500 208.5000
## 2 PBO F 165.5556 170.1111 175.7778 168.1111 171.2222
## 3 ACT M 163.1667 179.5000 195.5833 204.0833 207.6667
## 4 PBO M 178.3333 165.5556 173.3333 193.6667 194.2222
## Treatment Gender D0 D1 D2 D3 D4
## 1 ACT F 42.18306 33.10994 58.21006 41.65762 57.10892
## 2 PBO F 41.52743 39.13261 47.25404 33.43443 45.34252
## 3 ACT M 42.34669 34.51350 40.14623 28.40441 28.04001
## 4 PBO M 45.00556 45.53051 43.68352 55.37824 47.35182
## 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.0000000
## 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.0000000
## 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.0000000
## 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
#gls
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
summary(m1 <- gls(Distance ~ Time + Gender*Treatment + Time:Treatment,
weights = varIdent(form = ~ 1 | Time),
correlation = corSymm(form = ~ 1 | PID),
data = dtaL))
## Generalized least squares fit by REML
## Model: Distance ~ Time + Gender * Treatment + Time:Treatment
## Data: dtaL
## AIC BIC logLik
## 1749.506 1835.414 -847.7529
##
## Correlation Structure: General
## Formula: ~1 | PID
## Parameter estimate(s):
## Correlation:
## 1 2 3 4
## 2 0.878
## 3 0.781 0.718
## 4 0.770 0.635 0.673
## 5 0.746 0.736 0.715 0.837
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Time
## Parameter estimates:
## D0 D1 D2 D3 D4
## 1.0000000 0.8915599 1.0795057 0.9486381 1.0287427
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 170.03789 13.440748 12.650925 0.0000
## TimeD1 15.00000 4.607348 3.255669 0.0014
## TimeD2 30.45000 6.661777 4.570852 0.0000
## TimeD3 31.10000 6.377944 4.876180 0.0000
## TimeD4 37.95000 6.955426 5.456172 0.0000
## GenderM 0.02018 15.640507 0.001290 0.9990
## TreatmentPBO 1.79485 18.675206 0.096109 0.9235
## GenderM:TreatmentPBO 0.20322 22.484646 0.009038 0.9928
## TimeD1:TreatmentPBO -19.11111 6.694322 -2.854824 0.0048
## TimeD2:TreatmentPBO -27.83889 9.679338 -2.876115 0.0045
## TimeD3:TreatmentPBO -22.15556 9.266937 -2.390817 0.0179
## TimeD4:TreatmentPBO -27.17222 10.106000 -2.688722 0.0079
##
## Correlation:
## (Intr) TimeD1 TimeD2 TimeD3 TimeD4 GendrM TrtPBO GM:TPB
## TimeD1 -0.325
## TimeD2 -0.162 0.195
## TimeD3 -0.291 0.076 0.252
## TimeD4 -0.230 0.360 0.366 0.665
## GenderM -0.698 0.000 0.000 0.000 0.000
## TreatmentPBO -0.720 0.234 0.117 0.210 0.165 0.503
## GenderM:TreatmentPBO 0.486 0.000 0.000 0.000 0.000 -0.696 -0.660
## TimeD1:TreatmentPBO 0.223 -0.688 -0.134 -0.052 -0.248 0.000 -0.340 0.000
## TimeD2:TreatmentPBO 0.112 -0.134 -0.688 -0.174 -0.252 0.000 -0.170 0.000
## TimeD3:TreatmentPBO 0.200 -0.052 -0.174 -0.688 -0.457 0.000 -0.305 0.000
## TimeD4:TreatmentPBO 0.158 -0.248 -0.252 -0.457 -0.688 0.000 -0.240 0.000
## TD1:TP TD2:TP TD3:TP
## TimeD1
## TimeD2
## TimeD3
## TimeD4
## GenderM
## TreatmentPBO
## GenderM:TreatmentPBO
## TimeD1:TreatmentPBO
## TimeD2:TreatmentPBO 0.195
## TimeD3:TreatmentPBO 0.076 0.252
## TimeD4:TreatmentPBO 0.360 0.366 0.665
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.8562836 -0.6190432 -0.1082266 0.5224223 3.5519928
##
## Residual standard error: 43.03221
## Degrees of freedom: 190 total; 178 residual