[R] poly() in lm() leads to wrong coefficients (but correct residuals)

Markus Jäntti markus.jantti at iki.fi
Wed Jun 29 18:43:56 CEST 2005


On Wed, 2005-06-29 at 18:19 +0200, Andreas Neumann wrote:
> Dear all,
> 
> I am using poly() in lm() in the following form.
> 
> 1> DelsDPWOS.lm3 <- lm(DelsPDWOS[,1] ~ poly(DelsPDWOS[,4],3))
> 
> 2> DelsDPWOS.I.lm3 <- lm(DelsPDWOS[,1] ~ poly(I(DelsPDWOS[,4]),3))
> 
> 3> DelsDPWOS.2.lm3 <-
>    lm(DelsPDWOS[,1]~DelsPDWOS[,4]+I(DelsPDWOS[,4]^2)+I(DelsPDWOS[,4]^3))
> 
> 1 and 2 lead to identical but wrong results. 3 is correct. Surprisingly
> (to me) the residuals are the same (correct) in all cases although the
> coefficients of 1 (and 2) are wrong and do not in any way match the
> stated regression polynomial. (summaries below)
> 
> QUESTION:
> Is there a correct way to use poly() in lm() since version 3 becomes quite
> tedious if used more often especially with higher order polynomials?
> 

The coefficients  using 1 and 2 are not incorrect.  
poly() defines orthogonal polynomials, whereas your
DelsPDWOS[,4]+I(DelsPDWOS[,4]^2)+I(DelsPDWOS[,4]^3 
contruct defines an ordinary polynomial. 

You should be able to recover version 3 coefficients from 1 and 2.
See help(poly)

> x <- runif(10)
> x
 [1] 0.1878 0.2415 0.5834 0.6556 0.4112 0.3399 0.8144 0.1134 0.7360
0.0463
> model.matrix(~ poly(x, 2))
   (Intercept) poly(x, 2)1 poly(x, 2)2
1            1    -0.27648     -0.0452
2            1    -0.21052     -0.1899
3            1     0.20937     -0.2708
4            1     0.29799     -0.1021
5            1    -0.00212     -0.4117
6            1    -0.08970     -0.3621
7            1     0.49297      0.4968
8            1    -0.36790      0.2148
9            1     0.39672      0.1620
10           1    -0.45033      0.5082
attr(,"assign")
[1] 0 1 1
> model.matrix(~ x + I(x^2))
   (Intercept)      x  I(x^2)
1            1 0.1878 0.03528
2            1 0.2415 0.05834
3            1 0.5834 0.34040
4            1 0.6556 0.42982
5            1 0.4112 0.16911
6            1 0.3399 0.11554
7            1 0.8144 0.66320
8            1 0.1134 0.01286
9            1 0.7360 0.54169
10           1 0.0463 0.00214
attr(,"assign")
[1] 0 1 2
>


Regards,

-- 
Markus Jantti
Abo Akademi University
markus.jantti at iki.fi
http://www.iki.fi/~mjantti




More information about the R-help mailing list