# Recovering raw coefficients and variances from orthogonal polynomial regression

It seems that if I have a regression model such as $y_i \sim \beta_0 + \beta_1 x_i+\beta_2 x_i^2 +\beta_3 x_i^3$ I can either fit a raw polynomial and get unreliable results or fit an orthogonal polynomial and get coefficients that don't have a direct physical interpretation (e.g. I cannot use them to find the locations of the extrema on the original scale). Seems like I should be able to have the best of both worlds and be able to transform the fitted orthogonal coefficients and their variances back to the raw scale. I've taken a graduate course in applied linear regression (using Kutner, 5ed) and I looked through the polynomial regression chapter in Draper (3ed, referred to by Kutner) but found no discussion of how to do this. The help text for the poly() function in R does not. Nor have I found anything in my web searching, including here. Is reconstructing raw coefficients (and obtaining their variances) from coefficients fitted to an orthogonal polynomial...

1. impossible to do and I'm wasting my time.
2. maybe possible but not known how in the general case.
3. possible but not discussed because "who would want to?"
4. possible but not discussed because "it's obvious".

If the answer is 3 or 4, I would be very grateful if someone would have the patience to explain how to do this or point to a source that does so. If it's 1 or 2, I'd still be curious to know what the obstacle is. Thank you very much for reading this, and I apologize in advance if I'm overlooking something obvious.

• I don't understand your points. x, x$^2$ and x$^3$ are not orthogonal. Hence they are correlated and the regression parameters could be unstable, but it is not automatically the case that they are unreliable. Conversion to orthognonal polynomials may be more reliable. But what makes the coefficient of the original powers of x any more interpretable than the coefficients of the orthogonal polynomials? If x is the only variable as in the model y=a + bx then ∆y=yi-yi-1=b∆x and b is interpretable as the change in y per unit change in x. But with powers involved such interpretation is lost. – Michael R. Chernick Jul 7 '12 at 22:48
• I used a model with just x as the variable for simplicity, but in reality I'm comparing curves between treatment groups. So, depending on which terms are significant and their magnitude, I can interpret them-- for example an upward/downward overall shift, or a greater/lesser initial slope. Also, as my question says, a natural comparison to make between curves is the location of the maxima/minima, which is easier to interpret if it's on the original scale. So, your vote is for choice 3, I take it? – f1r3br4nd Jul 8 '12 at 4:44
• No I haven't figured out whether or not it is possible yet. I just did understand why you want to do it. – Michael R. Chernick Jul 8 '12 at 21:33
• Well, note that the model fit with orthogonal polynomials will have the exact same fit (i.e. the same $R^2$, same fitted values, etc.) as the model fit with the raw polynomial terms. So, if you're looking to relate this back to the original data, you can look at the coefficients for the raw terms but use the orthogonal polynomials to do inference for the individual terms in a way that "accounts for" the dependency between them. – Macro Jul 9 '12 at 12:27
• As it turns out, cubic splines and B-splines are in a class all by themselves, and are the best of two worlds. – Carl Nov 27 '16 at 1:35

Yes, it's possible.

Let $z_1, z_2, z_3$ be the non-constant parts of the orthogonal polynomials computed from the $x_i$. (Each is a column vector.) Regressing these against the $x_i$ must give a perfect fit. You can perform this with the software even when it does not document its procedures to compute orthogonal polynomials. The regression of $z_j$ yields coefficients $\gamma_{ij}$ for which

$$z_{ij} = \gamma_{j0} + x_i\gamma_{j1} + x_i^2\gamma_{j2} + x_i^3\gamma_{j3}.$$

The result is a $4\times 4$ matrix $\Gamma$ that, upon right multiplication, converts the design matrix $X=\pmatrix{1;&x;&x^2;&x^3}$ into $$Z=\pmatrix{1;&z_1;&z_2;&z_3} = X\Gamma.\tag{1}$$

After fitting the model

$$\mathbb{E}(Y) = Z\beta$$

and obtaining estimated coefficients $\hat\beta$ (a four-element column vector), you may substitute $(1)$ to obtain

$$\hat Y = Z\hat\beta = (X\Gamma)\hat\beta = X(\Gamma\hat\beta).$$

Therefore $\Gamma\hat\beta$ is the estimated coefficient vector for the model in terms of the original (raw, un-orthogonalized) powers of $x$.

The following R code illustrates these procedures and tests them with synthetic data.

n <- 10        # Number of observations
d <- 3         # Degree
#
# Synthesize a regressor, its powers, and orthogonal polynomials thereof.
#
x <- rnorm(n)
x.p <- outer(x, 0:d, ^); colnames(x.p) <- c("Intercept", paste0("x.", 1:d))
z <- poly(x, d)
#
# Compute the orthogonal polynomials in terms of the powers via OLS.
#
xform <- lm(cbind(1, z) ~ x.p-1)
gamma <- coef(xform)
#
# Verify the transformation: all components should be tiny, certainly
# infinitesimal compared to 1.
#
if (!all.equal(as.vector(1 + crossprod(x.p %*% gamma - cbind(1,z)) - 1),
rep(0, (d+1)^2)))
warning("Transformation is inaccurate.")
#
# Fit the model with orthogonal polynomials.
#
y <- x + rnorm(n)
fit <- lm(y ~ z)
#summary(fit)
#
# As a check, fit the model with raw powers.
#
fit.p <- lm(y ~ .-1, data.frame(x.p))
#summary(fit.p)
#
# Compare the results.
#
(rbind(Computed=as.vector(gamma %*% coef(fit)), Fit=coef(fit.p)))

if (!all.equal(as.vector(gamma %*% coef(fit)), as.vector(coef(fit.p))))
warning("Results were not the same.")
• Nice answer, especially $\Gamma$. In 8th line of code, do you have (1 + cross... - 1). Why we need +1 - 1.? – pauljohn32 Nov 29 '17 at 7:44
• @Paul Adding $1$ turns all components smaller than about $10^{-16}$ into true zeros due to floating point imprecision and subtracting $1$ cancels out that effect. This implements the test announced in the preceding code comment block. – whuber Nov 29 '17 at 15:50
• Two years later... @whuber, is it possible to expand this to the 95% CIs of the coefficients as well? – user2602640 Sep 21 '18 at 18:49
• @user2602640 Yes. You need to extract the variance-covariance matrix of the coefficients (use vcov in R) to convert variances computed in one basis to variances in the new basis, and then compute the CIs manually in the usual way. – whuber Sep 21 '18 at 19:00
• @whuber I followed your comment about halfway through, then lost you entirely... any chance you would take pity on a mathematically-challenged biologist and write it out in code? – user2602640 Sep 21 '18 at 19:03

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

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...

1. 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