Raw polynomials
To get ordinary polynomials as in the question use raw = TRUE
. Unfortunately there is an undesirable aspect with ordinary polynomials in regression. If we fit a quadratic, say, and then a cubic the lower order coefficients of the cubic are all different than for the quadratic, i.e. 56.900099702, -0.466189630, 0.001230536 for the quadratic vs. 6.068478e+01, -5.688501e-01, 2.079011e-03 after refitting with a cubic below.
library(ISLR)
fm2raw <- lm(mpg ~ poly(horsepower, 2, raw = TRUE), Auto)
cbind(coef(fm2raw))
## [,1]
## (Intercept) 56.900099702
## poly(horsepower, 2, raw = TRUE)1 -0.466189630
## poly(horsepower, 2, raw = TRUE)2 0.001230536
fm3raw <- lm(mpg ~ poly(horsepower, 3, raw = TRUE), Auto)
cbind(coef(fm3raw))
## [,1]
## (Intercept) 6.068478e+01
## poly(horsepower, 3, raw = TRUE)1 -5.688501e-01
## poly(horsepower, 3, raw = TRUE)2 2.079011e-03
## poly(horsepower, 3, raw = TRUE)3 -2.146626e-06
Orthogonal polynonials
What we would really like is to add the cubic term in such a way that the lower order coefficients that were fit using the quadratic stay the same after refitting with a cubic fit. To do this take linear combinations of the columns of poly(horsepower, 2, raw = TRUE)
and do the same with poly(horsepower, 3, raw = TRUE)
such that the columns in the quadratic fit are orthogonal to each other and similarly for the cubic fit. That is sufficient to guarantee that the lower order coefficients won't change when we add higher order coefficients. Note how the first three coefficients are now the same in the two sets below (whereas above they differ). That is, in both cases below the 3 lower order coefficients are 23.44592, -120.13774 and 44.08953 .
fm2 <- lm(mpg ~ poly(horsepower, 2), Auto)
cbind(coef(fm2))
## [,1]
## (Intercept) 23.44592
## poly(horsepower, 2)1 -120.13774
## poly(horsepower, 2)2 44.08953
fm3 <- lm(mpg ~ poly(horsepower, 3), Auto)
cbind(coef(fm3))
## [,1]
## (Intercept) 23.445918
## poly(horsepower, 3)1 -120.137744
## poly(horsepower, 3)2 44.089528
## poly(horsepower, 3)3 -3.948849
Same predictions
Importantly, since the columns of poly(horsepwer, 2)
are just linear combinations of the columnns of poly(horsepower, 2, raw = TRUE)
the two quadratic models (orthogonal and raw) represent the same models (i.e. they give the same predictions) and only differ in parameterization. For example, the fitted values are the same:
all.equal(fitted(fm2), fitted(fm2raw))
## [1] TRUE
This would also be true of the raw and orthogonal cubic models.
Orthogonality
We can verify that the polynomials do have orthogonal columns which are also orthogonal to the intercept:
nr <- nrow(Auto)
e <- rep(1, nr) / sqrt(nr) # constant vector of unit length
p <- cbind(e, poly(Auto$horsepower, 2))
zapsmall(crossprod(p))
## e 1 2
## e 1 0 0
## 1 0 1 0
## 2 0 0 1
Residual sum of squares
Another nice property of orthogonal polynomials is that due to the fact that poly
produces a matrix whose columns have unit length and are mutually orthogonal (and also orthogonal to the intercept column) the reduction in residual sum of squares due to the adding the cubic term is simply the square of the length of the projection of the response vector on the cubic column of the model matrix.
# these three give the same result
# 1. squared length of projection of y, i.e. Auto$mpg, on cubic term column
crossprod(model.matrix(fm3)[, 4], Auto$mpg)^2
## [,1]
## [1,] 15.5934
# 2. difference in sums of squares
deviance(fm2) - deviance(fm3)
## [1] 15.5934
# 3. difference in sums of squares from anova
anova(fm2, fm3)
##
## Analysis of Variance Table
##
## Model 1: mpg ~ poly(horsepower, 2)
## Model 2: mpg ~ poly(horsepower, 3)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 389 7442.0
## 2 388 7426.4 1 15.593 0.8147 0.3673 <-- note Sum of Sq value
poly(horsepower,2)
generates?poly(horsepower, degree=2, raw=TRUE)
; you're passing 2 as the wrong argument, andraw
defaults to FALSE.poly
generate the same output as the explicit formula, but I'd still like to know the actual form of the "orthogonal polynomial" thatpoly
generates without that parameter. Also, according to the manual, I am passing 2 as degree: " Although formally 'degree' should be named (as it follows '...'), an unnamed second argument of length 1 will be interpreted as the degree."...
, still it's good practice to name arguments after it.