1. Univariate & Factorial ANOVA

1.1. One-way ANOVA

One-way ANOVA (Analysis of Variance) determines whether at least one group mean significantly differs from the others across three or more independent groups.

$H_0:$ All group means are equal ($\mu_1=\mu_2=\cdots=\mu_k$)

$H_1:$ At least one group mean differs from the others

Why not multiple $t$-tests?

Running multiple independent $t$-tests inflates the Family-Wise Error Rate (FWER). Given an alpha level, the probability of making at least one false rejection across $k$ comparisons is:

$FWER=P(\text{making at least one false rejection across $k$ comparisons})= 1 - P(\text{not making a false rejection on any of the $k$ comparisons}) = 1 - (1-\alpha)^k$

For 3 pairwise comparisons at $\alpha=5\%$, FWER rises to $14.3\%$; for 10 comparisons, it exceeds $40\%$. ANOVA solves this problem by testing all groups simultaneously in a single test, keeping the FWER $\alpha = 5\%$.

Partitioning variance

Total variability is broken down into explained and unexplained components: \[SS_{\text{total}}=SS_{\text{between}}+SS_{\text{within}}\] where:

$SS_{\text{between}}$ (Treatment): Variability of group means around the grand mean
$SS_{\text{within}}$ (Error): Variability of individual scores around their respective group means

Assumptions

Before fitting a one-way ANOVA, the following assumptions should be verified:

Assumption	How to check
Independence of observations	Research design (each participant in one group only)
Normality of residuals within each group	Q-Q plots; Shapiro-Wilk test
Homogeneity of variances (homoskedasticity)	Levene’s test; boxplots
No extreme outliers	Boxplots; z-scores

ANOVA is relatively robust to mild violations of normality, especially with balanced designs (equal group sizes) and moderate-to-large sample sizes ($n \ge 15$ per group). It is more sensitive to violations of homogeneity of variances, particularly with unequal group sizes.

What if assumptions are violated?

Non-normal residuals with sufficient sample size ($n \ge 15$ per group): ANOVA is robust; proceed with caution
Non-normal residuals with small samples: Use the non-parametric alternative Kruskal-Wallis test
Unequal variances: Use Welch’s ANOVA; oneway.test(count ~ spray, data = InsectSprays, var.equal = FALSE) which does not assume homoskedasticity

Post-hoc tests

If the ANOVA test is significant, Tukey’s Honestly Significant Difference (HSD) test performs all pairwise comparisons while controlling FWER.

Example

Using InsectSprays, we test if the mean count of insects differs significantly across six different agricultural spray treatments (A through F).

# Load dataset
data(InsectSprays)
head(InsectSprays)

  count spray
1    10     A
2     7     A
3    20     A
4    14     A
5    14     A
6    12     A

tail(InsectSprays)

   count spray
67    13     F
68    10     F
69    26     F
70    26     F
71    24     F
72    13     F

str(InsectSprays)

'data.frame':   72 obs. of  2 variables:
 $ count: num  10 7 20 14 14 12 10 23 17 20 ...
 $ spray: Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...

# Fit one-way ANOVA
fit_oneway <- aov(count ~ spray, data = InsectSprays)
summary(fit_oneway)

            Df Sum Sq Mean Sq F value Pr(>F)    
spray        5   2669   533.8    34.7 <2e-16 ***
Residuals   66   1015    15.4                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Post-hoc pairwise comparisons
TukeyHSD(fit_oneway)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = count ~ spray, data = InsectSprays)

$spray
           diff        lwr       upr     p adj
B-A   0.8333333  -3.866075  5.532742 0.9951810
C-A -12.4166667 -17.116075 -7.717258 0.0000000
D-A  -9.5833333 -14.282742 -4.883925 0.0000014
E-A -11.0000000 -15.699409 -6.300591 0.0000000
F-A   2.1666667  -2.532742  6.866075 0.7542147
C-B -13.2500000 -17.949409 -8.550591 0.0000000
D-B -10.4166667 -15.116075 -5.717258 0.0000002
E-B -11.8333333 -16.532742 -7.133925 0.0000000
F-B   1.3333333  -3.366075  6.032742 0.9603075
D-C   2.8333333  -1.866075  7.532742 0.4920707
E-C   1.4166667  -3.282742  6.116075 0.9488669
F-C  14.5833333   9.883925 19.282742 0.0000000
E-D  -1.4166667  -6.116075  3.282742 0.9488669
F-D  11.7500000   7.050591 16.449409 0.0000000
F-E  13.1666667   8.467258 17.866075 0.0000000

# Assumption Diagnostics

# Shapiro-Wilk test (H0: Data is normally distributed)
shapiro.test(residuals(fit_oneway))


    Shapiro-Wilk normality test

data:  residuals(fit_oneway)
W = 0.96006, p-value = 0.02226

# If p < 0.05, we reject normality.

# Q-Q plot for visual inspection
library(car)
qqPlot(residuals(fit_oneway),
       main = "Q-Q Plot of Residuals",
       xlab = "Theoretical Quantiles",
       ylab = "Standardized Residuals")

[1] 69 70

