Analysis of Covariance (ANCOVA)

Analysis of Covariance (ANCOVA) extends the general linear model to include both categorical independent variables (factors) and continuous independent variables (covariates). ANCOVA combines the explanatory power of ANOVA for group comparisons with the precision-enhancing capability of regression through covariate adjustment. By statistically controlling for the linear influence of one or more covariates, ANCOVA reduces error variance, increases statistical power, and enables more precise comparisons of treatment means adjusted for confounding variables.

Mathematical Framework

Homogeneity of Regression Slopes

A critical assumption of ANCOVA is that the regression slopes relating the covariate to the dependent variable are equal across all groups.

Adjusted Means

ANCOVA produces adjusted (least-squares) means that account for the covariate, providing fairer comparisons when groups have different covariate distributions.

Effect Size: Partial Eta-Squared

Assumptions and Diagnostics

Python Implementation

import numpy as np
import pandas as pd
from scipy import stats
import statsmodels.api as sm
from statsmodels.formula.api import ols
import matplotlib.pyplot as plt

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

n_per_group = 25
n_groups = 3

# Group-specific true effects
group_effects = {'A': 5.0, 'B': 8.0, 'C': 3.0}
beta_true = 0.75  # True common slope

# Generate covariate (pretest) - different distributions per group
covariate_means = {'A': 75, 'B': 65, 'C': 70}
data_list = []

for group, mean_x in covariate_means.items():
    X = np.random.normal(mean_x, 10, n_per_group)
    Y = 50 + group_effects[group] + beta_true * X + np.random.normal(0, 5, n_per_group)
    
    df = pd.DataFrame({
        'Pretest': X,
        'Posttest': Y,
        'Group': group
    })
    data_list.append(df)

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

print("Data Summary:")
print(data.groupby('Group').agg(['mean', 'std', 'count']))

# Test homogeneity of slopes (interaction model)
model_interaction = ols('Posttest ~ C(Group) * Pretest', data=data).fit()
print("\nInteraction Model (Test Homogeneity of Slopes):")
print(f"Interaction F-statistic: {model_interaction.fvalue:.3f}")
print(f"Interaction p-value: {model_interaction.f_pvalue:.4f}")

# ANCOVA model (assuming parallel slopes)
model_ancova = ols('Posttest ~ C(Group) + Pretest', data=data).fit()
print("\nANCOVA Model Summary:")
print(model_ancova.summary())

# Extract ANCOVA table
from statsmodels.stats.anova import anova_lm
ancova_table = anova_lm(model_ancova, typ=2)
print("\nANCOVA Table (Type II):")
print(ancova_table)

# Compute partial eta-squared
ss_group = ancova_table.loc['C(Group)', 'sum_sq']
ss_covariate = ancova_table.loc['Pretest', 'sum_sq']
ss_error = ancova_table.loc['Residual', 'sum_sq']

eta2_partial_group = ss_group / (ss_group + ss_error)
eta2_partial_covariate = ss_covariate / (ss_covariate + ss_error)

print(f"\nPartial eta-squared (Group): {eta2_partial_group:.4f}")
print(f"Partial eta-squared (Pretest): {eta2_partial_covariate:.4f}")

# Compute adjusted means
group_means = data.groupby('Group')['Posttest'].mean()
group_cov_means = data.groupby('Group')['Pretest'].mean()
overall_cov_mean = data['Pretest'].mean()
beta_hat = model_ancova.params['Pretest']

adjusted_means = {}
for group in ['A', 'B', 'C']:
    adj_mean = group_means[group] - beta_hat * (group_cov_means[group] - overall_cov_mean)
    adjusted_means[group] = adj_mean

print("\nUnadjusted vs Adjusted Means:")
print(f"{'Group':<8} {'Unadjusted':<12} {'Adjusted':<12} {'Pretest Mean':<14}")
for group in ['A', 'B', 'C']:
    print(f"{group:<8} {group_means[group]:<12.3f} {adjusted_means[group]:<12.3f} "
          f"{group_cov_means[group]:<14.1f}")

# Standard error of adjusted means
ms_error = ancova_table.loc['Residual', 'mean_sq']
n_g = data.groupby('Group').size()

se_adjusted = {}
for group in ['A', 'B', 'C']:
    se = np.sqrt(ms_error * (1/n_g[group] + 
                            (group_cov_means[group] - overall_cov_mean)**2 / 
                            ancova_table.loc['Pretest', 'sum_sq']))
    se_adjusted[group] = se

print("\nAdjusted Means with 95% CI:")
for group in ['A', 'B', 'C']:
    ci_low = adjusted_means[group] - 1.96 * se_adjusted[group]
    ci_high = adjusted_means[group] + 1.96 * se_adjusted[group]
    print(f"  Group {group}: {adjusted_means[group]:.3f} [{ci_low:.3f}, {ci_high:.3f}]")

