R help - May 2004 - Orthogonal Polynomial Regression Parameter Estimation

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

You have sent this to both the S and R lists under different names:  who
are you and which system are you interested in?

You use poly(), as you did, but the internals of poly() are not the same 
on the two systems.  Do read the documentation to find out which 
orthogonal polynomials are used.

Hint: The S-PLUS description says

   Returns a matrix of orthonormal polynomials, which represents a basis
   for polynomial regression.

Notice the difference?  Orthogonal polynomials are not uniquely defined
(nor are orthonormal ones, BTW).


On Thu, 6 May 2004, WilD KID wrote:
> 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
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
> 
> 
-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

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.

Dear Prof Brian Ripley,

	First of all, thanks for your kind reply. 
	Actually, my subscribed e-mail account in both the
lists are different (due to space considerations) -
that is why e-mail ID's were different, not the name -
"Mohammad Ehsanul Karim". Also, i am interested in any
system that solves my problem (i use R 1.7.1 and
S-plus 4.0 version).

	Anyway, can you give me any suggestions how should i
solve my problem (in any given system)? Yes - it is
possible to solve my problem on a step-by-step basis -
but i thought that some functions might be available
to solve this Orthogonal Polynomial Regression problem
directly. If lm(Y~poly(X,degree = 6)) does not work -
then is there any other way - or i am using it in an
improper way?

	Thank you for your time and effort.

_______________________

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

Dear Douglas Bates,

Thank you very much for your kind reply. 

_________________________

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

> Draper and Smith (1981) are using unnormalized
orthogonal 
> polynomials whereas the R function poly(X, degree 6) 
> produces normalized columns.

Dear all,

I came up with some sort of solution(!), which is :

--------------------------------------------
>  orth.fit<-lm(Y~poly(X,degree = 6))
>  co<-coef(orth.fit)
>  s<-c(sum(p0^2), sum(p1^2), sum(p2^2), sum(p3^2),
sum(p4^2), sum(p5^2), sum(p6^2)) # these pi's were
previously defined
>  nor<-co/sqrt(ss)
> 
co1<-c(alpha0,alpha1,alpha2,alpha3,alpha4,alpha5,alpha6)
# from previous results
>  cbind(co,ss,nor,co1)
                               co   ss           nor  
        co1
(Intercept)           1.285000000    8  0.4543161069 
1.2850000000
poly(X, degree = 6)1  0.529260490  168  0.0408333333 
0.0408333333
poly(X, degree = 6)2 -0.316321867  168 -0.0244047619
-0.0244047619
poly(X, degree = 6)3 -0.221564684  264 -0.0136363636
-0.0136363636
poly(X, degree = 6)4  0.054795962  616  0.0022077922 
0.0022077922
poly(X, degree = 6)5  0.062910192 2184  0.0013461538 
0.0013461538
poly(X, degree = 6)6  0.006154575  264  0.0003787879 
0.0003787879
--------------------------------------------

Note that all values of "nor" (modified R result) and
"co1" (Draper and Smith result) matches, except the
intercept term! But i guess i can live with it:p
However, if any one has any more suggestions, please
let me know.

Thanks to all who helped me here (Douglas Bates and
Prof Brian Ripley) and elsewhere (including Dale
McLerran of comp.soft-sys.sas and Nick Cox of Stata
list).

Thank you very much for all your support.

_______________________

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