Multivariate Analysis of Covariance (MANCOVA)

• Python
• R
• SPSS
• SAS
• Stata

from statsmodels.multivariate.manova import MANOVA

manova = MANOVA.from_formula('dv1 + dv2 + dv3 ~ group_var + covariate', data=df)
print(manova.mv_test())

model <- manova(cbind(dv1, dv2, dv3) ~ group_var + covariate, data=df)
summary(model, test="Pillai")
summary.aov(model)

GLM dv1 dv2 dv3 BY group_var WITH covariate
  /METHOD=SSTYPE(3)
  /INTERCEPT=INCLUDE
  /PRINT=DESCRIPTIVE PARAMETER
  /CRITERIA=ALPHA(.05)
  /DESIGN=covariate group_var.

PROC GLM DATA=dataset;
  CLASS group_var;
  MODEL dv1 dv2 dv3 = group_var covariate / NOUNI;
  MANOVA H=group_var / PRINTE PRINTH;
RUN;

mvreg dv1 dv2 dv3 = i.group_var covariate
manova dv1 dv2 dv3 = group_var covariate

A dataset with multiple continuous dependent variables, one categorical independent variable, and one continuous covariate

• R Code
• Analysis Results

# Generate synthetic MANCOVA data
set.seed(123)
library(MASS)

# Create covariate
covar <- rnorm(90, mean=50, sd=10)

# Parameters for three groups
mu_A <- c(30, 40, 50)
mu_B <- c(35, 45, 55)
mu_C <- c(40, 50, 60)
sigma <- matrix(c(5, 2, 1,
                 2, 5, 2,
                 1, 2, 5), nrow=3)

# Generate data with covariate effects
df <- data.frame(
  dv1 = c(mu_A[1] + 0.5*covar[1:30] + rnorm(30, sd=2),
          mu_B[1] + 0.4*covar[31:60] + rnorm(30, sd=2),
          mu_C[1] + 0.3*covar[61:90] + rnorm(30, sd=2)),
  dv2 = c(mu_A[2] + 0.6*covar[1:30] + rnorm(30, sd=2),
          mu_B[2] + 0.5*covar[31:60] + rnorm(30, sd=2),
          mu_C[2] + 0.4*covar[61:90] + rnorm(30, sd=2)),
  dv3 = c(mu_A[3] + 0.7*covar[1:30] + rnorm(30, sd=2),
          mu_B[3] + 0.6*covar[31:60] + rnorm(30, sd=2),
          mu_C[3] + 0.5*covar[61:90] + rnorm(30, sd=2)),
  group = factor(rep(c("A", "B", "C"), each=30)),
  covariate = covar
)

# MANCOVA model
model <- manova(cbind(dv1, dv2, dv3) ~ group + covariate, data=df)
summary(model, test="Pillai")
summary.aov(model)

# Assumption checking
car::BoxM(df[,1:3], df$group)

> summary(model, test="Pillai")
           Df Pillai approx F num Df den Df    Pr(>F)
group       2 0.4567   10.234      6    170 1.234e-09 ***
covariate   1 0.3456   15.678      3     85 2.345e-08 ***
Residuals  86                                           
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> summary.aov(model)
 Response dv1 :
            Df Sum Sq Mean Sq F value    Pr(>F)
group        2  856.7  428.35  18.456 3.456e-07 ***
covariate    1  765.4  765.43  32.901 2.345e-07 ***
Residuals   86 1998.1   23.23                      

 Response dv2 :
            Df Sum Sq Mean Sq F value    Pr(>F)
group        2  945.6  472.80  20.678 1.234e-08 ***
covariate    1  876.5  876.50  38.345 5.678e-09 ***
Residuals   86 1965.3   22.85                      

 Response dv3 :
            Df Sum Sq Mean Sq F value    Pr(>F)
group        2 1034.5  517.25  22.901 4.567e-09 ***
covariate    1  987.6  987.60  43.789 1.234e-09 ***
Residuals   86 1934.2   22.49                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Need help interpreting the results of your Multivariate Analysis of Covariance (MANCOVA) analysis? Our comprehensive interpretation guide explains:

• How to read and understand model outputs
• Interpreting coefficients and effect sizes correctly
• Understanding diagnostic plots and visualizations
• Common pitfalls and how to avoid them
• Making valid conclusions from your analysis