# Levene's test
leveneTest(count ~ spray, data = InsectSprays)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value   Pr(>F)   
group  5  3.8214 0.004223 **
      66                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# If p < 0.05, variances are significantly different (heteroscedasticity).

1.2. Two-way (Factorial) ANOVA

A two-way ANOVA examines the individual and interactive effects of two categorical factors on a single continuous dependent variable, partitioning variance into three components: \[SS_{\text{total}}=SS_{\text{A}}+SS_{\text{B}}+SS_{\text{AB}}+SS_{\text{within}}\]

Main effect of A: Does the dependent variable differ across levels of factor A, averaging across all levels of factor B?
Main effect of B: Does the dependent variable differ across levels of factor B, averaging across all levels of factor A?
Interaction effect AB: Does the effect of factor A on the dependent variable depend on the level of factor B?

Core rule: Always evaluate the interaction term ($SS_{\text{AB}}$) first. If the interaction is statistically significant, main effects cannot be interpreted independently because the effect of Factor A depends on the specific level of Factor B. In this case, interpret the simple effects (the effect of Factor A at each level of Factor B) rather than main effects.

Example

Using ToothGrowth, we evaluate the main effects and interaction effect of Delivery Method (supp) and Dose Level (dose) on tooth length (len).

# Load dataset
data(ToothGrowth)
head(ToothGrowth)

   len supp dose
1  4.2   VC  0.5
2 11.5   VC  0.5
3  7.3   VC  0.5
4  5.8   VC  0.5
5  6.4   VC  0.5
6 10.0   VC  0.5

tail(ToothGrowth)

    len supp dose
55 24.8   OJ    2
56 30.9   OJ    2
57 26.4   OJ    2
58 27.3   OJ    2
59 29.4   OJ    2
60 23.0   OJ    2

str(ToothGrowth)

'data.frame':   60 obs. of  3 variables:
 $ len : num  4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...
 $ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...
 $ dose: num  0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...

ToothGrowth$dose <- as.factor(ToothGrowth$dose)

# Fit Two-Way ANOVA with interaction
fit_twoway <- aov(len ~ supp * dose, data = ToothGrowth)
summary(fit_twoway)

            Df Sum Sq Mean Sq F value   Pr(>F)    
supp         1  205.4   205.4  15.572 0.000231 ***
dose         2 2426.4  1213.2  92.000  < 2e-16 ***
supp:dose    2  108.3    54.2   4.107 0.021860 *  
Residuals   54  712.1    13.2                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Interaction Plot
interaction.plot(x.factor = ToothGrowth$dose,
                 trace.factor = ToothGrowth$supp,
                 response = ToothGrowth$len,
                 fun = mean,
                 type = "b",
                 col = c("blue", "red"),
                 xlab = "Dose Level",
                 ylab = "Mean Tooth Length",
                 trace.label = "Supplement")

The interaction plot is the most informative visualisation for a factorial design. Parallel lines indicate no interaction; non-parallel (crossing or diverging) lines suggest an interaction.

Because the interaction is significant, we focus on simple effects – the effect of Dose Level within each level of Supplement – rather than interpreting main effects in isolation.

Check assumptions

# Shapiro-Wilk test
shapiro.test(residuals(fit_twoway))


    Shapiro-Wilk normality test

data:  residuals(fit_twoway)
W = 0.98499, p-value = 0.6694

# Levene's test
leveneTest(len ~ supp * dose, data = ToothGrowth)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  5  1.7086 0.1484
      54

Simple effects analysis

When the interaction is significant, we decompose it by examining the effect of one factor at each level of the other:

library(emmeans)
emm_twoway <- emmeans(fit_twoway, ~ dose | supp)  
pairs(emm_twoway, adjust = "bonferroni")

supp = OJ:
 contrast        estimate   SE df t.ratio p.value
 dose0.5 - dose1    -9.47 1.62 54  -5.831 <0.0001
 dose0.5 - dose2   -12.83 1.62 54  -7.900 <0.0001
 dose1 - dose2      -3.36 1.62 54  -2.069  0.1301

supp = VC:
 contrast        estimate   SE df t.ratio p.value
 dose0.5 - dose1    -8.79 1.62 54  -5.413 <0.0001
 dose0.5 - dose2   -18.16 1.62 54 -11.182 <0.0001
 dose1 - dose2      -9.37 1.62 54  -5.770 <0.0001

P value adjustment: bonferroni method for 3 tests

The simple effects analysis reveals that for subjects receiving Vitamin C (VC), increasing dose levels consistently yielded significant increases in tooth growth across all pairwise comparisons (all $p < 0.0001$). However, for subjects receiving Orange Juice (OJ), increasing the dose from 0.5 mg to 1.0 mg and from 0.5 mg to 2.0 mg resulted in significant growth ($p < 0.0001$), but increasing the dose from 1.0 mg to 2.0 mg did not produce a statistically significant increase ($p = 0.1301$).

2. Repeated Measures ANOVA

Repeated measures ANOVA is used when the same participants are measured under three or more conditions. Because measurements come from the same individuals, they are correlated – the repeated measures design removes between-subject variability, increasing statistical power.

The Sphericity assumption

