Practical ANCOVA Example with R

1. Raw Data

Here we have a sample dataset of 20 students who were randomly assigned to either a traditional or interactive teaching method. The data includes students' prior GPAs and final exam scores.

2. Initial Analysis in R

In this step, we load the data into R, summarize the key statistics, and create a scatterplot to visualize the relationship between prior GPA and final exam scores by teaching method:

# Load required packages
library(tidyverse)
library(car) # for assumption testing
library(effectsize) # for effect size
library(emmeans) # for adjusted means

# Read and examine the data
data <- read_csv("teaching_methods.csv")

# Get summary statistics by group
data |>
  group_by(teaching_method) |>
  summarise(
    n = n(),
    mean_gpa = mean(prior_gpa),
    sd_gpa = sd(prior_gpa),
    mean_exam = mean(final_exam),
    sd_exam = sd(final_exam)
  )

# A tibble: 2 × 6
  teaching_method     n mean_gpa sd_gpa mean_exam sd_exam
  <chr>           <int>    <dbl>  <dbl>     <dbl>   <dbl>
1 Interactive        10     3.35  0.303      84.6    4.72
2 Traditional        10     3.25  0.303      78.6    4.84

Visualizing the data can help identify patterns and relationships between variables. Let's create a scatterplot with a regression line:

ggplot(data, aes(x = prior_gpa, y = final_exam, color = teaching_method)) +
  geom_point(size = 3, alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE) +
  theme_minimal() +
  labs(
    title = "Final Exam Scores by Prior GPA and Teaching Method",
    x = "Prior GPA",
    y = "Final Exam Score",
    color = "Teaching Method"
  )

Scatterplot with regression line for final exam scores by prior GPA and teaching method.

The scatterplot shows a positive linear relationship between prior GPA and final exam scores, with a clear distinction between teaching methods. Interactive method students tend to have higher exam scores. The regression lines are roughly parallel, indicating homogeneity of regression slopes. Next, we'll proceed with the ANCOVA analysis and check assumptions.

3. Assumption Testing

# Test linearity and homogeneity of regression slopes
model_interaction <- aov(final_exam ~ teaching_method * prior_gpa, data = data)
summary(model_interaction)

                          Df Sum Sq Mean Sq F value   Pr(>F)    
teaching_method            1  180.0   180.0 395.341 1.05e-12 ***
prior_gpa                  1  403.4   403.4 886.043 1.94e-15 ***
teaching_method:prior_gpa  1    0.1     0.1   0.213    0.651    
Residuals                 16    7.3     0.5                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Test normality of residuals
model_initial <- aov(final_exam ~ teaching_method + prior_gpa, data = data)
shapiro.test(residuals(model_initial))

	Shapiro-Wilk normality test

data:  residuals(model_initial)
W = 0.98373, p-value = 0.9728

# Check homogeneity of variances
leveneTest(final_exam ~ teaching_method, data = data)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  1  0.0338 0.8562
      18    

4. Fitting ANCOVA Model

With our assumptions met, we can proceed with the ANCOVA analysis to compare teaching methods while controlling for prior GPA:

# Fit ANCOVA model
ancova_model <- aov(final_exam ~ teaching_method + prior_gpa, data = data)
summary(ancova_model)

                Df Sum Sq Mean Sq F value   Pr(>F)    
teaching_method  1  180.0   180.0   414.5 2.24e-13 ***
prior_gpa        1  403.4   403.4   929.1 2.80e-16 ***
Residuals       17    7.4     0.4                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The ANCOVA model shows a significant effect of teaching method (p < .001) and prior GPA (p < .001) on final exam scores. The adjusted means will provide a clearer picture of the group differences while controlling for GPA.

5. Effect Sizes

Let's calculate effect sizes to quantify the practical significance of our results:

# Get effect sizes
eta_squared(ancova_model)

# Effect Size for ANOVA (Type I)

Parameter       | Eta2 (partial) |       95% CI
-----------------------------------------------
teaching_method |           0.96 | [0.92, 1.00]
prior_gpa       |           0.98 | [0.96, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

The effect sizes show a large impact of teaching method (η² = 0.96) and prior GPA (η² = 0.98) on final exam scores. This indicates a substantial practical significance of both variables.

6. Adjusted Means

Finally, let's calculate the adjusted means for each teaching method based on prior GPA:

# Calculate adjusted means
adjusted_means <- emmeans(ancova_model, "teaching_method")
pairs(adjusted_means)
confint(adjusted_means)

  contrast                  estimate    SE df t.ratio p.value
 Interactive - Traditional     4.44 0.299 17  14.831  <.0001

# Confidence Intervals
 teaching_method emmean   SE df lower.CL upper.CL
 Interactive      83.82 0.21 17    83.38    84.26
 Traditional      79.38 0.21 17    78.94    79.82

Confidence level used: 0.95

The adjusted means show a significant difference between teaching methods, with the interactive method scoring 4.44 points higher than the traditional method. The confidence intervals confirm the reliability of this difference.

7. Interpretation and Reporting

8. Practical Implications