[R] Orthogonal Polynomial Regression Parameter Estimation

[Douglas Bates]

WilD KID <wildscop at yahoo.com> writes:

> Dear all,
> 
> Can any one tell me how can i perform Orthogonal
> Polynomial Regression parameter estimation in R?
> 
> --------------------------------------------
> 
> Here is an "Orthogonal Polynomial" Regression problem
> collected from Draper, Smith(1981), page 269. Note
> that only value of alpha0 (intercept term) and signs
> of each estimate match with the result obtained from
> coef(orth.fit). What went wrong?

Draper and Smith (1981) are using unnormalized orthogonal polynomials
whereas the R function poly(X, degree = 6) produces normalized
columns.  That is, up to rounding error, the crossproduct of the
columns of the matrix produced by poly(X, degree = 6) is the identity.

The columns given in Draper and Smith are a set of orthogonal columns
but there are an infinite number of such sets and they each produce
different estimates of the coefficients.  The normalized orthogonal
polynomial columns are unique up to sign changes so the magnitudes of
the coefficients are well defined.

Generally the predictions from a linear model are well defined but the
coefficients are not necessarily well defined.  A better test would be
to see if the predictions from the poly model and the Draper and Smith
model were similar.  It would be best to check using all.equal to
allow for minor differences caused by numerical imprecision.

> Z<-1957:1964
> X<-Z-1956
> poly(X, degree = 6)
               1           2          3          4          5           6
[1,] -0.54006172  0.54006172 -0.4308202  0.2820380 -0.1497862  0.06154575
[2,] -0.38575837  0.07715167  0.3077287 -0.5237849  0.4921546 -0.30772873
[3,] -0.23145502 -0.23145502  0.4308202 -0.1208734 -0.3637664  0.55391171
[4,] -0.07715167 -0.38575837  0.1846372  0.3626203 -0.3209704 -0.30772873
[5,]  0.07715167 -0.38575837 -0.1846372  0.3626203  0.3209704 -0.30772873
[6,]  0.23145502 -0.23145502 -0.4308202 -0.1208734  0.3637664  0.55391171
[7,]  0.38575837  0.07715167 -0.3077287 -0.5237849 -0.4921546 -0.30772873
[8,]  0.54006172  0.54006172  0.4308202  0.2820380  0.1497862  0.06154575
attr(,"degree")
[1] 1 2 3 4 5 6
attr(,"coefs")
attr(,"coefs")$alpha
[1] 4.5 4.5 4.5 4.5 4.5 4.5

attr(,"coefs")$norm2
[1]    1.000    8.000   42.000  168.000  594.000 1810.286 4457.143 7854.545

attr(,"class")
[1] "poly"   "matrix"
> zapsmall(crossprod(poly(X, degree = 6)))
  1 2 3 4 5 6
1 1 0 0 0 0 0
2 0 1 0 0 0 0
3 0 0 1 0 0 0
4 0 0 0 1 0 0
5 0 0 0 0 1 0
6 0 0 0 0 0 1
> 
> And using the solution procedure given in Draper,
> Smith(1981) is -
> 
> --------------------------------------------
> 
> The following values are coefficients of 0-6th order
> (for n=8) polynomial collected from Pearson, Hartley
> (1958) table, page 212:
> 
> > p0<-rep(1,8)
> > p1<-c(-7,-5,-3,-1,1,3,5,7)
> > p2<-c(7,1,-3,-5,-5,-3,1,7)
> > p3<-c(-7,5,7,3,-3,-7,-5,7)
> > p4<-c(7,-13,-3,9,9,-3,-13,7)
> > p5<-c(-7,23,-17,-15,15,17,-23,7)
> > p6<-c(1,-5,9,-5,-5,9,-5,1)
> 
> Now, the estimated parameters of the orthogonal
> polynomial is calculated by the following formula:
> 
> > alpha0<-sum(Y*p0)/sum(p0^2);
> alpha1<-sum(Y*p1)/sum(p1^2);
> alpha2<-sum(Y*p2)/sum(p2^2);
> alpha3<-sum(Y*p3)/sum(p3^2);
> alpha4<-sum(Y*p4)/sum(p4^2);
> alpha5<-sum(Y*p5)/sum(p5^2);
> alpha6<-sum(Y*p6)/sum(p6^2)
> > alpha0;alpha1;alpha2;alpha3;alpha4;alpha5;alpha6
> [1] 1.285
> [1] 0.04083333
> [1] -0.02440476
> [1] -0.01363636
> [1] 0.002207792
> [1] 0.001346154
> [1] 0.0003787879
> 
> --------------------------------------------
> 
> Any response / help / comment / suggestion / idea /
> web-link / replies will be greatly appreciated.