# MANOVA Test Statistics with R

December 8, 2016
By

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Multiple tests of significance can be employed when performing MANOVA. The most well known and widely used MANOVA test statistics are Wilk’s $\Lambda$, Pillai, Lawley-Hotelling, and Roy’s test. Unlike ANOVA in which only one dependent variable is examined, several tests are often utilized in MANOVA due to its multidimensional nature. Each MANOVA test statistic is different and can lead to different conclusions depending on how the data and mean vectors lie, thus understanding how the test statistics are performed and how they differ from one another can be valuable when determining the results of MANOVA.

This post will introduce and explore the three MANOVA test statistics using the rootstock data from the previous MANOVA post. The rootstock data were obtained from the companion FTP site of the book Methods of Multivariate Analysis by Alvin Rencher.

The data is loaded, and the columns are named as before.

root <- read.table('ROOT.DAT', col.names = c('Tree.Number', 'Trunk.Girth.4.Years', 'Ext.Growth.4.Years', 'Trunk.Girth.15.Years', 'Weight.Above.Ground.15.Years'))
root$Tree.Number <- as.factor(root$Tree.Number)


Perform MANOVA on the rootstock data.

root.manova <- manova(cbind(root$Trunk.Girth.4.Years, root$Ext.Growth.4.Years, root$Trunk.Girth.15.Years, root$Weight.Above.Ground.15.Years) ~ Tree.Number,
data = root)
root.summary <- summary(root.manova)
root.summary
##             Df Pillai approx F num Df den Df    Pr(>F)
## Tree.Number  5 1.3055   4.0697     20    168 1.983e-07 ***
## Residuals   42
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


As before, $H_0$ is rejected, and it is concluded there are significant differences in the group mean vectors. The default test statistic used in the summary() call is the Pillai statistic. The same call can be made with summary(root.manova, test='Pillai'). The manova() function provides four multivariate tests by setting the test argument to either Pillai, Wilks, Hotelling-Lawley, or Roy.

The next section examines calculating the Pillai and other MANOVA test statistics and any departures from the above conclusion that may result.

###### Calculating MANOVA Test Statistics

The test statistics rely on the error $E$ and hypothesis $H$ matrices as explored in the previous post, which can be calculated with the following code.

n <- dim(root) / length(unique(root$Tree.Number)) total.means <- colMeans(root[,2:5]) root.group <- split(root[,2:5], root$Tree.Number)
root.means <- sapply(root.group, function(x) {
apply(x, 2, mean)
}, simplify = 'data.frame')

H = matrix(data = 0, nrow = 4, ncol = 4)
for (i in 1:dim(H)) {
for (j in 1:i) {
H[i,j] <- n * sum((root.means[i,] - total.means[i]) * (root.means[j,] - total.means[j]))
H[j,i] <- n * sum((root.means[j,] - total.means[j]) * (root.means[i,] - total.means[i]))
}
}

E = matrix(data = 0, nrow = 4, ncol = 4)
for (i in 1:dim(E)) {
for (j in 1:i) {
b <- c()
for (k in root.group) {
a <- sum((k[,i] - mean(k[,i])) * (k[,j] - mean(k[,j])))
b <- append(b, a)
}
E[i,j] <- sum(b)
E[j,i] <- sum(b)
}
}

H # Hypothesis Matrix
##            [,1]      [,2]      [,3]     [,4]
## [1,] 0.07356042 0.5373852 0.3322646 0.208470
## [2,] 0.53738521 4.1996619 2.3553885 1.637108
## [3,] 0.33226458 2.3553885 6.1139354 3.781044
## [4,] 0.20847000 1.6371084 3.7810438 2.493091

E # Error matrix
##           [,1]      [,2]      [,3]     [,4]
## [1,] 0.3199875  1.696564 0.5540875 0.217140
## [2,] 1.6965637 12.142790 4.3636125 2.110214
## [3,] 0.5540875  4.363613 4.2908125 2.481656
## [4,] 0.2171400  2.110214 2.4816563 1.722525


The same matrices can also be extracted from the previous summary call.

root.summary$SS ##$Tree.Number
##            [,1]      [,2]      [,3]     [,4]
## [1,] 0.07356042 0.5373852 0.3322646 0.208470
## [2,] 0.53738521 4.1996619 2.3553885 1.637108
## [3,] 0.33226458 2.3553885 6.1139354 3.781044
## [4,] 0.20847000 1.6371084 3.7810437 2.493091
##
## Tree.Number
##    1.305472

###### Wilk’s $\Lambda$

Wilk’s $\Lambda$, one of the more commonly used MANOVA test statistics, compares the $E$ matrix to the total $E + H$ matrix. Wilk’s $\Lambda$ is calculated using the determinants of those two matrices.

