Lecture 14: Review

# Placeholder for review image

Lecture 15: ANCOVA Overview

Overview

• Analysis of covariance (ANCOVA):

  • Introduction to ANCOVA

  • The ANCOVA workflow in R

  • When to use ANCOVA

  • Building the ANCOVA model in R

  • Assumptions of ANCOVA

    • Homogeneity of slopes (The critical test!)

  • Interpreting results & “Adjusted Means”

  • Handling violations (heterogeneous slopes)

Introduction to ANCOVA

When to Use ANCOVA

• Common Applications of ANCOVA
  • Increasing statistical power

    • The covariate (e.g., body size, initial temperature) explains some of the “noise” (residual error) in the response variable.

    • By “removing” this variation, we reduce the residual error, making it easier to detect the true effect of our categorical factor.

    • This results in a more powerful test of treatment effects.

  • Adjusting for confounding variables

    • Imagine your treatment groups happen to differ in their average covariate value (e.g., Group A was tested in a warmer room than Group B).

    • ANCOVA allows you to separate the treatment effect from the covariate effect.

  • Testing equality of regression lines
    • Do treatments have the same relationship (slope) with a continuous variable?

ANCOVA Example: Cricket Chirping

ANCOVA Linear Model: Conceptual Framework

• The ANCOVA Model
  • Key concept: Dfference between “unadjusted” group means and “adjusted” means
• In this visualization:
  • Group Means (asterisks): The raw average X and Y for each group.
  • Group A and Group B have different mean X values. Comparing directly misleading
  • Adjusted Means (triangles): ANCOVA statistically “moves” each group mean along regression line to single, common X value (the overall mean X).
  • ANCOVA compares heights (the Y values) of adjusted means.
  • Adjustment “levels the playing field,” estimating each group’s mean if all groups had same average covariate value

ANCOVA Model Visualization

Mathematical Model for ANCOVA

• Single-factor ANCOVA with
  • factor A (p levels, i = 1 to p), a continuous covariate (x), and response variable (y): \(Y_{ij} = \mu + \alpha_i + \beta(X_{ij} - \bar{X}) + \varepsilon_{ij}\)
• Where:
  • \(Y_{ij}\) = response value for observation j in level i of factor A
  • \(\mu\) = overall mean
  • \(\alpha_i\) = effect of level i of factor A (what we test in ANOVA)
  • \(\beta\) = common regression slope relating Y to X (what we test in Regression)
  • \(X_{ij}\) = covariate value for observation j in level i of factor A
  • \(\bar{X}\) = mean value of covariate
  • \(\varepsilon_{ij}\) = error term

ANCOVA Parameters Interpretation

Building the ANCOVA Model in R

• Remember lm(y ~ predictor).
  • ANCOVA combines the predictors from ANOVA and regression.
• Simple Linear Regression:
  • Predict y using a continuous variable x.
    • model_reg <- lm(y ~ x, data = mydata)
• One-Way ANOVA:
  • We want to compare y among groups in a categorical factor A.
    • model_aov <- lm(y ~ A, data = mydata)
• ANCOVA (Additive Model):
  • Compare y among groups in A, after accounting for effect of x.
  • Add terms together
    • model_ancova <- lm(y ~ A + x, data = mydata)
  • Model assumes slope of y ~ x relationship is same for all groups in A.

Building the ANCOVA Model in R

• ANCOVA (Interaction Model):
  • What if the slopes are different for each group?
  • Test for interaction between factor and covariate.
    • model_int <- lm(y ~ A * x, data = mydata)
  • Shorthand for: lm(y ~ A + x + A:x, data = mydata)
  • The A:x term tests if the slopes are different.

ANCOVA in R: The 3-Step Workflow

• Step 1: Test for Homogeneity of Slopes - Always Start Here

  • Fit interaction model first - tests key assumption that all groups share same slope.

  # The asterisk * tests for both main effects AND the interaction
  model_int <- lm(y ~ A * x, data = mydata)

• Step 2: Examine the Interaction Term

  • ANOVA table for model_int and find the p-value for interaction term (A:x).

      Anova(model_int, type = "III")

  • If p > 0.05 (e.g., p = 0.4) - interaction IS NOT significant - slopes are homogeneous (parallel).
  • If p < 0.05 (e.g., p = 0.02): - interaction IS significant - slopes heterogeneous (NOT parallel).
  • ** standard ANCOVA model (y ~ A + x) inappropriate**
  • Analysis stops here and interpret interaction (often most interesting result!)

ANCOVA in R: The 3-Step Workflow

• Step 3: Run the ANCOVA (Additive Model)

  • If interaction not significant, run the simpler additive model the final ANCOVA model.

  # We use the + sign, which means "additive"
  model_ancova <- lm(y ~ A + x, data = mydata)
  # Get the final ANCOVA table
  Anova(model_ancova, type = "III")
  # Get the adjusted means
  library(emmeans)
  emmeans(model_ancova, "A")