Repeated measures ANOVA requires an additional assumption not needed in independent ANOVA: sphericity – the assumption that variances of the differences between all pairs of repeated conditions are equal.

Sphericity is tested using Mauchly’s test:

$p > 0.05$: Sphericity is not violated – proceed normally
$p \le 0.05$: Sphericity is violated – apply a correction to the degrees of freedom

The most common correction is the Greenhouse-Geisser correction ($\widehat{\epsilon}$), which reduces the degrees of freedom to account for the violation. A more liberal alternative is the Huynh-Feldt correction. When sphericity is severely violated ($\widehat{\epsilon} < 0.75$), the Greenhouse-Geisser correction is preferred.

Example

Using ChickWeight, we track weight changes across time points (Time: 0, 2, 4, 6, 8, 10 days) for individual subjects (Chick) fed Diet 1.

# Filter for complete balanced subset (Diets across early timepoints)
data(ChickWeight)
cw_sub <- subset(ChickWeight, Time %in% c(0, 2, 4, 6, 8, 10))
head(cw_sub)

  weight Time Chick Diet
1     42    0     1    1
2     51    2     1    1
3     59    4     1    1
4     64    6     1    1
5     76    8     1    1
6     93   10     1    1

tail(cw_sub)

    weight Time Chick Diet
567     41    0    50    4
568     54    2    50    4
569     67    4    50    4
570     84    6    50    4
571    105    8    50    4
572    122   10    50    4

str(cw_sub)

Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame':  296 obs. of  4 variables:
 $ weight: num  42 51 59 64 76 93 40 49 58 72 ...
 $ Time  : num  0 2 4 6 8 10 0 2 4 6 ...
 $ Chick : Ord.factor w/ 50 levels "18"<"16"<"15"<..: 15 15 15 15 15 15 17 17 17 17 ...
 $ Diet  : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...

cw_sub$Time <- as.factor(cw_sub$Time)

# Fit Repeated measures ANOVA
# Error(Chick/Time): Specifies that Chick is a random effect and Time is nested within Chick
fit_rm <- aov(weight ~ Time + Error(Chick/Time), data = cw_sub)
summary(fit_rm)


Error: Chick
          Df Sum Sq Mean Sq F value Pr(>F)  
Time       1   2251  2250.7   4.749 0.0343 *
Residuals 48  22747   473.9                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: Chick:Time
           Df Sum Sq Mean Sq F value Pr(>F)    
Time        5 159410   31882   334.2 <2e-16 ***
Residuals 241  22993      95                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The rstatix package provides a more convenient interface that also performs Mauchly’s test and applies corrections automatically:

library(rstatix)

anova_test(weight ~ Time + Error(Chick/Time), data = cw_sub)

ANOVA Table (type III tests)

$ANOVA
  Effect DFn DFd      F         p p<.05   ges
1   Time   5 240 334.02 6.51e-106     * 0.777

$`Mauchly's Test for Sphericity`
  Effect        W       p p<.05
1   Time 1.75e-05 3.1e-98     *

$`Sphericity Corrections`
  Effect   GGe      DF[GG]    p[GG] p[GG]<.05   HFe      DF[HF]    p[HF]
1   Time 0.221 1.11, 53.08 1.86e-25         * 0.223 1.11, 53.42 1.33e-25
  p[HF]<.05
1         *

Post-hoc pairwise comparisons

emm_rm <- emmeans(fit_rm, ~ Time)  
pairs(emm_rm, adjust = "bonferroni")

 contrast       estimate   SE  df t.ratio p.value
 Time0 - Time2     -8.16 1.95 241  -4.177  0.0006
 Time0 - Time4    -18.73 1.97 241  -9.517 <0.0001
 Time0 - Time6    -33.08 1.97 241 -16.805 <0.0001
 Time0 - Time8    -50.02 1.97 241 -25.411 <0.0001
 Time0 - Time10   -66.61 1.97 241 -33.840 <0.0001
 Time2 - Time4    -10.57 1.97 241  -5.371 <0.0001
 Time2 - Time6    -24.92 1.97 241 -12.660 <0.0001
 Time2 - Time8    -41.86 1.97 241 -21.265 <0.0001
 Time2 - Time10   -58.45 1.97 241 -29.694 <0.0001
 Time4 - Time6    -14.35 1.97 241  -7.270 <0.0001
 Time4 - Time8    -31.29 1.97 241 -15.854 <0.0001
 Time4 - Time10   -47.88 1.97 241 -24.262 <0.0001
 Time6 - Time8    -16.94 1.97 241  -8.584 <0.0001
 Time6 - Time10   -33.53 1.97 241 -16.992 <0.0001
 Time8 - Time10   -16.59 1.97 241  -8.408 <0.0001

P value adjustment: bonferroni method for 15 tests

3. MANOVA

MANOVA (Multivariate ANOVA) extends ANOVA to handle two or more correlated dependent variables simultaneously, preventing the inflation of Type I error that occurs when running multiple separate ANOVAs.

Note: Do not use MANOVA simply because we have multiple outcomes – if the outcomes are unrelated, separate ANOVAs with an appropriate correction (e.g., Bonferroni) may be more interpretable.

Multivariate test statistics

MANOVA produces several test statistics, each with slightly different properties:

