---
title: "Week 8: Contrasts and Estimable Functions"
format:
html:
toc: true
toc-depth: 3
number-sections: true
code-fold: false
code-tools: true
bibliography: ../references.bib
---
::: {.callout-note icon=false}
## Learning Objectives
By the end of this week, you will be able to:
1. **Define** **linear contrasts** and **estimable functions** in **overparameterized models**
2. **Determine** which functions are **estimable** using the **row space criterion**
3. **Construct** and **test** hypotheses involving **single** and **multiple contrasts**
4. **Create** **orthogonal contrast sets** and understand their role in **partitioning sums of squares**
5. **Compute** **standard errors** and **test statistics** for estimable functions
:::
## Why Contrasts Matter in Animal Breeding
After fitting a one-way ANOVA model (Week 7), we can test whether there are overall differences among treatment groups. But in animal breeding and genetics, we often have **specific questions** that go beyond "Are the groups different?"
**Common questions in livestock science:**
- **Breed comparisons**: Is Angus different from Hereford? Are British breeds different from Continental breeds?
- **Genetic lines**: Is Line A superior to the average of Lines B, C, and D?
- **Sire evaluation**: Which sires produce significantly different offspring performance?
- **Treatment effects**: Is a new diet formulation different from the control?
- **Selection strategies**: Do selected lines differ from unselected control lines?
Each of these questions can be formulated as a **contrast** — a specific linear combination of treatment means. Contrasts allow us to test targeted hypotheses rather than just omnibus "any difference" tests.
::: {.callout-note}
## Connection to Week 7
In Week 7, we learned about one-way ANOVA and the F-test for overall treatment effects. The F-test tells us **if** there are differences among groups, but not **which** groups differ or **how** they differ.
Contrasts provide the tools to answer these specific questions. They are the natural next step after establishing that groups differ overall.
:::
::: {.callout-tip}
## Why Geneticists Care About Specific Comparisons
In genetic improvement programs, we rarely care about "Are all breeds different?" Instead, we want to know:
- Which crosses produce the best offspring? (comparing specific crosses)
- Is a new selection line better than the base population? (selected vs. control)
- Which sire should we use for specific mating goals? (pairwise sire comparisons)
- Do different production systems affect a trait? (system A vs. B, or A vs. average of others)
Contrasts are the statistical tool that translates breeding questions into testable hypotheses.
:::
### The Estimability Problem
In Week 7, we encountered an important issue: when using the **effects model** for one-way ANOVA:
$$
y_{ij} = \mu + \alpha_i + e_{ij}
$$
the design matrix $\mathbf{X}$ is **not full rank**. This means:
- The normal equations have infinitely many solutions
- Individual parameters ($\mu$, $\alpha_1$, $\alpha_2$, etc.) are not uniquely estimable
- Different constraints (or different generalized inverses) give different estimates
However, **differences** between treatment effects (contrasts like $\alpha_1 - \alpha_2$) **are** uniquely estimable. The concept of **estimability** helps us understand which functions of parameters have unique, interpretable estimates.
This week, we'll develop the mathematical theory of estimable functions and show how to identify, estimate, and test them.
## Mathematical Theory
### Linear Combinations and Contrasts
A **linear combination** of parameters is any expression of the form:
$$
\psi = c_1 \beta_1 + c_2 \beta_2 + \cdots + c_p \beta_p = \mathbf{c}'\boldsymbol{\beta}
$$ {#eq-linear-combination}
where $\mathbf{c}$ is a $p \times 1$ vector of constants and $\boldsymbol{\beta}$ is the $p \times 1$ vector of model parameters.
**Examples** (for one-way ANOVA with cell means $\mu_1, \mu_2, \mu_3, \mu_4$):
1. **Pairwise comparison**: $\mu_1 - \mu_2$, where $\mathbf{c} = [1, -1, 0, 0]'$
2. **One vs. average of others**: $\mu_1 - \frac{1}{3}(\mu_2 + \mu_3 + \mu_4)$, where $\mathbf{c} = [1, -\frac{1}{3}, -\frac{1}{3}, -\frac{1}{3}]'$
3. **Group comparison**: $\frac{\mu_1 + \mu_2}{2} - \frac{\mu_3 + \mu_4}{2}$, where $\mathbf{c} = [\frac{1}{2}, \frac{1}{2}, -\frac{1}{2}, -\frac{1}{2}]'$
A **contrast** is a special type of linear combination where the coefficients sum to zero:
$$
\sum_{i=1}^{p} c_i = 0 \quad \text{or} \quad \mathbf{c}'\mathbf{1} = 0
$$
where $\mathbf{1}$ is a vector of ones.
::: {.callout-note}
## Why "Contrast"?
The term "contrast" comes from the idea of comparing or contrasting groups. When $\sum c_i = 0$, we're comparing weighted averages of some groups against weighted averages of others — no single group stands alone.
In the context of treatment means (cell means model), all pairwise comparisons and group comparisons are contrasts.
:::
### Estimable Functions
::: {.callout-important}
## Definition: Estimable Function
A linear function $\mathbf{c}'\boldsymbol{\beta}$ is **estimable** if and only if there exists a vector $\mathbf{a}$ such that:
$$
\mathbf{c}' = \mathbf{a}'\mathbf{X}
$$
Equivalently, $\mathbf{c}'$ is in the **row space** of $\mathbf{X}$.
:::
**Intuition**: A function is estimable if it can be expressed as a linear combination of the rows of the design matrix. This ensures that the function can be uniquely estimated from the data, regardless of which generalized inverse we use.
#### Key Theorems
**Theorem 1**: If $\mathbf{X}$ is **full rank** (rank $= p$), then **all** linear functions $\mathbf{c}'\boldsymbol{\beta}$ are estimable.
**Theorem 2**: If $\mathbf{X}$ is **not full rank** (rank $< p$), then only certain linear functions are estimable.
**Theorem 3**: For the **cell means model** in one-way ANOVA, all contrasts $\sum c_i \mu_i$ (where $\sum c_i = 0$) are estimable.
**Theorem 4**: For the **effects model** with sum-to-zero constraint:
- Individual parameters $\mu$ and $\alpha_i$ are **NOT estimable**
- But differences $\alpha_i - \alpha_j$ **ARE estimable**
### Estimation and Variance
::: {.callout-important}
## Key Property: Uniqueness
If $\mathbf{c}'\boldsymbol{\beta}$ is estimable, then its estimate $\mathbf{c}'\mathbf{b}$ is **unique** — it does not depend on which generalized inverse $(\mathbf{X}'\mathbf{X})^{-}$ is used to solve the normal equations.
:::
For an estimable function $\mathbf{c}'\boldsymbol{\beta}$, the variance of its estimate is:
$$
\text{Var}(\mathbf{c}'\mathbf{b}) = \mathbf{c}'(\mathbf{X}'\mathbf{X})^{-}\mathbf{c} \, \sigma^2
$$ {#eq-var-contrast}
We estimate $\sigma^2$ from the residuals:
$$
\hat{\sigma}^2 = \text{MSE} = \frac{\text{SSE}}{n - r(\mathbf{X})}
$$
### Hypothesis Testing
#### Single Contrast
To test $H_0: \mathbf{c}'\boldsymbol{\beta} = 0$ vs. $H_a: \mathbf{c}'\boldsymbol{\beta} \neq 0$:
**Test statistic**:
$$
t = \frac{\mathbf{c}'\mathbf{b}}{\text{se}(\mathbf{c}'\mathbf{b})} = \frac{\mathbf{c}'\mathbf{b}}{\sqrt{\mathbf{c}'(\mathbf{X}'\mathbf{X})^{-}\mathbf{c} \, \hat{\sigma}^2}}
$$ {#eq-t-test-contrast}
Under $H_0$, $t \sim t(n - r(\mathbf{X}))$.
#### Multiple Contrasts
To test $H_0: \mathbf{C}\boldsymbol{\beta} = \mathbf{0}$ where $\mathbf{C}$ is a $q \times p$ matrix:
**F-test statistic**:
$$
F = \frac{\mathbf{b}'\mathbf{C}'[\mathbf{C}(\mathbf{X}'\mathbf{X})^{-}\mathbf{C}']^{-1}\mathbf{C}\mathbf{b} / q}{\text{MSE}}
$$ {#eq-f-test-contrasts}
Under $H_0$, $F \sim F(q, n - r(\mathbf{X}))$.
### Orthogonal Contrasts
Two estimable functions $\mathbf{c}_1'\boldsymbol{\beta}$ and $\mathbf{c}_2'\boldsymbol{\beta}$ are **orthogonal** if:
$$
\mathbf{c}_1'(\mathbf{X}'\mathbf{X})^{-}\mathbf{c}_2 = 0
$$
**For balanced designs**, this simplifies to $\mathbf{c}_1'\mathbf{c}_2 = 0$.
::: {.callout-important}
## Property: Orthogonal Contrasts Partition Variation
For a set of $g-1$ **orthogonal contrasts** in a one-way ANOVA with $g$ groups:
$$
\text{SS(Treatments)} = \text{SS}(\mathbf{c}_1) + \text{SS}(\mathbf{c}_2) + \cdots + \text{SS}(\mathbf{c}_{g-1})
$$
where $\text{SS}(\mathbf{c}_j) = \frac{(\mathbf{c}_j'\mathbf{b})^2}{\mathbf{c}_j'(\mathbf{X}'\mathbf{X})^{-}\mathbf{c}_j}$
:::
## Small Numerical Example: Beef Cattle Carcass Weight
We have data on carcass weights (kg) for 12 beef steers from 4 breeds (3 steers per breed).
```{r}
#| label: load-beef-data
#| message: false
# Load data
beef <- read.csv("data/beef_breeds_small.csv")
print(beef)
# Summary by breed
library(dplyr)
beef_summary <- beef %>%
group_by(breed) %>%
summarise(n = n(), mean = mean(carcass_weight_kg), sd = sd(carcass_weight_kg))
print(beef_summary)
```
### Fit One-Way ANOVA
```{r}
#| label: beef-anova
#| message: false
# Design matrix (cell means model)
y <- beef$carcass_weight_kg
n <- length(y)
X <- model.matrix(~ breed - 1, data = beef)
colnames(X) <- levels(beef$breed)
# Solve normal equations
XtX <- t(X) %*% X
Xty <- t(X) %*% y
b <- solve(XtX) %*% Xty
cat("Breed means:\n")
print(b)
# ANOVA components
y_hat <- X %*% b
residuals <- y - y_hat
SST <- sum((y - mean(y))^2)
SSE <- sum(residuals^2)
SSM <- SST - SSE
MSE <- SSE / (n - ncol(X))
cat("\nMSE =", MSE, "\n")
```
### Define and Test Contrasts
```{r}
#| label: beef-contrasts
#| message: false
# Three orthogonal contrasts
c1 <- c(1, 0, -1, 0) # Angus vs Hereford
c2 <- c(1, -1, 1, -1) # British vs Continental
c3 <- c(0, 1, 0, -1) # Charolais vs Simmental
# Check orthogonality
cat("Orthogonality checks:\n")
cat("c1' c2 =", sum(c1 * c2), "\n")
cat("c1' c3 =", sum(c1 * c3), "\n")
cat("c2' c3 =", sum(c2 * c3), "\n")
# Estimate contrasts
est1 <- sum(c1 * b)
est2 <- sum(c2 * b)
est3 <- sum(c3 * b)
# Standard errors (balanced design)
n_per_group <- 3
se1 <- sqrt(sum(c1^2) / n_per_group * MSE)
se2 <- sqrt(sum(c2^2) / n_per_group * MSE)
se3 <- sqrt(sum(c3^2) / n_per_group * MSE)
# t-statistics
t1 <- est1 / se1
t2 <- est2 / se2
t3 <- est3 / se3
# Results
results <- data.frame(
Contrast = c("Angus - Hereford", "British - Continental", "Charolais - Simmental"),
Estimate = round(c(est1, est2, est3), 2),
SE = round(c(se1, se2, se3), 3),
t = round(c(t1, t2, t3), 3),
p = round(2 * (1 - pt(abs(c(t1, t2, t3)), 8)), 4)
)
print(results)
```
::: {.callout-note}
## Biological Interpretation
- **Angus vs. Hereford**: No significant difference (both British breeds)
- **British vs. Continental**: Highly significant! Continental breeds ~39 kg heavier
- **Charolais vs. Simmental**: No significant difference (both Continental)
This pattern matches breeding history: Continental breeds selected for growth and carcass weight.
:::
## Realistic Application: Layer Strain Egg Production
```{r}
#| label: layers-data
#| message: false
# Load layer strain data
layers <- read.csv("data/layer_strains_egg.csv")
# Summary
layers_summary <- layers %>%
group_by(strain) %>%
summarise(n = n(), mean = mean(eggs_per_month), sd = sd(eggs_per_month))
print(layers_summary)
```
### Fit Model and Define Contrasts
```{r}
#| label: layers-contrasts
#| message: false
# Fit ANOVA
y_l <- layers$eggs_per_month
X_l <- model.matrix(~ strain - 1, data = layers)
b_l <- solve(t(X_l) %*% X_l) %*% t(X_l) %*% y_l
resid_l <- y_l - X_l %*% b_l
MSE_l <- sum(resid_l^2) / (length(y_l) - ncol(X_l))
# Contrasts (alphabetical order: HyLineW36, ISABrown, LohmannBrown, RhodeIslandRed, WhiteLeghorn)
c1_l <- c(-0.25, -0.25, -0.25, 1, -0.25) # Heritage vs Commercial
c2_l <- c(0.5, -0.5, -0.5, 0, 0.5) # White vs Brown
c3_l <- c(-1, 0, 0, 0, 1) # WL vs HyLine
c4_l <- c(0, 1, -1, 0, 0) # ISA vs Lohmann
# Test contrasts
n_per <- 10
test_contrast <- function(c_vec) {
est <- sum(c_vec * b_l)
se <- sqrt(sum(c_vec^2) / n_per * MSE_l)
t_stat <- est / se
p_val <- 2 * (1 - pt(abs(t_stat), 45))
return(c(est = est, se = se, t = t_stat, p = p_val))
}
r1 <- test_contrast(c1_l)
r2 <- test_contrast(c2_l)
r3 <- test_contrast(c3_l)
r4 <- test_contrast(c4_l)
layer_results <- data.frame(
Contrast = c("Heritage vs Commercial", "White vs Brown", "WL vs HyLine", "ISA vs Lohmann"),
rbind(r1, r2, r3, r4)
)
colnames(layer_results) <- c("Contrast", "Estimate", "SE", "t", "p")
layer_results[,2:5] <- round(layer_results[,2:5], 4)
print(layer_results)
```
::: {.callout-important}
## Commercial Implications
- **Heritage vs Commercial**: Significant! Modern strains produce 5+ more eggs/month
- **White vs Brown**: Significant difference favoring white-egg strains
- **Within types**: HyLine outperforms Leghorn; ISA and Lohmann similar
:::
## Building Contrast Testing Tools
```{r}
#| label: custom-functions
#| message: false
# Function 1: Check estimability
is_estimable <- function(c, X, tol = 1e-10) {
library(MASS)
XtX_ginv <- ginv(t(X) %*% X)
projection <- t(X) %*% X %*% XtX_ginv %*% c
return(all(abs(projection - c) < tol))
}
# Function 2: Test single contrast
test_single_contrast <- function(c, b, X, sigma2, method = "balanced") {
n <- nrow(X)
p <- ncol(X)
estimate <- sum(c * b)
if (method == "balanced") {
variance <- sum(c^2) / (n / p) * sigma2
} else {
library(MASS)
variance <- t(c) %*% ginv(t(X) %*% X) %*% c * sigma2
}
se <- sqrt(variance)
df <- n - qr(X)$rank
t_stat <- estimate / se
p_value <- 2 * (1 - pt(abs(t_stat), df))
return(list(estimate = estimate, se = se, t = t_stat, df = df, p_value = p_value))
}
# Function 3: Check orthogonality
check_orthogonality <- function(C, X = NULL, method = "simple") {
if (method == "simple") {
return(C %*% t(C))
} else {
library(MASS)
XtX_ginv <- ginv(t(X) %*% X)
return(C %*% XtX_ginv %*% t(C))
}
}
# Test functions
cat("Testing is_estimable():\n")
cat("c1 in beef data:", is_estimable(c1, X), "\n\n")
cat("Testing test_single_contrast():\n")
result <- test_single_contrast(c1, b, X, MSE, "balanced")
cat("Estimate:", round(result$estimate, 2), "\n")
cat("SE:", round(result$se, 3), "\n")
cat("t:", round(result$t, 3), "\n")
cat("p:", round(result$p_value, 4), "\n\n")
cat("Testing check_orthogonality():\n")
C_beef <- rbind(c1, c2, c3)
ortho_mat <- check_orthogonality(C_beef)
cat("Off-diagonal elements:", ortho_mat[row(ortho_mat) != col(ortho_mat)], "\n")
```
## Exercises
### Conceptual
1. Prove that in a one-way ANOVA effects model with sum-to-zero constraint, the contrast $\alpha_i - \alpha_j$ is estimable even though individual $\alpha_i$ are not.
2. Explain why orthogonal contrasts partition the treatment sum of squares in balanced designs.
3. In genetic evaluation, explain the difference between testing "Sire A is different from Sire B" versus stating "Sire A has breeding value +10 kg".
### Computational
1. **Hand Calculations**: Yorkshire (11, 12, 11), Landrace (10, 11, 10), Duroc (9, 10, 9)
- Test Yorkshire = Landrace
- Test Yorkshire = (Landrace + Duroc)/2
- Construct 95% CIs
2. **Non-orthogonal Contrasts**: Create contrasts that are NOT orthogonal and show their SSs don't sum to SS(Treatments).
3. **Verify Estimability**: Use different generalized inverses and show estimable functions remain constant.
### Applied
1. **Broiler Strains**: Design 4 meaningful contrasts comparing fast-growing vs. heritage breeds.
2. **Sheep Breeds**: Test if Dorset and Suffolk (meat breeds) are similar and together different from wool breeds.
## Summary
✓ **Defined** linear contrasts as $\mathbf{c}'\boldsymbol{\beta}$ where $\sum c_i = 0$
✓ **Determined** estimability using the row space criterion
✓ **Computed** estimates and standard errors for contrasts
✓ **Tested** contrasts using t-tests and F-tests
✓ **Constructed** orthogonal contrast sets
✓ **Interpreted** results in livestock breeding contexts
✓ **Built** custom R functions for contrast testing
**Key concepts**:
- Estimability depends on design matrix structure
- Contrasts are estimable even when individual parameters are not
- Orthogonal contrasts partition variation in balanced designs
- Generalized inverses give different solutions, but estimable functions are unique
- Contrasts answer specific research questions
## Looking Ahead
- **Week 9**: Contrasts in factorial models with interactions
- **Week 10**: Adjusted means as estimable functions in ANCOVA
- **Week 12**: Estimability in unbalanced designs
---
**Previous**: [Week 7: One-Way ANOVA](../Week07_ANOVA_OneWay/Week07_ANOVA_OneWay.qmd)
**Next**: [Week 9: Two-Way ANOVA and Factorial Models](../Week09_TwoWayANOVA/Week09_TwoWayANOVA.qmd)