data(Skulls)
library(car) # for Anova
# make shorter labels for epochs
Skulls$epoch <- factor(Skulls$epoch, labels=sub("c","",levels(Skulls$epoch)))
# longer variable labels
vlab <- c("maxBreadth", "basibHeight", "basialLength", "nasalHeight")
# fit manova model
sk.mod <- lm(cbind(mb, bh, bl, nh) ~ epoch, data=Skulls)
Anova(sk.mod)
summary(Anova(sk.mod))
# test trends over epochs
print(linearHypothesis(sk.mod, "epoch.L"), SSP=FALSE) # linear component
print(linearHypothesis(sk.mod, "epoch.Q"), SSP=FALSE) # quadratic component
# typical scatterplots are not very informative
scatterplot(mb ~ bh|epoch, data=Skulls,
ellipse = list(levels=0.68),
smooth=FALSE,
legend = list(coords="topright"),
xlab=vlab[2], ylab=vlab[1])
scatterplot(mb ~ bl|epoch, data=Skulls,
ellipse = list(levels=0.68),
smooth=FALSE,
legend = list(coords="topright"),
xlab=vlab[3], ylab=vlab[1])
# HE plots
heplot(sk.mod,
hypotheses=list(Lin="epoch.L", Quad="epoch.Q"),
xlab=vlab[1], ylab=vlab[2])
pairs(sk.mod,
hypotheses=list(Lin="epoch.L", Quad="epoch.Q"),
var.labels=vlab)
# 3D plot shows that nearly all of hypothesis variation is linear!
if (FALSE) {
heplot3d(sk.mod, hypotheses=list(Lin="epoch.L", Quad="epoch.Q"), col=c("pink", "blue"))
# view in canonical space
if (require(candisc)) {
sk.can <- candisc(sk.mod)
sk.can
heplot(sk.can)
heplot3d(sk.can)
}
}
Run the code above in your browser using DataLab