Why are my prediction intervals a little too narrow?

Question

I'm trying to make a little function to make population-level predictions from a model that includes random effects, which in turn was fit using a package that doesn't support interval="prediction". (As far as I can tell, the reason why the package doesn't do so is that prediction is fraught when trying to predict within one of the model's levels -- but I just want population-level predictions).

So anyway, testing my code and comparing against the in-built code in predict.lm, my intervals are a little too narrow. Anybody able to help me figure out why?

set.seed(13)
n=100
x1 = runif(n)
x2 = runif(n) 
X = cbind(x1,x2)

y = x1 + x2 + rnorm(n)
m = lm(y~x1+x2)
summary(m)

d = predict(m,interval = "prediction")
auto = d[,3]-d[,1]

getpred = function(vals){
    qt(.975,n-(2+1))*sqrt(var(m$resid)* (1+ t(vals) %*% solve(t(X)%*%X) %*%vals))
    }
man = apply(X,1,getpred)

plot(density(auto),xlim=c(2.1,2.3),ylim = c(0,30))
lines(density(man),col="red")

Thanks in advance!

Answer 1

You forgot the intercept and did not calculate the residual variance correctly.

Here is how I would do it:

#get the design matrix
Designmat <- model.matrix(eval(eval(m$call)[-2]), as.data.frame(X))
#this is the same as cbind(1,X) for your model

#calculate the standard errors
predvar <- diag(Designmat %*% vcov(m) %*% t(Designmat))
SE <- sqrt(predvar) 
SE2 <- sqrt(predvar+summary(m)$sigma^2) 

tfrac <- qt(0.975, m$df.residual)

all.equal(auto,tfrac*SE2)
#TRUE

And here is your code fixed (though it is less general this way):

getpred <- function(vals){
  qt(.975,n-(2+1))*sqrt(summary(m)$sigma^2* (1+ t(vals) %*% solve(t(cbind(1,X))%*%cbind(1,X)) %*%vals))
}
man <- apply(cbind(1,X),1,getpred)
all.equal(unname(auto),man)
#TRUE

You can find code for calculation of prediction intervals for mixed-effects models in the r-sig-mixed models FAQ. I also put an example (with some important warnings) on StackOverflow.