Statistic	Robustness	Recommendation
Pillai’s trace	Most robust to violations of multivariate normality and homogeneity of covariance matrices	Preferred for most applications
Wilks’ $\lambda$	Traditional default; performs well when assumptions are met	Report alongside Pillai when assumptions met
Hotelling’s trace	Similar to Wilks’ $\lambda$, but more sensitive to severe departures from normality	Rarely preferred

Assumptions of MANOVA

MANOVA extends the ANOVA assumptions to the multivariate case:

Assumption	How to check
Multivariate normality	Q-Q plots per group per variable; Mardia’s test; `rstatix::mshapiro_test()`
Homogeneity of covariance matrices	Box’s M test; `rstatix::box_m()` – Note: very sensitive, treat cautiously
Independence of observations	Research design
No extreme multivariate outliers	Mahalanobis distance
No multicollinearity between DVs	Correlation matrix; if $r > 0.90$, consider combining or dropping variables

Box’s M test

Box’s M test is extremely sensitive – with moderate sample sizes it almost always rejects the null hypothesis of equal covariance matrices even when the violation is minor. Treat a significant Box’s M as a warning rather than an absolute contraindication. If Pillai’s trace is used as the test statistic, MANOVA remains robust to moderate violations.

Example

Using UScereal from the MASS package, we evaluate whether Display Shelf (shelf) simultaneously affects two correlated nutritional variables: Calories (calories) and Carbohydrates (carbo).

library(MASS)
data(UScereal)
head(UScereal)

                          mfr calories   protein      fat   sodium     fibre
100% Bran                   N 212.1212 12.121212 3.030303 393.9394 30.303030
All-Bran                    K 212.1212 12.121212 3.030303 787.8788 27.272727
All-Bran with Extra Fiber   K 100.0000  8.000000 0.000000 280.0000 28.000000
Apple Cinnamon Cheerios     G 146.6667  2.666667 2.666667 240.0000  2.000000
Apple Jacks                 K 110.0000  2.000000 0.000000 125.0000  1.000000
Basic 4                     G 173.3333  4.000000 2.666667 280.0000  2.666667
                             carbo   sugars shelf potassium vitamins
100% Bran                 15.15152 18.18182     3 848.48485 enriched
All-Bran                  21.21212 15.15151     3 969.69697 enriched
All-Bran with Extra Fiber 16.00000  0.00000     3 660.00000 enriched
Apple Cinnamon Cheerios   14.00000 13.33333     1  93.33333 enriched
Apple Jacks               11.00000 14.00000     2  30.00000 enriched
Basic 4                   24.00000 10.66667     3 133.33333 enriched

tail(UScereal)

                    mfr calories  protein      fat   sodium    fibre    carbo
Total Whole Grain     G 100.0000 3.000000 1.000000 200.0000 3.000000 16.00000
Triples               G 146.6667 2.666667 1.333333 333.3333 0.000000 28.00000
Trix                  G 110.0000 1.000000 1.000000 140.0000 0.000000 13.00000
Wheat Chex            R 149.2537 4.477612 1.492537 343.2836 4.477612 25.37313
Wheaties              G 100.0000 3.000000 1.000000 200.0000 3.000000 17.00000
Wheaties Honey Gold   G 146.6667 2.666667 1.333333 266.6667 1.333333 21.33333
                       sugars shelf potassium vitamins
Total Whole Grain    3.000000     3  110.0000     100%
Triples              4.000000     3   80.0000 enriched
Trix                12.000000     2   25.0000 enriched
Wheat Chex           4.477612     1  171.6418 enriched
Wheaties             3.000000     1  110.0000 enriched
Wheaties Honey Gold 10.666667     1   80.0000 enriched

str(UScereal)

'data.frame':   65 obs. of  11 variables:
 $ mfr      : Factor w/ 6 levels "G","K","N","P",..: 3 2 2 1 2 1 6 4 5 1 ...
 $ calories : num  212 212 100 147 110 ...
 $ protein  : num  12.12 12.12 8 2.67 2 ...
 $ fat      : num  3.03 3.03 0 2.67 0 ...
 $ sodium   : num  394 788 280 240 125 ...
 $ fibre    : num  30.3 27.3 28 2 1 ...
 $ carbo    : num  15.2 21.2 16 14 11 ...
 $ sugars   : num  18.2 15.2 0 13.3 14 ...
 $ shelf    : int  3 3 3 1 2 3 1 3 2 1 ...
 $ potassium: num  848.5 969.7 660 93.3 30 ...
 $ vitamins : Factor w/ 3 levels "100%","enriched",..: 2 2 2 2 2 2 2 2 2 2 ...

# Convert shelf factor
UScereal$shelf <- as.factor(UScereal$shelf)

# Response matrix of joint dependent variables
dv_manova <- cbind(UScereal$calories, UScereal$carbo)

# Fit MANOVA
fit_manova <- manova(dv_manova ~ shelf, data = UScereal)

# Test multivariate significance via Pillai's Trace and Wilks' Lambda
summary(fit_manova, test = "Pillai")

          Df  Pillai approx F num Df den Df    Pr(>F)    
