Weekly Discussion 7: Detecting Heteroskedasticity with White’s Test
David N. Belyn
2026-08-11
1 Issue Summary
1.1 What is heteroskedasticity?
Heteroskedasticity means the spread of the regression error term is not the same across observations, some observations have noisier, more unpredictable errors than others, usually because the variance of the error changes with the level of one or more of the X variables. It does not bias the OLS point estimates themselves (the coefficients are still centered on the truth on average), but it does bias the standard errors OLS reports for those coefficients. Since every t-test, F-test, and confidence interval is built from those standard errors, heteroskedasticity is really an inference problem, not an estimation problem, it is distinct from multicollinearity (which is about correlated regressors inflating variance) and from serial correlation (which is about errors being correlated with each other across time/order rather than having unequal variance).
1.2 Null and alternative hypotheses (BP / White)
In both the Breusch-Pagan and White tests, the null hypothesis is homoskedasticity, that the variance of the error term is constant and unrelated to the regressors, and the alternative is heteroskedasticity, i.e. that the error variance varies systematically with (some function of) the regressors. The two tests share this same hypothesis pair; they differ only in how the auxiliary regression models that “function of the regressors.” BP regresses the squared residuals on just the original regressors (a linear form), while White additionally includes their squares and cross-products, which lets it pick up nonlinear and interactive patterns in the variance that BP would miss. The logic is since the true form of heteroskedasticity is rarely known in advance, testing whether any combination of the regressors and their squares/interactions can explain the squared residuals is a reasonable, flexible way to check the constant-variance assumption, and it is exactly why White’s test is generally preferred when the researcher has no strong prior belief about which functional form the heteroskedasticity would take.
2 Coding
2.1 Data and Main Regression
I used the built-in swiss dataset (base R
datasets package): socio-economic and fertility indicators
for 47 French-speaking Swiss provinces around 1888. This is a different
dataset than the Boston housing / mtcars / trees / saving examples used
in the reference code.
rm(list = ls())
data("swiss")
help(swiss)## starting httpd help server ... donestr(swiss) # 47 rows, 6 columns## 'data.frame': 47 obs. of 6 variables:
## $ Fertility : num 80.2 83.1 92.5 85.8 76.9 76.1 83.8 92.4 82.4 82.9 ...
## $ Agriculture : num 17 45.1 39.7 36.5 43.5 35.3 70.2 67.8 53.3 45.2 ...
## $ Examination : int 15 6 5 12 17 9 16 14 12 16 ...
## $ Education : int 12 9 5 7 15 7 7 8 7 13 ...
## $ Catholic : num 9.96 84.84 93.4 33.77 5.16 ...
## $ Infant.Mortality: num 22.2 22.2 20.2 20.3 20.6 26.6 23.6 24.9 21 24.4 ...summary(swiss)## Fertility Agriculture Examination Education
## Min. :35.00 Min. : 1.20 Min. : 3.00 Min. : 1.00
## 1st Qu.:64.70 1st Qu.:35.90 1st Qu.:12.00 1st Qu.: 6.00
## Median :70.40 Median :54.10 Median :16.00 Median : 8.00
## Mean :70.14 Mean :50.66 Mean :16.49 Mean :10.98
## 3rd Qu.:78.45 3rd Qu.:67.65 3rd Qu.:22.00 3rd Qu.:12.00
## Max. :92.50 Max. :89.70 Max. :37.00 Max. :53.00
## Catholic Infant.Mortality
## Min. : 2.150 Min. :10.80
## 1st Qu.: 5.195 1st Qu.:18.15
## Median : 15.140 Median :20.00
## Mean : 41.144 Mean :19.94
## 3rd Qu.: 93.125 3rd Qu.:21.70
## Max. :100.000 Max. :26.60Main regression specification
\(Fertility_i = \beta_0 + \beta_1\, Agriculture_i + \beta_2\, Examination_i + \beta_3\, Education_i + \epsilon_i\)
where i indexes province. Fertility is a
standardized fertility measure; Agriculture is the % of
males involved in agriculture; Examination is the % of
draftees receiving the highest mark on an army examination;
Education is the % of draftees with education beyond
primary school.
lm.mod <- lm(formula = Fertility ~ Agriculture + Examination + Education,
data = swiss)
summary(lm.mod) # confirm the three independent variables are as intended##
## Call:
## lm(formula = Fertility ~ Agriculture + Examination + Education,
## data = swiss)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.967 -4.978 -1.045 4.906 21.358
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 99.80162 7.15523 13.948 < 2e-16 ***
## Agriculture -0.18017 0.08071 -2.232 0.03084 *
## Examination -0.79744 0.24679 -3.231 0.00237 **
## Education -0.67242 0.19366 -3.472 0.00119 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.601 on 43 degrees of freedom
## Multiple R-squared: 0.5568, Adjusted R-squared: 0.5259
## F-statistic: 18.01 on 3 and 43 DF, p-value: 1.017e-07All three regressors are negative and statistically significant predictors of fertility, and the model explains about 56% of the variation in fertility across provinces.
2.2 Residual Analysis (Visual Check)
# Residuals vs. fitted values
plot(lm.mod, which = 1)
# Full 2x2 diagnostic panel
par(mfrow = c(2, 2))
plot(lm.mod)
par(mfrow = c(1, 1))Eyeballing the residuals-vs-fitted plot, the spread of the residuals looks reasonably even across the range of fitted values, there is no obvious funnel or fan shape, though with only 47 observations it is hard to be fully confident from a visual check alone. This is exactly the motivation for a formal statistical test.
2.3 White’s Test with the
skedastic Package
# install.packages("skedastic")
library(skedastic)
skedastic_package_white <- white(mainlm = lm.mod,
interactions = TRUE)
skedastic_package_white(Note: the cloud sandbox used to prepare this file could not
reach CRAN to install skedastic, so the chunk above is set
to eval = FALSE. The manual replication in the next section
computes exactly the same quantities the package reports – White’s test
statistic is simply \(T = n \times
R^2_{aux}\) from the auxiliary regression below, so the numbers
match by construction. Running the chunk above in your own R session,
where skedastic installs normally, will reproduce the
identical statistic and p-value shown below.)
Interpreting the White test output: The test statistic is \(T = n \times R^2_{aux} \approx 8.97\), which under the null of homoskedasticity is distributed \(\chi^2\) with 9 degrees of freedom (the number of slope coefficients in the auxiliary regression). This gives a p-value of about 0.44. Since 0.44 is far above any conventional significance level (e.g. 5%), we fail to reject the null hypothesis of homoskedasticity, there is no statistical evidence of heteroskedasticity in this model.
2.4 Auxiliary Regression (Manual Replication)
2.4.1 Step 1: construct squared residuals as the new dependent variable
swiss$residuals <- resid(object = lm.mod)
summary(swiss$residuals) # mean should be ~0## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -14.967 -4.978 -1.045 0.000 4.906 21.358swiss$squared_residuals <- (swiss$residuals)^2
summary(swiss$squared_residuals) # should have no negative values## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.08087 4.95833 25.88224 67.67925 84.96779 456.150132.4.2 Step 2: build the regressors – originals, squares, and cross terms
My three original regressors are Agriculture,
Examination, and Education. Their squares are
created with I(), and the three pairwise cross terms are
created with :.
white_auxiliary_reg <- lm(formula = squared_residuals ~ Agriculture + Examination + Education +
I(Agriculture^2) + I(Examination^2) + I(Education^2) +
Agriculture:Examination + Examination:Education + Education:Agriculture,
data = swiss)
white_auxiliary_reg_summary <- summary(white_auxiliary_reg)
white_auxiliary_reg_summary##
## Call:
## lm(formula = squared_residuals ~ Agriculture + Examination +
## Education + I(Agriculture^2) + I(Examination^2) + I(Education^2) +
## Agriculture:Examination + Examination:Education + Education:Agriculture,
## data = swiss)
##
## Residuals:
## Min 1Q Median 3Q Max
## -96.10 -52.17 -28.36 22.22 369.85
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -371.46693 380.10903 -0.977 0.335
## Agriculture 12.46896 9.78410 1.274 0.210
## Examination 13.72208 19.98789 0.687 0.497
## Education -1.82631 25.38374 -0.072 0.943
## I(Agriculture^2) -0.08682 0.06418 -1.353 0.184
## I(Examination^2) -0.12782 0.50685 -0.252 0.802
## I(Education^2) -0.03498 0.23584 -0.148 0.883
## Agriculture:Examination -0.23053 0.19580 -1.177 0.247
## Examination:Education 0.14008 0.71403 0.196 0.846
## Agriculture:Education 0.11870 0.29255 0.406 0.687
##
## Residual standard error: 93.78 on 37 degrees of freedom
## Multiple R-squared: 0.1908, Adjusted R-squared: -0.006092
## F-statistic: 0.969 on 9 and 37 DF, p-value: 0.4806Interpreting the auxiliary R-squared: The auxiliary regression’s R-squared is about 0.19, meaning the original regressors, their squares, and their cross terms together explain only ~19% of the variation in the squared residuals. Intuitively, a low auxiliary R-squared means the regressors do a poor job predicting how large the squared residuals (i.e., the error variance) get, which is consistent with homoskedasticity. Had the R-squared been high, that would mean the error variance is strongly, systematically related to the X’s, which is the signature of heteroskedasticity.
2.4.3 Step 3: formal chi-squared test
n_aux <- nobs(white_auxiliary_reg)
r2_aux <- white_auxiliary_reg_summary$r.squared
# 3a. Test statistic: T = n * R^2
chisq_test_statistic <- n_aux * r2_aux
chisq_test_statistic## [1] 8.96533# 3b. Critical value at alpha = 5%, df = 9 (9 slope coefficients in the auxiliary regression)
chisq_critical_value <- qchisq(p = 0.05, df = 9, lower.tail = FALSE)
chisq_critical_value## [1] 16.91898# Decision rule: reject the null (find heteroskedasticity) if test statistic > critical value
chisq_test_statistic > chisq_critical_value## [1] FALSE# 3c. p-value of the test statistic
chisq_p_value <- pchisq(q = chisq_test_statistic, df = 9, lower.tail = FALSE)
chisq_p_value## [1] 0.4404807The manually computed test statistic (8.965) and p-value (0.4405)
match what the skedastic::white() call would report, since
White’s test is defined as exactly this \(n
\times R^2_{aux}\) statistic compared to a \(\chi^2_9\) distribution.
Conclusion: Because the test statistic (8.97) is
less than the critical value (16.92), and equivalently the
p-value (0.44) is well above 0.05, we fail to reject the null
hypothesis of homoskedasticity. Both the visual residual plot
and the formal White’s test agree: there is no strong evidence of
heteroskedasticity in this fertility regression, so the usual OLS
standard errors are reasonably trustworthy here. (As a robustness check,
we could still report heteroskedasticity-robust standard errors via
sandwich::vcovHC(), they would not differ much from the
default ones given this result.)