Set Up Your Project and Load Libraries

knitr::opts_chunk$set(echo = TRUE) 
options(digits = 3)

## Load the libraries we will be using
pacman::p_load(tidyverse, Rfast, heplots, 
               GGally, MVN, car, rstatix)

## Changing the default theme to black white
theme_set(theme_bw())

theme_update(
  axis.title = element_text(size = 10),
  axis.text = element_text(size = 10),
  plot.title = element_text(size = 14, hjust = 0.5),
  plot.subtitle = element_text(size = 10, hjust = 0.5),
  plot.caption = element_text(size = 8)
)
  

# Reading in the data:
turkey <- readxl::read_excel("turkey.xlsx")

# Let's change Trt to a factor and adding the name of the level
# and change Day to a factor as well
turkey <- 
  turkey |> 
  mutate(
    Salt = factor(Trt, 
                  labels = c("Control", "STP3", "STP5",
                             "SASMP3", "SASMP5")),
    Day = as_factor(Day)) |> 
  dplyr::select(-Trt)


# skim creates the basic summary stats for each variable
skimr::skim(turkey)

Data summary
Name	turkey
Number of rows	25
Number of columns	7
_______________________
Column type frequency:
factor	2
numeric	5
________________________
Group variables	None

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
Day	0	1	FALSE	5	1: 5, 2: 5, 3: 5, 4: 5
Salt	0	1	FALSE	5	Con: 5, STP: 5, STP: 5, SAS: 5

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
Cooking_Loss	1	34.20	2.86	26.80	33.00	34.20	36.20	38.50	▂▁▇▅▇
Ph	1	6.61	0.23	6.20	6.49	6.64	6.72	7.34	▃▇▇▁▁
Moist	1	60.22	1.67	56.25	59.57	60.40	61.37	63.37	▂▁▆▇▁
Fat	1	10.00	0.62	9.18	9.52	10.00	10.18	11.72	▆▇▃▁▂
Hex	1	0.77	0.43	0.20	0.40	0.73	1.00	1.58	▇▇▆▂▅

Skim shows us that the experiment is balanced (same n for each treatment)

Preliminary Analysis

For our pre-MANOVA analysis, we want to look for

Outliers
High Correlations
Compare the variables across groups

With summaries
Graphically

We also should calculate the total sample size (N), the number of variables (p), the number of groups (k), and the individual group sizes ($n_1$ to $n_k$)

# Total sample size (N)
N <- nrow(turkey)

# number of treatments (k) 
k <- n_distinct(turkey$Salt)

# Number of response variables (p)
p <- turkey |> select(Cooking_Loss:Hex) |> ncol()

# Since the design is balance, the individual sample size is N/k
n_groups <- 
  turkey |> 
  count(Salt) |> 
  dplyr::select(n) |> 
  unlist()

c("N" = N, "k" = k, "p" = p, n_groups)

##  N  k  p n1 n2 n3 n4 n5 
## 25  5  5  5  5  5  5  5

Preliminary Analysis: Plots

########### Plotting the data #################
# Let's start by pivoting the turkey data into a "long" format
turkey_long <- 
  turkey|> 
  # Creating 1 column for all of the variable values 
  # and another column indicating which variable the value is measuring
  pivot_longer(
    cols = where(is.numeric),
    values_to = "measure",
    names_to = "variable"
  ) |> 
  
  # changing the order of variable groups to match the order of the columns
  mutate(variable = as_factor(variable))


## Side-by-side boxplots per group by variable:
gg_turkey_box <- 
  ggplot(
    data = turkey_long,
    mapping = aes(x = measure)
  ) + 
  geom_boxplot(
    mapping = aes(fill = Salt)
  ) + 
  facet_wrap(
    facets = vars(variable), 
    scales = "free"
  ) +
  theme(legend.position = "top") + 
  labs(fill = NULL)
gg_turkey_box

# We can use ggpairs() to make a scatterplot matrix, 
# But they're small and hard to read
turkey |> 
  dplyr::select(Cooking_Loss:Hex) |> 
  ggpairs(ggplot2::aes(color=turkey$Salt))

## Correlation plot
ggcorr(
  data = turkey,
  low = "red",
  high = "blue",
  mid = "grey90",
  label = T,
  label_round = 2
)