shelf      2 0.35303   6.6449      4    124 7.044e-05 ***
Residuals 62                                             
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(fit_manova, test = "Wilks")

          Df   Wilks approx F num Df den Df    Pr(>F)    
shelf      2 0.67746   6.5559      4    122 8.196e-05 ***
Residuals 62                                             
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Follow-up univariate ANOVAs

A significant MANOVA tells us the groups differ on the combination of outcomes. Follow-up univariate ANOVAs identify which individual outcomes drive the multivariate effect:

# Univariate follow-up ANOVA tests
aov_summary <- summary.aov(fit_manova)
aov_summary

 Response 1 :
            Df Sum Sq Mean Sq F value    Pr(>F)    
shelf        2  50435 25217.6  7.8623 0.0009054 ***
Residuals   62 198860  3207.4                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response 2 :
            Df Sum Sq Mean Sq F value   Pr(>F)   
shelf        2  839.9  419.96  6.9437 0.001901 **
Residuals   62 3749.8   60.48                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correction for multiple comparisons across the univariate ANOVAs

pvals <- sapply(aov_summary, function(tab) tab[["Pr(>F)"]][1])
round(p.adjust(pvals, method = "bonferroni"), 4)

 Response 1  Response 2 
     0.0018      0.0038

For each significant univariate ANOVA, run post-hoc comparisons as you would for a standard one-way ANOVA.

How do we run a two-way MANOVA with interactions?

Real designs often have two factors. The same manova() call handles them, we just put both factors on the right side joined with * for an interaction.

Example: We want to see how Gender (Factor 1) and Treatment Type (Factor 2) impact both Blood Pressure and Heart Rate (two distinct dependent variables).

# Two-way MANOVA on mtcars
# Two factors: cylinders and transmission
# Two outcomes: mpg and hp
tw_manova <- manova(cbind(mpg, hp) ~ factor(cyl) * factor(am), data = mtcars)
summary(tw_manova)

                       Df  Pillai approx F num Df den Df    Pr(>F)    
factor(cyl)             2 1.05064  14.3868      4     52 5.656e-08 ***
factor(am)              1 0.46991  11.0808      2     25 0.0003584 ***
factor(cyl):factor(am)  2 0.40390   3.2897      4     52 0.0177291 *  
Residuals              26                                             
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The main effect of cyl is significant: cylinder count moves the (mpg, hp) joint outcome. The main effect of am (transmission) is also significant. The interaction effect is also significant at the $5\%$ level, which means the effect of transmission on the (mpg, hp) bundle depends on cylinder count. In other words, the gap between automatic and manual cars is not the same in 4-cylinder, 6-cylinder, and 8-cylinder cars.

Note: Two-way MANOVA in base R uses Type I (sequential) sums of squares. When our design is unbalanced, the order of factors in the formula affects the main-effect tests. For Type III tests (each effect adjusted for all others), use car::Manova().

4. ANCOVA

ANCOVA (Analysis of Covariance) combines ANOVA with linear regression. It tests whether group means on a dependent variable differ after controlling for one or more continuous covariates.

The goal of ANCOVA is to answer one critical question, for example, if all students had started with the exact same baseline scores, would the intervention (blended, flipped, tutoring) they received still make a significant difference in their final academic performance? By statistically controlling for each student’s prior math and reading ability, ANCOVA isolates the intervention’s true effect on each outcome independently, removing the noise caused by pre-existing differences between students before the instruction even began.

ANCOVA serves two related purposes:

Reducing error variance: By explaining some of the within-group variability via the covariate, ANCOVA increases statistical power
Removing confounding: If groups differ on a variable that also predicts the outcome (a confound), ANCOVA statistically equates the groups on that variable before testing the group effect

Key assumptions of ANCOVA

ANCOVA requires all ANOVA assumptions plus three additional ones:

Assumption	Description	How to check
Linear relationship between covariate and DV	The covariate should correlate linearly with the outcome	Scatterplot; correlation
Homogeneity of regression slopes	The relationship between covariate and DV must be the same in all groups (no group × covariate interaction)	Test the interaction term: `aov(DV ~ group * covariate)`
Independence of covariate and group factor	The covariate should not itself be affected by the group manipulation (especially in experiments)	Design consideration

Homogeneity of regression slopes assumption

If the slopes relating the covariate to the outcome differ across groups, the covariate adjusts scores differently in different groups – making the adjusted means uninterpretable. Always test this assumption by including the group × covariate interaction in the model. If the interaction is significant, ANCOVA assumptions are violated and a different approach (e.g., moderated regression) is needed.

Example

Using UScereal, we test if Manufacturer (mfr) impacts cereal calories (calories), while controlling for continuous Protein content (protein) as a baseline covariate.

# Subset major manufacturers (G = General Mills, K = Kellogg, P = Post)
cereal_sub <- subset(UScereal, mfr %in% c("G", "K", "P"))
cereal_sub$mfr <- factor(cereal_sub$mfr)

# ANCOVA Model (Covariate entered first and evaluated via Type III SS)
fit_ancova <- aov(calories ~ protein + mfr, data = cereal_sub)
car::Anova(fit_ancova, type = "III")

Anova Table (Type III tests)

