Linear regression with violation of heteroskedasticity with small samples

TLDR: In small samples, the wild bootstrap implemented in the R hcci package is a good bet when heteroskedasticity is a concern.

Today while teaching the multiple regression lab, I showed the class the standardized residuals versus standardized predictor plot SPSS lets you produce. It is the plot we typically use to assess homoskedasticity. The sample size for the analysis was 44. I mentioned how the regression slopes are fine under heteroskedasticity, but inference \((t, SE, p\)-value) may be problematic. Supplemental R material I created included how to use the sandwich package to obtain heteroskedasticity-consistent standard errors (HCSEs). And after applying HC3 (or any of the HCs from 1 to 5), a regression coefficient was no longer statistically significant at \(\alpha=.05\)

I mentioned to the class that some folks would recommend applying HCSEs by default. After class, I tried to learn about the difference between the different HCs. The following papers were helpful: Zeileis (2004),¹ Long & Ervin (2000),² Cribari-Neto, Souza & Vasconcellos (2007),³ and Hausman & Palmer (2012)⁴. The documentation for the sandwich package was a big help. The Hausman & Palmer (H&P) paper is probably best if you’re only going to read one of the papers, and it can also serve as a short handy reference for dealing with heteroskedasticity at small sample sizes.

I learned that HCSEs can be problematic (H&P Table 1). Additionally, the Wild Bootstrap does a good job of maintaining the nominal error rate in small samples (n=40) under homoskedasticity, moderate heteroskedasticity and severe heteroskedasticity (H&P Table 1). It is also statistically powerful (H&P Fig. 1 & 2). The good thing is the hcci package contains a function called Pboot() which performs the wild bootstrap to correct for heteroskedasticity.

As far as I see, the function has one limitation: when you perform your regression, you cannot use the optional dataframe argument in lm(). Here’s an example with this dataset:

library(hcci)
atlschools <- read.csv("./atlschools.csv")

You can not pass the dataframe to the Pboot function so the next few lines are required prior to calling lm().

ppc <- atlschools$PPC # per-pupil costs
ptr_c <- scale(atlschools$PTR, scale = FALSE) # pupil/teacher ratio
mts_c_10 <- scale(atlschools$MTS, scale = FALSE) / 10 # monthly teacher salary

coef(summary(fit.0 <- lm(ppc ~ ptr_c + mts_c_10)))
             Estimate Std. Error   t value     Pr(>|t|)
(Intercept) 67.884318  1.1526357 58.894861 3.017231e-41
ptr_c       -2.798285  0.3685282 -7.593138 2.427617e-09
mts_c_10     2.477010  0.8167532  3.032752 4.190607e-03

Pboot(model = fit.0, J = 1000, K = 100)

$beta
[1] 67.884318 -2.798285  2.477010

$ci_lower_simple
[1] 65.5454924 -3.7301276 -0.0653991

$ci_upper_simple
[1] 70.221038 -1.904783  4.969260

The CI of monthly teacher salary includes 0, evidence to suggest we cannot distinguish its slope from 0. The inference at \(\alpha=.05\) is different from OLS.