[R] Error in predict.lm() for models without intercept?

[isabella at ghement.ca]

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.  

Many thanks, 

Isabella 


Isabella R. Ghement, Ph.D. 
Ghement Statistical Consulting Company Ltd. 
301-7031 Blundell Road, Richmond, B.C., Canada, V6Y 1J5 
Tel: 604-767-1250 
Fax: 604-270-3922 
E-mail: isabella at ghement.ca 
Web: www.ghement.ca