Response: calories
            Sum Sq Df F value    Pr(>F)    
(Intercept)  93824  1 44.9048 2.113e-08 ***
protein     104967  1 50.2382 5.393e-09 ***
mfr           5726  2  1.3703    0.2638    
Residuals   100291 48                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Test homogeneity of regression slopes

library(ggplot2)

ggplot(cereal_sub, aes(x = protein, y = calories, color = mfr, shape = mfr)) +
  geom_point(alpha = 0.7, size = 2.5) +
  geom_smooth(method = "lm", se = TRUE, alpha = 0.15) +
  labs(title = "Parallel regression lines support homogeneity of slopes", x = "Protein", y = "Calories")

No parallel lines suggest that the relationship between Protein content and cereal calories is not similar across manufactures – violating the homogeneity of slopes assumption.

# Include the interaction to test homogeneity of slopes  
slopes_test <- aov(calories ~ protein * mfr, data = cereal_sub)
car::Anova(slopes_test, type = "III")

Anova Table (Type III tests)

Response: calories
            Sum Sq Df F value    Pr(>F)    
(Intercept)  35158  1  24.760  9.53e-06 ***
protein      20278  1  14.280 0.0004516 ***
mfr          12189  2   4.292 0.0195398 *  
protein:mfr  34972  2  12.315  5.21e-05 ***
Residuals    65319 46                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A significant interaction ($p < 0.05$) confirms that the slopes are not homogeneous – ANCOVA is not appropriate.

5. MANCOVA

MANCOVA extends MANOVA by incorporating one or more continuous covariates into a multivariate model.

MANOVA examines whether groups differ in their mean vectors across multiple dependent variables. MANCOVA examines whether groups differ in their mean vectors after statistically adjusting for the effects of continuous covariates.

Consider a company comparing the performance of its employees across three training programmes. The company measures two dependent variables:

Sales revenue generated per month
Customer satisfaction ratings

However, employees have different levels of prior experience. Experience may influence both sales revenue and customer satisfaction.

A MANOVA could compare the training programmes using the two dependent variables, but it would not explicitly adjust for differences in prior experience. MANCOVA incorporates prior experience as a covariate.

The research question becomes: Do the training programmes differ in their multivariate performance outcomes after controlling for prior work experience?

Why use MANCOVA?

MANCOVA is useful when:

There are multiple related dependent variables
The primary objective is to compare groups
Continuous covariates may be associated with the dependent variables
Researchers want to compare groups after adjusting for those covariates

Example: Comparing training programmes

Suppose a company compares three employee training programmes. The dependent variables are:

Monthly sales revenue
Customer satisfaction
Employee productivity

The covariate is years of previous work experience.

The research question is: Do the training programmes differ in their combined performance outcomes after adjusting for previous work experience?

set.seed(123)

# Sample size
n <- 90

# Training programme
Training <- factor(rep(c("A", "B", "C"), each = 30))

# Continuous covariate
Experience <- runif(n, min = 1, max = 15)

# Generate dependent variables
Sales <- 100 + 5 * Experience + ifelse(Training == "B", 10, 0) +
  ifelse(Training == "C", 18, 0) + rnorm(n, 0, 12)

Satisfaction <- 60 + 1.2 * Experience + ifelse(Training == "B", 4, 0) +
  ifelse(Training == "C", 7, 0) + rnorm(n, 0, 5)

Productivity <- 50 + 2 * Experience + ifelse(Training == "B", 5, 0) +
  ifelse(Training == "C", 9, 0) + rnorm(n, 0, 8)

# Create dataset
training_data <- data.frame(Training, Experience, Sales, Satisfaction, Productivity)
head(training_data)

  Training Experience    Sales Satisfaction Productivity
1        A   5.026085 111.6531     71.68799     60.53017
2        A  12.036272 155.3467     67.14033     68.43577
3        A   6.725677 128.0285     71.77055     57.71361
4        A  13.362244 176.1708     85.58021     83.80169
5        A  14.166542 169.8323     69.78038     70.20834
6        A   1.637791 111.2288     65.47427     68.91793

The dataset contains three dependent variables, one categorical factor, and one continuous covariate.

Fit the MANCOVA model

# Fit multivariate linear model
mancova_model <- manova(cbind(Sales, Satisfaction, Productivity) ~ Training + Experience, data = training_data)

# Display multivariate tests
summary(mancova_model, test = "Pillai")

           Df  Pillai approx F num Df den Df    Pr(>F)    
Training    2 0.42108    7.556      6    170 3.393e-07 ***
Experience  1 0.85673  167.441      3     84 < 2.2e-16 ***
Residuals  86                                             
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The model evaluates the multivariate effects of Training and Experience.

The test for Training addresses whether the training programmes differ in their combined outcomes after adjusting for Experience.

The test for Experience examines whether the covariate is associated with the multivariate response, conditional on Training

As Training is significant, a suitable interpretation is: After adjusting for previous work experience, the training programmes differ in their combined sales, satisfaction, and productivity outcomes.

The specific dependent variables responsible for the multivariate result should be investigated using appropriate follow-up analyses. When follow-up pairwise comparisons are conducted, we should consider appropriate adjustments for multiple comparisons.

