Robust Prediction Intervals for LM

robust_predict.lm.R

predict.robust <- function(model, data, robust_vcov = NULL, level = 0.95, 
                           interval = "prediction"){
  # adapted from
  # https://stackoverflow.com/questions/38109501/how-does-predict-lm-compute-confidence-interval-and-prediction-interval
  # model is an lm object from r
  # data is the dataset to predict from
  # robust_vcov must be a robust vcov matrix created by V <- sandwich::vcovHC(model, ...)
  # level = the % of the confidence interval, default is 95%
  # interval = either "prediction" or "confidence" - prediction includes uncertainty about the model itself
  if(is.null(robust_vcov)){
    robust_vcov <- vcov(model)
  }
  if(is.null(data)){
     data <- model$model
   }
  fit <- as.numeric(model.matrix(model, data=data) %*% model$coefficients)
  se2 <- unname(rowSums((model.matrix(model) %*% robust_vcov) * model.matrix(model)))
  alpha <- qt((1-level)/2, df = model$df.residual)
  if(interval == "confidence"){
    upr <- fit + alpha*sqrt(se2)
    lwr <- fit - alpha*sqrt(se2)
  } else if(interval == "prediction"){
    sigma2 <- sum(model$residuals ^ 2) / model$df.residual
    upr <- fit + alpha*sqrt(se2+sigma2)
    lwr <- fit - alpha*sqrt(se2+sigma2)
  }
  preds <- cbind(fit, lwr, upr)
  return(list(fit = preds, se.fit = sqrt(se2)))
}