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

isabella at ghement.ca isabella at ghement.ca
Wed Sep 17 00:51:39 CEST 2014

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