MANOVA Test in R: Multivariate Analysis of Variance

Author

Affiliation

Mtro. Santiago Ríos

What is MANOVA test?

In the situation where there multiple response variables you can test them simultaneously using a multivariate analysis of variance (MANOVA). This article describes how to compute manova in R.

For example, we may conduct an experiment where we give two treatments (A and B) to two groups of mice, and we are interested in the weight and height of mice. In that case, the weight and height of mice are two dependent variables, and our hypothesis is that both together are affected by the difference in treatment. A multivariate analysis of variance could be used to test this hypothesis.

Assumptions of MANOVA

MANOVA can be used in certain conditions:

The dependent variables should be normally distribute within groups. The R function mshapiro.test( )[in the mvnormtest package] can be used to perform the Shapiro-Wilk test for multivariate normality. This is useful in the case of MANOVA, which assumes multivariate normality.
Homogeneity of variances across the range of predictors.
Linearity between all pairs of dependent variables, all pairs of covariates, and all dependent variable-covariate pairs in each cell

Interpretation of MANOVA

If the global multivariate test is significant, we conclude that the corresponding effect (treatment) is significant. In that case, the next question is to determine if the treatment affects only the weight, only the height or both. In other words, we want to identify the specific dependent variables that contributed to the significant global effect.

To answer this question, we can use one-way ANOVA (or univariate ANOVA) to examine separately each dependent variable.

Compute MANOVA in R

Import your data into R

Prepare your data as specified here: [url=/wiki/best-practices-for-preparing-your-data-set-for-r]Best practices for preparing your data set for R[/url]
Save your data in an external .txt tab or .csv files
Import your data into R as follow:

<span># If .txt tab file, use this</span>
<span>my_data</span> <span>&lt;-</span> <span>read.delim</span><span>(</span><span>file.choose</span><span>(</span><span>)</span><span>)</span>
<span># Or, if .csv file, use this</span>
<span>my_data</span> <span>&lt;-</span> <span>read.csv</span><span>(</span><span>file.choose</span><span>(</span><span>)</span><span>)</span>

Here, we’ll use iris data set:

<span># Store the data in the variable my_data</span>
<span>my_data</span> <span>&lt;-</span> <span>iris</span>

Check your data

The R code below display a random sample of our data using the function sample_n()[in dplyr package]. First, install dplyr if you don’t have it:

<span>install.packages</span><span>(</span><span>"dplyr"</span><span>)</span>

<span># Show a random sample</span>
<span>set.seed</span><span>(</span><span>1234</span><span>)</span>
<span>dplyr</span><span>:</span><span>:</span><span>sample_n</span><span>(</span><span>my_data</span>, <span>10</span><span>)</span>

    Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
94           5.0         2.3          3.3         1.0 versicolor
91           5.5         2.6          4.4         1.2 versicolor
93           5.8         2.6          4.0         1.2 versicolor
127          6.2         2.8          4.8         1.8  virginica
150          5.9         3.0          5.1         1.8  virginica
2            4.9         3.0          1.4         0.2     setosa
34           5.5         4.2          1.4         0.2     setosa
96           5.7         3.0          4.2         1.2 versicolor
74           6.1         2.8          4.7         1.2 versicolor
98           6.2         2.9          4.3         1.3 versicolor

Question: We want to know if there is any significant difference, in sepal and petal length, between the different species.

Compute MANOVA test

The function manova() can be used as follow:

<span>sepl</span> <span>&lt;-</span> <span>iris</span><span>$</span><span>Sepal.Length</span>
<span>petl</span> <span>&lt;-</span> <span>iris</span><span>$</span><span>Petal.Length</span>
<span># MANOVA test</span>
<span>res.man</span> <span>&lt;-</span> <span>manova</span><span>(</span><span>cbind</span><span>(</span><span>Sepal.Length</span>, <span>Petal.Length</span><span>)</span> <span>~</span> <span>Species</span>, <span>data</span> <span>=</span> <span>iris</span><span>)</span>
<span>summary</span><span>(</span><span>res.man</span><span>)</span>

           Df Pillai approx F num Df den Df    Pr(&gt;F)    
Species     2 0.9885   71.829      4    294 &lt; 2.2e-16 ***
Residuals 147                                            
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

<span># Look to see which differ</span>
<span>summary.aov</span><span>(</span><span>res.man</span><span>)</span>

 Response Sepal.Length :
             Df Sum Sq Mean Sq F value    Pr(&gt;F)    
Species       2 63.212  31.606  119.26 &lt; 2.2e-16 ***
Residuals   147 38.956   0.265                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 Response Petal.Length :
             Df Sum Sq Mean Sq F value    Pr(&gt;F)    
Species       2 437.10 218.551  1180.2 &lt; 2.2e-16 ***
Residuals   147  27.22   0.185                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

From the output above, it can be seen that the two variables are highly significantly different among Species.

Infos

This analysis has been performed using R software (ver. 3.2.4).

Enjoyed this article? I’d be very grateful if you’d help it spread by emailing it to a friend, or sharing it on Twitter, Facebook or Linked In.

Show me some love with the like buttons below… Thank you and please don’t forget to share and comment below!!

Recommended for You!

Get involved :
Click to follow us on Facebook :
Comment this article by clicking on “Discussion” button (top-right position of this page)

library(ggpubr)
# Data preparation
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Transform `dose` into factor variable
df <- ToothGrowth
df$dose <- as.factor(df$dose)
# Add a random grouping variable
df$group <- factor(rep(c("grp1", "grp2"), 30))
head(df, 3)
#>    len supp dose group
#> 1  4.2   VC  0.5  grp1
#> 2 11.5   VC  0.5  grp2
#> 3  7.3   VC  0.5  grp1


# Two groups by x position
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# Create a box plot
# Add 10% spaces between the p-value labels and the plot border
bxp <- ggpubr::ggboxplot(
  df, x = "dose", y = "len",
  color = "supp", palette = c("#00AFBB", "#E7B800")
) +
 scale_y_continuous(expand = expansion(mult = c(0.05, 0.10)))


# Add p-values onto the box plots
# label can be "p.format"  or "p.adj.format"
bxp + geom_pwc(
  aes(group = supp), tip.length = 0,
  method = "t_test", label = "p.format"
)

# Complex cases
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# 1. Add p-values of OJ vs VC at each dose group
bxp.complex <- bxp +
  geom_pwc(
    aes(group = supp), tip.length = 0,
    method = "t_test", label = "p.adj.format",
    p.adjust.method = "bonferroni", p.adjust.by = "panel"
  )
# 2. Add pairwise comparisons between dose levels
# Nudge up the brackets by 20% of the total height
bxp.complex <- bxp.complex +
  geom_pwc(
    method = "t_test", label = "p.adj.format",
    p.adjust.method = "bonferroni",
    bracket.nudge.y = 0.2
  )
# 3. Display the plot
bxp.complex

# Three groups by x position
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# Simple plots
#_____________________________________

# Box plots with p-values
bxp <- ggboxplot(
  df, x = "supp", y = "len", fill = "dose",
  palette = "npg"
)
bxp +
  geom_pwc(
    aes(group = dose), tip.length = 0,
    method = "t_test", label = "p.adj.format",
    bracket.nudge.y = -0.08
  ) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.1)))

# Bar plots with p-values
bp <- ggbarplot(
  df, x = "supp", y = "len", fill = "dose",
  palette = "npg", add = "mean_sd",
  position = position_dodge(0.8)
)
bp +
  geom_pwc(
    aes(group = dose), tip.length = 0,
    method = "t_test", label = "p.adj.format",
    bracket.nudge.y = -0.08
  ) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.1)))