# Assumption diagnostics
# 1. Linearity - scatter plots
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
colors = {'A': 'blue', 'B': 'red', 'C': 'green'}

for i, group in enumerate(['A', 'B', 'C']):
    mask = data['Group'] == group
    axes[i].scatter(data[mask]['Pretest'], data[mask]['Posttest'], 
                   c=colors[group], alpha=0.6, s=50)
    
    # Add regression line
    x_range = np.linspace(data[mask]['Pretest'].min(), 
                         data[mask]['Pretest'].max(), 100)
    y_pred = model_ancova.params['Intercept'] + model_ancova.params[f'C(Group)[T.{group}]'] + \
             beta_hat * x_range
    axes[i].plot(x_range, y_pred, 'k-', linewidth=2)
    
    axes[i].set_xlabel('Pretest (Covariate)')
    axes[i].set_ylabel('Posttest')
    axes[i].set_title(f'Group {group}')
    axes[i].grid(True, alpha=0.3)

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

# 2. Homogeneity of slopes visualization
fig, ax = plt.subplots(figsize=(8, 6))

for group in ['A', 'B', 'C']:
    mask = data['Group'] == group
    ax.scatter(data[mask]['Pretest'], data[mask]['Posttest'], 
              c=colors[group], alpha=0.6, label=group, s=50)
    
    # Separate regression lines
    z = np.polyfit(data[mask]['Pretest'], data[mask]['Posttest'], 1)
    x_line = np.linspace(data[mask]['Pretest'].min(), 
                        data[mask]['Pretest'].max(), 100)
    ax.plot(x_line, np.polyval(z, x_line), color=colors[group], 
           linestyle='--', linewidth=2)

# Common regression line
x_common = np.linspace(data['Pretest'].min(), data['Pretest'].max(), 100)
y_common = model_ancova.params['Intercept'] + beta_hat * x_common
ax.plot(x_common, y_common, 'k-', linewidth=2, label='Common slope (ANCOVA)')

ax.set_xlabel('Pretest (Covariate)')
ax.set_ylabel('Posttest')
ax.set_title('Test of Homogeneity of Regression Slopes')
ax.legend()
ax.grid(True, alpha=0.3)
plt.savefig('ancova_slopes.png', dpi=150)
plt.show()

# 3. Residual diagnostics
residuals = model_ancova.resid
fitted = model_ancova.fittedvalues

fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# Residuals vs Fitted
axes[0].scatter(fitted, residuals, alpha=0.6)
axes[0].axhline(0, color='red', linestyle='--')
axes[0].set_xlabel('Fitted Values')
axes[0].set_ylabel('Residuals')
axes[0].set_title('Residuals vs Fitted')
axes[0].grid(True, alpha=0.3)

# Q-Q plot
stats.probplot(residuals, dist="norm", plot=axes[1])
axes[1].set_title('Normal Q-Q Plot')
axes[1].grid(True, alpha=0.3)

# Residuals by group
for group in ['A', 'B', 'C']:
    mask = data['Group'] == group
    group_residuals = residuals[mask]
    axes[2].boxplot(group_residuals, positions=[list(['A', 'B', 'C']).index(group)],
                   widths=0.5)

axes[2].set_xticklabels(['A', 'B', 'C'])
axes[2].set_xlabel('Group')
axes[2].set_ylabel('Residuals')
axes[2].set_title('Residuals by Group')
axes[2].axhline(0, color='red', linestyle='--')
axes[2].grid(True, alpha=0.3)

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

# Levene's test for homoscedasticity
levene_stat, levene_p = stats.levene(*[data[data['Group'] == g]['Posttest'] 
                                       for g in ['A', 'B', 'C']])
print(f"\nLevene's test for homoscedasticity: F = {levene_stat:.3f}, p = {levene_p:.4f}")

# Shapiro-Wilk test for normality of residuals
shapiro_stat, shapiro_p = stats.shapiro(residuals)
print(f"Shapiro-Wilk test for normality: W = {shapiro_stat:.3f}, p = {shapiro_p:.4f}")

# Comparison: ANOVA vs ANCOVA
from scipy.stats import f_oneway

# Separate ANOVAs (one-way)
anova_groups = [data[data['Group'] == g]['Posttest'] for g in ['A', 'B', 'C']]
f_anova, p_anova = f_oneway(*anova_groups)

print("\nComparison: One-way ANOVA vs ANCOVA")
print(f"One-way ANOVA: F({n_groups-1}, {len(data)-n_groups}) = {f_anova:.3f}, p = {p_anova:.4f}")
print(f"ANCOVA: F({n_groups-1}, {len(data)-n_groups-1}) = {ancova_table.loc['C(Group)', 'F']:.3f}, "
      f"p = {ancova_table.loc['C(Group)', 'PR(>F)']:.4f}")
print(f"Error reduction: {(1 - ss_error/((len(data)-n_groups)*np.var(data['Posttest'], ddof=1)))*100:.1f}%")