0

I am analyzing a dataset in which the variance of the error term in the regression is not constant for all observations. For this, I re-built the model, estimating heteroskedasticity-robust (Huber-White) standard errors using the coeftest function. Now, I want to use these new results for a prediction with predict() function.

The dataset looks like the following but with multiple X: 

set.seed(123)
x <- rep(c(10, 15, 20, 25), each = 25)
e <- c()
e[1:25] <- rnorm(25, sd = 10)
e[26:50] <- rnorm(25, sd = 15)
e[51:75] <- rnorm(25, sd = 20)
e[76:100] <- rnorm(25, sd = 25)
y <- 720 - 3.3 * x + e
model <- lm(y ~ x)
library(lmtest)
library(sandwich)
coeftest(model, vcov=vcovHC(model, "HC1"))

I found the following solution for the issue on the internet:

predict.rob <- function(x,vcov,newdata){
    if(missing(newdata)){ newdata <- x$model }
    tt <- terms(x)
    Terms <- delete.response(tt)
    m.mat <- model.matrix(Terms,data=newdata)
    m.coef <- x$coef
    fit <- as.vector(m.mat %*% x$coef)
    se.fit <- sqrt(diag(m.mat%*%vcov%*%t(m.mat)))
    return(list(fit=fit,se.fit=se.fit))}

The remaining problem is that my regression has more than 1 regressor.

Is there any way to addapt this resolution to multiple (7) explanatory variables?

Thanks in advance!

Celeste
  • 1
  • Welcome StackOverflow! In order to get a high quality answer, it may help to include details on what exactly you actually tried to do to solve this problem, and what ended up happening. For example, have you tried using your adaptation of the internet solution? What happened? – ctt May 08 '19 at 18:21

1 Answers1

0

I'm not sure but coeftest function is only performing a test. You can't use directly its result for your prediction. Maybe, you can in a way specifify to predict.lm the covaraince via vcovHC(model, "HC1"). I hope it will help you a bit.

Rémi Coulaud
  • 1,684
  • 1
  • 8
  • 19