Analysis of Variance for ANCOVA: Partitioning

ANOVA Table for ANCOVA

The ANOVA table for a single-factor ANCOVA has these components: :::{.panel column=“screen”}

# Create a tibble with ANOVA table information
anova_table <- tibble(
  Source = c("Factor A (adjusted)", "Covariate", "Residual", "Total"),
  df = c("(p-1)", "1", "n-p-1", "n-1"),
  `Sum of Squares` = c("SS_A(adj)", "SS_Covariate", "SS_Residual", "SS_Total"),
  `Mean Square` = c("MS_A(adj) = SS_A(adj)/(p-1)", "MS_Covariate = SS_Covariate/1", 
                   "MS_Residual = SS_Residual/(n-p-1)", ""),
  `F-ratio` = c("MS_A(adj)/MS_Residual", "MS_Covariate/MS_Residual", "", ""),
  `Expected MS` = c("σ² + n∑α²/(p-1)", "σ² + β²∑(X-X̄)²", "σ²", "")
)

# Create
anova_table
## # A tibble: 4 × 6
##   Source            df    `Sum of Squares` `Mean Square` `F-ratio` `Expected MS`
##   <chr>             <chr> <chr>            <chr>         <chr>     <chr>        
## 1 Factor A (adjust… (p-1) SS_A(adj)        "MS_A(adj) =… "MS_A(ad… "σ² + n∑α²/(…
## 2 Covariate         1     SS_Covariate     "MS_Covariat… "MS_Cova… "σ² + β²∑(X-…
## 3 Residual          n-p-1 SS_Residual      "MS_Residual… ""        "σ²"         
## 4 Total             n-1   SS_Total         ""            ""        ""

:::

Visualization of Homogeneity of Slopes

Parallel vs. Non-Parallel Slopes

Partridge Example: Data Overview

ANCOVA on Longevity of Male Fruitflies

Partridge Example: Step 1 - Test Homogeneity

• Testing Homogeneity of Slopes
  • First, fit interaction model.
  • The p-value for the interaction term (THORAX:treatment) is xl.
  • Since this value is > 0.05, we conclude the slopes are homogeneous.
  • We can proceed to Step 3 and fit the simpler additive model.

# Step 1: Fit model with interaction term
homo_slopes_model <- lm(LONGEV ~ THORAX * treatment, data = partridge)
# Step 2: Test significance of interaction using Type III SS
Anova(homo_slopes_model, type = "III")
## Anova Table (Type III tests)
## 
## Response: LONGEV
##                   Sum Sq  Df F value    Pr(>F)    
## (Intercept)        755.6   1  6.6320   0.01128 *  
## THORAX            3486.3   1 30.5999 2.017e-07 ***
## treatment           36.9   4  0.0810   0.98805    
## THORAX:treatment    42.5   4  0.0933   0.98441    
## Residuals        13102.1 115                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Extract the p-value for interaction
p_interaction <- Anova(homo_slopes_model, type = "III")["THORAX:treatment", "Pr(>F)"]

Partridge Example: Step 3 ANCOVA Analysis

• Since the interaction was not significant, we re-run the model without it.
  • Interpretation:
    • THORAX: covariate highly significant
    • thorax length explains significant amount of variance in longevity
    • treatment: treatment effect, after adjusting for thorax length, is highly significant

# Step 3: Fit the ANCOVA model (additive model)
ancova_model <- lm(LONGEV ~ THORAX + treatment, data = partridge)
# View the final ANCOVA table
Anova(ancova_model, type = "III")
## Anova Table (Type III tests)
## 
## Response: LONGEV
##              Sum Sq  Df F value    Pr(>F)    
## (Intercept)  2234.9   1  20.233 1.605e-05 ***
## THORAX      13168.9   1 119.219 < 2.2e-16 ***
## treatment    9611.5   4  21.753 1.719e-13 ***
## Residuals   13144.7 119                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Partridge Example: Getting Adjusted Means

• The ANOVA table tells us the treatments are different, but not which ones.

  • Use emmeans package to get the adjusted means / perform pairwise comparisons

  • The “estimated marginal mean” (or adjusted mean)

  • predicted longevity for each group, assuming an average thorax length

  • This is the “level playing field” comparison.

# Get adjusted means using emmeans
library(emmeans)
adjusted_means <- emmeans(ancova_model, "treatment")
adjusted_means
##  treatment                       emmean   SE  df lower.CL upper.CL
##  No females                        65.4 2.11 119     61.3     69.6
##  One virgin female daily           61.5 2.11 119     57.3     65.7
##  Eight virgin females daily        64.2 2.10 119     60.0     68.3
##  One inseminated female daily      54.5 2.11 119     50.3     58.7
##  Eight inseminated females daily   41.6 2.12 119     37.4     45.8
## 
## Confidence level used: 0.95

