**R – Aaron Schlegel**, and kindly contributed to R-bloggers)

Multiple tests of significance can be employed when performing MANOVA. The most well known and widely used MANOVA test statistics are Wilk’s , 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, 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 and hypothesis 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 and defined as:

Where represents the th nonzero eigenvalue of . 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

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

Wilk’s can also be written in terms of the eigenvalues .

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.

Where and degrees of freedom and is equal to:

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 -statistic, is denoted and is defined as:

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:

Where , , and are defined as:

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 .

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 . Thus the F-statistic in Roy’s setting can be computed as:

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 , there are instances in which one test may accept 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.

**leave a comment**for the author, please follow the link and comment on their blog:

**R – Aaron Schlegel**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...