predict() still not yet "safe".. (PR#1840)

maechler@stat.math.ethz.ch maechler@stat.math.ethz.ch
Fri, 26 Jul 2002 20:12:16 +0200 (MET DST)


Executive Summary: This only concerns  predict() at new points
	  when the formula contained things like poly() or ns() {in the RHS}.

Because I hadn't have time to fix this yet, I submit it.
[R version 1.5.1 and at least 1.4.1, probably all versions at all:]
For these, R 1.5.x does not differ from 1.4.1 : Inspite of R
1.5.0's "NEWS" entry and the ?Safepredict help page, 
predict() doesn't (always?) work correctly for (some cases of)
   predict(lm(y ~  poly(....)), newdata = .....)
or
   predict(lm(y ~  ns(....)),   newdata = .....)

Here is a simple (longer than necessary) example 

#### This modified from the last part of  example(cars) {of R 1.5.1}:

data(cars)
## An example of polynomial regression
## -- to show the bug, the linear example is sufficient
cars.1 <- lm(dist ~ poly(speed, degree = 1), data = cars)
cars1  <- lm(dist ~      speed,              data = cars)# the `same'
all.equal(predict(cars.1), predict(cars1))# the fitted values are ok

cars.4 <- lm(dist ~ poly(speed, degree = 4), data = cars)
newd <- data.frame(speed = d <- seq(0, 25, length = 51)) # shorter than in ?cars
pr1d <- predict(cars.1, newd)
pr4d <- predict(cars.4, newd)

plot(cars, xlab = "Speed (mph)", ylab = "Stopping distance (ft)",
     las = 1, xlim = c(0, 25))
abline(cars1, col = "blue")# the correct one
## in R 1.5.0 and 1.5.1, these look wrong (particularly to the left!):
lines(d, pr1d)
lines(d, pr4d, col = "red")

## or see by simple prediction at one point
predict(cars1,  data.frame(speed = 4)) # -1.84946
predict(cars.1, data.frame(speed = 4))
## --> Error in poly(speed, degree = 1) :
##          degree must be less than number of points


-- 
Martin Maechler <maechler@stat.math.ethz.ch>	http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum  LEO C16	Leonhardstr. 27
ETH (Federal Inst. Technology)	8092 Zurich	SWITZERLAND
phone: x-41-1-632-3408		fax: ...-1228			<><

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._