R devel - Jan 2002 - Bug in predict(newdata=x) with poly() (PR#1258)

Bug in predict.lm & poly

The predict function doesn't work when used with poly and newdata.

For example, I'd expect the following code to work, and plot a fitted
cubic to the nearly straight line:

 x <- 1:10
 y <- x + rnorm(10)/100
 plot(x,y)
 fit <- lm(y ~ poly(x,3))
 newx <- seq(1,10,len=100)
 lines(newx,predict(fit,newdata=data.frame(x=newx)))

However, the plotted line is way off the data.  The problem is that
poly(x,3) doesn't agree with poly(newx,3) at equal values of x, so
predict applies the fitted coefficients to the wrong values.

To fix this we need a predict.poly procedure analogous to the
predict.ns() and predict.bs() procedures. This doesn't look easy to
write.

Duncan Murdoch

--please do not edit the information below--

Version:
 platform = i386-pc-mingw32
 arch = x86
 os = Win32
 system = x86, Win32
 status = 
 major = 1
 minor = 4.0
 year = 2001
 month = 12
 day = 19
 language = R

Windows 98 4.10 (build 1998)  

Search Path:
 .GlobalEnv, package:ctest, Autoloads, package:base


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To:
r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

This is long-term known in S: there have been examples in V&R for years.
The S solution is predict.gam.

predict.poly is easy. The problem (as for bs and ns) is getting it used by
predict.lm.  Thomas Yee (I think) has made a solution available, but the
consensus was that it could be done more elegantly at a lower level.  We
already preserve information on contrasts and sets of levels: we need to
preserve similar information for poly() and spline terms.

I noticed this was still open last week, and put it on my (private) TODO
list.

On Sat, 12 Jan 2002 murdoch@stats.uwo.ca wrote:
> Bug in predict.lm & poly
>
> The predict function doesn't work when used with poly and newdata.
>
> For example, I'd expect the following code to work, and plot a fitted
> cubic to the nearly straight line:
>
>  x <- 1:10
>  y <- x + rnorm(10)/100
>  plot(x,y)
>  fit <- lm(y ~ poly(x,3))
>  newx <- seq(1,10,len=100)
>  lines(newx,predict(fit,newdata=data.frame(x=newx)))
>
> However, the plotted line is way off the data.  The problem is that
> poly(x,3) doesn't agree with poly(newx,3) at equal values of x, so
> predict applies the fitted coefficients to the wrong values.
>
> To fix this we need a predict.poly procedure analogous to the
> predict.ns() and predict.bs() procedures. This doesn't look easy to
> write.
>
> Duncan Murdoch
>
> --please do not edit the information below--
>
> Version:
>  platform = i386-pc-mingw32
>  arch = x86
>  os = Win32
>  system = x86, Win32
>  status >  major = 1
>  minor = 4.0
>  year = 2001
>  month = 12
>  day = 19
>  language = R
>
> Windows 98 4.10 (build 1998)
>
> Search Path:
>  .GlobalEnv, package:ctest, Autoloads, package:base
>
>
>
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> r-devel mailing list -- Read
http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
> Send "info", "help", or "[un]subscribe"
> (in the "body", not the subject !)  To:
r-devel-request@stat.math.ethz.ch
>
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>
-- 
Brian D. Ripley,                  ripley@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 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To:
r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

On Sat, 12 Jan 2002 07:49:01 +0000 (GMT), you wrote:
>predict.poly is easy. The problem (as for bs and ns) is getting it used by
>predict.lm.  Thomas Yee (I think) has made a solution available, but the
>consensus was that it could be done more elegantly at a lower level.  We
>already preserve information on contrasts and sets of levels: we need to
>preserve similar information for poly() and spline terms.
Yes, I didn't notice that predict.bs wasn't being called by
predict.lm.
>I noticed this was still open last week, and put it on my (private) TODO
>list.
A short-term workaround would be to put code like this into
predict.lm, at the point where tt contains the terms of the fit, and
it has been determined that newdata is not missing:

  if (any(c('poly','bs','ns') %in% flatten(tt)))
     warning("newdata not supported with 'poly', 'ns', or
'bs'")

I don't know if there's already a function like flatten in R; here's
a
simple one:

flatten <- function (expr) 
{
  if (is.expression(expr) | is.call(expr) | is.list(expr)) {
    unlist(lapply(as.list(expr), flatten))
  } else {
    as.character(expr)
  }
}

There may be other functions besides poly, bs and ns that return
bases; I don't know.

Duncan Murdoch
 


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To:
r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._