Homework_1019_2
yuwenchen
10/19/2020
Exercise 1
The new drug Novafylline is developed to reduce the symptoms of intermittent claudication, which is a medical term for pain caused by too little blood flow, usually during exercise. Patients were randomly assigned to receive either Novafylline (active treatment) or a placebo in a double-blind study. A total of 38 patients (20 for treatment, 18 for placebo) underwent treadmill testing at baseline and at each of 4 monthly, follow-up visits. Patients were stratified by gender (17 females, 21 males). The primary measurement of efficacy is the walking distance on a treadmill until discontinuation due to claudication pain. Is there any distinction in exercise tolerance profiles between patients who receive Novafylline and those on placebo?
Source: Glenn Walker, G., & Shostak, J.(2010). Common Statistical Methods for Clinical Research with SAS Examples.
- Column 1: Treatment ID
- Column 2: Patient ID
- Column 3: Gender ID
- Column 4: Walking distance at baseline (month 0)
- Column 5: Walking distance at month 1
- Column 6: Walking distance at month 2
- Column 7: Walking distance at month 3
- Column 8: Walking distance at month 4
pacman::p_load("tidyverse","ggplot2","nlme","lme4", "afex","emmeans", "car")
dta_1<- read.csv("claudication.csv")
head(dta_1)## 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# wide form to long form
dtaL <- gather(data=dta_1, key=Time, value=Distance, D0:D4, factor_key=TRUE)
pd<-position_dodge(.1)
ggplot(dtaL, aes(Time, Distance, group=Treatment, col=Treatment))+
stat_summary(fun.y = mean, geom="line", alpha=.5,position = pd )+
stat_summary(fun.y= mean, geom="point", alpha=.5, position=pd)+
stat_summary(fun.data = mean_se, geom="errorbar",width=rel(.3), position = pd)+
geom_point(alpha=.3, position = pd)+
facet_grid(~Gender)+
theme_bw()+
theme(legend.position = c(.1, .8),
legend.background = element_rect(fill="white"))+
labs(x="Time(in month)",
y="Distance(in meters)")
# group mean
aggregate(cbind(D0,D1,D2,D3,D4)~Gender+Treatment, data=dta_1, 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.2222#group standard deviation
aggregate(cbind(D0,D1,D2,D3,D4)~Gender+Treatment, data=dta_1, 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# Female cases in ACT group repeated data correlation
cor(subset(dta_1, 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# Male cases in ACT group repeated data correlation
cor(subset(dta_1, 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#Female in PBO group repeated data correlation
cor(subset(dta_1, 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#Male in PBO group repeated data correlation
cor(subset(dta_1, 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.0000000m0_1<- gls(Distance ~ Time+ Time:Treatment+ Gender * Treatment , data=dtaL,
na.action=na.exclude,
weights = varIdent(form = ~1|PID),
correlation=corAR1(form = ~1|Time),
control=glsControl(opt='optim'))
summary(m0_1)## Generalized least squares fit by REML
## Model: Distance ~ Time + Time:Treatment + Gender * Treatment
## Data: dtaL
## AIC BIC logLik
## 1844.734 2007.005 -871.3668
##
## Correlation Structure: AR(1)
## Formula: ~1 | Time
## Parameter estimate(s):
## Phi
## 0.09611659
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | PID
## Parameter estimates:
## P101 P105 P109 P112 P117 P122 P123 P128
## 1.0000000 2.3100494 1.4860408 2.9762357 0.3706410 1.6900885 1.0376094 1.1800073
## P129 P132 P134 P136 P103 P107 P108 P113
## 0.6671421 0.7708146 1.8592185 2.6152716 0.8360770 4.9556513 0.8437202 0.8767682
## P115 P119 P124 P126 P104 P111 P114 P118
## 1.8701633 1.4966598 2.3960079 0.0214408 0.6225646 1.5492560 1.0861280 4.2638816
## P125 P127 P131 P133 P135 P102 P106 P110
## 0.6140227 1.8970965 0.5128498 1.2607998 3.3175320 1.1837631 1.2849888 2.3898384
## P116 P120 P121 P130 P137 P138
## 1.7955836 3.3151959 0.8419682 2.4160621 0.4500130 1.6842631
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 174.90612 0.463409 377.4337 0.0000
## TimeD1 22.02014 0.654884 33.6245 0.0000
## TimeD2 19.99438 0.654884 30.5312 0.0000
## TimeD3 7.08639 0.654884 10.8208 0.0000
## TimeD4 17.98459 0.654884 27.4623 0.0000
## GenderM 8.79758 2.626616 3.3494 0.0010
## TreatmentPBO -12.41750 5.717959 -2.1717 0.0312
## TimeD1:TreatmentPBO -36.49755 7.210934 -5.0614 0.0000
## TimeD2:TreatmentPBO -30.77176 7.210934 -4.2674 0.0000
## TimeD3:TreatmentPBO -0.42953 7.210934 -0.0596 0.9526
## TimeD4:TreatmentPBO -14.51739 7.210934 -2.0132 0.0456
## TreatmentPBO:GenderM 16.10332 5.286480 3.0461 0.0027
##
## Correlation:
## (Intr) TimeD1 TimeD2 TimeD3 TimeD4 GendrM TrtPBO TD1:TP
## TimeD1 -0.707
## TimeD2 -0.707 0.500
## TimeD3 -0.707 0.500 0.500
## TimeD4 -0.707 0.500 0.500 0.500
## GenderM -0.035 0.000 0.000 0.000 0.000
## TreatmentPBO -0.055 0.035 0.035 0.035 0.035 0.003
## TimeD1:TreatmentPBO 0.039 -0.055 -0.028 -0.028 -0.028 0.000 -0.631
## TimeD2:TreatmentPBO 0.039 -0.028 -0.055 -0.028 -0.028 0.000 -0.631 0.500
## TimeD3:TreatmentPBO 0.039 -0.028 -0.028 -0.055 -0.028 0.000 -0.631 0.500
## TimeD4:TreatmentPBO 0.039 -0.028 -0.028 -0.028 -0.055 0.000 -0.631 0.500
## TreatmentPBO:GenderM 0.030 0.000 0.000 0.000 0.000 -0.499 -0.394 0.000
## TD2:TP TD3:TP TD4:TP
## TimeD1
## TimeD2
## TimeD3
## TimeD4
## GenderM
## TreatmentPBO
## TimeD1:TreatmentPBO
## TimeD2:TreatmentPBO
## TimeD3:TreatmentPBO 0.500
## TimeD4:TreatmentPBO 0.500 0.500
## TreatmentPBO:GenderM 0.000 0.000 0.000
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.869631964 -0.913090362 -0.009671972 0.739278035 2.118939040
##
## Residual standard error: 21.85608
## Degrees of freedom: 190 total; 178 residualExercise 2
Patients were randomized to receive one of two daily doses (low or high) of a new treatment for Alzheimer’s disease or a placebo in a study. Each patient returned to the clinic every 2 months for 1 year for assessment of disease progression based on cognitive measurements on the Alzheimer’s Disease Assessment Scale (ADAS). This test evaluates memory, language, and praxis function, and is based on the sum of scores from an 11-item scale, with a potential range of 0 to 70, higher scores indicative of greater disease severity. The primary goal is to determine if the rate of disease progression is slowed with active treatment compared with a placebo. Is there a difference in response profiles over time among the three groups?
Source: Glenn Walker, G., & Shostak, J.(2010). Common Statistical Methods for Clinical Research with SAS Examples.
- Column 1: Treatment ID (L = Low dose, H = High dose, P = Placebo)
- Column 2: Patient ID
- Column 3: ADAS score at month 2
- Column 4: ADAS score at month 4
- Column 5: ADAS score at month 6
- Column 6: ADAS score at month 8
- Column 7: ADAS score at month 10
- Column 8: ADAS score at month 12
dta_2<-read.table("adas.txt", h=T,na.strings='.')
head(dta_2)## 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 57dta_2 <- dta_2%>% mutate(Baseline= dta_2[,3])
colnames(dta_2)<-c("Treatment", "PID", "2","4","6","8","10","12","Baseline")
dta_2L<-gather(dta_2, key=Time, value=Score, "2":"12")
dta_2L$Time_f <- as.factor(as.numeric(dta_2L$Time) - 2)
dta_2L$Score<-as.numeric(dta_2L$Score)
tail(dta_2L)## Treatment PID Baseline Time Score Time_f
## 475 P 64 21 12 27 10
## 476 P 65 26 12 32 10
## 477 P 70 53 12 NA 10
## 478 P 71 32 12 35 10
## 479 P 74 NA 12 54 10
## 480 P 77 24 12 30 10pd<-position_dodge(.2)
p <- ggplot(dta_2L, aes(Time_f, Score,
group=Treatment, shape=Treatment, linetype=Treatment)) +
stat_summary(fun=mean, geom="line", position = pd) +
stat_summary(fun=mean, geom="point", position = pd) +
stat_summary(fun.data=mean_se, geom="errorbar", width=rel(0.3), position = pd) +
scale_shape_manual(values=c(1, 2,5)) +
scale_y_continuous(breaks=seq(20, 55,5))+
labs(x="Time (months since baseline)",
y="Mean ADAS Score",
linetype="Treatment", shape="Treatment") +
theme_minimal()+
theme(legend.position=c(.2, .9))
suppressWarnings(suppressMessages(ggplot2:::print.ggplot(p)))
dta_2w <- read.table('adas.txt', h=T, na.strings='.')
dta_2w <- dta_2w%>% mutate(Baseline= dta_2[,3])
dta_2w %>% dplyr::group_by(Treatment) %>% dplyr::select(starts_with("adas")) %>%
furniture::table1(digits=2, total=FALSE, test=FALSE, output="html")| H | L | P | |
|---|---|---|---|
| n = 18 | n = 17 | n = 21 | |
| adas02 | |||
| 30.22 (8.56) | 32.82 (7.51) | 29.62 (6.78) | |
| adas04 | |||
| 32.44 (7.33) | 35.12 (8.64) | 33.43 (6.42) | |
| adas06 | |||
| 33.56 (8.56) | 37.65 (9.09) | 34.86 (6.43) | |
| adas08 | |||
| 33.33 (9.06) | 36.41 (11.09) | 36.38 (6.34) | |
| adas10 | |||
| 33.72 (7.44) | 38.35 (11.19) | 37.57 (7.03) | |
| adas12 | |||
| 34.28 (7.10) | 39.18 (10.77) | 39.05 (8.80) |
knitr::kable(cor(dta_2w[,2:7], use="pair"))| PID | adas02 | adas04 | adas06 | adas08 | adas10 | |
|---|---|---|---|---|---|---|
| PID | 1.0000000 | 0.1548313 | 0.1591352 | 0.1329467 | 0.1810840 | 0.0416975 |
| adas02 | 0.1548313 | 1.0000000 | 0.8954000 | 0.7544220 | 0.7740050 | 0.7006653 |
| adas04 | 0.1591352 | 0.8954000 | 1.0000000 | 0.8144181 | 0.7936916 | 0.7006978 |
| adas06 | 0.1329467 | 0.7544220 | 0.8144181 | 1.0000000 | 0.8628966 | 0.6716075 |
| adas08 | 0.1810840 | 0.7740050 | 0.7936916 | 0.8628966 | 1.0000000 | 0.8472388 |
| adas10 | 0.0416975 | 0.7006653 | 0.7006978 | 0.6716075 | 0.8472388 | 1.0000000 |
scatterplotMatrix(~ adas02 + adas04 + adas06 + adas08 + adas10 + adas12,
data=dta_2w, col=8,
smooth=FALSE, regLine=FALSE,
ellipse=list(levels=c(.975),
fill=TRUE,
fill.alpha=.1))
Repeated measures ANOVA
m_aov <- afex::aov_car(Score ~ Treatment*Time_f + Error(PID/Time_f), data=dta_2L)
knitr::kable(nice(m_aov))| Effect | df | MSE | F | ges | p.value |
|---|---|---|---|---|---|
| Treatment | 2, 53 | 328.01 | 1.10 | .032 | .342 |
| Time_f | 3.27, 173.16 | 25.22 | 18.57 *** | .066 | <.001 |
| Treatment:Time_f | 6.53, 173.16 | 25.22 | 1.35 | .010 | .234 |
emmeans(m_aov, ~ Time_f | Treatment)## Treatment = H:
## Time_f emmean SE df lower.CL upper.CL
## X0 30.2 1.93 82.6 26.4 34.1
## X2 32.4 1.93 82.6 28.6 36.3
## X4 33.5 1.93 82.6 29.7 37.4
## X6 33.3 1.93 82.6 29.5 37.2
## X8 33.7 1.93 82.6 29.9 37.6
## X10 34.3 1.93 82.6 30.4 38.1
##
## Treatment = L:
## Time_f emmean SE df lower.CL upper.CL
## X0 32.8 1.96 83.8 28.9 36.7
## X2 35.1 1.96 83.8 31.2 39.0
## X4 37.6 1.96 83.8 33.7 41.5
## X6 36.4 1.96 83.8 32.5 40.3
## X8 38.3 1.96 83.8 34.4 42.2
## X10 39.2 1.96 83.8 35.3 43.1
##
## Treatment = P:
## Time_f emmean SE df lower.CL upper.CL
## X0 29.6 1.87 79.5 25.9 33.3
## X2 33.4 1.87 79.5 29.7 37.1
## X4 34.8 1.87 79.5 31.1 38.6
## X6 36.4 1.87 79.5 32.7 40.1
## X8 37.6 1.87 79.5 33.8 41.3
## X10 39.0 1.87 79.5 35.3 42.7
##
## Warning: EMMs are biased unless design is perfectly balanced
## Confidence level used: 0.95lattice::xyplot(Score ~ Baseline| Time_f, groups=Treatment, data=dta_2L,
type=c("p","r","g"), layout=c(6,1),
xlab="Baseline ADAS score",
ylab="ADAS score",
auto.key = T)
Comparing adjusted means: analysis of covariance
m_ancova <- gls(Score ~ Baseline + Time_f, data=dta_2L, na.action=na.exclude, correlation=corSymm(form= ~ 1 | Treatment), weights=varIdent(form= ~ 1 | Time_f), control=glsControl(opt=‘optim’))