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