2  Week 2: Linear Algebra Essentials

NoteLearning Objectives

By the end of this week, you will be able to:

1. Define and compute the rank of a matrix
2. Understand linear independence and its implications for solving equations
3. Compute regular inverses and generalized inverses of matrices
4. Solve systems of linear equations using matrix methods
5. Recognize when a system has unique, infinite, or no solutions

2.1 Why Linear Algebra Matters for Linear Models

Last week, we worked with a simple model where \(\mathbf{X}'\mathbf{X}\) was just a scalar (the sample size \(n\)). Inverting was trivial: \((n)^{-1} = 1/n\).

But what happens when:

• We have multiple effects (breeds, herds, years)?
• Some combinations are missing (not all breeds on all farms)?
• We have more parameters than we can uniquely estimate?

In these cases, \(\mathbf{X}'\mathbf{X}\) becomes a larger matrix that may not be full rank - meaning it doesn’t have a regular inverse. Understanding matrix rank and generalized inverses is essential for:

• Building ANOVA models
• Handling unbalanced data
• Understanding estimability
• Working toward mixed models and BLUP

2.2 Vector Spaces and Linear Independence

2.2.1 Vectors as Points in Space

A vector \(\mathbf{v} = \begin{bmatrix} v_1 \\ v_2 \\ \vdots \\ v_n \end{bmatrix}\) can be thought of as a point in \(n\)-dimensional space.

For example, \(\mathbf{v} = \begin{bmatrix} 2 \\ 3 \end{bmatrix}\) is a point in 2-D space.

2.2.2 Linear Combinations

A linear combination of vectors \(\mathbf{v}_1, \mathbf{v}_2, \ldots, \mathbf{v}_k\) is:

\[ c_1\mathbf{v}_1 + c_2\mathbf{v}_2 + \cdots + c_k\mathbf{v}_k \]

where \(c_1, c_2, \ldots, c_k\) are scalars.

Example:

v1 <- c(1, 0)
v2 <- c(0, 1)

# Linear combination: 2*v1 + 3*v2
result <- 2*v1 + 3*v2
print("2*v1 + 3*v2 =")

[1] "2*v1 + 3*v2 ="

print(result)

[1] 2 3

2.2.3 Linear Independence

Vectors \(\mathbf{v}_1, \mathbf{v}_2, \ldots, \mathbf{v}_k\) are linearly independent if the only solution to:

\[ c_1\mathbf{v}_1 + c_2\mathbf{v}_2 + \cdots + c_k\mathbf{v}_k = \mathbf{0} \]

is \(c_1 = c_2 = \cdots = c_k = 0\).

In other words, no vector can be written as a linear combination of the others.

Example of linear independence:

# These vectors are linearly independent
v1 <- c(1, 0)
v2 <- c(0, 1)

# No combination (except 0*v1 + 0*v2) gives the zero vector
# They point in different directions

Example of linear dependence:

# These vectors are linearly dependent
v1 <- c(1, 2)
v2 <- c(2, 4)  # v2 = 2*v1

# v2 can be written as a multiple of v1
# Therefore: 2*v1 - 1*v2 = 0 (non-trivial combination)
print("2*v1 - v2 =")

[1] "2*v1 - v2 ="

print(2*v1 - v2)

[1] 0 0

NoteConnection to Linear Models

In a design matrix \(\mathbf{X}\):

• Each column represents an effect (e.g., breed, herd)
• Columns are linearly independent → all effects are estimable
• Columns are linearly dependent → we have rank deficiency

Rank deficiency means we can’t uniquely estimate all parameters!

2.3 Rank of a Matrix

2.3.1 Definition

The rank of a matrix \(\mathbf{A}\), denoted \(r(\mathbf{A})\), is the number of linearly independent rows (or columns).

Key properties:

• \(r(\mathbf{A}) \leq \min(m, n)\) for an \(m \times n\) matrix
• \(r(\mathbf{A}) = r(\mathbf{A}')\)
• \(r(\mathbf{AB}) \leq \min(r(\mathbf{A}), r(\mathbf{B}))\)
• \(r(\mathbf{A}'\mathbf{A}) = r(\mathbf{A})\)

2.3.2 Full Rank vs. Rank Deficient

A matrix \(\mathbf{A}\) (\(m \times n\)) is:

• Full rank if \(r(\mathbf{A}) = \min(m, n)\)
• Rank deficient if \(r(\mathbf{A}) < \min(m, n)\)

For a square matrix (\(n \times n\)):

• Full rank means \(r(\mathbf{A}) = n\)
• Rank deficient means \(r(\mathbf{A}) < n\)

2.3.3 Computing Rank in R

library(MASS)

# Full rank matrix
A_full <- matrix(c(1, 2, 3, 4), nrow = 2, ncol = 2)
print("Full rank matrix:")

[1] "Full rank matrix:"

print(A_full)

     [,1] [,2]
[1,]    1    3
[2,]    2    4

print(paste("Rank:", qr(A_full)$rank))

[1] "Rank: 2"

# Rank deficient matrix
A_deficient <- matrix(c(1, 2, 2, 4), nrow = 2, ncol = 2)
print("Rank deficient matrix:")

[1] "Rank deficient matrix:"

print(A_deficient)

     [,1] [,2]
[1,]    1    2
[2,]    2    4

print(paste("Rank:", qr(A_deficient)$rank))

[1] "Rank: 1"

# Note: second column is 2 times first column

2.3.4 Small Example: Design Matrix Rank

Consider a beef feedlot study with 4 steers in 2 pens:

# Data: 4 steers, 2 pens (2 steers each)
# Pen: 1, 1, 2, 2

# Design matrix (effects model: overall mean + pen effect)
X <- matrix(c(
  1, 1, 0,  # Steer 1: in pen 1
  1, 1, 0,  # Steer 2: in pen 1
  1, 0, 1,  # Steer 3: in pen 2
  1, 0, 1   # Steer 4: in pen 2
), nrow = 4, ncol = 3, byrow = TRUE)

colnames(X) <- c("mu", "pen1", "pen2")
print("Design matrix X:")

[1] "Design matrix X:"

print(X)

     mu pen1 pen2
[1,]  1    1    0
[2,]  1    1    0
[3,]  1    0    1
[4,]  1    0    1

print(paste("Number of columns (parameters):", ncol(X)))

[1] "Number of columns (parameters): 3"

print(paste("Rank of X:", qr(X)$rank))

[1] "Rank of X: 2"

# X'X
XtX <- t(X) %*% X
print("X'X:")

[1] "X'X:"

print(XtX)

     mu pen1 pen2
mu    4    2    2
pen1  2    2    0
pen2  2    0    2

print(paste("Rank of X'X:", qr(XtX)$rank))

[1] "Rank of X'X: 2"

WarningRank Deficiency in Design Matrices

Problem: We have 3 columns but rank = 2. Why?

The first column (overall mean) equals the sum of the second and third columns (pen effects):

\[ \text{Column 1} = \text{Column 2} + \text{Column 3} \]

This is called linear dependence or singularity.

Consequences: - \(\mathbf{X}'\mathbf{X}\) is not invertible (no regular inverse exists) - Individual parameters (\(\mu\), \(\alpha_1\), \(\alpha_2\)) cannot be uniquely estimated - We need constraints or generalized inverses to obtain solutions

Solution: Use generalized inverses (next section) or apply constraints to make the design matrix full rank.

2.4 Matrix Inverses

2.4.1 Regular Inverse

For a square matrix \(\mathbf{A}\) (\(n \times n\)), the inverse \(\mathbf{A}^{-1}\) satisfies:

\[ \mathbf{A}\mathbf{A}^{-1} = \mathbf{A}^{-1}\mathbf{A} = \mathbf{I} \]

Requirements for \(\mathbf{A}^{-1}\) to exist:

1. \(\mathbf{A}\) must be square (\(n \times n\))
2. \(\mathbf{A}\) must be full rank: \(r(\mathbf{A}) = n\)
3. Equivalently: \(\det(\mathbf{A}) \neq 0\) (determinant non-zero)

2.4.2 Computing Regular Inverse

# Full rank 2x2 matrix
A <- matrix(c(4, 3, 3, 2), nrow = 2, ncol = 2)
print("Matrix A:")

[1] "Matrix A:"

print(A)

     [,1] [,2]
[1,]    4    3
[2,]    3    2

# Compute inverse
A_inv <- solve(A)
print("A inverse:")

[1] "A inverse:"

print(A_inv)

     [,1] [,2]
[1,]   -2    3
[2,]    3   -4

# Verify: A * A^-1 = I
I_check <- A %*% A_inv
print("A * A^-1 (should be identity):")

[1] "A * A^-1 (should be identity):"

print(round(I_check, 10))

     [,1] [,2]
[1,]    1    0
[2,]    0    1

2.4.3 When Inverse Doesn’t Exist

# Rank deficient matrix
B <- matrix(c(1, 2, 2, 4), nrow = 2, ncol = 2)
print("Rank deficient matrix B:")

[1] "Rank deficient matrix B:"

print(B)

     [,1] [,2]
[1,]    1    2
[2,]    2    4

print(paste("Rank:", qr(B)$rank))

[1] "Rank: 1"

print(paste("Determinant:", det(B)))

[1] "Determinant: 0"

# Try to invert - this will fail
try({
  B_inv <- solve(B)
}, silent = FALSE)

Error in solve.default(B) : 
  Lapack routine dgesv: system is exactly singular: U[2,2] = 0

The error occurs because B is singular (not full rank). We need a generalized inverse instead!

2.5 Generalized Inverses

2.5.1 Definition

A generalized inverse (g-inverse) of matrix \(\mathbf{A}\) is any matrix \(\mathbf{A}^{-}\) that satisfies:

\[ \mathbf{A}\mathbf{A}^{-}\mathbf{A} = \mathbf{A} \tag{2.1}\]

Key facts:

• Generalized inverses exist for any matrix (even non-square, rank deficient)
• Generalized inverses are not unique (many possible g-inverses)
• If \(\mathbf{A}\) is full rank, then \(\mathbf{A}^{-} = \mathbf{A}^{-1}\) (unique)

WarningNOTATION EXTENSION

We now extend our notation from Week 1 to include generalized inverses:

• \(\mathbf{A}^{-}\) = generalized inverse (g-inverse) of matrix \(\mathbf{A}\)
  • Satisfies: \(\mathbf{A}\mathbf{A}^{-}\mathbf{A} = \mathbf{A}\) (Equation 12.5)
  • Not unique (multiple g-inverses exist)
  • Works for any matrix (even rank deficient)
• \(\mathbf{A}^{+}\) = Moore-Penrose inverse (special g-inverse)
  • Unique generalized inverse
  • Computed with ginv() in R

Key distinction from Week 1: - Regular inverse \(\mathbf{A}^{-1}\): requires square, full rank matrix - Generalized inverse \(\mathbf{A}^{-}\): works for any matrix, including rank deficient

This notation will be essential for handling overparameterized models and missing cells in ANOVA (Weeks 7-12).

2.5.2 Moore-Penrose Inverse

The Moore-Penrose inverse \(\mathbf{A}^{+}\) is a special g-inverse that satisfies four conditions:

1. \(\mathbf{A}\mathbf{A}^{+}\mathbf{A} = \mathbf{A}\)
2. \(\mathbf{A}^{+}\mathbf{A}\mathbf{A}^{+} = \mathbf{A}^{+}\)
3. \((\mathbf{A}\mathbf{A}^{+})' = \mathbf{A}\mathbf{A}^{+}\) (symmetric)
4. \((\mathbf{A}^{+}\mathbf{A})' = \mathbf{A}^{+}\mathbf{A}\) (symmetric)

The Moore-Penrose inverse is unique and computed in R using ginv() from the MASS package (Penrose 1955; Moore 1920).

2.5.3 Computing Generalized Inverse

library(MASS)

# Rank deficient matrix
B <- matrix(c(1, 2, 2, 4), nrow = 2, ncol = 2)
print("Rank deficient matrix B:")

[1] "Rank deficient matrix B:"

print(B)

     [,1] [,2]
[1,]    1    2
[2,]    2    4

print(paste("Rank:", qr(B)$rank))

[1] "Rank: 1"

# Compute Moore-Penrose generalized inverse
B_ginv <- ginv(B)
print("Generalized inverse B^-:")

[1] "Generalized inverse B^-:"

print(B_ginv)

     [,1] [,2]
[1,] 0.04 0.08
[2,] 0.08 0.16

# Verify the defining property: B * B^- * B = B
check <- B %*% B_ginv %*% B
print("B * B^- * B (should equal B):")

[1] "B * B^- * B (should equal B):"

print(round(check, 10))

     [,1] [,2]
[1,]    1    2
[2,]    2    4

# Note: B * B^- ≠ I (identity)
BB_ginv <- B %*% B_ginv
print("B * B^- (NOT identity):")

[1] "B * B^- (NOT identity):"

print(round(BB_ginv, 10))

     [,1] [,2]
[1,]  0.2  0.4
[2,]  0.4  0.8

ImportantKey Difference

Regular inverse: \(\mathbf{A}\mathbf{A}^{-1} = \mathbf{I}\)

Generalized inverse: \(\mathbf{A}\mathbf{A}^{-}\mathbf{A} = \mathbf{A}\) (but \(\mathbf{A}\mathbf{A}^{-} \neq \mathbf{I}\) in general)

2.5.4 Application to Normal Equations

For rank deficient \(\mathbf{X}'\mathbf{X}\):

\[ \mathbf{X}'\mathbf{X}\mathbf{b} = \mathbf{X}'\mathbf{y} \tag{2.2}\]

One solution is:

\[ \mathbf{b} = (\mathbf{X}'\mathbf{X})^{-}\mathbf{X}'\mathbf{y} \tag{2.3}\]

Important: This gives a solution, but not necessarily a unique solution. Different g-inverses give different \(\mathbf{b}\) vectors.

However, certain functions of \(\mathbf{b}\) (called estimable functions) are unique regardless of which g-inverse is used. We’ll explore this in Week 8.

2.6 Solving Systems of Linear Equations

2.6.1 General Form

A system of linear equations can be written as:

\[ \mathbf{A}\mathbf{x} = \mathbf{b} \tag{2.4}\]

Where:

• \(\mathbf{A}\) is an \(m \times n\) coefficient matrix
• \(\mathbf{x}\) is an \(n \times 1\) vector of unknowns
• \(\mathbf{b}\) is an \(m \times 1\) vector of constants

2.6.2 Three Cases

Case 1: Unique Solution

• \(\mathbf{A}\) is square and full rank
• Solution: \(\mathbf{x} = \mathbf{A}^{-1}\mathbf{b}\)

A <- matrix(c(2, 1, 1, 3), nrow = 2, ncol = 2)
b <- c(8, 7)

print("System: Ax = b")

[1] "System: Ax = b"

print("A:")

[1] "A:"

print(A)

     [,1] [,2]
[1,]    2    1
[2,]    1    3

print("b:")

[1] "b:"

print(b)

[1] 8 7

# Unique solution
x <- solve(A) %*% b
print("Solution x:")

[1] "Solution x:"

print(x)

     [,1]
[1,]  3.4
[2,]  1.2

# Verify
print("Ax (should equal b):")

[1] "Ax (should equal b):"

print(A %*% x)

     [,1]
[1,]    8
[2,]    7

Case 2: Infinite Solutions

• More unknowns than equations, or
• Rows/columns linearly dependent

# System: x1 + 2*x2 = 5
#         2*x1 + 4*x2 = 10  (second equation is 2× first)

A <- matrix(c(1, 2, 2, 4), nrow = 2, ncol = 2, byrow = TRUE)
b <- c(5, 10)

print("System (second equation redundant):")

[1] "System (second equation redundant):"

print(A)

     [,1] [,2]
[1,]    1    2
[2,]    2    4

print(paste("Rank:", qr(A)$rank))

[1] "Rank: 1"

# One solution using generalized inverse
A_ginv <- ginv(A)
x_particular <- A_ginv %*% b
print("One particular solution:")

[1] "One particular solution:"

print(x_particular)

     [,1]
[1,]    1
[2,]    2

# But infinitely many solutions exist!
# Any x = [5 - 2*t, t] for any t works

Case 3: No Solution (Inconsistent)

• Equations are contradictory

# System: x1 + 2*x2 = 5
#         2*x1 + 4*x2 = 12  (contradicts first equation!)

A <- matrix(c(1, 2, 2, 4), nrow = 2, ncol = 2, byrow = TRUE)
b_inconsistent <- c(5, 12)

print("Inconsistent system:")

[1] "Inconsistent system:"

print(A)

     [,1] [,2]
[1,]    1    2
[2,]    2    4

print("b:")

[1] "b:"

print(b_inconsistent)

[1]  5 12

# ginv() gives "least squares" solution (minimizes ||Ax - b||)
x_ls <- ginv(A) %*% b_inconsistent
print("Least squares solution (doesn't satisfy exactly):")

[1] "Least squares solution (doesn't satisfy exactly):"

print(x_ls)

     [,1]
[1,] 1.16
[2,] 2.32

print("Ax (does NOT equal b):")

[1] "Ax (does NOT equal b):"

print(A %*% x_ls)

     [,1]
[1,]  5.8
[2,] 11.6

2.6.3 Connection to Animal Breeding

In genetic evaluation:

• Case 1 (unique solution): Balanced experiments with all factor combinations
• Case 2 (infinite solutions): Missing cells, confounded effects → need constraints
• Case 3 (no solution): Usually doesn’t occur with real data (every observation provides information)

We almost always encounter Case 2 in animal breeding, which is why understanding generalized inverses and estimability is crucial.

2.7 Small Numerical Example: Pig Litter Size

Let’s work through a complete example with rank deficiency.

2.7.1 Problem Setup

Three breeds of pigs with unequal replication:

• Yorkshire: 3 litters
• Landrace: 2 litters
• Duroc: 1 litter

Litter sizes: Yorkshire (11, 12, 10), Landrace (10, 11), Duroc (9)

2.7.2 Data

# Litter size data
litter_size <- c(11, 12, 10, 10, 11, 9)
breed <- factor(c("Yorkshire", "Yorkshire", "Yorkshire",
                  "Landrace", "Landrace", "Duroc"))

data_pig <- data.frame(breed = breed, litter_size = litter_size)
print(data_pig)

      breed litter_size
1 Yorkshire          11
2 Yorkshire          12
3 Yorkshire          10
4  Landrace          10
5  Landrace          11
6     Duroc           9

2.7.3 Cell Means Model (Full Rank)

Model: \(y_{ij} = \mu_i + e_{ij}\) where \(i\) = breed, \(j\) = observation within breed

# Design matrix: indicator for each breed
X_cell <- model.matrix(~ breed - 1, data = data_pig)
print("Cell means design matrix:")

[1] "Cell means design matrix:"

print(X_cell)

  breedDuroc breedLandrace breedYorkshire
1          0             0              1
2          0             0              1
3          0             0              1
4          0             1              0
5          0             1              0
6          1             0              0
attr(,"assign")
[1] 1 1 1
attr(,"contrasts")
attr(,"contrasts")$breed
[1] "contr.treatment"

print(paste("Rank of X:", qr(X_cell)$rank))

[1] "Rank of X: 3"

# This is full rank! (3 columns, rank = 3)

# Solve using normal equations
y <- litter_size
XtX <- t(X_cell) %*% X_cell
Xty <- t(X_cell) %*% y

print("X'X:")

[1] "X'X:"

print(XtX)

               breedDuroc breedLandrace breedYorkshire
breedDuroc              1             0              0
breedLandrace           0             2              0
breedYorkshire          0             0              3

print("X'y:")

[1] "X'y:"

print(Xty)

               [,1]
breedDuroc        9
breedLandrace    21
breedYorkshire   33

# Invert (regular inverse exists)
b_cell <- solve(XtX) %*% Xty
print("Breed means (cell means model):")

[1] "Breed means (cell means model):"

print(b_cell)

               [,1]
breedDuroc      9.0
breedLandrace  10.5
breedYorkshire 11.0

# These are simply the breed averages!
print("Verify - breed averages:")

[1] "Verify - breed averages:"

print(tapply(litter_size, breed, mean))

    Duroc  Landrace Yorkshire 
      9.0      10.5      11.0 

2.7.4 Effects Model (Rank Deficient)

Model: \(y_{ij} = \mu + \alpha_i + e_{ij}\) where \(\mu\) = overall mean, \(\alpha_i\) = breed effect

# Design matrix: intercept + breed effects
X_effects <- model.matrix(~ breed, data = data_pig)
print("Effects model design matrix:")

[1] "Effects model design matrix:"

print(X_effects)

  (Intercept) breedLandrace breedYorkshire
1           1             0              1
2           1             0              1
3           1             0              1
4           1             1              0
5           1             1              0
6           1             0              0
attr(,"assign")
[1] 0 1 1
attr(,"contrasts")
attr(,"contrasts")$breed
[1] "contr.treatment"

print(paste("Number of columns:", ncol(X_effects)))

[1] "Number of columns: 3"

print(paste("Rank of X:", qr(X_effects)$rank))

[1] "Rank of X: 3"

# Rank deficient! (4 columns but rank = 3)

# X'X is singular
XtX_effects <- t(X_effects) %*% X_effects
print("X'X (effects model):")

[1] "X'X (effects model):"

print(XtX_effects)

               (Intercept) breedLandrace breedYorkshire
(Intercept)              6             2              3
breedLandrace            2             2              0
breedYorkshire           3             0              3

print(paste("Rank of X'X:", qr(XtX_effects)$rank))

[1] "Rank of X'X: 3"

print(paste("Determinant of X'X:", det(XtX_effects)))

[1] "Determinant of X'X: 6"

# Cannot use regular inverse - use generalized inverse
XtX_ginv <- ginv(XtX_effects)
Xty_effects <- t(X_effects) %*% y

b_effects <- XtX_ginv %*% Xty_effects
print("Solution using generalized inverse:")

[1] "Solution using generalized inverse:"

print(b_effects)

     [,1]
[1,]  9.0
[2,]  1.5
[3,]  2.0

Note: The solution depends on which generalized inverse we use. R’s ginv() sets one parameter to zero by default.

Let’s verify estimable functions are consistent:

# Contrast: Yorkshire - Duroc
# In cell means model: mu_Yorkshire - mu_Duroc
contrast_cell <- b_cell[1] - b_cell[3]

# In effects model: alpha_Yorkshire - alpha_Duroc
# (mu cancels out)
contrast_effects <- b_effects[2] - b_effects[3]

print("Yorkshire - Duroc difference:")

[1] "Yorkshire - Duroc difference:"

print(paste("Cell means model:", round(contrast_cell, 4)))

[1] "Cell means model: -2"

print(paste("Effects model:", round(contrast_effects, 4)))

[1] "Effects model: -0.5"

print(paste("Match:", all.equal(contrast_cell, contrast_effects)))

[1] "Match: Mean relative difference: 0.75"

The contrast (difference between breeds) is the same regardless of parameterization! This is an estimable function.

2.8 Realistic Livestock Application

2.8.1 Scenario

A sheep researcher measures fleece weight (kg) for 3 breeds across 2 farms, but not all breed × farm combinations are present (unbalanced design).

2.8.2 Data Structure

set.seed(456)

# Farm-Breed combinations (missing cells)
farm <- factor(c(rep("A", 5), rep("B", 4)))
breed <- factor(c(rep("Merino", 3), rep("Suffolk", 2),  # Farm A
                  rep("Suffolk", 2), rep("Romney", 2)))  # Farm B

# Note: No Merino on Farm B, No Romney on Farm A

fleece_weight <- c(
  # Farm A: Merino
  4.8, 5.1, 4.9,
  # Farm A: Suffolk
  5.5, 5.8,
  # Farm B: Suffolk
  5.2, 5.4,
  # Farm B: Romney
  5.9, 6.1
)

sheep_data <- data.frame(
  farm = farm,
  breed = breed,
  fleece_weight = fleece_weight
)

print("Sheep fleece weight data:")

[1] "Sheep fleece weight data:"

print(sheep_data)

  farm   breed fleece_weight
1    A  Merino           4.8
2    A  Merino           5.1
3    A  Merino           4.9
4    A Suffolk           5.5
5    A Suffolk           5.8
6    B Suffolk           5.2
7    B Suffolk           5.4
8    B  Romney           5.9
9    B  Romney           6.1

# Check structure
table(sheep_data$farm, sheep_data$breed)

   
    Merino Romney Suffolk
  A      3      0       2
  B      0      2       2

2.8.3 Analysis

# Full model: farm + breed + farm×breed
# But missing cells cause rank deficiency

X_full <- model.matrix(~ farm * breed, data = sheep_data)
print("Design matrix (first 5 rows):")

[1] "Design matrix (first 5 rows):"

print(head(X_full, 5))

  (Intercept) farmB breedRomney breedSuffolk farmB:breedRomney
1           1     0           0            0                 0
2           1     0           0            0                 0
3           1     0           0            0                 0
4           1     0           0            1                 0
5           1     0           0            1                 0
  farmB:breedSuffolk
1                  0
2                  0
3                  0
4                  0
5                  0

print(paste("Number of columns:", ncol(X_full)))

[1] "Number of columns: 6"

print(paste("Rank:", qr(X_full)$rank))

[1] "Rank: 4"

# Rank deficient due to missing cells

# Solve using generalized inverse
y_sheep <- fleece_weight
XtX <- t(X_full) %*% X_full
Xty <- t(X_full) %*% y_sheep

print("X'X:")

[1] "X'X:"

print(XtX)

                   (Intercept) farmB breedRomney breedSuffolk farmB:breedRomney
(Intercept)                  9     4           2            4                 2
farmB                        4     4           2            2                 2
breedRomney                  2     2           2            0                 2
breedSuffolk                 4     2           0            4                 0
farmB:breedRomney            2     2           2            0                 2
farmB:breedSuffolk           2     2           0            2                 0
                   farmB:breedSuffolk
(Intercept)                         2
farmB                               2
breedRomney                         0
breedSuffolk                        2
farmB:breedRomney                   0
farmB:breedSuffolk                  2

print(paste("Rank of X'X:", qr(XtX)$rank))

[1] "Rank of X'X: 4"

# Generalized inverse solution
XtX_ginv <- ginv(XtX)
b_sheep <- XtX_ginv %*% Xty

print("Parameter estimates:")

[1] "Parameter estimates:"

print(b_sheep)

            [,1]
[1,]  4.93333333
[2,]  0.07333333
[3,]  0.49666667
[4,]  0.71666667
[5,]  0.49666667
[6,] -0.42333333

# Estimable contrasts within farm
# Suffolk at Farm A vs Suffolk at Farm B
# (these cells both exist, so difference is estimable)

# Verify with lm()
model_sheep <- lm(fleece_weight ~ farm * breed, data = sheep_data)
print("lm() coefficients:")

[1] "lm() coefficients:"

print(coef(model_sheep))

       (Intercept)              farmB        breedRomney       breedSuffolk 
         4.9333333         -0.3500000          1.4166667          0.7166667 
 farmB:breedRomney farmB:breedSuffolk 
                NA                 NA 

TipPractical Advice

When dealing with missing cells or rank deficiency:

1. Use cell means models when possible (always full rank)
2. For effects models, use generalized inverse (ginv())
3. Only interpret estimable functions (contrasts, differences)
4. Week 8 will formalize estimability criteria

2.9 Summary

This week covered essential linear algebra for linear models:

2.9.1 Key Concepts

NoteWhat We Learned

1. Linear independence: Vectors that can’t be written as combinations of each other
2. Rank: Number of linearly independent rows/columns
   • Full rank: \(r(\mathbf{A}) = \min(m, n)\)
   • Rank deficient: \(r(\mathbf{A}) < \min(m, n)\)
3. Regular inverse \(\mathbf{A}^{-1}\):
   • Requires square, full rank matrix
   • \(\mathbf{A}\mathbf{A}^{-1} = \mathbf{I}\)
4. Generalized inverse \(\mathbf{A}^{-}\):
   • Works for any matrix
   • Not unique
   • \(\mathbf{A}\mathbf{A}^{-}\mathbf{A} = \mathbf{A}\) (Equation 12.5)
5. Moore-Penrose inverse \(\mathbf{A}^{+}\):
   • Unique generalized inverse
   • Computed with ginv() in R
6. Solving \(\mathbf{A}\mathbf{x} = \mathbf{b}\) (Equation 2.4):
   • Full rank → unique solution
   • Rank deficient → infinite solutions (use g-inverse)
   • Inconsistent → no exact solution

2.9.2 Implications for Linear Models

• Cell means models are always full rank
• Effects models are often rank deficient (overparameterized)
• Generalized inverses let us solve normal equations
• Only estimable functions have unique values
• Understanding rank is key to understanding what we can estimate

2.9.3 Looking Ahead

Week 3: Build design matrices for different model types

Week 4-5: Regression models (usually full rank)

Week 7-8: ANOVA models (often rank deficient) and estimable functions

Week 12: Deep dive into non-full rank models and constraints

2.10 R Functions Reference

• qr(A)$rank - Compute rank of matrix A
• solve(A) - Regular inverse of A
• ginv(A) - Moore-Penrose generalized inverse (MASS package)
• det(A) - Determinant of A
• model.matrix() - Create design matrix from formula

2.11 Additional Resources

2.11.1 Key References

• Penrose (1955) - Original paper on generalized inverses
• Searle (1971) - Linear models textbook with extensive matrix algebra
• Searle and Gruber (2006) - Updated edition

2.11.2 Practice

Work through Exercise Set 2 to solidify these concepts!

---

Previous: Week 1: Course Overview

Next: Week 3: Building the Design Matrix Framework

Moore, Eliakim H. 1920. “On the Reciprocal of the General Algebraic Matrix.” Bulletin of the American Mathematical Society 26 (9): 394–95.

Penrose, Roger. 1955. “A Generalized Inverse for Matrices.” Mathematical Proceedings of the Cambridge Philosophical Society 51 (3): 406–13. https://doi.org/10.1017/S0305004100030401.

Searle, Shayle R. 1971. Linear Models. New York: John Wiley & Sons.

Searle, Shayle R, and Marvin HJ Gruber. 2006. Linear Models. 2nd ed. John Wiley & Sons. https://doi.org/10.1002/9780470316962.