Multivariate Analysis of Variance (MANOVA) is an ANOVA with two or more continuous dependent (or response) variables.

The one-way MANOVA tests simultaneously statistical differences for multiple response variables by one grouping variables.

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.

Null hypothesis: the multivariate means of weight and height in mice group A and B are equal.

Alternative hypothesis: the multivariate means of weight and height in mice group A and B are not equal when subjected to different treatments.

The procedure of MANOVA can be summarized as:

Loading packages

setwd("D:/Class Materials & Work/Summer 2020 practice/MANOVA")

library(tidyverse) #data manipulation and visualization
library(ggpubr) #for plot
library(rstatix) #for the analysis
library(car) #for the analysis
library(broom) #for data frame coercing
library(GGally) #scatterplot matrix

Data preparation

We’ll use the built-in R dataset iris. First, we will filter for variables of interest.

iris2 <- iris %>%
  select(Sepal.Length, Petal.Length, Species) %>%
  add_column(id = 1:nrow(iris), .before = 1)

str(iris2)
## 'data.frame':    150 obs. of  4 variables:
##  $ id          : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

Visualization

We will create a boxplot shows mean of sepal and petal length grouped by plant species.

ggboxplot(
  iris2, x = "Species", y = c("Sepal.Length", "Petal.Length"), 
  merge = TRUE, palette = "jco")

Summary statistics

Compute summary statistics (mean, SD) by groups (species) for each outcome variable (sepal and petal length):

iris2 %>%
  group_by(Species) %>%
  get_summary_stats(Sepal.Length, Petal.Length, type = "mean_sd")
## # A tibble: 6 x 5
##   Species    variable         n  mean    sd
##   <fct>      <chr>        <dbl> <dbl> <dbl>
## 1 setosa     Petal.Length    50  1.46 0.174
## 2 setosa     Sepal.Length    50  5.01 0.352
## 3 versicolor Petal.Length    50  4.26 0.47 
## 4 versicolor Sepal.Length    50  5.94 0.516
## 5 virginica  Petal.Length    50  5.55 0.552
## 6 virginica  Sepal.Length    50  6.59 0.636

Assumptions and preliminary tests

MANOVA makes the following assumptions about the data:

Check sample size assumption

iris2 %>%
  group_by(Species) %>%
  summarise(N = n())
## # A tibble: 3 x 2
##   Species        N
##   <fct>      <int>
## 1 setosa        50
## 2 versicolor    50
## 3 virginica     50

The result indicates 50 data point per group, which is adequate for the analysis.

Identify univariate outliers

Univariate outliers can be easily identified using box plot methods with rstatix::identify_outliers().

Group the data by Species and then, identify outliers in the Sepal.Length variable:

iris2 %>%
  group_by(Species) %>%
  identify_outliers(Sepal.Length)
## # A tibble: 1 x 6
##   Species      id Sepal.Length Petal.Length is.outlier is.extreme
##   <fct>     <int>        <dbl>        <dbl> <lgl>      <lgl>     
## 1 virginica   107          4.9          4.5 TRUE       FALSE

Group the data by Species and then, identify outliers in the Petal.Length variable:

iris2 %>%
  group_by(Species) %>%
  identify_outliers(Petal.Length)
## # A tibble: 5 x 6
##   Species       id Sepal.Length Petal.Length is.outlier is.extreme
##   <fct>      <int>        <dbl>        <dbl> <lgl>      <lgl>     
## 1 setosa        14          4.3          1.1 TRUE       FALSE     
## 2 setosa        23          4.6          1   TRUE       FALSE     
## 3 setosa        25          4.8          1.9 TRUE       FALSE     
## 4 setosa        45          5.1          1.9 TRUE       FALSE     
## 5 versicolor    99          5.1          3   TRUE       FALSE

There were no univariate extreme outliers in the Sepal.Length and Petal.length variable, as assessed by box plot methods.

Should we have one, it could be because of data entry error, measurement error, or unique data point. Extreme outlier can be included if the researchers believe that results will not be affected. This can be double-checked by running the analysis with- and without extreme outliers.

It is recommended for the researcher to report of decisions about the detected outlier.

Detect multivariate outliers

Multivariate outliers are data points that have an unusual combination of dependent variable value.

