Claudication Example
Szu-Yu Chen
2020-10-26
1 Data Management
Load data
pacman::p_load(tidyr,tidyverse, nlme, car, afex, emmeans, lattice)
# Input data
dta <- read.csv("C:/Users/ASUS/Desktop/data/claudication.csv", header = T)
str(dta)## '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 ...head(dta)## 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 205Reshape Data (wide to long)
dta <- mutate(dta, Treatment = factor(Treatment),
Gender = factor(Gender),
PID = factor(PID))
dtaL <- gather(dta, Month, Distance, D0:D4, factor_key=TRUE)
head(dtaL)## 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 140ggplot(dtaL, aes(Month, Distance,
group = Treatment,
shape = Treatment,
linetype = Treatment,
color = Treatment)) +
geom_point() +
stat_summary(fun = mean, geom = "line") +
stat_summary(fun = mean, geom = "point") +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.3) +
scale_shape_manual(values = c(1, 2)) +
facet_wrap( ~ Gender)+
labs(x = "TIME (in months)",
y = "Walking distance (in meters)",
linetype = "Treatment", shape = "Treatment") +
theme(legend.justification = c(1, 1),
legend.position = c(1, 1),
legend.background = element_rect(fill = "white",color = "black"))
2 Analysis & Output
2.1 Mean, SD, Correlation matrix
Mean & SD
dta0 <- dta %>% select(-2)
# Mean
aggregate(.~ Treatment + Gender, data = dta0, 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# SD
aggregate(.~ Treatment + Gender, data = dta0, sd, na.rm=T)## 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.35182Female & Male Correlation matrix (Treatment)
cor(subset(dta, 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(dta, 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.0000000Female & Male Correlation matrix (Placebo)
cor(subset(dta, 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(dta, 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.00000002.2 GLS
2.2.1 m1
M1: (Distance ~ Month + Gender*Treatment + Month:Treatment), Covariance Matrix: CS
summary(m1 <- gls(Distance ~ Month + Gender*Treatment + Month:Treatment,
weights = varIdent(form = ~ 1 | Month),
correlation = corSymm(form = ~ 1 | PID),
data = dtaL))## Generalized least squares fit by REML
## Model: Distance ~ Month + Gender * Treatment + Month: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 | Month
## 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
## MonthD1 15.00000 4.607348 3.255669 0.0014
## MonthD2 30.45000 6.661777 4.570852 0.0000
## MonthD3 31.10000 6.377944 4.876180 0.0000
## MonthD4 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
## MonthD1:TreatmentPBO -19.11111 6.694322 -2.854824 0.0048
## MonthD2:TreatmentPBO -27.83889 9.679338 -2.876115 0.0045
## MonthD3:TreatmentPBO -22.15556 9.266937 -2.390817 0.0179
## MonthD4:TreatmentPBO -27.17222 10.106000 -2.688722 0.0079
##
## Correlation:
## (Intr) MnthD1 MnthD2 MnthD3 MnthD4 GendrM TrtPBO GM:TPB
## MonthD1 -0.325
## MonthD2 -0.162 0.195
## MonthD3 -0.291 0.076 0.252
## MonthD4 -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
## MonthD1:TreatmentPBO 0.223 -0.688 -0.134 -0.052 -0.248 0.000 -0.340 0.000
## MonthD2:TreatmentPBO 0.112 -0.134 -0.688 -0.174 -0.252 0.000 -0.170 0.000
## MonthD3:TreatmentPBO 0.200 -0.052 -0.174 -0.688 -0.457 0.000 -0.305 0.000
## MonthD4:TreatmentPBO 0.158 -0.248 -0.252 -0.457 -0.688 0.000 -0.240 0.000
## MD1:TP MD2:TP MD3:TP
## MonthD1
## MonthD2
## MonthD3
## MonthD4
## GenderM
## TreatmentPBO
## GenderM:TreatmentPBO
## MonthD1:TreatmentPBO
## MonthD2:TreatmentPBO 0.195
## MonthD3:TreatmentPBO 0.076 0.252
## MonthD4: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