Run multiple heteroscedasticity diagnostics

Source: man/runHeteroTests.Rd · Help: ?runHeteroTests

Convenience wrapper that executes several heteroscedasticity tests on a fitted linear model.

Usage

runHeteroTests(
  model,
  data = NULL,
  tests = c("white", "koenker"),
  use_cache = TRUE,
  chunk_threshold_mb = Inf,
  chunk_size = 10000,
  progress = interactive()
)

Arguments

model

a fitted lm, or a glm with the Gaussian family and the identity link, a model formula, a tidymodels workflow or a parsnip model fit. The tests refuse a glm of any other family or link since 0.15.0, the deprecated registry names "ncv" and "student_bp" and the substitutes of the fallback (Details) included: the call stops with the error of the first test asked for, without running a substitute and with the class heteroTests_glm_family_error since 0.18.0, and with grouped data each test is returned as failed, with the issue "glm_family". A workflow, or a parsnip fit, is tested on the fit of its engine, which names its response ..y and its predictors as the workflow's preprocessing made them: the data the workflow was fitted to are mapped to those variables by that preprocessing (hardhat::forge()), the outcome as the preprocessing transforms it among them, or the outcome of a fit of parsnip::fit_xy() to ..y, and once their rows are matched to the observations of the fit, every column of its model frame compared, its model frame is tested, as without data. Before 0.17.0 they were refused. Data in another order than the workflow was fitted in are refused, and so are data that cannot be mapped, such as data without a column the preprocessing needs, grouped data, and the data of a workflow whose recipe skips on new data a step that changed the outcome or a predictor at the fit, such as step_log(y, skip = TRUE), which the message names. The case weights of a workflow are compared with those of the fit, on the rows its preprocessing keeps; the data of a parsnip fit of fit_xy() with case weights, whose column the fit does not record, are refused. Each result keeps what generics::augment() needs to map the data in the same way, the blueprint of the preprocessing, the number of rows it kept and the column of the case weights, or the outcome of a parsnip fit, as attr(, "tidymodels"), the blueprint as the workflow keeps it: a recipe step keeps the environment it was written in, so that a saved result carries what that environment holds, as saving the workflow does. The results are otherwise those of the fit of the engine.

data

optional data frame used to fit model. If omitted, model.frame(model) is used, which serves only the tests that find what they need there (see their data argument). With data grouped by dplyr::group_by() the model is refitted within each group by evaluating its own call on that group, so every variable it uses, including weights, offset and subset, must be a column of data. A fit returned by fitWLS() is refused with grouped data since 0.17.0: its call does not hold the weights it estimated, so each group was tested as the unweighted fit of its rows; fit fitWLS() within each group, or pass the data ungrouped. So is a MASS::rlm() fit, which every test refuses. A test that refuses the model of a group as the model context refuses a model, such as a multivariate lm or a penalised mgcv::gam() fit given a test that needs least squares, stops the grouped run with that refusal and its class since 0.17.0, as the run on ungrouped data does; before, the group recorded it as a failed test, whose statistic was NA. The failures of a test for another reason are recorded in its group, as before.

tests

character vector of test names to run: any of the names listDiagnostics() returns, which include those added with registerDiagnostic. On a fit returned by fitWLS() the resampling tests "wild_bootstrap" and "white_bootstrap" estimate the weights again in each replicate, as performWildBootstrapTest() and performWhiteTestBootstrap() do by default since 0.15.0; to hold them fixed, call those functions with reestimate_weights = FALSE.

use_cache

logical: reuse the results of an earlier call on the same inputs? Requires the digest package and defaults to TRUE. Only deterministic tests are cached, under a key that holds the test and the parts of the model and the data it reads (see cachedTest, which also says what the key leaves out): the resampling tests "wild_bootstrap", "white_bootstrap", "rank_permutation" and "spatial_hetero", and diagnostics added with registerDiagnostic, run on every call, and a call that runs one, or falls back on one, is not stored as a whole. A seeded call returns what it returns with use_cache = FALSE, identical(), and leaves the random number generator where it leaves it. Each test is kept in the cache of cachedTest, and each call whose tests are all deterministic in a second cache (clearAnalysisCache); each holds at most getOption("heteroTests.cache_size", 100) results and drops the one used least recently, and options(heteroTests.cache_size = 0) turns both off. The checks the tests make first have a cache of their own, bounded by getOption("heteroTests.validation_cache_size", 1000). Before 0.16.0 the resampling tests were cached too, and two fits with the same coefficients and residuals shared their results.

chunk_threshold_mb

size of data in megabytes above which "white", "breusch_pagan" and "koenker" are run by their streaming forms. The default, Inf, never streams: give a finite value to ask for them, such as 100, or 0 for data of any size. The result of a streaming test has attr(, "extras")$streaming set to TRUE. Up to 0.15.0 the default was 100. A warning about the size of the data is given above half the threshold, or above 50 MB with the default.

chunk_size

integer number of observations processed per chunk when the streaming tests run (chunk_threshold_mb).

progress

logical flag controlling whether textual progress bars are displayed when running multiple diagnostics.

Value

A hetero_test_suite: a named list of htest objects, one for each test asked for, which prints as a table of the tests, and whose summary() is laid out as that of an lm (see summary.hetero_test_suite). With grouped data, a hetero_grouped_suite holding one such suite for each group. Since 0.16.0 a suite whose model the runner fitted itself, from a formula, to the data of a survey design or within a group, says so in attr(, "model_origin") ("formula", "survey" or "group"; see hetero_test_suite).