In MANOVA, mahalanobis distance is oftentimes used to check the assumption. The distance tells us how far an observation is from the center of the data cluster while considering cluster shape.

The function rstatix::mahalanobis_distance() can be easily used to compute the Mahalanobis distance and to flag multivariate outliers.

This metric needs to be calculated by groups:

# Compute distance by groups and filter outliers
# Use -id to omit the id column in the computation

iris2 %>%
  group_by(Species) %>%
  mahalanobis_distance(-id) %>%
  filter(is.outlier == TRUE) %>%
  as.data.frame()
## [1] id           Sepal.Length Petal.Length mahal.dist   is.outlier  
## <0 rows> (or 0-length row.names)

There were no multivariate outliers in the data, as assessed by Mahalanobis distance (p > 0.001).

If you have multivariate outliers, you could consider running MANOVA before and after removing the outlier to check whether or not their presence alter the results. You should report your final decision.

Check univariate normality assumption

The normality assumption can be checked by computing Shapiro-Wilk test for each outcome variable at each level of the grouping variable. If the data is normally distributed, the p-value should be greater than 0.05 (fail to reject the Null).

iris2 %>%
  group_by(Species) %>%
  shapiro_test(Sepal.Length, Petal.Length) %>%
  arrange(variable)
## # A tibble: 6 x 4
##   Species    variable     statistic      p
##   <fct>      <chr>            <dbl>  <dbl>
## 1 setosa     Petal.Length     0.955 0.0548
## 2 versicolor Petal.Length     0.966 0.158 
## 3 virginica  Petal.Length     0.962 0.110 
## 4 setosa     Sepal.Length     0.978 0.460 
## 5 versicolor Sepal.Length     0.978 0.465 
## 6 virginica  Sepal.Length     0.971 0.258

Sepal.Length and Petal.length were normally distributed for each Species groups, as assessed by Shapiro-Wilk’s test (p > 0.05).

You can also create QQ plot for each group. QQ plot draws the correlation between a given data and the normal distribution.

# QQ plot of Sepal.Length
ggqqplot(iris2, "Sepal.Length", facet.by = "Species",
         ylab = "Sepal Length", ggtheme = theme_bw())

# QQ plot of Petal.Length
ggqqplot(iris2, "Petal.Length", facet.by = "Species",
         ylab = "Petal Length", ggtheme = theme_bw())

All the points fall approximately along the reference line, for each group. So we can assume normality of the data.

Note that, if your sample size is greater than 50, the normal QQ plot is preferred because at larger sample sizes the Shapiro-Wilk test becomes very sensitive even to a minor deviation from normality; Thus, making the test more likely to be significant (normality rejected).

In the situation where the assumptions are not met, you could consider running MANOVA on the data after transforming the outcome variables. You can also perform the test regardless as MANOVA is fairly robust to deviations from normality.

Multivariate normality

iris2 %>%
  select(Sepal.Length, Petal.Length) %>%
  mshapiro_test()
## # A tibble: 1 x 2
##   statistic p.value
##       <dbl>   <dbl>
## 1     0.995   0.855

The test is not significant (p > 0.05), so we can assume multivariate normality.

Identify multicollinearity

Ideally, the correlation between the outcome variables should be moderate. A correlation above 0.9 is an indication of multicollinearity, which is problematic for MANOVA.

In other hand, if the correlation is too low, you should consider running separate one-way ANOVA for each outcome variable.

Compute pairwise Pearson correlation coefficients between the outcome variable with rstatix::cor_test(), or cor_mat() in case of more than two dependent variables.

iris2 %>% cor_test(Sepal.Length, Petal.Length)
## # A tibble: 1 x 8
##   var1         var2           cor statistic        p conf.low conf.high method 
##   <chr>        <chr>        <dbl>     <dbl>    <dbl>    <dbl>     <dbl> <chr>  
## 1 Sepal.Length Petal.Length  0.87      21.6 1.04e-47    0.827     0.906 Pearson

There was no multicollinearity, as assessed by Pearson correlation (r = 0.87, p < 0.0001).

In the situation, where you have multicollinearity, you could consider removing one of the outcome variables that is highly correlated.

Check linearity assumption

The pairwise relationship between the dependent variables should be linear for each group. This can be checked with GGally::ggpairs() by creating scatter plot matrix. In our case, we only have one pair of Sepal.Length and Petal.Length.

