Practical — Linear Algebra & Linear Models (STS-106)

Topics Covered

Objectives

1. Knowing the computational procedures and also their implementation using R.
2. Finding the inverse of a matrix using various methods.
3. Applying transformations on matrices.
4. Applying matrix operations to given data sets — determinant, eigenvalues, eigenvectors, transformations.
5. Summarizing the properties of data sets based on matrix operations.

The Fourteen Practicals

All fourteen are worked below, each set out as 1. Problem, 2. Aim, 3. Formula, 4. Calculation, 5. Result, with the R program (written without packages) and its real output inside the Calculation.

Practical 1 — Inverse by the Partition Method

1. Problem

Find the inverse of \(A = \begin{pmatrix} 2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{pmatrix}\) by the partition method.

2. Aim

To invert a matrix by partitioning it into blocks and inverting the smaller pieces.

3. Formula

\(S\) is the Schur complement of \(A_{11}\).

Applying it:

1. Partition, taking \(A_{11}\) as the leading \(2 \times 2\).
2. Invert \(A_{11}\).
3. Form the Schur complement and its inverse.
4. The off-diagonal blocks, then the leading block.
5. Assemble and check \(AA^{-1} = I\).

4. Calculation

since \(\tfrac23 + \tfrac1{12} = \tfrac34\) and \(-\tfrac13 + \tfrac1{12} = -\tfrac14\).

A <- matrix(c(2,1,1, 1,2,1, 1,1,2), 3, 3, byrow = TRUE)
A11 <- A[1:2, 1:2]; A12 <- A[1:2, 3, drop = FALSE]; A21 <- A[3, 1:2, drop = FALSE]; A22 <- A[3, 3, drop = FALSE]
A11i <- solve(A11); S <- A22 - A21 %*% A11i %*% A12; Si <- solve(S)
B12 <- -A11i %*% A12 %*% Si; B21 <- -Si %*% A21 %*% A11i
B11 <- A11i + A11i %*% A12 %*% Si %*% A21 %*% A11i
Ainv <- rbind(cbind(B11, B12), cbind(B21, Si)); Ainv * 4
round(max(abs(A %*% Ainv - diag(3))), 10)   # 0 -- A times its inverse is I

     [,1] [,2] [,3]
[1,]    3   -1   -1
[2,]   -1    3   -1
[3,]   -1   -1    3
[1] 0

5. Result

the matrix obtained in Unit 2 by Cayley–Hamilton and again by spectral decomposition (Practical 7): three independent routes, one answer. \(\checkmark\)

Practical 2 — Solution of Linear Equations by the Sweep-Out Method

1. Problem

