R help - Sep 2014 - Error in predict.lm() for models without intercept?

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

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, isabella at 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