PSYCH 710 RLab 10
Higher-Order Within-Subjects ANOVA
visual development
In this section we will analyze data from a fictitious experiment
that measured visual contrast sensitivity in six male and six female
children at 30, 36, 42, and 48 months of age. Each child was tested at
each age in two conditions that differed in brightness. The
between-subjects factor (gender) and the two within-subject
factors (age and brightness) are fixed,
whereas subjects is random. The data are in the file
visDevelopment.rda, which contains the data frames
visDevWide and visDev. The data frame
visDev also contains a numeric variable, csf,
which contains the dependent variable.
options(contrasts=c("contr.sum","contr.poly"))
load(url("http://pnb.mcmaster.ca/bennett/psy710/datasets/visDevelopment.rda") )
summary(visDev)
## subj gender age brightness csf
## s1 : 8 female:48 a30:24 high:48 Min. : 68.55
## s2 : 8 male :48 a36:24 low :48 1st Qu.: 91.08
## s3 : 8 a42:24 Median :102.34
## s4 : 8 a48:24 Mean :102.49
## s5 : 8 3rd Qu.:113.00
## s6 : 8 Max. :133.00
## (Other):48
sapply(visDevWide,class)
## subj gender high.30 high.36 high.40 high.48 low.30 low.36 low.42 low.48
## "factor" "factor" "integer" "integer" "integer" "integer" "numeric" "numeric" "numeric" "numeric"Questions
- Use ANOVA to evaluate the effects of
gender,ageandbrightnessoncsf. Conduct two analyses: a split-plot analysis that includes the between-subjects variable and a within-subjects analysis that does not. Where necessary, correct the \(p\) values for deviations from sphericity. Do all of the \(p\) have to be adjusted? Why or why not?
options(contrasts=c("contr.sum","contr.poly"))
# without between-subj factor:
vis.aov.00 <- aov(csf~age*brightness+Error(subj/(age*brightness)),data=visDev)
summary(vis.aov.00) # assumes sphericity
##
## Error: subj
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 11 9875 897.7
##
## Error: subj:age
## Df Sum Sq Mean Sq F value Pr(>F)
## age 3 1186 395.2 3.028 0.0432
## Residuals 33 4307 130.5
##
## Error: subj:brightness
## Df Sum Sq Mean Sq F value Pr(>F)
## brightness 1 2912.1 2912.1 62.79 7.15e-06
## Residuals 11 510.2 46.4
##
## Error: subj:age:brightness
## Df Sum Sq Mean Sq F value Pr(>F)
## age:brightness 3 9.7 3.242 0.117 0.95
## Residuals 33 918.3 27.828
# with between-subj factor:
vis.aov.01 <- aov(csf~gender*age*brightness+Error(subj/(age*brightness)),data=visDev)
summary(vis.aov.01) # assumes sphericity
##
## Error: subj
## Df Sum Sq Mean Sq F value Pr(>F)
## gender 1 61 60.6 0.062 0.809
## Residuals 10 9815 981.5
##
## Error: subj:age
## Df Sum Sq Mean Sq F value Pr(>F)
## age 3 1186 395.2 2.773 0.0586
## gender:age 3 32 10.5 0.074 0.9735
## Residuals 30 4275 142.5
##
## Error: subj:brightness
## Df Sum Sq Mean Sq F value Pr(>F)
## brightness 1 2912.1 2912.1 83.381 3.63e-06
## gender:brightness 1 160.9 160.9 4.607 0.0574
## Residuals 10 349.2 34.9
##
## Error: subj:age:brightness
## Df Sum Sq Mean Sq F value Pr(>F)
## age:brightness 3 9.7 3.242 0.117 0.950
## gender:age:brightness 3 85.8 28.592 1.030 0.393
## Residuals 30 832.5 27.751# following anovas correct p values:
library(afex)
vis.aov_car.00 <- aov_car(csf~age*brightness+Error(subj/(age*brightness)),data=visDev)
summary(vis.aov_car.00)
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 1008450 1 9875.2 11 1123.3191 1.994e-12
## age 1186 3 4306.7 33 3.0284 0.04315
## brightness 2912 1 510.2 11 62.7904 7.152e-06
## age:brightness 10 3 918.3 33 0.1165 0.94979
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## age 0.45067 0.17276
## age:brightness 0.75116 0.73488
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## age 0.65748 0.06981
## age:brightness 0.84113 0.92771
##
## HF eps Pr(>F[HF])
## age 0.8001121 0.05707083
## age:brightness 1.1121152 0.94978986
vis.aov_car.01 <- aov_car(csf~gender*age*brightness+Error(subj/(age*brightness)),data=visDev)
summary(vis.aov_car.01)
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 1008450 1 9814.5 10 1027.5067 2.055e-11
## gender 61 1 9814.5 10 0.0618 0.80876
## age 1186 3 4275.1 30 2.7734 0.05855
## gender:age 32 3 4275.1 30 0.0740 0.97350
## brightness 2912 1 349.2 10 83.3814 3.632e-06
## gender:brightness 161 1 349.2 10 4.6072 0.05741
## age:brightness 10 3 832.5 30 0.1168 0.94953
## gender:age:brightness 86 3 832.5 30 1.0303 0.39321
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## age 0.44446 0.21810
## gender:age 0.44446 0.21810
## age:brightness 0.81124 0.87367
## gender:age:brightness 0.81124 0.87367
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## age 0.65179 0.08801
## gender:age 0.65179 0.92563
## age:brightness 0.88017 0.93371
## gender:age:brightness 0.88017 0.38797
##
## HF eps Pr(>F[HF])
## age 0.8083701 0.07321342
## gender:age 0.8083701 0.95360925
## age:brightness 1.2249865 0.94953271
## gender:age:brightness 1.2249865 0.39320879## We don't have to worry about age because it has 1 degree of freedom, so sphericity MUST be valid. We could use the unadjusted p values for age and age:brightness because the Mauchly test is not significant in either case, but I would still use the adjusted p values because they do a (slightly) better job of controlling Type I error.- Use
lmerin thelmerTestpackage to evaluate the effects ofgender,age, andbrightness(and interactions) with a mixed model. Treatsubjas a random factor and all the other factors as fixed. Evaluate the fixed effects with \(F\) tests that assume sphericity andchi-squaretests that do not assume spericity.
library(lmerTest)
vis.lmer.01 <- lmer(csf~gender*age*brightness + (1|subj),data=visDev)
anova(vis.lmer.01) # assumes sphericity/independence
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## gender 4.81 4.81 1 10 0.0618 0.808756
## age 1185.65 395.22 3 70 5.0698 0.003103
## brightness 2912.06 2912.06 1 70 37.3557 4.965e-08
## gender:age 31.63 10.54 3 70 0.1352 0.938703
## gender:brightness 160.91 160.91 1 70 2.0641 0.155259
## age:brightness 9.73 3.24 3 70 0.0416 0.988604
## gender:age:brightness 85.78 28.59 3 70 0.3668 0.777193
# following is a short-cut... instead of creating multiple nested models
# the command Anova does this for us...
library(car)
Anova(vis.lmer.01,type="III") # chi-sqaure test does not assume sphericity
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: csf
## Chisq Df Pr(>Chisq)
## (Intercept) 1027.5067 1 < 2.2e-16
## gender 0.0618 1 0.803726
## age 15.2095 3 0.001646
## brightness 37.3557 1 9.843e-10
## gender:age 0.4057 3 0.939056
## gender:brightness 2.0641 1 0.150805
## age:brightness 0.1248 3 0.988706
## gender:age:brightness 1.1003 3 0.776996
# Note: You could try to use lme in the nlme package to fit a model that allowed
# the correlations among within-subject variables to vary. However, that model
# would fail to converge on a stable set of estimated parameters because the model
# has a LOT of within-subject levels but not many subjects. In other words, lme doesn't
# have enough data to estimate all of the correlations and so it gives up.
# library(nlme)
# this fails to converge...
# vis.lme.03 <- lme(fixed = csf~gender*age*brightness,
# data=visDev,
# random = ~1|subj,
# correlation=corSymm(form=~1|subj))- Check the residuals of your split-plot model to see if they are distributed normally
shapiro.test(residuals(vis.aov_car.01)) # aov model
##
## Shapiro-Wilk normality test
##
## data: residuals(vis.aov_car.01)
## W = 0.97709, p-value = 0.09089
shapiro.test(residuals(vis.lmer.01)) # lmer model
##
## Shapiro-Wilk normality test
##
## data: residuals(vis.lmer.01)
## W = 0.98895, p-value = 0.6097op <- par(no.readonly = TRUE)
par(mfrow=c(1,2))
qqnorm(residuals(vis.aov_car.01)); qqline(residuals(vis.aov_car.01))
qqnorm(residuals(vis.lmer.01)); qqline(residuals(vis.lmer.01))
Fig 1. QQ Plots
## You might wonder why we are not checking for constant variance across conditions with some sort of Bartlett test. The reason is that the Mauchly test already examined the variance-covariance structure of the residuals and the departure from sphericity was not significant (and in any case we adjusted for the departure from perfect sphericity).- Estimate association strength and effect size for each fixed effect in the split-plot model.
# using aov_car object:
library(effectsize)
omega_squared(vis.aov_car.01)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Omega2 (partial) | 95% CI
## -------------------------------------------------------
## gender | 0.00 | [0.00, 1.00]
## age | 0.05 | [0.00, 1.00]
## gender:age | 0.00 | [0.00, 1.00]
## brightness | 0.20 | [0.00, 1.00]
## gender:brightness | 0.01 | [0.00, 1.00]
## age:brightness | 0.00 | [0.00, 1.00]
## gender:age:brightness | 2.15e-04 | [0.00, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
cohens_f(vis.aov_car.01)
## # Effect Size for ANOVA (Type III)
##
## Parameter | Cohen's f (partial) | 95% CI
## ---------------------------------------------------------
## gender | 0.08 | [0.00, Inf]
## age | 0.53 | [0.00, Inf]
## gender:age | 0.09 | [0.00, Inf]
## brightness | 2.89 | [1.68, Inf]
## gender:brightness | 0.68 | [0.00, Inf]
## age:brightness | 0.11 | [0.00, Inf]
## gender:age:brightness | 0.32 | [0.00, Inf]
##
## - One-sided CIs: upper bound fixed at [Inf].- Evaluate the linear trend of
csfacross age. Is the overall linear trend (ignoring brightness and gender) significant? Does the linear trend differ between genders? Does it differ between brightness conditions?
# create weights
names(visDevWide)
## [1] "subj" "gender" "high.30" "high.36" "high.40" "high.48" "low.30" "low.36" "low.42" "low.48"
ymat <- as.matrix(visDevWide[,c(-1,-2)]) # extract data
ymat[1:2,]
## high.30 high.36 high.40 high.48 low.30 low.36 low.42 low.48
## [1,] 108 96 110 122 89.38860 97.27518 104.7083 117.4466
## [2,] 103 117 127 133 88.54495 116.88224 103.4521 121.1088
# linear trend weights for trend within 1 brightness condition:
lw <- contr.poly(n=4)[,1]
# linear trend weights for OVERALL test:
linWeights <- c(lw,lw) # need to cancat 2 sets (1 for each brightness condition )
ymat[1:2,]
## high.30 high.36 high.40 high.48 low.30 low.36 low.42 low.48
## [1,] 108 96 110 122 89.38860 97.27518 104.7083 117.4466
## [2,] 103 117 127 133 88.54495 116.88224 103.4521 121.1088
linWeights
## [1] -0.6708204 -0.2236068 0.2236068 0.6708204 -0.6708204 -0.2236068 0.2236068 0.6708204# is overall trend significant?
visDevWide$linScores <- ymat %*% linWeights
t.test(visDevWide$linScores,alternative="two.sided") # significant
##
## One Sample t-test
##
## data: visDevWide$linScores
## t = 2.2275, df = 11, p-value = 0.04773
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.1661178 27.7268848
## sample estimates:
## mean of x
## 13.9465
# does trend differ between genders?
t.test(linScores~gender,data=visDevWide) # not significant
##
## Welch Two Sample t-test
##
## data: linScores by gender
## t = -0.14498, df = 8.6774, p-value = 0.888
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
## -31.74933 27.94513
## sample estimates:
## mean in group female mean in group male
## 12.99545 14.89755
# does trend differ between brightness conditions?
# get data for each brightness condition:
yHigh <- ymat[,1:4]
yLow <- ymat[,5:8]
# calculate trend scores:
visDevWide$linHigh <- yHigh %*% lw
visDevWide$linLow <- yLow %*% lw
# compare 2 sets of scores:
with(visDevWide,t.test(linHigh,linLow,paired=TRUE) ) # this is a paired t test
##
## Paired t-test
##
## data: linHigh and linLow
## t = -0.25694, df = 11, p-value = 0.802
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -5.070876 4.010689
## sample estimates:
## mean difference
## -0.5300935# using emmeans:
library(emmeans)
# is overall trend significant?
# compare next 2 results...
visDev.emm.00 <- emmeans(vis.aov_car.00,specs="age") # model >>lacks<< "gender"
# contrast(visDev.emm.00,method=list(lw))
contrast(visDev.emm.00,method="poly")
## contrast estimate SE df t.ratio p.value
## linear 31.19 14.00 11 2.228 0.0477
## quadratic -1.44 3.88 11 -0.371 0.7175
## cubic 2.27 7.41 11 0.306 0.7654
##
## Results are averaged over the levels of: brightness
visDev.emm.01 <- emmeans(vis.aov_car.01,specs="age") # model >>includes<< "gender"
contrast(visDev.emm.01,method="poly")
## contrast estimate SE df t.ratio p.value
## linear 31.19 14.67 10 2.126 0.0594
## quadratic -1.44 4.05 10 -0.356 0.7293
## cubic 2.27 7.69 10 0.295 0.7742
##
## Results are averaged over the levels of: gender, brightness
# note that t tests differ slightly
# reason is that emm objects differ... the >>model<< matters
# because they will have different error terms
# does trend differ between genders?
visDev.emm.02 <- emmeans(vis.aov_car.01,specs="age",by="gender")
contrast(visDev.emm.02,interaction=list(age=list(linTrend=lw),gender=list(gender=c(1,-1))),by=NULL)
## age_custom gender_custom estimate SE df t.ratio p.value
## linTrend gender -0.951 6.56 10 -0.145 0.8876
##
## Results are averaged over the levels of: brightness
# does trend differ between brightnesses?
visDev.emm.03 <- emmeans(vis.aov_car.01,specs="age",by="brightness")
contrast(visDev.emm.03,interaction=list(age=list(linTrend=lw),brightness=list(brightness=c(1,-1))),by=NULL)
## age_custom brightness_custom estimate SE df t.ratio p.value
## linTrend brightness -0.53 2.07 10 -0.256 0.8031
##
## Results are averaged over the levels of: gender