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.

Load packages and the data

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")

fit manova model

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"))

Visualizing covariance matrices

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)