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)

Box’s M test

skulls.boxm <- boxM(skulls.mod)
skulls.boxm
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  Y
## Chi-Sq (approx.) = 46, df = 40, p-value = 0.2
summary(skulls.boxm)
## Summary for Box's M-test of Equality of Covariance Matrices
## 
## Chi-Sq:   45.67 
## df:   40 
## p-value: 0.2 
## 
## log of Covariance determinants:
## 4000BC 3300BC 1850BC  200BC  AD150 pooled 
##  11.48  11.03  11.22  10.47  12.15  11.61 
## 
## Eigenvalues:
##   4000BC 3300BC 1850BC  200BC AD150 pooled
## 1 34.823 28.399  26.93 35.906 37.33 29.462
## 2 30.334 22.031  19.07 16.956 28.60 21.441
## 3 18.436 15.596  13.21 13.675 14.80 18.779
## 4  4.952  6.336  11.01  4.229 12.00  9.246
## 
## Statistics based on eigenvalues:
##              4000BC    3300BC    1850BC     200BC     AD150    pooled
## product   96431.296 61827.946 74729.073 35211.658 1.897e+05 1.097e+05
## sum          88.545    72.362    70.228    70.767 9.273e+01 7.893e+01
## precision     3.146     3.305     3.905     2.523 4.703e+00 4.132e+00
## max          34.823    28.399    26.931    35.906 3.733e+01 2.946e+01
plot(skulls.boxm, gplabel="Epoch")

plot other eigenvalue measures

skulls.stats <- summary(skulls.boxm, quiet=TRUE)
eigs <- skulls.stats$eigs
eigstats <- skulls.stats$eigstats

matplot(1:4, log(skulls.stats$eigs), type="b",
    lwd=c(rep(1,5), 3), pch=c(1:5, 'p'),
    lty=c(rep(5,5), 1),
    xlab = "Dimension", ylab="log Eigenvalue")

 op <- par(mfrow=c(2,2), mar=c(5,4,1,1))
 plot(skulls.boxm, which="product", gplabel="Epoch", xlim=c(10,14))
 plot(skulls.boxm, which="sum", gplabel="Epoch")
 plot(skulls.boxm, which="precision", gplabel="Epoch")
 plot(skulls.boxm, which="max", gplabel="Epoch")

 par(op)