Wild bootstrap heteroscedasticity test¶

Source: R/modern_diagnostics.R · Help: ?performWildBootstrapTest

Calibrates the Breusch-Pagan Lagrange multiplier statistic with a null-imposed wild bootstrap, providing a p-value that is accurate in small samples and under non-normal errors without relying on the asymptotic \(\chi^2\) approximation.

Usage¶

performWildBootstrapTest(
  model,
  data = NULL,
  B = 999,
  distribution = c("rademacher", "mammen"),
  progress = interactive(),
  reestimate_weights = NULL
)

Arguments¶

model: A fitted stats::lm object describing the mean structure whose residual variance is to be assessed.
data: A base::data.frame (or compatible object) containing the variables referenced in model. The data must include all observations used to fit model and should not contain unresolved missing values. When NULL, the default, the model frame of model is used. It holds the terms of the model rather than its variables (log(x) rather than x), so a model that transforms a variable needs data, as does a model with an offset() term, and so may a fit whose weights or offset argument is a column of the data, which the model frame holds as (weights) or (offset).
B: Integer number of bootstrap replications. Defaults to 999, which is the default of every resampling test in the package since 0.15.0; it was 499 before. The p-value is (1 + #) / (B_eff + 1), so a test at level alpha is exact only when alpha * (B + 1) is an integer, which 999 is at the 5% and 1% levels.
distribution: Multiplier distribution used for the wild perturbation. Supported options are "rademacher" and "mammen".
progress: Logical flag controlling whether progress should be reported when running the bootstrap loop.
reestimate_weights: Whether each replicate estimates the weights of a fitWLS() fit again from the variance function fitWLS() recorded on it. NULL, the default, does so for a fitWLS() fit and holds the weights fixed for any other model. FALSE holds them fixed, as every call did before 0.15.0. TRUE is an error for a model that fitWLS() did not return, or on which it fell back to equal weights. Each replicate regenerates the response and holds the design fixed, so a variance regressor that uses a variable of the response which the right-hand side of the model does not use, such as y or x for I(y / x) ~ z, or that is the response term itself, as ~ . gives, cannot be estimated again: NULL then holds the weights fixed, with a warning, and TRUE is an error. Giving the response as a column of the data, as in r ~ z, says which variables stay. See Weighted fits.

Value¶

An object of class htest containing the observed Breusch-Pagan statistic, bootstrap p-value, and additional details under the bootstrap element.

Details¶

The observed statistic is the usual \(n R^2\) from the auxiliary regression of the squared residuals on the regressors. The bootstrap reference distribution, however, must be generated under the homoscedastic null — otherwise it inherits the heteroscedasticity the test is looking for and the procedure loses all power. Each bootstrap sample is therefore built from leverage-standardised, mean-centred residuals \(r_i = e_i/\sqrt{1-h_{ii}}\) (which have a common variance under \(H_0\)) resampled i.i.d. and perturbed by a wild multiplier \(v_i\), so that \(y_i^* = \hat y_i + r^*_i v_i\). The statistic is recomputed on each refit and the p-value is \((1 + \#\{T_b^* \ge T\})/(B + 1)\). In the package's validation study of the multiplier, with B = 999, the test held its level under \(H_0\) with heavy-tailed \(t_5\) errors but not with skewed ones: with centred exponential errors and 30, 50 and 100 observations it rejected a true null hypothesis 8.2%, 7.3% and 6.1% of the time at the 5% level with the Rademacher multiplier, and 7.8%, 7.0% and 5.9% with the Mammen multiplier. Under the alternative the observed statistic is extreme relative to the homoscedastic reference, giving power.

Each replicate is refitted on the design matrix of the fit, with its offset, whether that was given as an argument or written in the formula. The formula is not evaluated again, so a transformed regressor such as log(x), or a transformed response, is refitted on the scale of the fit.

Weighted fits¶

A fit with weights \(w_i\) is ordinary least squares on \(\sqrt{w_i}\, y_i\) and \(\sqrt{w_i}\, x_i\), whose errors have constant variance under the null hypothesis, and the bootstrap works in that regression. The statistic is computed from its residuals, the Pearson residuals \(\sqrt{w_i}\, e_i\), with the auxiliary regressors on their original scale; the residuals that are resampled are the Pearson residuals divided by \(\sqrt{1 - h_{ii}}\), with the leverages of that regression, and centred; and each replicate refits it. With equal weights the result is that of the unweighted fit.

Weights given to the fit are known, and they are held fixed in every replicate. On a fit returned by fitWLS() they were estimated from the same data, and since 0.15.0 they are estimated again in each replicate by default (reestimate_weights). The replicate is fitted as the model fitWLS() was given: by least squares, with the weights that model had, if any. The variance function fitWLS() recorded, with its form, is then fitted again to the residuals, on the regressors fitWLS() evaluated, and the replicate is refitted with the weights it gives. The statistic is computed from the Pearson residuals of that refit, so the reference distribution includes the estimation the observed statistic went through. A replicate whose variance function cannot be fitted, where fitWLS() would fall back to equal weights, counts as a failed replicate: the p-value is computed from the replicates that remain, and a warning gives the number that failed.

Holding the weights of a fitWLS() fit fixed, as every call did before 0.15.0 and reestimate_weights = FALSE still does, leaves the estimation out of the replicates, and the test inherits the error of the tests that take the weights as known. In inst/validation/weighted-fits-size.R, the data have standard deviation \(\exp(0.3 x)\), the variance function is correctly specified, the errors are normal and B = 499. At the 5% level the test then rejected 0.086 of the time with 150 observations and 0.105 with 600 when it held the weights fixed, and 0.0605 and 0.0515 when it estimated them again. With the true weights it rejected 0.054 and 0.048. At 2000 replications 0.0605 is 2.2 standard errors above 0.05, and outside the release gate's interval, 0.042 to 0.058. performVarianceFormTest() allows for the estimation too.

References¶

Davidson, R., & Flachaire, E. (2008). The wild bootstrap, tamed at last. Journal of Econometrics, 146(1), 162–169.

Godfrey, L. G. (2006). Tests for regression models with heteroskedasticity of unknown form. Computational Statistics & Data Analysis, 50(10), 2715–2733.

Wu, C. F. J. (1986). Jackknife, bootstrap and other resampling methods in regression analysis. The Annals of Statistics, 14(4), 1261–1295.

Examples¶

data(mtcars)
model <- lm(mpg ~ wt + qsec, data = mtcars)
if (interactive()) {
  set.seed(123)
  performWildBootstrapTest(model, mtcars, B = 199, progress = FALSE)
}