A result whose data.name is "model" reports the formula of the model instead, as Koenker's test and the others do. Up to 0.15.0 five tests reported "model" here, the name of the argument of the function the registry calls them with: White's test, its bootstrap and its form without cross products, the studentized Breusch-Pagan statistic and Szroeter's test. Any other data.name, such as that of a diagnostic added with registerDiagnostic, is kept. Called directly, a test still names the object it is given: performWhiteTest(m, mtcars)$data.name is "m". runHeteroTestsParallel and cachedTest, which call the tests as the runner does, do the same since 0.18.0.

Details

By default runs performWhiteTest and performKoenkerTest, the studentized Breusch-Pagan test, but additional tests can be requested. Custom diagnostics registered via registerDiagnostic are also available.

Up to 0.14.0 the default was c("white", "breusch_pagan"), whose second test is the classical Breusch-Pagan statistic of performBPTest. That statistic divides by \(2 \hat\sigma^4\), the variance of the squared errors only when they are normal. In the full sweep of inst/validation/, with 150 observations and 400 replications, it rejected a true null hypothesis 28.75% of the time at the 5% level when the errors had a \(t_5\) distribution, against 4.75% for Koenker's statistic. The classical test is still run by name, tests = "breusch_pagan".

When a test stops with an error, the substitute registered for it is run, if one is and it succeeds, and its result is returned under the name of the test asked for: "white" falls back on White's test without cross products, the classical Breusch-Pagan, Koenker and Cook-Weisberg tests, "koenker" on the Cook-Weisberg test. The method of the result names the test that ran, attr(result, "diagnostic") its registry name and attr(result, "extras")$status is "fallback". A weighted fit whose weights are a column of the data, called without data, is such a case: White's and Koenker's tests stop, and both names hold the Cook-Weisberg test. Pass data for such a fit. A refusal that every substitute would make too runs none, and the call stops with it and its class: of the model by the model context, of data whose rows are not those of the fit (heteroTests_data_rows_error; since 0.18.0 also data that lack rows of it, which were tried on every substitute and called a data-quality issue), and, since 0.18.0, of a glm that is not a least-squares fit (heteroTests_glm_family_error). These two refusals, of data that lack rows and of such a glm, are no longer simpleErrors since 0.18.0; catch them by their class or as an error. Nor is a diagnostic registered under the name of a substitute, as by registerDiagnostic("breusch_pagan", f), run after such a refusal. An auxiliary regression that the design forces to fit the squared residuals exactly, as on a paired design, runs no substitute either when the regressors alone fit them, or when that is not known: refused by White's test since 0.16.0 and, since 0.18.0, by the Breusch-Pagan, Koenker, studentized Breusch-Pagan and wild bootstrap tests (see Exact fit in performKoenkerTest()), where "breusch_pagan", "koenker", "student_bp" and "wild_bootstrap" returned the statistic of that fit. The call stops with its message, as an error without its class, as on White's refusal. A substitute that refuses such a fit ends the search, so that "high_dimensional" on fewer than 25 observations and "spatial_hetero" without listw, which fall back on the Breusch-Pagan test, no longer return its statistic of the fit, and stop with their own refusal. White's refusal of as many terms as observations when the regressors alone do not fit them does not end it: the next substitute is tried, as before.

The runner runs the exact tests unless chunk_threshold_mb asks for the streaming ones, performWhiteTestStreaming(), performBPTestStreaming() and performKoenkerTestStreaming(). Up to 0.15.0 it switched to them by itself on data above 100 MB, and performWhiteTest(), performBPTest() and performKoenkerTest() did so on their own. Streaming saves memory mainly in White's test, whose auxiliary design is several times wider than the model matrix; the auxiliary design of the Breusch-Pagan and Koenker tests is the model matrix itself. And until 0.17.0 the streaming White test stopped on any model with a factor or a poly() term, whose squares are collinear with the other terms of its auxiliary regression, so that above 100 MB such a model got a fallback under the name "white" where smaller data got White's test; since 0.17.0 it drops those terms, as performWhiteTest() does.

References

White, H. (1980). A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity. Econometrica, 48(4), 817–838. Breusch, T. S., & Pagan, A. R. (1979). A Simple Test for Heteroscedasticity and Random Coefficient Variation. Econometrica, 47(5), 1287–1294. Koenker, R. (1981). A note on studentizing a test for heteroscedasticity. Journal of Econometrics, 17(1), 107–112. doi:10.1016/0304-4076(81)90062-2

Examples

data(mtcars)
m <- lm(mpg ~ wt + qsec, data = mtcars)
# Default: White and Koenker
runHeteroTests(m, mtcars)
# Laid out as summary.lm() is, with the significance codes
summary(runHeteroTests(m, mtcars))
# Specify additional diagnostics, or the classical Breusch-Pagan test
runHeteroTests(m, mtcars, tests = c("white", "koenker", "cook_weisberg"))
runHeteroTests(m, mtcars, tests = "breusch_pagan")
custom <- function(model, data) list(stat = 1)
registerDiagnostic("custom", custom)
runHeteroTests(m, mtcars, tests = c("white", "custom"))