[R] Orthogonal Polynomial Regression Parameter Estimation

WilD KID wildscop at yahoo.com
Thu May 6 13:09:13 CEST 2004


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?

--------------------------------------------

Data:

> Z<-1957:1964
> Y<-c(0.93,0.99,1.11,1.33,1.52,1.60,1.47,1.33)
> X<-Z-1956
> X
[1] 1 2 3 4 5 6 7 8

--------------------------------------------

Using lm function to get orthogonal polynomial, we
get-

> orth.fit<-lm(Y~poly(X,degree =6))
> coef(orth.fit)

(Intercept)            poly(X, degree = 6)1   poly(X,
degree = 6)2 
1.285000000                     0.529260490          
-0.316321867 

poly(X, degree = 6)3   poly(X, degree = 6)4   poly(X,
degree = 6)5 
        -0.221564684            0.054795962           
0.062910192 

poly(X, degree = 6)6 
         0.006154575 

--------------------------------------------

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. 

Thanks in advance for your time.

_______________________

Mohammad Ehsanul Karim <wildscop at yahoo.com> 
Institute of Statistical Research and Training 
University of Dhaka, Dhaka- 1000, Bangladesh




More information about the R-help mailing list