$$\Lambda = \frac{\left| E \right|}{\left| E + H \right|}$$

Wilk’s $\Lambda$ can also be written in terms of the $E^{-1}H$ eigenvalues $(\lambda_1, \lambda_2, \cdots, \lambda_s)$.

$$\Lambda = \prod^s_{i=1} \frac{1}{1 + \lambda_i}$$
det(E) / det(E + H)
##  0.1540077

prod(1 / (1 + e1h.eigen$values)) ##  0.1540077  summary(root.manova, 'Wilks')$stats[,2]
## Tree.Number
##   0.1540077


An approximate F-statistic can be calculated.

$$F = \frac{1 - \Lambda^{(1/t)}}{\Lambda^{(1/t)}} \frac{df_2}{df_1}$$

Where $df_1$ and $df_2$ degrees of freedom and $t$ is equal to:

$$df_1 = pv_H, \qquad df_2 = wt - \frac{1}{2}(pv_H - 2), \qquad w = v_E + v_H - \frac{1}{2}(p + v_H + 1)$$
$$t = \sqrt{\frac{p^2v^2_H - 4}{p^2 + v^2_H - 5}}$$

k <- length(unique(root$Tree.Number)) p <- length(root[,2:5]) vh <- k - 1 ve <- dim(root) - k t <- sqrt((p^2 * vh^2 - 4) / (p^2 + vh^2 -5)) df1 <- p * vh df2 <- (ve + vh - .5 * (p + vh + 1)) * t - .5 * (p * vh - 2) f <- (1 - (det(E) / det(E + H))^(1/t)) / (det(E) / det(E + H))^(1/t) * df2 / df1 f ##  4.936888  summary(root.manova, 'Wilks')$stats[,3]
## Tree.Number
##    4.936888

###### Lawley-Hotelling Statistic

The Lawley-Hotelling statistic, also known as Hotelling’s generalized $T^2$-statistic, is denoted $U^{(s)}$ and is defined as:

$$U^{(s)} = tr(E^{-1}H) = \sum^s_{i=1} \lambda_i$$
sum(diag(solve(E) %*% H))
##  2.921368

sum(e1h.eigen$values) ##  2.921368  summary(root.manova, 'Hotelling-Lawley')$stats[,2]
## Tree.Number
##    2.921368


The approximate F-statistic in the Lawley-Hotelling setting is defined as:

$$F = \frac{2 (sN + 1)U^{(s)}}{s^2(2m + s + 1)}$$

Where $s$, $m$, and $N$ are defined as:

$$s = min(p, V_h) \qquad m = \frac{1}{2} (\left| V_h - p \right| - 1) \qquad N = \frac{1}{2} (V_E - p - 1)$$
s <- min(vh, p)
m <- .5 * (abs(vh - p) - 1)
N <- .5 * (ve - p - 1)

lawley.hotelling.f <- (2 * (s * N + 1) * sum(diag(solve(E) %*% H))) / (s^2 * (2 * m + s + 1))
lawley.hotelling.f
##  5.477566

summary(root.manova, 'Hotelling-Lawley')$stats[,3] ## Tree.Number ## 5.477566  ###### Roy’s Test Roy’s test statistic is the largest eigenvalue of the matrix $E^{-1}H$. roy.stat <- e1h.eigen$values
roy.stat
##  1.875671

summary(root.manova, 'Roy')$stats[,2] ## Tree.Number ## 1.875671  Roy’s test differs from the previous three in that the F-statistic is found by maximizing the spread of the transformed means. The maximum occurs at the eigenvector corresponding to the largest (first) eigenvalue of $E^{-1}H$. Thus the F-statistic in Roy’s setting can be computed as: $$\frac{k(n - 1)}{k - 1} \lambda_1$$ roy.f <- k * (n - 1) / (k - 1) * e1h.eigen$values
roy.f
##  15.75564

summary(root.manova, 'Roy')\$stats[,3]
## Tree.Number
##    15.75564

?summary.manova

###### Comparison of the MANOVA Test Statistics

According to the documentation in ?summary.manova, “Wilks’ statistic is most popular in the literature, but the default Pillai-Bartlett statistic is recommended by Hand and Taylor (1987).” Although in this example each test rejected $H_0$, there are instances in which one test may accept $H_0$ while the others determine rejection. Rencher recommends not using Roy’s test in any situation unless there is collinearity amongst the dependent variables (Rencher, n.d., pp. 177–177).

###### References

Andrews, D. F., and Herzberg, A. M. (1985), Data, New York: Springer-Verlag.

Rencher, A. (n.d.). Methods of Multivariate Analysis (2nd ed.). Brigham Young University: John Wiley & Sons, Inc.

The post MANOVA Test Statistics with R appeared first on Aaron Schlegel.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.