Multivariate Analysis of Variance (MANOVA)

Multivariate Analysis of Variance (MANOVA) extends the univariate ANOVA framework to situations where multiple dependent variables are measured simultaneously. MANOVA tests whether group means differ across multiple correlated outcome variables, controlling the overall Type I error rate while exploiting the correlation structure among dependent variables. The method provides greater statistical power than conducting separate ANOVAs when dependent variables are correlated and the hypothesis involves simultaneous group differences.

Mathematical Foundation

Test Statistics

MANOVA test statistics are functions of the eigenvalues of $\mathbf{E}^{-1}\mathbf{H}$ or related matrices. Let $s = \min(p, g-1)$ and $\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_s$ be the eigenvalues of $\mathbf{E}^{-1}\mathbf{H}$.

Assumptions and Diagnostics

Post-Hoc Tests and Discriminant Analysis

Python Implementation

import numpy as np
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

# Generate MANOVA data
np.random.seed(42)

# Parameters
n_per_group = 30
n_groups = 3
n_vars = 2

# Group means
means = np.array([
    [10, 8],    # Group 1
    [12, 10],   # Group 2
    [14, 12]    # Group 3
])

# Covariance matrix
Sigma = np.array([
    [4.0, 2.5],
    [2.5, 3.5]
])

# Generate data
data_list = []
for g in range(n_groups):
    Y = np.random.multivariate_normal(means[g], Sigma, n_per_group)
    group_df = pd.DataFrame(Y, columns=['Anxiety', 'Depression'])
    group_df['Group'] = f'Treatment_{g+1}'
    data_list.append(group_df)

data = pd.concat(data_list, ignore_index=True)

# Compute H and E matrices
grand_mean = data[['Anxiety', 'Depression']].mean().values

H = np.zeros((n_vars, n_vars))
E = np.zeros((n_vars, n_vars))

for g in range(n_groups):
    group_data = data[data['Group'] == f'Treatment_{g+1}'][['Anxiety', 'Depression']].values
    n_g = len(group_data)
    group_mean = group_data.mean(axis=0)
    
    # Between-groups matrix
    H += n_g * np.outer(group_mean - grand_mean, group_mean - grand_mean)
    
    # Within-groups matrix
    for i in range(n_g):
        E += np.outer(group_data[i] - group_mean, group_data[i] - group_mean)

print("Hypothesis (Between-groups) matrix H:")
print(np.round(H, 3))
print("\nError (Within-groups) matrix E:")
print(np.round(E, 3))

# Compute eigenvalues of E^{-1}H
E_inv = np.linalg.inv(E)
E_inv_H = E_inv @ H
eigenvalues = np.linalg.eigvals(E_inv_H)
eigenvalues = np.sort(eigenvalues)[::-1]

print(f"\nEigenvalues of E^{{-1}}H: {np.round(eigenvalues, 4)}")

# MANOVA test statistics
s = min(n_vars, n_groups - 1)

# Wilks' Lambda
Lambda = np.prod(1 / (1 + eigenvalues[:s]))
print(f"\nWilks' Lambda: {Lambda:.4f}")

# Pillai's Trace
Pillai = np.sum(eigenvalues[:s] / (1 + eigenvalues[:s]))
print(f"Pillai's Trace: {Pillai:.4f}")

# Hotelling-Lawley Trace
HL_trace = np.sum(eigenvalues[:s])
print(f"Hotelling-Lawley Trace: {HL_trace:.4f}")

# Roy's Largest Root
Roy = eigenvalues[0] / (1 + eigenvalues[0])
print(f"Roy's Largest Root: {Roy:.4f}")

# Rao's approximation for Wilks' Lambda
p = n_vars
g = n_groups
n = len(data)

t = np.sqrt((p**2 * (g-1)**2 - 4) / (p**2 + (g-1)**2 - 5))
df1 = p * (g - 1)
df2 = 4 + (p * (g-1) + 2) * t - (p**2 * (g-1)**2) / (2*t)

F_approx = ((1 - Lambda**(1/t)) / Lambda**(1/t)) * (df2 / df1)
p_value = 1 - stats.f.cdf(F_approx, df1, df2)

