### MANCOVA with covariate adjusted posthocs
#install.packages("jmv")
library(jmv)
### Read the data
Wk03b<- read.csv("Class03_homework_2018b.csv")
head(Wk03b)
## ï..const ID Condition PSS DmSEff SIMS WM
## 1 1 1 -1 52.57735 66.10522 41.86421 59.43669
## 2 1 2 -1 27.95034 51.07199 58.17004 37.15727
## 3 1 3 -1 31.96956 55.64821 60.25409 35.50487
## 4 1 4 -1 38.03308 57.91234 46.41795 45.38161
## 5 1 5 -1 60.94577 62.95876 40.40288 43.78198
## 6 1 6 -1 59.59354 65.72283 33.15095 49.70716
#We will run it with and without the covariate, so we can compare
Wk03b$Condition<-as.factor(Wk03b$Condition)
Wk03b$Condition<-as.numeric(Wk03b$Condition)
head(Wk03b)
## ï..const ID Condition PSS DmSEff SIMS WM
## 1 1 1 1 52.57735 66.10522 41.86421 59.43669
## 2 1 2 1 27.95034 51.07199 58.17004 37.15727
## 3 1 3 1 31.96956 55.64821 60.25409 35.50487
## 4 1 4 1 38.03308 57.91234 46.41795 45.38161
## 5 1 5 1 60.94577 62.95876 40.40288 43.78198
## 6 1 6 1 59.59354 65.72283 33.15095 49.70716
options(contrasts = c("contr.helmert", "contr.poly"))
#Without covariate
mancova(data=Wk03b, deps = vars(PSS,DmSEff,SIMS), factors = Condition,
multivar = list("pillai", "wilks", "hotel", "roy"), boxM = TRUE,
shapiro = TRUE, qqPlot = TRUE)
##
## MANCOVA
##
## Multivariate Tests
## ---------------------------------------------------------------------------------------
## value F df1 df2 p
## ---------------------------------------------------------------------------------------
## Condition Pillai's Trace 0.3207947 5.476469 6 172 0.0000328
## Wilks' Lambda 0.6819781 5.976014 6 170 0.0000109
## Hotelling's Trace 0.4622568 6.471595 6 168 0.0000037
## Roy's Largest Root 0.4532870 12.99423 3 86 0.0000004
## ---------------------------------------------------------------------------------------
##
##
## Univariate Tests
## ----------------------------------------------------------------------------------------------------
## Dependent Variable Sum of Squares df Mean Square F p
## ----------------------------------------------------------------------------------------------------
## Condition PSS 1675.2892 2 837.64460 9.815231 0.0001433
## DmSEff 3074.5832 2 1537.29159 11.757021 0.0000302
## SIMS 223.6476 2 111.82382 1.085769 0.3421726
## Residuals PSS 7424.6931 87 85.34130
## DmSEff 11375.7017 87 130.75519
## SIMS 8960.1658 87 102.99041
## ----------------------------------------------------------------------------------------------------
##
##
## ASSUMPTION CHECKS
##
## Box's Homogeneity of Covariance Matrices Test
## ---------------------------------------------
## <U+03C7>² df p
## ---------------------------------------------
## 13.01614 12 0.3678687
## ---------------------------------------------
##
##
## Shapiro-Wilk Multivariate Normality Test
## ----------------------------------------
## W p
## ----------------------------------------
## 0.9801197 0.1841270
## ----------------------------------------
#With covariate
mancova(data=Wk03b, deps = vars(PSS,DmSEff,SIMS), factors = Condition,
covs = vars(WM),
multivar = list("pillai", "wilks", "hotel", "roy"), boxM = TRUE,
shapiro = TRUE, qqPlot = TRUE)
##
## MANCOVA
##
## Multivariate Tests
## -----------------------------------------------------------------------------------------
## value F df1 df2 p
## -----------------------------------------------------------------------------------------
## Condition Pillai's Trace 0.3576326 6.169706 6 170 0.0000071
## Wilks' Lambda 0.6468316 6.814681 6 168 0.0000017
## Hotelling's Trace 0.5390959 7.457493 6 166 0.0000004
## Roy's Largest Root 0.5259743 14.90260 3 85 < .0000001
##
## WM Pillai's Trace 0.4562081 23.490283 3 84 < .0000001
## Wilks' Lambda 0.5437919 23.490283 3 84 < .0000001
## Hotelling's Trace 0.8389387 23.490283 3 84 < .0000001
## Roy's Largest Root 0.8389387 23.49028 3 84 < .0000001
## -----------------------------------------------------------------------------------------
##
##
## Univariate Tests
## -------------------------------------------------------------------------------------------------------
## Dependent Variable Sum of Squares df Mean Square F p
## -------------------------------------------------------------------------------------------------------
## Condition PSS 1675.28920 2 837.64460 16.26289597 0.0000010
## DmSEff 3074.58319 2 1537.29159 11.93373652 0.0000267
## SIMS 223.64763 2 111.82382 1.07450786 0.3460049
## WM PSS 2995.13545 1 2995.13545 58.15064777 < .0000001
## DmSEff 297.27064 1 297.27064 2.30766205 0.1324055
## SIMS 10.16319 1 10.16319 0.09765741 0.7554156
## Residuals PSS 4429.55768 86 51.50648
## DmSEff 11078.43103 86 128.81897
## SIMS 8950.00257 86 104.06980
## -------------------------------------------------------------------------------------------------------
##
##
## ASSUMPTION CHECKS
##
## Box's Homogeneity of Covariance Matrices Test
## ---------------------------------------------
## <U+03C7>² df p
## ---------------------------------------------
## 13.01614 12 0.3678687
## ---------------------------------------------
##
##
## Shapiro-Wilk Multivariate Normality Test
## ----------------------------------------
## W p
## ----------------------------------------
## 0.9801197 0.1841270
## ----------------------------------------
#this is a second way to do the ANCOVA; If you put the covariate second,
#it yields identical results to the
#jmv package, and has the same discrepancies from SPSS
#If you put the covariate first, as below, the F-approximations for your IV
#match SPSS (!), but the df are still off, and the covariate now has
#different estimates from SPSS
#In the Univariates, "group' effects match SPSS perfectly, though covariate
#is off
attach(Wk03b)
Y <- cbind(PSS,DmSEff,SIMS)
man2 <- manova(Y ~ WM+Condition)
summary(man2, test="Pillai")
## Df Pillai approx F num Df den Df Pr(>F)
## WM 1 0.50036 28.3738 3 85 8.280e-13 ***
## Condition 1 0.23207 8.5622 3 85 4.957e-05 ***
## Residuals 87
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(man2, test="Wilks")
## Df Wilks approx F num Df den Df Pr(>F)
## WM 1 0.49964 28.3738 3 85 8.280e-13 ***
## Condition 1 0.76793 8.5622 3 85 4.957e-05 ***
## Residuals 87
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(man2, test="Hotelling-Lawley")
## Df Hotelling-Lawley approx F num Df den Df Pr(>F)
## WM 1 1.00143 28.3738 3 85 8.280e-13 ***
## Condition 1 0.30219 8.5622 3 85 4.957e-05 ***
## Residuals 87
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(man2, test="Roy")
## Df Roy approx F num Df den Df Pr(>F)
## WM 1 1.00143 28.3738 3 85 8.280e-13 ***
## Condition 1 0.30219 8.5622 3 85 4.957e-05 ***
## Residuals 87
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary.aov(man2) # for univariate statistics.
## Response PSS :
## Df Sum Sq Mean Sq F value Pr(>F)
## WM 1 4373.1 4373.1 83.4330 2.386e-14 ***
## Condition 1 166.9 166.9 3.1835 0.07787 .
## Residuals 87 4560.0 52.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response DmSEff :
## Df Sum Sq Mean Sq F value Pr(>F)
## WM 1 71.3 71.3 0.5599 0.4563
## Condition 1 3293.0 3293.0 25.8424 2.101e-06 ***
## Residuals 87 11086.0 127.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response SIMS :
## Df Sum Sq Mean Sq F value Pr(>F)
## WM 1 9.9 9.929 0.0964 0.7570
## Condition 1 211.3 211.343 2.0515 0.1556
## Residuals 87 8962.5 103.018
#Both sets of the above results look a little different from SPSS in the univariates
#Specifically, without covariate, the analysis matches SPSS
#With covariate, the *covariate* effects match SPSS, but the group doesn't
#Degrees of freedom are also much higher than SPSS (in both approaches if the cov is 2nd)
#If the cov is 1st, see above for greater congruence in the second approach
#A better way to get the univariates is to do them as three separate ANCOVAS
library(multcomp) #adds glht
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
##
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
##
## geyser
library(psych) #adds describBy
##
## Attaching package: 'psych'
## The following object is masked from 'package:jmv':
##
## pca
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
## The following object is masked from 'package:car':
##
## logit
#install.packages("effects")
library(effects) #adds covariate adjusted means
## Registered S3 methods overwritten by 'lme4':
## method from
## cooks.distance.influence.merMod car
## influence.merMod car
## dfbeta.influence.merMod car
## dfbetas.influence.merMod car
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
Wk03b$Condition<-as.factor(Wk03b$Condition)
#PSS
# NOTE: covariate goes first!! NOTE: there is
# no interaction
#without covariate
res1a <- aov(PSS ~ Condition, data = Wk03b)
summary(res1a)
## Df Sum Sq Mean Sq F value Pr(>F)
## Condition 2 1675 837.6 9.815 0.000143 ***
## Residuals 87 7425 85.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#with covariate
res1b <- aov(PSS ~ WM + Condition, data = Wk03b)
summary(res1b)
## Df Sum Sq Mean Sq F value Pr(>F)
## WM 1 4373 4373 84.903 1.8e-14 ***
## Condition 2 297 149 2.886 0.0612 .
## Residuals 86 4430 52
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Unadjusted Posthoc
TukeyHSD(res1a) #does not produce correct differences!
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = PSS ~ Condition, data = Wk03b)
##
## $Condition
## diff lwr upr p adj
## 2-1 3.335946 -2.400457 9.072348 0.3522848
## 3-1 10.263012 4.621485 15.904539 0.0001136
## 3-2 6.927066 1.236323 12.617810 0.0128944
#To get the raw (unadjusted) means by group
describeBy(Wk03b$PSS, Wk03b$Condition)
##
## Descriptive statistics by group
## group: 1
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 30 46.29 8.62 46 46.58 9.6 27.95 60.95 33 -0.27 -0.79 1.57
## ------------------------------------------------------------
## group: 2
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 29 49.63 9.07 49.86 49.67 12.45 33.66 65.66 32 -0.03 -1.38
## se
## X1 1.68
## ------------------------------------------------------------
## group: 3
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 31 56.56 9.94 55.37 55.65 7.68 42.37 85.33 42.97 0.84 0.58 1.79
#Adjusted posthoc
# uses the general linear hypothesis function
# looks odd-- uses the glht function. Just
# replace the output object (res1), the
# grouping variable (group)
posthoc1 <- glht(res1b, linfct = mcp(Condition = "Tukey"))
summary(posthoc1)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = PSS ~ WM + Condition, data = Wk03b)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 2 - 1 == 0 -0.807 1.946 -0.415 0.9096
## 3 - 1 == 0 3.550 2.038 1.742 0.1956
## 3 - 2 == 0 4.357 1.885 2.312 0.0594 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
confint(posthoc1)
##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = PSS ~ WM + Condition, data = Wk03b)
##
## Quantile = 2.3839
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 2 - 1 == 0 -0.8070 -5.4469 3.8328
## 3 - 1 == 0 3.5496 -1.3088 8.4080
## 3 - 2 == 0 4.3566 -0.1358 8.8491
#To get the covariate adjusted means by group
effect("Condition", res1b)
##
## Condition effect
## Condition
## 1 2 3
## 49.94049 49.13344 53.49009
#DmSEff
# NOTE: covariate goes first!! NOTE: there is
# no interaction
#without covariate
res2a <- aov(DmSEff ~ Condition, data = Wk03b)
summary(res2a)
## Df Sum Sq Mean Sq F value Pr(>F)
## Condition 2 3075 1537.3 11.76 3.02e-05 ***
## Residuals 87 11376 130.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#with covariate
res2b <- aov(DmSEff ~ WM + Condition, data = Wk03b)
summary(res2b)
## Df Sum Sq Mean Sq F value Pr(>F)
## WM 1 71 71.3 0.554 0.459
## Condition 2 3301 1650.3 12.811 1.35e-05 ***
## Residuals 86 11078 128.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Unadjusted Posthoc
TukeyHSD(res2a) #does not produce correct differences!
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = DmSEff ~ Condition, data = Wk03b)
##
## $Condition
## diff lwr upr p adj
## 2-1 -7.471860 -14.57237 -0.3713516 0.0368012
## 3-1 -14.197936 -21.18101 -7.2148636 0.0000160
## 3-2 -6.726075 -13.77007 0.3179167 0.0644563
#To get the raw (unadjusted) means by group
describeBy(Wk03b$DmSEff, Wk03b$Condition)
##
## Descriptive statistics by group
## group: 1
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 30 56.57 9.91 57.94 57.17 11.11 36.66 74.63 37.98 -0.38 -0.69
## se
## X1 1.81
## ------------------------------------------------------------
## group: 2
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 29 49.1 12.95 50.37 49.7 15.47 14.29 68.73 54.45 -0.5 -0.24 2.4
## ------------------------------------------------------------
## group: 3
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 31 42.37 11.31 40.46 40.87 10.8 28.98 71.3 42.31 0.95 0.23 2.03
#Adjusted posthoc
# uses the general linear hypothesis function
# looks odd-- uses the glht function. Just
# replace the output object (res1), the
# grouping variable (group)
posthoc2 <- glht(res2b, linfct = mcp(Condition = "Tukey"))
summary(posthoc2)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = DmSEff ~ WM + Condition, data = Wk03b)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 2 - 1 == 0 -8.777 3.078 -2.852 0.0149 *
## 3 - 1 == 0 -16.313 3.223 -5.061 <1e-04 ***
## 3 - 2 == 0 -7.536 2.980 -2.529 0.0350 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
confint(posthoc2)
##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = DmSEff ~ WM + Condition, data = Wk03b)
##
## Quantile = 2.3853
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 2 - 1 == 0 -8.7771 -16.1190 -1.4352
## 3 - 1 == 0 -16.3129 -24.0007 -8.6252
## 3 - 2 == 0 -7.5359 -14.6445 -0.4273
#To get the covariate adjusted means by group
effect("Condition", res2b)
##
## Condition effect
## Condition
## 1 2 3
## 57.71756 48.94048 41.40462
#SIMS
# NOTE: covariate goes first!! NOTE: there is
# no interaction
#without covariate
res3a <- aov(SIMS ~ Condition, data = Wk03b)
summary(res3a)
## Df Sum Sq Mean Sq F value Pr(>F)
## Condition 2 224 111.8 1.086 0.342
## Residuals 87 8960 103.0
#with covariate
res3b <- aov(SIMS ~ WM + Condition, data = Wk03b)
summary(res3b)
## Df Sum Sq Mean Sq F value Pr(>F)
## WM 1 10 9.93 0.095 0.758
## Condition 2 224 111.94 1.076 0.346
## Residuals 86 8950 104.07
#Unadjusted Posthoc
TukeyHSD(res3a) #does not produce correct differences!
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = SIMS ~ Condition, data = Wk03b)
##
## $Condition
## diff lwr upr p adj
## 2-1 1.003100 -5.298609 7.304810 0.9237348
## 3-1 3.698509 -2.498975 9.895994 0.3336074
## 3-2 2.695409 -3.556142 8.946960 0.5612504
#To get the raw (unadjusted) means by group
describeBy(Wk03b$SIMS, Wk03b$Condition)
##
## Descriptive statistics by group
## group: 1
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 30 47.59 9.48 45.99 47.49 8.95 30.21 71.82 41.6 0.28 -0.2 1.73
## ------------------------------------------------------------
## group: 2
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 29 48.6 10.81 50.51 49.14 12.08 23.13 65.03 41.9 -0.44 -0.62
## se
## X1 2.01
## ------------------------------------------------------------
## group: 3
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 31 51.29 10.14 51.69 50.94 10.53 33.22 74.79 41.57 0.23 -0.62
## se
## X1 1.82
#Adjusted posthoc
# uses the general linear hypothesis function
# looks odd-- uses the glht function. Just
# replace the output object (res1), the
# grouping variable (group)
posthoc3 <- glht(res3b, linfct = mcp(Condition = "Tukey"))
summary(posthoc3)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = SIMS ~ WM + Condition, data = Wk03b)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 2 - 1 == 0 1.244 2.767 0.450 0.895
## 3 - 1 == 0 4.090 2.897 1.412 0.339
## 3 - 2 == 0 2.845 2.679 1.062 0.540
## (Adjusted p values reported -- single-step method)
confint(posthoc3)
##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = SIMS ~ WM + Condition, data = Wk03b)
##
## Quantile = 2.3842
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 2 - 1 == 0 1.2444 -5.3516 7.8405
## 3 - 1 == 0 4.0896 -2.8172 10.9963
## 3 - 2 == 0 2.8451 -3.5413 9.2316
#To get the covariate adjusted means by group
effect("Condition", res3b)
##
## Condition effect
## Condition
## 1 2 3
## 47.38118 48.62561 51.47075