This script reproduces all of the analysis and graphs for the MANOVA of the Wine 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(Wine, package="candisc")

Fit the MANOVA model

Wine.mod <- lm(as.matrix(Wine[, -1]) ~ Cultivar, data=Wine)
Anova(Wine.mod)
## 
## Type II MANOVA Tests: Pillai test statistic
##          Df test stat approx F num Df den Df Pr(>F)    
## Cultivar  2      1.71     73.2     26    328 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

HE plots

heplot(Wine.mod, fill=TRUE, fill.alpha=.1)

#pairs(Wine.mod, fill=TRUE, fill.alpha=.1)
pairs(Wine.mod, fill=TRUE, fill.alpha=.1, variables=1:5)

univariate F tests

There is no simple way to get a compact summary of the univariate F-tests from a MLM. This code extracts the result of summary() and makes them into a data frame.

ss <- summary(Wine.mod)
UniStats <- as.data.frame(matrix(0, nrow=length(ss), 5))
for (i in 1:length(ss)) {
    UniStats[i,1] <- ss[[i]]$r.squared
    f <- ss[[i]]$fstatistic
    UniStats[i,2:4] <- f
    UniStats[i,5] <- pf(f[1], f[2], f[3], lower.tail=FALSE)
}

rownames(UniStats) <- sub("Response ", "", names(ss))
UniStats$stars <- c(gtools:::stars.pval(UniStats[,5]))
UniStats[,5] <- format.pval(UniStats[,5], eps=0.001)
colnames(UniStats) <- c("R^2", "F", "df1", "df2", "Pr (>F)", "")
UniStats
##                   R^2      F df1 df2 Pr (>F)    
## Alcohol        0.6069 135.08   2 175  <0.001 ***
## MalicAcid      0.2969  36.94   2 175  <0.001 ***
## Ash            0.1321  13.31   2 175  <0.001 ***
## AlcAsh         0.2902  35.77   2 175  <0.001 ***
## Mg             0.1244  12.43   2 175  <0.001 ***
## Phenols        0.5172  93.73   2 175  <0.001 ***
## Flav           0.7278 233.93   2 175  <0.001 ***
## NonFlavPhenols 0.2396  27.58   2 175  <0.001 ***
## Proa           0.2570  30.27   2 175  <0.001 ***
## Color          0.5797 120.66   2 175  <0.001 ***
## Hue            0.5366 101.32   2 175  <0.001 ***
## OD             0.6847 189.97   2 175  <0.001 ***
## Proline        0.7038 207.92   2 175  <0.001 ***

view in canonical space

Wine.can <- candisc(Wine.mod)
Wine.can
## 
## Canonical Discriminant Analysis for Cultivar:
## 
##   CanRsq Eigenvalue Difference Percent Cumulative
## 1  0.901       9.08       4.95    68.7       68.7
## 2  0.805       4.13       4.95    31.3      100.0
## 
## Test of H0: The canonical correlations in the 
## current row and all that follow are zero
## 
##   LR test stat approx F num Df den Df Pr(> F)    
## 1       0.0193      539      4    348  <2e-16 ***
## 2       0.1950      722      1    175  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(Wine.can, ellipse=TRUE, var.lwd=2)

heplot(Wine.can, var.lwd=2)

## Vector scale factor set to  20.73
# need to use extremely small p-value to preserve resolution
heplot(Wine.can, var.lwd=2, alpha=0.00001, fill=TRUE, fill.alpha=.1, var.col="black")

## Vector scale factor set to  10.32

effect ordered pairs.mlm plot

The HE pairs plot is confusing with so many responses, but can be simplified by rearranging the variables in the order derived from the canonical discriminant plot.

struc <- Wine.can$structure
angles <- ifelse( struc[,1] > 0, 
                  atan(struc[,2]/struc[,1]), 
                  atan(struc[,2]/struc[,1]) + pi)
#angles[order(angles)]
colnames(Wine[,-1])[order(angles)]
##  [1] "AlcAsh"         "NonFlavPhenols" "MalicAcid"      "Color"         
##  [5] "Ash"            "Alcohol"        "Mg"             "Proline"       
##  [9] "Phenols"        "Flav"           "Proa"           "OD"            
## [13] "Hue"

This computation is now done by the function varOrder() in the candisc package.

(ord <- varOrder(Wine.mod))
##  [1]  4  8  2 10  3  1  5 13  6  7  9 12 11
varOrder(Wine.mod, names=TRUE)
##  [1] "AlcAsh"         "NonFlavPhenols" "MalicAcid"      "Color"         
##  [5] "Ash"            "Alcohol"        "Mg"             "Proline"       
##  [9] "Phenols"        "Flav"           "Proa"           "OD"            
## [13] "Hue"
ord <- varOrder(Wine.mod)
pairs(Wine.mod, variables=ord, fill=TRUE, fill.alpha=.1, var.cex=1.5)