pairwise_scatter <- iris2 %>%
  select(Sepal.Length, Petal.Length, Species) %>%
  group_by(Species) %>%
  doo(~ggpairs(.) + theme_bw(), result = "plots")


pairwise_scatter
## # A tibble: 3 x 2
##   Species    plots 
##   <fct>      <list>
## 1 setosa     <gg>  
## 2 versicolor <gg>  
## 3 virginica  <gg>
pairwise_scatter$plots
## [[1]]

## 
## [[2]]

## 
## [[3]]

Relationships between the two dependent variables are linear as indicated by the scatter plot. In the situation, where you detect non-linear relationships, You can either transform or remove the concerned DV, or run the analysis at the risk of losing some statistical power.

Check the homogeneity of covariances assumption

This assumption can be checked with rstatix ::box_m for Box’s M-test.

box_m(iris2[, c("Sepal.Length", "Petal.Length")], iris2$Species)
## # A tibble: 1 x 4
##   statistic  p.value parameter method                                           
##       <dbl>    <dbl>     <dbl> <chr>                                            
## 1      58.4 9.62e-11         6 Box's M-test for Homogeneity of Covariance Matri~

The result shows statistical significance, which indicates the violation of homogeneity of covariances assumption.

Note that, if you have balanced design (i.e., groups with similar sizes), you don’t need to worry too much about the violation of this assumption.

However, having an unbalanced design is problematic. Possible solutions include: 1) transforming the dependent variables; 2) running the test anyway, but using Pillai’s multivariate statistic instead of Wilks’ statistic.

Check the homogneity of variance assumption

For each DV, the one-way MANOVA assumes that there are equal variances between groups. This can be checked through Levene’s test with rstatix::levene_test().

iris2 %>% 
  gather(key = "variable", value = "value", Sepal.Length, Petal.Length) %>%
  group_by(variable) %>%
  levene_test(value ~ Species)
## # A tibble: 2 x 5
##   variable       df1   df2 statistic            p
##   <chr>        <int> <int>     <dbl>        <dbl>
## 1 Petal.Length     2   147     19.5  0.0000000313
## 2 Sepal.Length     2   147      6.35 0.00226

Both DVs have no homogeneity of variances.

If you do not have homogeneity of variances, you can try to transform the outcome (dependent) variable to correct for the unequal variances.

Alternatively, you can continue, but accept a lower level of statistical significance (alpha level) for your MANOVA result.

Additionally, any follow-up univariate ANOVAs will need to be corrected for this violation (i.e., you will need to use different post-hoc tests).

Computation

There are four different types of multivariate statistics that can be used for computing MANOVA, “Pillai”, “Wilks”, “Hotelling-Lawley”, and “Roy”.

The most commonly recommended multivariate statistic to use is Wilks’ Lambda.

However, Pillai’s Trace is more robust and is recommended when you have unbalanced design and also have a statistically significant Box’s M result (as in our example)

Note that, “Pillai” is the default in the car::Manova() function.

iris_model <- lm(cbind(Sepal.Length, Petal.Length) ~ Species, iris2)

Manova(iris_model, test.statistic = "Pillai")
## 
## Type II MANOVA Tests: Pillai test statistic
##         Df test stat approx F num Df den Df    Pr(>F)    
## Species  2    0.9885   71.829      4    294 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

There was a statistically significant difference between the Species on the combined dependent variables (Sepal.Length and Petal.Length), F(4, 294) = 71.829, p < 0.0001.

Therefore, the null hypothesis is rejected. There is at least one pair of difference in the multivariate mean of Sepal.Length and Petal.Length between the provided plant species.

Post-hoc test is recommended.

Post-hoc tests

A statistically significant one-way MANOVA can be followed up by univariate one-way ANOVA examining, separately, each dependent variable. The goal is to identify the specific dependent variables that contributed to the significant overall effect. Then, pairwise comparison can be done to narrow down the variable level that differs from the rest.

The procedure is as follows:

Note that, there are different R function to compute one-way ANOVA depending whether the assumptions are met or not:

We will try all three functions for result comparison:

# Group the data by variable
grouped.data <- iris2 %>%
  gather(key = "variable", value = "value", Sepal.Length, Petal.Length) %>%
  group_by(variable)

