The goal of this chapter 5 tutorial is to express the matrix arithmetic in R. Many of the operations in multiple regression analysis have R functions fully capable of performing these operations internally, releaving the analyst to focus on the details. For pedagogical reasons, however, this tutorial will follow the examples in the book, so that if the interested reader needed to perform the requisite matrix operations, the details of those methods will be laid out.
# The data set (from chapter 6) to be used in this chapter
df <- structure(list(x1 = c(68.5, 45.2, 91.3, 47.8, 46.9, 66.1, 49.5,
52, 48.9, 38.4, 87.9, 72.8, 88.4, 42.9, 52.5, 85.7, 41.3, 51.7,
89.6, 82.7, 52.3), x2 = c(16.7, 16.8, 18.2, 16.3, 17.3, 18.2,
15.9, 17.2, 16.6, 16, 18.3, 17.1, 17.4, 15.8, 17.8, 18.4, 16.5,
16.3, 18.1, 19.1, 16), y = c(174.4, 164.4, 244.2, 154.6, 181.6,
207.5, 152.8, 163.2, 145.4, 137.2, 241.9, 191.1, 232, 145.3,
161.1, 209.7, 146.4, 144, 232.6, 224.1, 166.5)), .Names = c("x1",
"x2", "y"), class = "data.frame", row.names = c(NA, -21L))
The linear model object is a simple interface for the R user, but there are many facilities to interact with the matricies involved in the regression examples. Some of these will be detailed below.
(A = c(4, 7, 10)) # Vector in R is not the type of entities we want.
## [1] 4 7 10
as.matrix(A) # Matrix of R vector is our basic column vector.
## [,1]
## [1,] 4
## [2,] 7
## [3,] 10
t(A) # Transpose of a matrix. Not the same as R vector! A is first coerced to column matrix.
## [,1] [,2] [,3]
## [1,] 4 7 10
fit <- lm(y ~ x1 + x2, data = df)
(X = model.matrix(fit)) # (5.6) Note: carries an attribute "assign"
## (Intercept) x1 x2
## 1 1 68.5 16.7
## 2 1 45.2 16.8
## 3 1 91.3 18.2
## 4 1 47.8 16.3
## 5 1 46.9 17.3
## 6 1 66.1 18.2
## 7 1 49.5 15.9
## 8 1 52.0 17.2
## 9 1 48.9 16.6
## 10 1 38.4 16.0
## 11 1 87.9 18.3
## 12 1 72.8 17.1
## 13 1 88.4 17.4
## 14 1 42.9 15.8
## 15 1 52.5 17.8
## 16 1 85.7 18.4
## 17 1 41.3 16.5
## 18 1 51.7 16.3
## 19 1 89.6 18.1
## 20 1 82.7 19.1
## 21 1 52.3 16.0
## attr(,"assign")
## [1] 0 1 2
t(X) # (5.7)
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
## (Intercept) 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
## x1 68.5 45.2 91.3 47.8 46.9 66.1 49.5 52.0 48.9 38.4 87.9 72.8 88.4 42.9 52.5 85.7 41.3 51.7 89.6 82.7 52.3
## x2 16.7 16.8 18.2 16.3 17.3 18.2 15.9 17.2 16.6 16.0 18.3 17.1 17.4 15.8 17.8 18.4 16.5 16.3 18.1 19.1 16.0
## attr(,"assign")
## [1] 0 1 2
(Y = model.frame(fit)[1]) # (5.4) Note: this returns a data.frame object
## y
## 1 174.4
## 2 164.4
## 3 244.2
## 4 154.6
## 5 181.6
## 6 207.5
## 7 152.8
## 8 163.2
## 9 145.4
## 10 137.2
## 11 241.9
## 12 191.1
## 13 232.0
## 14 145.3
## 15 161.1
## 16 209.7
## 17 146.4
## 18 144.0
## 19 232.6
## 20 224.1
## 21 166.5
t(Y) # (5.5) Note: Can still transpose frame.
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## y 174.4 164.4 244.2 154.6 181.6 207.5 152.8 163.2 145.4 137.2 241.9 191.1 232 145.3 161.1 209.7 146.4 144 232.6 224.1
## 21
## y 166.5
(e = as.matrix(residuals(fit))) # (5.10)
## [,1]
## 1 -12.7841
## 2 10.1706
## 3 9.8037
## 4 1.2715
## 5 20.2151
## 6 9.7586
## 7 0.7449
## 8 -4.6666
## 9 -12.3382
## 10 0.3540
## 11 11.5126
## 12 -6.0849
## 13 9.3143
## 14 3.7816
## 15 -13.1132
## 16 -18.4239
## 17 0.6530
## 18 -15.0013
## 19 1.6130
## 20 -6.2161
## 21 9.4356
(A = matrix(1:6, ncol = 2)) # 3x2 matrix
## [,1] [,2]
## [1,] 1 4
## [2,] 2 5
## [3,] 3 6
(B = matrix(c(1:3, 2:4), ncol = 2)) # 3x2 matrix
## [,1] [,2]
## [1,] 1 2
## [2,] 2 3
## [3,] 3 4
A + B
## [,1] [,2]
## [1,] 2 6
## [2,] 4 8
## [3,] 6 10
A - B # 3x2 (5.8)
## [,1] [,2]
## [1,] 0 2
## [2,] 0 2
## [3,] 0 2
(A = matrix(c(2, 9, 7, 3), ncol = 2)) # 2x2
## [,1] [,2]
## [1,] 2 7
## [2,] 9 3
4 * A # 2x2 (5.11)
## [,1] [,2]
## [1,] 8 28
## [2,] 36 12
(A = matrix(c(2, 4, 5, 1), ncol = 2)) # 2x2
## [,1] [,2]
## [1,] 2 5
## [2,] 4 1
(B = matrix(c(4, 5, 6, 8), ncol = 2)) # 2x2
## [,1] [,2]
## [1,] 4 6
## [2,] 5 8
A %*% B # 2x2
## [,1] [,2]
## [1,] 33 52
## [2,] 21 32
(A = matrix(c(1, 3, 4, 0, 5, 8), ncol = 3, byrow = TRUE)) ## 2x3
## [,1] [,2] [,3]
## [1,] 1 3 4
## [2,] 0 5 8
(B = matrix(c(3, 5, 2))) # 3x1
## [,1]
## [1,] 3
## [2,] 5
## [3,] 2
A %*% B # 2x1 (= 2x3 %*% 3x1)
## [,1]
## [1,] 26
## [2,] 41
(beta = matrix(coef(fit))) # coefficients from model
## [,1]
## [1,] -68.857
## [2,] 1.455
## [3,] 9.366
cbind (
X %*% beta, # fitted values
fitted(fit) # wrapper to get fitted values from object
);
## [,1] [,2]
## 1 187.2 187.2
## 2 154.2 154.2
## 3 234.4 234.4
## 4 153.3 153.3
## 5 161.4 161.4
## 6 197.7 197.7
## 7 152.1 152.1
## 8 167.9 167.9
## 9 157.7 157.7
## 10 136.8 136.8
## 11 230.4 230.4
## 12 197.2 197.2
## 13 222.7 222.7
## 14 141.5 141.5
## 15 174.2 174.2
## 16 228.1 228.1
## 17 145.7 145.7
## 18 159.0 159.0
## 19 231.0 231.0
## 20 230.3 230.3
## 21 157.1 157.1
class(Y)
## [1] "data.frame"
class(t(Y)) # See difference
## [1] "matrix"
t(Y) %*% Y # Error produced, Y not a matrix
## Error: requires numeric/complex matrix/vector arguments
t(Y) %*% as.matrix(Y) # (5.13) SSy = sum of squares of Y
## y
## y 721072
t(X) %*% X # (5.14)
## (Intercept) x1 x2
## (Intercept) 21 1302 360
## x1 1302 87708 22609
## x2 360 22609 6190
t(X) %*% as.matrix(Y) # (5.15)
## y
## (Intercept) 3820
## x1 249643
## x2 66073
solve(t(X) %*% X) # (5.24) See Inverse below
## (Intercept) x1 x2
## (Intercept) 29.72892 0.0721835 -1.99255
## x1 0.07218 0.0003702 -0.00555
## x2 -1.99255 -0.0055499 0.13631
t( t(X) %*% X ) # Symmetric matrix A = t(A)
## (Intercept) x1 x2
## (Intercept) 21 1302 360
## x1 1302 87708 22609
## x2 360 22609 6190
(I = diag(3)) # Diagonal matrix of size 3
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1
I %*% (t(X) %*% X) # Also defines an identity matrix
## (Intercept) x1 x2
## [1,] 21 1302 360
## [2,] 1302 87708 22609
## [3,] 360 22609 6190
2*I # Scalar matrix
## [,1] [,2] [,3]
## [1,] 2 0 0
## [2,] 0 2 0
## [3,] 0 0 2
(J = matrix(1, 4, 4)) # (5.18)
## [,1] [,2] [,3] [,4]
## [1,] 1 1 1 1
## [2,] 1 1 1 1
## [3,] 1 1 1 1
## [4,] 1 1 1 1
matrix(1, 4) %*% t(matrix(1, 4)) # Also equals J
## [,1] [,2] [,3] [,4]
## [1,] 1 1 1 1
## [2,] 1 1 1 1
## [3,] 1 1 1 1
## [4,] 1 1 1 1
(A = matrix(c(1, 2, 5, 1, 2, 2, 10, 6, 3, 4, 15, 1), ncol = 4, byrow = T))
## [,1] [,2] [,3] [,4]
## [1,] 1 2 5 1
## [2,] 2 2 10 6
## [3,] 3 4 15 1
qr(A)$rank # QR decomposition returns rank value
## [1] 3
5*A[,1] + 0*A[,2] - 1*A[,3] + 0*A[,4] # equals the 0 vector)
## [1] 0 0 0
(A = matrix(c(2, 3, 4, 1), ncol = 2))
## [,1] [,2]
## [1,] 2 4
## [2,] 3 1
solve(A) # Inverse of A
## [,1] [,2]
## [1,] -0.1 0.4
## [2,] 0.3 -0.2
round(solve(A) %*% A) # Round to remove unnecessary precision
## [,1] [,2]
## [1,] 1 0
## [2,] 0 1
A %*% solve(A) # Rounding unnecssary; first term
## [,1] [,2]
## [1,] 1 0
## [2,] 0 1
# in product determines precision?
det(A) # (5.23b)
## [1] -10
(A = matrix(ncol = 2, sample(seq(100), 4))) # Randomly defined 2x2 matrices
## [,1] [,2]
## [1,] 44 78
## [2,] 8 25
(B = matrix(ncol = 2, sample(seq(100), 4))) # from 0 to 100 with equal prob.
## [,1] [,2]
## [1,] 50 39
## [2,] 56 6
(C = matrix(ncol = 2, sample(seq(100), 4)))
## [,1] [,2]
## [1,] 83 84
## [2,] 72 81
(k = sample(0:100, 1))
## [1] 16
A + B == B + A # (5.25)
## [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUE
(A + B) + C == A + (B + C) # (5.26)
## [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUE
(A %*% B) %*% C == A %*% (B %*% C) # (5.27)
## [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUE
C %*% (A + B) == C %*% A + C %*% B # (5.28)
## [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUE
k * (A + B) == k*A + k*B # (5.29)
## [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUE
t(t(A)) == A # (5.30)
## [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUE
t(A + B) == t(A) + t(B) # (5.31)
## [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUE
t(A %*% B) == t(B) %*% t(A) # (5.32)
## [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUE
t(A %*% B %*% C) == t(C) %*% t(B) %*% t(A) # (5.33)
## [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUE
solve(A %*% B)
## [,1] [,2]
## [1,] -0.0005152 0.002435
## [2,] 0.0020072 -0.007324
solve(B) %*% solve(A) # (5.34) Check manually, precision causes FALSE comparisons.
## [,1] [,2]
## [1,] -0.0005152 0.002435
## [2,] 0.0020072 -0.007324
solve(A %*% B %*% C)
## [,1] [,2]
## [1,] -0.0003116 0.001204
## [2,] 0.0003018 -0.001160
solve(C) %*% solve(B) %*% solve(A) # (5.33) Same here.
## [,1] [,2]
## [1,] -0.0003116 0.001204
## [2,] 0.0003018 -0.001160
solve(solve(A))
## [,1] [,2]
## [1,] 44 78
## [2,] 8 25
A # (5.36) and here.
## [,1] [,2]
## [1,] 44 78
## [2,] 8 25
solve(t(A))
## [,1] [,2]
## [1,] 0.05252 -0.01681
## [2,] -0.16387 0.09244
t(solve(A)) # (5.37) and here.
## [,1] [,2]
## [1,] 0.05252 -0.01681
## [2,] -0.16387 0.09244
options(scipen = 3, digits = 4) # Ease up on scientific notation
Y = as.matrix(Y) # Recall Y = model.frame(fit)[1]
vcov(fit) # (5.40) see methods(class = "lm")
## (Intercept) x1 x2
## (Intercept) 3602.035 8.74594 -241.4230
## x1 8.746 0.04485 -0.6724
## x2 -241.423 -0.67244 16.5158
diag(vcov(fit)) # Variances of coefficients
## (Intercept) x1 x2
## 3602.03467 0.04485 16.51576
cbind(Y, X, e)
## y (Intercept) x1 x2
## 1 174.4 1 68.5 16.7 -12.7841
## 2 164.4 1 45.2 16.8 10.1706
## 3 244.2 1 91.3 18.2 9.8037
## 4 154.6 1 47.8 16.3 1.2715
## 5 181.6 1 46.9 17.3 20.2151
## 6 207.5 1 66.1 18.2 9.7586
## 7 152.8 1 49.5 15.9 0.7449
## 8 163.2 1 52.0 17.2 -4.6666
## 9 145.4 1 48.9 16.6 -12.3382
## 10 137.2 1 38.4 16.0 0.3540
## 11 241.9 1 87.9 18.3 11.5126
## 12 191.1 1 72.8 17.1 -6.0849
## 13 232.0 1 88.4 17.4 9.3143
## 14 145.3 1 42.9 15.8 3.7816
## 15 161.1 1 52.5 17.8 -13.1132
## 16 209.7 1 85.7 18.4 -18.4239
## 17 146.4 1 41.3 16.5 0.6530
## 18 144.0 1 51.7 16.3 -15.0013
## 19 232.6 1 89.6 18.1 1.6130
## 20 224.1 1 82.7 19.1 -6.2161
## 21 166.5 1 52.3 16.0 9.4356
c(beta = beta) # (5.52) R indexes by 1, not 0
## beta1 beta2 beta3
## -68.857 1.455 9.366
cbind(X %*% beta + e, Y) # (5.53)
## y
## 1 174.4 174.4
## 2 164.4 164.4
## 3 244.2 244.2
## 4 154.6 154.6
## 5 181.6 181.6
## 6 207.5 207.5
## 7 152.8 152.8
## 8 163.2 163.2
## 9 145.4 145.4
## 10 137.2 137.2
## 11 241.9 241.9
## 12 191.1 191.1
## 13 232.0 232.0
## 14 145.3 145.3
## 15 161.1 161.1
## 16 209.7 209.7
## 17 146.4 146.4
## 18 144.0 144.0
## 19 232.6 232.6
## 20 224.1 224.1
## 21 166.5 166.5
cbind(t(X) %*% X %*% beta, t(X) %*% Y) # (5.59)
## y
## (Intercept) 3820 3820
## x1 249643 249643
## x2 66073 66073
cbind(beta,solve(t(X) %*% X) %*% t(X) %*% Y) # (5.60)
## y
## (Intercept) -68.857 -68.857
## x1 1.455 1.455
## x2 9.366 9.366
fitted(fit) # (5.70)
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 187.2 154.2 234.4 153.3 161.4 197.7 152.1 167.9 157.7 136.8 230.4 197.2 222.7 141.5 174.2 228.1 145.7 159.0 231.0 230.3
## 21
## 157.1
cbind(fitted(fit), X %*% beta) # (5.71)
## [,1] [,2]
## 1 187.2 187.2
## 2 154.2 154.2
## 3 234.4 234.4
## 4 153.3 153.3
## 5 161.4 161.4
## 6 197.7 197.7
## 7 152.1 152.1
## 8 167.9 167.9
## 9 157.7 157.7
## 10 136.8 136.8
## 11 230.4 230.4
## 12 197.2 197.2
## 13 222.7 222.7
## 14 141.5 141.5
## 15 174.2 174.2
## 16 228.1 228.1
## 17 145.7 145.7
## 18 159.0 159.0
## 19 231.0 231.0
## 20 230.3 230.3
## 21 157.1 157.1
(H = X %*% solve(t(X) %*% X) %*% t(X)) # (5.73a) Big! 21x21 matrix
## 1 2 3 4 5 6 7 8 9 10 11
## 1 0.121759 -0.0010405 0.087994 0.059755 -0.040951 -0.034400 0.106546 -0.006547 0.036197 0.04300 0.0618474
## 2 -0.001041 0.1043505 -0.029695 0.069806 0.120306 0.079249 0.043813 0.093597 0.079033 0.09646 -0.0103354
## 3 0.087994 -0.0296952 0.173747 -0.007564 -0.030446 0.048451 0.008251 -0.003248 -0.007617 -0.04878 0.1550009
## 4 0.059755 0.0698056 -0.007564 0.086271 0.050822 0.007197 0.099665 0.051432 0.074833 0.10257 -0.0091703
## 5 -0.040951 0.1203057 -0.030446 0.050822 0.161974 0.132569 -0.002307 0.118450 0.075305 0.08003 0.0020808
## 6 -0.034400 0.0792490 0.048451 0.007197 0.132569 0.158232 -0.048789 0.098206 0.038840 0.01171 0.0754077
## 7 0.106546 0.0438125 0.008251 0.099665 -0.002307 -0.048789 0.143487 0.019230 0.072174 0.10837 -0.0094384
## 8 -0.006547 0.0935966 -0.003248 0.051432 0.118450 0.098206 0.019230 0.091578 0.066022 0.07026 0.0167796
## 9 0.036197 0.0790329 -0.007617 0.074833 0.075305 0.038840 0.072174 0.066022 0.072449 0.09252 -0.0014679
## 10 0.043000 0.0964607 -0.048779 0.102569 0.080030 0.011712 0.108369 0.070258 0.092519 0.13254 -0.0430880
## 11 0.061847 -0.0103354 0.155001 -0.009170 0.002081 0.075408 -0.009438 0.016780 -0.001468 -0.04309 0.1456715
## 12 0.104110 -0.0009865 0.102011 0.042846 -0.026635 -0.004552 0.076305 0.001499 0.027795 0.02280 0.0810663
## 13 0.150977 -0.0544453 0.174053 0.022915 -0.095950 -0.036076 0.081635 -0.042288 -0.001321 -0.02206 0.1345664
## 14 0.084123 0.0676846 -0.022724 0.107129 0.029854 -0.032182 0.138541 0.039462 0.084461 0.12668 -0.0316936
## 15 -0.061920 0.1194005 -0.011806 0.029554 0.178413 0.168898 -0.039599 0.127601 0.064387 0.05424 0.0268151
## 16 0.041529 0.0038365 0.142222 -0.011479 0.026845 0.097137 -0.024412 0.031976 0.002469 -0.04028 0.1401323
## 17 0.008918 0.1072281 -0.043564 0.082883 0.113935 0.059804 0.064965 0.090280 0.086579 0.11323 -0.0268815
## 18 0.078697 0.0529452 0.011827 0.083986 0.025593 -0.009793 0.108493 0.035731 0.067644 0.09321 0.0031477
## 19 0.089371 -0.0270067 0.167135 -0.002971 -0.029982 0.043712 0.014397 -0.002743 -0.004364 -0.04223 0.1482222
## 20 -0.040476 0.0494325 0.114421 -0.034905 0.119983 0.195222 -0.101157 0.088430 0.007169 -0.05037 0.1405225
## 21 0.110511 0.0363685 0.020332 0.094427 -0.009887 -0.048842 0.139831 0.014297 0.066894 0.09918 0.0008147
## 12 13 14 15 16 17 18 19 20 21
## 1 0.1041098 0.150977 0.084123 -0.061920 0.041529 0.008918 0.078697 0.089371 -0.040476 0.1105115
## 2 -0.0009865 -0.054445 0.067685 0.119401 0.003837 0.107228 0.052945 -0.027007 0.049433 0.0363685
## 3 0.1020114 0.174053 -0.022724 -0.011806 0.142222 -0.043564 0.011827 0.167135 0.114421 0.0203321
## 4 0.0428458 0.022915 0.107129 0.029554 -0.011479 0.082883 0.083986 -0.002971 -0.034905 0.0944274
## 5 -0.0266353 -0.095950 0.029854 0.178413 0.026845 0.113935 0.025593 -0.029982 0.119983 -0.0098868
## 6 -0.0045520 -0.036076 -0.032182 0.168898 0.097137 0.059804 -0.009793 0.043712 0.195222 -0.0488420
## 7 0.0763047 0.081635 0.138541 -0.039599 -0.024412 0.064965 0.108493 0.014397 -0.101157 0.1398314
## 8 0.0014986 -0.042288 0.039462 0.127601 0.031976 0.090280 0.035731 -0.002743 0.088430 0.0142971
## 9 0.0277948 -0.001321 0.084461 0.064387 0.002469 0.086579 0.067644 -0.004364 0.007169 0.0668938
## 10 0.0227985 -0.022058 0.126682 0.054238 -0.040277 0.113232 0.093208 -0.042229 -0.050366 0.0991774
## 11 0.0810663 0.134566 -0.031694 0.026815 0.140132 -0.026882 0.003148 0.148222 0.140523 0.0008147
## 12 0.0960232 0.142288 0.054963 -0.035792 0.065196 0.002224 0.059338 0.101390 0.006537 0.0815775
## 13 0.1422884 0.238960 0.037737 -0.104935 0.105086 -0.053557 0.055435 0.171014 0.002118 0.0938465
## 14 0.0549628 0.037737 0.143758 -0.006512 -0.040213 0.089302 0.108593 -0.015668 -0.095195 0.1318987
## 15 -0.0357922 -0.104935 -0.006512 0.209459 0.056831 0.104645 0.001588 -0.013856 0.178028 -0.0454366
## 16 0.0651956 0.105086 -0.040213 0.056831 0.140190 -0.015121 -0.004502 0.135187 0.162777 -0.0154087
## 17 0.0022241 -0.053557 0.089302 0.104645 -0.015121 0.115017 0.066885 -0.039327 0.016337 0.0562161
## 18 0.0593377 0.055435 0.108593 0.001588 -0.004502 0.066885 0.087332 0.016130 -0.047410 0.1051329
## 19 0.1013901 0.171014 -0.015668 -0.013856 0.135187 -0.039327 0.016130 0.161069 0.104671 0.0258499
## 20 0.0065368 0.002118 -0.095195 0.178028 0.162777 0.016337 -0.047410 0.104671 0.278797 -0.0949343
## 21 0.0815775 0.093846 0.131899 -0.045437 -0.015409 0.056216 0.105133 0.025850 -0.094934 0.1373330
cbind(fitted(fit), H %*% Y) # (5.73)
## y
## 1 187.2 187.2
## 2 154.2 154.2
## 3 234.4 234.4
## 4 153.3 153.3
## 5 161.4 161.4
## 6 197.7 197.7
## 7 152.1 152.1
## 8 167.9 167.9
## 9 157.7 157.7
## 10 136.8 136.8
## 11 230.4 230.4
## 12 197.2 197.2
## 13 222.7 222.7
## 14 141.5 141.5
## 15 174.2 174.2
## 16 228.1 228.1
## 17 145.7 145.7
## 18 159.0 159.0
## 19 231.0 231.0
## 20 230.3 230.3
## 21 157.1 157.1
round(H %*% H - H, 12) == 0 # (5.74) Difference to high precision is 0. H is idempotent
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
## 1 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 2 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 3 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 4 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 5 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 6 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 7 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 8 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 9 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 10 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 11 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 12 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 13 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 14 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 15 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 16 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 17 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 18 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 19 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 20 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 21 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
cbind(hatvalues(fit), diag(H)) # Useful in later chapters.
## [,1] [,2]
## 1 0.12176 0.12176
## 2 0.10435 0.10435
## 3 0.17375 0.17375
## 4 0.08627 0.08627
## 5 0.16197 0.16197
## 6 0.15823 0.15823
## 7 0.14349 0.14349
## 8 0.09158 0.09158
## 9 0.07245 0.07245
## 10 0.13254 0.13254
## 11 0.14567 0.14567
## 12 0.09602 0.09602
## 13 0.23896 0.23896
## 14 0.14376 0.14376
## 15 0.20946 0.20946
## 16 0.14019 0.14019
## 17 0.11502 0.11502
## 18 0.08733 0.08733
## 19 0.16107 0.16107
## 20 0.27880 0.27880
## 21 0.13733 0.13733
cbind(
e, # (5.75)
Y - X %*% beta, # (5.76)
Y - H %*% Y,
(diag(nrow(Y)) - H) %*% Y # (5.78)
);
## y y y
## 1 -12.7841 -12.7841 -12.7841 -12.7841
## 2 10.1706 10.1706 10.1706 10.1706
## 3 9.8037 9.8037 9.8037 9.8037
## 4 1.2715 1.2715 1.2715 1.2715
## 5 20.2151 20.2151 20.2151 20.2151
## 6 9.7586 9.7586 9.7586 9.7586
## 7 0.7449 0.7449 0.7449 0.7449
## 8 -4.6666 -4.6666 -4.6666 -4.6666
## 9 -12.3382 -12.3382 -12.3382 -12.3382
## 10 0.3540 0.3540 0.3540 0.3540
## 11 11.5126 11.5126 11.5126 11.5126
## 12 -6.0849 -6.0849 -6.0849 -6.0849
## 13 9.3143 9.3143 9.3143 9.3143
## 14 3.7816 3.7816 3.7816 3.7816
## 15 -13.1132 -13.1132 -13.1132 -13.1132
## 16 -18.4239 -18.4239 -18.4239 -18.4239
## 17 0.6530 0.6530 0.6530 0.6530
## 18 -15.0013 -15.0013 -15.0013 -15.0013
## 19 1.6130 1.6130 1.6130 1.6130
## 20 -6.2161 -6.2161 -6.2161 -6.2161
## 21 9.4356 9.4356 9.4356 9.4356
The notation can get messy with so many ways to say the same thing. In R anova output, the error sum of squares, sum of squares error, or residual sum of squares is denoted by RSS. The authors use SSE. The authors refer to the regression sum of squares (not to be confused with RSS here!) by SSR. In R anova output it is always in a “Sum Sq” column of the output table.
n = nrow(Y)
J = matrix(1, n, n)
I = diag(n)
anova(lm(y ~ 1, df), fit) # SSTO = SSE for reduced model
## Analysis of Variance Table
##
## Model 1: y ~ 1
## Model 2: y ~ x1 + x2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 20 26196
## 2 18 2181 2 24015 99.1 1.9e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(SSTO = t(Y) %*% Y - (1/n) * t(Y) %*% J %*% Y) # (5.83)
## y
## y 26196
sum(anova(fit)[1:2, "Sum Sq"]) # SSR in R is sum of SSR of all predictors in model
## [1] 24015
anova(fit)
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 23372 23372 192.90 4.6e-11 ***
## x2 1 643 643 5.31 0.033 *
## Residuals 18 2181 121
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(SSE = t(e) %*% e) # (5.84)
## [,1]
## [1,] 2181
(SSR = t(beta) %*% t(X) %*% Y - (1/n) * t(Y) %*% J %*% Y) # (5.85)
## y
## [1,] 24015
c(SSTO = SSTO, SSE = SSE, SSR = SSR) # In summary
## SSTO SSE SSR
## 26196 2181 24015
anova(lm(y ~ 1, df), fit) # All included here
## Analysis of Variance Table
##
## Model 1: y ~ 1
## Model 2: y ~ x1 + x2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 20 26196
## 2 18 2181 2 24015 99.1 1.9e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(SSTO = t(Y) %*% (I - (1/n) * J) %*% Y) # (5.89a)
## y
## y 26196
(SSE = t(Y) %*% (I - H) %*% Y) # (5.89b)
## y
## y 2181
(SSR = t(Y) %*% (H - (1/n) * J) %*% Y) # (5.89b)
## y
## y 24015
(MSE = as.vector(SSE / (n-3)))
## [1] 121.2
(se = MSE * solve(t(X) %*% X)) # (5.93)
## (Intercept) x1 x2
## (Intercept) 3602.035 8.74594 -241.4230
## x1 8.746 0.04485 -0.6724
## x2 -241.423 -0.67244 16.5158
vcov(fit) # Convenient!
## (Intercept) x1 x2
## (Intercept) 3602.035 8.74594 -241.4230
## x1 8.746 0.04485 -0.6724
## x2 -241.423 -0.67244 16.5158
sqrt(diag(se)) # Standard error for coefs.
## (Intercept) x1 x2
## 60.0170 0.2118 4.0640
sqrt(diag(vcov(fit))) # Much easier!
## (Intercept) x1 x2
## 60.0170 0.2118 4.0640
summary(fit) # See something familiar?
##
## Call:
## lm(formula = y ~ x1 + x2, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.424 -6.216 0.745 9.436 20.215
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -68.857 60.017 -1.15 0.266
## x1 1.455 0.212 6.87 0.000002 ***
## x2 9.366 4.064 2.30 0.033 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11 on 18 degrees of freedom
## Multiple R-squared: 0.917, Adjusted R-squared: 0.907
## F-statistic: 99.1 on 2 and 18 DF, p-value: 1.92e-10
(h = matrix(c(1, 34, 13))) # (5.95) x1 = 34, x2 = 13
## [,1]
## [1,] 1
## [2,] 34
## [3,] 13
(t(h) %*% beta) # (5.96)
## [,1]
## [1,] 102.3
(sh = MSE * (t(h) %*% solve(t(X) %*% X) %*% h)) # (5.98)
## [,1]
## [1,] 168.3
(sh = sqrt(sh))
## [,1]
## [1,] 12.97
# Compare with the below output for $fit and $se.fit
predict(fit, newdata = data.frame(x1 = 34, x2 = 13), se.fit = TRUE)
## $fit
## 1
## 102.3
##
## $se.fit
## [1] 12.97
##
## $df
## [1] 18
##
## $residual.scale
## [1] 11.01