1 data input

# input data
dta <- read.csv("claudication.csv", h = T)

# see the data structure
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 ...
# show the first six data
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 205

2 Wide to long format

library(dplyr)
## 
## 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 ...
head(dtaL)
##   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

3 line chart by gender

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"))

4 Mean & sd

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
#sd
aggregate(.~ Treatment + Gender, data = dta1_nid, 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.35182

5 Corr

#ACT and Female
cor(subset(dta1, 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.0000000
#ACT and Male
cor(subset(dta1, 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.0000000
#PBO and Female
cor(subset(dta1, 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.0000000
#PBO and Male
cor(subset(dta1, 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

#gls

library(nlme)
## 
## 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