<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/dealing-with-heteroskedasticity-regression-with-robust-standard-errors-using-r/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.brodrigues.co/blog/2018-07-08-rob_stderr/"> Econometrics and Free Software</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<div>
<p><a rel="nofollow" href="https://cran.r-project.org/web/packages/sandwich/index.html" target="_blank"><br>
<img src="https://i1.wp.com/www.brodrigues.co/blog/2018-07-08-rob_stderr/img/bread-breakfast-bun-5678.jpg?w=450" data-recalc-dims="1"></a></p>
</div>
<p>First of all, is it heteros<strong>k</strong>edasticity or heteros<strong>c</strong>edasticity? According to<br>
<a rel="nofollow" href="https://www.jstor.org/stable/1911250" target="_blank">McCulloch (1985)</a>,<br>
heteros<strong>k</strong>edasticity is the proper spelling, because when transliterating Greek words, scientists<br>
use the Latin letter k in place of the Greek letter κ (kappa). κ sometimes is transliterated as<br>
the Latin letter c, but only when these words entered the English language through French, such<br>
as scepter.</p>
<p>Now that this is out of the way, we can get to the meat of this blogpost (foreshadowing pun).<br>
A random variable is said to be heteroskedastic, if its variance is not constant. For example,<br>
the variability of expenditures may increase with income. Richer families may spend a similar<br>
amount on groceries as poorer people, but some rich families will sometimes buy expensive<br>
items such as lobster. The variability of expenditures for rich families is thus quite large.<br>
However, the expenditures on food of poorer families, who cannot afford lobster, will not vary much.<br>
Heteroskedasticity can also appear when data is clustered; for example, variability of<br>
expenditures on food may vary from city to city, but is quite constant within a city.</p>
<p>To illustrate this, let’s first load all the packages needed for this blog post:</p>
<pre><code>library(robustbase)
library(tidyverse)
library(sandwich)
library(lmtest)
library(modelr)
library(broom)</code></pre>
<p>First, let’s load and prepare the data:</p>
<pre><code>data(&quot;education&quot;)

education &lt;- 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)</code></pre>
<p>I will be using the <code>education</code> data sat from the <code>{robustbase}</code> package. I renamed some columns<br>
and changed the values of the <code>Region</code> column. Now, let’s do a scatterplot of per capita expenditures<br>
on per capita income:</p>
<pre><code>ggplot(education, aes(per_capita_income, per_capita_exp)) + 
    geom_point() +
    theme_dark()</code></pre>
<p><img src="https://i2.wp.com/www.brodrigues.co/blog/2018-07-08-rob_stderr/blog/2018-07-08-rob_stderr_files/figure-html/unnamed-chunk-4-1.png?w=450" data-recalc-dims="1"></p>
<p>It would seem that, as income increases, variability of expenditures increases too. Let’s look<br>
at the same plot by <code>region</code>:</p>
<pre><code>ggplot(education, aes(per_capita_income, per_capita_exp)) + 
    geom_point() + 
    facet_wrap(~region) + 
    theme_dark()</code></pre>
<p><img src="https://i0.wp.com/www.brodrigues.co/blog/2018-07-08-rob_stderr/blog/2018-07-08-rob_stderr_files/figure-html/unnamed-chunk-5-1.png?w=450" data-recalc-dims="1"></p>
<p>I don’t think this shows much; it would seem that observations might be clustered, but there are<br>
not enough observations to draw any conclusion from this plot (in any case, drawing conclusions<br>
from plots only is dangerous).</p>
<p>Let’s first run a good ol’ linear regression:</p>
<pre><code>lmfit &lt;- lm(per_capita_exp ~ region + residents + young_residents + per_capita_income, data = education)

summary(lmfit)</code></pre>
<pre><code>## 
## Call:
## lm(formula = per_capita_exp ~ region + residents + young_residents + 
##     per_capita_income, data = education)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -77.963 -25.499  -2.214  17.618  89.106 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|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</code></pre>
<p>Let’s test for heteroskedasticity using the Breusch-Pagan test that you can find in the <code>{lmtest}</code><br>
package:</p>
<pre><code>bptest(lmfit)</code></pre>
<pre><code>## 
##  studentized Breusch-Pagan test
## 
## data:  lmfit
## BP = 17.921, df = 6, p-value = 0.006432</code></pre>
<p>This test shows that we can reject the null that the variance of the residuals is constant,<br>
thus heteroskedacity is present. To get the correct standard errors, we can use the <code>vcovHC()</code><br>
function from the <code>{sandwich}</code> package (hence the choice for the header picture of this post):</p>
<pre><code>lmfit %>% 
    vcovHC() %>% 
    diag() %>% 
    sqrt()</code></pre>
<pre><code>##       (Intercept)   regionnortheast       regionsouth        regionwest 
##      311.31088691       25.30778221       23.56106307       24.12258706 
##         residents   young_residents per_capita_income 
##        0.09184368        0.68829667        0.02999882</code></pre>
<p>By default <code>vcovHC()</code> estimates a heteroskedasticity consistent (HC) variance covariance<br>
matrix for the parameters. There are several ways to estimate such a HC matrix, and by default<br>
<code>vcovHC()</code> estimates the “HC3” one. You can refer to <a rel="nofollow" href="https://www.jstatsoft.org/article/view/v011i10" target="_blank">Zeileis (2004)</a><br>
for more details.</p>
<p>We see that the standard errors are much larger than before! The intercept and <code>regionwest</code> variables<br>
are not statistically significant anymore.</p>
<p>You can achieve the same in one single step:</p>
<pre><code>coeftest(lmfit, vcov = vcovHC(lmfit))</code></pre>
<pre><code>## 
## 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</code></pre>
<p>It’s is also easy to change the estimation method for the variance-covariance matrix:</p>
<pre><code>coeftest(lmfit, vcov = vcovHC(lmfit, type = &quot;HC0&quot;))</code></pre>
<pre><code>## 
## 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</code></pre>
<p>As I wrote above, by default, the <code>type</code> argument is equal to “HC3”.</p>
<p>Another way of dealing with heteroskedasticity is to use the <code>lmrob()</code> function from the<br>
<code>{robustbase}</code> package. This package is quite interesting, and offers quite a lot of functions<br>
for robust linear, and nonlinear, regression models. Running a robust linear regression<br>
is just the same as with <code>lm()</code>:</p>
<pre><code>lmrobfit &lt;- lmrob(per_capita_exp ~ region + residents + young_residents + per_capita_income, 
                  data = education)

summary(lmrobfit)</code></pre>
<pre><code>## 
## Call:
## lmrob(formula = per_capita_exp ~ region + residents + young_residents + per_capita_income, 
##     data = education)
##  \--> 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 ( &lt; 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)</code></pre>
<p>This however, gives you different estimates than when fitting a linear regression model.<br>
The estimates should be the same, only the standard errors should be different. This is because<br>
the estimation method is different, and is also robust to outliers (at least that’s my understanding,<br>
I haven’t read the theoretical papers behind the package yet).</p>
<p>Finally, it is also possible to bootstrap the standard errors. For this I will use the<br>
<code>bootstrap()</code> function from the <code>{modelr}</code> package:</p>
<pre><code>resamples &lt;- 100

boot_education &lt;- education %>% 
 modelr::bootstrap(resamples)</code></pre>
<p>Let’s take a look at the <code>boot_education</code> object:</p>
<pre><code>boot_education</code></pre>
<pre><code>## # A tibble: 100 x 2
##    strap          .id  
##    </code></pre>