Using model with heteroskedasticity for predictions?

I am running a model with ordinary least squares regression, and am using robust standard errors (RSE's) because diagnostics tests indicated heteroskedasticity of the model. I'm a bit limited in the kinds of regression I can use because the data is non-normally distributed (I tried a few transformations but nothing really helped).

I will also be graphing the relationship between predicted values of y and values of different independent variables. I know the incorrect standard error issue affects interval estimates, but the coefficient values used in the regression equation remain the same after adjusting standard errors, so I am not sure if it would affect the predicted y values as well. My question is, since the model itself is still heteroskedastistic, would it still be "valid" to use the model for predictions? If not, is there a way I could adjust the actual model in R so that I could use it with the predict() function?

If it helps to have an idea of what I'm working with, here is the model without adjusted standard errors:

Call:
lm(formula = ortho ~ forb + sfdist + year, data = insect.ortho)

Residuals:
Min      1Q  Median      3Q     Max
-5.2720 -1.5416 -0.6649  0.7500 18.1564

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  3.140313   0.367324   8.549 3.76e-15 ***
forb         0.218135   0.051850   4.207 3.96e-05 ***
sfdist      -0.001098   0.000373  -2.943  0.00365 **
year        -0.762830   0.401225  -1.901  0.05876 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.603 on 193 degrees of freedom Multiple R-squared: 0.149, Adjusted R-squared: 0.1358 F-statistic: 11.27 on 3 and 193 DF, p-value: 7.621e-07

And after using robust standard errors:

t test of coefficients:

Estimate  Std. Error t value  Pr(>|t|)
(Intercept)  3.14031323  0.37927174  8.2799 2.017e-14 ***
forb         0.21813498  0.06895898  3.1633  0.001813 **
sfdist      -0.00109765  0.00033704 -3.2567  0.001331 **
year        -0.76282962  0.39007561 -1.9556  0.051956 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
• Personally, I would inspect the heteroskedasticity and my assumptions regarding the model error term a bit more. What is the reason for it? Should you possibly use a GLM or transform the DV? I might even try to fit a variance structure. – Roland Sep 19 '17 at 6:40

The prediction will not be altered in any way by using het-robust standard errors. It remains the same and is still valid.

It's the interval around that prediction (and any hypothesis tests about coefficients or the predictions) that will changed by the choice of of whether to use the het-robust errors or not. In general, if you have heteroskedasticity and use the non-het-robust errors, your intervals will be too small.

Is it possible you haven't chosen the correct transformation when you were testing them? I would suggest performing a Box-Cox test test on the data to see what the algorithm recommends (if you haven't already). The 'caret' package has a useful function called preProcess() which let's you run the Box-Cox, then fit it to your data using predict(). On a side note, when you perform a Breusch–Pagan test, what p-value are you getting (is it borderline?).

• Hi- I've tried Box-Cox transformation on the model, and log and square root transformations on the response variable. With Box-Cox, I get the message: "Lambda could not be estimated; no transformation is applied". For square root I get a TINY p-value when I perform Shapiro-Wilk on the response variable, and for log I get something not much better Interestingly, I get a "normal" value on the B-P test (0.3367) but a very small p-value with the NCV test (0.0015), and a graph of residuals vs. fitted values indicates heteroskedasticity so that's what I'm going with. – zugunru Sep 21 '17 at 15:18

I would not use the model to make predictions because one of the assumptions of regression models is homoscedasticity. This problem might be better solved using something like a GLM.

• The assumption of homoskedasticity can usually be relaxed with regression quite easily. – Dimitriy V. Masterov Sep 18 '17 at 22:50
• (-1) This is not correct. Linear regression is still a consistent estimator of the conditional mean even when errors are heteroskedastic. Also, the phrase "This problem might be better solved using gradient descent instead of regression" makes no sense to me. That's an apples an oranges comparison. Linear regression is a statistical/ML model, gradient descent is an optimization procedure. – Matthew Drury Sep 19 '17 at 0:07
• Thank you for the replies! @MatthewDrury are you saying you think the model would be ok to use for predictions? Dimitriy Masterov, I'm not sure if you are saying I don't need to assume homoskedasticity (in which case I'm not sure what type of regression you mean) or that I can adjust the model itself, which is what I'm asking for technical advice on. I am posing this question because I've found plenty on how to find the robust standard errors, but nothing on how I could actually apply them to the model in R. Thank you – zugunru Sep 19 '17 at 17:20
• Yup, homoskedacity does not infringe on the predictive capabilities of linear regression. There may be other, better models to use, due to either better capturing non-linear effects (ala logistic regression), or using the data more efficiently, but the concerns that lead to the homoskedactiy "assumption" are non-predictive in nature. – Matthew Drury Sep 19 '17 at 18:31
• Thanks for the comments. I saw the t-test and assumed statistical inference was in order, not model prediction (should have read all the text!). Yes, Matthew is right--heteroscedasticity will not influence a model's predictive accuracy; homoscedasticity is only necessary for consistency and unbiasedness. – Jesse Lawson Sep 20 '17 at 3:20