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. A glm is rejected: each replicate is refitted with least squares, which would discard its family and link.

data

Data frame used when fitting model, holding the rows it was fitted to under their row names. 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.

Neither strategy accepts a glm: each replicate is refitted with least squares, so the family and link would be discarded and the replicates would describe a different model from the one supplied.

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 with an offset the model was given as an argument, lm(..., offset = o); an offset written in the formula, offset(o), is evaluated again from the data of the replicate. 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 model is then refitted to the rows of data, and data is refused when that refit does not give the residuals of the fit to a relative 1e-8: rows in another order taken by position, values that differ from those the model was fitted to, or a variable the formula finds outside data, which keeps the order of the fit. Pairs resampling of a model with such a variable, for example offset(o) with o in the environment of the formula, is refused too; give such an offset as an argument instead.

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.

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.

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
}