Partridge Example: Pairwise Comparisons

• Pairwise Comparisons of Adjusted Means to see which groups differ

  • Interpretation:

    • The pairs() output gives p-values for all combinations.

    • The plot shows these comparisons visually. Red arrows connect any groups that are NOT significantly different.

    • Here, “One virgin female” and “One inseminated female” are not different from each other, but all other comparisons are significant.

# Pairwise comparisons of adjusted means
pairs(adjusted_means, adjust = "tukey")
##  contrast                                                       estimate   SE
##  No females - One virgin female daily                               3.93 3.00
##  No females - Eight virgin females daily                            1.28 2.98
##  No females - One inseminated female daily                         10.95 3.00
##  No females - Eight inseminated females daily                      23.88 2.97
##  One virgin female daily - Eight virgin females daily              -2.65 2.98
##  One virgin female daily - One inseminated female daily             7.02 2.97
##  One virgin female daily - Eight inseminated females daily         19.95 3.01
##  Eight virgin females daily - One inseminated female daily          9.67 2.98
##  Eight virgin females daily - Eight inseminated females daily      22.60 2.99
##  One inseminated female daily - Eight inseminated females daily    12.93 3.01
##   df t.ratio p.value
##  119   1.311  0.6849
##  119   0.428  0.9929
##  119   3.650  0.0035
##  119   8.031  <.0001
##  119  -0.891  0.8996
##  119   2.361  0.1336
##  119   6.636  <.0001
##  119   3.249  0.0129
##  119   7.560  <.0001
##  119   4.298  0.0003
## 
## P value adjustment: tukey method for comparing a family of 5 estimates

Partridge Example: Pairwise Comparisons

• Pairwise Comparisons of Adjusted Means to see which groups differ

  • Interpretation:

    • The pairs() output gives p-values for all combinations.

    • The plot shows these comparisons visually. Red arrows connect any groups that are NOT significantly different.

    • Here, “One virgin female” and “One inseminated female” are not different from each other, but all other comparisons are significant.

# Plot adjusted means with confidence intervals
plot(adjusted_means, comparisons = TRUE) + 
  labs(title = "Adjusted Means by Treatment",
       x = "Adjusted Longevity (days)",
       y = "Treatment") +
  theme_minimal()

Visualizing ANCOVA Results

Visualization Options for ANCOVA

What if Slopes are HETEROGENEOUS? (p < 0.05)

Interpretation of Heterogeneous Slopes

• With heterogeneous slopes, interpretation shifts from “adjusted means” to “regions of significance”:
  • Initial > Low Food when cube root body volume > 2.95
    For large urchins, initial sample has wider sutures than low food urchins
  • High Food > Initial when cube root body volume > 1.81
    For most urchins, high food treatment results in wider sutures than initial samples.
  • High Food > Low Food when cube root body volume > 2.07
    For most medium to large urchins, high food results in wider sutures than low food.
• Biological interpretation is food regime affects suture width differently depending on urchin size

Checking Assumptions in R

• For assumptions 2 and 3, we check the final model (the additive one, if appropriate).
• Interpretation:
  • Residuals vs. Fitted: The red line should be roughly flat at 0. Looks good.
  • Normal Q-Q: Points should fall along the dashed line. Looks good.
  • Scale-Location: Red line should be flat (tests homogeneity of variance). Looks good.
  • Residuals vs. Leverage: No points should be in the top/bottom right (high leverage/influence). Looks good.
• Conclusion: The assumptions are met.

Robust ANCOVA Approaches

# Example of rank-transformed ANCOVA
partridge$rank_LONGEV <- rank(partridge$LONGEV)
partridge$rank_THORAX <- rank(partridge$THORAX)

# Fit ANCOVA on ranked data
rank_ancova <- lm(rank_LONGEV ~ 
                  rank_THORAX + treatment, 
                data = partridge)

Anova(rank_ancova, type = "III")
## Anova Table (Type III tests)
## 
## Response: rank_LONGEV
##             Sum Sq  Df F value    Pr(>F)    
## (Intercept)  31251   1  63.616 1.045e-12 ***
## rank_THORAX  54314   1 110.564 < 2.2e-16 ***
## treatment    40447   4  20.584 6.536e-13 ***
## Residuals    58458 119                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Compare p-values for treatment effect
cat("P-value (parametric):", 
   Anova(ancova_model, type="III")["treatment", "Pr(>F)"],  "\n")
## P-value (parametric): 1.719359e-13
cat("P-value (rank-based):", 
    Anova(rank_ancova, type="III")["treatment", "Pr(>F)"])
## P-value (rank-based): 6.536035e-13