Important assumption: Homogeneity of regression slopes

A MANCOVA commonly assumes that the relationship between each dependent variable and the covariate is consistent across groups. This can be investigated by adding the interaction between Training and Experience:

# Fit multivariate linear model
mancova_interaction <- manova(cbind(Sales, Satisfaction, Productivity) ~ Training * Experience, data = training_data)

summary(mancova_interaction, test = "Pillai")

                    Df  Pillai approx F num Df den Df    Pr(>F)    
Training             2 0.42196    7.398      6    166  4.99e-07 ***
Experience           1 0.85688  163.648      3     82 < 2.2e-16 ***
Training:Experience  2 0.03715    0.524      6    166    0.7898    
Residuals           84                                             
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The interaction term evaluates whether the covariate’s multivariate relationship with the outcomes differs across training programmes.

If the interaction is significant, the common-slope MANCOVA model may not be appropriate without further consideration. Group-specific slopes or conditional comparisons may be needed.

The interaction model is a diagnostic extension; it should not automatically replace the main-effects model without considering the research question.

6. Repeated-Measures MANOVA (Profile Analysis)

Repeated-measures MANOVA (Profile analysis) evaluates multivariate profiles across groups by testing three specific hypotheses:

Parallelism: Do the groups follow the same pattern of change?
Level: Do the groups differ in their overall mean response?
Flatness: Does the mean response remain constant across the repeated occasions?

Example: Comparing quarterly financial performance metrics across two market regions

A company wants to compare the quarterly financial performance of two market regions:

North
South

Financial performance is recorded for four quarters.

The research questions are:

Do North and South follow the same pattern of financial performance across quarters?
Do the regions differ in their overall average performance?
Does financial performance change across quarters?

set.seed(123)
region <- factor(rep(c("North", "South"), each = 10))
Q1 <- rnorm(20, 100, 10)
Q2 <- Q1 + rnorm(20, 5, 3)
Q3 <- Q2 + rnorm(20, 8, 4)
Q4 <- Q3 + rnorm(20, 12, 5)

fin_data <- data.frame(region, Q1, Q2, Q3, Q4)
head(fin_data)

  region        Q1        Q2       Q3       Q4
1  North  94.39524  96.19177 101.4129 115.3111
2  North  97.69823 102.04430 109.2126 118.7010
3  North 115.58708 117.50907 120.4475 130.7814
4  North 100.70508 103.51841 120.1942 127.1014
5  North 101.29288 104.41776 117.2496 123.8907
6  North 117.15065 117.09057 120.5981 134.1158

# Summary statistics by region
aggregate(cbind(Q1, Q2, Q3, Q4) ~ region, data = fin_data, FUN = mean)

  region       Q1       Q2       Q3       Q4
1  North 100.7463 104.4726 112.4377 125.0531
2  South 102.0862 108.0524 116.9391 127.1245

Visualize the profiles

# Calculate mean quarterly performance by region
quarter_means <- aggregate(cbind(Q1, Q2, Q3, Q4) ~ region, data = fin_data, FUN = mean)

# Convert the quarterly means into a matrix
profile_matrix <- t(as.matrix(quarter_means[, c("Q1", "Q2", "Q3", "Q4")]))

# Plot the mean profiles
matplot(1:4, profile_matrix, type = "b", pch = c(16, 17), lty = 1, xaxt = "n",
        xlab = "Quarter", ylab = "Mean Financial Performance",
        main = "Quarterly Financial Performance by Region")

axis(1, at = 1:4, labels = c("Q1", "Q2", "Q3", "Q4"))

legend("topleft", legend = quarter_means$region, pch = c(16, 17), lty = 1, bty = "n")

The profile plot helps us assess whether the regions appear to have similar patterns of change, different average levels, or stable performance.

Visual inspection should be followed by formal statistical testing.

Conducting repeated-measures MANOVA using `car` package

# Define response matrix for repeated measures
ymatrix <- as.matrix(fin_data[, c("Q1", "Q2", "Q3", "Q4")])

# Fit multivariate repeated measures model
library(car)
profile_model <- Anova(lm(ymatrix ~ region),
                       idata = data.frame(Quarter = factor(c("Q1", "Q2", "Q3", "Q4"))),
                       idesign = ~Quarter, type = "III")
summary(profile_model, multivariate = TRUE)


Type III Repeated Measures MANOVA Tests:

------------------------------------------
 
Term: (Intercept) 

 Response transformation matrix:
   (Intercept)
Q1           1
Q2           1
Q3           1
Q4           1

Sum of squares and products for the hypothesis:
            (Intercept)
(Intercept)     1959919

Multivariate Tests: (Intercept)
                 Df test stat approx F num Df den Df     Pr(>F)    
Pillai            1   0.98518 1196.322      1     18 < 2.22e-16 ***
Wilks             1   0.01482 1196.322      1     18 < 2.22e-16 ***
Hotelling-Lawley  1  66.46231 1196.322      1     18 < 2.22e-16 ***
Roy               1  66.46231 1196.322      1     18 < 2.22e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------
 
Term: region 

 Response transformation matrix:
   (Intercept)
