Bootstrap White test¶

Source: R/additional_tests.R · Help: ?performWhiteTestBootstrap

Approximates the finite-sample distribution of White's LM statistic by resampling the fitted residuals, refitting the model, and recalculating the test statistic across many bootstrap replications.

Usage¶

performWhiteTestBootstrap(
  model,
  data = NULL,
  B = 999,
  parallel = FALSE,
  n_cores = NULL,
  reestimate_weights = NULL
)

Arguments¶

model: A fitted stats::lm object representing the mean specification to be diagnosed.
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.
B: Integer scalar giving the number of bootstrap replications. Defaults to 999. The p-value is (1 + #) / (B_eff + 1), where B_eff is the number of replicates whose statistic is finite (B when every replicate succeeds), so a test at level alpha is exact only when alpha * (B_eff + 1) is an integer: with 999 it is at the 5% and 1% levels, and with 1000, the default before 0.15.0, at neither.
parallel: Logical scalar; if TRUE the replicates are computed in parallel, by forked workers (parallel::mclapply()) where the platform can fork and on a socket cluster (parallel::makeCluster()) elsewhere. See Parallel replicates.
n_cores: Optional integer number of worker processes when parallel = TRUE. Defaults to parallel::detectCores() - 1. The result does not depend on it.
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¶

A htest object reporting both the asymptotic White statistic and its bootstrap p-value. The returned object also stores the simulated test statistics in boot_statistics for further inspection.

Details¶

For a fitted model \(\hat{y} = X\hat{\beta}\), the algorithm proceeds as follows:

Compute White's LM statistic \(n R^2\) from the auxiliary regression on squares and cross-products of the regressors.
Generate \(B\) bootstrap samples by resampling the centred residuals with replacement, forming \(y^{*(b)} = \hat{y} + \hat{e}^{*(b)}\).
Refit the model to each bootstrap sample and recompute White's statistic \(T^{*(b)}\) using the same auxiliary specification.
Estimate the bootstrap p-value as \(\hat{p} = (1 + \sum_{b} I\{T^{*(b)} \ge T_{\text{obs}}\}) / (B_{\text{eff}} + 1)\), where \(b\) runs over the \(B_{\text{eff}}\) replicates whose statistic is finite, so the p-value is never zero.

The bootstrap distribution offers improved size control for moderate sample sizes or high-dimensional designs where the chi-squared approximation may be inaccurate. When parallel = TRUE the resampling step exploits available CPU cores to reduce computation time.

Parallel replicates¶

With parallel = TRUE replicate i draws its random numbers from stream i of a L'Ecuyer-CMRG generator (parallel::nextRNGStream()), seeded by one number drawn from the session's generator. A call made after set.seed() therefore returns the same result each time, whatever the number of workers and whether they are forked or run on a socket cluster. It does not return the result of the same call with parallel = FALSE, which draws every replicate from the session's generator, as it did before 0.15.0: the two are bootstraps of the same data with different draws, and each is valid. Before 0.15.0 a seeded call with parallel = TRUE gave a different result each time, and on Windows more than one worker was an error.

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 is carried out in that regression. What is resampled is its residuals, the Pearson residuals \(\sqrt{w_i}\, e_i\). When the model has an intercept they are centred, as the residuals of an unweighted fit with an intercept are, because the weighted regression has no constant column. Each draw is divided by \(\sqrt{w_i}\) to return it to the scale of the observation it is added to, and each replicate is refitted with the weights, so performWhiteTest() computes its statistic from the Pearson residuals as it does for the fit supplied. With equal weights the result is that of the unweighted fit.

Weighted or not, each replicate is refitted on the stored design of the fit with its total offset, the offset argument and the offset() terms (see The design of the replicates). A transformed response, such as log(y) in log(y) ~ x, is regenerated on the scale the model uses.

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, the variance function fitWLS() recorded is fitted again to its residuals, on the same regressors and with the same form, and the replicate is refitted with the weights it gives, so that performWhiteTest() sees a fitWLS() fit of the replicate. A replicate whose variance function cannot be fitted, where fitWLS() would fall back to equal weights, counts as a failed replicate, and the p-value is computed from the replicates that remain; a warning gives the number that failed.

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, there are 150 observations and B = 499. At the 5% level the test then rejected 0.064 of the time when it held the weights of fitWLS() fixed, as every call did before 0.15.0 and reestimate_weights = FALSE still does. It rejected 0.0415 when it estimated them again, just below the release gate's interval, 0.042 to 0.058, and 0.054 with the true weights. The section on weighted fits in performWildBootstrapTest() has the larger error of Koenker's statistic.

The design of the replicates¶

Since 0.17.0 the observed statistic and every replicate are computed on the design stored with the fit, the design performWhiteTest() reads: each replicate, the fitted values plus the errors drawn, is refitted by least squares on that design, with the weights and the total offset of the fit, and its statistic is White's on the same design, without the columns its refit aliased. Nothing is evaluated again from data or the formula. Up to 0.16.0 both were computed on a design rebuilt from the rows of data, each replicate refitted from the formula: cut(x, quantile(x)) computed other breaks on the rows of the fit and stopped the test ("factor ... has new levels"); a term computed from all the rows of the data, such as I(x > median(x)), took other values on a fit that dropped rows; a fit whose weights or offset are a column of the data, or with an offset() term or a transformed variable, stopped without data; and an mgcv::gam() fit, whose residuals carry no names, given its data in another order, had each draw added to the fitted value of another row. The random numbers are drawn as before, in the same order, so a call made after set.seed() returns what it returned wherever the stored design is the design rebuilt from data.

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 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 (before, the refit from the formula stopped on the smooth term), and it tests a gam that is that fit, as one whose smooth terms are all fixed (s(x, fx = TRUE)), on that design, with the value and the replicates of the lm() of it; a gam without smooth terms is tested as the Gaussian glm it is (decisions 37 and 38 of the design of the model context). It also 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 a gam 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. Which model each test accepts, and why, is set out in test_results, under Supported models.

References¶

Efron, B., & Tibshirani, R. J. (1993). An Introduction to the Bootstrap. Chapman & Hall/CRC.

Davidson, R., & MacKinnon, J. G. (2006). The power of bootstrap and asymptotic tests. Journal of Econometrics, 133(2), 421–441. doi:10.1016/j.jeconom.2005.02.002

Examples¶

# The bootstrap needs at least 50 observations, so mtcars (32) is too small.
set.seed(42)
n <- 120
sim <- data.frame(x = runif(n, 1, 5))
sim$y <- 1 + 2 * sim$x + rnorm(n, sd = 0.3 + 0.6 * sim$x)
mod <- lm(y ~ x, data = sim)
performWhiteTestBootstrap(mod, sim, B = 199)


# More replications give a smoother p-value estimate
performWhiteTestBootstrap(mod, sim, B = 499)

# parallel = TRUE runs the replicates on n_cores workers, forked where the
# platform can fork and on a socket cluster elsewhere. Two is the most
# R CMD check --as-cran allows, and the result does not depend on it.
set.seed(42)
performWhiteTestBootstrap(mod, sim, B = 199, parallel = TRUE, n_cores = 2)