When I run your code (in the single predictor case) I get exactly what I
would expect. In particular the standard errors are indeed proportional
to the (absolute) value of x, and the standard error is indeed 0 at x = 0.
The proportionality constant is exactly what it should be, explicitly
s/sqrt(sum(x^2)) where "s" is the estimate of sigma (obtained via
summary(mod)$sigma).
So all is in harmony in the universe.
It is ***VERY*** unlikely that there is any error in predict.lm(), which
has been around and in heavy use for a long, long time.
I don't know what you are seeing to make you think that there is an
error, but I am not seeing anything untoward.
cheers,
Rolf Turner
On 17/09/14 10:51, isabe...@ghement.ca wrote:
Hi everyone,
It appears my R code didn't come through the first time (thanks for letting me
know, Ista). Here is my message again:
Could there be an error in the predict() function in R for models without
intercepts which include one or more predictors?
When using the predict() function to obtain the standard errors of the fitted
values produced by a linear model (lm), the
behaviour of the standard errors seems to mirror that of models which include
an intercept (which should not happen).
Here is an attempt to produce standard errors (and hence confidence bands) for
the fitted values in a linear model with a
single predictor and no intercept using R code:
## simulate a response y and two predictors x and z
x <- rnorm(100,mean=0, sd=1)
z <- runif(100,min=-1,max=1)
y <- 1*x + 2*z + rnorm(100,mean=0, sd=1)
## fit a linear model with no intercept but with one predictor
mod <- lm(y ~ 0 + x)
## compute confidence bands (i.e., fitted values +/- 1.96 standard errors of
fitted values)
conf.band.x <- predict(mod,newdata=data.frame(x =
seq(from=ceiling(min(x)),to=floor(max(x)),by=0.01)),
interval="confidence")
## display confidence bands
conf.band.x <- data.frame(lwr=conf.band.x[,"lwr"],
fit=conf.band.x[,"fit"],
upr=conf.band.x[,"upr"])
matplot(x=seq(from=ceiling(min(x)),to=floor(max(x)),by=0.01), y=conf.band.x, type="l",
xlab="x", ylab="y")
abline(v=mean(x),lty=3,col="magenta")
title("Effect of x on y")
According to statistical theory, in a model with no intercept and one
predictor, the standard errors should be directly
proportional to the value of x at which they are evaluated. If x=0, the
standard errors should also be zero. If x increases,
the standard errors should also increase. The resulting plot produced by
matplot shows this is not the case - the standard
errors appear to increase as one moves away from the average value of x. We
would expect this behaviour if the model included
an intercept, which is not the case here.
Here is some R code for looking at standard errors of fitted values when the
model includes no intercept and two predictors x
and z. In this code, the value of the predictor z is set to its average level.
## linear model with no intercept but with two predictors
mod <- lm(y ~ 0 + x + z)
conf.band.x <- predict(mod,newdata=data.frame(x =
seq(from=ceiling(min(x)),to=floor(max(x)),by=0.01),
z = mean(z)),
interval="confidence")
conf.band.x <- data.frame(lwr=conf.band.x[,"lwr"],
fit=conf.band.x[,"fit"],
upr=conf.band.x[,"upr"])
matplot(x=seq(from=ceiling(min(x)),to=floor(max(x)),by=0.01), y=conf.band.x, type="l",
xlab="x", ylab="y")
abline(v=mean(x),lty=3,col="magenta")
title("Partial Effect of x on y (obtained by setting z to its average level)")
Again, the standard errors seem to behave as though they would come from a
model with an intercept (given that they are
flaring up as one moves away from the average value of the predictor x).
I would very much appreciate any clarifications or suggestions for how to fix
this problem.
If the problem is confirmed, it appears to also carry over to the effects
package in R, which constructs plots similar to the
ones produced by matplot above by relying on the predict() function.
--
Rolf Turner
Technical Editor ANZJS
______________________________________________
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.
