In Weeks 1 and 2, we established the computational foundations and reviewed the linear algebra essentials needed for linear models. Now we address the critical question: How do we get from raw data to the matrix equations we solve?
The design matrix (denoted \(\mathbf{X}\)) is the bridge between:
Understanding how to construct \(\mathbf{X}\) is fundamental to applying linear models in animal breeding and genetics.
In animal breeding applications, we constantly work with different types of predictors:
The way we code these predictors into the design matrix determines:
The general linear model expresses the relationship between observations and parameters as:
\[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{e} \]
where:
Always verify matrix dimensions for compatibility:
The model equation is dimensionally consistent.
Taking expected values of both sides:
\[ E(\mathbf{y}) = E(\mathbf{X}\boldsymbol{\beta} + \mathbf{e}) = \mathbf{X}\boldsymbol{\beta} + E(\mathbf{e}) = \mathbf{X}\boldsymbol{\beta} \]
This assumes \(E(\mathbf{e}) = \mathbf{0}\) (errors have mean zero, discussed below).
For the \(i\)-th observation:
\[ y_i = \sum_{j=1}^{p} x_{ij}\beta_j + e_i = x_{i1}\beta_1 + x_{i2}\beta_2 + \cdots + x_{ip}\beta_p + e_i \]
where \(x_{ij}\) is the element in row \(i\), column \(j\) of \(\mathbf{X}\).
The structure of \(\mathbf{X}\) depends on the type of predictors in your model. Let’s examine the three main types.
When all predictors are continuous variables, the design matrix includes:
Model: \(y_i = \beta_0 + \beta_1 x_i + e_i\)
For \(n=4\) observations:
# Data: broiler weight (kg) vs. age (days)
age <- c(21, 28, 35, 42)
weight <- c(0.5, 0.9, 1.4, 1.9)
# Design matrix: first column is intercept, second is predictor
X <- cbind(1, age)
print(X) age
[1,] 1 21
[2,] 1 28
[3,] 1 35
[4,] 1 42# Check dimensions
cat("\nDimensions: n =", nrow(X), "observations, p =", ncol(X), "parameters\n")
Dimensions: n = 4 observations, p = 2 parametersThe design matrix is:
\[ \mathbf{X} = \begin{bmatrix} 1 & 21 \\ 1 & 28 \\ 1 & 35 \\ 1 & 42 \end{bmatrix}_{4 \times 2} \]
The parameter vector is:
\[ \boldsymbol{\beta} = \begin{bmatrix} \beta_0 \\ \beta_1 \end{bmatrix}_{2 \times 1} \]
With multiple continuous predictors: \(y_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \cdots + \beta_k x_{ki} + e_i\)
# Example: Predict lamb weaning weight from birth weight and dam age
birth_wt <- c(4.5, 4.0, 4.8, 4.2, 4.6)
dam_age <- c(3, 5, 4, 6, 4)
wean_wt <- c(28, 24, 30, 26, 29)
# Design matrix: intercept, birth weight, dam age
X_mult <- cbind(1, birth_wt, dam_age)
colnames(X_mult) <- c("Intercept", "BirthWt", "DamAge")
print(X_mult) Intercept BirthWt DamAge
[1,] 1 4.5 3
[2,] 1 4.0 5
[3,] 1 4.8 4
[4,] 1 4.2 6
[5,] 1 4.6 4\[ \mathbf{X} = \begin{bmatrix} 1 & 4.5 & 3 \\ 1 & 4.0 & 5 \\ 1 & 4.8 & 4 \\ 1 & 4.2 & 6 \\ 1 & 4.6 & 4 \end{bmatrix}_{5 \times 3}, \quad \boldsymbol{\beta} = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \end{bmatrix}_{3 \times 1} \]
Categorical predictors (factors) require indicator variables (also called dummy variables). There are two primary coding schemes.
The cell means model estimates a separate mean for each group, with no explicit intercept.
Model: \(y_{ij} = \mu_i + e_{ij}\)
where \(\mu_i\) is the mean for group \(i\) (\(i = 1, \ldots, g\) groups).
The design matrix has one column per group, with indicators for group membership:
\[ x_{ij} = \begin{cases} 1 & \text{if observation } j \text{ is in group } i \\ 0 & \text{otherwise} \end{cases} \]
# Load the pig data
pig_data <- data.frame(
breed = c("Yorkshire", "Yorkshire", "Landrace", "Landrace", "Duroc", "Duroc"),
litter_size = c(11, 12, 10, 11, 9, 10)
)
print(pig_data) breed litter_size
1 Yorkshire 11
2 Yorkshire 12
3 Landrace 10
4 Landrace 11
5 Duroc 9
6 Duroc 10# Cell means design matrix: one column per breed
# Yorkshire = column 1, Landrace = column 2, Duroc = column 3
X_cell <- model.matrix(~ breed - 1, data = pig_data) # -1 removes intercept
print(X_cell) breedDuroc breedLandrace breedYorkshire
1 0 0 1
2 0 0 1
3 0 1 0
4 0 1 0
5 1 0 0
6 1 0 0
attr(,"assign")
[1] 1 1 1
attr(,"contrasts")
attr(,"contrasts")$breed
[1] "contr.treatment"# Check rank
cat("\nRank of X:", qr(X_cell)$rank, "\n")
Rank of X: 3 cat("Number of parameters:", ncol(X_cell), "\n")Number of parameters: 3 The cell means design matrix is:
\[ \mathbf{X}_{\text{cell}} = \begin{bmatrix} 1 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 1 \end{bmatrix}_{6 \times 3} \]
Parameter vector:
\[ \boldsymbol{\beta}_{\text{cell}} = \begin{bmatrix} \mu_{\text{Duroc}} \\ \mu_{\text{Landrace}} \\ \mu_{\text{Yorkshire}} \end{bmatrix}_{3 \times 1} \]
The effects model decomposes each observation into an overall mean plus group-specific deviations.
Model: \(y_{ij} = \mu + \alpha_i + e_{ij}\)
where:
Problem: This model is overparameterized without constraints.
With \(g\) groups, we have \(g+1\) parameters (\(\mu\) and \(\alpha_1, \ldots, \alpha_g\)), but only \(g\) distinct means. We need a constraint to make parameters identifiable.
Common constraint: \(\sum_{i=1}^{g} \alpha_i = 0\) (sum-to-zero constraint)
With this constraint, \(\alpha_i\) represents the deviation of group \(i\) from the overall mean.
# Effects model design matrix: intercept + group indicators
# R uses reference cell coding by default (first level = 0)
X_effects <- model.matrix(~ breed, data = pig_data)
print(X_effects) (Intercept) breedLandrace breedYorkshire
1 1 0 1
2 1 0 1
3 1 1 0
4 1 1 0
5 1 0 0
6 1 0 0
attr(,"assign")
[1] 0 1 1
attr(,"contrasts")
attr(,"contrasts")$breed
[1] "contr.treatment"# Check rank
cat("\nRank of X:", qr(X_effects)$rank, "\n")
Rank of X: 3 cat("Number of parameters:", ncol(X_effects), "\n")Number of parameters: 3 The effects model design matrix (without constraint applied yet):
\[ \mathbf{X}_{\text{effects}} = \begin{bmatrix} 1 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 0 \\ 1 & 0 & 1 \\ 1 & 0 & 1 \end{bmatrix}_{6 \times 3} \]
The full effects model matrix (with all \(g+1\) parameters) would be:
\[ \mathbf{X}_{\text{full}} = \begin{bmatrix} 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 \end{bmatrix}_{6 \times 4} \]
This has rank \(r(\mathbf{X}) = 3 < 4\), making it rank deficient. Individual parameters (\(\mu\), \(\alpha_i\)) are not uniquely estimable, but contrasts like \(\alpha_i - \alpha_j\) are estimable.
We’ll address non-full rank models in detail in Week 12.
R’s default for categorical variables uses reference cell coding (also called treatment coding):
With this constraint:
# Compute group means
group_means <- tapply(pig_data$litter_size, pig_data$breed, mean)
print(group_means) Duroc Landrace Yorkshire
9.5 10.5 11.5 # Interpretation with Duroc as reference (first alphabetically)
cat("\nWith reference cell coding (Duroc as baseline):\n")
With reference cell coding (Duroc as baseline):cat("μ (intercept) = Duroc mean =", group_means["Duroc"], "\n")μ (intercept) = Duroc mean = 9.5 cat("α_Landrace = Landrace - Duroc =", group_means["Landrace"] - group_means["Duroc"], "\n")α_Landrace = Landrace - Duroc = 1 cat("α_Yorkshire = Yorkshire - Duroc =", group_means["Yorkshire"] - group_means["Duroc"], "\n")α_Yorkshire = Yorkshire - Duroc = 2 When models include both categorical and continuous predictors, we have Analysis of Covariance (ANCOVA).
Model: \(y_{ij} = \mu + \alpha_i + \beta(x_{ij} - \bar{x}) + e_{ij}\)
The design matrix combines indicator columns (for groups) and continuous columns (for covariates).
# Data: egg production (eggs/month) by strain, with body weight covariate
egg_data <- data.frame(
strain = rep(c("Strain1", "Strain2", "Strain3"), each = 3),
body_weight = c(1.8, 2.0, 1.9, 1.7, 1.8, 1.9, 2.1, 2.3, 2.2),
eggs = c(24, 26, 25, 22, 23, 24, 26, 28, 27)
)
# ANCOVA design matrix: strain indicators + centered body weight
egg_data$bw_centered <- egg_data$body_weight - mean(egg_data$body_weight)
X_ancova <- model.matrix(~ strain + bw_centered, data = egg_data)
print(X_ancova) (Intercept) strainStrain2 strainStrain3 bw_centered
1 1 0 0 -0.16666667
2 1 0 0 0.03333333
3 1 0 0 -0.06666667
4 1 1 0 -0.26666667
5 1 1 0 -0.16666667
6 1 1 0 -0.06666667
7 1 0 1 0.13333333
8 1 0 1 0.33333333
9 1 0 1 0.23333333
attr(,"assign")
[1] 0 1 1 2
attr(,"contrasts")
attr(,"contrasts")$strain
[1] "contr.treatment"We often center continuous covariates by subtracting the mean: \(x_{ij}^* = x_{ij} - \bar{x}\)
Benefits:
Note: Centering does NOT change \(R^2\), residuals, or fitted values—only the interpretation of \(\beta_0\).
For the general linear model \(\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{e}\), we make specific assumptions about the error term \(\mathbf{e}\).
Under these assumptions, ordinary least squares (OLS) estimators are BLUE (Best Linear Unbiased Estimators):
The assumption \(\text{Var}(\mathbf{e}) = \sigma^2\mathbf{I}\) can be written explicitly:
\[ \text{Var}(\mathbf{e}) = \begin{bmatrix} \sigma^2 & 0 & 0 & \cdots & 0 \\ 0 & \sigma^2 & 0 & \cdots & 0 \\ 0 & 0 & \sigma^2 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & \sigma^2 \end{bmatrix} = \sigma^2 \mathbf{I} \]
For hypothesis testing and confidence intervals, we often add:
This implies \(\mathbf{y} \sim N(\mathbf{X}\boldsymbol{\beta}, \sigma^2\mathbf{I})\)
For large samples, the Central Limit Theorem provides approximate normality even if errors aren’t exactly normal.
| Violation | Consequence | Solution |
|---|---|---|
| Non-linearity | Biased estimates | Transform variables, add polynomial terms |
| \(E(\mathbf{e}) \neq \mathbf{0}\) | Biased intercept | Check model specification |
| Heteroscedasticity | Inefficient estimates, incorrect SEs | Weighted least squares, robust SEs |
| Correlation (e.g., time series) | Inefficient estimates, incorrect SEs | Generalized least squares, mixed models |
| Non-normality | Invalid inference (test statistics) | Transformations, bootstrapping, larger samples |
We’ll explore diagnostics to check these assumptions in Week 11.
Let’s work through a complete small example, showing both cell means and effects model formulations.
# Read the data
pig_data <- read.csv("data/pig_litter_breeds.csv")
print(pig_data) sow_id breed litter_size
1 1 Yorkshire 11
2 2 Yorkshire 12
3 3 Landrace 10
4 4 Landrace 11
5 5 Duroc 9
6 6 Duroc 10# Summary statistics by breed
cat("\nSummary by breed:\n")
Summary by breed:by(pig_data$litter_size, pig_data$breed, function(x) {
cat(sprintf(" n = %d, mean = %.2f, sd = %.2f\n", length(x), mean(x), sd(x)))
}) n = 2, mean = 9.50, sd = 0.71
n = 2, mean = 10.50, sd = 0.71
n = 2, mean = 11.50, sd = 0.71pig_data$breed: Duroc
NULL
------------------------------------------------------------
pig_data$breed: Landrace
NULL
------------------------------------------------------------
pig_data$breed: Yorkshire
NULLModel: \(y_{ij} = \mu_i + e_{ij}\) where \(i \in \{\text{Duroc, Landrace, Yorkshire}\}\)
# Construct design matrix manually
n <- nrow(pig_data)
X_cell <- matrix(0, nrow = n, ncol = 3)
X_cell[pig_data$breed == "Duroc", 1] <- 1
X_cell[pig_data$breed == "Landrace", 2] <- 1
X_cell[pig_data$breed == "Yorkshire", 3] <- 1
colnames(X_cell) <- c("Duroc", "Landrace", "Yorkshire")
cat("Cell Means Design Matrix:\n")Cell Means Design Matrix:print(X_cell) Duroc Landrace Yorkshire
[1,] 0 0 1
[2,] 0 0 1
[3,] 0 1 0
[4,] 0 1 0
[5,] 1 0 0
[6,] 1 0 0# Response vector
y <- pig_data$litter_size
# Normal equations: X'X b = X'y
XtX <- t(X_cell) %*% X_cell
Xty <- t(X_cell) %*% y
cat("\nX'X (3×3):\n")
X'X (3×3):print(XtX) Duroc Landrace Yorkshire
Duroc 2 0 0
Landrace 0 2 0
Yorkshire 0 0 2cat("\nX'y (3×1):\n")
X'y (3×1):print(Xty) [,1]
Duroc 19
Landrace 21
Yorkshire 23# Solve for estimates
b_cell <- solve(XtX) %*% Xty
cat("\nParameter estimates:\n")
Parameter estimates:print(b_cell) [,1]
Duroc 9.5
Landrace 10.5
Yorkshire 11.5# These are just the group means!
cat("\nVerify these equal group means:\n")
Verify these equal group means:print(tapply(y, pig_data$breed, mean)) Duroc Landrace Yorkshire
9.5 10.5 11.5 Interpretation:
In the cell means model with balanced data, the estimates are simply the group means. The math gives us:
\[ \hat{\mu}_i = \frac{\sum_{j=1}^{n_i} y_{ij}}{n_i} = \bar{y}_{i\cdot} \]
This is intuitive: our best estimate of a group’s mean is the average of observations in that group!
Model: \(y_{ij} = \mu + \alpha_i + e_{ij}\) with \(\alpha_{\text{Duroc}} = 0\) (Duroc as reference)
# Use R's default reference cell coding
X_effects <- model.matrix(~ breed, data = pig_data)
cat("Effects Model Design Matrix:\n")Effects Model Design Matrix:print(X_effects) (Intercept) breedLandrace breedYorkshire
1 1 0 1
2 1 0 1
3 1 1 0
4 1 1 0
5 1 0 0
6 1 0 0
attr(,"assign")
[1] 0 1 1
attr(,"contrasts")
attr(,"contrasts")$breed
[1] "contr.treatment"# Note: Duroc is reference (alphabetically first), so no column for it
# breedLandrace and breedYorkshire are deviations from Duroc
# Solve using lm() for comparison
fit_effects <- lm(litter_size ~ breed, data = pig_data)
cat("\nParameter estimates (reference cell coding):\n")
Parameter estimates (reference cell coding):print(coef(fit_effects)) (Intercept) breedLandrace breedYorkshire
9.5 1.0 2.0 # Manual calculation
XtX_eff <- t(X_effects) %*% X_effects
Xty_eff <- t(X_effects) %*% y
b_effects <- solve(XtX_eff) %*% Xty_eff
cat("\nManual calculation:\n")
Manual calculation:print(b_effects) [,1]
(Intercept) 9.5
breedLandrace 1.0
breedYorkshire 2.0Interpretation (with Duroc as reference):
# Cell means model: estimates are group means
cat("Cell Means Model Estimates:\n")Cell Means Model Estimates:cat("Duroc:", b_cell[1], "\n")Duroc: 9.5 cat("Landrace:", b_cell[2], "\n")Landrace: 10.5 cat("Yorkshire:", b_cell[3], "\n\n")Yorkshire: 11.5 # Effects model: reconstruct group means
cat("Effects Model - Reconstructed Group Means:\n")Effects Model - Reconstructed Group Means:cat("Duroc (reference):", b_effects[1], "\n")Duroc (reference): 9.5 cat("Landrace:", b_effects[1] + b_effects[2], "\n")Landrace: 10.5 cat("Yorkshire:", b_effects[1] + b_effects[3], "\n\n")Yorkshire: 11.5 # Both give identical fitted values
fitted_cell <- X_cell %*% b_cell
fitted_effects <- X_effects %*% b_effects
cat("Fitted values match:", all.equal(fitted_cell[,1], fitted_effects[,1]), "\n")Fitted values match: names for current but not for target Different coding schemes (cell means vs. effects model) give different parameter estimates, but:
The choice of coding affects parameter interpretation, not model fit.
Now let’s analyze a larger dataset: body weight (kg) for 20 broiler chickens at 42 days of age, by sex.
# Load data
broiler_data <- read.csv("data/broiler_bodyweight_sex.csv")
cat("Data structure:\n")Data structure:str(broiler_data)'data.frame': 20 obs. of 3 variables:
$ bird_id : int 1 2 3 4 5 6 7 8 9 10 ...
$ sex : chr "Male" "Male" "Male" "Male" ...
$ body_weight_kg: num 2.85 2.92 2.78 2.88 2.95 2.82 2.9 2.87 2.93 2.8 ...cat("\nFirst few rows:\n")
First few rows:head(broiler_data)| bird_id | sex | body_weight_kg |
|---|---|---|
| 1 | Male | 2.85 |
| 2 | Male | 2.92 |
| 3 | Male | 2.78 |
| 4 | Male | 2.88 |
| 5 | Male | 2.95 |
| 6 | Male | 2.82 |
# Summary statistics by sex
cat("\nSummary statistics by sex:\n")
Summary statistics by sex:by(broiler_data$body_weight_kg, broiler_data$sex, function(x) {
cat(sprintf(" n = %d, mean = %.3f, sd = %.3f, min = %.3f, max = %.3f\n",
length(x), mean(x), sd(x), min(x), max(x)))
}) n = 10, mean = 2.496, sd = 0.032, min = 2.450, max = 2.550
n = 10, mean = 2.870, sd = 0.057, min = 2.780, max = 2.950broiler_data$sex: Female
NULL
------------------------------------------------------------
broiler_data$sex: Male
NULL# Box plot
boxplot(body_weight_kg ~ sex, data = broiler_data,
main = "Broiler Body Weight by Sex",
xlab = "Sex",
ylab = "Body Weight (kg)",
col = c("lightblue", "lightpink"))
# Add means
means <- tapply(broiler_data$body_weight_kg, broiler_data$sex, mean)
points(1:2, means, pch = 19, col = "red", cex = 1.5)
legend("topright", legend = "Group mean", pch = 19, col = "red")
Observations:
# Cell means model: estimate mean for each sex
X_cell_broiler <- model.matrix(~ sex - 1, data = broiler_data)
y_broiler <- broiler_data$body_weight_kg
# Solve normal equations
XtX_broiler <- t(X_cell_broiler) %*% X_cell_broiler
Xty_broiler <- t(X_cell_broiler) %*% y_broiler
cat("X'X:\n")X'X:print(XtX_broiler) sexFemale sexMale
sexFemale 10 0
sexMale 0 10cat("\nX'y:\n")
X'y:print(Xty_broiler) [,1]
sexFemale 24.96
sexMale 28.70b_cell_broiler <- solve(XtX_broiler) %*% Xty_broiler
cat("\nEstimated means (kg):\n")
Estimated means (kg):print(b_cell_broiler) [,1]
sexFemale 2.496
sexMale 2.870# Compute residuals and fit statistics
fitted_broiler <- X_cell_broiler %*% b_cell_broiler
residuals_broiler <- y_broiler - fitted_broiler
SSE_broiler <- sum(residuals_broiler^2)
SST_broiler <- sum((y_broiler - mean(y_broiler))^2)
R2_broiler <- 1 - SSE_broiler / SST_broiler
cat("\nModel fit:\n")
Model fit:cat("SSE =", round(SSE_broiler, 4), "\n")SSE = 0.0386 cat("SST =", round(SST_broiler, 4), "\n")SST = 0.738 cat("R² =", round(R2_broiler, 4), "\n")R² = 0.9476 # Fit using lm() with effects model
fit_broiler <- lm(body_weight_kg ~ sex, data = broiler_data)
cat("Effects model summary:\n")Effects model summary:summary(fit_broiler)
Call:
lm(formula = body_weight_kg ~ sex, data = broiler_data)
Residuals:
Min 1Q Median 3Q Max
-0.0900 -0.0285 0.0020 0.0310 0.0800
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.49600 0.01465 170.36 < 2e-16 ***
sexMale 0.37400 0.02072 18.05 5.62e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.04633 on 18 degrees of freedom
Multiple R-squared: 0.9476, Adjusted R-squared: 0.9447
F-statistic: 325.8 on 1 and 18 DF, p-value: 5.617e-13# Compare with manual calculation
X_eff_broiler <- model.matrix(fit_broiler)
b_eff_broiler <- solve(t(X_eff_broiler) %*% X_eff_broiler) %*%
t(X_eff_broiler) %*% y_broiler
cat("\nManual calculation of effects model:\n")
Manual calculation of effects model:print(b_eff_broiler) [,1]
(Intercept) 2.496
sexMale 0.374With Female as the reference category (alphabetically first):
The lm() output shows:
This provides strong evidence that male broilers weigh more than females at 42 days of age.
We’ll cover the details of hypothesis testing in Weeks 5-6.
# The estimates should match simple group means
group_means_broiler <- tapply(y_broiler, broiler_data$sex, mean)
cat("Group means from data:\n")Group means from data:print(group_means_broiler)Female Male
2.496 2.870 cat("\nDifference (Male - Female):\n")
Difference (Male - Female):print(group_means_broiler["Male"] - group_means_broiler["Female"]) Male
0.374 cat("\nThese match our effects model estimates:\n")
These match our effects model estimates:cat("Intercept (Female mean):", coef(fit_broiler)[1], "\n")Intercept (Female mean): 2.496 cat("sexMale (difference):", coef(fit_broiler)[2], "\n")sexMale (difference): 0.374 R provides several ways to construct design matrices. Understanding both manual and automatic approaches deepens your understanding.
For full control and learning, build \(\mathbf{X}\) manually:
# Example: 2 groups, 4 observations
group <- c("A", "A", "B", "B")
y_vals <- c(5, 7, 3, 4)
# Cell means model: create indicator columns
n <- length(y_vals)
X_manual <- matrix(0, nrow = n, ncol = 2)
X_manual[group == "A", 1] <- 1
X_manual[group == "B", 2] <- 1
colnames(X_manual) <- c("GroupA", "GroupB")
cat("Manually constructed design matrix:\n")Manually constructed design matrix:print(X_manual) GroupA GroupB
[1,] 1 0
[2,] 1 0
[3,] 0 1
[4,] 0 1model.matrix()The model.matrix() function automatically creates design matrices from formulas:
# Create a data frame
example_data <- data.frame(group = group, y = y_vals)
# Cell means model (no intercept)
X_cell_auto <- model.matrix(~ group - 1, data = example_data)
cat("Cell means (~ group - 1):\n")Cell means (~ group - 1):print(X_cell_auto) groupA groupB
1 1 0
2 1 0
3 0 1
4 0 1
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"# Effects model (with intercept)
X_effects_auto <- model.matrix(~ group, data = example_data)
cat("\nEffects model (~ group):\n")
Effects model (~ group):print(X_effects_auto) (Intercept) groupB
1 1 0
2 1 0
3 1 1
4 1 1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"~ group: Effects model with intercept (reference cell coding)~ group - 1: Cell means model (no intercept, one column per level)~ 0 + group: Equivalent to ~ group - 1~ x1 + x2: Multiple predictors (intercept + x1 + x2)~ group + x: ANCOVA (categorical + continuous)~ group * x: Includes main effects and interactionAlways verify your design matrix:
# Function to check design matrix properties
check_design <- function(X, name = "X") {
cat("\n=== Design Matrix Check:", name, "===\n")
cat("Dimensions:", nrow(X), "×", ncol(X), "\n")
cat("Rank:", qr(X)$rank, "\n")
cat("Full rank?", qr(X)$rank == ncol(X), "\n")
# Check for linear dependencies
if (qr(X)$rank < ncol(X)) {
cat("WARNING: Matrix is rank deficient!\n")
cat("Number of parameters:", ncol(X), "\n")
cat("Effective rank:", qr(X)$rank, "\n")
}
}
# Check our broiler design matrices
check_design(X_cell_broiler, "Cell Means (Broiler)")
=== Design Matrix Check: Cell Means (Broiler) ===
Dimensions: 20 × 2
Rank: 2
Full rank? TRUE check_design(model.matrix(fit_broiler), "Effects Model (Broiler)")
=== Design Matrix Check: Effects Model (Broiler) ===
Dimensions: 20 × 2
Rank: 2
Full rank? TRUE Design Matrix is the Bridge: \(\mathbf{X}\) connects raw data to the model \(\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{e}\)
Predictor Types:
Coding Schemes:
Gauss-Markov Assumptions:
Different Coding, Same Fit: Coding schemes change parameter interpretation but not fitted values, residuals, or \(R^2\)
Always Check: Verify dimensions, rank, and linear independence of \(\mathbf{X}\)
Next week (Week 4: Simple Linear Regression), we’ll:
The design matrix concepts from this week provide the foundation for all subsequent work with linear models.
| Function | Purpose | Example |
|---|---|---|
model.matrix() |
Create design matrix from formula | model.matrix(~ breed, data) |
cbind() |
Combine vectors/matrices by columns | cbind(1, x1, x2) |
qr()$rank |
Compute matrix rank | qr(X)$rank |
solve() |
Matrix inverse | solve(XtX) |
t() |
Matrix transpose | t(X) |
%*% |
Matrix multiplication | X %*% beta |
Previous: Week 2: Linear Algebra Essentials Next: Week 4: Simple Linear Regression