Fitting a linear model, lm(y ~ X), by least squares can be thought of geometrically as the orthogonal projection of
y on the column space of X. This function is designed to allow exploration of projections
and orthogonality.
Proj(y, X, list = FALSE)a vector, treated as a one-column matrix
a vector or matrix. Number of rows of y and X must match
logical; if FALSE, return just the projected vector; otherwise returns a list
the projection of y on X (if list=FALSE) or a list with elements y and P
The projection is defined as \(P y\) where \(P = X (X'X)^- X'\) and \(X^-\) is a generalized inverse.
Other vector diagrams:
arc(),
arrows3d(),
circle3d(),
corner(),
plot.regvec3d(),
pointOnLine(),
regvec3d(),
vectors3d(),
vectors()
X <- matrix( c(1, 1, 1, 1, 1, -1, 1, -1), 4,2, byrow=TRUE)
y <- 1:4
Proj(y, X[,1]) # project y on unit vector
#> [1] 2.5 2.5 2.5 2.5
Proj(y, X[,2])
#> [1] -1 -1 1 1
Proj(y, X)
#> [1] 1.5 1.5 3.5 3.5
# project unit vector on line between two points
y <- c(1,1)
p1 <- c(0,0)
p2 <- c(1,0)
Proj(y, cbind(p1, p2))
#> [1] 1 0
# orthogonal complements
y <- 1:4
yp <-Proj(y, X, list=TRUE)
yp$y
#> [1] 1.5 1.5 3.5 3.5
P <- yp$P
IP <- diag(4) - P
yc <- c(IP %*% y)
crossprod(yp$y, yc)
#> [,1]
#> [1,] 0
# P is idempotent: P P = P
P %*% P
#> [,1] [,2] [,3] [,4]
#> [1,] 0.5 0.5 0.0 0.0
#> [2,] 0.5 0.5 0.0 0.0
#> [3,] 0.0 0.0 0.5 0.5
#> [4,] 0.0 0.0 0.5 0.5
all.equal(P, P %*% P)
#> [1] TRUE