Q1           1
Q2           1
Q3           1
Q4           1

Sum of squares and products for the hypothesis:
            (Intercept)
(Intercept)    660.3864

Multivariate Tests: region
                 Df test stat  approx F num Df den Df  Pr(>F)
Pillai            1 0.0219037 0.4030956      1     18 0.53348
Wilks             1 0.9780963 0.4030956      1     18 0.53348
Hotelling-Lawley  1 0.0223942 0.4030956      1     18 0.53348
Roy               1 0.0223942 0.4030956      1     18 0.53348

------------------------------------------
 
Term: Quarter 

 Response transformation matrix:
   Quarter1 Quarter2 Quarter3
Q1        1        0        0
Q2        0        1        0
Q3        0        0        1
Q4       -1       -1       -1

Sum of squares and products for the hypothesis:
         Quarter1 Quarter2 Quarter3
Quarter1 5908.244 5002.491 3066.415
Quarter2 5002.491 4235.593 2596.323
Quarter3 3066.415 2596.323 1591.488

Multivariate Tests: Quarter
                 Df test stat approx F num Df den Df     Pr(>F)    
Pillai            1  0.875078 37.35982      3     16 1.8668e-07 ***
Wilks             1  0.124922 37.35982      3     16 1.8668e-07 ***
Hotelling-Lawley  1  7.004966 37.35982      3     16 1.8668e-07 ***
Roy               1  7.004966 37.35982      3     16 1.8668e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------
 
Term: region:Quarter 

 Response transformation matrix:
   Quarter1 Quarter2 Quarter3
Q1        1        0        0
Q2        0        1        0
Q3        0        0        1
Q4       -1       -1       -1

Sum of squares and products for the hypothesis:
          Quarter1  Quarter2  Quarter3
Quarter1  2.674795 -5.516291 -8.886646
Quarter2 -5.516291 11.376373 18.327133
Quarter3 -8.886646 18.327133 29.524682

Multivariate Tests: region:Quarter
                 Df test stat approx F num Df den Df   Pr(>F)  
Pillai            1 0.3518786 2.895578      3     16 0.067495 .
Wilks             1 0.6481214 2.895578      3     16 0.067495 .
Hotelling-Lawley  1 0.5429209 2.895578      3     16 0.067495 .
Roy               1 0.5429209 2.895578      3     16 0.067495 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Univariate Type III Repeated-Measures ANOVA Assuming Sphericity

               Sum Sq num Df Error SS den Df   F value Pr(>F)    
(Intercept)    489980      1   7372.3     18 1196.3215 <2e-16 ***
region            165      1   7372.3     18    0.4031 0.5335    
Quarter          3469      3    665.3     54   93.8547 <2e-16 ***
region:Quarter     31      3    665.3     54    0.8312 0.4826    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Mauchly Tests for Sphericity

               Test statistic    p-value
Quarter               0.19929 6.0022e-05
region:Quarter        0.19929 6.0022e-05


Greenhouse-Geisser and Huynh-Feldt Corrections
 for Departure from Sphericity

                GG eps Pr(>F[GG])    
Quarter        0.54826  1.003e-12 ***
region:Quarter 0.54826     0.4245    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                  HF eps   Pr(>F[HF])
Quarter        0.5961531 1.195465e-13
region:Quarter 0.5961531 4.328032e-01

Interpret the output

The three effects correspond to the profile-analysis hypotheses.

region: Equality of levels
Quarter: Flatness
region:Quarter: Parallelism

At $\alpha=0.05$:

A significant Region × Quarter interaction provides evidence that the regions follow different patterns of change.
A significant Region effect provides evidence of an overall difference in average performance.
A significant Quarter effect provides evidence that average performance changes across quarters.

These effects should be interpreted together rather than as three unrelated tests.

Analysis of Variance (ANOVA) Models & Their Multivariate Extensions in R

Asanka Gunawardana

2026-10-25

1. Univariate & Factorial ANOVA

1.1. One-way ANOVA

Why not multiple \(t\)-tests?

Partitioning variance

Assumptions

What if assumptions are violated?

Post-hoc tests

Example

1.2. Two-way (Factorial) ANOVA

Example

Check assumptions

Simple effects analysis

2. Repeated Measures ANOVA

The Sphericity assumption

Example

Post-hoc pairwise comparisons

3. MANOVA

Multivariate test statistics

Assumptions of MANOVA

Box’s M test

Example

Follow-up univariate ANOVAs

Correction for multiple comparisons across the univariate ANOVAs

How do we run a two-way MANOVA with interactions?

4. ANCOVA

Key assumptions of ANCOVA

Homogeneity of regression slopes assumption

Example

Test homogeneity of regression slopes

5. MANCOVA

Why use MANCOVA?

Example: Comparing training programmes

Fit the MANCOVA model

Important assumption: Homogeneity of regression slopes

6. Repeated-Measures MANOVA (Profile Analysis)

Example: Comparing quarterly financial performance metrics across two market regions

Visualize the profiles

Conducting repeated-measures MANOVA using car package

Interpret the output

Analysis of Variance (ANOVA) Models & Their Multivariate Extensions in `R`

Conducting repeated-measures MANOVA using `car` package