print(f"\nRao's F approximation: F({df1:.1f}, {df2:.1f}) = {F_approx:.3f}, p = {p_value:.4f}")

# Box's M test for homogeneity of covariance matrices
def box_m_test(data, group_col, dep_vars):
    """Perform Box's M test for equality of covariance matrices."""
    groups = data[group_col].unique()
    g = len(groups)
    p = len(dep_vars)
    n = len(data)
    
    # Compute pooled and group covariance matrices
    S_pooled = np.zeros((p, p))
    df_total = 0
    
    log_dets = []
    df_groups = []
    
    for group in groups:
        group_data = data[data[group_col] == group][dep_vars].values
        n_g = len(group_data)
        S_g = np.cov(group_data, rowvar=False) * (n_g - 1)
        
        S_pooled += S_g
        df_total += n_g - 1
        
        log_dets.append(np.log(np.linalg.det(S_g / (n_g - 1))))
        df_groups.append(n_g - 1)
    
    S_pooled /= df_total
    log_det_pooled = np.log(np.linalg.det(S_pooled))
    
    # Box's M statistic
    M = (df_total) * log_det_pooled - sum([(df_groups[i]) * log_dets[i] for i in range(g)])
    
    # Correction factor
    sum_inv_df = sum([1/df_groups[i] for i in range(g)])
    C = (2 * p**2 + 3 * p - 1) / (6 * (p + 1) * (g - 1)) * (sum_inv_df - 1/df_total)
    
    # Chi-square approximation
    chi2 = M * (1 - C)
    df_chi2 = p * (p + 1) * (g - 1) / 2
    p_value = 1 - stats.chi2.cdf(chi2, df_chi2)
    
    return M, chi2, df_chi2, p_value

M_stat, chi2_stat, df_m, p_m = box_m_test(data, 'Group', ['Anxiety', 'Depression'])
print(f"\nBox's M test: M = {M_stat:.3f}, Chi2({df_m}) = {chi2_stat:.3f}, p = {p_m:.4f}")

# Discriminant Analysis
lda = LinearDiscriminantAnalysis()
X = data[['Anxiety', 'Depression']].values
y = data['Group'].values
lda.fit(X, y)

# Structure coefficients (correlations with discriminant functions)
disc_scores = lda.transform(X)
structure_corr = np.corrcoef(X.T, disc_scores.T)[:p, p:]
print(f"\nStructure coefficients:\n{np.round(structure_corr, 3)}")

# Visualize discriminant space
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Scatter plot in original space
colors = ['blue', 'red', 'green']
for g, group in enumerate(data['Group'].unique()):
    mask = data['Group'] == group
    axes[0].scatter(data[mask]['Anxiety'], data[mask]['Depression'], 
                   c=colors[g], alpha=0.6, label=group, s=50)

axes[0].set_xlabel('Anxiety')
axes[0].set_ylabel('Depression')
axes[0].set_title('Original Variable Space')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Discriminant scores
for g, group in enumerate(data['Group'].unique()):
    mask = data['Group'] == group
    axes[1].hist(disc_scores[mask, 0], bins=10, alpha=0.5, 
                color=colors[g], label=group, density=True)

axes[1].set_xlabel('Discriminant Score (Function 1)')
axes[1].set_ylabel('Density')
axes[1].set_title('Distribution of Discriminant Scores')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('manova_analysis.png', dpi=150)
plt.show()

# Effect sizes
# Partial eta-squared from Wilks' Lambda
eta2_partial = 1 - Lambda**(1/s)
print(f"\nPartial eta-squared (from Wilks' Lambda): {eta2_partial:.4f}")

# Compare with separate ANOVAs (for illustration)
print("\nSeparate one-way ANOVAs (for comparison):")
for var in ['Anxiety', 'Depression']:
    groups_data = [data[data['Group'] == g][var].values 
                   for g in data['Group'].unique()]
    f_stat, p_val = stats.f_oneway(*groups_data)
    print(f"  {var}: F = {f_stat:.3f}, p = {p_val:.4f}")