First of all, is it heteroskedasticity or heteroscedasticity? According to
McCulloch (1985),
heteroskedasticity is the proper spelling, because when transliterating Greek words, scientists
use the Latin letter k in place of the Greek letter κ (kappa). κ sometimes is transliterated as
the Latin letter c, but only when these words entered the English language through French, such
as scepter.
Now that this is out of the way, we can get to the meat of this blogpost (foreshadowing pun).
A random variable is said to be heteroskedastic, if its variance is not constant. For example,
the variability of expenditures may increase with income. Richer families may spend a similar
amount on groceries as poorer people, but some rich families will sometimes buy expensive
items such as lobster. The variability of expenditures for rich families is thus quite large.
However, the expenditures on food of poorer families, who cannot afford lobster, will not vary much.
Heteroskedasticity can also appear when data is clustered; for example, variability of
expenditures on food may vary from city to city, but is quite constant within a city.
To illustrate this, let’s first load all the packages needed for this blog post:
library(robustbase)
library(tidyverse)
library(sandwich)
library(lmtest)
library(modelr)
library(broom)
First, let’s load and prepare the data:
data(&quot;education&quot;)
education % rename(residents = X1,
per_capita_income = X2,
young_residents = X3,
per_capita_exp = Y,
state = State) %__% mutate(region = case_when(
Region == 1 ~ &quot;northeast&quot;,
Region == 2 ~ &quot;northcenter&quot;,
Region == 3 ~ &quot;south&quot;,
Region == 4 ~ &quot;west&quot;
)) %__% select(-Region)
I will be using the education data sat from the {robustbase} package. I renamed some columns
and changed the values of the Region column. Now, let’s do a scatterplot of per capita expenditures
on per capita income:
ggplot(education, aes(per_capita_income, per_capita_exp)) + geom_point() +
theme_dark()
It would seem that, as income increases, variability of expenditures increases too. Let’s look
at the same plot by region:
ggplot(education, aes(per_capita_income, per_capita_exp)) + geom_point() + facet_wrap(~region) + theme_dark()
I don’t think this shows much; it would seem that observations might be clustered, but there are
not enough observations to draw any conclusion from this plot (in any case, drawing conclusions
from plots only is dangerous).
Let’s first run a good ol’ linear regression:
lmfit |t|) ## (Intercept) -467.40283 142.57669 -3.278 0.002073 ** ## regionnortheast 15.72741 18.16260 0.866 0.391338 ## regionsouth 7.08742 17.29950 0.410 0.684068 ## regionwest 34.32416 17.49460 1.962 0.056258 . ## residents -0.03456 0.05319 -0.650 0.519325 ## young_residents 1.30146 0.35717 3.644 0.000719 ***
## per_capita_income 0.07204 0.01305 5.520 1.82e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## ## Residual standard error: 39.88 on 43 degrees of freedom
## Multiple R-squared: 0.6292, Adjusted R-squared: 0.5774 ## F-statistic: 12.16 on 6 and 43 DF, p-value: 6.025e-08
Let’s test for heteroskedasticity using the Breusch-Pagan test that you can find in the {lmtest}
package:
bptest(lmfit)
## ## studentized Breusch-Pagan test
## ## data: lmfit
## BP = 17.921, df = 6, p-value = 0.006432
This test shows that we can reject the null that the variance of the residuals is constant,
thus heteroskedacity is present. To get the correct standard errors, we can use the vcovHC()
function from the {sandwich} package (hence the choice for the header picture of this post):
lmfit %__% vcovHC() %__% diag() %__% sqrt()
## (Intercept) regionnortheast regionsouth regionwest ## 311.31088691 25.30778221 23.56106307 24.12258706 ## residents young_residents per_capita_income ## 0.09184368 0.68829667 0.02999882
By default vcovHC() estimates a heteroskedasticity consistent (HC) variance covariance
matrix for the parameters. There are several ways to estimate such a HC matrix, and by default
vcovHC() estimates the “HC3” one. You can refer to Zeileis (2004)
for more details.
We see that the standard errors are much larger than before! The intercept and regionwest variables
are not statistically significant anymore.
You can achieve the same in one single step:
coeftest(lmfit, vcov = vcovHC(lmfit))
## ## t test of coefficients:
## ## Estimate Std. Error t value Pr(__|t|) ## (Intercept) -467.402827 311.310887 -1.5014 0.14056 ## regionnortheast 15.727405 25.307782 0.6214 0.53759 ## regionsouth 7.087424 23.561063 0.3008 0.76501 ## regionwest 34.324157 24.122587 1.4229 0.16198 ## residents -0.034558 0.091844 -0.3763 0.70857 ## young_residents 1.301458 0.688297 1.8908 0.06540 .
## per_capita_income 0.072036 0.029999 2.4013 0.02073 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
It’s is also easy to change the estimation method for the variance-covariance matrix:
coeftest(lmfit, vcov = vcovHC(lmfit, type = &quot;HC0&quot;))
## ## t test of coefficients:
## ## Estimate Std. Error t value Pr(__|t|) ## (Intercept) -467.402827 172.577569 -2.7084 0.009666 ** ## regionnortheast 15.727405 20.488148 0.7676 0.446899 ## regionsouth 7.087424 17.755889 0.3992 0.691752 ## regionwest 34.324157 19.308578 1.7777 0.082532 . ## residents -0.034558 0.054145 -0.6382 0.526703 ## young_residents 1.301458 0.387743 3.3565 0.001659 ** ## per_capita_income 0.072036 0.016638 4.3296 8.773e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
As I wrote above, by default, the type argument is equal to “HC3”.
Another way of dealing with heteroskedasticity is to use the lmrob() function from the
{robustbase} package. This package is quite interesting, and offers quite a lot of functions
for robust linear, and nonlinear, regression models. Running a robust linear regression
is just the same as with lm():
lmrobfit method = &quot;MM&quot;
## Residuals:
## Min 1Q Median 3Q Max ## -57.074 -14.803 -0.853 24.154 174.279 ## ## Coefficients:
## Estimate Std. Error t value Pr(__|t|) ## (Intercept) -156.37169 132.73828 -1.178 0.24526 ## regionnortheast 20.64576 26.45378 0.780 0.43940 ## regionsouth 10.79694 29.42747 0.367 0.71549 ## regionwest 45.22589 33.07950 1.367 0.17867 ## residents 0.03406 0.04412 0.772 0.44435 ## young_residents 0.57896 0.25512 2.269 0.02832 * ## per_capita_income 0.04328 0.01442 3.000 0.00447 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## ## Robust residual standard error: 26.4 ## Multiple R-squared: 0.6235, Adjusted R-squared: 0.571 ## Convergence in 24 IRWLS iterations
## ## Robustness weights: ## observation 50 is an outlier with |weight| = 0 ( __ 0.002); ## 7 weights are ~= 1. The remaining 42 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.05827 0.85200 0.93870 0.85250 0.98700 0.99790 ## Algorithmic parameters: ## tuning.chi bb tuning.psi refine.tol ## 1.548e+00 5.000e-01 4.685e+00 1.000e-07 ## rel.tol scale.tol solve.tol eps.outlier ## 1.000e-07 1.000e-10 1.000e-07 2.000e-03 ## eps.x warn.limit.reject warn.limit.meanrw ## 1.071e-08 5.000e-01 5.000e-01 ## nResample max.it best.r.s k.fast.s k.max ## 500 50 2 1 200 ## maxit.scale trace.lev mts compute.rd fast.s.large.n ## 200 0 1000 0 2000 ## psi subsampling cov ## &quot;bisquare&quot; &quot;nonsingular&quot; &quot;.vcov.avar1&quot; ## compute.outlier.stats ## &quot;SM&quot; ## seed : int(0)
This however, gives you different estimates than when fitting a linear regression model.
The estimates should be the same, only the standard errors should be different. This is because
the estimation method is different, and is also robust to outliers (at least that’s my understanding,
I haven’t read the theoretical papers behind the package yet).
Finally, it is also possible to bootstrap the standard errors. For this I will use the
bootstrap() function from the {modelr} package:
resamples % broom::tidy() %__% full_join(r.std.error) %__% select(term, estimate, std.error, r.std.error)
## Joining, by = &quot;term&quot;
## term estimate std.error r.std.error
## 1 (Intercept) -467.40282655 142.57668653 226.08735720
## 2 regionnortheast 15.72740519 18.16259519 19.65120904
## 3 regionsouth 7.08742365 17.29949951 19.05153934
## 4 regionwest 34.32415663 17.49459866 18.72767470
## 5 residents -0.03455770 0.05318811 0.06285984
## 6 young_residents 1.30145825 0.35716636 0.50928879
## 7 per_capita_income 0.07203552 0.01305066 0.02109277
As you see, using the whole bootstrapping procedure is longer than simply using either one of
the first two methods. However, this procedure is very flexible and can thus be adapted to a very
large range of situations. Either way, in the case of heteroskedasticity, you can see that
results vary a lot depending on the procedure you use, so I would advise to use them all as
robustness tests and discuss the differences.
If you found this blog post useful, you might want to follow me on twitter
for blog post updates.