# Heteroscedasticity in linear regression, there is a a pattern. What to do?

I'm modelling the behaviour of two variables with a linear regression. Since I saw (and believe) there is a multiplicative behaviour I transformed the dependent and independent variables taking the log.

fit <- lm(log_GDP ~ log_population, data=df)
par(mfrow = c(2, 2))
plot(fit)

From the second graph it seems that the residuals are not Normally distributed. Moreover, from the third graph we can see that there is no homogenity.

This is confirmed by:

library(gvlma )
gvmodel <- gvlma(fit)
summary(gvmodel)

Value   p-value                   Decision
Global Stat        63.543 5.216e-13 Assumptions NOT satisfied!
Skewness            8.617 3.330e-03 Assumptions NOT satisfied!
Kurtosis           49.484 2.000e-12 Assumptions NOT satisfied!
Link Function       3.163 7.534e-02    Assumptions acceptable.
Heteroscedasticity  2.279 1.311e-01    Assumptions acceptable.

Since I believe it is possible to remove some outliers to better fit this behaviour, I removed some of them

mod <- lm(log_GDP ~ log_population, data=df[-c(190,59, 228, 47, 221,48),])
gvlma <- gvlma(mod)
summary(gvlma)

Value p-value                Decision
Global Stat        5.36084 0.25223 Assumptions acceptable.
Skewness           2.86748 0.09039 Assumptions acceptable.
Kurtosis           0.82106 0.36487 Assumptions acceptable.
Link Function      1.65349 0.19848 Assumptions acceptable.
Heteroscedasticity 0.01882 0.89089 Assumptions acceptable.

Great! Let's see them:

From the third graph I can still see a pattern. Homoscedasticity!! Should I do something (and what?)? Are the assumptions valid?

If everything is all right, why should I specify Instrument Variables?

• "Homoscedasticity" (meaning a dispersion of residuals that is approximately constant regardless of the fitted values) is a correct interpretation of these plots and it's exactly what you want to see. What's the problem, then? – whuber Dec 19 '15 at 17:42
• @whuber I don't like so much the third graph (fitted values vs standardized residuals) because there is a decreasing pattern... I'm worried it causes problems – marcodena Dec 19 '15 at 17:50
• That is primarily because the number of data points decreases, too. As they decrease, the range of the data will tend to decrease. You need to account for that factor when reading such plots. I suspect you meant to write "Heteroscedasticity!" in your post. – whuber Dec 19 '15 at 19:46

First off, one should remember that homo/heteroskedasticity does NOT play a role in the important part of the model. As a model-builder you would like to know if the estimates are unbiased or consistent, and the assumptions you make about the variance of the residuals does NOT alter this.

The good thing about heteroskedasticityis is that is it an easy problem to fix. You have a couple of options, one is to compute heteroscedasticity consistent (HC) standard errors, in R the package “sandwich” can do this for you. Another choice is to use weighted least squares, I would not advice you to do so - since you would need to make more, not trivial, functional form assumptions (essentially build a model for the variance). Finally you can opt for the non-parametric way and bootstrap the standard errors.

Regardless of if opt for HC or bootstrapped standard errors the usual way of inference is completely valid, you simple use these new standard errors in place of the old ones.

I see that you use R, so here is an example:

set.seed(42)
x <- rnorm(10000)
u <- rnorm(10000)
y <- 50 + x + u
reg1 <- lm(y ~ x)

# usual standard errors:
U_se    <- sqrt(diag(vcov(reg1)))

# HC consistent standard errors:
require(sandwich)
HC_vcov <- vcovHC(reg1, type = "HC0")
HC_se   <- sqrt(diag(HC_vcov))

# Test the coefficents:
require(lmtest)
coeftest(reg1, vcov = HC_vcov)

# Compare:
require(texreg)
screenreg(l = list(reg1, reg1),
override.se = list(U_se, HC_se),
digits = 4)

===========================================
Model 1         Model 2
-------------------------------------------
(Intercept)     50.0007 ***     50.0007 ***
(0.0101)        (0.0101)
x                0.9927 ***      0.9927 ***
(0.0100)        (0.0101)
-------------------------------------------
R^2              0.4956          0.4956
As you can see the estimates are unchanged, also the $R^2$ and other statistics that does not depend on the estimates themselves are unchanged. In fact, since this model was generated under the CLM assumptions, there is not a whole lot of difference between the robust and usual OLS standard errors, generally one would expect that HC standard errors are larger. Finally note that it is always valid to use HC standard errors, even if the model is homoskedlastic - but in this case you loose efficency.