### Repeated measures ANOVA (doubly multivariate mixed between-within), multivariate approach
### 1. Read the data
setwd("H:/data/User/Classes/multiv_13/YTs/YTSplit/Week 03 wmv - 2016/Wk03.R")
getwd()
## [1] "H:/data/User/Classes/multiv_13/YTs/YTSplit/Week 03 wmv - 2016/Wk03.R"
Wk03a<- read.csv("Class03_homework_2018a.csv")
### 2. Complex mixed between-within with multivariate effect estimation - omnibus
library("car")
## Loading required package: carData
#I define the "data" for the within-subjects design as follows:
# This syntax is in lieu of reshape/restructure that we did in CLP 6528 to
# make the data set long
# It is the R equivalent of the "Define" repeated measures window before
# We actually set up the analysis.
domain <- factor(rep(c("cogment", "phys", "pain"), each=3),
levels=c("cogment", "phys", "pain"))
task <- ordered(rep(1:3, 3))
idata <- data.frame(domain,task)
idata
## domain task
## 1 cogment 1
## 2 cogment 2
## 3 cogment 3
## 4 phys 1
## 5 phys 2
## 6 phys 3
## 7 pain 1
## 8 pain 2
## 9 pain 3
Wk03a$Condition<-factor(Wk03a$Condition)
options(contrasts = c("contr.helmert", "contr.poly"))
mod.1 <- lm(cbind(MOCA,BDI,STAI,WalkSpeed,CalfStrength,V02,PainDis,PainCat,PainIntens) ~ Condition,
data=Wk03a)
(av.1 <- Anova(mod.1, idata=idata, idesign=~domain*task, type=3))
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.99535 42341 1 198 <2e-16 ***
## Condition 1 0.67564 412 1 198 <2e-16 ***
## domain 1 0.02173 2 2 197 0.1149
## Condition:domain 1 0.00856 1 2 197 0.4286
## task 1 0.00171 0 2 197 0.8448
## Condition:task 1 0.01906 2 2 197 0.1502
## domain:task 1 0.02502 1 4 195 0.2908
## Condition:domain:task 1 0.72301 127 4 195 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Like SPSS, when multivariate is true, we get both multivariate and univariate
# Repeated Measures tests
summary(av.1, multivariate=TRUE)
## Warning in summary.Anova.mlm(av.1, multivariate = TRUE): HF eps > 1 treated as 1
##
## Type III Repeated Measures MANOVA Tests:
##
## ------------------------------------------
##
## Term: (Intercept)
##
## Response transformation matrix:
## (Intercept)
## MOCA 1
## BDI 1
## STAI 1
## WalkSpeed 1
## CalfStrength 1
## V02 1
## PainDis 1
## PainCat 1
## PainIntens 1
##
## Sum of squares and products for the hypothesis:
## (Intercept)
## (Intercept) 40998264
##
## Multivariate Tests: (Intercept)
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.99535 42341.4 1 198 < 2.22e-16 ***
## Wilks 1 0.00465 42341.4 1 198 < 2.22e-16 ***
## Hotelling-Lawley 1 213.84548 42341.4 1 198 < 2.22e-16 ***
## Roy 1 213.84548 42341.4 1 198 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: Condition
##
## Response transformation matrix:
## (Intercept)
## MOCA 1
## BDI 1
## STAI 1
## WalkSpeed 1
## CalfStrength 1
## V02 1
## PainDis 1
## PainCat 1
## PainIntens 1
##
## Sum of squares and products for the hypothesis:
## (Intercept)
## (Intercept) 399344
##
## Multivariate Tests: Condition
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.6756368 412.4269 1 198 < 2.22e-16 ***
## Wilks 1 0.3243632 412.4269 1 198 < 2.22e-16 ***
## Hotelling-Lawley 1 2.0829641 412.4269 1 198 < 2.22e-16 ***
## Roy 1 2.0829641 412.4269 1 198 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: domain
##
## Response transformation matrix:
## domain1 domain2
## MOCA 1 0
## BDI 1 0
## STAI 1 0
## WalkSpeed 0 1
## CalfStrength 0 1
## V02 0 1
## PainDis -1 -1
## PainCat -1 -1
## PainIntens -1 -1
##
## Sum of squares and products for the hypothesis:
## domain1 domain2
## domain1 229.3897 775.7309
## domain2 775.7309 2623.3025
##
## Multivariate Tests: domain
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.0217284 2.187789 2 197 0.11488
## Wilks 1 0.9782716 2.187789 2 197 0.11488
## Hotelling-Lawley 1 0.0222111 2.187789 2 197 0.11488
## Roy 1 0.0222111 2.187789 2 197 0.11488
##
## ------------------------------------------
##
## Term: Condition:domain
##
## Response transformation matrix:
## domain1 domain2
## MOCA 1 0
## BDI 1 0
## STAI 1 0
## WalkSpeed 0 1
## CalfStrength 0 1
## V02 0 1
## PainDis -1 -1
## PainCat -1 -1
## PainIntens -1 -1
##
## Sum of squares and products for the hypothesis:
## domain1 domain2
## domain1 35.81157 -151.8094
## domain2 -151.80937 643.5374
##
## Multivariate Tests: Condition:domain
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.0085633 0.8507679 2 197 0.42865
## Wilks 1 0.9914367 0.8507679 2 197 0.42865
## Hotelling-Lawley 1 0.0086372 0.8507679 2 197 0.42865
## Roy 1 0.0086372 0.8507679 2 197 0.42865
##
## ------------------------------------------
##
## Term: task
##
## Response transformation matrix:
## task.L task.Q
## MOCA -7.071068e-01 0.4082483
## BDI -7.850462e-17 -0.8164966
## STAI 7.071068e-01 0.4082483
## WalkSpeed -7.071068e-01 0.4082483
## CalfStrength -7.850462e-17 -0.8164966
## V02 7.071068e-01 0.4082483
## PainDis -7.071068e-01 0.4082483
## PainCat -7.850462e-17 -0.8164966
## PainIntens 7.071068e-01 0.4082483
##
## Sum of squares and products for the hypothesis:
## task.L task.Q
## task.L 53.14105 51.43377
## task.Q 51.43377 49.78135
##
## Multivariate Tests: task
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.0017106 0.168785 2 197 0.84481
## Wilks 1 0.9982894 0.168785 2 197 0.84481
## Hotelling-Lawley 1 0.0017136 0.168785 2 197 0.84481
## Roy 1 0.0017136 0.168785 2 197 0.84481
##
## ------------------------------------------
##
## Term: Condition:task
##
## Response transformation matrix:
## task.L task.Q
## MOCA -7.071068e-01 0.4082483
## BDI -7.850462e-17 -0.8164966
## STAI 7.071068e-01 0.4082483
## WalkSpeed -7.071068e-01 0.4082483
## CalfStrength -7.850462e-17 -0.8164966
## V02 7.071068e-01 0.4082483
## PainDis -7.071068e-01 0.4082483
## PainCat -7.850462e-17 -0.8164966
## PainIntens 7.071068e-01 0.4082483
##
## Sum of squares and products for the hypothesis:
## task.L task.Q
## task.L 0.7613683 -27.70725
## task.Q -27.7072508 1008.30539
##
## Multivariate Tests: Condition:task
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.0190639 1.914287 2 197 0.15018
## Wilks 1 0.9809361 1.914287 2 197 0.15018
## Hotelling-Lawley 1 0.0194344 1.914287 2 197 0.15018
## Roy 1 0.0194344 1.914287 2 197 0.15018
##
## ------------------------------------------
##
## Term: domain:task
##
## Response transformation matrix:
## domain1:task.L domain2:task.L domain1:task.Q domain2:task.Q
## MOCA -7.071068e-01 0.000000e+00 0.4082483 0.0000000
## BDI -7.850462e-17 0.000000e+00 -0.8164966 0.0000000
## STAI 7.071068e-01 0.000000e+00 0.4082483 0.0000000
## WalkSpeed 0.000000e+00 -7.071068e-01 0.0000000 0.4082483
## CalfStrength 0.000000e+00 -7.850462e-17 0.0000000 -0.8164966
## V02 0.000000e+00 7.071068e-01 0.0000000 0.4082483
## PainDis 7.071068e-01 7.071068e-01 -0.4082483 -0.4082483
## PainCat 7.850462e-17 7.850462e-17 0.8164966 0.8164966
## PainIntens -7.071068e-01 -7.071068e-01 -0.4082483 -0.4082483
##
## Sum of squares and products for the hypothesis:
## domain1:task.L domain2:task.L domain1:task.Q domain2:task.Q
## domain1:task.L 11.80999 60.80139 -49.68966 -88.41242
## domain2:task.L 60.80139 313.02400 -255.81742 -455.17396
## domain1:task.Q -49.68966 -255.81742 209.06561 371.98883
## domain2:task.Q -88.41242 -455.17396 371.98883 661.87683
##
## Multivariate Tests: domain:task
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.0250180 1.250921 4 195 0.29083
## Wilks 1 0.9749820 1.250921 4 195 0.29083
## Hotelling-Lawley 1 0.0256599 1.250921 4 195 0.29083
## Roy 1 0.0256599 1.250921 4 195 0.29083
##
## ------------------------------------------
##
## Term: Condition:domain:task
##
## Response transformation matrix:
## domain1:task.L domain2:task.L domain1:task.Q domain2:task.Q
## MOCA -7.071068e-01 0.000000e+00 0.4082483 0.0000000
## BDI -7.850462e-17 0.000000e+00 -0.8164966 0.0000000
## STAI 7.071068e-01 0.000000e+00 0.4082483 0.0000000
## WalkSpeed 0.000000e+00 -7.071068e-01 0.0000000 0.4082483
## CalfStrength 0.000000e+00 -7.850462e-17 0.0000000 -0.8164966
## V02 0.000000e+00 7.071068e-01 0.0000000 0.4082483
## PainDis 7.071068e-01 7.071068e-01 -0.4082483 -0.4082483
## PainCat 7.850462e-17 7.850462e-17 0.8164966 0.8164966
## PainIntens -7.071068e-01 -7.071068e-01 -0.4082483 -0.4082483
##
## Sum of squares and products for the hypothesis:
## domain1:task.L domain2:task.L domain1:task.Q domain2:task.Q
## domain1:task.L 171.1624 -1237.890 3447.108 2806.134
## domain2:task.L -1237.8901 8952.737 -24930.366 -20294.672
## domain1:task.Q 3447.1085 -24930.366 69422.701 56513.850
## domain2:task.Q 2806.1336 -20294.672 56513.850 46005.344
##
## Multivariate Tests: Condition:domain:task
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.7230087 127.2483 4 195 < 2.22e-16 ***
## Wilks 1 0.2769913 127.2483 4 195 < 2.22e-16 ***
## Hotelling-Lawley 1 2.6102214 127.2483 4 195 < 2.22e-16 ***
## Roy 1 2.6102214 127.2483 4 195 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 4555363 1 21302 198 42341.4045 <2e-16 ***
## Condition 44372 1 21302 198 412.4269 <2e-16 ***
## domain 462 2 42337 396 2.1585 0.1169
## Condition:domain 185 2 42337 396 0.8638 0.4223
## task 34 2 39109 396 0.1737 0.8406
## Condition:task 336 2 39109 396 1.7029 0.1835
## domain:task 509 4 79831 792 1.2616 0.2836
## Condition:domain:task 46184 4 79831 792 114.5470 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## domain 0.99899 0.90516
## Condition:domain 0.99899 0.90516
## task 0.98622 0.25494
## Condition:task 0.98622 0.25494
## domain:task 0.96307 0.59653
## Condition:domain:task 0.96307 0.59653
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## domain 0.99899 0.1169
## Condition:domain 0.99899 0.4223
## task 0.98641 0.8378
## Condition:task 0.98641 0.1839
## domain:task 0.98149 0.2839
## Condition:domain:task 0.98149 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## domain 1.0091681 1.168501e-01
## Condition:domain 1.0091681 4.223473e-01
## task 0.9962656 8.398436e-01
## Condition:task 0.9962656 1.836109e-01
## domain:task 1.0038254 2.835536e-01
## Condition:domain:task 1.0038254 4.546606e-77
### 3. Complex mixed between-within with multivariate effect estimation - univariates
#Given the three-way, we will do followup Condition by Task ANOVAS
#For each Domain
#cogment -- interaction is significant and will require followup
task_cogment <- factor(rep(c("MOCA","BDI","STAI"), each=1),
levels=c("MOCA","BDI","STAI"))
idata_cogment <- data.frame(task_cogment)
idata_cogment
## task_cogment
## 1 MOCA
## 2 BDI
## 3 STAI
# we have to specify this format of internal contrast (for how R calculates
# the sums of squares) here, as well as type 3 SS below, to get results
# to match SPSS and SAS. I repeat this below for each univariate
options(contrasts = c("contr.helmert", "contr.poly"))
mod.cogment <- lm(cbind(MOCA,BDI,STAI) ~ Condition,
data=Wk03a)
(av.cogment <- Anova(mod.cogment, idata=idata_cogment, idesign=~task_cogment, type=3))
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.98603 13971.0 1 198 <2e-16 ***
## Condition 1 0.43693 153.6 1 198 <2e-16 ***
## task_cogment 1 0.00015 0.0 2 197 0.9857
## Condition:task_cogment 1 0.48816 93.9 2 197 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(av.cogment, multivariate=TRUE)
##
## Type III Repeated Measures MANOVA Tests:
##
## ------------------------------------------
##
## Term: (Intercept)
##
## Response transformation matrix:
## (Intercept)
## MOCA 1
## BDI 1
## STAI 1
##
## Sum of squares and products for the hypothesis:
## (Intercept)
## (Intercept) 4525635
##
## Multivariate Tests: (Intercept)
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.98603 13971.02 1 198 < 2.22e-16 ***
## Wilks 1 0.01397 13971.02 1 198 < 2.22e-16 ***
## Hotelling-Lawley 1 70.56072 13971.02 1 198 < 2.22e-16 ***
## Roy 1 70.56072 13971.02 1 198 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: Condition
##
## Response transformation matrix:
## (Intercept)
## MOCA 1
## BDI 1
## STAI 1
##
## Sum of squares and products for the hypothesis:
## (Intercept)
## (Intercept) 49769.64
##
## Multivariate Tests: Condition
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.4369292 153.6431 1 198 < 2.22e-16 ***
## Wilks 1 0.5630708 153.6431 1 198 < 2.22e-16 ***
## Hotelling-Lawley 1 0.7759755 153.6431 1 198 < 2.22e-16 ***
## Roy 1 0.7759755 153.6431 1 198 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: task_cogment
##
## Response transformation matrix:
## task_cogment1 task_cogment2
## MOCA 1 0
## BDI 0 1
## STAI -1 -1
##
## Sum of squares and products for the hypothesis:
## task_cogment1 task_cogment2
## task_cogment1 2.768401 -1.240718
## task_cogment2 -1.240718 0.556054
##
## Multivariate Tests: task_cogment
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.0001457 0.01435734 2 197 0.98575
## Wilks 1 0.9998543 0.01435734 2 197 0.98575
## Hotelling-Lawley 1 0.0001458 0.01435734 2 197 0.98575
## Roy 1 0.0001458 0.01435734 2 197 0.98575
##
## ------------------------------------------
##
## Term: Condition:task_cogment
##
## Response transformation matrix:
## task_cogment1 task_cogment2
## MOCA 1 0
## BDI 0 1
## STAI -1 -1
##
## Sum of squares and products for the hypothesis:
## task_cogment1 task_cogment2
## task_cogment1 3195.316 9541.467
## task_cogment2 9541.467 28491.574
##
## Multivariate Tests: Condition:task_cogment
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.4881629 93.94406 2 197 < 2.22e-16 ***
## Wilks 1 0.5118371 93.94406 2 197 < 2.22e-16 ***
## Hotelling-Lawley 1 0.9537468 93.94406 2 197 < 2.22e-16 ***
## Roy 1 0.9537468 93.94406 2 197 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 1508545 1 21379 198 13971.0220 <2e-16 ***
## Condition 16590 1 21379 198 153.6431 <2e-16 ***
## task_cogment 3 2 38078 396 0.0158 0.9843
## Condition:task_cogment 14764 2 38078 396 76.7678 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## task_cogment 0.96472 0.02909
## Condition:task_cogment 0.96472 0.02909
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## task_cogment 0.96593 0.9823
## Condition:task_cogment 0.96593 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## task_cogment 0.9752702 9.828453e-01
## Condition:task_cogment 0.9752702 2.963558e-28
#phys -- interaction is significant and will require followup
task_phys <- factor(rep(c("WalkSpeed", "CalfStrength","V02"), each=1),
levels=c("WalkSpeed", "CalfStrength","V02"))
idata_phys <- data.frame(task_phys)
idata_phys
## task_phys
## 1 WalkSpeed
## 2 CalfStrength
## 3 V02
# we have to specify this format of internal contrast (for how R calculates
# the sums of squares) here, as well as type 3 SS below, to get results
# to match SPSS and SAS. I repeat this below for each univariate
options(contrasts = c("contr.helmert", "contr.poly"))
mod.phys <- lm(cbind(WalkSpeed, CalfStrength,V02) ~ Condition,
data=Wk03a)
(av.phys <- Anova(mod.phys, idata=idata_phys, idesign=~task_phys, type=3))
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.98617 14119.2 1 198 < 2.2e-16 ***
## Condition 1 0.35902 110.9 1 198 < 2.2e-16 ***
## task_phys 1 0.01267 1.3 2 197 0.2847
## Condition:task_phys 1 0.29328 40.9 2 197 1.416e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(av.phys, multivariate=TRUE)
## Warning in summary.Anova.mlm(av.phys, multivariate = TRUE): HF eps > 1 treated
## as 1
##
## Type III Repeated Measures MANOVA Tests:
##
## ------------------------------------------
##
## Term: (Intercept)
##
## Response transformation matrix:
## (Intercept)
## WalkSpeed 1
## CalfStrength 1
## V02 1
##
## Sum of squares and products for the hypothesis:
## (Intercept)
## (Intercept) 4680414
##
## Multivariate Tests: (Intercept)
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.98617 14119.22 1 198 < 2.22e-16 ***
## Wilks 1 0.01383 14119.22 1 198 < 2.22e-16 ***
## Hotelling-Lawley 1 71.30917 14119.22 1 198 < 2.22e-16 ***
## Roy 1 71.30917 14119.22 1 198 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: Condition
##
## Response transformation matrix:
## (Intercept)
## WalkSpeed 1
## CalfStrength 1
## V02 1
##
## Sum of squares and products for the hypothesis:
## (Intercept)
## (Intercept) 36763.76
##
## Multivariate Tests: Condition
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.3590237 110.9038 1 198 < 2.22e-16 ***
## Wilks 1 0.6409763 110.9038 1 198 < 2.22e-16 ***
## Hotelling-Lawley 1 0.5601200 110.9038 1 198 < 2.22e-16 ***
## Roy 1 0.5601200 110.9038 1 198 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: task_phys
##
## Response transformation matrix:
## task_phys1 task_phys2
## WalkSpeed 1 0
## CalfStrength 0 1
## V02 -1 -1
##
## Sum of squares and products for the hypothesis:
## task_phys1 task_phys2
## task_phys1 342.14129 -55.012367
## task_phys2 -55.01237 8.845353
##
## Multivariate Tests: task_phys
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.0126736 1.264378 2 197 0.2847
## Wilks 1 0.9873264 1.264378 2 197 0.2847
## Hotelling-Lawley 1 0.0128363 1.264378 2 197 0.2847
## Roy 1 0.0128363 1.264378 2 197 0.2847
##
## ------------------------------------------
##
## Term: Condition:task_phys
##
## Response transformation matrix:
## task_phys1 task_phys2
## WalkSpeed 1 0
## CalfStrength 0 1
## V02 -1 -1
##
## Sum of squares and products for the hypothesis:
## task_phys1 task_phys2
## task_phys1 9175.001 -3125.838
## task_phys2 -3125.838 1064.944
##
## Multivariate Tests: Condition:task_phys
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.2932765 40.87558 2 197 1.4161e-15 ***
## Wilks 1 0.7067235 40.87558 2 197 1.4161e-15 ***
## Hotelling-Lawley 1 0.4149805 40.87558 2 197 1.4161e-15 ***
## Roy 1 0.4149805 40.87558 2 197 1.4161e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 1560138 1 21879 198 14119.2158 <2e-16 ***
## Condition 12255 1 21879 198 110.9038 <2e-16 ***
## task_phys 271 2 42110 396 1.2727 0.2812
## Condition:task_phys 8911 2 42110 396 41.8972 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## task_phys 0.99405 0.5558
## Condition:task_phys 0.99405 0.5558
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## task_phys 0.99409 0.2811
## Condition:task_phys 0.99409 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## task_phys 1.004143 2.812257e-01
## Condition:task_phys 1.004143 3.123973e-17
#pain -- interaction is significant and will require followup
head(Wk03a)
## ï..const ID Condition MOCA BDI STAI WalkSpeed CalfStrength
## 1 1 1 -1 42.29156 36.06104 60.28816 44.05973 51.29967
## 2 1 2 -1 29.04060 60.75548 31.91447 46.59351 41.14846
## 3 1 3 -1 57.59534 50.02032 37.15589 35.49083 47.94157
## 4 1 4 -1 33.71176 52.26377 28.57453 22.67779 55.12031
## 5 1 5 -1 54.47724 41.30908 50.76746 35.36213 42.69576
## 6 1 6 -1 48.87171 46.27396 36.27426 51.92100 46.44171
## V02 PainDis PainCat PainIntens
## 1 45.54295 56.23380 36.08248 40.76668
## 2 39.85025 52.15265 33.64651 60.19698
## 3 33.87827 47.50576 39.05012 53.52304
## 4 43.51130 63.64303 21.64790 47.95669
## 5 67.13485 51.19453 48.65334 27.84717
## 6 62.29283 57.53720 24.14119 47.63323
task_pain <- factor(rep(c("PainDis", "PainCat", "PainIntens"), each=1),
levels=c("PainDis", "PainCat", "PainIntens"))
idata_pain <- data.frame(task_pain)
idata_pain
## task_pain
## 1 PainDis
## 2 PainCat
## 3 PainIntens
# we have to specify this format of internal contrast (for how R calculates
# the sums of squares) here, as well as type 3 SS below, to get results
# to match SPSS and SAS. I repeat this below for each univariate
options(contrasts = c("contr.helmert", "contr.poly"))
mod.pain <- lm(cbind(PainDis, PainCat, PainIntens) ~ Condition,
data=Wk03a)
(av.pain <- Anova(mod.pain, idata=idata_pain, idesign=~task_pain, type=3))
##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 0.98648 14447.3 1 198 <2e-16 ***
## Condition 1 0.43531 152.6 1 198 <2e-16 ***
## task_pain 1 0.01685 1.7 2 197 0.1875
## Condition:task_pain 1 0.58959 141.5 2 197 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(av.pain, multivariate=TRUE)
##
## Type III Repeated Measures MANOVA Tests:
##
## ------------------------------------------
##
## Term: (Intercept)
##
## Response transformation matrix:
## (Intercept)
## PainDis 1
## PainCat 1
## PainIntens 1
##
## Sum of squares and products for the hypothesis:
## (Intercept)
## (Intercept) 4461424
##
## Multivariate Tests: (Intercept)
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.98648 14447.31 1 198 < 2.22e-16 ***
## Wilks 1 0.01352 14447.31 1 198 < 2.22e-16 ***
## Hotelling-Lawley 1 72.96620 14447.31 1 198 < 2.22e-16 ***
## Roy 1 72.96620 14447.31 1 198 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: Condition
##
## Response transformation matrix:
## (Intercept)
## PainDis 1
## PainCat 1
## PainIntens 1
##
## Sum of squares and products for the hypothesis:
## (Intercept)
## (Intercept) 47135.37
##
## Multivariate Tests: Condition
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.4353138 152.6372 1 198 < 2.22e-16 ***
## Wilks 1 0.5646862 152.6372 1 198 < 2.22e-16 ***
## Hotelling-Lawley 1 0.7708949 152.6372 1 198 < 2.22e-16 ***
## Roy 1 0.7708949 152.6372 1 198 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: task_pain
##
## Response transformation matrix:
## task_pain1 task_pain2
## PainDis 1 0
## PainCat 0 1
## PainIntens -1 -1
##
## Sum of squares and products for the hypothesis:
## task_pain1 task_pain2
## task_pain1 42.56113 -104.5413
## task_pain2 -104.54129 256.7808
##
## Multivariate Tests: task_pain
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.0168489 1.688054 2 197 0.18754
## Wilks 1 0.9831511 1.688054 2 197 0.18754
## Hotelling-Lawley 1 0.0171376 1.688054 2 197 0.18754
## Roy 1 0.0171376 1.688054 2 197 0.18754
##
## ------------------------------------------
##
## Term: Condition:task_pain
##
## Response transformation matrix:
## task_pain1 task_pain2
## PainDis 1 0
## PainCat 0 1
## PainIntens -1 -1
##
## Sum of squares and products for the hypothesis:
## task_pain1 task_pain2
## task_pain1 1445.909 -6203.972
## task_pain2 -6203.972 26619.429
##
## Multivariate Tests: Condition:task_pain
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.5895947 141.5066 2 197 < 2.22e-16 ***
## Wilks 1 0.4104053 141.5066 2 197 < 2.22e-16 ***
## Hotelling-Lawley 1 1.4366154 141.5066 2 197 < 2.22e-16 ***
## Roy 1 1.4366154 141.5066 2 197 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 1487141 1 20381 198 14447.3069 <2e-16 ***
## Condition 15712 1 20381 198 152.6372 <2e-16 ***
## task_pain 269 2 38752 396 1.3757 0.2539
## Condition:task_pain 22846 2 38752 396 116.7316 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## task_pain 0.96267 0.023568
## Condition:task_pain 0.96267 0.023568
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## task_pain 0.96401 0.2539
## Condition:task_pain 0.96401 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## task_pain 0.9733047 2.538584e-01
## Condition:task_pain 0.9733047 1.428276e-39
### 4. Complex mixed between-within with multivariate effect estimation - posthocs
# Now, the "idata" solution we used above, unfortunately doesn't *name*
# the repeated measures factor
# This means we don't have a factor that we can easily pull into any
# post-hoc approach.
# That means we have to revert to the long data set we just created
# And rerun the between-within interaction analysis in a different way
# so we can grab posthocs
#First we have to create three "melted" long data sets, one for each outcome
(Wk03a_cogment <- Wk03a[, c(2,3,4,5,6)])
## ID Condition MOCA BDI STAI
## 1 1 -1 42.29156 36.06104 60.288157
## 2 2 -1 29.04060 60.75548 31.914466
## 3 3 -1 57.59534 50.02032 37.155894
<truncated for length>
(Wk03a_phys <- Wk03a[, c(2,3,7,8,9)])
## ID Condition WalkSpeed CalfStrength V02
## 1 1 -1 44.05973 51.29967 45.54295
## 2 2 -1 46.59351 41.14846 39.85025
## 3 3 -1 35.49083 47.94157 33.87827
<truncated for length>
(Wk03a_pain <- Wk03a[, c(2,3,10,11,12)])
## ID Condition PainDis PainCat PainIntens
## 1 1 -1 56.23380 36.08248 40.76668
## 2 2 -1 52.15265 33.64651 60.19698
## 3 3 -1 47.50576 39.05012 53.52304
<truncated for length>
library(reshape2)
melt.cogment<-melt(Wk03a_cogment, id=c("ID","Condition"))
head(melt.cogment)
## ID Condition variable value
## 1 1 -1 MOCA 42.29156
## 2 2 -1 MOCA 29.04060
## 3 3 -1 MOCA 57.59534
## 4 4 -1 MOCA 33.71176
## 5 5 -1 MOCA 54.47724
## 6 6 -1 MOCA 48.87171
melt.cogment$task<-melt.cogment$variable
melt.cogment$score<-melt.cogment$value
head(melt.cogment)
## ID Condition variable value task score
## 1 1 -1 MOCA 42.29156 MOCA 42.29156
## 2 2 -1 MOCA 29.04060 MOCA 29.04060
## 3 3 -1 MOCA 57.59534 MOCA 57.59534
## 4 4 -1 MOCA 33.71176 MOCA 33.71176
## 5 5 -1 MOCA 54.47724 MOCA 54.47724
## 6 6 -1 MOCA 48.87171 MOCA 48.87171
melt.cogment$variable<-NULL
melt.cogment$value<-NULL
head(melt.cogment)
## ID Condition task score
## 1 1 -1 MOCA 42.29156
## 2 2 -1 MOCA 29.04060
## 3 3 -1 MOCA 57.59534
## 4 4 -1 MOCA 33.71176
## 5 5 -1 MOCA 54.47724
## 6 6 -1 MOCA 48.87171
melt.phys<-melt(Wk03a_phys, id=c("ID","Condition"))
head(melt.phys)
## ID Condition variable value
## 1 1 -1 WalkSpeed 44.05973
## 2 2 -1 WalkSpeed 46.59351
## 3 3 -1 WalkSpeed 35.49083
## 4 4 -1 WalkSpeed 22.67779
## 5 5 -1 WalkSpeed 35.36213
## 6 6 -1 WalkSpeed 51.92100
melt.phys$task<-melt.phys$variable
melt.phys$score<-melt.phys$value
head(melt.phys)
## ID Condition variable value task score
## 1 1 -1 WalkSpeed 44.05973 WalkSpeed 44.05973
## 2 2 -1 WalkSpeed 46.59351 WalkSpeed 46.59351
## 3 3 -1 WalkSpeed 35.49083 WalkSpeed 35.49083
## 4 4 -1 WalkSpeed 22.67779 WalkSpeed 22.67779
## 5 5 -1 WalkSpeed 35.36213 WalkSpeed 35.36213
## 6 6 -1 WalkSpeed 51.92100 WalkSpeed 51.92100
melt.phys$variable<-NULL
melt.phys$value<-NULL
head(melt.phys)
## ID Condition task score
## 1 1 -1 WalkSpeed 44.05973
## 2 2 -1 WalkSpeed 46.59351
## 3 3 -1 WalkSpeed 35.49083
## 4 4 -1 WalkSpeed 22.67779
## 5 5 -1 WalkSpeed 35.36213
## 6 6 -1 WalkSpeed 51.92100
melt.pain<-melt(Wk03a_pain, id=c("ID","Condition"))
head(melt.pain)
## ID Condition variable value
## 1 1 -1 PainDis 56.23380
## 2 2 -1 PainDis 52.15265
## 3 3 -1 PainDis 47.50576
## 4 4 -1 PainDis 63.64303
## 5 5 -1 PainDis 51.19453
## 6 6 -1 PainDis 57.53720
melt.pain$task<-melt.pain$variable
melt.pain$score<-melt.pain$value
head(melt.pain)
## ID Condition variable value task score
## 1 1 -1 PainDis 56.23380 PainDis 56.23380
## 2 2 -1 PainDis 52.15265 PainDis 52.15265
## 3 3 -1 PainDis 47.50576 PainDis 47.50576
## 4 4 -1 PainDis 63.64303 PainDis 63.64303
## 5 5 -1 PainDis 51.19453 PainDis 51.19453
## 6 6 -1 PainDis 57.53720 PainDis 57.53720
melt.pain$variable<-NULL
melt.pain$value<-NULL
head(melt.pain)
## ID Condition task score
## 1 1 -1 PainDis 56.23380
## 2 2 -1 PainDis 52.15265
## 3 3 -1 PainDis 47.50576
## 4 4 -1 PainDis 63.64303
## 5 5 -1 PainDis 51.19453
## 6 6 -1 PainDis 57.53720
##### Build model #####
library(nlme)
head(melt.cogment)
## ID Condition task score
## 1 1 -1 MOCA 42.29156
## 2 2 -1 MOCA 29.04060
## 3 3 -1 MOCA 57.59534
## 4 4 -1 MOCA 33.71176
## 5 5 -1 MOCA 54.47724
## 6 6 -1 MOCA 48.87171
model1 = lme(score ~ Condition*task,
random = ~1| ID ,
data=melt.cogment,
method="REML")
# The produces same f-statistics and p-values as idata above
anova.lme(model1)
## numDF denDF F-value p-value
## (Intercept) 1 396 13971.017 <.0001
## Condition 1 198 153.643 <.0001
## task 2 396 0.016 0.9843
## Condition:task 2 396 76.768 <.0001
model2 = lme(score ~ Condition*task,
random = ~1| ID ,
data=melt.phys,
method="REML")
# The produces same f-statistics and p-values as idata above
anova.lme(model2)
## numDF denDF F-value p-value
## (Intercept) 1 396 14119.225 <.0001
## Condition 1 198 110.904 <.0001
## task 2 396 1.273 0.2812
## Condition:task 2 396 41.897 <.0001
model3 = lme(score ~ Condition*task,
random = ~1| ID ,
data=melt.pain,
method="REML")
# The produces same f-statistics and p-values as idata above
anova.lme(model3)
## numDF denDF F-value p-value
## (Intercept) 1 396 14447.305 <.0001
## Condition 1 198 152.637 <.0001
## task 2 396 1.376 0.2539
## Condition:task 2 396 116.732 <.0001
##### Post-hoc tests #####
### These are kind of nice, because they are close to what
### SPSS does with simple effects decomposition, INCLUDING
### imposing homogeneity of variance on standard errors!
library(multcompView)
library(lsmeans)
## Loading required package: emmeans
## The 'lsmeans' package is now basically a front end for 'emmeans'.
## Users are encouraged to switch the rest of the way.
## See help('transition') for more information, including how to
## convert old 'lsmeans' objects and scripts to work with 'emmeans'.
#cogment
(leastsquare1 = lsmeans::lsmeans(model1,
pairwise ~ Condition*task,
adjust="none") ) ### Tukey-adjusted comparisons
## Warning in sweep(X, 1, sqrt(weights), "*"): STATS is longer than the extent of
## 'dim(x)[MARGIN]'
## $lsmeans
## Condition task lsmean SE df lower.CL upper.CL
## -1 MOCA 43.5 1 198 41.5 45.4
## 1 MOCA 56.6 1 198 54.6 58.6
## -1 BDI 51.6 1 198 49.6 53.6
## 1 BDI 48.9 1 198 46.9 50.8
## -1 STAI 39.6 1 198 37.6 41.6
## 1 STAI 60.7 1 198 58.8 62.7
##
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## (-1 MOCA) - 1 MOCA -13.14 1.41 198 -9.290 <.0001
## (-1 MOCA) - (-1 BDI) -8.11 1.39 396 -5.847 <.0001
## (-1 MOCA) - 1 BDI -5.38 1.41 198 -3.800 0.0002
## (-1 MOCA) - (-1 STAI) 3.88 1.39 396 2.797 0.0054
## (-1 MOCA) - 1 STAI -17.26 1.41 198 -12.198 <.0001
## 1 MOCA - (-1 BDI) 5.04 1.41 198 3.559 0.0005
## 1 MOCA - 1 BDI 7.77 1.39 396 5.602 <.0001
## 1 MOCA - (-1 STAI) 17.02 1.41 198 12.032 <.0001
## 1 MOCA - 1 STAI -4.11 1.39 396 -2.967 0.0032
## (-1 BDI) - 1 BDI 2.73 1.41 198 1.931 0.0549
## (-1 BDI) - (-1 STAI) 11.99 1.39 396 8.645 <.0001
## (-1 BDI) - 1 STAI -9.15 1.41 198 -6.467 <.0001
## 1 BDI - (-1 STAI) 9.26 1.41 198 6.541 <.0001
## 1 BDI - 1 STAI -11.88 1.39 396 -8.569 <.0001
## (-1 STAI) - 1 STAI -21.14 1.41 198 -14.940 <.0001
##
## Degrees-of-freedom method: containment
#lsmeans adjust options include
#"tukey", "scheffe", "sidak", "bonferroni", "dunnettx", "mvt", and "none"
#These are tricky to look at, but they have the information we need
#Some contrasts compare different treatments at pretest or at posttest
#These will give us the equivalent of SPSS simple effects tables
leastsquare1$contrasts
## contrast estimate SE df t.ratio p.value
## (-1 MOCA) - 1 MOCA -13.14 1.41 198 -9.290 <.0001
## (-1 MOCA) - (-1 BDI) -8.11 1.39 396 -5.847 <.0001
## (-1 MOCA) - 1 BDI -5.38 1.41 198 -3.800 0.0002
## (-1 MOCA) - (-1 STAI) 3.88 1.39 396 2.797 0.0054
## (-1 MOCA) - 1 STAI -17.26 1.41 198 -12.198 <.0001
## 1 MOCA - (-1 BDI) 5.04 1.41 198 3.559 0.0005
## 1 MOCA - 1 BDI 7.77 1.39 396 5.602 <.0001
## 1 MOCA - (-1 STAI) 17.02 1.41 198 12.032 <.0001
## 1 MOCA - 1 STAI -4.11 1.39 396 -2.967 0.0032
## (-1 BDI) - 1 BDI 2.73 1.41 198 1.931 0.0549
## (-1 BDI) - (-1 STAI) 11.99 1.39 396 8.645 <.0001
## (-1 BDI) - 1 STAI -9.15 1.41 198 -6.467 <.0001
## 1 BDI - (-1 STAI) 9.26 1.41 198 6.541 <.0001
## 1 BDI - 1 STAI -11.88 1.39 396 -8.569 <.0001
## (-1 STAI) - 1 STAI -21.14 1.41 198 -14.940 <.0001
##
## Degrees-of-freedom method: containment
#To slice the interaction in the opposite direction
(leastsquare1b = lsmeans::lsmeans(model1,
pairwise ~ task*Condition,
adjust="none") ) ### Tukey-adjusted comparisons
## Warning in sweep(X, 1, sqrt(weights), "*"): STATS is longer than the extent of
## 'dim(x)[MARGIN]'
## $lsmeans
## task Condition lsmean SE df lower.CL upper.CL
## MOCA -1 43.5 1 198 41.5 45.4
## BDI -1 51.6 1 198 49.6 53.6
## STAI -1 39.6 1 198 37.6 41.6
## MOCA 1 56.6 1 198 54.6 58.6
## BDI 1 48.9 1 198 46.9 50.8
## STAI 1 60.7 1 198 58.8 62.7
##
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## (MOCA -1) - (BDI -1) -8.11 1.39 396 -5.847 <.0001
## (MOCA -1) - (STAI -1) 3.88 1.39 396 2.797 0.0054
## (MOCA -1) - MOCA 1 -13.14 1.41 198 -9.290 <.0001
## (MOCA -1) - BDI 1 -5.38 1.41 198 -3.800 0.0002
## (MOCA -1) - STAI 1 -17.26 1.41 198 -12.198 <.0001
## (BDI -1) - (STAI -1) 11.99 1.39 396 8.645 <.0001
## (BDI -1) - MOCA 1 -5.04 1.41 198 -3.559 0.0005
## (BDI -1) - BDI 1 2.73 1.41 198 1.931 0.0549
## (BDI -1) - STAI 1 -9.15 1.41 198 -6.467 <.0001
## (STAI -1) - MOCA 1 -17.02 1.41 198 -12.032 <.0001
## (STAI -1) - BDI 1 -9.26 1.41 198 -6.541 <.0001
## (STAI -1) - STAI 1 -21.14 1.41 198 -14.940 <.0001
## MOCA 1 - BDI 1 7.77 1.39 396 5.602 <.0001
## MOCA 1 - STAI 1 -4.11 1.39 396 -2.967 0.0032
## BDI 1 - STAI 1 -11.88 1.39 396 -8.569 <.0001
##
## Degrees-of-freedom method: containment
leastsquare1b$contrasts
## contrast estimate SE df t.ratio p.value
## (MOCA -1) - (BDI -1) -8.11 1.39 396 -5.847 <.0001
## (MOCA -1) - (STAI -1) 3.88 1.39 396 2.797 0.0054
## (MOCA -1) - MOCA 1 -13.14 1.41 198 -9.290 <.0001
## (MOCA -1) - BDI 1 -5.38 1.41 198 -3.800 0.0002
## (MOCA -1) - STAI 1 -17.26 1.41 198 -12.198 <.0001
## (BDI -1) - (STAI -1) 11.99 1.39 396 8.645 <.0001
## (BDI -1) - MOCA 1 -5.04 1.41 198 -3.559 0.0005
## (BDI -1) - BDI 1 2.73 1.41 198 1.931 0.0549
## (BDI -1) - STAI 1 -9.15 1.41 198 -6.467 <.0001
## (STAI -1) - MOCA 1 -17.02 1.41 198 -12.032 <.0001
## (STAI -1) - BDI 1 -9.26 1.41 198 -6.541 <.0001
## (STAI -1) - STAI 1 -21.14 1.41 198 -14.940 <.0001
## MOCA 1 - BDI 1 7.77 1.39 396 5.602 <.0001
## MOCA 1 - STAI 1 -4.11 1.39 396 -2.967 0.0032
## BDI 1 - STAI 1 -11.88 1.39 396 -8.569 <.0001
##
## Degrees-of-freedom method: containment
#phys
(leastsquare2 = lsmeans::lsmeans(model2,
pairwise ~ Condition*task,
adjust="none") ) ### Tukey-adjusted comparisons
## Warning in sweep(X, 1, sqrt(weights), "*"): STATS is longer than the extent of
## 'dim(x)[MARGIN]'
## $lsmeans
## Condition task lsmean SE df lower.CL upper.CL
## -1 WalkSpeed 42.1 1.04 198 40.1 44.2
## 1 WalkSpeed 61.7 1.04 198 59.7 63.8
## -1 CalfStrength 49.7 1.04 198 47.6 51.7
## 1 CalfStrength 51.1 1.04 198 49.1 53.2
## -1 V02 47.6 1.04 198 45.5 49.6
## 1 V02 53.7 1.04 198 51.6 55.7
##
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## (-1 WalkSpeed) - 1 WalkSpeed -19.61 1.47 198 -13.358 <.0001
## (-1 WalkSpeed) - (-1 CalfStrength) -7.56 1.46 396 -5.186 <.0001
## (-1 WalkSpeed) - 1 CalfStrength -9.01 1.47 198 -6.138 <.0001
## (-1 WalkSpeed) - (-1 V02) -5.47 1.46 396 -3.748 0.0002
## (-1 WalkSpeed) - 1 V02 -11.53 1.47 198 -7.853 <.0001
## 1 WalkSpeed - (-1 CalfStrength) 12.05 1.47 198 8.206 <.0001
## 1 WalkSpeed - 1 CalfStrength 10.60 1.46 396 7.268 <.0001
## 1 WalkSpeed - (-1 V02) 14.14 1.47 198 9.635 <.0001
## 1 WalkSpeed - 1 V02 8.08 1.46 396 5.541 <.0001
## (-1 CalfStrength) - 1 CalfStrength -1.45 1.47 198 -0.985 0.3256
## (-1 CalfStrength) - (-1 V02) 2.10 1.46 396 1.438 0.1512
## (-1 CalfStrength) - 1 V02 -3.96 1.47 198 -2.701 0.0075
## 1 CalfStrength - (-1 V02) 3.54 1.47 198 2.414 0.0167
## 1 CalfStrength - 1 V02 -2.52 1.46 396 -1.727 0.0850
## (-1 V02) - 1 V02 -6.06 1.47 198 -4.130 0.0001
##
## Degrees-of-freedom method: containment
leastsquare2$contrasts
## contrast estimate SE df t.ratio p.value
## (-1 WalkSpeed) - 1 WalkSpeed -19.61 1.47 198 -13.358 <.0001
## (-1 WalkSpeed) - (-1 CalfStrength) -7.56 1.46 396 -5.186 <.0001
## (-1 WalkSpeed) - 1 CalfStrength -9.01 1.47 198 -6.138 <.0001
## (-1 WalkSpeed) - (-1 V02) -5.47 1.46 396 -3.748 0.0002
## (-1 WalkSpeed) - 1 V02 -11.53 1.47 198 -7.853 <.0001
## 1 WalkSpeed - (-1 CalfStrength) 12.05 1.47 198 8.206 <.0001
## 1 WalkSpeed - 1 CalfStrength 10.60 1.46 396 7.268 <.0001
## 1 WalkSpeed - (-1 V02) 14.14 1.47 198 9.635 <.0001
## 1 WalkSpeed - 1 V02 8.08 1.46 396 5.541 <.0001
## (-1 CalfStrength) - 1 CalfStrength -1.45 1.47 198 -0.985 0.3256
## (-1 CalfStrength) - (-1 V02) 2.10 1.46 396 1.438 0.1512
## (-1 CalfStrength) - 1 V02 -3.96 1.47 198 -2.701 0.0075
## 1 CalfStrength - (-1 V02) 3.54 1.47 198 2.414 0.0167
## 1 CalfStrength - 1 V02 -2.52 1.46 396 -1.727 0.0850
## (-1 V02) - 1 V02 -6.06 1.47 198 -4.130 0.0001
##
## Degrees-of-freedom method: containment
#To slice the interaction in the opposite direction
(leastsquare2b = lsmeans::lsmeans(model2,
pairwise ~ task*Condition,
adjust="none") ) ### Tukey-adjusted comparisons
## Warning in sweep(X, 1, sqrt(weights), "*"): STATS is longer than the extent of
## 'dim(x)[MARGIN]'
## $lsmeans
## task Condition lsmean SE df lower.CL upper.CL
## WalkSpeed -1 42.1 1.04 198 40.1 44.2
## CalfStrength -1 49.7 1.04 198 47.6 51.7
## V02 -1 47.6 1.04 198 45.5 49.6
## WalkSpeed 1 61.7 1.04 198 59.7 63.8
## CalfStrength 1 51.1 1.04 198 49.1 53.2
## V02 1 53.7 1.04 198 51.6 55.7
##
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## (WalkSpeed -1) - (CalfStrength -1) -7.56 1.46 396 -5.186 <.0001
## (WalkSpeed -1) - (V02 -1) -5.47 1.46 396 -3.748 0.0002
## (WalkSpeed -1) - WalkSpeed 1 -19.61 1.47 198 -13.358 <.0001
## (WalkSpeed -1) - CalfStrength 1 -9.01 1.47 198 -6.138 <.0001
## (WalkSpeed -1) - V02 1 -11.53 1.47 198 -7.853 <.0001
## (CalfStrength -1) - (V02 -1) 2.10 1.46 396 1.438 0.1512
## (CalfStrength -1) - WalkSpeed 1 -12.05 1.47 198 -8.206 <.0001
## (CalfStrength -1) - CalfStrength 1 -1.45 1.47 198 -0.985 0.3256
## (CalfStrength -1) - V02 1 -3.96 1.47 198 -2.701 0.0075
## (V02 -1) - WalkSpeed 1 -14.14 1.47 198 -9.635 <.0001
## (V02 -1) - CalfStrength 1 -3.54 1.47 198 -2.414 0.0167
## (V02 -1) - V02 1 -6.06 1.47 198 -4.130 0.0001
## WalkSpeed 1 - CalfStrength 1 10.60 1.46 396 7.268 <.0001
## WalkSpeed 1 - V02 1 8.08 1.46 396 5.541 <.0001
## CalfStrength 1 - V02 1 -2.52 1.46 396 -1.727 0.0850
##
## Degrees-of-freedom method: containment
leastsquare2b$contrasts
## contrast estimate SE df t.ratio p.value
## (WalkSpeed -1) - (CalfStrength -1) -7.56 1.46 396 -5.186 <.0001
## (WalkSpeed -1) - (V02 -1) -5.47 1.46 396 -3.748 0.0002
## (WalkSpeed -1) - WalkSpeed 1 -19.61 1.47 198 -13.358 <.0001
## (WalkSpeed -1) - CalfStrength 1 -9.01 1.47 198 -6.138 <.0001
## (WalkSpeed -1) - V02 1 -11.53 1.47 198 -7.853 <.0001
## (CalfStrength -1) - (V02 -1) 2.10 1.46 396 1.438 0.1512
## (CalfStrength -1) - WalkSpeed 1 -12.05 1.47 198 -8.206 <.0001
## (CalfStrength -1) - CalfStrength 1 -1.45 1.47 198 -0.985 0.3256
## (CalfStrength -1) - V02 1 -3.96 1.47 198 -2.701 0.0075
## (V02 -1) - WalkSpeed 1 -14.14 1.47 198 -9.635 <.0001
## (V02 -1) - CalfStrength 1 -3.54 1.47 198 -2.414 0.0167
## (V02 -1) - V02 1 -6.06 1.47 198 -4.130 0.0001
## WalkSpeed 1 - CalfStrength 1 10.60 1.46 396 7.268 <.0001
## WalkSpeed 1 - V02 1 8.08 1.46 396 5.541 <.0001
## CalfStrength 1 - V02 1 -2.52 1.46 396 -1.727 0.0850
##
## Degrees-of-freedom method: containment
#pain
(leastsquare3 = lsmeans::lsmeans(model3,
pairwise ~ Condition*task,
adjust="none") ) ### Tukey-adjusted comparisons
## Warning in sweep(X, 1, sqrt(weights), "*"): STATS is longer than the extent of
## 'dim(x)[MARGIN]'
## $lsmeans
## Condition task lsmean SE df lower.CL upper.CL
## -1 PainDis 49.6 0.998 198 47.7 51.6
## 1 PainDis 48.6 0.998 198 46.6 50.5
## -1 PainCat 37.0 0.998 198 35.0 39.0
## 1 PainCat 64.4 0.998 198 62.4 66.4
## -1 PainIntens 47.4 0.998 198 45.4 49.4
## 1 PainIntens 51.7 0.998 198 49.8 53.7
##
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## (-1 PainDis) - 1 PainDis 1.04 1.41 198 0.738 0.4612
## (-1 PainDis) - (-1 PainCat) 12.63 1.40 396 9.029 <.0001
## (-1 PainDis) - 1 PainCat -14.78 1.41 198 -10.473 <.0001
## (-1 PainDis) - (-1 PainIntens) 2.23 1.40 396 1.592 0.1121
## (-1 PainDis) - 1 PainIntens -2.11 1.41 198 -1.494 0.1367
## 1 PainDis - (-1 PainCat) 11.59 1.41 198 8.213 <.0001
## 1 PainDis - 1 PainCat -15.82 1.40 396 -11.308 <.0001
## 1 PainDis - (-1 PainIntens) 1.19 1.41 198 0.840 0.4017
## 1 PainDis - 1 PainIntens -3.15 1.40 396 -2.252 0.0249
## (-1 PainCat) - 1 PainCat -27.41 1.41 198 -19.425 <.0001
## (-1 PainCat) - (-1 PainIntens) -10.40 1.40 396 -7.437 <.0001
## (-1 PainCat) - 1 PainIntens -14.74 1.41 198 -10.446 <.0001
## 1 PainCat - (-1 PainIntens) 17.01 1.41 198 12.052 <.0001
## 1 PainCat - 1 PainIntens 12.67 1.40 396 9.056 <.0001
## (-1 PainIntens) - 1 PainIntens -4.34 1.41 198 -3.073 0.0024
##
## Degrees-of-freedom method: containment
leastsquare3$contrasts
## contrast estimate SE df t.ratio p.value
## (-1 PainDis) - 1 PainDis 1.04 1.41 198 0.738 0.4612
## (-1 PainDis) - (-1 PainCat) 12.63 1.40 396 9.029 <.0001
## (-1 PainDis) - 1 PainCat -14.78 1.41 198 -10.473 <.0001
## (-1 PainDis) - (-1 PainIntens) 2.23 1.40 396 1.592 0.1121
## (-1 PainDis) - 1 PainIntens -2.11 1.41 198 -1.494 0.1367
## 1 PainDis - (-1 PainCat) 11.59 1.41 198 8.213 <.0001
## 1 PainDis - 1 PainCat -15.82 1.40 396 -11.308 <.0001
## 1 PainDis - (-1 PainIntens) 1.19 1.41 198 0.840 0.4017
## 1 PainDis - 1 PainIntens -3.15 1.40 396 -2.252 0.0249
## (-1 PainCat) - 1 PainCat -27.41 1.41 198 -19.425 <.0001
## (-1 PainCat) - (-1 PainIntens) -10.40 1.40 396 -7.437 <.0001
## (-1 PainCat) - 1 PainIntens -14.74 1.41 198 -10.446 <.0001
## 1 PainCat - (-1 PainIntens) 17.01 1.41 198 12.052 <.0001
## 1 PainCat - 1 PainIntens 12.67 1.40 396 9.056 <.0001
## (-1 PainIntens) - 1 PainIntens -4.34 1.41 198 -3.073 0.0024
##
## Degrees-of-freedom method: containment
#To slice the interaction in the opposite direction
(leastsquare3b = lsmeans::lsmeans(model3,
pairwise ~ task*Condition,
adjust="none") ) ### Tukey-adjusted comparisons
## Warning in sweep(X, 1, sqrt(weights), "*"): STATS is longer than the extent of
## 'dim(x)[MARGIN]'
## $lsmeans
## task Condition lsmean SE df lower.CL upper.CL
## PainDis -1 49.6 0.998 198 47.7 51.6
## PainCat -1 37.0 0.998 198 35.0 39.0
## PainIntens -1 47.4 0.998 198 45.4 49.4
## PainDis 1 48.6 0.998 198 46.6 50.5
## PainCat 1 64.4 0.998 198 62.4 66.4
## PainIntens 1 51.7 0.998 198 49.8 53.7
##
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## (PainDis -1) - (PainCat -1) 12.63 1.40 396 9.029 <.0001
## (PainDis -1) - (PainIntens -1) 2.23 1.40 396 1.592 0.1121
## (PainDis -1) - PainDis 1 1.04 1.41 198 0.738 0.4612
## (PainDis -1) - PainCat 1 -14.78 1.41 198 -10.473 <.0001
## (PainDis -1) - PainIntens 1 -2.11 1.41 198 -1.494 0.1367
## (PainCat -1) - (PainIntens -1) -10.40 1.40 396 -7.437 <.0001
## (PainCat -1) - PainDis 1 -11.59 1.41 198 -8.213 <.0001
## (PainCat -1) - PainCat 1 -27.41 1.41 198 -19.425 <.0001
## (PainCat -1) - PainIntens 1 -14.74 1.41 198 -10.446 <.0001
## (PainIntens -1) - PainDis 1 -1.19 1.41 198 -0.840 0.4017
## (PainIntens -1) - PainCat 1 -17.01 1.41 198 -12.052 <.0001
## (PainIntens -1) - PainIntens 1 -4.34 1.41 198 -3.073 0.0024
## PainDis 1 - PainCat 1 -15.82 1.40 396 -11.308 <.0001
## PainDis 1 - PainIntens 1 -3.15 1.40 396 -2.252 0.0249
## PainCat 1 - PainIntens 1 12.67 1.40 396 9.056 <.0001
##
## Degrees-of-freedom method: containment
leastsquare3b$contrasts
## contrast estimate SE df t.ratio p.value
## (PainDis -1) - (PainCat -1) 12.63 1.40 396 9.029 <.0001
## (PainDis -1) - (PainIntens -1) 2.23 1.40 396 1.592 0.1121
## (PainDis -1) - PainDis 1 1.04 1.41 198 0.738 0.4612
## (PainDis -1) - PainCat 1 -14.78 1.41 198 -10.473 <.0001
## (PainDis -1) - PainIntens 1 -2.11 1.41 198 -1.494 0.1367
## (PainCat -1) - (PainIntens -1) -10.40 1.40 396 -7.437 <.0001
## (PainCat -1) - PainDis 1 -11.59 1.41 198 -8.213 <.0001
## (PainCat -1) - PainCat 1 -27.41 1.41 198 -19.425 <.0001
## (PainCat -1) - PainIntens 1 -14.74 1.41 198 -10.446 <.0001
## (PainIntens -1) - PainDis 1 -1.19 1.41 198 -0.840 0.4017
## (PainIntens -1) - PainCat 1 -17.01 1.41 198 -12.052 <.0001
## (PainIntens -1) - PainIntens 1 -4.34 1.41 198 -3.073 0.0024
## PainDis 1 - PainCat 1 -15.82 1.40 396 -11.308 <.0001
## PainDis 1 - PainIntens 1 -3.15 1.40 396 -2.252 0.0249
## PainCat 1 - PainIntens 1 12.67 1.40 396 9.056 <.0001
##
## Degrees-of-freedom method: containment
### 5. Basic interaction plot, no error bars
------------------------------------------------------------------------------
### Begin with simple interaction plots
### Note this has no error bars;
### What follows is a more complex approach from CLP 6528 which has a detailed
### helper function to get correct repeated measures error
### bars for within factors
require(graphics)
## [1] 1
head(melt.cogment)
## ID Condition task score
## 1 1 -1 MOCA 42.29156
## 2 2 -1 MOCA 29.04060
## 3 3 -1 MOCA 57.59534
## 4 4 -1 MOCA 33.71176
## 5 5 -1 MOCA 54.47724
## 6 6 -1 MOCA 48.87171
with(melt.cogment, {
interaction.plot(task, Condition, score, fixed = TRUE)})
with(melt.phys, {
interaction.plot(task, Condition, score, fixed = TRUE)})
with(melt.pain, {
interaction.plot(task, Condition, score, fixed = TRUE)})
### 5. Rich interaction plot, with error bars
# First, collapse the data using summarySEwithin
#(defined here
# both of the helper functions below must be entered before the function is called here).
############################################################################################
# Just run everything between these helper function bars
# Begin helper function
############################################################################################
#Helper Functions
## Gives count, mean, standard deviation, standard error of the mean, and confidence interval (default 95%).
## data: a data frame.
## measurevar: the name of a column that contains the variable to be summariezed
## groupvars: a vector containing names of columns that contain grouping variables
## na.rm: a boolean that indicates whether to ignore NA's
## conf.interval: the percent range of the confidence interval (default is 95%)
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
conf.interval=.95, .drop=TRUE) {
library(plyr)
# New version of length which can handle NA's: if na.rm==T, don't count them
length2 <- function (x, na.rm=FALSE) {
if (na.rm) sum(!is.na(x))
else length(x)
}
# This does the summary. For each group's data frame, return a vector with
# N, mean, and sd
datac <- ddply(data, groupvars, .drop=.drop,
.fun = function(xx, col) {
c(N = length2(xx[[col]], na.rm=na.rm),
mean = mean (xx[[col]], na.rm=na.rm),
sd = sd (xx[[col]], na.rm=na.rm)
)
},
measurevar
)
# Rename the "mean" column
datac <- rename(datac, c("mean" = measurevar))
datac$se <- datac$sd / sqrt(datac$N) # Calculate standard error of the mean
# Confidence interval multiplier for standard error
# Calculate t-statistic for confidence interval:
# e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
ciMult <- qt(conf.interval/2 + .5, datac$N-1)
datac$ci <- datac$se * ciMult
return(datac)
}
## Norms the data within specified groups in a data frame; it normalizes each
## subject (identified by idvar) so that they have the same mean, within each group
## specified by betweenvars.
## data: a data frame.
## idvar: the name of a column that identifies each subject (or matched subjects)
## measurevar: the name of a column that contains the variable to be summariezed
## betweenvars: a vector containing names of columns that are between-subjects variables
## na.rm: a boolean that indicates whether to ignore NA's
normDataWithin <- function(data=NULL, idvar, measurevar, betweenvars=NULL,
na.rm=FALSE, .drop=TRUE) {
library(plyr)
# Measure var on left, idvar + between vars on right of formula.
data.subjMean <- ddply(data, c(idvar, betweenvars), .drop=.drop,
.fun = function(xx, col, na.rm) {
c(subjMean = mean(xx[,col], na.rm=na.rm))
},
measurevar,
na.rm
)
# Put the subject means with original data
data <- merge(data, data.subjMean)
# Get the normalized data in a new column
measureNormedVar <- paste(measurevar, "_norm", sep="")
data[,measureNormedVar] <- data[,measurevar] - data[,"subjMean"] +
mean(data[,measurevar], na.rm=na.rm)
# Remove this subject mean column
data$subjMean <- NULL
return(data)
}
## Summarizes data, handling within-subjects variables by removing inter-subject variability.
## It will still work if there are no within-S variables.
## Gives count, un-normed mean, normed mean (with same between-group mean),
## standard deviation, standard error of the mean, and confidence interval.
## If there are within-subject variables, calculate adjusted values using method from Morey (2008).
## data: a data frame.
## measurevar: the name of a column that contains the variable to be summariezed
## betweenvars: a vector containing names of columns that are between-subjects variables
## withinvars: a vector containing names of columns that are within-subjects variables
## idvar: the name of a column that identifies each subject (or matched subjects)
## na.rm: a boolean that indicates whether to ignore NA's
## conf.interval: the percent range of the confidence interval (default is 95%)
summarySEwithin <- function(data=NULL, measurevar, betweenvars=NULL, withinvars=NULL,
idvar=NULL, na.rm=FALSE, conf.interval=.95, .drop=TRUE) {
# Ensure that the betweenvars and withinvars are factors
factorvars <- vapply(data[, c(betweenvars, withinvars), drop=FALSE],
FUN=is.factor, FUN.VALUE=logical(1))
if (!all(factorvars)) {
nonfactorvars <- names(factorvars)[!factorvars]
message("Automatically converting the following non-factors to factors: ",
paste(nonfactorvars, collapse = ", "))
data[nonfactorvars] <- lapply(data[nonfactorvars], factor)
}
# Get the means from the un-normed data
datac <- summarySE(data, measurevar, groupvars=c(betweenvars, withinvars),
na.rm=na.rm, conf.interval=conf.interval, .drop=.drop)
# Drop all the unused columns (these will be calculated with normed data)
datac$sd <- NULL
datac$se <- NULL
datac$ci <- NULL
# Norm each subject's data
ndata <- normDataWithin(data, idvar, measurevar, betweenvars, na.rm, .drop=.drop)
# This is the name of the new column
measurevar_n <- paste(measurevar, "_norm", sep="")
# Collapse the normed data - now we can treat between and within vars the same
ndatac <- summarySE(ndata, measurevar_n, groupvars=c(betweenvars, withinvars),
na.rm=na.rm, conf.interval=conf.interval, .drop=.drop)
# Apply correction from Morey (2008) to the standard error and confidence interval
# Get the product of the number of conditions of within-S variables
nWithinGroups <- prod(vapply(ndatac[,withinvars, drop=FALSE], FUN=nlevels,
FUN.VALUE=numeric(1)))
correctionFactor <- sqrt( nWithinGroups / (nWithinGroups-1) )
# Apply the correction factor
ndatac$sd <- ndatac$sd * correctionFactor
ndatac$se <- ndatac$se * correctionFactor
ndatac$ci <- ndatac$ci * correctionFactor
# Combine the un-normed means with the normed results
merge(datac, ndatac)
}
############################################################################################
# End helper function
############################################################################################
##Plot cogment
head(melt.cogment)
## ID Condition task score
## 1 1 -1 MOCA 42.29156
## 2 2 -1 MOCA 29.04060
## 3 3 -1 MOCA 57.59534
## 4 4 -1 MOCA 33.71176
## 5 5 -1 MOCA 54.47724
## 6 6 -1 MOCA 48.87171
#This creates a summary data set with means, stds, cis for each combination of conditions
cogrepc <- summarySEwithin(melt.cogment, measurevar="score", betweenvars="Condition",
withinvars="task",
idvar="ID", na.rm=FALSE, conf.interval=.95)
cogrepc
## Condition task N score score_norm sd se ci
## 1 -1 BDI 100 51.58296 56.84127 9.774080 0.9774080 1.939390
## 2 -1 MOCA 100 43.47408 48.73239 11.338775 1.1338775 2.249859
## 3 -1 STAI 100 39.59466 44.85297 9.828430 0.9828430 1.950174
## 4 1 BDI 100 48.85020 43.59189 8.441122 0.8441122 1.674902
## 5 1 MOCA 100 56.61832 51.36001 9.931792 0.9931792 1.970683
## 6 1 STAI 100 60.73304 55.47473 9.292674 0.9292674 1.843868
#The next section adds some more meaningful labels for plotting
library(expss)
##
## Use 'expss_output_viewer()' to display tables in the RStudio Viewer.
## To return to the console output, use 'expss_output_default()'.
##
## Attaching package: 'expss'
## The following object is masked from 'package:plyr':
##
## split_labels
## The following object is masked from 'package:car':
##
## recode
cogrepc = apply_labels(cogrepc,
task = "Task")
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:expss':
##
## between, compute, contains, first, last, na_if, recode, vars
## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## The following object is masked from 'package:nlme':
##
## collapse
## The following object is masked from 'package:car':
##
## recode
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
cogrepc$task2<-recode(cogrepc$task, PainDis = "Pain Disability", PainCat = "Pain Catastrophizing",
PainIntens = "Pain Intensity")
cogrepc$Condition2<-as.character(cogrepc$Condition)
cogrepc$Condition2<-recode(cogrepc$Condition2, "-1" = "marijuana", "1" = "opiate")
#Now we do the plot
#All those terms involving Gender ensure that separate lines and a legend are produced
#by gender
#scale_color_discrete "name" gives a title to your legend (not very intuitive)
library(ggplot2)
##
## Attaching package: 'ggplot2'
## The following object is masked from 'package:expss':
##
## vars
# Make the graph with the 95% confidence interval
ggplot(cogrepc, aes(x=task2, y=score, group=Condition2, colour=Condition2)) +
geom_line(aes(x=task2, y=score, group=Condition2, colour=Condition2)) +
scale_colour_discrete(name = "Condition") +
geom_errorbar(width=.1, aes(ymin=score-ci, ymax=score+ci)) +
geom_point(shape=21, size=3, fill="white") +
labs(x = "Task", y="Task Score") +
ylim(35,65) +
theme_bw()
# We should also do this for phys and pain, but will stop here for time.
# R may require you to shut down and reopen for each diagram, because the helper function gets
# "stuck" in memory with the last time it was run and sometimes can't be updated.