# Explicit formulation of predict() output of an orthogonal polynomial regression?

In R, I am trying to reproduce the predict() output value of the orthogonal polynomial regression below. Based on my understanding of polynomial regression, I get 0.03869436 which is different from 0.05947406. Can anyone please help me by providing the explicit formulation of the predict output as a function of the fit model coefficients and p variable?

> q0 <- c(0.200,0.100,0.050,0.025)
> p0 <- c(0.325,0.409,0.477,0.534)
> p <- 0.4612118
> fit <- lm(q0 ~ poly(p0,3))
> predict(fit, newdata = list(p0 = p))
1
0.05947406
> # which is not the same as f(p) = (beta_0) + (beta_1)*p + (beta_2)*(p^2) + (beta_3)*(p^3) below
> as.numeric(fit$coefficients + fit$coefficients*p + fit$coefficients*(p^2) + fit$coefficients*(p^3))
 0.03869436


My internet searches have not turned up anything yet. Thank you.

• Have you looked at the R help page for poly()? – whuber Sep 20 '11 at 22:48

This does not handle orthogonal polynomials and implements ordinary (not orthogonal which IMHO are not worth the effort with modern stat computing algorithms) polynomials and restricted cubic splines (natural splines), but the R rms package is useful in a similar context.

require(rms)
f <- ols(y ~ x1 + rcs(x2,4) + pol(x3,3)) # spline for x2 with 4 knots; cubic for x3
Function(f)   # derive R function containing algebraic expression for predictions
latex(f)      # compose LaTeX markup for pretty algebraic form of fitted model


This looks to be an acceptable solution to my problem:

> fit <- lm(q0 ~ p0 + I(p0^2) + I(p0^3))
> as.numeric(fit$coefficients+fit$coefficients*(p) + fit$coefficients*(p^2) + fit$coefficients*(p^3))
 0.05947406


Thanks to all to considered the problem! (Also, thank you whuber for your suggestion.)

• I endorse the suggestion from @whuber to read the help page for poly(). Compare the results of coef(lm(q0 ~ poly(p0, 3, raw = FALSE))) (the default) and coef(lm(q0 ~ poly(p0, 3, raw = TRUE))) – user20637 Feb 25 '16 at 12:39