Just a potentially useful additions to whuber's answer. Looking at the code for poly, you can deduce the linear map yourself. Let $\vec h_{m:n} =(h_m, h_{m + 1}, \dots, h_n)^\top$, let negative indices be zero by definition, and undefined $\Gamma$ entries be zero. Then we can find that if we disregard the scaling here then the map to the orthogonal polynomial is given by
$$\begin{align*}
z_0 &= 1 = \gamma_{0,0:0}\cdot 1\\
z_1 &= x - \alpha_1 \\
&= \underbrace{(\vec\gamma_{0,-1:0} - \alpha_1
\vec\gamma_{0,0:1})^\top}_{\vec\gamma_{1,0:1}^\top}(1 , x) \\
z_2 &= (x - \alpha_2)z_1 - \frac{\sigma_2}{\sigma_1} z_0 \\
&= x^2 + (\gamma_{10} -\alpha_2)x
- \alpha_2\gamma_{10} - \frac{\sigma_2}{\sigma_1} \\
&= (1, x, x^2)\underbrace{
(\vec\gamma_{1,-1:1}-\alpha_2\vec\gamma_{1,0:2}
-\frac{\sigma_2}{\sigma_1}\vec\gamma_{0,0:2})}_{
\vec\gamma_{2,0:2}^\top}\\
z_3 &= (x - \alpha_3)z_2 - \frac{\sigma_3}{\sigma_2} z_1 \\
&= x^3 + (\gamma_{21}-\alpha_3)x^2
+ (\gamma_{20}-\alpha_3\gamma_{21})x
-\alpha_3\gamma_{20}
-\frac{\sigma_3}{\sigma_2} x
- \frac{\sigma_3}{\sigma_2}\gamma_{01} \\
&= (1, x, x^2, x^3)\underbrace{(\vec\gamma_{2,-1:2}-\alpha_3\vec\gamma_{2,0:3}
-\frac{\sigma_3}{\sigma_2}\vec\gamma_{1,0:3})}_{
\vec\gamma_{3,0:3}}\\
\vdots\, &= \,\vdots
\end{align*}$$
Thus, we can compute the $\Gamma$ matrix with this code
get_poly_orth_map <- function(object){
stopifnot(inherits(object, "poly"))
sigs <- attr(object, "coefs")$norm2
alpha <- attr(object, "coefs")$alpha
nc <- length(alpha) + 1L
Gamma <- matrix(0., nc, nc)
Gamma[1, 1] <- 1
if(nc > 1){
Gamma[ , 2] <- -alpha[1] * Gamma[, 1]
Gamma[2, 2] <- 1
}
if(nc > 2)
for(i in 3:nc){
i_m1 <- i - 1L
Gamma[, i] <- c(0, Gamma[-nc, i_m1]) - alpha[i_m1] * Gamma[, i_m1] -
sigs[i] / sigs[i_m1] * Gamma[, i - 2L]
}
tmp <- sigs[-1]
tmp[1] <- 1
Gamma / rep(sqrt(tmp), each = nc)
}
and confirm that this gives the right matrix
# from whuber's answer
set.seed(1)
lm_method <- function(d, n = d * 4){
x <- rnorm(n, mean = 2)
x_p <- outer(x, 1:d, `^`)
colnames(x_p) <- paste0("x", 1:d)
poly_obj <- poly(x, d)
list(poly_obj = poly_obj, gamma = coef(lm(cbind(1, poly_obj) ~ x_p)))
}
# check that we get the same with different degrees
for(d in 1:10){
dat <- lm_method(d)
stopifnot(all.equal(
dat$gamma, get_poly_orth_map(dat$poly_obj), check.attributes = FALSE))
}
Is reconstructing raw coefficients (and obtaining their variances) from coefficients fitted to an orthogonal polynomial...
- impossible to do and I'm wasting my time.
As whuber shows, it is not. As another addition, here is an example to get the standard errors of the estimates as mentioned in the comments
# from `help(cars)`
fm <- lm(dist ~ speed + I(speed^2) + I(speed^3), data = cars)
summary(fm)
#R> Call:
#R> lm(formula = dist ~ speed + I(speed^2) + I(speed^3), data = cars)
#R>
#R> Residuals:
#R> Min 1Q Median 3Q Max
#R> -26.67 -9.60 -2.23 7.08 44.69
#R>
#R> Coefficients:
#R> Estimate Std. Error t value Pr(>|t|)
#R> (Intercept) -19.5050 28.4053 -0.69 0.50
#R> speed 6.8011 6.8011 1.00 0.32
#R> I(speed^2) -0.3497 0.4999 -0.70 0.49
#R> I(speed^3) 0.0103 0.0113 0.91 0.37
#R>
#R> Residual standard error: 15.2 on 46 degrees of freedom
#R> Multiple R-squared: 0.673, Adjusted R-squared: 0.652
#R> F-statistic: 31.6 on 3 and 46 DF, p-value: 3.07e-11
fp <- lm(dist ~ poly(speed, 3), data = cars)
gamma <- get_poly_orth_map(poly(cars$speed, 3))
drop(gamma %*% coef(fp))
#R> [1] -19.5050 6.8011 -0.3497 0.0103
sqrt(diag(tcrossprod(gamma %*% vcov(fp), gamma)))
#R> [1] 28.4053 6.8011 0.4999 0.0113
