This script reproduces all of the analysis and graphs for the MANCOVA of the Rohwer data in the paper and also includes other analyses not described there. It is set up as an R script that can be “compiled” to HTML, Word, or PDF using knitr::knit()
. This is most convenient within R Studio via the File -> Compile Notebook
option.
library(heplots)
data(Rohwer, package="heplots")
SES
is the group factor, representing Hi
and Lo
socioeconomic status
rohwer.mod <- lm(cbind(SAT, PPVT, Raven) ~ SES + n + s + ns + na + ss,
data=Rohwer)
Anova(rohwer.mod)
##
## Type II MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## SES 1 0.379 12.18 3 60 2.5e-06 ***
## n 1 0.040 0.84 3 60 0.4773
## s 1 0.093 2.04 3 60 0.1173
## ns 1 0.193 4.78 3 60 0.0047 **
## na 1 0.231 6.02 3 60 0.0012 **
## ss 1 0.050 1.05 3 60 0.3770
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pairs(rohwer.mod,
hypotheses=list("Regr" = c("n", "s", "ns", "na", "ss")),
fill=TRUE, fill.alpha=0.1)
rohwer.mod1 <- lm(cbind(SAT, PPVT, Raven) ~ SES * (n + s + ns + na + ss),
data=Rohwer)
Anova(rohwer.mod1)
##
## Type II MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## SES 1 0.391 11.78 3 55 4.5e-06 ***
## n 1 0.079 1.57 3 55 0.20638
## s 1 0.125 2.62 3 55 0.05952 .
## ns 1 0.254 6.25 3 55 0.00100 ***
## na 1 0.307 8.11 3 55 0.00015 ***
## ss 1 0.060 1.17 3 55 0.32813
## SES:n 1 0.072 1.43 3 55 0.24417
## SES:s 1 0.099 2.02 3 55 0.12117
## SES:ns 1 0.118 2.44 3 55 0.07383 .
## SES:na 1 0.148 3.18 3 55 0.03081 *
## SES:ss 1 0.057 1.12 3 55 0.35094
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
linearHypothesis()
takes a vector of terms to test simulataneously. Here we use grep()
to find all interaction terms, which have a “:” in their name.
coefs <- rownames(coef(rohwer.mod1)) # store coefficient names in a vector
print(linearHypothesis(rohwer.mod1, # only test for interaction effects
coefs[grep(":", coefs)]), SSP=FALSE)
##
## Multivariate Tests:
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 5 0.4179 1.845 15 171.0 0.03209 *
## Wilks 5 0.6236 1.894 15 152.2 0.02769 *
## Hotelling-Lawley 5 0.5387 1.927 15 161.0 0.02396 *
## Roy 5 0.3846 4.385 5 57.0 0.00191 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
This also allows us to see possible differences in within-group error variances and covariances
rohwer.ses1 <- lm(cbind(SAT, PPVT, Raven) ~ n + s + ns + na + ss,
data = Rohwer, subset = SES == "Hi")
rohwer.ses2 <- lm(cbind(SAT, PPVT, Raven) ~ n + s + ns + na + ss,
data = Rohwer, subset = SES == "Lo")
For this to work correctly, it is necessary to set the x,y axis limits the explicitly to be the same in both plots
## Low SES students:
heplot(rohwer.ses2, col = c("red", rep("black",5), "blue"),
hypotheses = list("B=0, Low SES" = c("n", "s", "ns", "na", "ss")),
level = 0.5, cex = 1.25,
fill = c(TRUE, FALSE), fill.alpha = 0.05,
xlim = c(-15, 110), ylim = c(40, 110),
xlab = "Student Achievement Test",
ylab = "Peabody Picture Vocabulary Test",
label.pos = c(1, rep(NULL, 5), 1))
## High SES students:
heplot(rohwer.ses1, col = c("red", rep("black", 5), "blue"),
hypotheses = list("B=0, High SES" = c("n", "s", "ns", "na", "ss")),
level = 0.5, cex = 1.25,
add = TRUE, # place both plots on same graphic
error = TRUE, # error ellipse is not drawn by default with add = TRUE
fill = c(TRUE, FALSE), fill.alpha = 0.05,
xlim = c(-15, 110), ylim = c(40, 110))