This script reproduces all of the analysis and graphs for the MANOVA of the Skulls
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)
library(car)
library(candisc)
data(Skulls, package="heplots")
# make shorter labels for epochs
Skulls$epoch <- ordered(Skulls$epoch, labels=sub("c","",levels(Skulls$epoch)))
# make longervariable labels
vlab <- c("maxBreadth", "basibHeight", "basialLength", "nasalHeight")
skulls.mod <- lm(cbind(mb, bh, bl, nh) ~ epoch, data=Skulls)
Manova(skulls.mod)
##
## Type II MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## epoch 4 0.353 3.51 16 580 4.7e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# test trends over epochs
linearHypothesis(skulls.mod, "epoch.L") # linear component
##
## Sum of squares and products for the hypothesis:
## mb bh bl nh
## mb 486.4 -264.85 -606.1 129.88
## bh -264.9 144.21 330.0 -70.72
## bl -606.1 330.03 755.3 -161.84
## nh 129.9 -70.72 -161.8 34.68
##
## Sum of squares and products for error:
## mb bh bl nh
## mb 3061.067 5.333 11.47 291.3
## bh 5.333 3405.267 754.00 412.5
## bl 11.467 754.000 3505.97 164.3
## nh 291.300 412.533 164.33 1472.1
##
## Multivariate Tests:
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.2914 14.6 4 142 5.2e-10 ***
## Wilks 1 0.7086 14.6 4 142 5.2e-10 ***
## Hotelling-Lawley 1 0.4112 14.6 4 142 5.2e-10 ***
## Roy 1 0.4112 14.6 4 142 5.2e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#linearHypothesis(skulls.mod, "epoch.Q") # quadratic component
# test all nonlinear terms
print(linearHypothesis(skulls.mod, c("epoch.Q", "epoch.C", "epoch^4")), SSP=FALSE)
##
## Multivariate Tests:
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 3 0.0682 0.8373 12 432 0.612
## Wilks 3 0.9330 0.8326 12 376 0.617
## Hotelling-Lawley 3 0.0706 0.8279 12 422 0.622
## Roy 3 0.0452 1.6268 4 144 0.171
heplot(skulls.mod, hypotheses=list(Lin="epoch.L", Quad="epoch.Q"),
xlab=vlab[1], ylab=vlab[2])
heplot(skulls.mod,
hypotheses=list(Lin="epoch.L", NonLin=c("epoch.Q", "epoch.C", "epoch^4") ),
xlab=vlab[1], ylab=vlab[2])
pairs(skulls.mod, var.labels=vlab, hypotheses=list(Lin="epoch.L"))
covEllipses(Skulls[,-1], Skulls$epoch,
fill=c(rep(FALSE, 5), TRUE), label.pos=1:4,
xlab=vlab[1], ylab=vlab[2])
covEllipses(Skulls[,-1], Skulls$epoch,
fill=c(rep(FALSE, 5), TRUE), label.pos=1:4,
xlab=vlab[1], ylab=vlab[2], center=TRUE)
All pairwise plots
covEllipses(Skulls[,-1], Skulls$epoch, variables=1:4, vlabels=vlab,
fill=c(rep(FALSE, 5), TRUE), label.pos=1:4,
center=TRUE, fill.alpha=.1)
covEllipses(Skulls[,-1], Skulls$epoch, variables=1:4, vlabels=vlab,
fill=c(rep(FALSE, 5), TRUE), label.pos=1:4)