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)[1] / 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)[1]) {
  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)[1]) {
  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
## 
## $Residuals
##           [,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.363612 4.2908125 2.481656
## [4,] 0.2171400  2.110214 2.4816562 1.722525

Pillai Test Statistic

The Pillai test statistic is denoted as \(V^{(s)}\) and defined as:

\( V^{(s)} = tr[(E + H0)^{-1} H] = \sum^s_{i=1} \frac{\lambda_i}{1 + \lambda_i} \)

Where \(\lambda_i\) represents the \(i\)th nonzero eigenvalue of \(E^{-1}H\). The Pillai statistic can be computed with either of the below.

e1h.eigen <- eigen(solve(E) %*% H)

sum(diag(solve(E+H) %*% H))
## [1] 1.305472


sum(e1h.eigen$values / (1 + e1h.eigen$values))
## [1] 1.305472


summary(root.manova, 'Pillai')$stats[,2][1]
## 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)
## [1] 0.1540077


prod(1 / (1 + e1h.eigen$values))
## [1] 0.1540077


summary(root.manova, 'Wilks')$stats[,2][1]
## 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)[1] - 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
## [1] 4.936888


summary(root.manova, 'Wilks')$stats[,3][1]
## 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))
## [1] 2.921368


sum(e1h.eigen$values)
## [1] 2.921368


summary(root.manova, 'Hotelling-Lawley')$stats[,2][1]
## 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
## [1] 5.477566


summary(root.manova, 'Hotelling-Lawley')$stats[,3][1]
## 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[1]
roy.stat
## [1] 1.875671


summary(root.manova, 'Roy')$stats[,2][1]
## 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[1]
roy.f
## [1] 15.75564


summary(root.manova, 'Roy')$stats[,3][1]
## 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.