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

The data frame model was fitted to, or NULL, the default. The response, design, weights and offset are read from the fit itself, so a model that transforms its variables, has an offset() term, or takes its weights or offset from a column of the data needs no data, and the model frame of model is accepted as data (since 0.17.0; before, such models needed data). Data that are given are matched to the observations of the fit: rows with row names of their own by name, in any order; rows that only number themselves (a tibble, a data.table, a data frame whose row names were reset) by the row numbers of the fit, or in the order of the fit when they hold one row per observation. The rows the fit dropped may be there or not. Data whose rows cannot be shown to hold the observations of the fit are refused; see test_results. The test reads nothing else from data: its replicates are refitted on the stored design.

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. Its parameter holds B and df, the degrees of freedom of performKoenkerTest(), the rank of the auxiliary regression less its intercept, which the p-value does not use. Up to 0.17.0 df was the number of regressors, one more on a model without an intercept whose columns span the constant, such as lm(y ~ 0 + group + x1). The statistic, for the data and for each replicate, is that of performKoenkerTest(), computed on its auxiliary regression since 0.18.0. Up to 0.17.0 that regression was fitted by lm() from a formula, which dropped a regressor far from zero that a model without an intercept or a Gaussian glm keeps: on lm(y ~ 0 + I(x1 + 1e8) + x2) the statistic was 2.372, that of x2 alone, where Koenker's test gives 13.92, and the p-value was that of the wrong statistic. On other designs the statistic and the replicates are the ones they were. A model whose regressors are all constant, such as lm(y ~ 0 + one), or with an intercept alone, lm(y ~ 1), is refused, since 0.18.0, as Koenker's test refuses it: its statistic is 0 in every replicate, and the test returned 0 with a p-value of 1. Since 0.18.0 it also refuses, as performKoenkerTest() does (see its Exact fit), data whose auxiliary regression fits the squared residuals exactly, as on a paired design, where every replicate does too: on 30 subjects of two observations it gave 60, the number of observations, with a p-value of 1. The refusal comes before the replicates are drawn. Only a fit the design forces is refused; one the data make exact is bootstrapped, as up to 0.17.0.

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)\). The test holds its level with Gaussian and heavy-tailed \(t_5\) errors, and with skewed errors only from about 200 observations; see Size. 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.

Size

In inst/validation/wild-multiplier-size.R and inst/validation/wild-no-multiplier-size.R, with B = 999, 5000 replications per cell and the design lm(y ~ x1 + x2) with x1 uniform on (1, 5) and x2 standard normal, the test rejected a true null hypothesis at the 5% level as often as this. With B = 999 a test whose level is exact rejects 4.9% of the time, and the release gate's interval is 0.042 to 0.058.llrrrr Errors Multiplier n = 30 n = 50 n = 100 n = 200

Gaussian Rademacher 0.0442 0.0476 0.0482 0.0448

Gaussian Mammen 0.0430 0.0508 0.0488 0.0470

\(t_5\) Rademacher 0.0472 0.0444 0.0492 0.0438

\(t_5\) Mammen 0.0460 0.0454 0.0476 0.0460

centred exponential Rademacher 0.0820 0.0734 0.0608 0.0526

centred exponential Mammen 0.0782 0.0698 0.0586 0.0534

With skewed errors, rexp(n) - 1, either multiplier needs about 200 observations to hold the level: both are outside the gate at 30, 50 and 100 and inside it at 200. The excess comes from the multiplier: without it the rate falls, by 1.7 points at 30 and 50 observations. Plausibly, the Rademacher multiplier makes the bootstrap errors symmetric, so the reference distribution loses the skewness of the residuals, on which the squared residuals of the refit depend, and the Mammen multiplier keeps their third moment but doubles their fourth. A plain residual bootstrap, the same resampling without a multiplier (\(v_i = 1\)), which already imposes the null hypothesis, rejected 0.0648, 0.0562, 0.0482 and 0.0466 of the time with the same errors at 30, 50, 100 and 200 observations, holding the level from 50, and 0.042 to 0.047 with Gaussian and \(t_5\) errors. It costs a little power in small samples: against \(\sigma_i^2 = \exp(0.4 x_{1i})\) with Gaussian errors it rejected 0.180, 0.354 and 0.744 of the time at 30, 50 and 100 observations, where the Rademacher multiplier rejected 0.194, 0.361 and 0.745, about 1.4 points less at 30. It will be available as distribution = "none" in 1.0.0, with "rademacher" kept as the default.

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.

Supported models

It accepts an lm and a glm with the Gaussian family and the identity link, unweighted or with positive weights, a fitWLS() fit and, since 0.16.0, a rank-deficient fit, which it tests on the estimable part of its design: it gives the value of the same model fitted without the aliased columns. It refuses a glm of another family or link, a survey::svyglm() fit, a fit with a zero weight, a coefficient that is not finite for another reason, and data that lack rows or variables of the fit. Since 0.17.0 it refuses a fit whose coefficients, residuals or weights were changed after it was fitted, from which it computed the value of the fit before the edit, and a fitWLS() fit of an mgcv::gam() fit with smooth terms, whose terms do not rebuild the design of its coefficients, on which it computed a value on the variables of the formula. Since 0.17.0 it refuses an mgcv::gam() fit with penalised smooth terms, which is not the least-squares fit of its design, the basis of its smooths, that each replicate fits again, and it tests a gam that is that fit, as one whose smooth terms are all fixed (s(x, fx = TRUE)) or one without smooth terms, on that design, its leverages those of the design, with the value and the replicates of the lm() of it (decisions 37 and 38 of the design of the model context). Before, it stopped on every gam ("lm object does not have a proper 'qr' component"), reading its leverages with hatvalues(). Which model each test accepts, and why, is set out in test_results, under Supported models.

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)
}