Here is an example of the numerical instability that poly() can help with.
Imagine a 100-minute experiment where you chose to record the times as POSIXt
objects (stored as seconds since 1970). Using poly() gives the right predicted
values, but using I(x^2) does not:
> d <- data.frame(Date=as.POSIXct("2013-04-01")+60*(1:100),
> yQuadratic=(function(x)5*x-4*x^2)( (1:100)/50 ))
> d$NumDate <- as.numeric(d$Date)
> newD <- d[c(10, 80, 95),]
> cbind(newD,
> PQuadratic=predict(lm(yQuadratic~NumDate+I(NumDate^2),data=d),newdata=newD))
Date yQuadratic NumDate PQuadratic
10 2013-04-01 00:10:00 0.84 1364800200 2.1312
80 2013-04-01 01:20:00 -2.24 1364804400 -2.1808
95 2013-04-01 01:35:00 -4.94 1364805300 -3.1048
Warning message:
In predict.lm(lm(yQuadratic ~ NumDate + I(NumDate^2), data = d), :
prediction from a rank-deficient fit may be misleading
> cbind(newD,
> PQuadratic=predict(lm(yQuadratic~poly(NumDate,2),data=d),newdata=newD))
Date yQuadratic NumDate PQuadratic
10 2013-04-01 00:10:00 0.84 1364800200 0.84
80 2013-04-01 01:20:00 -2.24 1364804400 -2.24
95 2013-04-01 01:35:00 -4.94 1364805300 -4.94
Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com
> -----Original Message-----
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
> Behalf
> Of Paul Johnson
> Sent: Monday, April 01, 2013 10:21 AM
> To: R-help
> Subject: [R] example to demonstrate benefits of poly in regression?
>
> Here's my little discussion example for a quadratic regression:
>
> http://pj.freefaculty.org/R/WorkingExamples/regression-quadratic-1.R
>
> Students press me to know the benefits of poly() over the more obvious
> regression formulas.
>
> I think I understand the theory on why poly() should be more numerically
> stable, but I'm having trouble writing down an example that proves the
> benefit of this.
>
> I am beginning to suspect that because R's lm uses a QR decomposition under
> the hood, the instability that we used to see when we inverted (X'X) is no
> longer so serious as it used to be. When the columns in X are
>
> x1 x2 x3 x3squared
>
> then the QR decomposition is doing about the same thing as we would do if
> we replaced x3 and x3squared by the poly(x3, 2) results.
>
> Or not. I'm not arguing that, I'm thinking out loud, hoping you'll correct
> me.
>
> pj
>
>
> --
> Paul E. Johnson
> Professor, Political Science Assoc. Director
> 1541 Lilac Lane, Room 504 Center for Research Methods
> University of Kansas University of Kansas
> http://pj.freefaculty.org http://quant.ku.edu
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
