Math scores: HE plot examples

Load packages and the data

library(car)
library(heplots)
library(candisc)
library(MASS)       # for lda()
data(mathscore, package="heplots")

Fit the MLM

math.mod <- lm(cbind(BM, WP) ~ group, data=mathscore)
Anova(math.mod)

## 
## Type II MANOVA Tests: Pillai test statistic
##       Df test stat approx F num Df den Df  Pr(>F)    
## group  1     0.865     28.9      2      9 0.00012 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Compare with univariate ANOVAs

Univariate tests treat the BM and WP responses as independent & don’t take correlation into account.

Anova(lm(BM ~ group, data=mathscore))

## Anova Table (Type II tests)
## 
## Response: BM
##           Sum Sq Df F value Pr(>F)  
## group       1302  1    4.24  0.066 .
## Residuals   3071 10                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(lm(WP ~ group, data=mathscore))

## Anova Table (Type II tests)
## 
## Response: WP
##           Sum Sq Df F value Pr(>F)   
## group       4408  1    10.4  0.009 **
## Residuals   4217 10                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

HE plot

heplot() returns an object containing useful information about the plot. I save this here, to be able to use the plot limits later so the scales align.

hep <- heplot(math.mod, fill=TRUE, cex=2, cex.lab=1.8,
       xlab="Basic math", ylab="Word problems")

Multivariate tests: H & E matrices

math.aov <- Anova(math.mod)

Extract the H & E matrices from the object

(H <- math.aov$SSP)

## $group
##       BM    WP
## BM  1302 -2396
## WP -2396  4408

(E <- math.aov$SSPE)

##      BM   WP
## BM 3071 2808
## WP 2808 4217

Compare with direct calculation

resids <- residuals(math.mod)
crossprod(resids)            # E matrix

##      BM   WP
## BM 3071 2808
## WP 2808 4217

fit <- fitted(math.mod)
ybar <- colMeans(mathscore[,2:3])
n <- nrow(mathscore)
crossprod(fit) - n * outer(ybar, ybar)   # H matrix

##       BM    WP
## BM  1302 -2396
## WP -2396  4408

Visualize data ellipses

Set colors and point symbols

cols <- c("darkgreen", "blue", "red")
pch <- ifelse(mathscore$group==1, 15, 16)

The E matrix is the data ellipse of the residuals – the pooled within group covariance matrix

covEllipses(resids, mathscore$group, 
            pooled=TRUE, cex=2,
            xlab="Basic math", ylab="Word problems",
            fill = c(F, F, TRUE), fill.alpha=0.1,
            col = cols,
            asp=1,
            cex.lab=1.5)
points(resids, pch=pch, col=cols, cex=1.25)

fitted values: H matrix

covEllipses(fit, mathscore$group, 
            pooled=FALSE, cex=2,
            xlab="Basic math", ylab="Word problems",
            fill = c(F, F, TRUE), fill.alpha=0.1,
            col = cols,
            asp=1,
            xlim = hep$xlim, ylim = hep$ylim,
            cex.lab=1.5)
points(jitter(fit), pch=pch, col=cols, cex=1.25)
car::dataEllipse(fit[,1], fit[,2], add=TRUE, level=0.68, col="black", lwd=3)

Discriminant analysis

mod.lda <- MASS::lda(group ~ ., mathscore)

mod.lda

## Call:
## lda(group ~ ., data = mathscore)
## 
## Prior probabilities of groups:
##   1   2 
## 0.5 0.5 
## 
## Group means:
##      BM     WP
## 1 178.3  83.33
## 2 157.5 121.67
## 
## Coefficients of linear discriminants:
##         LD1
## BM -0.08350
## WP  0.07527

names(mod.lda)

##  [1] "prior"   "counts"  "means"   "scaling" "lev"     "svd"     "N"      
##  [8] "call"    "terms"   "xlevels"

Run canonical discriminant analysis

mod.can <- candisc(math.mod)
mod.can

## 
## Canonical Discriminant Analysis for group:
## 
##   CanRsq Eigenvalue Difference Percent Cumulative
## 1  0.865       6.42                100        100
## 
## Test of H0: The canonical correlations in the 
## current row and all that follow are zero
## 
##   LR test stat approx F numDF denDF Pr(> F)
## 1        0.135              2

Show the canonical scores and structure coefficients

mod.can$structure

##       Can1
## BM -0.5867
## WP  0.7686

mod.can$scores

##    group    Can1
## 1      1 -2.7849
## 2      1 -1.8676
## 3      1 -2.7026
## 4      1 -1.3619
## 5      1 -1.7029
## 6      1 -3.4553
## 7      2  1.9783
## 8      2  1.7313
## 9      2  0.5552
## 10     2  2.7310
## 11     2  2.8957
## 12     2  3.9836

Plot the canonical scores and structure coefficients

plot(mod.can, var.lwd=3, points.1d = TRUE, pch=c(15,16), col=cols)