Bootstrap a test statistic¶

Source: R/bootstrap_methods.R · Help: ?bootstrapTestStatistic

Provides a consistent wrapper around bootstrap resampling for test statistics produced by heteroscedasticity diagnostics. The routine refits the supplied model on each bootstrap replicate and evaluates test_function on it, returning the bootstrap distribution, a percentile interval, and—for null-imposed resampling—a p-value.

Usage¶

bootstrapTestStatistic(
  test_function,
  model,
  data,
  B = 999,
  parallel = FALSE,
  ci_level = 0.95,
  resample = c("null", "pairs"),
  n_cores = NULL,
  progress = interactive(),
  ...,
  reestimate_weights = NULL
)

rbootstrap_test_statistic(
  test_function,
  model,
  data,
  B = 999,
  parallel = FALSE,
  ci_level = 0.95,
  resample = c("null", "pairs"),
  n_cores = NULL,
  progress = interactive(),
  ...,
  reestimate_weights = NULL
)

Arguments¶

test_function: A function that accepts (model, data, ...) and returns an object with a numeric statistic element (typically an htest result).
model: A fitted model of class lm, with or without weights, or, since 0.17.0, a glm with the Gaussian family and the identity link, the least-squares fit of lm() with its prior weights, whose replicates are those of that lm. Another glm is rejected: each replicate is refitted with least squares, which would discard its family and link. A Gaussian mgcv::gam() fit, a glm by class, is accepted since 0.17.0 when it is the least-squares fit of its design, as one whose smooth terms are all fixed (s(x, fx = TRUE)) or one without smooth terms, and its replicates are those of the lm() of that design; one with penalised smooth terms is refused, as every gam was before.
data: Data frame used when fitting model, holding the rows it was fitted to, or its model frame. Other rows are left out. See Details.
B: Integer, number of bootstrap replications. Defaults to 999. The p-value is (1 + #) / (B_eff + 1), so a test at level alpha is exact only when alpha * (B + 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, compute replicates in parallel when possible? See Parallel replicates.
ci_level: Confidence level for the percentile interval.
resample: Bootstrap strategy. "null" (the default) regenerates the response under homoscedasticity and yields a testable p_value; "pairs" samples rows of data with replacement and yields replicates and an interval but no p-value. See Details.
n_cores: Optional integer specifying the number of worker processes to use when parallel = TRUE. Defaults to parallel::detectCores() - 1. The result does not depend on it.
progress: Logical flag controlling whether a textual progress bar is displayed during bootstrap replication.
...: Additional arguments forwarded to test_function.
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. It comes after ..., so it has to be named in full. See Details.

Value¶

A list with components:

original: Result returned by test_function on the original data.
original_statistic: Numeric value of the test statistic from the original sample.
replicates: Numeric vector of bootstrap statistics.
ci: Percentile interval of the bootstrap distribution of the statistic, at ci_level. This summarises where the resampled statistic falls; it is not a confidence interval for a parameter and carries no coverage guarantee.
p_value: Empirical p-value from the null-imposed replicates, computed as \((1 + \#\{T^*_b \ge T\}) / (B_{\mathrm{eff}} + 1)\); NA when resample = "pairs".
B: Requested number of replications.
effective_samples: Number of non-missing bootstrap replications.
call: Matched call for reproducibility.

Details¶

The resample argument decides what the replicates mean.

With resample = "null" (the default) the response is regenerated as \(\hat y_i + e^*_i\), where the \(e^*_i\) are drawn with replacement from the model's residuals after dividing by \(\sqrt{1 - h_i}\) and centring. Those data satisfy homoscedasticity by construction, so the replicates estimate the null distribution and p_value is a test. The leverage correction matters: OLS residuals have variance \(\sigma^2 (1 - h_i)\), and resampling them unscaled under-disperses the regenerated errors, which inflates the rejection rate.

With resample = "pairs" rows are sampled with replacement. The replicates then follow the statistic's distribution under whatever variance structure the data actually have, which is useful for describing its variability but is not a null distribution: the replicates are centred on the observed statistic, so comparing the two returns roughly 0.5 whatever the data. Measured rejection under \(\sigma_i = x_i^2\) was 0%. p_value is therefore NA for pairs resampling, and was removed rather than kept as a number that cannot be interpreted.

A weighted lm is resampled on the scale where its null hypothesis is constant variance, the regression of \(\sqrt{w_i}\, y_i\) on \(\sqrt{w_i}\, x_i\). Null-imposed resampling draws the \(e^*_i\) from the Pearson residuals \(\sqrt{w_i}\, e_i\), divided by \(\sqrt{1 - h_i}\) with the leverages of that regression and centred, and divides each draw by \(\sqrt{w_i}\) before adding it to \(\hat y_i\). Pairs resampling draws each row with its weight. Either way each replicate is refitted with the weights, so test_function sees a weighted fit, as it did on the original data.

Weights given to the fit are known, and they are held fixed. On a fit returned by fitWLS() they were estimated from the same data, and since 0.15.0 each replicate estimates them again by default (reestimate_weights), so that test_function sees what fitWLS() would have returned for the replicate. The replicate is fitted as the model fitWLS() was given, without weights or with the weights that model had. The variance function fitWLS() recorded is fitted again to its residuals, with the same form and on the regressors fitWLS() evaluated, whose rows follow the rows drawn under pairs resampling, and the replicate is refitted with the weights it gives. A replicate whose variance function cannot be fitted, where fitWLS() would fall back to equal weights, fails as a replicate whose refit fails does, and the p-value and the interval are computed from the replicates that remain; a warning gives the number that failed. In inst/validation/weighted-fits-size.R, with the variance function correctly specified, normal errors, 150 observations, B = 499 and performKoenkerTest() as test_function, the rejection rate at the 5% level was 0.083 with the weights held fixed, as every call did before 0.15.0 and reestimate_weights = FALSE still does, 0.045 with them estimated again, and 0.045 with the true weights. See the section on weighted fits in performWildBootstrapTest().

Weighted or not, each replicate is refitted by least squares on the design stored with the fit, in the rows of the replicate, with the weights and the total offset of those rows, the offset argument and the offset() terms. Nothing is evaluated again from the formula or from data: a term computed from the rows it is given, such as cut(x, quantile(x)), the knots of splines::ns() or I(x > median(x)), keeps in every replicate the values the fit computed, and a variable the formula found outside data, such as offset(o) with o in its environment, follows the rows drawn. Under pairs resampling, a level of a factor or a value of a character regressor that no row drawn holds is left out of the replicate, as a refit of the formula leaves it out, and the contrasts of the fit are kept, but for a contrasts matrix of such a factor, which gives way to the default contrasts. test_function is given the replicate as an lm, whose model frame is that of the fit in the rows of the replicate with the response of the replicate, and the rows of data, regenerated or drawn. Up to 0.16.0 each replicate was refitted from the formula on the rows of data, and a fit whose model frame could not be rebuilt from them so was refused ("data does not reproduce the fit"), as with cut(x, quantile(x)), a term computed from all the rows on a fit that dropped some, or a model that transforms a variable given its model frame as data; so was pairs resampling of a model with a variable outside data. 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 the formula rebuilt.

The rows of data are matched to those of the fit by row name, so that each row is regenerated, or drawn, with its own fitted value, weight and offset. Rows the model was not fitted to, such as those it dropped for missing values, are left out, for the observed statistic as for the replicates; the rows that stay keep the order of data. data that lacks a row the model was fitted to is refused. Data without row names of its own, such as a tibble or a data.table, that has as many rows as the fit is taken in the order of the fit. The rows are then matched to the observations of the fit as the tests match them (see test_results), and data whose rows cannot be shown to hold them are refused, as are values that differ from those the model was fitted to. A model whose least-squares refit of its design does not give its residuals, such as an lm fitted with a smaller tol that keeps a column lm() would alias, is refused too ("data does not reproduce the fit"): its replicates would be those of another model. Since 0.17.0 a MASS::rlm() fit, an M-estimate, is refused first, as by every test (see Supported models); it was refused here before, with that message.

Null-imposed resampling additionally requires the response to be a plain variable: for log(y) ~ x the fitted values are on the log scale while the data column holds y, so regenerating one from the other and refitting would take the logarithm twice. That case raises an error; resample = "pairs" resamples rows and is unaffected by it.

The warning White's, Koenker's and the Breusch-Pagan tests give since 0.16.0 on a factor with more levels than observations per level is about their asymptotic p-value (see Fixed effects for small groups in performKoenkerTest()). The observed statistic and the replicates are computed here without it, so that performWhiteTestRobust() and performBPTestRobust() with bootstrap = TRUE give it once, for the asymptotic p-value they report. The result's original, the test's own result, carries that asymptotic p-value: on such a design it is unreliable, and the bootstrap's, p_value, is the one to use.

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 n_cores and whether the workers are forked (parallel::mclapply()) or run on a socket cluster (parallel::makeCluster(), where the platform cannot fork). 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.

Supported models¶

It accepts an lm, unweighted or with positive weights, a fitWLS() fit, since 0.16.0 a rank-deficient lm, and since 0.17.0 a glm with the Gaussian family and the identity link, weighted or not, bootstrapped as the lm it is; a replicate whose refit aliases a coefficient the fit estimated, as a pairs resample of an interaction can, fails. Since 0.17.0 it accepts an mgcv::gam() fit that is the least-squares fit of its design, as one whose smooth terms are all fixed (s(x, fx = TRUE)) or one without smooth terms, bootstrapped on that design, the basis of its smooths, as the lm() of it (decision 38 of the design of the model context); it refused every gam before. It refuses a glm of another family or link and a survey::svyglm() fit, as each replicate is refitted by least squares, an mgcv::gam() fit with penalised smooth terms, which is not the least-squares fit of its design, with the class of the refusals of the model context since 0.17.0 (decision 37), a fit with a zero weight, a MASS::rlm() fit, an M-estimate, with the message of every test since 0.17.0 (it was refused before as a fit the least-squares refit does not reproduce), a fit whose coefficients, residuals or weights were changed after it was fitted (since 0.17.0), a fitWLS() fit of an mgcv::gam() fit with smooth terms, whose terms do not rebuild the design of its coefficients (since 0.17.0), and data whose rows cannot be shown to hold the observations of the fit. Which model each test accepts, and why, is set out in test_results, under Supported models.

Deprecated name¶

rbootstrap_test_statistic() is deprecated since 0.14.0: use bootstrapTestStatistic(), the same function under the package's naming scheme. The old name still works, returns what it returned and warns.

Examples¶

data(mtcars)
model <- lm(mpg ~ wt + cyl, data = mtcars)
if (interactive()) {
  set.seed(123)
  boot <- bootstrapTestStatistic(performWhiteTest, model, mtcars, B = 100)
  boot$ci
}