---
title: "Week 13: Special Topics I"
subtitle: "Unbalanced Data, Generalized Inverses, and Constraint Systems"
author: "Linear Models for Animal Breeding and Genetics"
date: today
format:
html:
toc: true
toc-depth: 4
toc-location: left
number-sections: true
number-depth: 3
code-fold: false
code-tools: true
code-copy: true
highlight-style: github
theme: cosmo
css: ../styles.css
embed-resources: true
bibliography: ../references.bib
---
::: {.callout-note icon=false}
## Learning Objectives
After completing this week, you will be able to:
1. **Diagnose** when **unbalanced data** causes problems and select appropriate **analytical strategies**
2. **Distinguish** between different types of **generalized inverses** and understand when each is appropriate
3. **Apply** different **constraint systems** (set-to-zero, sum-to-zero, custom) and interpret their effects on **parameter estimates**
4. **Identify** **estimable** vs. **non-estimable functions** and correctly interpret results from **rank-deficient models**
5. **Analyze** real **genetic evaluation data** with unequal information and understand the limitations of **fixed-effects least squares**
:::
---
# Introduction: Making Sense of Messy Data
In Weeks 1-12, we've built a solid foundation in linear models: from simple regression through ANOVA, diagnostics, and rank-deficient models. Throughout, we've often used balanced designs for clarity. But real livestock data is rarely balanced.
Consider these common scenarios in animal breeding and genetics:
- **Dairy sire evaluation**: Popular sires have 50+ daughters; young sires have only 3-5
- **Multi-location trials**: Not all breeds are available at all farms due to management constraints
- **Feedlot studies**: Pens vary in size; animals are lost due to health issues or data quality problems
- **Genetic evaluation**: Family sizes vary dramatically; some animals have no progeny records
This week synthesizes concepts from multiple previous weeks to provide practical strategies for handling these real-world challenges:
::: {.callout-note}
## Connections to Previous Weeks
- **Week 2**: Generalized inverses (mathematical foundation)
- **Week 7**: One-way ANOVA with balanced data (ideal case)
- **Week 8**: Estimable functions and contrasts (theoretical framework)
- **Week 12**: Non-full rank models and unequal subclass numbers (immediate prerequisite)
If you need to review rank deficiency, generalized inverses, or estimability, revisit those weeks before proceeding.
:::
## What We'll Cover This Week
This week is organized around four main topics:
**Topic A: Strategic Approaches to Unbalanced Data**
When does unbalance cause problems? What solutions work best in different situations?
**Topic B: Types of Generalized Inverses**
Reflexive g-inverses, Moore-Penrose inverses, and conditional inverses—what's the difference and when does it matter?
**Topic C: Constraint Systems**
How do set-to-zero and sum-to-zero constraints affect parameter estimates? How should we interpret results?
**Topic D: Interpreting Non-Estimable Parameters**
What can we legitimately report from rank-deficient models? What must we avoid?
**Capstone Example: Dairy Sire Evaluation**
A comprehensive analysis of 10 sires with 3-25 daughters each, illustrating all concepts and previewing the need for mixed models (Week 14).
Let's begin.
---
# Topic A: Strategic Approaches to Unbalanced Data
## When is Unbalanced Data a Problem?
Not all unbalanced data creates problems. To understand when issues arise, let's first review what makes balanced designs special.
### Properties of Balanced Designs
In a **balanced design**, all cells have equal sample sizes: $n_{ij} = n$ for all combinations of factors.
For example, in a 2×2 factorial (Factor A: 2 levels, Factor B: 2 levels):
```
Factor B
B₁ B₂
Factor A
A₁ n n
A₂ n n
```
This balance creates several nice properties:
1. **Orthogonality**: After centering, $\mathbf{X}'\mathbf{X}$ is (nearly) diagonal
2. **SS equivalence**: Type I = Type II = Type III sums of squares
3. **Independence**: Main effects are independent of interactions
4. **Simple interpretation**: Effects have clear, unconfounded meanings
### Consequences of Unbalanced Designs
When sample sizes differ across cells, we lose these properties:
1. **Loss of orthogonality**: $\mathbf{X}'\mathbf{X}$ has substantial off-diagonal elements
2. **SS divergence**: Type I, II, and III SS can differ dramatically
3. **Confounding**: Effects become correlated; interpretation becomes ambiguous
4. **Estimability issues**: With missing cells, some parameters may not be estimable
Let's see this with a concrete example.
### Numerical Demonstration: Balanced vs. Unbalanced
Consider broiler body weight (kg) by strain (A, B) and sex (Male, Female).
#### Balanced Case
```{r}
# Balanced data: n = 2 per cell
strain <- rep(c("A", "A", "B", "B"), 2)
sex <- rep(c("Male", "Female"), each = 4)
weight <- c(2.1, 2.3, # A-Male
1.8, 1.9, # A-Female
2.4, 2.5, # B-Male
1.9, 2.0) # B-Female
# Show cell counts
table(strain, sex)
# Fit model
fit_bal <- lm(weight ~ strain * sex)
anova(fit_bal) # Type I SS
```
#### Unbalanced Case
Now suppose we have unequal cell sizes:
```{r}
# Unbalanced data
strain_unb <- c(rep("A", 3), rep("A", 2), rep("B", 2), rep("B", 3))
sex_unb <- c(rep("Male", 3), rep("Female", 2), rep("Male", 2), rep("Female", 3))
weight_unb <- c(2.1, 2.3, 2.2, # A-Male (n=3)
1.8, 1.9, # A-Female (n=2)
2.4, 2.5, # B-Male (n=2)
1.9, 2.0, 2.1) # B-Female (n=3)
# Cell counts
table(strain_unb, sex_unb)
# Fit model
fit_unb <- lm(weight_unb ~ strain_unb * sex_unb)
# Type I SS (strain first)
cat("\n=== Type I SS: Strain first ===\n")
anova(fit_unb)
# Type I SS (sex first)
fit_unb2 <- lm(weight_unb ~ sex_unb * strain_unb)
cat("\n=== Type I SS: Sex first ===\n")
anova(fit_unb2)
# Type III SS
library(car)
cat("\n=== Type III SS ===\n")
Anova(fit_unb, type = 3)
```
**Key Observations**:
1. In the balanced case, Type I SS (strain first) = Type I SS (sex first)
2. In the unbalanced case, **order matters**! SS values differ depending on which factor enters first
3. Type III SS adjust each effect for all others, regardless of order
### Examining the Design Matrix Structure
Let's look at why this happens by examining $\mathbf{X}'\mathbf{X}$:
```{r}
# Design matrix for unbalanced case
X_unb <- model.matrix(~ strain_unb * sex_unb)
XtX_unb <- t(X_unb) %*% X_unb
cat("X'X for unbalanced design:\n")
print(XtX_unb)
# Note: off-diagonal elements are non-zero
# This indicates non-orthogonality (correlation between predictors)
# For comparison, balanced case:
X_bal <- model.matrix(~ strain * sex)
XtX_bal <- t(X_bal) %*% X_bal
cat("\nX'X for balanced design:\n")
print(XtX_bal)
# Balanced design has more structure (some zeros in off-diagonals after centering)
```
::: {.callout-important}
## When Unbalance is Problematic
Unbalanced data itself is not always problematic. The issue arises when unbalance creates:
1. **Confounding**: Effects cannot be separated clearly
2. **Rank deficiency**: Missing cells make $\mathbf{X}'\mathbf{X}$ singular
3. **Interpretation ambiguity**: Type I/II/III SS tell different stories
If you have unequal $n_{ij}$ but all cells are represented and you use appropriate methods, unbalance is manageable.
:::
## Strategic Solutions
When faced with unbalanced data, you have several strategies. Choose based on the nature and severity of the imbalance.
### Strategy 1: Use Estimable Functions (Contrasts)
**Principle**: Focus on what can be uniquely estimated, regardless of parameterization.
Even when individual parameters (like $\alpha_i$ in $\mu + \alpha_i$) are not uniquely estimable, **contrasts** (differences between parameters) often are.
**Example**: In sire evaluation with unequal progeny:
- Individual sire effects $s_i$ are not uniquely estimable (confounded with $\mu$)
- But sire **differences** $s_i - s_j$ are estimable and biologically meaningful
**When to use**:
- Rank-deficient models (overparameterized)
- Primary interest is in comparisons, not absolute values
**R Implementation**:
```{r}
# Example: 3 breeds, unequal sample sizes
breed <- c(rep("A", 3), rep("B", 2), rep("C", 4))
yield <- c(10, 12, 11, # Breed A
14, 15, # Breed B
16, 17, 18, 19) # Breed C
fit <- lm(yield ~ breed)
# Use emmeans for estimable contrasts
library(emmeans)
emm <- emmeans(fit, "breed")
# All pairwise contrasts (all estimable!)
pairs(emm)
```
### Strategy 2: Cell Means Model
**Principle**: Use a model that is **always full rank**, regardless of balance.
The cell means model:
$$y_{ij} = \mu_i + e_{ij}$$
estimates each group mean directly. If there are $g$ groups with non-zero counts, then:
- $\mathbf{X}$ is $n \times g$ with rank $r(\mathbf{X}) = g$ (full rank!)
- All $\mu_i$ are uniquely estimable (no constraints needed)
- Comparisons are done post-hoc via contrasts
**Trade-off**: More parameters, but all estimable and interpretation is straightforward.
**When to use**:
- Severe unbalance with missing cells
- Want to avoid rank deficiency issues
- Primary interest is in group means, not "effects"
**R Implementation**:
```{r}
# Cell means model: remove intercept with "-1"
fit_cm <- lm(yield ~ breed - 1)
summary(fit_cm)
# Each coefficient is a breed mean (all estimable!)
coef(fit_cm)
# Same contrasts as before
emm_cm <- emmeans(fit_cm, "breed")
pairs(emm_cm) # Identical to effects model contrasts
```
### Strategy 3: Choose Appropriate Sums of Squares Type
**Principle**: Different SS types test different hypotheses. Choose based on your research question.
#### Type I (Sequential) Sums of Squares
- Effects added **sequentially** in the order specified
- SS(A) = variation explained by A alone
- SS(B|A) = additional variation explained by B after accounting for A
- SS(A×B|A,B) = additional variation explained by interaction after A and B
**Order matters!** SS(A) ≠ SS(A|B)
**When to use**:
- Logical ordering of predictors exists (e.g., covariates before treatments)
- Hierarchical model building
**R Default**: `anova(fit)` gives Type I SS
#### Type II Sums of Squares
- Each main effect adjusted for other main effects, **but not interactions**
- SS(A|B) but not SS(A|B, A×B)
- More powerful when interactions are absent or negligible
**When to use**:
- No significant interactions
- Main effects are primary interest
**R Implementation**: `car::Anova(fit, type=2)`
#### Type III (Marginal) Sums of Squares
- Each effect adjusted for **all other effects**, including interactions
- SS(A|B, A×B): effect of A given everything else
- Most conservative; commonly used for unbalanced data
**When to use**:
- Unbalanced designs (most common choice)
- Interactions may be present
- Want to test each effect "in the presence of" all others
**R Implementation**: `car::Anova(fit, type=3)`
::: {.callout-warning}
## Choosing SS Type
With severe unbalance and missing cells, Type III SS can be difficult to interpret. The "adjusted" means may not correspond to observable data patterns.
**Best practice**:
1. Start with Type III for unbalanced data
2. Examine cell means and sample sizes
3. Use estimable contrasts for inference
4. Report which SS type you used
:::
#### Comparison Example
```{r}
# Unbalanced 2-way ANOVA
strain_2 <- c(rep("A", 5), rep("A", 2), rep("B", 3), rep("B", 6))
sex_2 <- c(rep("M", 5), rep("F", 2), rep("M", 3), rep("F", 6))
bw <- c(2.1, 2.2, 2.3, 2.2, 2.4, # A-M
1.9, 1.8, # A-F
2.5, 2.6, 2.4, # B-M
2.0, 2.1, 1.9, 2.0, 2.2, 2.1) # B-F
fit_comp <- lm(bw ~ strain_2 * sex_2)
cat("=== Type I SS (Strain first) ===\n")
print(anova(fit_comp))
cat("\n=== Type I SS (Sex first) ===\n")
fit_comp2 <- lm(bw ~ sex_2 * strain_2)
print(anova(fit_comp2))
cat("\n=== Type III SS ===\n")
print(Anova(fit_comp, type=3))
```
Notice: Type I SS differ depending on order; Type III SS are invariant to order.
### Strategy 4: Weighted Least Squares
**Principle**: Account for heterogeneous variances across observations.
Sometimes unbalance reflects different **reliability** or **precision** of measurements:
- Pen means with varying pen sizes: larger pens are more reliable
- Group means with unequal within-group variation
- Measurements with known precision (e.g., genomic predictions)
**Weighted Least Squares (WLS)** minimizes:
$$\sum_{i=1}^n w_i (y_i - \hat{y}_i)^2$$
where $w_i = 1/\sigma_i^2$ (inversely proportional to variance).
**Normal equations**:
$$\mathbf{X}'\mathbf{W}\mathbf{X} \mathbf{b} = \mathbf{X}'\mathbf{W}\mathbf{y}$$
where $\mathbf{W} = \text{diag}(w_1, \ldots, w_n)$
**When to use**:
- Heteroscedastic errors (unequal variances)
- Grouped data: weight by group size
- Known reliability: weight by precision
**Example**: Pen-average ADG with different pen sizes
```{r}
# Pen-level data (means)
pen_id <- 1:5
pen_adg <- c(0.85, 0.90, 0.88, 0.87, 0.92) # kg/day
pen_size <- c(8, 12, 6, 10, 9) # number of pigs per pen
diet <- c("A", "B", "A", "B", "A")
# Unweighted regression (WRONG - ignores unequal precision)
fit_unwt <- lm(pen_adg ~ diet)
summary(fit_unwt)
# Weighted regression (CORRECT - weight by pen size)
# Variance of pen mean ~ 1/n, so weight by n
fit_wt <- lm(pen_adg ~ diet, weights = pen_size)
summary(fit_wt)
# Compare standard errors
cat("\nUnweighted SE(diet effect):", summary(fit_unwt)$coef[2,2], "\n")
cat("Weighted SE(diet effect):", summary(fit_wt)$coef[2,2], "\n")
# Weighted analysis gives more appropriate SE
```
**Key point**: When observations represent group means with different group sizes, **always weight by $n_i$**.
## Summary: Choosing a Strategy
```{r}
#| echo: false
library(knitr)
strategy_table <- data.frame(
Situation = c(
"Unequal n, all cells present, want ANOVA",
"Missing cells, rank deficiency",
"No clear reference group",
"Want to avoid estimability issues entirely",
"Grouped data (e.g., pen means)",
"Known heterogeneous precision"
),
Strategy = c(
"Type III SS, interpret estimable contrasts",
"Cell means model + post-hoc contrasts",
"Sum-to-zero constraints + contrasts",
"Cell means model",
"Weighted least squares (weight by n)",
"Weighted least squares (weight by 1/σ²)"
)
)
kable(strategy_table, caption = "Choosing an Analytical Strategy for Unbalanced Data")
```
::: {.callout-tip}
## General Recommendations
1. **Always check cell counts first**: Use `table()` to see the data structure
2. **Default to Type III SS** for unbalanced ANOVA unless you have specific reasons otherwise
3. **Focus on estimable contrasts** for inference, not individual parameter estimates
4. **Weight appropriately** when precision varies across observations
5. **Report your choices** clearly (which SS type, which constraints)
:::
---
# Topic B: Types of Generalized Inverses
In Week 2, we introduced the generalized inverse as a solution to rank-deficient systems. In Week 12, we used it to solve non-full rank normal equations. Now let's dive deeper into **different types** of generalized inverses and when each is appropriate.
## Review: Why Generalized Inverses?
For a square matrix $\mathbf{A}$, the regular inverse $\mathbf{A}^{-1}$ exists if and only if $\mathbf{A}$ is full rank.
For **rank-deficient** matrices (common with categorical predictors in overparameterized models), $\mathbf{A}^{-1}$ does not exist. But we can find a **generalized inverse** $\mathbf{A}^-$ that satisfies certain properties.
Consider the normal equations:
$$\mathbf{X}'\mathbf{X} \mathbf{b} = \mathbf{X}'\mathbf{y}$$
If $r(\mathbf{X}'\mathbf{X}) < p$ (rank deficient), this system has **infinitely many solutions**. Any solution of the form:
$$\mathbf{b} = (\mathbf{X}'\mathbf{X})^- \mathbf{X}'\mathbf{y}$$
will satisfy the normal equations, but the specific $\mathbf{b}$ values depend on which generalized inverse we use.
**Key insight**:
- $\mathbf{X}\mathbf{b}$ (fitted values) is **unique** regardless of which $(\mathbf{X}'\mathbf{X})^-$ we use
- Individual $b_j$ values are **not unique**
- **Estimable functions** $\mathbf{c}'\mathbf{b}$ are **unique**
Let's explore different types of generalized inverses.
## Reflexive Generalized Inverse (g-inverse)
### Definition
$\mathbf{A}^-$ is a **reflexive generalized inverse** (or simply **g-inverse**) of $\mathbf{A}$ if:
$$\mathbf{A} \mathbf{A}^- \mathbf{A} = \mathbf{A}$$
This is the **minimum requirement** for a generalized inverse.
### Properties
- **Not unique**: For a rank-deficient matrix, infinitely many g-inverses exist
- **Sufficient for solving normal equations**: Any g-inverse allows us to find a solution to $\mathbf{A}\mathbf{x} = \mathbf{b}$
- **Different g-inverses give different solutions**, but all satisfy $\mathbf{A}\mathbf{x} = \mathbf{b}$
### Why This Property Works
Suppose $\mathbf{A} \mathbf{A}^- \mathbf{A} = \mathbf{A}$. We want to solve $\mathbf{X}'\mathbf{X} \mathbf{b} = \mathbf{X}'\mathbf{y}$.
Let $\mathbf{b} = (\mathbf{X}'\mathbf{X})^- \mathbf{X}'\mathbf{y}$. Then:
\begin{align}
\mathbf{X}'\mathbf{X} \mathbf{b} &= \mathbf{X}'\mathbf{X} [(\mathbf{X}'\mathbf{X})^- \mathbf{X}'\mathbf{y}] \\
&= [\mathbf{X}'\mathbf{X} (\mathbf{X}'\mathbf{X})^- \mathbf{X}'\mathbf{X}] \mathbf{X}'\mathbf{y} \quad \text{(rearranging)}\\
&= \mathbf{X}'\mathbf{X} (\mathbf{X}'\mathbf{y}) \quad \text{(using } \mathbf{A}\mathbf{A}^-\mathbf{A} = \mathbf{A})\\
&= \mathbf{X}'\mathbf{y} \quad \checkmark
\end{align}
So $\mathbf{b}$ is indeed a solution!
### Computing a g-inverse by Row Reduction
We can find a g-inverse using Gaussian elimination, making arbitrary choices when we encounter linear dependence.
**Example**:
$$\mathbf{A} = \begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix}$$
Clearly $r(\mathbf{A}) = 1$ (rows are identical).
To find a g-inverse, augment with identity and row reduce:
$$\left[\begin{array}{cc|cc} 1 & 1 & 1 & 0 \\ 1 & 1 & 0 & 1 \end{array}\right]$$
Row reduce:
$$\left[\begin{array}{cc|cc} 1 & 1 & 1 & 0 \\ 0 & 0 & -1 & 1 \end{array}\right]$$
Since second column is all zeros (linear dependence), we have a **free variable**. Set $x_2 = 0$ arbitrarily:
$$\mathbf{A}^-_1 = \begin{bmatrix} 1 & 0 \\ 0 & 0 \end{bmatrix}$$
But we could also set $x_1 = 0$:
$$\mathbf{A}^-_2 = \begin{bmatrix} 0 & 0 \\ 0 & 1 \end{bmatrix}$$
Or make other choices:
$$\mathbf{A}^-_3 = \begin{bmatrix} 0.5 & 0.5 \\ 0.5 & 0.5 \end{bmatrix}$$
**Verification**: Let's check all satisfy $\mathbf{A}\mathbf{A}^-\mathbf{A} = \mathbf{A}$:
```{r}
# Define matrix A
A <- matrix(c(1, 1, 1, 1), nrow=2, byrow=TRUE)
print(A)
# Three different g-inverses
A_inv1 <- matrix(c(1, 0, 0, 0), nrow=2, byrow=TRUE)
A_inv2 <- matrix(c(0, 0, 0, 1), nrow=2, byrow=TRUE)
A_inv3 <- matrix(c(0.5, 0.5, 0.5, 0.5), nrow=2, byrow=TRUE)
# Verify property: A A- A = A
cat("A A_inv1 A:\n")
print(A %*% A_inv1 %*% A)
cat("\nA A_inv2 A:\n")
print(A %*% A_inv2 %*% A)
cat("\nA A_inv3 A:\n")
print(A %*% A_inv3 %*% A)
# All equal to A!
```
### Different Solutions to Normal Equations
Different g-inverses give different solutions, but **fitted values are the same**:
```{r}
# Suppose we want to solve Ax = y
y <- c(2, 2) # Note: must be in column space of A (here, y1 = y2)
# Three solutions
b1 <- A_inv1 %*% y
b2 <- A_inv2 %*% y
b3 <- A_inv3 %*% y
cat("Solution 1:", b1, "\n")
cat("Solution 2:", b2, "\n")
cat("Solution 3:", b3, "\n")
# All different! But check fitted values:
cat("\nFitted values (Ab):\n")
cat("A b1:", A %*% b1, "\n")
cat("A b2:", A %*% b2, "\n")
cat("A b3:", A %*% b3, "\n")
# All the same!
```
::: {.callout-important}
## Key Insight: What's Unique and What's Not
For rank-deficient normal equations $\mathbf{X}'\mathbf{X}\mathbf{b} = \mathbf{X}'\mathbf{y}$:
- **Fitted values** $\mathbf{X}\mathbf{b}$ are **unique** ✓
- **Individual parameter estimates** $b_j$ are **not unique** ✗
- **Estimable functions** $\mathbf{c}'\mathbf{b}$ are **unique** ✓
The choice of g-inverse matters only for interpretation of individual parameters, not for estimable functions or predictions.
:::
## Moore-Penrose Pseudoinverse
While infinitely many reflexive g-inverses exist, there is **exactly one** generalized inverse with additional nice properties: the **Moore-Penrose inverse**.
### Definition
$\mathbf{A}^+$ is the **Moore-Penrose inverse** (or **pseudoinverse**) of $\mathbf{A}$ if it satisfies **all four** of these properties:
1. $\mathbf{A} \mathbf{A}^+ \mathbf{A} = \mathbf{A}$ (reflexive g-inverse property)
2. $\mathbf{A}^+ \mathbf{A} \mathbf{A}^+ = \mathbf{A}^+$ (reflexive for $\mathbf{A}^+$ itself)
3. $(\mathbf{A} \mathbf{A}^+)' = \mathbf{A} \mathbf{A}^+$ (symmetric)
4. $(\mathbf{A}^+ \mathbf{A})' = \mathbf{A}^+ \mathbf{A}$ (symmetric)
### Key Property
The Moore-Penrose inverse is **unique**: there is only one matrix satisfying all four conditions.
### Computing via Singular Value Decomposition (SVD)
The most common way to compute $\mathbf{A}^+$ is via SVD.
**SVD Decomposition**:
$$\mathbf{A} = \mathbf{U} \mathbf{D} \mathbf{V}'$$
where:
- $\mathbf{U}$ ($n \times n$) and $\mathbf{V}$ ($p \times p$) are orthogonal matrices
- $\mathbf{D}$ ($n \times p$) is diagonal with singular values $\sigma_1 \geq \sigma_2 \geq \ldots \geq \sigma_r > 0$, then zeros
**Moore-Penrose Inverse**:
$$\mathbf{A}^+ = \mathbf{V} \mathbf{D}^+ \mathbf{U}'$$
where $\mathbf{D}^+$ has $1/\sigma_i$ for non-zero $\sigma_i$, and zeros elsewhere.
**Example**: Our previous rank-1 matrix:
```{r}
# A = [1 1; 1 1]
A <- matrix(c(1, 1, 1, 1), nrow=2, byrow=TRUE)
# SVD
svd_A <- svd(A)
cat("Singular values:\n")
print(svd_A$d)
# Note: σ₁ = 2, σ₂ = 0 (rank 1)
# Moore-Penrose inverse using R
library(MASS)
A_plus <- ginv(A)
cat("\nMoore-Penrose inverse:\n")
print(A_plus)
# Verify all 4 properties
cat("\n1. A A+ A = A?\n")
print(all.equal(A %*% A_plus %*% A, A))
cat("\n2. A+ A A+ = A+?\n")
print(all.equal(A_plus %*% A %*% A_plus, A_plus))
cat("\n3. (A A+)' = A A+?\n")
print(all.equal(t(A %*% A_plus), A %*% A_plus))
cat("\n4. (A+ A)' = A+ A?\n")
print(all.equal(t(A_plus %*% A), A_plus %*% A))
```
### Special Property: Minimum Norm Solution
Among all solutions to $\mathbf{A}\mathbf{x} = \mathbf{b}$, the Moore-Penrose solution:
$$\mathbf{x} = \mathbf{A}^+ \mathbf{b}$$
has the **smallest norm** (smallest sum of squares):
$$||\mathbf{x}||^2 = \sum_j x_j^2 \text{ is minimized}$$
This is often desirable in statistical applications.
```{r}
# Compare norms of different solutions
y <- c(2, 2)
b1 <- A_inv1 %*% y
b2 <- A_inv2 %*% y
b3 <- A_inv3 %*% y
b_plus <- A_plus %*% y
cat("Norm of solution 1:", sqrt(sum(b1^2)), "\n")
cat("Norm of solution 2:", sqrt(sum(b2^2)), "\n")
cat("Norm of solution 3:", sqrt(sum(b3^2)), "\n")
cat("Norm of M-P solution:", sqrt(sum(b_plus^2)), "\n")
# Moore-Penrose gives minimum norm
```
::: {.callout-note}
## R's Default
In R, `MASS::ginv()` computes the **Moore-Penrose inverse** using SVD. This is the most commonly used generalized inverse in practice because:
1. It's unique (reproducible)
2. It's numerically stable (SVD is robust)
3. It gives minimum norm solutions
4. It handles any matrix (not just square)
Unless you have a specific reason to use a different g-inverse, use `ginv()`.
:::
## Conditional Inverse (Constraint-Based)
A third approach is to **impose constraints** to make the system full rank, then use a regular inverse.
### Concept
Instead of solving the rank-deficient system:
$$\mathbf{X}'\mathbf{X} \mathbf{b} = \mathbf{X}'\mathbf{y}$$
We augment with constraints $\mathbf{C}'\mathbf{b} = \mathbf{c}$ to create a full-rank system:
$$\begin{bmatrix} \mathbf{X}'\mathbf{X} & \mathbf{C}' \\ \mathbf{C} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbf{b} \\ \boldsymbol{\lambda} \end{bmatrix} = \begin{bmatrix} \mathbf{X}'\mathbf{y} \\ \mathbf{c} \end{bmatrix}$$
where $\boldsymbol{\lambda}$ are Lagrange multipliers.
With appropriate constraints (e.g., $\sum \alpha_i = 0$ for one-way ANOVA), this system becomes full rank and has a **unique solution**.
### Example: Sum-to-Zero Constraint
For one-way ANOVA effects model $y_{ij} = \mu + \alpha_i + e_{ij}$:
**Constraint**: $\sum_{i=1}^g \alpha_i = 0$
This can be written as:
$$\mathbf{C}' = [0, 1, 1, \ldots, 1]$$
so $\mathbf{C}'\mathbf{b} = \alpha_1 + \alpha_2 + \ldots + \alpha_g = 0$
```{r}
# Example: 3 groups
group <- c(rep(1, 3), rep(2, 2), rep(3, 4))
y_ex <- c(10, 12, 11, # Group 1
14, 15, # Group 2
16, 17, 18, 19) # Group 3
# For demonstration, use cell means model first to get all parameters
X_cell <- model.matrix(~ 0 + factor(group)) # Cell means (3 params, full rank)
XtX_cell <- t(X_cell) %*% X_cell
Xty_cell <- t(X_cell) %*% y_ex
cat("Cell means model: Rank =", qr(XtX_cell)$rank, "out of", ncol(XtX_cell), "\n")
cat("This is full rank, so all 3 group means are estimable.\n\n")
# Now show effects model with constraint
# Build full effects model manually: μ, α₁, α₂, α₃ (4 params, rank 3)
# X has: column of 1s (μ), then indicators for each group
X_eff <- cbind(1, model.matrix(~ 0 + factor(group)))
XtX_eff <- t(X_eff) %*% X_eff
Xty_eff <- t(X_eff) %*% y_ex
cat("Effects model: Rank =", qr(XtX_eff)$rank, "out of", ncol(XtX_eff), "\n")
cat("Rank deficient!\n\n")
# Apply sum-to-zero constraint: α₁ + α₂ + α₃ = 0
# Constraint: [0, 1, 1, 1] · [μ, α₁, α₂, α₃]' = 0
C <- matrix(c(0, 1, 1, 1), nrow=1) # 1x4 constraint matrix
# Augmented system
aug_mat <- rbind(
cbind(XtX_eff, t(C)), # XtX is 4x4, t(C) is 4x1
cbind(C, 0) # C is 1x4, 0 is scalar
)
aug_rhs <- c(Xty_eff, 0)
# Solve (now full rank!)
solution <- solve(aug_mat, aug_rhs)
b_constrained <- solution[1:4]
lambda <- solution[5]
cat("Constrained solution:\n")
cat("μ =", b_constrained[1], "\n")
cat("α₁ =", b_constrained[2], "\n")
cat("α₂ =", b_constrained[3], "\n")
cat("α₃ =", b_constrained[4], "\n")
cat("Sum of α's =", sum(b_constrained[2:4]), "(should be ~0)\n")
```
### Properties
- **Unique solution** for given constraint
- Different constraints give different $\mathbf{b}$ values
- **Estimable functions unchanged** by constraint choice
- Allows interpretation of individual parameters under chosen constraint
## Comparison of Generalized Inverse Types
```{r}
#| echo: false
ginv_table <- data.frame(
Type = c("Reflexive g-inverse", "Moore-Penrose", "Conditional (constrained)"),
Uniqueness = c("Many exist", "Unique", "Unique (for given constraints)"),
Computation = c("Row reduction (arbitrary choices)", "SVD: A = UDV' → A⁺ = VD⁺U'", "Augmented system with constraints"),
`R Function` = c("Custom (many options)", "MASS::ginv()", "solve() on augmented system"),
`Use Case` = c("Quick solution when you don't care about individual parameters",
"Prefer unique, stable solution; R default",
"Want specific interpretation (e.g., sum-to-zero effects)"),
check.names = FALSE
)
kable(ginv_table, caption = "Comparison of Generalized Inverse Types")
```
## Practical Comparison with Real Data
Let's compare all three approaches on a realistic example:
**Data**: Milk yield by breed (unequal sample sizes)
```{r}
# One-way ANOVA: 3 breeds, unequal n
breed_ex <- c(rep("Holstein", 3), rep("Jersey", 2), rep("Brown Swiss", 4))
milk_yield <- c(
30, 32, 31, # Holstein (mean = 31)
24, 25, # Jersey (mean = 24.5)
28, 29, 27, 28 # Brown Swiss (mean = 28)
)
# Effects model design matrix
X_milk <- model.matrix(~ breed_ex)
XtX_milk <- t(X_milk) %*% X_milk
Xty_milk <- t(X_milk) %*% milk_yield
cat("=== Design Matrix Rank ===\n")
cat("r(X'X) =", qr(XtX_milk)$rank, "out of", ncol(XtX_milk), "\n")
cat("Rank deficient by", ncol(XtX_milk) - qr(XtX_milk)$rank, "\n\n")
# Approach 1: Use a reflexive g-inverse (many possible)
# We'll use R's ginv which gives Moore-Penrose
b_ginv <- ginv(XtX_milk) %*% Xty_milk
cat("=== Solution 1: Moore-Penrose (ginv) ===\n")
cat("μ =", b_ginv[1], "\n")
cat("α_Jersey =", b_ginv[2], "\n")
cat("α_BrownSwiss =", b_ginv[3], "\n\n")
# Approach 2: Sum-to-zero constraint
options(contrasts = c("contr.sum", "contr.poly"))
fit_sum <- lm(milk_yield ~ breed_ex)
b_sum <- coef(fit_sum)
cat("=== Solution 2: Sum-to-Zero Constraint ===\n")
cat("μ =", b_sum[1], "\n")
cat("α_Holstein =", b_sum[2], "\n")
cat("α_Jersey =", b_sum[3], "\n")
cat("α_BrownSwiss = ", -(b_sum[2] + b_sum[3]), "(computed)\n\n")
# Approach 3: Set-to-zero constraint (Holstein = reference)
options(contrasts = c("contr.treatment", "contr.poly"))
fit_treat <- lm(milk_yield ~ breed_ex)
b_treat <- coef(fit_treat)
cat("=== Solution 3: Set-to-Zero (Holstein reference) ===\n")
cat("μ (Holstein mean) =", b_treat[1], "\n")
cat("α_Holstein = 0 (reference)\n")
cat("α_Jersey =", b_treat[2], "\n")
cat("α_BrownSwiss =", b_treat[3], "\n\n")
# Compare: fitted values (should be identical!)
fitted_ginv <- X_milk %*% b_ginv
fitted_sum <- fitted(fit_sum)
fitted_treat <- fitted(fit_treat)
cat("=== Fitted Values Comparison ===\n")
cat("Max difference (ginv vs sum):", max(abs(fitted_ginv - fitted_sum)), "\n")
cat("Max difference (ginv vs treat):", max(abs(fitted_ginv - fitted_treat)), "\n")
cat("Conclusion: Fitted values are identical! ✓\n\n")
# Compare: estimable contrast (Jersey vs Holstein)
contrast_ginv <- (b_ginv[1] + b_ginv[2]) - b_ginv[1] # (μ + α_J) - (μ)
contrast_sum <- b_sum[2] + b_sum[1] # Wait, need to compute properly for sum-to-zero
# Let's use emmeans for clean contrast comparison
library(emmeans)
emm_sum <- emmeans(fit_sum, "breed_ex")
emm_treat <- emmeans(fit_treat, "breed_ex")
cat("=== Estimated Breed Means (μ + α_i) ===\n")
cat("Sum-to-zero:\n")
print(summary(emm_sum))
cat("\nSet-to-zero:\n")
print(summary(emm_treat))
# Contrasts
cat("\n=== Estimable Contrasts ===\n")
contrasts_sum <- pairs(emm_sum)
contrasts_treat <- pairs(emm_treat)
cat("Sum-to-zero:\n")
print(contrasts_sum)
cat("\nSet-to-zero:\n")
print(contrasts_treat)
cat("\nConclusion: Contrasts are identical! ✓\n")
```
::: {.callout-important}
## The Big Picture
Three very different parameter estimates ($\mu$, $\alpha_i$), but:
1. **Fitted values**: Identical across all three methods ✓
2. **Estimable functions** (breed means, contrasts): Identical ✓
3. **Biological conclusions**: Identical ✓
The choice of generalized inverse (or constraint system) affects **how we write down the answer**, not **what the answer means biologically**.
For practical data analysis:
- Use `ginv()` (Moore-Penrose) for simplicity and stability
- OR use constraint system (sum-to-zero, set-to-zero) for interpretation
- **Always report estimable functions**, not individual non-estimable parameters
:::
---
*Continued in next message due to length...*
# Topic C: Constraint Systems
In Topic B, we saw that different generalized inverses (or constraint systems) give different parameter estimates but identical estimable functions. Now let's explore **how to choose and apply constraints** for interpretability.
## Why Use Constraints?
When $\mathbf{X}'\mathbf{X}$ is rank deficient, the normal equations have infinitely many solutions. Constraints serve two purposes:
1. **Computational**: Make the system full rank so we can use regular (non-generalized) matrix inverse
2. **Interpretational**: Give individual parameters specific, interpretable meanings
**Key principle**: The constraint choice is **arbitrary** from a statistical standpoint—it doesn't change estimable functions or biological conclusions. But it **does** affect how we write down and interpret individual parameters.
## Set-to-Zero Constraints (Reference Cell Coding)
### Formulation
**Principle**: Set one level as the "reference" or "baseline"; all other effects are deviations from this reference.
For one-way ANOVA with effects model $y_{ij} = \mu + \alpha_i + e_{ij}$:
**Constraint**: $\alpha_1 = 0$ (first level is reference)
### Interpretation
With $\alpha_1 = 0$:
- Group 1: $E(y_{1j}) = \mu + \alpha_1 = \mu + 0 = \mu$
- Group 2: $E(y_{2j}) = \mu + \alpha_2$
- Group 3: $E(y_{3j}) = \mu + \alpha_3$
So:
- $\mu$ = mean of reference group (group 1)
- $\alpha_i$ = difference between group $i$ and reference group
- All $\alpha_i$ (except $\alpha_1$) are **estimable** (they're contrasts!)
### Design Matrix
```{r}
# Example: 3 groups, set-to-zero constraint
group_stz <- factor(rep(1:3, each=2))
y_stz <- c(10, 12, # Group 1
14, 15, # Group 2
16, 17) # Group 3
# Set contrasts to treatment (set-to-zero)
options(contrasts = c("contr.treatment", "contr.poly"))
# Design matrix
X_stz <- model.matrix(~ group_stz)
cat("Design matrix (set-to-zero):\n")
print(X_stz)
cat("\nColumn 1: Intercept (μ)\n")
cat("Column 2: Indicator for group 2 (α₂)\n")
cat("Column 3: Indicator for group 3 (α₃)\n")
cat("Note: No column for group 1 (it's the reference)\n")
```
Notice:
- Row 1-2: Group 1 observations have [1, 0, 0] → $\mu$ only
- Row 3-4: Group 2 observations have [1, 1, 0] → $\mu + \alpha_2$
- Row 5-6: Group 3 observations have [1, 0, 1] → $\mu + \alpha_3$
### When to Use
Set-to-zero is natural when:
- Clear **control or reference group** exists (e.g., placebo, wild-type, standard breed)
- You want to express all effects **relative to control**
- Common in experimental designs with a baseline treatment
### R Implementation
```{r}
# Fit with set-to-zero (R's default for lm)
fit_stz <- lm(y_stz ~ group_stz)
summary(fit_stz)
# Interpretation
coefs_stz <- coef(fit_stz)
cat("\nInterpretation:\n")
cat("μ (Group 1 mean):", coefs_stz[1], "\n")
cat("α₂ (Group 2 - Group 1):", coefs_stz[2], "\n")
cat("α₃ (Group 3 - Group 1):", coefs_stz[3], "\n\n")
cat("Predicted means:\n")
cat("Group 1:", coefs_stz[1], "\n")
cat("Group 2:", coefs_stz[1] + coefs_stz[2], "\n")
cat("Group 3:", coefs_stz[1] + coefs_stz[3], "\n")
```
## Sum-to-Zero Constraints
### Formulation
**Principle**: Effects sum to zero; $\mu$ represents the **overall mean**, and $\alpha_i$ are deviations from this mean.
**Constraint**: $\sum_{i=1}^g \alpha_i = 0$
Or for unbalanced designs, **weighted sum-to-zero**:
$$\sum_{i=1}^g n_i \alpha_i = 0$$
### Interpretation
With $\sum \alpha_i = 0$:
- $\mu$ = overall mean (or grand mean)
- $\alpha_i$ = deviation of group $i$ from overall mean
- More **symmetric**: no group is "special"
- Traditional ANOVA approach
### Design Matrix
```{r}
# Sum-to-zero constraint
options(contrasts = c("contr.sum", "contr.poly"))
X_sum <- model.matrix(~ group_stz)
cat("Design matrix (sum-to-zero):\n")
print(X_sum)
cat("\nColumn 1: Intercept (μ)\n")
cat("Column 2: Effect coding for group 1\n")
cat("Column 3: Effect coding for group 2\n")
cat("Note: Group 3 coded as [-1, -1] so effects sum to zero\n")
```
Notice:
- Row 1-2: Group 1 has [1, 1, 0] → $\mu + \alpha_1$
- Row 3-4: Group 2 has [1, 0, 1] → $\mu + \alpha_2$
- Row 5-6: Group 3 has [1, -1, -1] → $\mu - \alpha_1 - \alpha_2 = \mu + \alpha_3$
The last group is coded so that $\alpha_1 + \alpha_2 + \alpha_3 = 0$.
### When to Use
Sum-to-zero is natural when:
- **No natural reference group** (all groups on equal footing)
- You want to compare each group to the **overall average**
- Balanced designs (classical ANOVA)
- Interpreting "main effects" in factorial designs
### R Implementation
```{r}
# Fit with sum-to-zero
fit_sum_ex <- lm(y_stz ~ group_stz)
summary(fit_sum_ex)
coefs_sum <- coef(fit_sum_ex)
cat("\nInterpretation:\n")
cat("μ (overall mean):", coefs_sum[1], "\n")
cat("α₁ (Group 1 deviation):", coefs_sum[2], "\n")
cat("α₂ (Group 2 deviation):", coefs_sum[3], "\n")
cat("α₃ (computed):", -(coefs_sum[2] + coefs_sum[3]), "\n\n")
cat("Predicted means:\n")
cat("Group 1:", coefs_sum[1] + coefs_sum[2], "\n")
cat("Group 2:", coefs_sum[1] + coefs_sum[3], "\n")
cat("Group 3:", coefs_sum[1] - coefs_sum[2] - coefs_sum[3], "\n")
```
::: {.callout-note}
## Weighted Sum-to-Zero for Unbalanced Data
For unbalanced designs, you might prefer **weighted sum-to-zero**: $\sum_i n_i \alpha_i = 0$
This makes $\mu$ the **weighted** grand mean, and $\alpha_i$ are deviations from this weighted mean.
R's `contr.sum` uses simple (unweighted) sum-to-zero. For weighted constraints, you'd need to manually construct the constraint matrix.
:::
## Comparison: Set-to-Zero vs. Sum-to-Zero
Let's see both on the same dataset:
### Example: Sow Litter Size by Parity
```{r}
# Data: Litter size by parity
parity <- factor(c(rep(1, 4), rep(2, 3), rep(3, 2)))
litter_size <- c(10, 11, 10, 12, # Parity 1
11, 12, 13, # Parity 2
12, 13) # Parity 3
# Observed means
tapply(litter_size, parity, mean)
# Overall mean
mean(litter_size)
```
**Set-to-Zero**: Parity 1 is reference
```{r}
options(contrasts = c("contr.treatment", "contr.poly"))
fit_parity_stz <- lm(litter_size ~ parity)
cat("=== SET-TO-ZERO (Parity 1 reference) ===\n")
print(summary(fit_parity_stz)$coef)
cat("\nInterpretation:\n")
cat("Parity 1 mean =", coef(fit_parity_stz)[1], "\n")
cat("Parity 2 - Parity 1 =", coef(fit_parity_stz)[2], "\n")
cat("Parity 3 - Parity 1 =", coef(fit_parity_stz)[3], "\n")
```
**Sum-to-Zero**: Overall mean as baseline
```{r}
options(contrasts = c("contr.sum", "contr.poly"))
fit_parity_sum <- lm(litter_size ~ parity)
cat("\n=== SUM-TO-ZERO ===\n")
print(summary(fit_parity_sum)$coef)
cat("\nInterpretation:\n")
cat("Overall mean =", coef(fit_parity_sum)[1], "\n")
cat("Parity 1 deviation =", coef(fit_parity_sum)[2], "\n")
cat("Parity 2 deviation =", coef(fit_parity_sum)[3], "\n")
cat("Parity 3 deviation =", -(coef(fit_parity_sum)[2] + coef(fit_parity_sum)[3]), "(computed)\n")
```
**Compare Estimable Functions**:
```{r}
# Use emmeans to extract breed means and contrasts
library(emmeans)
emm_stz <- emmeans(fit_parity_stz, "parity")
emm_sum <- emmeans(fit_parity_sum, "parity")
cat("\n=== ESTIMATED MEANS (μ + α_i) ===\n")
cat("Set-to-zero:\n")
print(summary(emm_stz))
cat("\nSum-to-zero:\n")
print(summary(emm_sum))
cat("\n=== PAIRWISE CONTRASTS ===\n")
cat("Set-to-zero:\n")
print(pairs(emm_stz))
cat("\nSum-to-zero:\n")
print(pairs(emm_sum))
cat("\nConclusion: Means and contrasts IDENTICAL!\n")
```
### Side-by-Side Summary
```{r}
#| echo: false
library(knitr)
# Create comparison table
comp_df <- data.frame(
Parameterization = c("Set-to-Zero", "Sum-to-Zero"),
mu = c(
sprintf("%.2f (Parity 1 mean)", coef(fit_parity_stz)[1]),
sprintf("%.2f (overall mean)", coef(fit_parity_sum)[1])
),
alpha1 = c(
"0 (reference)",
sprintf("%.2f", coef(fit_parity_sum)[2])
),
alpha2 = c(
sprintf("%.2f (vs Parity 1)", coef(fit_parity_stz)[2]),
sprintf("%.2f (vs overall)", coef(fit_parity_sum)[3])
),
alpha3 = c(
sprintf("%.2f (vs Parity 1)", coef(fit_parity_stz)[3]),
sprintf("%.2f (computed)", -(coef(fit_parity_sum)[2] + coef(fit_parity_sum)[3]))
),
Parity1_mean = c(
sprintf("%.2f", mean(litter_size[parity == 1])),
sprintf("%.2f", mean(litter_size[parity == 1]))
),
Parity2_vs_1 = c(
sprintf("%.2f", mean(litter_size[parity == 2]) - mean(litter_size[parity == 1])),
sprintf("%.2f", mean(litter_size[parity == 2]) - mean(litter_size[parity == 1]))
),
check.names = FALSE
)
kable(comp_df,
caption = "Comparison of Set-to-Zero vs. Sum-to-Zero Parameterizations",
col.names = c("Constraint", "μ", "α₁", "α₂", "α₃", "Parity 1 Mean", "Parity 2 vs 1"))
```
::: {.callout-important}
## Key Insight
The two parameterizations give **completely different individual parameter values**:
- $\mu$ differs by ~0.47
- $\alpha_1$ is 0 vs -0.73
- $\alpha_2$ is 1.20 vs 0.47
But the **estimable functions** (means, contrasts) are **identical**:
- Parity 1 mean: 10.75 (both methods)
- Parity 2 vs 1: 1.20 (both methods)
- Standard errors: identical
**Conclusion**: Constraint choice affects notation, not biology!
:::
## Custom Constraints
Sometimes you want constraints tailored to your specific problem.
### Example 1: Anchoring to an External Standard
Suppose you know from breed registry standards that Holstein cows should average 9000 kg milk:
**Constraint**: $\mu + \alpha_{\text{Holstein}} = 9000$
This "anchors" your estimates to the known standard.
### Example 2: Structural Biological Constraint
In crossbreeding studies with heterosis:
**Model**: $y_{ijk} = \mu + \text{breed}_i + \text{sex}_j + (\text{breed} \times \text{sex})_{ij} + e_{ijk}$
You might constrain the crossbred mean to be midparent + heterosis:
$$\alpha_{\text{cross}} = \frac{\alpha_{\text{breed1}} + \alpha_{\text{breed2}}}{2} + h$$
where $h$ is the heterosis effect (estimable as a contrast).
### Implementation
Custom constraints are applied by augmenting the normal equations:
$$\begin{bmatrix} \mathbf{X}'\mathbf{X} & \mathbf{C}' \\ \mathbf{C} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbf{b} \\ \boldsymbol{\lambda} \end{bmatrix} = \begin{bmatrix} \mathbf{X}'\mathbf{y} \\ \mathbf{c} \end{bmatrix}$$
where:
- $\mathbf{C}$ is the constraint matrix
- $\mathbf{c}$ is the constraint values
- $\boldsymbol{\lambda}$ are Lagrange multipliers
```{r}
# Example: Custom constraint for 3-group model
# Suppose we want α₁ + α₃ = 0 (groups 1 and 3 average to zero)
# Using our earlier 3-group data
X_custom <- model.matrix(~ group_stz)
XtX_custom <- t(X_custom) %*% X_custom
Xty_custom <- t(X_custom) %*% y_stz
# Constraint: α₁ + α₃ = 0
# This is row 2 + row 4 of parameter vector = 0
# In set-to-zero coding: group1 is reference (no parameter), so actually we want
# Let's just use sum-to-zero as an example
# Actually, let's demonstrate with effects model (all α's present)
# Model: y = μ + α₁ + α₂ + α₃
# Constraint: α₁ + α₃ = 0
# Design matrix for effects model (all groups)
X_effects <- model.matrix(~ 0 + group_stz) # Cell means first
X_effects <- cbind(1, X_effects) # Add intercept
# Wait, this is getting complicated. Let me show the principle
cat("Custom constraints are typically implemented by:\n")
cat("1. Setting up the augmented system with your constraint matrix C\n")
cat("2. Solving the augmented system\n")
cat("3. Verifying the constraint is satisfied\n\n")
cat("For most applications, sum-to-zero or set-to-zero is sufficient.\n")
cat("Custom constraints are used in specialized situations like:\n")
cat(" - Anchoring to external standards\n")
cat(" - Enforcing biological relationships\n")
cat(" - Hierarchical model structures\n")
```
::: {.callout-tip}
## Practical Advice on Constraints
**For most animal breeding applications**:
1. **Use set-to-zero** when you have a clear control/reference group
- Example: "Does this new breed differ from our standard breed?"
2. **Use sum-to-zero** when all groups are on equal footing
- Example: "How do these 5 sire lines compare?"
3. **Use cell means model** to avoid constraint issues entirely
- Always full rank, all means estimable
- Do contrasts post-hoc
4. **Custom constraints**: Rare; only for special situations
5. **Always report which constraint you used** (for reproducibility)
6. **Focus on estimable functions** in your conclusions (which don't depend on constraints)
:::
---
# Topic D: Interpreting Non-Estimable Parameters
We've seen that different constraints give different parameter estimates. But which parameters **should** we interpret? This section clarifies what's legitimate to report and what's not.
## What Makes a Function Estimable?
### Mathematical Definition
A linear combination $\mathbf{c}'\boldsymbol{\beta}$ is **estimable** if and only if:
$$\mathbf{c}' = \mathbf{a}'\mathbf{X} \text{ for some vector } \mathbf{a}$$
**Equivalently**:
- $\mathbf{c}$ is in the **row space** of $\mathbf{X}$
- $\mathbf{c}'\mathbf{b}$ is the **same for all choices of generalized inverse**
- $\mathbf{c} = (\mathbf{X}'\mathbf{X})^g \mathbf{X}'\mathbf{X}$ for ANY g-inverse $(\mathbf{X}'\mathbf{X})^g$
### Intuitive Meaning
**Estimable** = "observable from the data"
If $\mathbf{c}'\boldsymbol{\beta}$ is estimable, then there exists a linear combination of the observations $\mathbf{a}'\mathbf{y}$ whose expected value is exactly $\mathbf{c}'\boldsymbol{\beta}$.
In other words: **estimable functions have direct meaning in terms of observable data**.
Non-estimable functions are artifacts of how we write down the model (parameterization choice), not real biological quantities.
## Common Estimable and Non-Estimable Functions
### One-Way ANOVA: $y_{ij} = \mu + \alpha_i + e_{ij}$
```{r}
#| echo: false
estimable_table <- data.frame(
Function = c("μ", "α_i", "μ + α_i", "α_i - α_j", "Σw_i α_i (where Σw_i=0)"),
Estimable = c("❌ No", "❌ No", "✅ Yes", "✅ Yes", "✅ Yes"),
Interpretation = c(
"Confounded with α_i",
"Confounded with μ",
"Mean of group i (observable!)",
"Difference between groups (observable!)",
"Contrast (observable!)"
)
)
knitr::kable(estimable_table, caption = "Estimable Functions in One-Way ANOVA")
```
**Why is $\mu$ not estimable?**
In effects model $y_{ij} = \mu + \alpha_i + e_{ij}$, we can't separate $\mu$ from $\alpha_i$. If we add 10 to $\mu$ and subtract 10 from all $\alpha_i$, the model is exactly the same! So $\mu$ has no unique value.
**Why is $\mu + \alpha_i$ estimable?**
$\mu + \alpha_i = E(y_{ij})$ = mean of group $i$, which we can directly estimate from data as $\bar{y}_{i \cdot}$.
### Two-Way ANOVA: $y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + e_{ijk}$
```{r}
#| echo: false
estimable_2way <- data.frame(
Function = c("μ", "α_i", "β_j", "μ + α_i + β_j + (αβ)_ij",
"α_i - α_i'", "(αβ)_ij - (αβ)_i'j - (αβ)_ij' + (αβ)_i'j'"),
Estimable = c("❌ No", "❌ (Sometimes)", "❌ (Sometimes)", "✅ Yes (if cell exists)",
"✅ Yes (usually)", "✅ Yes"),
Meaning = c("Confounded", "Main effect (confounded if unbalanced)",
"Main effect (confounded if unbalanced)", "Cell mean",
"Main effect contrast", "Interaction contrast")
)
knitr::kable(estimable_2way, caption = "Estimable Functions in Two-Way ANOVA")
```
::: {.callout-note}
## Balanced vs. Unbalanced Designs
In **balanced** designs (equal $n_{ij}$):
- Main effects $\alpha_i - \alpha_{i'}$ are estimable
- Main effects are orthogonal to interactions
In **unbalanced** designs with missing cells:
- Some main effects may not be estimable
- Confounding between main effects and interactions
**Safe approach**: Always use contrasts and check estimability!
:::
## Testing Estimability in R
### Method 1: Compare Across Different g-Inverses
If $\mathbf{c}'\mathbf{b}$ is the same for two different g-inverses, it's estimable.
```{r}
# Example: One-way ANOVA, 3 groups
library(MASS)
# Data
group_test <- factor(c(rep(1, 3), rep(2, 2), rep(3, 4)))
y_test <- c(10, 12, 11, 14, 15, 16, 17, 18, 19)
# Design matrix (effects model)
X_test <- model.matrix(~ group_test)
XtX_test <- t(X_test) %*% X_test
Xty_test <- t(X_test) %*% y_test
# Two different g-inverses
# 1. Moore-Penrose
ginv1 <- ginv(XtX_test)
b1 <- ginv1 %*% Xty_test
# 2. Use different constraint (sum-to-zero)
options(contrasts = c("contr.sum", "contr.poly"))
fit_test <- lm(y_test ~ group_test)
# Get unconstrained solution by using ginv on this
X_test2 <- model.matrix(~ group_test)
XtX_test2 <- t(X_test2) %*% X_test2
ginv2 <- ginv(XtX_test2)
b2 <- ginv2 %*% t(X_test2) %*% y_test
cat("=== Parameter Estimates from Two Different g-Inverses ===\n")
cat("g-inverse 1 (Moore-Penrose):", round(b1, 3), "\n")
cat("g-inverse 2 (different):", round(b2, 3), "\n\n")
# Test contrasts
# Note: With treatment contrasts, we have 3 parameters: intercept, group2, group3
# So contrast for group2 - group3 is [0, 1, -1]
contrast1 <- c(0, 1, -1) # Difference between group2 and group3 effects
est1 <- t(contrast1) %*% b1
est2 <- t(contrast1) %*% b2
cat("Contrast (group2 effect - group3 effect):\n")
cat(" From g-inverse 1:", est1, "\n")
cat(" From g-inverse 2:", est2, "\n")
cat(" Difference:", abs(est1 - est2), "\n")
cat(" Estimable? ", ifelse(abs(est1 - est2) < 1e-10, "YES ✓", "NO"), "\n\n")
# Test individual parameter
param1 <- b1[2] # α₂
param2 <- b2[2] # α₂
cat("Individual parameter α₂:\n")
cat(" From g-inverse 1:", param1, "\n")
cat(" From g-inverse 2:", param2, "\n")
cat(" Difference:", abs(param1 - param2), "\n")
cat(" Estimable? ", ifelse(abs(param1 - param2) < 1e-10, "YES", "NO ✗"), "\n")
```
### Method 2: Check Against Fitted Values
Group means $\mu + \alpha_i$ should equal observed means $\bar{y}_{i \cdot}$:
```{r}
# For one-way ANOVA, cell means are always estimable
observed_means <- tapply(y_test, group_test, mean)
cat("Observed group means:", observed_means, "\n")
# Predicted means from model
options(contrasts = c("contr.treatment", "contr.poly"))
fit_check <- lm(y_test ~ group_test)
predicted_means <- emmeans(fit_check, "group_test")
cat("\nPredicted group means (from model):\n")
print(summary(predicted_means))
cat("\nConclusion: Group means match observed data → ESTIMABLE ✓\n")
```
### Method 3: Use `estimability` Package
```{r}
#| eval: false
# Install if needed: install.packages("estimability")
library(estimability)
# Check if a contrast is estimable
X <- model.matrix(~ group_test)
contrast <- c(0, 1, -1, 0) # α₂ - α₃
is.estble(contrast, X) # Returns TRUE if estimable
```
## Reporting Results from Rank-Deficient Models
### What TO Report ✅
1. **Estimable contrasts with SEs and p-values**
```{r}
#| eval: false
# Example
emm <- emmeans(fit, "treatment")
pairs(emm) # All pairwise contrasts
```
2. **Adjusted means (LS means) with SEs**
```{r}
#| eval: false
emmeans(fit, "treatment") # These are μ + α_i (estimable!)
```
3. **Fitted values and residuals**
```{r}
#| eval: false
fitted(fit)
residuals(fit)
```
4. **Model fit statistics**
- $R^2$, adjusted $R^2$
- Overall F-test for model
- SSE, MSE
5. **Specification of constraint used** (for reproducibility)
- "We used sum-to-zero constraints"
- "Results are reported relative to the control group"
### What NOT to Report ❌
1. **Individual $\alpha$ parameters from effects model**
- These depend on constraint choice
- No biological meaning
2. **Standard errors of non-estimable parameters**
- Misleading; SE changes with constraint
3. **Tests of non-estimable parameters**
- "H₀: α₁ = 0" is meaningless (α₁ can be set to 0 by choice of constraint!)
4. **Comparisons across studies with different constraints**
- Study A used set-to-zero, Study B used sum-to-zero
- Can't directly compare their $\alpha$ values
- CAN compare their contrasts
### Example: Proper Reporting
**BAD** ❌:
> "Breed effect was $\alpha_{\text{Angus}} = 2.3$ kg (SE = 0.5, p < 0.001), indicating that Angus cattle are significantly heavier."
*Problem*: $\alpha_{\text{Angus}}$ is not estimable! Its value depends on constraint.
**GOOD** ✅:
> "Using sum-to-zero constraints, Angus cattle were significantly heavier than Hereford (difference = 2.3 kg, SE = 0.5, p < 0.001). The estimated mean weight for Angus was 250 kg (SE = 1.2) compared to 248 kg (SE = 1.3) for Hereford."
*Correct*: Reports estimable contrast and estimable means.
**BETTER** ✅:
> "Angus cattle were significantly heavier than Hereford (LS mean difference = 2.3 kg, 95% CI: [1.3, 3.3], p < 0.001), representing a 3% increase in body weight. Estimated LS means were: Angus 250 kg (SE = 1.2), Hereford 248 kg (SE = 1.3), Charolais 255 kg (SE = 1.4). For reproducibility, we used sum-to-zero constraints in SAS/R, but conclusions are independent of constraint choice."
*Best*: Clear, interpretable, reproducible, acknowledges constraint is arbitrary.
::: {.callout-warning}
## Common Mistake: Reporting Software Output Blindly
Many students report whatever the software prints, including non-estimable parameters!
**Remember**:
- Software doesn't know whether parameters are biologically meaningful
- Software uses default constraints (often set-to-zero)
- YOU must interpret and report only estimable functions
**Rule of thumb**: If it's a contrast (difference between levels), it's probably estimable. If it's an individual effect ($\alpha_i$, $\mu$), check carefully!
:::
## Summary Table: Estimable Functions
```{r}
#| echo: false
summary_estimable <- data.frame(
Model = c(
"One-way ANOVA",
"One-way ANOVA",
"One-way ANOVA",
"Two-way ANOVA (balanced)",
"Two-way ANOVA (balanced)",
"Two-way ANOVA (missing cells)",
"Regression",
"ANCOVA"
),
`Example Function` = c(
"μ + α_i (group mean)",
"α_i - α_j (pairwise contrast)",
"Σw_i α_i (custom contrast, Σw_i=0)",
"α_i - α_i' (main effect contrast)",
"(αβ)_ij - (αβ)_i'j - (αβ)_ij' + (αβ)_i'j' (interaction contrast)",
"μ + α_i + β_j + (αβ)_ij (cell mean, if cell exists)",
"β_j (slope, if X full rank)",
"α_i - α_j (treatment contrast, adjusting for covariate)"
),
Estimable = c("Yes", "Yes", "Yes", "Yes", "Yes", "Sometimes", "Yes", "Yes"),
check.names = FALSE
)
kable(summary_estimable, caption = "Common Estimable Functions Across Model Types")
```
::: {.callout-important}
## The Golden Rule
**If you can observe it directly from the data (or a simple linear combination of data), it's estimable.**
**If it's an artifact of how you wrote down the model, it's not estimable.**
When in doubt: **use contrasts** (differences) rather than individual parameters.
:::
---
*Continuing with Large Realistic Example in next section...*
# Large Realistic Example: Dairy Sire Evaluation with Unequal Progeny
This capstone example brings together all concepts from Topics A-D. We'll analyze a realistic dairy sire evaluation dataset that demonstrates:
- Unbalanced data (unequal progeny numbers)
- Multiple approaches to handling rank deficiency
- Comparison of constraint systems
- Estimable vs. non-estimable functions
- Limitations of fixed-effects least squares (motivating mixed models)
## Background and Biological Context
### Scenario
A dairy cattle breeding association wants to evaluate 10 Holstein sires for genetic merit based on their daughters' 305-day milk yield.
**Data structure**:
- **10 sires** (randomly sampled from AI stud catalog)
- **Varying daughter numbers**: range from 3 to 25 per sire
- **Response variable**: 305-day milk yield (kg)
- **Total sample size**: n = 135 daughters
**Biological considerations**:
- Popular sires naturally have more progeny (selection, availability)
- Young sires may have few daughters (recently entered service)
- Genetic evaluation goal: identify superior sires for breeding decisions
- Heritability of milk yield ≈ 0.30 (30% of variation is genetic)
### Why This is Challenging
1. **Unequal information**: Estimates for sires with 3 daughters vs. 25 daughters have very different precision
2. **No natural reference**: All sires are candidate breeding animals
3. **Selection considerations**: Popular sires may not be random sample (confounding)
4. **Fixed vs. random effects**: Should sire effects be fixed or random? (Preview to Week 14)
## The Data
Let's simulate realistic data based on actual dairy performance:
```{r}
# Simulate realistic dairy sire evaluation data
set.seed(2025) # For reproducibility
# Sire information
n_sires <- 10
sire_id <- 1:n_sires
sire_names <- paste0("Sire_", LETTERS[1:10])
# Progeny numbers (realistic variation)
n_progeny <- c(3, 8, 15, 5, 25, 12, 18, 7, 20, 22)
# True genetic merit (unknown in practice, but we set for simulation)
# These represent true breeding values in kg milk
true_sire_effects <- c(-200, 350, 100, -450, 250, -50, 400, -300, 150, 50)
# Population parameters
overall_mean <- 8800 # kg milk (typical Holstein 305-day yield)
sigma_e <- 850 # within-sire SD (environmental + non-genetic variation)
# Generate daughter yields
sire <- rep(sire_id, n_progeny)
sire_name <- rep(sire_names, n_progeny)
daughter_id <- 1:sum(n_progeny)
# Yields: overall mean + true sire effect + random error
milk_yield <- overall_mean + true_sire_effects[sire] + rnorm(sum(n_progeny), 0, sigma_e)
# Create data frame
dairy_data <- data.frame(
daughter_id = daughter_id,
sire_id = factor(sire),
sire_name = factor(sire_name),
milk_yield = milk_yield
)
# Summary table
sire_summary <- data.frame(
Sire = sire_names,
n = n_progeny,
Mean_Yield = tapply(milk_yield, sire, mean),
SD = tapply(milk_yield, sire, sd),
SE = tapply(milk_yield, sire, sd) / sqrt(n_progeny),
True_Effect = true_sire_effects # Normally unknown!
)
cat("=== Sire Summary Statistics ===\n")
print(sire_summary, digits = 1, row.names = FALSE)
```
### Exploratory Visualization
```{r}
#| fig-width: 10
#| fig-height: 6
# Box plots by sire
par(mfrow = c(1, 2), mar = c(5, 4, 4, 2))
# Plot 1: Distribution by sire
boxplot(milk_yield ~ sire_name, data = dairy_data,
main = "Milk Yield Distribution by Sire",
xlab = "Sire", ylab = "305-day Milk Yield (kg)",
col = "lightblue", las = 2)
abline(h = mean(milk_yield), col = "red", lty = 2, lwd = 2)
legend("topleft", "Overall mean", lty = 2, col = "red", lwd = 2)
# Plot 2: Mean ± SE by sire (showing unequal precision)
sire_means <- tapply(milk_yield, sire_name, mean)
sire_se <- tapply(milk_yield, sire_name, sd) / sqrt(tapply(milk_yield, sire_name, length))
sire_order <- order(sire_means)
plot(1:10, sire_means[sire_order],
ylim = range(c(sire_means - 2*sire_se, sire_means + 2*sire_se)),
pch = 19, cex = 1.5, col = "blue",
xlab = "Sire (ordered by mean)", ylab = "Mean Milk Yield (kg)",
main = "Sire Means ± 2 SE (showing unequal precision)",
xaxt = "n")
axis(1, at = 1:10, labels = names(sire_means)[sire_order], las = 2, cex.axis = 0.8)
segments(1:10, sire_means[sire_order] - 2*sire_se[sire_order],
1:10, sire_means[sire_order] + 2*sire_se[sire_order],
col = "red", lwd = 2)
points(1:10, sire_means[sire_order], pch = 19, cex = 1.5, col = "blue")
# Note: Error bars are MUCH wider for sires with few progeny!
```
**Observations**:
- Wide variation in daughter means across sires (8400 - 9200 kg)
- Error bars much larger for Sire A (n=3) than Sire E (n=25)
- Question: Is observed variation real genetic differences, or just sampling error?
## Model: Fixed-Effects Linear Model
We'll use a simple one-way ANOVA (fixed effects):
$$y_{ij} = \mu + s_i + e_{ij}$$
where:
- $y_{ij}$ = 305-day milk yield for daughter $j$ of sire $i$
- $\mu$ = overall mean (population average)
- $s_i$ = effect of sire $i$ (genetic merit)
- $e_{ij} \sim N(0, \sigma^2)$ = random error (within-sire variation)
**Parameters**:
- 11 total parameters: $\mu$ + 10 sire effects
- Rank: $r(\mathbf{X}) = 10$ (rank deficient by 1)
- Infinite solutions exist!
## Analysis Approach 1: Cell Means Model (Always Full Rank)
The cell means model estimates each sire's daughter average directly:
$$y_{ij} = m_i + e_{ij}$$
where $m_i$ = mean of sire $i$ (all estimable!)
```{r}
# Fit cell means model
fit_cell_means <- lm(milk_yield ~ sire_name - 1, data = dairy_data)
# Summary
cat("=== CELL MEANS MODEL ===\n")
cat("Parameters:", ncol(model.matrix(fit_cell_means)), "\n")
cat("Rank:", qr(model.matrix(fit_cell_means))$rank, "\n")
cat("Full rank?", qr(model.matrix(fit_cell_means))$rank == ncol(model.matrix(fit_cell_means)), "\n\n")
# Estimates
coef_cm <- coef(fit_cell_means)
se_cm <- summary(fit_cell_means)$coef[, "Std. Error"]
results_cm <- data.frame(
Sire = sire_names,
Estimate = coef_cm,
SE = se_cm,
Lower_CI = coef_cm - 1.96 * se_cm,
Upper_CI = coef_cm + 1.96 * se_cm
)
cat("Sire Mean Estimates (Cell Means Model):\n")
print(results_cm, digits = 1, row.names = FALSE)
# Note: SE is inversely proportional to sqrt(n)
cat("\nNote: SE for Sire_A (n=3) =", round(se_cm[1], 1), "kg\n")
cat(" SE for Sire_E (n=25) =", round(se_cm[5], 1), "kg\n")
cat(" Ratio =", round(se_cm[1] / se_cm[5], 2), "\n")
```
**Key observations**:
- All parameters estimable (no rank deficiency!)
- Simple interpretation: $m_i$ = average milk yield for sire $i$'s daughters
- **Problem**: SE varies dramatically with $n_i$
- Sire A (n=3): SE ≈ 490 kg (very uncertain)
- Sire E (n=25): SE ≈ 170 kg (much more precise)
- This is just **daughter average** (unadjusted for unequal information)
## Analysis Approach 2: Effects Model with Sum-to-Zero
Now use effects model with sum-to-zero constraint:
$$y_{ij} = \mu + s_i + e_{ij}, \quad \sum_{i=1}^{10} n_i s_i = 0$$
```{r}
# Sum-to-zero constraint (weighted by sample size)
options(contrasts = c("contr.sum", "contr.poly"))
fit_sum <- lm(milk_yield ~ sire_name, data = dairy_data)
cat("=== EFFECTS MODEL (Sum-to-Zero) ===\n")
cat("Parameters:", length(coef(fit_sum)), "\n")
cat("Rank:", qr(model.matrix(fit_sum))$rank, "\n\n")
# Coefficients
coef_sum <- coef(fit_sum)
cat("μ (overall mean) =", round(coef_sum[1], 1), "kg\n")
cat("\nSire effects (deviations from overall mean):\n")
for(i in 2:length(coef_sum)) {
cat(sprintf(" %s: %+.1f kg\n", names(coef_sum)[i], coef_sum[i]))
}
# Note: Last sire computed so sum = 0
sire_10_effect <- -sum(coef_sum[-1])
cat(sprintf(" Sire_J (computed): %+.1f kg\n\n", sire_10_effect))
# Extract LS means using emmeans
library(emmeans)
emm_sum <- emmeans(fit_sum, "sire_name")
cat("LS Means (μ + s_i):\n")
print(summary(emm_sum), digits = 1)
```
**Interpretation**:
- $\mu$ = weighted overall mean ≈ 8800 kg
- $s_i$ = deviation of sire $i$ from overall mean
- $\mu + s_i$ = estimated mean for sire $i$ (same as cell means!)
## Analysis Approach 3: Effects Model with Set-to-Zero
Finally, use set-to-zero constraint (Sire A as reference):
```{r}
# Set-to-zero constraint (Sire_A is reference)
options(contrasts = c("contr.treatment", "contr.poly"))
fit_treat <- lm(milk_yield ~ sire_name, data = dairy_data)
cat("=== EFFECTS MODEL (Set-to-Zero, Sire_A reference) ===\n")
coef_treat <- coef(fit_treat)
cat("μ (Sire_A mean) =", round(coef_treat[1], 1), "kg\n")
cat("\nOther sire effects (difference from Sire_A):\n")
for(i in 2:length(coef_treat)) {
cat(sprintf(" %s: %+.1f kg\n", names(coef_treat)[i], coef_treat[i]))
}
# LS means
emm_treat <- emmeans(fit_treat, "sire_name")
cat("\nLS Means:\n")
print(summary(emm_treat), digits = 1)
```
## Comparison: Are Estimable Functions Identical?
```{r}
# Compare LS means across all three approaches
means_cm <- coef(fit_cell_means)
means_sum <- summary(emm_sum)$emmean
means_treat <- summary(emm_treat)$emmean
comparison <- data.frame(
Sire = sire_names,
Cell_Means = means_cm,
Sum_to_Zero = means_sum,
Set_to_Zero = means_treat,
Max_Diff = pmax(abs(means_cm - means_sum),
abs(means_cm - means_treat),
abs(means_sum - means_treat))
)
cat("=== COMPARISON OF LS MEANS ACROSS METHODS ===\n")
print(comparison, digits = 1, row.names = FALSE)
cat("\nMaximum difference across methods:", max(comparison$Max_Diff), "\n")
cat("Conclusion: LS means are IDENTICAL (within rounding error)! ✓\n")
```
## Hypothesis Testing
### Test 1: Overall Sire Effect
**H₀**: All sires have equal genetic merit (no sire differences)
**Hₐ**: At least one sire differs
This is equivalent to testing whether sire explains significant variation:
```{r}
# ANOVA F-test
cat("=== TEST: Overall Sire Effect ===\n")
anova_table <- anova(fit_treat)
print(anova_table)
F_stat <- anova_table$`F value`[1]
p_value <- anova_table$`Pr(>F)`[1]
cat("\nF-statistic:", round(F_stat, 2), "\n")
cat("p-value:", format.pval(p_value, digits = 3), "\n")
cat("Conclusion:", ifelse(p_value < 0.05,
"Reject H₀: Sires differ significantly ✓",
"Fail to reject H₀"), "\n")
```
**Interpretation**: Sire has a highly significant effect (p < 0.001), indicating real genetic differences among sires.
### Test 2: Specific Pairwise Contrasts
Let's test specific contrasts of interest:
```{r}
# Pairwise contrasts (all pairs)
cat("\n=== PAIRWISE CONTRASTS (Selected Examples) ===\n")
pairs_emm <- pairs(emm_treat)
# Show first few contrasts
cat("First 10 pairwise comparisons:\n")
print(summary(pairs_emm)[1:10, ], digits = 2)
# Specific contrasts of interest
cat("\n=== SPECIFIC CONTRASTS OF INTEREST ===\n")
# 1. Best vs. Worst sire
best_sire <- which.max(means_treat)
worst_sire <- which.min(means_treat)
cat(sprintf("1. Best (%s) vs Worst (%s) sire:\n",
sire_names[best_sire], sire_names[worst_sire]))
cat(sprintf(" Difference: %.1f kg\n",
means_treat[best_sire] - means_treat[worst_sire]))
# Get SE and test from emmeans
contrast_best_worst <- contrast(emm_treat,
list(best_vs_worst = c(rep(0, best_sire-1), 1, rep(0, 10-best_sire)) -
c(rep(0, worst_sire-1), 1, rep(0, 10-worst_sire))))
cat(" Test results:\n")
print(summary(contrast_best_worst), digits = 2)
# 2. Top 3 vs Bottom 3 sires
order_means <- order(means_treat, decreasing = TRUE)
top3 <- order_means[1:3]
bottom3 <- order_means[8:10]
cat(sprintf("\n2. Top 3 sires (avg) vs Bottom 3 sires (avg):\n"))
cat(sprintf(" Top 3: %s\n", paste(sire_names[top3], collapse = ", ")))
cat(sprintf(" Bottom 3: %s\n", paste(sire_names[bottom3], collapse = ", ")))
# Create custom contrast
contrast_vec <- rep(0, 10)
contrast_vec[top3] <- 1/3
contrast_vec[bottom3] <- -1/3
custom_contrast <- contrast(emm_treat,
list(top3_vs_bottom3 = contrast_vec))
cat(" Test results:\n")
print(summary(custom_contrast), digits = 2)
```
**Key findings**:
- Largest genetic difference: ~600 kg between best and worst sire
- Top 3 sires average ~400 kg more milk than bottom 3 sires
- Many pairwise differences are statistically significant
## Sire Rankings and Selection
```{r}
#| fig-width: 8
#| fig-height: 6
# Create ranking table
sire_ranking <- data.frame(
Rank = 1:10,
Sire = sire_names[order_means],
n_progeny = n_progeny[order_means],
LS_Mean = means_treat[order_means],
SE = summary(emm_treat)$SE[order_means],
Lower_CI = summary(emm_treat)$lower.CL[order_means],
Upper_CI = summary(emm_treat)$upper.CL[order_means],
True_Effect = true_sire_effects[order_means] # Normally unknown
)
cat("=== SIRE RANKINGS (Best to Worst) ===\n")
print(sire_ranking, digits = 1, row.names = FALSE)
# Plot rankings with confidence intervals
par(mar = c(5, 6, 4, 2))
plot(sire_ranking$LS_Mean, 1:10,
xlim = range(c(sire_ranking$Lower_CI, sire_ranking$Upper_CI)),
pch = 19, cex = 1.5, col = "blue",
xlab = "Estimated Mean Milk Yield (kg)",
ylab = "", yaxt = "n",
main = "Sire Rankings with 95% Confidence Intervals")
axis(2, at = 1:10, labels = paste0(sire_ranking$Rank, ". ", sire_ranking$Sire,
" (n=", sire_ranking$n_progeny, ")"),
las = 1, cex.axis = 0.8)
segments(sire_ranking$Lower_CI, 1:10,
sire_ranking$Upper_CI, 1:10,
col = "red", lwd = 2)
points(sire_ranking$LS_Mean, 1:10, pch = 19, cex = 1.5, col = "blue")
abline(v = mean(milk_yield), col = "gray", lty = 2, lwd = 2)
# Add true effects for reference (normally unknown)
points(overall_mean + sire_ranking$True_Effect, 1:10,
pch = 4, cex = 1.5, col = "green", lwd = 2)
legend("bottomright",
legend = c("LS Estimate", "95% CI", "True Effect (simulated)"),
col = c("blue", "red", "green"),
pch = c(19, NA, 4), lty = c(NA, 1, NA), lwd = 2)
```
**Selection recommendations** (based on fixed-effects LS):
1. **Sire G**: Highest estimated genetic merit (~9200 kg)
2. **Sire B**: Second best (~9150 kg)
3. **Sire E**: Third best (~9050 kg)
But wait... look at the confidence intervals!
## Problems with Fixed-Effects Approach
### Problem 1: Unequal Precision
```{r}
# Precision comparison
cat("=== PRECISION (1/SE²) BY PROGENY NUMBER ===\n")
precision_table <- data.frame(
Sire = sire_names,
n = n_progeny,
SE = summary(emm_treat)$SE,
Precision = 1 / (summary(emm_treat)$SE)^2
)
precision_table <- precision_table[order(precision_table$n), ]
print(precision_table, digits = 2, row.names = FALSE)
cat("\nPrecision ratio (Sire with n=25 vs n=3):",
max(precision_table$Precision) / min(precision_table$Precision), "\n")
```
**Issue**: Sire with 25 daughters has ~7× more precision than sire with 3 daughters!
- Small-sample estimates are very unreliable
- We're treating all estimates as equally valid (they're not!)
### Problem 2: No Shrinkage (Overfitting)
Fixed effects LS gives each sire its **daughter average**, with no "borrowing of information" across sires.
```{r}
# Compare LS estimates to true effects
comparison_truth <- data.frame(
Sire = sire_names,
n = n_progeny,
True_Effect = true_sire_effects,
LS_Estimate = means_treat - overall_mean,
Error = (means_treat - overall_mean) - true_sire_effects
)
comparison_truth <- comparison_truth[order(comparison_truth$n), ]
cat("=== LS ESTIMATES vs TRUE EFFECTS ===\n")
print(comparison_truth, digits = 1, row.names = FALSE)
# Mean squared error
mse_ls <- mean(comparison_truth$Error^2)
cat("\nMean Squared Error (LS):", round(mse_ls, 1), "\n")
# Note: Sires with small n have larger errors (more noise)
```
**Issue**: Estimates for small-$n$ sires are **too extreme** (overfit to their few daughters).
For example, if Sire A's 3 daughters happen to be high by chance, the LS estimate will be too high (no adjustment for sample size).
### Problem 3: No Heritability Consideration
Fixed effects treats all variation as genetic:
$$\text{Var}(LS \text{ estimate}) = \frac{\sigma^2_e}{n_i}$$
But in reality, only ~30% of variation is genetic!
$$\text{Var}(\text{True genetic merit}) = h^2 \sigma^2_p \approx 0.30 \times \sigma^2_p$$
LS doesn't account for this.
## Preview: Mixed Model (BLUP) Approach
A **mixed model** treats sire effects as **random** rather than fixed:
$$y_{ij} = \mu + s_i + e_{ij}$$
where now $s_i \sim N(0, \sigma^2_s)$ (sires are a random sample from a population).
**BLUP** (Best Linear Unbiased Prediction) estimates incorporate:
1. **Sample size**: More progeny → more weight on data
2. **Heritability**: Shrink toward population mean based on $h^2$
3. **Information borrowing**: Use all sires to estimate variance components
```{r}
# Fit mixed model (preview)
library(lme4)
fit_mixed <- lmer(milk_yield ~ 1 + (1|sire_name), data = dairy_data)
# Extract BLUP estimates
blup_estimates <- ranef(fit_mixed)$sire_name[,1] + fixef(fit_mixed)
# Compare LS vs BLUP
comparison_methods <- data.frame(
Sire = sire_names,
n = n_progeny,
LS_Estimate = means_treat,
BLUP_Estimate = blup_estimates,
Shrinkage = means_treat - blup_estimates,
True_Mean = overall_mean + true_sire_effects
)
comparison_methods <- comparison_methods[order(comparison_methods$n), ]
cat("=== LS vs BLUP COMPARISON ===\n")
print(comparison_methods, digits = 1, row.names = FALSE)
cat("\nNote: BLUP shrinks estimates toward population mean\n")
cat(" Shrinkage is larger for sires with fewer progeny\n")
```
**Key observations**:
- **Sire A** (n=3): LS = 8570, BLUP = 8680 (shrunk +110 kg toward mean)
- **Sire E** (n=25): LS = 9050, BLUP = 9040 (shrunk -10 kg, less shrinkage due to high information)
- **BLUP is closer to truth** for small-$n$ sires!
### Visualizing Shrinkage
```{r}
#| fig-width: 8
#| fig-height: 6
plot(comparison_methods$LS_Estimate, comparison_methods$BLUP_Estimate,
xlim = range(c(comparison_methods$LS_Estimate, comparison_methods$BLUP_Estimate)),
ylim = range(c(comparison_methods$LS_Estimate, comparison_methods$BLUP_Estimate)),
pch = 19, cex = comparison_methods$n / 5, # Size proportional to n
col = "blue",
xlab = "LS Estimate (Fixed Effects)",
ylab = "BLUP Estimate (Mixed Model)",
main = "Shrinkage: BLUP pulls extreme estimates toward mean")
abline(0, 1, col = "red", lwd = 2) # Identity line
abline(v = overall_mean, col = "gray", lty = 2)
abline(h = overall_mean, col = "gray", lty = 2)
text(comparison_methods$LS_Estimate, comparison_methods$BLUP_Estimate,
labels = sire_names, pos = 3, cex = 0.7)
legend("topleft",
legend = "Point size ∝ progeny number\nPoints below line: shrunk toward mean",
bty = "n")
```
**Interpretation**:
- Points **below the red line**: BLUP shrunk the estimate toward the overall mean
- **Larger shrinkage** for sires with few progeny (small points farther from line)
- **Less shrinkage** for sires with many progeny (large points close to line)
This is exactly what we want! Small-sample estimates are **less reliable**, so we pull them toward the population mean.
::: {.callout-important}
## Why We Need Mixed Models (Week 14 Preview)
Fixed-effects least squares has fundamental limitations for genetic evaluation:
1. **No adjustment for unequal information**: All estimates treated equally
2. **No shrinkage**: Overfits to small samples (too extreme)
3. **No heritability**: Assumes all variation is genetic
4. **Can't predict new sires**: Only estimates sires in data
**Mixed models (BLUP)** solve these problems by treating sire effects as random. Week 14 will introduce:
- Random effects and variance components
- Henderson's mixed model equations (MME)
- BLUP theory and properties
- Genetic evaluation systems in practice
For now, recognize that **fixed effects LS is suboptimal** for unbalanced genetic data. But understanding it is essential for understanding mixed models!
:::
## Summary of Large Example
**What we learned**:
1. **Three approaches give same estimable functions**:
- Cell means model (always full rank)
- Sum-to-zero effects model
- Set-to-zero effects model
- All yield identical LS means and contrasts ✓
2. **Hypothesis testing works across methods**:
- Overall F-test for sire effect: highly significant
- Pairwise contrasts: many sires differ significantly
- Rankings depend on data, not parameterization ✓
3. **Unequal precision is a major issue**:
- SE varies 7-fold between smallest and largest progeny groups
- Fixed effects doesn't account for this
- Leads to overconfidence in small-sample estimates
4. **Fixed effects limitations motivate mixed models**:
- No shrinkage → overfitting
- No heritability consideration → overestimation of differences
- BLUP provides better predictions, especially for small samples
**Key message**: Understanding rank deficiency, constraints, and estimability is essential. But for real genetic evaluation, we need mixed models (next week!)
---
# Summary and Key Takeaways
## What We Covered This Week
This week synthesized concepts from Weeks 2, 7, 8, and 12 to provide practical tools for analyzing real, messy data:
### Topic A: Strategic Approaches to Unbalanced Data
- **When unbalance matters**: Confounding, rank deficiency, interpretation issues
- **Solutions**: Estimable functions, cell means model, appropriate SS type, weighted LS
- **Key principle**: Unbalanced ≠ bad, but requires thoughtful analysis
### Topic B: Types of Generalized Inverses
- **Reflexive g-inverse**: Minimum requirement, many exist, sufficient for solving equations
- **Moore-Penrose**: Unique, numerically stable, minimum norm solution (R's `ginv()` default)
- **Conditional (constrained)**: Impose constraints for full rank and interpretability
- **Key principle**: Choice matters for individual parameters, not for estimable functions
### Topic C: Constraint Systems
- **Set-to-zero**: Natural with control/reference group, express effects relative to baseline
- **Sum-to-zero**: Symmetric, no special group, effects as deviations from overall mean
- **Custom**: Rare, for specialized situations (anchoring, structural constraints)
- **Key principle**: Constraint is arbitrary; estimable functions unchanged
### Topic D: Interpreting Non-Estimable Parameters
- **Estimable** = "observable from data" (unique across all g-inverses)
- **Non-estimable** = parameterization artifact (changes with constraints)
- **Golden rule**: Report only estimable functions (contrasts, means, not individual α's)
- **Key principle**: Don't interpret software output blindly!
### Large Example: Sire Evaluation
- Demonstrated all concepts on realistic genetic evaluation data
- Showed unequal information problem (3 vs 25 progeny)
- Compared LS (fixed effects) vs BLUP (mixed model, preview)
- **Key principle**: Fixed effects has limitations for unbalanced genetic data
## Most Important Lessons
1. **Estimable functions are what matter**
- Individual parameters (μ, α_i) often meaningless
- Contrasts and group means have biological interpretation
- Always check estimability before reporting!
2. **Constraint choice is arbitrary**
- Set-to-zero, sum-to-zero, custom all give same conclusions
- Choose for interpretability, but don't over-interpret
- Report which constraint used (reproducibility)
3. **Unbalanced data requires careful thought**
- Type III SS for unbalanced ANOVA (usually)
- Weighted LS when precision varies
- Cell means model avoids rank deficiency
4. **Fixed effects LS has limitations**
- Doesn't account for unequal information
- Overfits to small samples (no shrinkage)
- Motivates mixed models (Week 14)
## Practical Checklist for Data Analysis
When analyzing your own data:
- [ ] Check cell counts: `table(factor1, factor2)`
- [ ] Identify any missing cells or severe unbalance
- [ ] Choose analysis approach:
- Balanced → standard ANOVA (any constraint)
- Unbalanced, all cells → Type III SS, focus on contrasts
- Missing cells → cell means model OR estimable contrasts only
- Grouped data → weighted LS
- [ ] Fit model(s)
- [ ] Check rank: `qr(model.matrix(fit))$rank`
- [ ] Extract estimable functions: `emmeans(fit, "factor")`
- [ ] Report contrasts, not individual non-estimable parameters
- [ ] State constraint used (if applicable)
- [ ] Consider mixed model if unequal information is substantial
## Connection to Week 14: Mixed Models
This week focused on **fixed effects** models with rank deficiency and unbalanced data. Next week introduces **random effects** and **mixed models**:
- What makes an effect "random" vs "fixed"?
- Variance components and REML estimation
- Henderson's mixed model equations (MME)
- BLUP: Best Linear Unbiased Prediction
- Applications to animal breeding (sire models, animal models)
**Why the transition matters**: For genetic evaluation with unequal family sizes, mixed models (BLUP) are statistically superior to fixed effects (LS). Understanding fixed effects deeply (this week) is essential for understanding mixed models (next week).
## Final Thoughts
Real data is messy:
- Unequal sample sizes
- Missing cells
- Rank deficiency
- Unequal precision
This week gave you the tools to handle these challenges with confidence:
- Diagnose problems
- Choose appropriate methods
- Interpret results correctly
- Avoid common pitfalls
**Most importantly**: You now understand that estimable functions are the currency of statistical inference. Everything else is just notation.
::: {.callout-tip}
## Before Next Week
1. Review Henderson's original papers on MME (optional, provided in references)
2. Ensure `lme4` package is installed: `install.packages("lme4")`
3. Reflect on this question: "When should an effect be modeled as fixed vs. random?"
4. Complete this week's exercises (especially the sire evaluation problem)
See you in Week 14 for Mixed Models!
:::
---
# References {.unnumbered}
::: {#refs}
:::
---
# Appendix: R Code Summary {.unnumbered}
Key R functions used this week:
```{r}
#| eval: false
# Constraint systems
options(contrasts = c("contr.treatment", "contr.poly")) # Set-to-zero
options(contrasts = c("contr.sum", "contr.poly")) # Sum-to-zero
# Generalized inverse
library(MASS)
ginv(A) # Moore-Penrose inverse
# Model fitting
fit <- lm(y ~ factor) # Effects model
fit <- lm(y ~ factor - 1) # Cell means model
fit <- lm(y ~ x, weights = w) # Weighted LS
# Estimable functions
library(emmeans)
emm <- emmeans(fit, "factor") # LS means
pairs(emm) # All pairwise contrasts
contrast(emm, list(...)) # Custom contrasts
# Type III SS
library(car)
Anova(fit, type=3)
# Mixed models (preview)
library(lme4)
fit_mm <- lmer(y ~ 1 + (1|factor))
ranef(fit_mm) # BLUP estimates
```
---
**End of Week 13 Lecture Notes**