No major cause for concern, but there could be a couple of outliers (but hard to tell since the sample size per group is small)

From the first plot, there does appear to be a difference in salt treatments for all 5 variables, with the control and STP5 having the most differences across the 5 variables.

The correlation plot shows that there doesn’t appear to be linear dependence when ignoring the salt groups

Preliminary Analysis: Comparing Means

# Calculating the overall mean of each variable: Mu
overall_mean <- 
  turkey |> 
  dplyr::select(where(is.numeric)) |> 
  colMeans() 

# Finding the mean vector for each salt type and level: Mu_i
salt_means <- 
  turkey |> 
  summarise(
    .by = Salt,
    across(.cols = where(is.numeric),
           .fns = mean)
  )

# Comparing the overall mean with the individual means:
bind_rows(salt_means,
          overall_mean)

Fitting the MANOVA Model

# manova() needs a formula: cbind(variables) ~ groups
turkey_man <- 
  manova(
    cbind(Cooking_Loss,Ph,Moist,Fat,Hex) ~ Salt, 
    data = turkey
  )

# Use summary(manova_model) to get the results
summary(
  turkey_man, 
  test = "Wilks"
)

##           Df Wilks approx F num Df den Df Pr(>F)  
## Salt       4 0.172      1.9     20     54  0.033 *
## Residuals 20                                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Or alternatively we can use Manova() in the car package, just make sure type = 3 (important for imbalanced designs)

Manova(
  turkey_man, 
  test = "Wilks", 
  type = 3
)

## 
## Type III MANOVA Tests: Wilks test statistic
##             Df test stat approx F num Df den Df Pr(>F)    
## (Intercept)  1    0.0007     4652      5     16 <2e-16 ***
## Salt         4    0.1715        2     20     54  0.033 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Note: The summary() function defaults to Pillai’s test, which is more robust when the sample sizes aren’t balanced or the covariance matrices aren’t equal

summary(turkey_man)

##           Df Pillai approx F num Df den Df Pr(>F)
## Salt       4   1.09     1.42     20     76   0.14
## Residuals 20

Using Wilks’ statistic, we have some evidence of a difference, but using Pillai’s statistic we don’t have evidence of a difference across groups.

We should check the conditions!

Finding what the manova and summary objects offer

# Names function will let you learn what you can pull from the summary of the MANOVA model
names(turkey_man)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "contrasts"     "xlevels"       "call"          "terms"        
## [13] "model"

names(summary(turkey_man))

## [1] "row.names"   "SS"          "Eigenvalues" "stats"

We can get our E and H matrix from summary(manova_model)!

# Getting the treatment matrix (H) and the error matrix (E)
summary(turkey_man)$SS

## $Salt
##              Cooking_Loss     Ph  Moist    Fat    Hex
## Cooking_Loss       106.72 -2.856 -57.93 11.691 13.865
## Ph                  -2.86  0.138   1.40 -0.397 -0.501
## Moist              -57.93  1.397  34.85 -6.565 -6.812
## Fat                 11.69 -0.397  -6.57  1.849  2.018
## Hex                 13.87 -0.501  -6.81  2.018  2.488
## 
## $Residuals
##              Cooking_Loss      Ph  Moist   Fat     Hex
## Cooking_Loss       90.156 -0.2628 -11.80  0.34  3.6296
## Ph                 -0.263  1.1217   1.64 -0.27 -0.0348
## Moist             -11.802  1.6350  32.19 -4.61  1.9382
## Fat                 0.340 -0.2699  -4.61  7.29  1.2521
## Hex                 3.630 -0.0348   1.94  1.25  1.9677

# Getting just the H matrix
H <- summary(turkey_man)$SS$Salt

# Getting just the E matrix
E <- summary(turkey_man)$SS$Residual

The diagonal of the H and E matrices are the sums of squares treatment and sums of squares error of the individual ANOVAs, respectively.

We can calculate our test statistic from the $\mathbf{E}$ and $\mathbf{H}$

det(E)/det(E + H)

## [1] 0.172

summary(turkey_man,
        test = "Wilks")

##           Df Wilks approx F num Df den Df Pr(>F)  
## Salt       4 0.172      1.9     20     54  0.033 *
## Residuals 20                                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Testing the Assumptions