# Do welch one way anova test
grouped.data %>% welch_anova_test(value ~ Species)
## # A tibble: 2 x 8
##   variable     .y.       n statistic   DFn   DFd        p method     
## * <chr>        <chr> <int>     <dbl> <dbl> <dbl>    <dbl> <chr>      
## 1 Petal.Length value   150     1828.     2  78.1 2.69e-66 Welch ANOVA
## 2 Sepal.Length value   150      139.     2  92.2 1.51e-28 Welch ANOVA
# or do Kruskal-Wallis test
grouped.data %>% kruskal_test(value ~ Species)
## # A tibble: 2 x 7
##   variable     .y.       n statistic    df        p method        
## * <chr>        <chr> <int>     <dbl> <int>    <dbl> <chr>         
## 1 Petal.Length value   150     130.      2 4.80e-29 Kruskal-Wallis
## 2 Sepal.Length value   150      96.9     2 8.92e-22 Kruskal-Wallis
# or use aov()
grouped.data %>% anova_test(value ~ Species)
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## # A tibble: 2 x 8
##   variable     Effect    DFn   DFd     F        p `p<.05`   ges
##   <chr>        <chr>   <dbl> <dbl> <dbl>    <dbl> <chr>   <dbl>
## 1 Petal.Length Species     2   147 1180. 2.86e-91 *       0.941
## 2 Sepal.Length Species     2   147  119. 1.67e-31 *       0.619

From anova_test(), there is a difference in Sepal.Length (F(2, 147) = 119, p < 0.0001 ) and Petal.Length (F(2, 147) = 1180, p < 0.0001 ) between iris Species.

Note that, as we have two dependent variables, we need to apply Bonferroni multiple testing correction by decreasing the the level we declare statistical significance.

This is done by dividing classic alpha level (0.05) by the number of tests (or dependent variables, for our case, 2). This leads to a significance acceptance criteria of p < 0.025 rather than p < 0.05 because there are two dependent variables.

Compute multiple pairwise comparisons

A statistically significant univariate ANOVA can be followed up by multiple pairwise comparisons to determine which groups are different.

rstatix::tukey_hsd() can be used to compute Tukey post-hoc tests if the homogeneity of variance assumption is met.

If you had violated the assumption of homogeneity of variances, as in our example, you might prefer to run a Games-Howell post-hoc test. It’s also possible to use the function rstatix::pairwise_t_test() with the option pool.sd = FALSE and var.equal = FALSE.

pairwise <- iris2 %>%
  gather(key = "variables", value = "value", Sepal.Length, Petal.Length) %>%
  group_by(variables) %>%
  games_howell_test(value ~ Species) %>%
  select(-estimate, -conf.low, -conf.high) # Remove details

pairwise
## # A tibble: 6 x 6
##   variables    .y.   group1     group2        p.adj p.adj.signif
## * <chr>        <chr> <chr>      <chr>         <dbl> <chr>       
## 1 Petal.Length value setosa     versicolor 1.85e-11 ****        
## 2 Petal.Length value setosa     virginica  1.68e-11 ****        
## 3 Petal.Length value versicolor virginica  4.45e-10 ****        
## 4 Sepal.Length value setosa     versicolor 2.86e-10 ****        
## 5 Sepal.Length value setosa     virginica  0.       ****        
## 6 Sepal.Length value versicolor virginica  5.58e- 7 ****

All pairwise comparisons were significant for each of the outcome variable (Sepal.Length and Petal.Length).

Report

# Visualization: box plots with p-values

pairwise <- pairwise %>% add_xy_position(x = "Species")

test.label <- create_test_label(
  description = "MANOVA", statistic.text = quote(italic("F")),
  statistic = 71.83, p= "<0.0001", parameter = "4,294",
  type = "expression", detailed = TRUE
)
ggboxplot(
  iris2, x = "Species", y = c("Sepal.Length", "Petal.Length"), 
  merge = TRUE, palette = "jco"
) + 
  stat_pvalue_manual(
    pairwise, hide.ns = TRUE, tip.length = 0, 
    step.increase = 0.1, step.group.by = "variables",
    color = "variables"
  ) +
  labs(
    subtitle = test.label,
    caption = get_pwc_label(pairwise, type = "expression")
  )

Reference:https://www.datanovia.com/en/lessons/one-way-manova-in-r/#check-linearity-assumption