Solve \(A\mathbf{x} = \mathbf{b}\) by the sweep-out (Gauss–Jordan) method, where \(A = \begin{pmatrix} 2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{pmatrix}\) and \(\mathbf{b} = (7, 8, 9)'\).

2. Aim

To solve a system of linear equations by sweeping the augmented matrix to the identity.

3. Formula

At each step, divide the pivot row by the pivot, then subtract multiples of it from every other row so that its column becomes a unit vector. Augmenting with \(I\) instead of \(\mathbf{b}\) produces \(A^{-1}\) with identical work.

Applying it:

1. First pivot: \(R_1 \div 2\); \(R_2 - R_1\), \(R_3 - R_1\).
2. Second pivot: \(R_2 \div \tfrac32\); \(R_1 - \tfrac12R_2\), \(R_3 - \tfrac12R_2\).
3. Third pivot: \(R_3 \div \tfrac43\); \(R_1 - \tfrac13R_3\), \(R_2 - \tfrac13R_3\).
4. Read \(\mathbf{x}\) and substitute back.

4. Calculation

Substituting back: \(2(1)+2+3 = 7\), \(1+2(2)+3 = 8\), \(1+2+2(3) = 9\). \(\checkmark\)

A <- matrix(c(2,1,1, 1,2,1, 1,1,2), 3, 3, byrow = TRUE)
sweep <- function(M) {                 # Gauss-Jordan on an augmented matrix
  for (k in seq_len(nrow(M))) {
    M[k, ] <- M[k, ] / M[k, k]
    for (i in seq_len(nrow(M))[-k]) M[i, ] <- M[i, ] - M[i, k] * M[k, ]
  }
  M
}
sweep(cbind(A, c(7, 8, 9)))[, 4]       # the solution
sweep(cbind(A, diag(3)))[, 4:6] * 4    # and, augmenting with I, 4 times the inverse

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

5. Result

\(\mathbf{x} = (1, 2, 3)'\). The same sweep, applied to \((A \mid I)\), gives \(A^{-1}\) of Practical 1.

Practical 3 — Solution of Linear Equations by the Doolittle Method

1. Problem

Solve \(A\mathbf{x} = \mathbf{b}\) with \(A = \begin{pmatrix} 2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{pmatrix}\) and \(\mathbf{b} = (7, 8, 9)'\) by the Doolittle (LU) method, and find \(\det A\).

2. Aim

To factorise a matrix into unit lower and upper triangular factors once, and solve by forward and back substitution.

3. Formula

Applying it:

1. The first row of \(U\) and first column of \(L\).
2. The second row and column; then \(u_{33}\).
3. Check \(LU = A\) on one entry.
4. Forward substitution, then back substitution.

4. Calculation

\(u_{1j} = a_{1j}\): first row of \(U\) is \((2, 1, 1)\); \(\ell_{21} = \ell_{31} = \tfrac12\).

Check: row 3 of \(L\) against column 3 of \(U\): \(\tfrac12 + \tfrac16 + \tfrac43 = 2 = a_{33}\). \(\checkmark\)

A <- matrix(c(2,1,1, 1,2,1, 1,1,2), 3, 3, byrow = TRUE)
doolittle <- function(A) {             # Doolittle LU, written out rather than calling a library
  n <- nrow(A); L <- diag(n); U <- matrix(0, n, n)
  for (i in 1:n) {
    for (j in i:n) U[i, j] <- A[i, j] - sum(L[i, ] * U[, j])
    if (i < n) for (j in (i+1):n)
      L[j, i] <- (A[j, i] - sum(L[j, ] * U[, i])) / U[i, i]
  }
  list(L = L, U = U)
}
lu <- doolittle(A); lu$L; lu$U
max(abs(lu$L %*% lu$U - A))          # 0 -- the factorisation is exact
z <- forwardsolve(lu$L, c(7, 8, 9)); backsolve(lu$U, z)
prod(diag(lu$U))                     # det(A)

     [,1]      [,2] [,3]
[1,]  1.0 0.0000000    0
[2,]  0.5 1.0000000    0
[3,]  0.5 0.3333333    1
     [,1] [,2]     [,3]
[1,]    2  1.0 1.000000
[2,]    0  1.5 0.500000
[3,]    0  0.0 1.333333
[1] 0
[1] 1 2 3
[1] 4

5. Result

\(\mathbf{x} = (1, 2, 3)'\) and \(\det A = 1 \times 2 \times \tfrac32 \times \tfrac43 = 4\). Factorise once and any number of right-hand sides cost only the two substitutions — which is why this, and not Gauss–Jordan, is what numerical libraries do.

Practical 4 — Moore–Penrose Inverse by the Penrose Method

1. Problem

Find the Moore–Penrose inverse of \(A = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ 1 & 1 \end{pmatrix}\) and verify the four Penrose conditions.

2. Aim

To find the Moore–Penrose inverse of a matrix of full column rank and verify that it satisfies all four Penrose conditions.

3. Formula

Applying it:

1. Check the rank.
2. Form and invert \(A'A\); multiply by \(A'\).
3. Form \(A^{+}A\) and \(AA^{+}\); check the four conditions.

4. Calculation

The two columns are not multiples of one another, so \(\operatorname{rank}(A) = 2\), full column rank.

A <- matrix(c(1,0, 0,1, 1,1), 3, 2, byrow = TRUE)
Ap <- solve(t(A) %*% A) %*% t(A); Ap * 3
c(i = max(abs(A %*% Ap %*% A - A)), ii = max(abs(Ap %*% A %*% Ap - Ap)),
  iii = max(abs(t(A %*% Ap) - A %*% Ap)), iv = max(abs(t(Ap %*% A) - Ap %*% A)))   # all 0
sum(diag(A %*% Ap))                  # the trace of the projection = rank

     [,1] [,2] [,3]
[1,]    2   -1    1
[2,]   -1    2    1
  i  ii iii  iv 
  0   0   0   0 
[1] 2

5. Result

\(A^{+} = \frac13\begin{pmatrix} 2 & -1 & 1 \\ -1 & 2 & 1 \end{pmatrix}\), satisfying all four conditions. \(AA^{+}\) is not \(I_3\) and could not be (its rank is 2): it is the projection onto the column space of \(A\), with trace 2 = rank. Applied to \(\mathbf{y}\), \(A^{+}\) returns the least-squares coefficients (worked also in Unit 1, Example 1.3).

Practical 5 — Generalized Inverse of a Matrix

1. Problem

Find a generalized inverse of the singular matrix \(S = \begin{pmatrix} 1 & 1 & 0 \\ 1 & 1 & 0 \\ 0 & 0 & 2 \end{pmatrix}\), and check whether it is the Moore–Penrose inverse.

2. Aim

To construct a generalized inverse of a singular matrix from a non-singular submatrix, and to see that it need not be the Moore–Penrose inverse.

3. Formula

Here \(\det S = 0\) and \(\operatorname{rank}(S) = 2\).

Applying it:

1. Find a non-singular \(2 \times 2\) submatrix and invert it.
2. Place its entries back in the transposed positions, zeros elsewhere.
3. Verify \(SGS = S\); then test the other Penrose conditions.

4. Calculation

Rows 1 and 3, columns 1 and 3: \(\begin{pmatrix} 1 & 0 \\ 0 & 2 \end{pmatrix}\), \(\det = 2\), inverse \(\begin{pmatrix} 1 & 0 \\ 0 & \tfrac12 \end{pmatrix}\). Placing it back:

\(GSG = G\) also holds, so \(G\) is reflexive. But \((SG)' = \begin{pmatrix} 1 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix} \ne SG\): condition (iii) fails.

ginv_sub <- function(S, rows, cols) {   # a g-inverse from a non-singular submatrix
  G <- matrix(0, ncol(S), nrow(S))
  G[cols, rows] <- solve(S[rows, cols])
  G
}
S <- matrix(c(1,1,0, 1,1,0, 0,0,2), 3, 3, byrow = TRUE)
G <- ginv_sub(S, c(1, 3), c(1, 3)); G
max(abs(S %*% G %*% S - S))          # 0 -- condition (i) holds
max(abs(G %*% S %*% G - G))          # 0 -- and (ii): G is reflexive
isTRUE(all.equal(t(S %*% G), S %*% G))   # FALSE -- so NOT Moore-Penrose

     [,1] [,2] [,3]
[1,]    1    0  0.0
[2,]    0    0  0.0
[3,]    0    0  0.5
[1] 0
[1] 0
[1] FALSE

5. Result

\(G\) is a (reflexive) generalized inverse of \(S\) but not its Moore–Penrose inverse. There are many g-inverses — a different submatrix gives a different \(G\), and every one satisfies \(SGS = S\); the Moore–Penrose inverse is the single one that is also symmetric in both products. In Unit 4 that non-uniqueness is why \(\hat{\boldsymbol\beta}\) is not unique in a rank-deficient model while \(\boldsymbol\ell'\hat{\boldsymbol\beta}\) for estimable \(\boldsymbol\ell\) is.

Practical 6 — Characteristic Equation from Traces of Successive Powers

1. Problem

Find the characteristic equation of \(A = \begin{pmatrix} 2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{pmatrix}\) from the traces of its powers, its roots, and its eigenvectors.

2. Aim

To obtain the characteristic equation of a matrix from the traces of its successive powers, and its roots and vectors.

3. Formula

Checks: \(c_3 = \det A\); the roots sum to the trace and multiply to the determinant.

Applying it:

1. Compute \(A^{2}\), \(A^{3}\) and the three traces.
2. The coefficients; the equation.
3. Factorise; check.
4. The eigenvectors from \((A - \lambda I)\mathbf{x} = \mathbf{0}\).

4. Calculation

\(c_3 = 4 = \det A\); \(4 + 1 + 1 = 6\) and \(4 \times 1 \times 1 = 4\). \(\checkmark\) For \(\lambda = 4\): \(A - 4I\) has rank 2 and \(\mathbf{x} = (1, 1, 1)'\). For \(\lambda = 1\): \(A - I\) has rank 1, the single equation \(x_1 + x_2 + x_3 = 0\), and a basis \((-1, 1, 0)'\), \((-1, 0, 1)'\), both orthogonal to \((1, 1, 1)'\).

A <- matrix(c(2,1,1, 1,2,1, 1,1,2), 3, 3, byrow = TRUE)
tr <- function(M) sum(diag(M))       # characteristic equation from traces of powers
A2 <- A %*% A; A3 <- A2 %*% A
t1 <- tr(A); t2 <- tr(A2); t3 <- tr(A3)
c1 <- t1; c2 <- (c1*t1 - t2)/2; c3 <- (c2*t1 - c1*t2 + t3)/3
c(t1, t2, t3); c(c1, c2, c3)
Re(polyroot(c(-c3, c2, -c1, 1)))     # the roots

[1]  6 18 66
[1] 6 9 4
[1] 1 1 4

5. Result

The characteristic equation is \(\lambda^{3} - 6\lambda^{2} + 9\lambda - 4 = 0\), with roots 4 (vector \((1,1,1)'\)) and 1 twice (any vector with \(x_1 + x_2 + x_3 = 0\)) (worked also in Unit 2, Example 2.1).

Practical 7 — Spectral Decomposition of a Third-Order Matrix

1. Problem

Find the spectral decomposition of \(A = \begin{pmatrix} 2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{pmatrix}\) (roots 4 and 1, twice, from Practical 6), and from it \(A^{-1}\) and \(A^{1/2}\).

2. Aim

To write a symmetric matrix as a sum of its eigenvalues times the projections onto its eigenspaces, and use it for functions of the matrix.

3. Formula

Applying it:

1. The projection onto the \(\lambda = 4\) eigenspace, \(P_1 = \mathbf{v}\mathbf{v}'/\mathbf{v}'\mathbf{v}\).
2. \(P_2 = I - P_1\).
3. Check idempotence, orthogonality and the traces (the multiplicities).
4. Reassemble \(A\); then \(m = -1\) and \(m = \tfrac12\).

4. Calculation

\((I + J/3)^{2} = I + 2J/3 + J/3 = I + J = A\). \(\checkmark\)

A <- matrix(c(2,1,1, 1,2,1, 1,1,2), 3, 3, byrow = TRUE)
J <- matrix(1, 3, 3); P1 <- J / 3; P2 <- diag(3) - J / 3
round(max(abs(4 * P1 + P2 - A)), 10)   # 0 -- the decomposition reproduces A
round(c(max(abs(P1 %*% P1 - P1)), max(abs(P1 %*% P2))), 10)   # idempotent, orthogonal
round(max(abs((P1 / 4 + P2) - solve(A))), 10)   # 0 -- the inverse
R <- 2 * P1 + P2; R; round(max(abs(R %*% R - A)), 10)   # the square root, and its square

[1] 0
[1] 0 0
[1] 0
          [,1]      [,2]      [,3]
[1,] 1.3333333 0.3333333 0.3333333
[2,] 0.3333333 1.3333333 0.3333333
[3,] 0.3333333 0.3333333 1.3333333
[1] 0

5. Result

\(A = 4\cdot\frac{J}{3} + 1\cdot\left(I - \frac{J}{3}\right)\); hence \(A^{-1} = I - J/4\) (diagonal \(\tfrac34\), off-diagonal \(-\tfrac14\), as in Practical 1) and \(A^{1/2} = I + J/3\). Every function of \(A\) comes from one decomposition (worked also in Unit 2, Example 2.4).

Practical 8 — Simultaneous Reduction of a Pair of Quadratic Forms

1. Problem

Reduce the quadratic forms \(\mathbf{x}'A\mathbf{x}\) and \(\mathbf{x}'B\mathbf{x}\) simultaneously, where \(A = \begin{pmatrix} 5 & 2 \\ 2 & 2 \end{pmatrix}\) and \(B = \begin{pmatrix} 2 & 1 \\ 1 & 1 \end{pmatrix}\).

2. Aim

To reduce two quadratic forms, one positive definite, to sums of squares by one transformation.

3. Formula

\(B\) must be positive definite.

Applying it:

1. Check \(B\) (and \(A\)) positive definite by leading minors.
2. Form \(|A - \lambda B|\); factor out the common term.
3. Verify each root; write the canonical forms.

4. Calculation

\(B\): \(2 > 0\), \(\det B = 1 > 0\); \(A\): \(5 > 0\), \(\det A = 6 > 0\). \(\checkmark\)

A <- matrix(c(5, 2, 2, 2), 2); B <- matrix(c(2, 1, 1, 1), 2)
c(det(B), det(A))                    # both positive definite
eigen(solve(B) %*% A)$values         # the roots of |A - lambda B| = 0
c(det(A - 2 * B), det(A - 3 * B))    # both 0

[1] 1 6
[1] 3 2
[1] 0 0

5. Result

The roots are \(\lambda = 2, 3\), and the pair reduces to \(\mathbf{x}'B\mathbf{x} = y_1^{2} + y_2^{2}\), \(\mathbf{x}'A\mathbf{x} = 2y_1^{2} + 3y_2^{2}\). The ratio \(\mathbf{x}'A\mathbf{x}/\mathbf{x}'B\mathbf{x}\) therefore ranges over \([2, 3]\) and no further — in a discriminant problem, the smallest and largest achievable separation (worked also in Unit 3, Example 3.3).

Practical 9 — Orthonormal Basis by the Gram–Schmidt Process

1. Problem

Construct an orthonormal basis of \(\mathbb{R}^{3}\) from \(\mathbf{v}_1 = (1,1,0)'\), \(\mathbf{v}_2 = (1,0,1)'\), \(\mathbf{v}_3 = (0,1,1)'\).

2. Aim

To turn a set of independent vectors into an orthonormal basis of the same space.

3. Formula

Applying it:

1. Take \(\mathbf{u}_1 = \mathbf{v}_1\).
2. Subtract from \(\mathbf{v}_2\) its projection on \(\mathbf{u}_1\); check orthogonality.
3. Subtract from \(\mathbf{v}_3\) its projections on \(\mathbf{u}_1\) and \(\mathbf{u}_2\); check.
4. Normalise.

4. Calculation

gram_schmidt <- function(V) {          # columns of V in, orthonormal columns out
  E <- V
  for (k in seq_len(ncol(V))) {
    u <- V[, k]
    for (j in seq_len(k - 1)) u <- u - sum(V[, k] * E[, j]) * E[, j]
    E[, k] <- u / sqrt(sum(u^2))
  }
  E
}
V <- cbind(c(1, 1, 0), c(1, 0, 1), c(0, 1, 1))
E <- gram_schmidt(V); round(E, 6)
round(t(E) %*% E, 10)                # the identity: orthonormal

         [,1]      [,2]     [,3]
[1,] 0.707107  0.408248 -0.57735
[2,] 0.707107 -0.408248  0.57735
[3,] 0.000000  0.816497  0.57735
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

5. Result

\(\{\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3\}\) above is an orthonormal basis: \(E'E = I\). The starting vectors were independent but at \(60^{\circ}\) to one another; the output is a genuine set of coordinate axes for the same space — the construction that, as the Helmert matrix, makes the sample mean and variance independent (worked also in Unit 1, Example 1.1).

Practical 10 — The Variance–Covariance Matrix and Its Characteristics

1. Problem

Five units are measured on three variables:

Find the variance–covariance matrix, its trace and determinant, and the correlation matrix, and interpret them.

2. Aim

To compute the variance–covariance matrix of a data set and read the data's spread from its trace, determinant and correlations.

3. Formula

Applying it:

1. Compute each entry from the deviations.
2. Trace (total variance) and determinant (generalized variance); compare with \(\prod s_{jj}\).
3. Divide by the standard deviations for the correlation matrix.

4. Calculation

X <- cbind(c(2,4,6,8,10), c(1,3,2,5,4), c(5,4,6,5,8))
S <- cov(X); S
c(trace = sum(diag(S)), det = det(S), prod_diag = prod(diag(S)))
round(cov2cor(S), 6)

     [,1] [,2] [,3]
[1,] 10.0  4.0  3.5
[2,]  4.0  2.5  0.5
[3,]  3.5  0.5  2.3
    trace       det prod_diag 
   14.800     1.575    57.500 
       [,1]     [,2]     [,3]
[1,] 1.0000 0.800000 0.729800
[2,] 0.8000 1.000000 0.208514
[3,] 0.7298 0.208514 1.000000

5. Result

Total variance 14.8; generalized variance 1.575, against 57.5 for uncorrelated variables with these variances (Hadamard's inequality, \(1.575 \le 57.5\)): the correlations have collapsed the cloud to under 3% of that volume. \(x_1\) is strongly related to both others, but \(x_2\) and \(x_3\) are almost unrelated to each other (\(r = 0.2085\)): their apparent connection runs through \(x_1\), which is what the partial correlations of Practical 13 are designed to expose.

Practical 11 — Simple Linear Regression with Lack of Fit

1. Problem

Eight observations at four distinct \(x\) values, two replicates each:

Fit the simple linear regression; find \(R^{2}\), adjusted \(R^{2}\), a 95% confidence interval for the slope, and test for lack of fit using the pure error.

2. Aim

To fit a simple linear regression, measure its fit, and test whether a straight line is adequate by separating pure error from lack of fit.

3. Formula

Applying it:

1. The raw sums and the corrected sums of squares and products.
2. The coefficients; the analysis of variance; \(R^{2}\) and \(\bar R^{2}\).
3. The confidence interval for the slope.
4. Pure error from the replicates; lack of fit by subtraction; the \(F\) test.

4. Calculation

x <- c(1,1,2,2,3,3,4,4); y <- c(3,4,6,7,8,10,11,13)
fit <- lm(y ~ x); coef(fit)
c(R2 = summary(fit)$r.squared, adjR2 = summary(fit)$adj.r.squared)
confint(fit, "x")
anova(fit, lm(y ~ factor(x)))       # lack of fit: the line against the group means

(Intercept)           x 
       0.75        2.80 
       R2     adjR2 
0.9389222 0.9287425 
     2.5 %   97.5 %
x 2.086609 3.513391
Analysis of Variance Table

Model 1: y ~ x
Model 2: y ~ factor(x)
  Res.Df RSS Df Sum of Sq    F Pr(>F)
1      6 5.1                         
2      4 5.0  2       0.1 0.04 0.9612

5. Result

\(\hat y = 0.75 + 2.8x\), \(R^{2} = 0.939\), \(\bar R^{2} = 0.929\), and with 95% confidence the slope lies between 2.09 and 3.51. The lack-of-fit \(F = 0.04\) (\(p = 0.96\)) gives no evidence against a straight line. Without replication \(SSE\) could not have been split, and the residual would confound curvature with noise.

Practical 12 — Multiple Linear Regression

1. Problem

Regress \(y = 12, 18, 20, 28, 32\) on \(x_1\) and \(x_2\) of Practical 10 (\(x_1 = 2, 4, 6, 8, 10\); \(x_2 = 1, 3, 2, 5, 4\)). Find the estimates, \(R^{2}\), adjusted \(R^{2}\), the overall \(F\) test and the standard errors.

2. Aim

To fit a multiple linear regression by solving the normal equations in matrix form, and to test the model and each coefficient.

3. Formula

Applying it:

1. Form \(X'X\) and \(X'\mathbf{y}\); check \(\det(X'X) \ne 0\).
2. Solve for \(\hat{\boldsymbol\beta}\).
3. The analysis of variance, \(R^{2}\), \(\bar R^{2}\), \(F\).
4. Standard errors and \(t\) from the diagonal of \((X'X)^{-1}\).

4. Calculation

x1 <- c(2,4,6,8,10); x2 <- c(1,3,2,5,4); y <- c(12,18,20,28,32)
X <- cbind(1, x1, x2); XtX <- t(X) %*% X; XtX; det(XtX)
b <- solve(XtX, t(X) %*% y); round(drop(b), 6)
fit <- lm(y ~ x1 + x2); s <- summary(fit)
round(c(R2 = s$r.squared, adjR2 = s$adj.r.squared, F = unname(s$fstatistic[1])), 6)
round(coef(s)[-1, ], 6)

       x1  x2
    5  30  15
x1 30 220 106
x2 15 106  55
[1] 720
               x1       x2 
6.333333 2.055556 1.111111 
        R2      adjR2          F 
  0.993924   0.987847 163.571429 
   Estimate Std. Error  t value Pr(>|t|)
x1 2.055556   0.232406 8.844692 0.012543
x2 1.111111   0.464811 2.390457 0.139337

5. Result

\(\hat y = 6.333 + 2.056x_1 + 1.111x_2\), \(R^{2} = 0.994\), \(\bar R^{2} = 0.988\). The model as a whole is highly significant (\(p = 0.006\)) and \(x_1\) is individually significant (\(p = 0.013\)), but \(x_2\) is not (\(p = 0.139\)) — with only two error degrees of freedom there is very little power, the classic sign of a model over-fitted to its sample size, and why \(\bar R^{2}\) is reported beside \(R^{2}\).

Practical 13 — Simple, Partial and Multiple Correlation

1. Problem

For the data of Practical 12, find the simple correlations of \(y\), \(x_1\), \(x_2\), the partial correlations \(r_{01.2}\) and \(r_{02.1}\), and the multiple correlation \(R_{0.12}\).

2. Aim

To measure the relation between two variables with a third held fixed, and of one variable with two others jointly.

3. Formula

Subscript 0 is \(y\). Check: \(R_{0.12}^{2}\) equals the \(R^{2}\) of the regression of Practical 12.

Applying it:

1. The three simple correlations.
2. The two partial correlations.
3. The multiple correlation; cross-check with Practical 12.

4. Calculation

Numerator: \(0.988212 - 0.869626 \times 0.8 = 0.292511\). Denominator: \(\sqrt{(1 - 0.756249)(1 - 0.64)} = \sqrt{0.087750} = 0.296226\).

Check: \(\sqrt{R^{2}} = \sqrt{0.993924} = 0.996957\) (Practical 12). \(\checkmark\)

x1 <- c(2,4,6,8,10); x2 <- c(1,3,2,5,4); y <- c(12,18,20,28,32)
r01 <- cor(y, x1); r02 <- cor(y, x2); r12 <- cor(x1, x2)
round(c(r01 = r01, r02 = r02, r12 = r12), 6)
p012 <- (r01 - r02 * r12) / sqrt((1 - r02^2) * (1 - r12^2))
p021 <- (r02 - r01 * r12) / sqrt((1 - r01^2) * (1 - r12^2))
R012 <- sqrt((r01^2 + r02^2 - 2 * r01 * r02 * r12) / (1 - r12^2))
round(c(r01.2 = p012, r02.1 = p021, R0.12 = R012, sqrtR2 = sqrt(summary(lm(y ~ x1 + x2))$r.squared)), 6)

     r01      r02      r12 
0.988212 0.869626 0.800000 
   r01.2    r02.1    R0.12   sqrtR2 
0.987457 0.860663 0.996957 0.996957 

5. Result

Holding \(x_2\) fixed barely changes the relation between \(y\) and \(x_1\) (0.9882 → 0.9875), so \(x_1\) carries information of its own; holding \(x_1\) fixed also leaves \(r(y, x_2)\) high (0.8696 → 0.8607). The multiple correlation is 0.997. Contrast Practical 10, where \(x_2\) and \(x_3\) appeared related only through \(x_1\) — the case where a partial correlation collapses towards zero, and what the technique is for.

Practical 14 — Testing Multicollinearity

1. Problem

Three regressors are observed on five units:

Test the design for multicollinearity.

2. Aim

To detect multicollinearity among regressors by the correlation matrix, the variance inflation factors and the condition number.

3. Formula

\(R_j^{2}\) comes from regressing \(x_j\) on the other regressors. The usual warning lines are \(VIF > 10\) and \(\kappa\) above about 30.

Applying it:

1. Look before computing: \(x_2\) is exactly \(2x_1\) at the first four units and 11 rather than 10 at the fifth.
2. The correlation matrix — cheap, but blind to a relation among three or more regressors.
3. The variance inflation factors, which do catch multi-way relations.
4. The condition number, a single number for the whole design.

4. Calculation

The eigenvalues of the correlation matrix are 2.714426, 0.282689 and 0.002884, so \(\kappa = \sqrt{2.714426/0.002884} = 30.68\).

Z <- cbind(x1 = c(1,2,3,4,5), x2 = c(2,4,6,8,11), x3 = c(1,3,2,5,4))
round(cor(Z), 6)
vif_manual <- function(X) {          # VIF without the car package
  sapply(seq_len(ncol(X)), function(j)
    1 / (1 - summary(lm(X[, j] ~ X[, -j]))$r.squared))
}
round(vif_manual(Z), 4)
ev <- eigen(cor(Z))$values; round(ev, 6)
sqrt(max(ev) / min(ev))              # the condition number

         x1       x2       x3
x1 1.000000 0.995893 0.800000
x2 0.995893 1.000000 0.769554
x3 0.800000 0.769554 1.000000
[1] 184.2593 162.6667   3.7037
[1] 2.714426 0.282689 0.002884
[1] 30.67827

5. Result

\(x_1\) and \(x_2\) are severely collinear (\(VIF\) 184 and 163, far above 10); \(x_3\) is not part of the problem (\(VIF = 3.7\)). The condition number of the correlation matrix, 30.7, is just above the warning line of 30: one eigenvalue, 0.0029, is nearly zero — the near-exact relation \(x_2 \approx 2x_1\). The standard error of \(\hat\beta_1\) is \(\sqrt{184.26} = 13.6\) times what orthogonal regressors would give. The estimates stay unbiased and the fit stays good inside the data's range: multicollinearity damages the separation of effects, not the fit. The remedy is in the design — drop or combine one of the pair, or collect data that breaks the pattern — not a larger sample of the same shape (worked also in Unit 4, Example 4.2).

What the Practical Record Should Contain