We have 2 major assumptions that need checked:

Values within groups are MVN or that the Residuals are MVN
The covariance matrices are equal for each group

And since we are taking the inverse of a covariance matrix, we need to check for any potential linear dependencies!

Assumption: Residuals are MVN

We can get the residuals using manova_object$residuals. Our example: turkey_man$residuals

#### MVN for the residuals #####
mvn(turkey_man$res, 
    mvnTest="mardia", 
    multivariatePlot = "qq",
    desc = F,
    univariateTest = "SW")

## $multivariateNormality
##              Test        Statistic            p value Result
## 1 Mardia Skewness 49.1992665063625 0.0561892215861621    YES
## 2 Mardia Kurtosis 1.02644074059785  0.304683882760934    YES
## 3             MVN             <NA>               <NA>    YES
## 
## $univariateNormality
##           Test     Variable Statistic   p value Normality
## 1 Shapiro-Wilk Cooking_Loss     0.965    0.5280    YES   
## 2 Shapiro-Wilk      Ph          0.795    0.0002    NO    
## 3 Shapiro-Wilk    Moist         0.967    0.5812    YES   
## 4 Shapiro-Wilk     Fat          0.968    0.5909    YES   
## 5 Shapiro-Wilk     Hex          0.955    0.3178    YES

Small sample, but we still reject that Ph is univariate Normal. There is evidence that the complete data isn’t MVN. May want to remove Ph and rerun the analysis.

From the MV Q-Q plot, we can see that there is one point with a much larger Mahalanobis distance, indicating that there is an outlier that may be cause for concern.

Assumption: Box’s M test for Equal Covariance Matrices

All the tests from above have assumed that the covariance matrices are equal.

We can (try to) test it using Box’s M test: box_m() in rstatix package

?box_m

box_m(
  data = turkey |> 
    dplyr::select(where(is.numeric)),
  group = turkey$Salt
  ) #:(

## Warning in box_m(data = dplyr::select(turkey, where(is.numeric)), group =
## turkey$Salt): there are one or more levels with less observations than
## variables!

What does the warning mean?

Let’s calculate the determinant of the covariance matrix for the control cases

turkey |> 
  filter(Salt=="Control") |> 
  dplyr::select(where(is.numeric)) |> 
  var() |> 
  det()

## [1] 8.19e-19

If the number of cases in a group is less than or equal to the number of variables being tested, then the covariance matrix will always be singular.

Since Box’s M test depends on the determinant of the covariance matrix for each group, if any of the group sizes is less than or equal to p, we won’t be able to use the test!

While we can’t test if the covariance matrices are equal, we can test to see if the variances are all about the same for each salt type using Levene’s test and levene_test() in rstatix package:

turkey_long |> 
  group_by(variable) |> 
  levene_test(measure ~ Salt)

While we can’t check if correlations are the same across salt treatments, we don’t have any evidence that the variances are different!

So are our conditions met? Maybe, we get a little bit of mixed messages:

Mardia test says residuals are MVN but Ph is not Normal
Can’t test if the $\mathbf{\Sigma}_i$ for each salt group are equal, but the test of the variances came back insignificant

Safest (most conservative) to go with Pillai’s statistic and determine no strong evidence of a difference among salt treatments.

Post Hoc Analysis: Finding the differences

If we reject the null hypothesis, all we are claiming is that at least one group mean is different for at least one variable being tested.

Since we have 5 variables, we are essentially conducting 5 ANOVAs at once. Let’s perform each ANOVA individually!

But how can we do it quickly? By using anova_test() in rstatix package

To test to each variable, we need to put in “long” format. Thankfully we already did that earlier when we created turkey_long!

anova_results <- 
  turkey_long |> 
  group_by(variable) |> 
  anova_test(measure ~ Salt) |> 
  data.frame()

anova_results

If we had failed the equal variance test, we could use a more robust test: Kruskal test

kruskal_results <- 
  turkey_long |> 
  group_by(variable) |>  
  kruskal_test(measure ~ Salt) |> 
  data.frame()
kruskal_results

Comparing the results:

left_join(x = anova_results |> 
              dplyr::select(variable, `F`, p),
          y = kruskal_results |>
              dplyr::select(variable, statistic, p),
          by = "variable",
          suffix = c("_aov", "_kru")) |> 
  dplyr::select(variable, `F`, statistic, p_aov, p_kru) |> 
  arrange(p_aov)

Using Kruskal-Wallis test and using Bonferroni’s adjustment (sig = alpha/p), it looks like only Hex would be statistically significant.

But for sake of example, we will use all with a p-value below 0.05

Post Hoc 2: Comparing Significant Variables Across Groups

After we determine which variables are different across groups, we want to find which groups are significantly different from one another.

We can use Tukey’s Honest Significant Differences (HSD) to determine where the differences occur.

Similar to using anova_test(), we want to the data to be in the “long” format to use tukey_hsd() in the rstatix package

### Using tukey_hsd() in the rstatix package
turkey_tukey <- 
  turkey_long |> 
  filter(!(variable %in% c("Ph", "Fat"))) |> # Removing the two non-significant variables:
  group_by(variable) |> 
  tukey_hsd(measure ~ Salt) |> 
  filter(p.adj < 0.05)

turkey_tukey

If we failed the equal variance assumption, we can use a more robust test: Games-Howell test

turkey_long |> 
  filter(!(variable %in% c("Ph", "Fat"))) |> # Removing the two non-significant variables:
  group_by(variable) |>
  games_howell_test(measure ~ Salt) |> 
  filter(p.adj < 0.05)

# Comparing the results to the boxplots from earlier:
gg_turkey_box

NonParametric MANOVA: Randomization Test

If our conditions aren’t met, then we can’t rely on the p-value to give a reliable conclusion to the hypothesis test.

So what do we do instead? Randomization tests!

We’ll calculate the test statistic for a bunch of samples where there truly is no difference between the groups and see how often our test statistic occurred!

# Number of randomization tests to perform:
n_tests <- 10^4

# Getting the observed Wilks and Pillai statistics
Wilk_test_stat <- summary(turkey_man, 
                          test = "Wilks")$stats[1,3]


Pillai_test_stat <- summary(turkey_man, 
                            test = "Pillai")$stats[1,3]

# Creating a matrix to store the test statistics from the
# Randomization tests for both Wilks and Pillai's test
rand_test_stats <- data.frame(Wilks = rep(-1,n_tests),
                              Pillai = rep(-1,n_tests))


# Create a loop to scramble the Salt treatments and
# Find the test statistic of the MANOVA with no real groups

for (i in 1:n_tests){
  # Creating a data set with the 5 variables 
  # and a randomized salt treatment that has no effect
  ran_turkey <- 
    turkey |> 
    mutate(ran_salt = sample(Salt))
  
  # Calculating the manova for the randomization data
  ran_manova <- manova(cbind(Cooking_Loss, Ph, Moist, Fat, Hex) ~ ran_salt, 
                       data = ran_turkey)
  
  # Storing the test statistic from it
  rand_test_stats[i,1] <- summary(ran_manova, test="Wilks")$stats[1,3]
  rand_test_stats[i,2] <- summary(ran_manova)$stats[1,3]
}

rand_test_stats |> 
  pivot_longer(cols = everything(),
               names_to = "type",
               values_to = "statistic") |> 
  
  ggplot(mapping = aes(x = statistic)) +
  
  geom_histogram(mapping = aes(fill = type),
                 bins = 30,
                 color = "black",
                 show.legend = F) +
  
  facet_wrap(facets = ~ type,
             scales = "fixed") + 
  
  geom_vline(data = tibble(statisitic = c(Wilk_test_stat, Pillai_test_stat),
                           type = c("Wilks", "Pillai")),
             mapping = aes(xintercept = statisitic),
             color = "darkred") + 
  
  scale_y_continuous(expand = c(0, 0, 0.05, 0))

# Calculating the proportion of randomization test stats
# Greater than the actual test stat from the experiment
c(
  "Wilks p-value"  = mean(rand_test_stats$Wilks > Wilk_test_stat),
  "Pillai p-value" = mean(rand_test_stats$Pillai > Pillai_test_stat)
)

##  Wilks p-value Pillai p-value 
##         0.0266         0.1121

Chapter 6: One-way Manova

Jacob Martin

STAT 5230