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