Package {heteroTests}

Title:	Heteroscedasticity Diagnostics for Linear Models
Version:	0.16.0
Description:	Provides a unified set of heteroscedasticity diagnostics for linear-model workflows. It implements classical auxiliary-regression tests, including those of White (1980) <doi:10.2307/1912934>, Breusch and Pagan (1979) <doi:10.2307/1911963>, Koenker (1981) <doi:10.1016/0304-4076(81)90062-2>, Goldfeld and Quandt (1965) <doi:10.1080/01621459.1965.10480811> and Harvey (1976) <doi:10.2307/1913974>; the score test of Cook and Weisberg (1983) <doi:10.1093/biomet/70.1.1>; the ARCH test of Engle (1982) <doi:10.2307/1912773>; and group-wise tests of equal variance, including those of Bartlett (1937) <doi:10.1098/rspa.1937.0109>, Brown and Forsythe (1974) <doi:10.1080/01621459.1974.10482955> and Hartley (1950) <doi:10.2307/2332383>. Resampling and scalable variants, simulation utilities, diagnostic visualisation and remediation helpers share a consistent interface designed for reproducible statistical workflows and integration with common modelling tools.
License:	Apache License (≥ 2.0)
Encoding:	UTF-8
LazyData:	true
RoxygenNote:	7.3.3
Depends:	R (≥ 4.1)
Imports:	MASS, stats, ggplot2, curl, generics, scales, R6, parallel, SuppDists
Suggests:	quickcheck, testthat (≥ 3.0.0), styler, lintr, digest, covr, knitr, rmarkdown, rjtools, gridExtra, shiny, DT, plotly, htmlwidgets, readxl, Matrix, bench, lmtest, plm, withr, car, vartest, mgcv, quantreg, sandwich, spdep, broom, workflows, parsnip, recipes, survey, data.table, dtplyr, dplyr
VignetteBuilder:	knitr
Config/testthat/edition:	3
URL:	https://github.com/DiogoRibeiro7/heteroTests, https://diogoribeiro7.github.io/heteroTests/
BugReports:	https://github.com/DiogoRibeiro7/heteroTests/issues
NeedsCompilation:	no
Packaged:	2026-10-27 17:14:19 UTC; runner
Author:	Diogo Ribeiro [aut, cre]
Maintainer:	Diogo Ribeiro <dfr@esmad.ipp.pt>
Repository:	CRAN
Date/Publication:	2026-10-29 10:00:12 UTC

Diagnostics registry

Description

Maintains a registry of diagnostic functions that can be called by runHeteroTests() and runDiagnostics(). Users can register custom diagnostics via registerDiagnostic().

Usage

.diagnostic_registry

Format

An object of class environment of length 16.

Ecosystem integration helpers

Description

These utilities extend the package to work smoothly with broader R ecosystem tools such as tidymodels workflows, parsnip model objects, survey designs, data.table-backed pipelines and broom tidiers. They are exported only where user interaction is expected; internal helpers remain undocumented.

Usage

.ht_null

Format

An object of class NULL of length 0.

Create a diagnostic object

Description

Create a diagnostic object that exposes test(), plot() and summary() methods. Formulas are accepted for convenience.

Usage

HeteroDiagnostic(model, data = NULL)

Arguments

Value

An object of class HeteroDiagnostic.

Examples

 data(mtcars)
 m <- lm(mpg ~ wt + qsec, data = mtcars)
 d <- HeteroDiagnostic(m, mtcars)
 runDiagnostics(d$model, d$data)

Heteroscedasticity Test Factory

Description

Provides a registry of heteroscedasticity tests with metadata and convenience methods to run them with input validation. This is an alternative to the older environment-based registry used by runHeteroTests().

Deprecated

TestFactory and its shared instance test_factory are deprecated since 0.14.0: they are the package's internal machinery and leave the exports in 1.0.0. They are objects rather than functions, so they cannot warn when used; this page and NEWS.md are the notice. Use registerDiagnostic() to add a test and runHeteroTests() to run one.

Usage

The shared instance test_factory can be used to register new tests and execute them. See register and run_test methods for details.

Methods

Public methods

• TestFactory$register()

• TestFactory$get_available()

• TestFactory$run_test()

• TestFactory$clone()

Method register()

Usage

TestFactory$register(name, func, metadata = list())

Arguments

Returns

Invisibly returns the factory Get available tests

Method get_available()

Usage

TestFactory$get_available(data_type = NULL, min_n = NULL)

Arguments

Returns

Character vector of test names Run a registered test

Method run_test()

Usage

TestFactory$run_test(test_name, model, data, ...)

Arguments

Returns

Result of the test

Method clone()

The objects of this class are cloneable with this method.

Usage

TestFactory$clone(deep = FALSE)

Arguments

Additional heteroscedasticity tests

Description

These functions extend the package with more specialized tests for heteroscedasticity such as a studentized Breusch-Pagan test, a bootstrap version of White's test and the Szroeter ordered test.

Bibliography of Diagnostic Algorithms

Description

References for the statistical tests implemented in heteroTests.

References

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. Goldfeld, S. M., & Quandt, R. E. (1965). Some tests for homoscedasticity. Journal of the American Statistical Association, 60(310), 539–547. Levene, H. (1960). Robust tests for equality of variances. In Contributions to Probability and Statistics (pp. 278–292). Stanford University Press. Brown, M. B., & Forsythe, A. B. (1974). Robust tests for the equality of variances. Journal of the American Statistical Association, 69(346), 364–367. doi:10.1080/01621459.1974.10482955. Fligner, M. A., & Killeen, T. J. (1976). Distribution-free two-sample tests for scale. Journal of the American Statistical Association, 71(353), 210–213. doi:10.1080/01621459.1976.10481517. Bartlett, M. S. (1937). Properties of sufficiency and statistical tests. Proceedings of the Royal Society of London, 160(901), 268–282. doi:10.1098/rspa.1937.0109. Hartley, H. O. (1950). The maximum F-ratio as a short-cut test for heterogeneity of variance. Biometrika, 37(3/4), 308–312. doi:10.2307/2332383. Engle, R. F. (1982). Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation. Econometrica, 50(4), 987–1007. doi:10.2307/1912773. McLeod, A. I., & Li, W. K. (1983). Diagnostic checking ARMA time series models using squared-residual autocorrelations. Journal of Time Series Analysis, 4(4), 269–273. doi:10.1111/j.1467-9892.1983.tb00373.x. Cook, R. D., & Weisberg, S. (1983). Diagnostics for heteroscedasticity in regression. Biometrika, 70(1), 1–10. doi:10.1093/biomet/70.1.1. O'Brien, R. G. (1979). A general ANOVA method for robust tests of additive models for variances. Journal of the American Statistical Association, 74(368), 877–880. doi:10.1080/01621459.1979.10481047. Cameron, A. C., & Trivedi, P. K. (1990). The information matrix test and its applied alternative hypotheses. Working Paper, University of California, Davis. Spearman, C. (1904). The proof and measurement of association between two things. The American Journal of Psychology, 15(1), 72–101. Cameron, A. C. (1991). Regression based tests of heteroskedasticity in models where the variance depends on the mean (Working Paper No. WP-379). University of California, Davis. Greene, W. H. (2018). Econometric Analysis (8th ed.). Pearson. Wooldridge, J. M. (2020). Introductory Econometrics: A Modern Approach (7th ed.). Cengage Learning. Davidson, R., & MacKinnon, J. G. (2004). Econometric Theory and Methods. Oxford University Press. Godfrey, L. G. (1988). Misspecification Tests in Econometrics. Cambridge University Press. Fox, J., & Weisberg, S. (2019). An R Companion to Applied Regression (3rd ed.). Sage Publications. Kleiber, C., & Zeileis, A. (2008). Applied Econometrics with R. Springer-Verlag. Efron, B., & Tibshirani, R. J. (1993). An Introduction to the Bootstrap. Chapman & Hall.

Machine-learning residual analysis

Description

Fits a GAM using mgcv and compares its residuals to the input linear model.

Usage

analyzeMLResiduals(model, data = NULL)

Arguments

Value

A list with the GAM model, residuals and RMSE reduction.

See Also

compareModelDiagnostics

Examples

 data(mtcars)
 analyzeMLResiduals(mpg ~ wt + qsec, mtcars)

Automatically compare remedial models

Description

Fits several remedial models for heteroscedasticity and compares their performance using AIC and residual RMSE. Currently evaluates weighted least squares and robust regression against the original model.

Usage

autoCompareRemediations(model, data = NULL)

Arguments

Value

A list with components metrics, models, and best indicating the recommended method.

Examples

data(mtcars)
m <- lm(mpg ~ wt + qsec, mtcars)
autoCompareRemediations(m)

Automatic transformation helper

Description

Chooses between no transformation, log, square root and Box–Cox based on the Breusch-Pagan test statistic.

Usage

autoTransform(model)

Arguments

Value

A list with elements model (the transformed fit) and method indicating the chosen transformation.

Examples

data(mtcars)
m <- lm(mpg ~ wt + qsec, data = mtcars)
autoTransform(m)

Autoplot heteroscedasticity diagnostics

Description

Draws one horizontal bar per test. The length of a bar is the evidence against the null hypothesis, -\log_{10} p, on an axis labelled with the p-values themselves, so a smaller p-value gives a longer bar. A dashed line marks alpha; bars that pass it are drawn in the accent colour and the others in grey. The p-value of each test is printed beside its bar, so no reading depends on colour.

Usage

## S3 method for class 'hetero_test_suite'
autoplot(object, ..., alpha = 0.05)

## S3 method for class 'hetero_grouped_suite'
autoplot(object, ..., alpha = 0.05)

Arguments

Details

Bars stop at p = 10^{-4}: extra length beyond that would only compress the tests that are near the threshold. Labels give the p-value to three decimals, and p < 0.001 below that. A test that failed has no bar and is labelled as such.

For a grouped suite there is one panel per group.

Value

A ggplot object.

Earlier versions

Up to 0.12.0 the bars were the p-values themselves on a linear axis from 0 to 1, so a test that rejected had no visible bar, and the highlighting of significant tests was never applied to a suite of more than one test.

Examples

fit <- lm(mpg ~ wt + qsec, data = mtcars)
suite <- runHeteroTests(
  fit, mtcars,
  tests = c("white", "breusch_pagan", "koenker"),
  progress = FALSE
)
ggplot2::autoplot(suite)

Bootstrap a test statistic

Description

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

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.

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.

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.

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.

Supported models

It accepts an lm, unweighted or with positive weights, a fitWLS() fit and, since 0.16.0, a rank-deficient lm; a replicate whose refit aliases a coefficient the fit estimated, as a pairs resample of an interaction can, fails. It refuses every glm, as each replicate is refitted by least squares, a fit with a zero weight, and data that do not reproduce the fit. Which model each test accepts, and why, is set out in test_results, under Supported models.

Examples

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

Bootstrap utilities for heteroscedasticity diagnostics

Description

These helpers provide reusable bootstrap and power-estimation routines for the heteroscedasticity testing framework. They support the robust implementations introduced in robust_implementations.R but are written generically so additional tests can reuse them.

Bootstrap implementation for White test

Description

Bootstrap implementation for White test

Usage

bootstrap_white_test(model, data, B, parallel, cross_products)

Cache results of a heteroscedasticity test

Description

Runs a registered diagnostic, and returns the result of an earlier call on the same inputs instead when the diagnostic is deterministic.

Usage

cachedTest(test_name, model, data, ..., use_cache = TRUE)

Arguments

Details

What is cached. Only the deterministic tests the package registers: "white", "breusch_pagan", "koenker", "cook_weisberg", "ncv", "spread_level", "student_bp", "szroeter", "quantile_regression", "high_dimensional" and "variance_form".

What is not. A test that draws random numbers is never cached: "wild_bootstrap", "white_bootstrap", "rank_permutation" and "spatial_hetero" run on every call, so a seeded call returns what it returns without the cache, and leaves the random number generator where it leaves it without the cache. Nor is "box_m", which tests data rather than a model, a diagnostic added with registerDiagnostic() or test_factory$register(), or one registered again under a name of the package: the package cannot know whether it draws random numbers or reads anything besides its inputs. A result whose computation moved the random number generator is not stored either. Before 0.16.0 every name was cached, the resampling tests included: a call with another seed returned the replicates of the first, and the random numbers drawn after it were not those of a call without the cache.

The key. The test, its arguments, and the parts of the model and of the data the test reads: the class of the model, its family and link; from the model context of the fit, its model frame, response, design, residuals and fitted values, prior weights, total offset, rank and aliased coefficients, the names of its observations and its na.action; its coefficients, the weights it stores, its terms (with their predvars), call, contrasts and levels, residual degrees of freedom, whether it keeps its model frame, and its attributes, such as the variance function of a fitWLS() fit; and data, whole: its values, row names, class and attributes. The key is an xxhash64 digest of those values (from the digest package; without it nothing is cached). It holds values, not environments: a model refitted to the same data, in another environment or another session, has the same key and is answered from the cache, with the model of the call as attr(result, "model"), so that the result is identical() to the one the test computes. A model the tests refuse before reading it, such as a fit with model = FALSE whose data changed after the fit, has no key and is refused as without the cache; so has a fit with model = FALSE whose data are gone, whose model frame the tests rebuild from data, and it runs without the cache.

What the key leaves out. The key does not hold the environment of the formula: a test that evaluates the formula on data again, as the check of the rows of data does, reads the functions and variables of that environment as they are, so that after f is defined again a model y ~ f(x) is answered with the result computed with the old f, where without the cache the check may refuse data. Clear both caches, with clearTestCache() and clearAnalysisCache(), after changing one of them. Nor does it hold the options of the session: a result read from the cache gives none of the warnings and messages of the call that computed it but one (below), so a call made under options(warn = 2) after one that stored its result is answered with that result, where without the cache a warning of the test, turned into an error, stops it. The results of the tests otherwise do not depend on the options; the names of the estimates of the quantile regression test depended on options(digits) before 0.16.0.

Before 0.16.0 the key held only the test, the coefficients, the residuals, data and the arguments. Fits that share those shared one entry, and the second call returned the result of the first: lm(y ~ x + offset(o)) and lm(I(y - o) ~ x), which have the same coefficients and residuals and other fitted values; a fit and the same fit with its weights changed, which the tests refuse; or a fit with model = FALSE before and after its data changed, which they refuse too.

The bound. The cache holds at most getOption("heteroTests.cache_size", 100) results, and drops the one used least recently to store another; options(heteroTests.cache_size = 0) turns it off, together with the cache of runHeteroTests(). A stored result keeps no copy of the model. Before 0.16.0 the cache grew without bound, by about 25 KB per call on a new model, and kept about 0.3 KB per call after clearTestCache(). The checks the tests make first keep their results in a cache of their own, of at most getOption("heteroTests.validation_cache_size", 1000) results, which 0 turns off; it holds no check that drew random numbers.

A result read from the cache is returned without running the test, so the call gives none of the warnings and messages of the call that computed it, but one: since 0.16.0 the warning White's, Koenker's and the Breusch-Pagan tests give on a factor with more levels than observations per level, of class heteroTests_fixed_effects_warning, is stored with the result and given again on every call that reads it, as by runHeteroTests() (see Fixed effects for small groups in performKoenkerTest()). The call that computes the result gives it once the result is stored, so that under options(warn = 2) the warning stops the call after the result is stored.

Value

The diagnostic result, either retrieved from cache or freshly computed.

See Also

clearTestCache() empties the cache; runHeteroTests() uses it with use_cache = TRUE.

Validate data argument (deprecated)

Description

Ensures that data is a data.frame. Used internally for input validation across the package.

Usage

checkData(data)

Arguments

Details

Deprecated since 0.14.0: an internal helper of the package, which leaves the exports in 1.0.0. It still works and warns.

Value

Invisible data if valid, otherwise an error is thrown.

Validate model argument

Description

Checks that model is of class lm or glm. With enhanced = TRUE it also runs the soft checks of the deprecated checkModelEnhanced(): it warns when the residual degrees of freedom fall below six, when a near-perfect fit is detected (R^2 > 0.999), or, when data is given, when a variance inflation factor exceeds 10 (Belsley et al., 1980).

Usage

checkModel(model, data = NULL, enhanced = FALSE)

Arguments

Details

It does not check whether a test accepts model: it accepts a binomial or a Poisson glm, a survey::svyglm() fit and a fit with a zero weight, which the tests refuse, each with its own message. Which test accepts which model is set out in test_results, under Supported models.

The R-squared is that of the least-squares fit: for a glm with the Gaussian family and the identity link, that of the equivalent lm, which summary() of a glm does not report; for any other glm it is not checked. Before 0.16.0, enhanced = TRUE stopped on any glm with "missing value where TRUE/FALSE needed".

Value

Invisible model if valid, otherwise an error is thrown.

References

Belsley, D. A., Kuh, E., & Welsch, R. E. (1980). Regression Diagnostics: Identifying Influential Data and Sources of Collinearity. Wiley.

Enhanced model diagnostics (deprecated)

Description

Deprecated since 0.14.0: use checkModel(model, data, enhanced = TRUE), which runs the same checks. This name still works and warns.

Usage

checkModelEnhanced(model, data = NULL)

Arguments

Details

Extends checkModel() with additional soft checks that flag potential numerical issues before computationally expensive diagnostics are executed. The helper mirrors the warnings used throughout the package and retains the original return semantics for backward compatibility.

The function warns when residual degrees of freedom fall below six, when a near-perfect fit is detected (R^2 > 0.999), or when variance inflation factors (VIFs) exceed 10. The latter threshold follows the regression diagnostics literature (Belsley et al., 1980).

Value

Invisibly returns model after issuing any relevant warnings.

References

Belsley, D. A., Kuh, E., & Welsch, R. E. (1980). Regression Diagnostics: Identifying Influential Data and Sources of Collinearity. Wiley.

See Also

checkModel(), performVIFDiagnostic(), validateTestInputs()

Examples

data(mtcars)
fit <- stats::lm(mpg ~ wt + hp, data = mtcars)
checkModel(fit, enhanced = TRUE)

Validate numeric vector

Description

Ensures x is numeric and non-empty. Used internally by helpers that expect numeric input.

Usage

checkNumericVector(x, name = "x")

Arguments

Value

Invisible x if valid, otherwise an error is thrown.

Check memory usage and warn for large datasets

Description

Check memory usage and warn for large datasets

Usage

check_memory_usage(data, threshold_mb = 100)

Arguments

Clear cached diagnostic run results

Description

Removes the results of whole calls that runHeteroTests stored with use_cache = TRUE, so that later calls run their tests again.

Usage

clearAnalysisCache()

Details

A call is stored only when every test it runs is deterministic, under a key that holds the tests, whether they are streamed and in chunks of what size, and every part of the model and the data they read (see cachedTest). The cache holds at most getOption("heteroTests.cache_size", 100) calls and drops the one used least recently; 0 turns it off. Each test of the call is also kept in the cache of cachedTest, which clearTestCache empties.

Value

NULL, invisibly.

See Also

runHeteroTests, clearTestCache

Clear cached test results

Description

Removes every result from the cache of cachedTest(), which runHeteroTests() also uses for each test with use_cache = TRUE, so that later calls compute them again.

Usage

clearTestCache()

Details

The cache is bounded and holds no test that draws random numbers (see cachedTest()), so clearing it is needed neither to bound the memory of a long session nor to draw new replicates with another seed. It is needed after a change the key does not see, such as a function or a variable of the environment of a model's formula; call clearAnalysisCache() then too, since the cache of whole calls of runHeteroTests() keeps the results of such a model as well.

Value

NULL, invisibly.

See Also

cachedTest() for invoking diagnostics with memoisation.

Reset cached validation results

Description

Clears the internal validation cache so that subsequent checks recompute their diagnostics. This helper is primarily used in automated tests where deterministic behaviour is required.

Usage

clearValidationCache()

Details

The cache holds at most getOption("heteroTests.validation_cache_size", 1000) results and drops the one used least recently; options(heteroTests.validation_cache_size = 0) turns it off. A check that draws random numbers, as the check of normality does on a variable longer than its sample_limit, is not stored.

Compare diagnostics across models

Description

Runs runDiagnostics for multiple models and collects test statistics. Diagnostics that cannot be computed (for example, due to validation failures) emit warnings and fill the corresponding entries with NA_real_. The original error messages are preserved on the returned data frame via the diagnostic_errors attribute.

Usage

compareModelDiagnostics(models, data = NULL, tests = c("white", "koenker"))

Arguments

Value

A data frame of statistics, one row per model. Failed diagnostics are represented by NA_real_ entries and the underlying error messages can be retrieved from the diagnostic_errors attribute.

See Also

analyzeMLResiduals

Examples

 data(mtcars)
 m1 <- lm(mpg ~ wt + qsec, mtcars)
 m2 <- lm(mpg ~ wt + hp, mtcars)
 compareModelDiagnostics(list(m1, m2))

Summarise heteroscedasticity test results

Description

Runs selected diagnostics for a single model and returns a tidy data frame of test statistics and p-values.

Usage

compareTestResults(model, data = NULL, tests = c("white", "koenker"))

Arguments

Value

Data frame with columns test, statistic, and p.value.

Examples

data(mtcars)
m <- lm(mpg ~ wt + qsec, mtcars)
compareTestResults(m)

Simulated dataset with collinearity and mild nonlinearity

Description

This dataset contains two predictors with correlation and a response including a quadratic term. It is useful for testing diagnostics beyond heteroscedasticity.

Usage

data(diagnostic_data)

Format

A data frame with 150 rows and 3 variables:

x1

first predictor

x2

second predictor, correlated with x1

y

response

Source

Simulated via set.seed(123).

Examples

data(diagnostic_data)
head(diagnostic_data)

Download example teaching dataset

Description

Provides convenient access to example datasets hosted online. Currently supports the 'tips' dataset used in various tutorials.

Usage

downloadTeachingData(
  name = "tips",
  destfile = tempfile(fileext = ".csv"),
  quiet = FALSE
)

Arguments

Value

Path to the downloaded dataset on success, or NULL invisibly if the download fails or no internet connection is available.

Examples

# Online example (requires internet)
if (curl::has_internet()) {
  data_path <- downloadTeachingData(quiet = TRUE)
  if (!is.null(data_path)) {
    # Use downloaded data
  }
}


# Offline example using built-in data
data(mtcars)
model <- lm(mpg ~ wt + hp, data = mtcars)
result <- performWhiteTest(model, mtcars)
print(result)

Robust regression wrapper

Description

Provides a simple interface to MASS::rlm for robust estimation.

Usage

fitRobust(model, data = NULL, ...)

Arguments

Value

An object of class rlm.

Examples

data(mtcars)
m <- fitRobust(mpg ~ wt + qsec, mtcars)
summary(m)

Weighted Least Squares wrapper

Description

Refit a model by feasible generalised least squares, weighting each observation by the inverse of an estimated error variance.

Usage

fitWLS(
  model,
  data = NULL,
  var_formula = NULL,
  form = c("exponential", "power")
)

Arguments

Details

The variance is modelled, not read off the residuals directly. A single squared residual is a one-degree-of-freedom estimate of \sigma_i^2 and far too noisy to invert: weighting by 1/e_i^2 hands almost all of the weight to whichever observations the initial fit happened to reproduce most closely. This function instead fits a variance function to the squared residuals and takes \hat\sigma_i^2 from its fitted values, the standard feasible-GLS recipe; weights are 1/\hat\sigma_i^2.

Value

A new lm object fitted with weights. The estimated variances are attached as the "variance_model" attribute, and the fitted variance function (its form, regressors and coefficients) as the "variance_function" attribute, which is what performVarianceFormTest() reads, and what the bootstrap tests read to estimate the weights again in each replicate. When model was fitted with weights, the attribute also holds them, as initial_weights.

Variance functions

form chooses the shape of the variance function and var_formula the variables z_i it depends on. Both forms are linear in the logarithm of the variance and carry their own constant.

"exponential"

\log \sigma_i^2 = \gamma_0 + z_i^\top \gamma, the multiplicative model of Harvey (1976), estimated by regressing \log e_i^2 on z_i. The log scale keeps the fitted variances positive without constraining the auxiliary regression. This is the default, and with var_formula = NULL it is what fitWLS() has computed since 0.9.0.

"power"

\log \sigma_i^2 = \gamma_0 + \sum_j \delta_j \log z_{ij}, that is \sigma_i^2 \propto \prod_j z_{ij}^{\delta_j}, the model of Park (1966). It is the exponential form in the logarithms of the regressors, which must be strictly positive. \delta = 2 is the textbook case of a standard deviation proportional to z.

Other log-linear shapes are written through var_formula, for example ~ x + I(x^2) or ~ sqrt(x). The additive model \sigma_i^2 = \alpha_0 + z_i^\top \alpha is not offered: nothing keeps its fitted variances positive, and in the validation study it was not usable in a material share of samples drawn from an additive model itself (inst/validation/README.md).

Residuals that are numerically zero are floored before the logarithm is taken, with a warning.

Weights estimated this way are consistent under a correctly specified variance function, so standard errors from the returned fit are usable. They were not before 0.9.0, when the weights were the raw inverse squared residuals: the weighted residual sum of squares then collapsed towards n regardless of the data, and nominal 95% intervals covered the truth about 10% of the time.

Whether the variance function is correctly specified can be tested: pass the result to performVarianceFormTest().

A test of constant variance on the result takes its weights as known, although they were estimated from the same data, so its reference distribution is only approximate: see the validation study in inst/validation/README.md. The bootstrap tests performWildBootstrapTest(), performWhiteTestBootstrap() and bootstrapTestStatistic() allow for the estimation. Since 0.15.0 they run fitWLS() again on each replicate by default: they fit it as model was fitted, fit the recorded variance function again to its residuals, on the recorded regressors, and refit it with the weights that gives. Their argument reestimate_weights = FALSE holds the weights fixed instead, as they did before 0.15.0.

If the variance function cannot be fitted, or yields no usable variation, the function falls back to equal weights, which reduces the result to the original OLS fit.

The weighted fit keeps the design and the offset of model, whether the offset was given as an argument, lm(..., offset = o), or written in the formula, offset(o). Its coefficients, fitted values and residuals are those lm() returns for the same call with the estimated weights added. On a model fitted with weights, the refit uses the estimated weights in their place, and weights() of the result returns them. Since 0.16.0 the response, the design and the offset are read from the fit, as the tests read them, and a fit with model = FALSE whose model frame, rebuilt from the data it was fitted to, no longer holds its response, because those data changed after the fit, is refused, where it was refitted to the response of the changed data; a fit with model = FALSE whose data are gone is refitted from data that give its response and fitted values.

A glm with the Gaussian family and the identity link is a least-squares fit, with the design and residuals of the lm, and since 0.16.0 the result for it is the lm returned for that lm, prior weights included, with the call of that lm: stats::lm(), without the arguments only glm() takes, so that update() fits an lm again. Before, it kept the class and the call of the glm without its family, which summary() could not read and every test refused. Any other glm is refused.

References

Harvey, A. C. (1976). Estimating regression models with multiplicative heteroscedasticity. Econometrica, 44(3), 461–465. doi:10.2307/1913974

Park, R. E. (1966). Estimation with heteroscedastic error terms. Econometrica, 34(4), 888. doi:10.2307/1910108

See Also

performVarianceFormTest() to test the variance function that was fitted.

Examples

data(mtcars)
m <- lm(mpg ~ wt + qsec, data = mtcars)
wls <- fitWLS(m)
summary(wls)

# Variance proportional to a power of one regressor
wls_power <- fitWLS(m, var_formula = ~ wt, form = "power")
attr(wls_power, "variance_function")$coefficients

Automated diagnostic report generation

Description

Creates a standalone report summarizing heteroscedasticity tests and diagnostic plots for a fitted model.

Usage

generateDiagnosticReport(
  model,
  data = NULL,
  output_format = "html",
  output_file = NULL,
  include_remediation = TRUE,
  include_theory = FALSE
)

Arguments

Value

Invisibly returns the path to the generated report.

Examples

if (requireNamespace("rmarkdown", quietly = TRUE) &&
    rmarkdown::pandoc_available()) {
  model <- lm(mpg ~ wt + hp, data = mtcars)
  generateDiagnosticReport(
    model, mtcars,
    output_file = file.path(tempdir(), "diagnostic_report.html")
  )
}

Intelligent recommendations for heteroscedasticity diagnostics

Description

These helpers assemble an end-to-end recommendation workflow that profiles the input dataset, selects appropriate heteroscedasticity diagnostics, interprets results with qualitative confidence labels, and suggests remediation strategies with decision-tree style guidance.

Usage

analyseDatasetCharacteristics(data, response = NULL)

suggestDiagnosticsForProfile(profile)

interpretHeteroTestResults(test_results, alpha = 0.05)

recommendRemediationStrategies(test_results, model, data, profile = NULL)

warnDiagnosticAssumptions(profile, interpretations, model, data, test_results)

buildDiagnosticDecisionTree(profile, recommendations)

composeDiagnosticNarrative(profile, interpretations, remediation,
  assumption_warnings, decision_tree)

generateHeteroRecommendations(model, data = NULL, test_results = NULL,
  alpha = 0.05, include_report = TRUE)

Arguments

Details

The recommendation engine analyses dataset characteristics such as sample size, missingness, and predictor skewness to propose diagnostics tailored to the context. It automatically interprets test outcomes, attaches qualitative confidence labels, warns when assumptions are at risk, and suggests remediation strategies matched to detected variance patterns. The generated decision tree and narrative target non-specialist audiences who need actionable guidance.

Value

analyseDatasetCharacteristics() returns a list describing dataset size, missingness, skewness, and risk factors. suggestDiagnosticsForProfile() returns a data frame of recommended diagnostics with rationales. interpretHeteroTestResults() yields a tidy data frame of interpretations and an overall summary. recommendRemediationStrategies() produces a structured remediation plan. warnDiagnosticAssumptions() returns warnings about assumption violations. buildDiagnosticDecisionTree() constructs a tidy decision tree. composeDiagnosticNarrative() outputs a plain-language character vector, and generateHeteroRecommendations() returns an object of class hetero_recommendation_report bundling all components.

See Also

runHeteroTests(), suggestRemediation(), generateDiagnosticReport()

Examples

model <- lm(mpg ~ wt + hp, data = mtcars)
recs <- generateHeteroRecommendations(model, mtcars)
print(recs)

Simulated dataset with heteroscedastic errors

Description

This dataset is generated from the model y_i = 1 + 2 x_i + \varepsilon_i, where \varepsilon_i \sim N(0, (0.5 + 2 x_i)^2). A fixed seed ensures repeatable values.

Usage

data(hetero_data)

Format

A data frame with 100 rows and 2 variables:

x

predictor

y

response

Source

Simulated with set.seed(42); runif() and rnorm().

Examples

data(hetero_data)
plot(hetero_data$x, hetero_data$y)

S3 containers for heteroscedasticity diagnostics

Description

Collections of heteroscedasticity test results returned by runHeteroTests. A hetero_test_suite is a simple named list of htest objects with additional metadata. When diagnostics are evaluated on grouped data, the returned object inherits from hetero_grouped_suite and contains one hetero_test_suite per group along with the grouping keys.

Details

These classes are primarily useful for method dispatch (e.g. the broom tidiers and autoplot() methods). Users typically interact with the objects via print() and summary() (see summary.hetero_test_suite), generics::tidy(), generics::glance(), generics::augment() and ggplot2::autoplot().

Attributes

A hetero_test_suite carries attr(, "tests"), the names asked for; attr(, "model"), the model tested; attr(, "data_source"), the class of the data, or NA; and, since 0.16.0, attr(, "model_origin") when runHeteroTests fitted the model itself: "formula" when it was given a formula, "group" for the refit of the model within a group of a grouped run, and "survey" for the fit to the data of a survey design (runSurveyHeteroTests). The call of such a model names variables of the runner rather than the caller's, so summary() shows its formula instead. A suite of a model you fitted has no model_origin. A hetero_grouped_suite carries attr(, "group_keys") and attr(, "tests").

Log a formatted message

Description

This helper wraps base::message() but prepends a log level for clearer diagnostics when running algorithms. Messages are filtered according to ht_set_log_level() and can optionally be captured with ht_enable_log_capture(). Intended for internal use.

Usage

ht_log(level = c("INFO", "WARN", "ERROR"), msg)

Arguments

Value

NULL, invoked for its side effect.

Control logging for heteroTests diagnostics

Description

These helpers manage the package-level logging subsystem used to debug complex diagnostic workflows. Users can adjust verbosity, enable log capture, inspect recent log entries and clear the stored history when finished.

Usage

ht_set_log_level(level = c("INFO", "WARN", "ERROR", "SILENT"))

ht_enable_log_capture(enabled = TRUE, max_entries = 1000L, sink = NULL)

ht_log_history()

ht_clear_log_history()

Arguments

Value

ht_set_log_level() returns the previous log level invisibly. ht_enable_log_capture() and ht_clear_log_history() return NULL. ht_log_history() returns a data frame with columns timestamp, level and message.

Deprecated

ht_enable_log_capture(), ht_log_history() and ht_clear_log_history() are deprecated since 0.14.0: of the logging functions only ht_set_log_level() stays exported after 1.0.0. They still work and warn.

Interpret White test results

Description

Interpret White test results

Usage

interpret_white_test(result, alpha = 0.05)

Interactive diagnostic dashboard

Description

Launch a production-ready Shiny dashboard for running heteroscedasticity diagnostics, exploring interactive plots, and experimenting with simulation studies.

Usage

launchDiagnosticDashboard(model, data)

Arguments

Details

Requires optional packages shiny, DT, and plotly. Excel uploads additionally depend on readxl and exporting interactive graphics relies on htmlwidgets. The dashboard guides users through data preparation, model fitting, diagnostic testing, interactive Plotly visualisations, and a real-time simulation lab for exploring heteroscedastic data-generating processes. Download buttons are provided for both results tables and interactive plots.

Value

A shiny.appobj that can be run with shiny::runApp().

Examples

if (interactive() && requireNamespace("shiny", quietly = TRUE) &&
    requireNamespace("DT", quietly = TRUE) &&
    requireNamespace("plotly", quietly = TRUE)) {
  mod <- lm(mpg ~ wt + hp, data = mtcars)
  launchDiagnosticDashboard(mod, mtcars)
}

List the registered diagnostics

Description

Returns the names runHeteroTests() accepts in tests: the diagnostics the package registers and any added with registerDiagnostic().

Usage

listDiagnostics()

Details

"box_m" is registered for runMultivariateTests(), which calls it with data and a grouping factor rather than with a model.

Value

A character vector of the registered names, sorted.

See Also

registerDiagnostic(), runHeteroTests()

Examples

listDiagnostics()

Standardized error and warning messages

Description

Provides a central location for common error and warning strings used across the package. The package signals problems through internal versions of std_error() and std_warning(); the exported names are deprecated since 0.14.0.

Modern heteroscedasticity diagnostics

Description

Collection of advanced heteroscedasticity tests including wild bootstrap, heteroscedasticity-consistent covariance LM statistics, quantile regression comparisons, rank-based permutation checks, high-dimensional projections and spatial correlation diagnostics.

Usage

modern_diagnostics

Format

An object of class NULL of length 0.

Perform Engle's ARCH LM test

Description

Regresses squared residuals on their lags to detect ARCH effects.

Usage

performArchLMTest(model, data = NULL, lags = 1)

Arguments

Details

The statistic is n R^2 from an auxiliary regression of e_t^2 on its lagged values, where R^2 is the coefficient of determination. Under the null of no ARCH effects it follows a chi-square distribution with degrees of freedom equal to the number of lags.

Value

An object of class htest containing the test statistic, p-value and degrees of freedom.

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.16.0 it is built on the model context, and also refuses a multivariate lm ("needs a model with one response"), a fit whose coefficients, residuals or weights were changed after the fit, and a fit with model = FALSE whose model frame, rebuilt from the data it was fitted to, no longer holds its response, because those data changed after the fit, with or without data. A fit with model = FALSE whose data are gone is tested on data that give its response and fitted values, and refused without them. The test reads nothing from data, and since 0.16.0 matches it to the observations of the fit as every test does: data with row names of their own in any order, data that number their rows (a tibble, a data.table) in the order of the fit; other numbered data, rows in front of those of the fit and repeated row names are refused, where it gave the value of the fit. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Engle, R. F. (1982). Autoregressive conditional heteroskedasticity with estimates of the variance of United Kingdom inflation. Econometrica, 50(4), 987–1007. doi:10.2307/1912773

Hamilton, J. D. (1994). Time Series Analysis. Princeton University Press.

Examples

 data(mtcars)
 m <- lm(mpg ~ wt + qsec, data = mtcars)
 performArchLMTest(m, lags = 2)

Breusch-Pagan LM test for random effects

Description

Lagrange Multiplier test for the presence of a random individual effect in panel data. This is not a test for heteroscedasticity and does not respond to one; use the auxiliary-regression diagnostics for that.

Usage

performBPRandomEffectsTest(model, data = NULL, id)

Arguments

Details

The statistic is Breusch and Pagan (1980) equation 5, LM = nT / (2(T-1)) [ \sum_i (\sum_t e_{it})^2 / \sum_{it} e_{it}^2 - 1 ]^2, which follows a chi-square distribution with one degree of freedom under the null of no individual effect. The bracketed ratio is close to one under the null and the statistic measures its squared departure from one. Before 0.11.0 the "- 1" and the square were absent and the scaling used T^2 rather than nT, which left the statistic sitting at the critical value: it rejected about a third of the time when no individual effect was present.

Value

An object of class htest.

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. The model is a pooled lm of the panel, with the units and the periods as columns of data; a plm::plm() fit is not an lm and is refused. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Breusch, T. S., & Pagan, A. R. (1980). The Lagrange Multiplier Test and Its Applications to Model Specification in Econometrics. Review of Economic Studies, 47(1), 239–253.

Examples

 df <- data.frame(id = rep(1:5, each = 4), time = rep(1:4, 5), x = runif(20), y = rnorm(20))
 m <- lm(y ~ x, data = df)
 performBPRandomEffectsTest(m, df, "id")

Perform Breusch-Pagan test for heteroscedasticity

Description

Implements the classical Breusch-Pagan (1979) test on a fitted linear model, in which the scaled squared residuals are regressed on the original regressors.

Usage

performBPTest(model, data = NULL)
performBreuschPaganTest(model, data = NULL)

Arguments

Details

The test statistic is half the explained sum of squares from regressing the scaled squared residuals e_i^2/\hat{\sigma}^2 - 1 (with \hat{\sigma}^2 = \sum e_i^2 / n) on the explanatory variables. Under the null hypothesis of homoscedasticity and normal disturbances it follows a chi-square distribution with degrees of freedom equal to the number of regressors, matching lmtest::bptest(..., studentize = FALSE). For the studentized n R^2 form that drops the normality assumption use performKoenkerTest.

The test is computed in memory whatever the size of the data. Up to 0.15.0 it switched to performBPTestStreaming() on data above 100 MB; call that function to stream.

Value

An object of class htest containing the test statistic, p-value and degrees of freedom.

Fixed effects for small groups

The Breusch-Pagan test, like White's and Koenker's tests, is not valid with fixed effects for small groups: a factor with many levels and few observations of each, such as subject in lm(y ~ x + subject) with each subject observed a few times. On homoscedastic data, with normal errors, it rejects too often, and the more often the more levels the factor has and the fewer observations each level has. The excess comes from the dummies of the factor in the auxiliary regression. The residuals of a level are taken from its own mean, so under normal errors the squared residuals of a level of T observations vary between levels T / (T - 1) times as much as the auxiliary regression assumes, and with G levels the statistic exceeds its degrees of freedom by about (G - 1) / (T - 1), as Koenker's did in the study below. Regressed on x alone, without the dummies, the squared residuals were rejected by this test in 4.3% to 5.6% of the draws (inst/validation/fixed-effects-dummies.R). The excess grows with G at a fixed T, while the spread of the chi-squared distribution grows as \sqrt{2G}, so with a fixed number of observations per level the size goes to 1 as the number of levels grows.

Since 0.16.0 the test warns when a term of the model that holds a factor or a character variable has more levels than observations per level: L > n / L, that is L^2 > n, with L the number of distinct combinations of the factors of the term among the n observations of the fit. The term is a factor, such as subject; an interaction of factors, counted by its cells, such as a:b of a crossed a * b or school:class of a nested school / class; or a slope for each level, x:subject. A variable left out of the formula, as in y ~ . - id, does not count. The warning names the term and gives L, n and the mean number of observations per level; the statistic, the degrees of freedom and the p-value are those the test computes without it. It has class heteroTests_fixed_effects_warning, which suppressWarnings(expr, classes = "heteroTests_fixed_effects_warning") muffles, and no other warning. It is given on every call, also when cachedTest() or runHeteroTests() answers from its cache, and runHeteroTests() and runDiagnostics() give it once for each test that gives it. performBPTestStreaming(), performBPTestRobust() with studentize = FALSE, the deprecated performBreuschPaganTest() and the registered "breusch_pagan" give it too.

On a balanced design of G subjects of T observations the rule is G \le T. In the study it kept the three tests within 2.5% to 7.5% at the 5% level where it holds, and put at least one of them outside where it does not, in 68 of the 69 designs where they computed; the exception is this test on 20 subjects of 20, at 8.0%. Unbalanced groups, a second regressor and a regressor correlated with the subject effects followed the same pattern, with the mean number of observations per level, n / G, in place of T. The rule is a warning, not a correction, and many panels fail it: on 100 firms observed over 10 years the three tests reject in 18% to 21% of the draws. Rejection at the 5% level, in percent, in 2,000 homoscedastic draws of lm(y ~ x + subject), with x, the subject effects and the errors standard normal:

A slope for each level counts too, but there the error goes the other way: with a slope for each subject, lm(y ~ x + x:subject), on 30 subjects of 3 and of 4 and on 100 of 8, the three tests rejected too rarely, White's in 0.7% to 2.2% of the draws, Koenker's in 0.1% to 0.8% and the Breusch-Pagan test in 0.4% to 2.8% (inst/validation/fixed-effects-variants-size.R). The warning therefore says that the p-value is unreliable: in the simulations usually far too small, and too large with a slope per level.

The study is inst/validation/fixed-effects-size.R, with its results in fixed-effects-size.csv (system.file("validation", package = "heteroTests")). Juhl and Sosa-Escudero (2014) derive a test of heteroskedasticity for this setting, after fixed effects estimation of a panel of many units observed a fixed number of times; it needs to be told which factor is the unit, and the package does not implement it.

Deprecated name

performBreuschPaganTest() is deprecated since 0.14.0. It was the same function as performBPTest() under a second name; it still works, returns what performBPTest() returns and warns. Use performBPTest().

Weighted fits

On a fit with weights the test is computed from the Pearson residuals \sqrt{w_i}\, e_i, so it tests whether the weights are adequate; see the section on weighted fits in performKoenkerTest(), which also covers weights estimated by fitWLS().

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. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Breusch, T. S., & Pagan, A. R. (1979). A simple test for heteroscedasticity and random coefficient variation. Econometrica, 47(5), 1287–1294. doi:10.2307/1911963

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

Juhl, T., & Sosa-Escudero, W. (2014). Testing for heteroskedasticity in fixed effects models. Journal of Econometrics, 178(Part 3), 484–494. doi:10.1016/j.jeconom.2013.07.005

Examples

 data(mtcars)
 m <- lm(mpg ~ wt + qsec, data = mtcars)
 performBPTest(m, mtcars)

Robust Breusch–Pagan test with bootstrap and effect sizes

Description

Enhances performBPTest() by optionally studentising residuals, resampling the test statistic via bootstrap, and reporting effect sizes together with power diagnostics.

Usage

performBPTestRobust(
  model,
  data = NULL,
  studentize = TRUE,
  bootstrap = FALSE,
  B = 999,
  ci_level = 0.95,
  parallel = FALSE,
  studentized
)

Arguments

Details

Offers expanded reporting for the Breusch–Pagan family of tests including bootstrap p-values and asymptotic confidence intervals. The optional studentize argument toggles between the classic and Koenker versions.

Value

An augmented htest object containing a robust_details element describing the additional diagnostics. Its studentized field records the value of studentize.

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. It accounts for an offset with studentize = FALSE, and since 0.16.0 with studentize = TRUE, the default, whose refit of the design carries the offset of the fit. It refuses a multivariate lm ("needs a model with one response"), and, since 0.16.0, an mgcv::gam() fit with smooth terms: a smooth term is usually penalised, and a penalised fit is not the least-squares fit of the design the test fits again. A later step settles which tests take such a fit; a gam without smooth terms is tested as the Gaussian glm it is. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

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.

See Also

performBPTest() for the baseline test and performKoenkerTest() for the standalone studentised variant.

Examples

data(mtcars)
mod <- lm(mpg ~ wt + qsec, data = mtcars)
performBPTestRobust(mod, mtcars, bootstrap = TRUE, B = 200)

Streaming Breusch-Pagan test for large datasets

Description

Computes the classical Breusch-Pagan statistic using chunked cross-products so that large datasets can be evaluated without materialising the full auxiliary regression in memory. The streamed result is algebraically identical to performBPTest.

Usage

performBPTestStreaming(model, data = NULL, chunk_size = 10000, progress = interactive())

Arguments

Details

Each chunk contributes to the cross-product matrices X'X and X'y for the auxiliary regression of squared residuals on the original regressors. The chunks are aggregated to recover the exact Breusch-Pagan statistic while keeping memory usage bounded. When Matrix is installed, sparse cross-products are used automatically for large chunks.

Since 0.16.0 the chunks are rows of the design of the fit, as in performWhiteTestStreaming(), whose details say what changed.

Value

A htest object mirroring performBPTest and reporting the chi-squared statistic, degrees of freedom, p-value, and metadata describing the chunked computation.

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, data that lack rows or variables of the fit, and, since 0.16.0, a multivariate lm, a fit changed after it was fitted and a fit with model = FALSE whose data are gone, unless data gives its response and fitted values. Which model each test accepts, and why, is set out in test_results, under Supported models.

See Also

performBPTest for the standard implementation.

Examples

data(mtcars)
mod <- lm(mpg ~ wt + qsec, data = mtcars)
performBPTestStreaming(mod, mtcars, chunk_size = 16, progress = FALSE)

Perform Bartlett's test for equality of variances

Description

Uses bartlett.test on model residuals grouped by a factor.

Usage

performBartlettTest(model, data = NULL, group)

Arguments

Details

The test statistic compares the pooled variance to the individual group variances. It is computed as \chi^2 = (N - k) \ln S_p^2 - \sum (n_i - 1) \ln s_i^2, where S_p^2 is the pooled variance and s_i^2 the group variances. Under the null it approximates a chi-square distribution with k-1 degrees of freedom.

Value

An object of class htest containing the chi-squared statistic, p-value and degrees of freedom.

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. It refuses a multivariate lm ("needs a model with one response"), since 0.16.0; before, it stopped with a coercion error. An mgcv::gam() fit is tested as before 0.16.0, with smooth terms too, until a later step settles which tests take such a fit; since 0.16.0 its observations are matched to the rows of data by the row names of its model frame. Before, they were matched by position, since the residuals of a gam carry no names, and data with the rows of the fit in another order gave a wrong value. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Bartlett, M. S. (1937). Properties of sufficiency and statistical tests. Proceedings of the Royal Society of London, 160(901), 268–282. doi:10.1098/rspa.1937.0109

Hartley, H. O. (1950). The maximum F-ratio as a short-cut test for heterogeneity of variance. Biometrika, 37(3/4), 308–312. doi:10.2307/2332383

Examples

 data(mtcars)
 mtcars$cyl <- factor(mtcars$cyl)
 m <- lm(mpg ~ wt, data = mtcars)
 performBartlettTest(m, mtcars, "cyl")

Box's M Test for Equality of Covariance Matrices

Description

Test whether multiple groups have equal covariance matrices, by Box's (1949) chi-square approximation to his M statistic.

Usage

performBoxMTest(data, group)

Arguments

Details

performBoxMTest(data, group) is the one test that does not take (model, data, ...): it tests the covariance matrices of the data themselves, not the residuals of a fitted model.

With g groups of n_i observations of p variables, \nu_i = n_i - 1, \nu = \sum_i \nu_i, S_i the covariance matrix of group i and S the pooled one, Box's statistic is

M = \nu \log|S| - \sum_i \nu_i \log|S_i|,

his correction factor is

c_1 = \left(\sum_i \frac{1}{\nu_i} - \frac{1}{\nu}\right) \frac{2p^2 + 3p - 1}{6 (p + 1) (g - 1)},

and M (1 - c_1) is referred to a chi-square distribution with f_1 = (g - 1) p (p + 1) / 2 degrees of freedom: Box's (1949) chi-square approximation. The test assumes that each group is multivariate normal. It gives the values of biotools::boxM(), heplots::boxM() and rstatix::box_m(): 140.9430 on 20 degrees of freedom for the example on iris.

Before 0.16.0 the correction factor was (g - 1) c_1: the same for two groups, whose values do not change, and too large for three or more, which made the statistic too small and the test conservative. The example on iris gave 135.2229.

The groups are the levels of group that have observations: a level with none, as subset() leaves on a factor, is dropped (since 0.16.0; before, it made the result NA). There must be at least two of them, and each group needs more observations than there are variables, n_i > p: a group of n_i \le p observations has a singular covariance matrix (a group of one has none), whose log-determinant is -\infty, and the statistic is not defined, whatever the data. Since 0.16.0 the test refuses such data with an error that names the groups; before, rounding made the statistic a finite number (133.74, p = 7.7e-19, for a group of three observations of four variables beside two species of iris), a group of one gave NA and a single group NaN. With more observations than variables a group's covariance matrix is still singular when a variable is constant in the group, or a linear combination of the others in it; the test refuses that too, naming the group and saying which, where it returned Inf (p = 0) or NaN with warnings. A group is checked by the rank of its centred data at the tolerance of qr(). A missing value in data or in group still gives NA.

Value

An object of class htest: statistic, named M, is Box's corrected statistic M (1 - c_1); parameter, named df, its degrees of freedom f_1; p.value the upper tail of the chi-square distribution; and method and data.name.

Size

In inst/validation/box-m-approximations-size.R, at 5000 replications per cell and the 5% level, with g = 2, 3 or 5 groups of 10 to 100 observations of p = 1 to 10 variables, the test held its level under Gaussian data in 60 of 67 cells, inside the release gate's interval, 0.042 to 0.058.

• It is liberal when the groups are small for the number of variables. With groups of 10 and 4 or 6 variables it rejected a true null hypothesis 0.055 to 0.106 of the time, outside the gate in four of those six cells (0.062 to 0.106, the most with five groups of 10 and six variables), and with groups of 20 and 10 variables 0.071 to 0.087. Box's F approximation, which the package does not compute, was closer to the level in each of those nine cells (0.051 to 0.070); in the other Gaussian cells the two differed by at most 0.005.

• It is not robust to heavy tails. With independent t_5 coordinates, three groups of 20, 50 or 100 and two or four variables, it rejected 0.26 to 0.48 of the time (0.18 to 0.46 with the factor of 0.15.0), more with larger groups: the test then detects the tails as well as unequal covariance matrices.

• It is conservative with light tails. On the data of the full sweep (inst/validation/full-sweep-size-power.R), 150 observations of (y, x) with x uniform and grouped by cut(x, 3), it rejected 0.017 of the time, as Box's F did. The sweep's 0.012, at 400 replications, was measured with the factor of 0.15.0, which gives 0.016 in that design at 5000 replications: the low size comes from the uniform, light-tailed data, not from the approximation.

When the groups are small for the number of variables, or the data may have heavy tails, a rejection is weak evidence of unequal covariance matrices.

Validation

tests/testthat/test-boxm.R checks the statistic against Box's approximation written out for two to five groups and one to four variables, and against the values of biotools::boxM() 4.3 and heplots::boxM() 1.8.6 on iris and on mtcars, and that two groups give the values of 0.15.0. inst/validation/box-m-approximations-size.R compares the chi-square and the F approximation with biotools, heplots and DFA.CANCOR on ten data sets, and measures the sizes above.

References

Box, G. E. P. (1949). A general distribution theory for a class of likelihood criteria. Biometrika, 36(3/4), 317–346. doi:10.1093/biomet/36.3-4.317

Examples

data(iris)
performBoxMTest(iris[,1:4], iris$Species)

# Two species: with two groups the values are those of 0.15.0
two <- droplevels(subset(iris, Species != "virginica"))
performBoxMTest(two[, 1:4], two$Species)

Perform Brown-Forsythe test for equality of variances

Description

Brown-Forsythe test using medians instead of means.

Usage

performBrownForsytheTest(model, data = NULL, group)

Arguments

Details

The absolute deviations from the group medians are compared across groups. The resulting statistic is an F ratio with k-1 and N-k degrees of freedom.

Value

An object of class htest containing the F statistic, p-value and degrees of freedom.

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. It refuses a multivariate lm ("needs a model with one response"), since 0.16.0; before, it stopped with a coercion error. An mgcv::gam() fit is tested as before 0.16.0, with smooth terms too, until a later step settles which tests take such a fit; since 0.16.0 its observations are matched to the rows of data by the row names of its model frame. Before, they were matched by position, since the residuals of a gam carry no names, and data with the rows of the fit in another order gave a wrong value. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Brown, M. B., & Forsythe, A. B. (1974). Robust tests for the equality of variances. Journal of the American Statistical Association, 69(346), 364–367.

Examples

 data(mtcars)
 mtcars$cyl <- factor(mtcars$cyl)
 m <- lm(mpg ~ wt, data = mtcars)
 performBrownForsytheTest(m, mtcars, "cyl")

Perform Cook-Weisberg test for heteroscedasticity

Description

The Cook-Weisberg (1983) score test, with the fitted values as the sole variance regressor by default, or against a variance model given by var_formula. Since 0.14.0 it is the one entry point for the test; performNCVTest is deprecated.

Usage

performCookWeisbergTest(model, data = NULL, var_formula = NULL)

Arguments

Details

Writing \hat{\sigma}^2 = \sum_i \hat{e}_i^2 / n for the maximum-likelihood error variance, the test regresses the scaled squared residuals \hat{e}_i^2 / \hat{\sigma}^2 on the fitted values and refers half the explained sum of squares to a chi-square distribution with one degree of freedom. This is the diagnostic returned by Stata's estat hettest with the fitted option, and it is exactly performNCVTest with its default variance model.

The chi-square reference distribution follows from the null variance of \hat{e}^2 / \sigma^2 being 2, which holds under normal errors. When that assumption is doubtful, prefer performKoenkerTest.

Value

An object of class htest containing the score statistic, its single degree of freedom and the p-value.

Validation

Reproduces car::ncvTest() to within 1e-8; see ‘tests/testthat/test-pass-a-reference.R’. Releases before 0.7.0 returned n R^2 from regressing the raw squared residuals on the fitted values, which is the studentized (Koenker) statistic rather than the Cook-Weisberg score test.

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.16.0 it is built on the model context, and also refuses a multivariate lm ("needs a model with one response"), a fit whose coefficients, residuals or weights were changed after the fit, and a fit with model = FALSE whose model frame, rebuilt from the data it was fitted to, no longer holds its response, because those data changed after the fit, with or without data. A fit with model = FALSE whose data are gone is tested on data that give its response and fitted values, and refused without them. Since 0.16.0 data is matched to the observations of the fit as by every test, with or without var_formula: without it the test reads nothing from data, and gave the value of the fit on numbered data in another order, rows in front of those of the fit and repeated row names, which are refused. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Cook, R. D., & Weisberg, S. (1983). Diagnostics for heteroscedasticity in regression. Biometrika, 70(1), 1–10. doi:10.1093/biomet/70.1.1

See Also

performNCVTest, performKoenkerTest, performBPTest

Examples

 data(mtcars)
 m <- lm(mpg ~ wt + qsec, data = mtcars)
 performCookWeisbergTest(m)

 # Score test against a specific variance model rather than the fitted values
 performCookWeisbergTest(m, var_formula = ~wt)

 # With the data, the variance model can use variables the mean model does not
 performCookWeisbergTest(m, mtcars, var_formula = ~ wt + hp)

Perform Davidian-Carroll test

Description

Regresses log(residual^2) on fitted values using a polynomial.

Usage

performDavidianCarrollTest(model, data = NULL, degree = 2)

Arguments

Details

A polynomial in the fitted values is fitted to \log(e^2). Joint significance of the polynomial terms is evaluated with an F statistic.

Value

An object of class htest.

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.16.0 it is built on the model context, and also refuses a multivariate lm ("needs a model with one response"), a fit whose coefficients, residuals or weights were changed after the fit, and a fit with model = FALSE whose model frame, rebuilt from the data it was fitted to, no longer holds its response, because those data changed after the fit, with or without data. A fit with model = FALSE whose data are gone is tested on data that give its response and fitted values, and refused without them. The test reads nothing from data, and since 0.16.0 matches it to the observations of the fit as every test does: data with row names of their own in any order, data that number their rows (a tibble, a data.table) in the order of the fit; other numbered data, rows in front of those of the fit and repeated row names are refused, where it gave the value of the fit. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Davidian, M., & Carroll, R. J. (1987). Variance function estimation. Journal of the American Statistical Association, 82(400), 1079–1091.

Examples

 data(mtcars)
 m <- lm(mpg ~ wt + qsec, data = mtcars)
 performDavidianCarrollTest(m)

Perform Fligner-Killeen test for homogeneity of variances

Description

The median-centred Fligner-Killeen test of homogeneity of variances, on the Pearson residuals of the fit, computed with fligner.test(): it ranks the absolute deviations of the residuals from their group medians and scores the ranks by quantiles of the normal distribution, which makes it robust to outliers and to errors that are not normal.

Usage

performFlignerKilleenTest(model, data = NULL, group)

Arguments

Details

The test is the median-centred version of the test of Fligner and Killeen (1976), the "F-K:med X^2" of Conover, Johnson and Johnson (1981), whose simulation study found it among the tests of homogeneity of variances most robust to departures from normality; it is the test fligner.test() computes. The residuals of each group are centred at the median of the group, the absolute values of the centred residuals are ranked over all n observations, and the rank r_i is given the score a_i = \Phi^{-1}\{(1 + r_i / (n + 1)) / 2\}. With \bar a_j the mean score of group j, of size n_j, \bar a the mean of all the scores and s_a^2 = \sum_i (a_i - \bar a)^2 / (n - 1), the statistic

X^2 = \sum_j n_j (\bar a_j - \bar a)^2 / s_a^2

is referred to a chi-squared distribution with k-1 degrees of freedom, for k groups. No logarithm is taken. The observations whose group is missing are left out before the statistic is computed.

Tied values are given the mean of their ranks. Since 0.16.0 the deviations are those of the residuals r_i = \sqrt{w_i}\,(y_i - x_i^\top \hat\beta - o_i) computed observation by observation, as the Gaussian glm computes them, and not the residuals an lm stores, which its QR decomposition rounds differently from one observation to the next. For an mgcv::gam() fit x_i is the row of its linear predictor matrix. Observations with the same values of the model's variables, weight and offset, such as repeated values of an integer response or of a response recorded with few digits, are then given the same residual, as in exact arithmetic, where their residuals agree to within rounding, even when a basis such as poly() gives them design rows that differ in their last bits; a term of the row position, such as seq_along(x), keeps them apart. This needs the model's variables: it is done when data is given, or when the model frame holds them, as it holds a variable the model uses under its own name (x, not poly(x, 3)). And absolute deviations closer than 8 \epsilon s_j, with \epsilon = .Machine$double.eps and s_j the largest |r_i| of group j (the smaller of the two for deviations of different groups), are ranked as tied: in a group of even size the two residuals in the middle lie at the same distance from the median of the group, half their difference, a tie that the centring breaks by up to \epsilon s_j. On data with repeated observations the statistic can differ much from before: on lm(extra ~ group, data = sleep) it was 0.154 and is 0.213, the value of exact ties. On continuous data it moves by much less (see NEWS). Where neither changes a rank, the value is that of fligner.test() on the residuals of the fit, as before. Two limits remain. A tie between different observations that exact decimal arithmetic would give, as when two pairs of observations of a regressor recorded with few digits differ by the same amounts, is restored only when rounding leaves it within the tolerance, and two forms of one model can then differ by several percent: on 70 observations that repeat 30 values of x recorded with one decimal, y ~ poly(x, 3) + g gives 1.21547, and y ~ x + I(x^2) + I(x^3) + g and y ~ poly(x, 3, raw = TRUE) + g give 1.30280, since poly() breaks in its last bits the ties between repeated values of x in different groups. And in a group whose residuals are far from zero, 10^8 times their spread, the tolerance may also tie two deviations that differ by a few units in the last place of the residuals.

Value

An object of class htest containing the chi-squared statistic, p-value and degrees of freedom.

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. It refuses a multivariate lm ("needs a model with one response"), since 0.16.0; before, it stopped with a coercion error. An mgcv::gam() fit is tested as before 0.16.0, with smooth terms too, until a later step settles which tests take such a fit; since 0.16.0 its observations are matched to the rows of data by the row names of its model frame. Before, they were matched by position, since the residuals of a gam carry no names, and data with the rows of the fit in another order gave a wrong value. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Conover, W. J., Johnson, M. E., & Johnson, M. M. (1981). A comparative study of tests for homogeneity of variances, with applications to the outer continental shelf bidding data. Technometrics, 23(4), 351–361. doi:10.1080/00401706.1981.10487680

Fligner, M. A., & Killeen, T. J. (1976). Distribution-free two-sample tests for scale. Journal of the American Statistical Association, 71(353), 210–213. doi:10.1080/01621459.1976.10481517

Examples

 data(mtcars)
 mtcars$cyl <- factor(mtcars$cyl)
 m <- lm(mpg ~ wt, data = mtcars)
 performFlignerKilleenTest(m, mtcars, "cyl")

Perform Goldfeld-Quandt test for heteroscedasticity

Description

Implements the Goldfeld-Quandt test on a fitted linear model with directional or two-sided alternatives.

Usage

performGQTest(model, data = NULL, order_by, fraction = 0.2,
  alternative = c("greater", "two.sided", "less"))

Arguments

Details

Observations are ordered by a suspected variance-driving variable. With the split point fixed at the sample midpoint, a central fraction is omitted and the model is fitted again on the lower and upper segments: each segment is fitted to its rows of the design, the response and the offset of the fit, as lmtest::gqtest() fits them to the rows of the design. The statistic is the residual mean square of segment 2 divided by the residual mean square of segment 1 and is therefore directional rather than being forced above one. Split arithmetic and p-value conventions match lmtest::gqtest(..., point = 0.5).

Before 0.16.0 the segments of an unweighted fit were refitted from the formula, which computed again, on the rows of each segment, a term that depends on the rows it is computed from: ns(x) and bs(x) placed their knots at quantiles of the segment, and I(x > median(x)) took its median, so the segments were fitted to columns the fit never used; poly(x) computed another basis of the same columns, which changed the last bits only. The refit also left out an offset given as an argument, and a weighted fit was refitted without its offset.

Value

An object of class htest containing the directional GQ statistic, p-value, degrees of freedom and alternative hypothesis.

Weighted fits

On a fit with weights each segment is refitted with its own weights and the segments are compared on the weighted residual sum of squares, \sum w_i e_i^2. That is the Goldfeld-Quandt statistic of the equivalent unweighted regression, so the test asks whether the weights are adequate. Releases before 0.12.0 refitted the segments without the weights.

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. It accepts a model with an offset, given as an argument or as an offset() term, weighted or not, and gives the value of the model fitted to the response minus the offset (since 0.16.0). It refuses a multivariate lm ("needs a model with one response"), and, since 0.16.0, an mgcv::gam() fit with smooth terms: a smooth term is usually penalised, and a penalised fit is not the least-squares fit of the design the test fits again. A later step settles which tests take such a fit; a gam without smooth terms is tested as the Gaussian glm it is. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Goldfeld, S. M., & Quandt, R. E. (1965). Some tests for homoscedasticity. Journal of the American Statistical Association, 60(310), 539–547. doi:10.1080/01621459.1965.10480811

Greene, W. H. (2018). Econometric Analysis (8th ed.). Pearson.

Examples

 data(mtcars)
 m <- lm(mpg ~ wt + qsec, data = mtcars)
 performGQTest(m, mtcars, order_by = "wt")
 performGQTest(m, mtcars, order_by = "wt", alternative = "two.sided")

Perform Glejser test for heteroscedasticity

Description

Regresses absolute residuals on a transformation of a suspected variable. By default, robust = TRUE since 0.15.0, the absolute residuals are first corrected for the asymmetry of the errors: this is the version of Im (2000) and of Machado and Santos Silva (2000), which keeps its level when the errors are asymmetric, from about 150 observations on (see the section on asymmetric errors for smaller samples). robust = FALSE gives Glejser's (1969) statistic, the default up to 0.14.0.

Usage

performGlejserTest(model, data = NULL, variable,
                   transformation = c("abs", "sqrt", "inverse", "inverse_sqrt"),
                   robust = TRUE)

Arguments

Details

A variety of transformations of the explanatory variable (e.g. absolute value, square root) can reveal a relationship between scale and the covariate. Significance of the slope parameter is assessed with a t test. For a weighted fit the residuals are the Pearson residuals.

Glejser's statistic, robust = FALSE, has its nominal level only when positive and negative errors are equally likely. The corrected statistic of the default does not need that, and is the one to use when the residuals may be skewed or their shape is unknown.

Either form conditions on one variable and one transformation. Scanning several transformations and reporting only the smallest p-value invalidates the nominal level. Simulated size and power of both statistics are recorded in inst/validation/glejser-skewed-errors.csv; those of Glejser's statistic in inst/validation/pass-a-size-power.csv, and those of the default in inst/validation/full-sweep-size-power.csv and inst/validation/weighted-fits-size.csv.

Value

An object of class htest containing the t statistic, p-value and degrees of freedom. The method element says which of the two statistics was computed.

Asymmetric errors

Glejser's statistic takes the absolute residuals for the absolute errors. The two differ by the estimation error of the mean equation,

|\hat{e}_i| \approx |e_i| - \mathrm{sign}(e_i)\, x_i^\top (\hat{\beta} - \beta),

and the second term averages out only when positive and negative errors are equally likely (Godfrey 1996). When they are not, and the variable tested is correlated with the regressors, the test does not in general have its nominal level, and a larger sample does not help. In the designs of inst/validation/glejser-skewed-errors.R Glejser's statistic, robust = FALSE, rejects a true null hypothesis about 13% of the time at the 5% level under exponential errors and about 17% under lognormal errors, at every sample size from 50 to 1000. For other error distributions the rate can also fall below the nominal level.

robust = TRUE replaces the regressand |\hat{e}_i| by

|\hat{e}_i| - \hat{m}\, \hat{e}_i, \qquad \hat{m} = \frac{1}{n} \sum_{i=1}^{n} \mathrm{sign}(\hat{e}_i),

where \hat{m} is the proportion of positive residuals less the proportion of negative ones. This is the correction of Im (2000). The regressand of Machado and Santos Silva (2000), \hat{e}_i\,[I(\hat{e}_i \ge 0) - \hat{\eta}] with \hat{\eta} the proportion of non-negative residuals, is half of it when no residual is exactly zero, and the two then give the same test. A residual that is zero up to rounding error is counted as zero, so that the statistic does not depend on the sign of a rounding error.

The auxiliary regression is otherwise the same. Both references state the statistic as n R^2 of that regression, referred to a chi-squared distribution with one degree of freedom. The t statistic reported here is the signed form of it, n R^2 = n t^2 / (t^2 + n - 2), and its p-value is taken from the t distribution with n - 2 degrees of freedom, as for Glejser's statistic. The two p-values agree as the sample grows.

In the same study, with 5000 replications in each cell, the corrected statistic is within simulation error of the 5% level from 150 observations on: all 27 of its cells with a true null hypothesis and 150 observations or more are inside the 99% interval around 5%, from 4.21% to 5.79%, at 4.4% to 5.7%. At 50 observations it rejects 5.7% to 6.6% of the time under skewed errors, and four of those six cells are above that interval, at 5.80% to 6.58%. Under symmetric errors \hat{m} is close to zero and the two statistics differ little: the corrected one rejects 0.5 to 0.7 percentage points less often at 50 observations and gives up one or two points of power there. For a weighted fit the correction is applied to the Pearson residuals.

The default

Up to 0.14.0 the default was robust = FALSE. Since 0.15.0 it is robust = TRUE, so a call that does not give robust returns a different statistic and p-value from those it returned before, and its method ends in "(robust to asymmetric errors)". robust = FALSE returns what such a call returned up to 0.14.0.

Glejser's statistic needs positive and negative errors to be equally likely, which is seldom known: in the study above it rejects a true null hypothesis 8.8% to 17.7% of the time at the 5% level when the errors are skewed, at every sample size from 50 to 1000, where the corrected statistic rejects 4.6% to 5.7% of the time from 150 observations on. The change costs some power when the errors are symmetric, 1.1 and 1.9 percentage points at 50 observations under Gaussian errors in that study, and the corrected statistic is above the 5% level at 50 observations under skewed errors, as above.

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. It refuses a multivariate lm ("needs a model with one response"), since 0.16.0; before, it stopped with a coercion error. An mgcv::gam() fit is tested as before 0.16.0, with smooth terms too, until a later step settles which tests take such a fit; since 0.16.0 its observations are matched to the rows of data by the row names of its model frame. Before, they were matched by position, since the residuals of a gam carry no names, and data with the rows of the fit in another order gave a wrong value. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Glejser, H. (1969). A new test for heteroskedasticity. Journal of the American Statistical Association, 64(325), 316–323.

Godfrey, L. G. (1996). Some results on the Glejser and Koenker tests for heteroskedasticity. Journal of Econometrics, 72(1–2), 275–299. doi:10.1016/0304-4076(94)01722-0

Im, K. S. (2000). Robustifying Glejser test of heteroskedasticity. Journal of Econometrics, 97(1), 179–188. doi:10.1016/S0304-4076(99)00061-5

Machado, J. A. F. and Santos Silva, J. M. C. (2000). Glejser's test revisited. Journal of Econometrics, 97(1), 189–202. doi:10.1016/S0304-4076(00)00016-6

See Also

performKoenkerTest for a test on the squared residuals that needs neither normal nor symmetric errors; performParkTest and performHarveyTest for related parametric tests.

Examples

 data(mtcars)
 m <- lm(mpg ~ wt + qsec, data = mtcars)
 # The default, which keeps its level under asymmetric errors from about
 # 150 observations
 performGlejserTest(m, mtcars, "wt")

 # Glejser's (1969) statistic, the default up to 0.14.0
 performGlejserTest(m, mtcars, "wt", robust = FALSE)

 # Examine an alternative transformation
 performGlejserTest(m, mtcars, "wt", transformation = "inverse")

Harrison-McCabe test for heteroscedasticity

Description

Orders the observations by a variable and tests whether the share of the residual sum of squares that falls in the first part of the ordered sample is too small, or too large, for errors of constant variance. The p-value is exact under normal errors, given the regressors.

Usage

performHarrisonMcCabeTest(
  model,
  data = NULL,
  order_by,
  point = 0.5,
  alternative = c("less", "two.sided", "greater")
)

Arguments

Details

With \hat{e}_{(1)}, \dots, \hat{e}_{(n)} the least-squares residuals in the order of order_by, the statistic of Harrison and McCabe (1979) is

b = \frac{\sum_{i=1}^{m} \hat{e}_{(i)}^2}{\sum_{i=1}^{n} \hat{e}_{(i)}^2},

the share of the residual sum of squares that falls in the first m observations. A variance that increases along the ordering makes b small, and one that decreases makes it large. It is the statistic of lmtest::hmctest().

The 1979 article was not consulted for this implementation. The definition follows Harrison (1981), by its first author, which restates the statistic with m the integer part of n / 2 and rejects for small b, and it agrees with lmtest::hmctest().

Under normal errors of constant variance, b is a ratio of quadratic forms in normal variables, and its distribution depends on the regressors. With M the residual-maker of the design in the ordered sample and D the diagonal matrix that selects the first m observations,

P(b \le b_0) = P\left(\sum_j (\lambda_j - b_0) z_j^2 \le 0\right),

where the z_j are independent standard normal and the \lambda_j are the n - k eigenvalues of M D M on the residual space, k being the number of coefficients. They are found from the k singular values of the first m rows of an orthonormal basis of the design, and the probability is evaluated by Imhof's (1961) formula with stats::integrate(), as Harrison (1981) does for the exact variant of the test. For a constant and a linear trend, the 5% critical values this gives at 10, 20, 30 and 40 observations are those of Harrison's (1981) Table 3 to the three decimals printed there. Under the null hypothesis the mean of b is (m - \sum_{i \le m} h_{(i)}) / (n - k), with h_{(i)} the leverages in the ordered sample.

The probability is computed to an absolute accuracy of about 2e-11. Smaller p-values are not resolved, and a statistic far in a tail can return a p-value of 0.

The p-value depends only on the space spanned by the regressors, so it does not change with the parametrisation of the model.

Value

An object of class htest with the share HMC as the statistic, the number of observations in the first part m and the sample size n as parameters, the p-value, and the alternative in words.

The rank

k is the rank of the fit, the number of its coefficients that are not aliased, decided at the tolerance the fit used: 1e-7 for an lm, or the tol it was given, and 1e-11 for a glm (since 0.16.0). Before, the rank was decided again at 1e-7, so a Gaussian glm, or an lm fitted with a smaller tol, whose design holds a column that is a combination of the others to within 1e-7 but not to within its own tolerance, was given a null distribution of a smaller rank than its residuals. With x3 = x1 + 1e-9 * rnorm(120) the glm has rank 4, and its p-value is 0.002273 where it was 0.002288, computed with rank 3. On every other fit the two ranks agree and the value is unchanged.

Direction of the alternative

In performGQTest() and performSzroeterTest() the statistic grows when the variance increases along the ordering, and in this test it shrinks. The alternative names the tail of each statistic, so the same hypothesis, a variance that increases with order_by, is alternative = "less" here and alternative = "greater" in those two functions. In all three it is the default, and the alternative element of the result states it in words.

Weighted fits

On a fit with weights the statistic is computed from the Pearson residuals \sqrt{w_i}\, \hat{e}_i, and the null distribution from the weighted design \sqrt{w_i}\, x_i. Both are those of the equivalent unweighted regression, so the test asks whether the weights are adequate. The weights of a Gaussian glm are its prior weights.

Heavy tails

The null distribution assumes normal errors. In the designs of inst/validation/harrison-mccabe-size-power.R the test rejects a true null hypothesis 11.7% to 16.3% of the time at the 5% level when the errors follow a t_5 distribution, and more often the larger the sample; performGQTest() and performSzroeterTest() do the same. When the tails may be heavy, use performGlejserTest() with the ordering variable as variable and robust = TRUE, which in the same study rejects 4.7% to 4.9% of the time under those errors. Its default transformation takes the absolute value of the variable, which preserves the ordering only when the variable is positive.

Reference implementations

The statistic agrees with lmtest::hmctest() to 1e-8. Its p-value does not: with simulate.p = TRUE it is simulated from residuals of a regression on a constant alone, whatever the model, and with simulate.p = FALSE it is NA. skedastic::harrison_mccabe() reports a different p-value as well; inst/validation/README.md says why.

Supported models

It accepts an lm and, since 0.16.0, 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. The Gaussian glm with prior weights w is the least-squares fit lm(..., weights = w): it has the same design, residuals and weights, and it states the normal errors of constant variance, given the weights, under which the null distribution above is exact, so the test gives the value of the lm when the two decide the rank alike (see the section on the rank). It refused every glm before 0.16.0. 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. It refuses a multivariate lm ("needs a model with one response"), since 0.16.0; before, it stopped with a coercion error. An mgcv::gam() fit is tested as before 0.16.0, with smooth terms too, until a later step settles which tests take such a fit; since 0.16.0 its observations are matched to the rows of data by the row names of its model frame. Before, they were matched by position, since the residuals of a gam carry no names, and data with the rows of the fit in another order gave a wrong value. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Harrison, M. J. (1981). A comparison of the bounds, beta-approximate, and exact variants of two tests for heteroscedasticity based on ordinary least squares residuals. The Economic and Social Review, 12(4), 235–252.

Harrison, M. J., & McCabe, B. P. M. (1979). A test for heteroscedasticity based on ordinary least squares residuals. Journal of the American Statistical Association, 74(366a), 494–499. doi:10.1080/01621459.1979.10482544

Imhof, J. P. (1961). Computing the distribution of quadratic forms in normal variables. Biometrika, 48(3-4), 419–426. doi:10.1093/biomet/48.3-4.419

See Also

performGQTest() and performSzroeterTest() for the other tests along an ordering of the observations; performGlejserTest() when the errors may be heavy-tailed.

Examples

data(mtcars)
mod <- lm(mpg ~ wt + qsec, data = mtcars)
performHarrisonMcCabeTest(mod, mtcars, order_by = "wt")

# Variance that increases with x
set.seed(7)
x <- runif(80, 1, 5)
d <- data.frame(x = x, y = 1 + 2 * x + rnorm(80, sd = exp(0.3 * x)))
performHarrisonMcCabeTest(lm(y ~ x, data = d), d, order_by = "x")

# Either direction, with a third of the sample in the first part
performHarrisonMcCabeTest(lm(y ~ x, data = d), d,
  order_by = "x",
  point = 1 / 3, alternative = "two.sided"
)

Perform Hartley's Fmax test

Description

Compares the maximum and minimum group residual variances using Hartley's maximum F-ratio distribution.

Usage

performHartleyFmaxTest(model, data = NULL, group)

Arguments

Details

The statistic is F_{max}=\max_j s_j^2/\min_j s_j^2. Under normality and equal group sizes its null distribution is the maximum F-ratio distribution for k independent mean squares with common degrees of freedom, evaluated with SuppDists::pmaxFratio(). If group sizes differ, the function warns and rounds the mean group size, then subtracts one, to obtain an approximate common integer degrees of freedom.

Value

An object of class htest containing the Fmax statistic, number of groups, reference degrees of freedom and p-value.

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. It refuses a multivariate lm ("needs a model with one response"), since 0.16.0; before, it stopped with a coercion error. An mgcv::gam() fit is tested as before 0.16.0, with smooth terms too, until a later step settles which tests take such a fit; since 0.16.0 its observations are matched to the rows of data by the row names of its model frame. Before, they were matched by position, since the residuals of a gam carry no names, and data with the rows of the fit in another order gave a wrong value. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Hartley, H. O. (1950). The maximum F-ratio as a short-cut test for heterogeneity of variance. Biometrika, 37(3/4), 308–312.

Examples

 set.seed(1701)
 n <- 20
 d <- data.frame(g = factor(rep(letters[1:3], each = n)), x = rnorm(3 * n))
 d$y <- 1 + d$x + rnorm(3 * n)
 m <- lm(y ~ x, data = d)
 performHartleyFmaxTest(m, d, "g")

Perform Harvey test for multiplicative heteroscedasticity

Description

Harvey's (1976) Lagrange multiplier test, which regresses log(residuals^2) on a set of variance regressors.

Usage

performHarveyTest(model, data = NULL, auxiliary = c("regressors", "fitted"),
                  studentize = FALSE)

Arguments

Details

Harvey (1976) models the error variance as \sigma_i^2 = \exp(z_i^\top \gamma) and tests \gamma = 0 through the auxiliary regression of \log \hat{e}_i^2 on z_i. Under the null the auxiliary error behaves like a centred \log \chi^2_1 variate with variance \pi^2 / 2 \approx 4.9348, so the classical statistic is \mathrm{ESS} / (\pi^2 / 2), asymptotically chi-square with q degrees of freedom.

Setting studentize = TRUE replaces the fixed constant \pi^2 / 2 with the auxiliary residual mean square and refers the overall F statistic to an F distribution. The two forms are asymptotically equivalent under normal errors; the studentized form is more reliable when normality is doubtful, because \pi^2 / 2 is the null variance of \log \hat{e}^2 only for Gaussian errors.

Degrees of freedom are taken from the realised rank of the auxiliary design, so a rank-deficient variance model is not credited with degrees of freedom for aliased columns.

Value

An object of class htest containing the test statistic, its degrees of freedom and the p-value.

Validation

The default statistic reproduces an independent reconstruction of Harvey (1976) to within 1e-8; see ‘tests/testthat/test-pass-a-reference.R’. Simulated size and power are recorded in ‘inst/validation/pass-a-size-power.csv’.

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.16.0 it is built on the model context, and also refuses a multivariate lm ("needs a model with one response"), a fit whose coefficients, residuals or weights were changed after the fit, and a fit with model = FALSE whose model frame, rebuilt from the data it was fitted to, no longer holds its response, because those data changed after the fit, with or without data. A fit with model = FALSE whose data are gone is tested on data that give its response and fitted values, and refused without them. The test reads nothing from data, and since 0.16.0 matches it to the observations of the fit as every test does: data with row names of their own in any order, data that number their rows (a tibble, a data.table) in the order of the fit; other numbered data, rows in front of those of the fit and repeated row names are refused, where it gave the value of the fit. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Harvey, A. C. (1976). Estimating regression models with multiplicative heteroscedasticity. Econometrica, 44(3), 461–465. doi:10.2307/1913974

Greene, W. H. (2018). Econometric Analysis (8th ed.). Pearson. Section 9.5 derives the \mathrm{ESS} / 4.9348 form of the statistic.

See Also

performParkTest, performGlejserTest, performBPTest

Examples

 data(mtcars)
 m <- lm(mpg ~ wt + qsec, data = mtcars)
 performHarveyTest(m)

 # Studentized variant, which does not assume normal errors
 performHarveyTest(m, studentize = TRUE)

 # Variance driven by the conditional mean rather than by the regressors
 performHarveyTest(m, auxiliary = "fitted")

High-dimensional projection heteroscedasticity test

Description

Applies principal component projections to the design matrix and evaluates a Breusch-Pagan style statistic on the reduced representation. The procedure is designed for settings where the number of predictors rivals or exceeds the sample size.

Usage

performHighDimensionalTest(
  model,
  data = NULL,
  variance_threshold = 0.9,
  max_components = NULL
)

Arguments

Value

An htest object summarising the chi-squared statistic of the auxiliary regression on principal component scores.

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. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Fan, J., & Lv, J. (2008). Sure independence screening for ultrahigh dimensional feature space. Journal of the Royal Statistical Society: Series B, 70(5), 849–911.

Examples

set.seed(123)
X <- matrix(rnorm(80 * 20), nrow = 80)
beta <- c(rep(1, 5), rep(0, 15))
y <- X %*% beta + rnorm(80, sd = 0.5 + 0.2 * scale(X[, 1]))
df <- as.data.frame(cbind(y = as.numeric(y), X))
model <- lm(y ~ ., data = df)
performHighDimensionalTest(model, df)

Identify influential observations

Description

Uses Cook's distance to flag influential observations.

Usage

performInfluenceDiagnostics(model, cutoff = NULL)

Arguments

Value

A list with Cook's distances and indices of influential points.

Examples

data(mtcars)
m <- lm(mpg ~ wt + qsec, data = mtcars)
performInfluenceDiagnostics(m)

Perform Koenker studentized Breusch-Pagan test

Description

Implementation of the Koenker version of the Breusch-Pagan test.

Usage

performKoenkerTest(model, data = NULL)

Arguments

Details

Squared residuals are regressed on the regressors but the statistic n R^2 uses a studentized form that is robust to non-normality.

The test is computed in memory whatever the size of the data. Up to 0.15.0 it switched to performKoenkerTestStreaming() on data above 100 MB; call that function to stream.

Value

An object of class htest containing the test statistic, p-value and degrees of freedom.

Fixed effects for small groups

Koenker's test, like White's and the Breusch-Pagan tests, is not valid with fixed effects for small groups: a factor with many levels and few observations of each, such as subject in lm(y ~ x + subject) with each subject observed a few times. On homoscedastic data it rejects too often, and the more often the more levels the factor has and the fewer observations each level has. The excess comes from the dummies of the factor in the auxiliary regression. The residuals of a level are taken from its own mean, so under normal errors the squared residuals of a level of T observations vary between levels T / (T - 1) times as much as the auxiliary regression assumes, and with G levels the statistic has mean its degrees of freedom plus (G - 1) / (T - 1), as it had in the study below, within 2.3% for T of 3 or more. Regressed on x alone, without the dummies, the squared residuals were rejected in 4.6% to 5.5% of the draws (inst/validation/fixed-effects-dummies.R). The excess grows with G at a fixed T, while the spread of the chi-squared distribution grows as \sqrt{2G}, so with a fixed number of observations per level the size goes to 1 as the number of levels grows.

Since 0.16.0 the test warns when a term of the model that holds a factor or a character variable has more levels than observations per level: L > n / L, that is L^2 > n, with L the number of distinct combinations of the factors of the term among the n observations of the fit. The term is a factor, such as subject; an interaction of factors, counted by its cells, such as a:b of a crossed a * b or school:class of a nested school / class; or a slope for each level, x:subject. A variable left out of the formula, as in y ~ . - id, does not count. The warning names the term and gives L, n and the mean number of observations per level; the statistic, the degrees of freedom and the p-value are those the test computes without it. It has class heteroTests_fixed_effects_warning, which suppressWarnings(expr, classes = "heteroTests_fixed_effects_warning") muffles, and no other warning. It is given on every call, also when cachedTest() or runHeteroTests() answers from its cache, and runHeteroTests() and runDiagnostics() give it once for each test that gives it. performKoenkerTestStreaming(), performBPTestRobust() with studentize = TRUE, the deprecated performStudentizedBPTest() and the registered "koenker" and "student_bp" give it too. performWildBootstrapTest(), which bootstraps the same statistic, held its level on these designs and does not give it.

On a balanced design of G subjects of T observations the rule is G \le T. In the study it kept the three tests within 2.5% to 7.5% at the 5% level where it holds, and put at least one of them outside where it does not, in 68 of the 69 designs where they computed; the exception is the Breusch-Pagan test on 20 subjects of 20, at 8.0%. Unbalanced groups, a second regressor and a regressor correlated with the subject effects followed the same pattern, with the mean number of observations per level, n / G, in place of T. The rule is a warning, not a correction, and many panels fail it: on 100 firms observed over 10 years the three tests reject in 18% to 21% of the draws. Rejection at the 5% level, in percent, in 2,000 homoscedastic draws of lm(y ~ x + subject), with x, the subject effects and the errors standard normal:

A slope for each level counts too, but there the error goes the other way: with a slope for each subject, lm(y ~ x + x:subject), on 30 subjects of 3 and of 4 and on 100 of 8, the three tests rejected too rarely, White's in 0.7% to 2.2% of the draws, Koenker's in 0.1% to 0.8% and the Breusch-Pagan test in 0.4% to 2.8% (inst/validation/fixed-effects-variants-size.R). The warning therefore says that the p-value is unreliable: in the simulations usually far too small, and too large with a slope per level.

The study is inst/validation/fixed-effects-size.R, with its results in fixed-effects-size.csv (system.file("validation", package = "heteroTests")). Juhl and Sosa-Escudero (2014) derive a test of heteroskedasticity for this setting, after fixed effects estimation of a panel of many units observed a fixed number of times; it needs to be told which factor is the unit, and the package does not implement it.

Weighted fits

lm(..., weights = w) states that the error variance is \sigma^2 / w_i. On such a fit the test is computed from the Pearson residuals \sqrt{w_i}\, e_i, the residuals of the equivalent unweighted regression, which have constant variance when the weights are right. The null hypothesis is then that the weights are adequate; the variance regressors stay on their original scale. Releases before 0.12.0 used the raw residuals, whose variance differs across observations by assumption, so a correctly weighted fit was rejected as often as an unweighted one. A fit with a zero weight is refused.

The reference distribution takes the weights as known. When they were estimated from the same data, as by fitWLS(), it is only approximate, and not more accurate in larger samples. With a correctly specified variance function and normal errors, the Koenker test rejected 8% of the time at the 5% level with 150 observations and 11% with 600, the White test 6% and 8%, and the Harvey test never, because it repeats the regression the weights were estimated from (inst/validation/weighted-fits-size.csv). performVarianceFormTest() accounts for the estimation, held its level in the same designs, and is the test to use on a fitWLS() fit.

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. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Koenker, R. (1981). A note on studentizing a test for heteroscedasticity. Journal of Econometrics, 17(1), 107–112.

Juhl, T., & Sosa-Escudero, W. (2014). Testing for heteroskedasticity in fixed effects models. Journal of Econometrics, 178(Part 3), 484–494. doi:10.1016/j.jeconom.2013.07.005

Examples

 data(mtcars)
 m <- lm(mpg ~ wt + qsec, data = mtcars)
 performKoenkerTest(m, mtcars)

Streaming Koenker studentized Breusch-Pagan test

Description

Evaluates Koenker's studentized variant of the Breusch-Pagan test (the n R^2 statistic from regressing the squared residuals on the regressors) using streamed cross-products so that large datasets can be processed without allocating the full auxiliary regression matrix.

Usage

performKoenkerTestStreaming(
  model,
  data = NULL,
  chunk_size = 10000,
  progress = interactive()
)

Arguments

Details

Chunks of the data contribute to the cross-product matrices for the regression of squared residuals on the model regressors. Aggregating these matrices yields the exact Koenker statistic while bounding memory usage. Sparse cross-products from Matrix are used automatically when available for large chunks.

Since 0.16.0 the chunks are rows of the design of the fit, as in performWhiteTestStreaming(), whose details say what changed.

Value

A htest object equivalent to performKoenkerTest with metadata recording the streaming configuration.

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, data that lack rows or variables of the fit, and, since 0.16.0, a multivariate lm, a fit changed after it was fitted and a fit with model = FALSE whose data are gone, unless data gives its response and fitted values. Which model each test accepts, and why, is set out in test_results, under Supported models.

See Also

performKoenkerTest for the standard implementation.

Examples

data(mtcars)
mod <- lm(mpg ~ wt + qsec, data = mtcars)
performKoenkerTestStreaming(mod, mtcars, chunk_size = 16, progress = FALSE)

Perform Levene's test for equality of variances

Description

Tests equality of variances across groups using residuals from a linear model.

Usage

performLeveneTest(model, data = NULL, group)

Arguments

Details

Absolute deviations from group means are analyzed via a one-way ANOVA. The resulting F statistic has k-1 and N-k degrees of freedom under the null of equal variances.

Value

An object of class htest containing the F statistic, p-value and degrees of freedom.

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. It refuses a multivariate lm ("needs a model with one response"), since 0.16.0; before, it stopped with a coercion error. An mgcv::gam() fit is tested as before 0.16.0, with smooth terms too, until a later step settles which tests take such a fit; since 0.16.0 its observations are matched to the rows of data by the row names of its model frame. Before, they were matched by position, since the residuals of a gam carry no names, and data with the rows of the fit in another order gave a wrong value. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Levene, H. (1960). Robust tests for equality of variances. In Contributions to Probability and Statistics (pp. 278–292). Stanford University Press.

Brown, M. B., & Forsythe, A. B. (1974). Robust tests for the equality of variances. Journal of the American Statistical Association, 69(346), 364–367. doi:10.1080/01621459.1974.10482955

Examples

 data(mtcars)
 mtcars$cyl <- factor(mtcars$cyl)
 m <- lm(mpg ~ wt, data = mtcars)
 performLeveneTest(m, mtcars, "cyl")

Perform McLeod-Li test

Description

Applies the Ljung-Box test to squared residuals.

Usage

performMcLeodLiTest(model, data = NULL, lags = 10)

Arguments

Details

Serial correlation in e_t^2 indicates conditional heteroscedasticity. The Ljung-Box statistic is compared to a chi-square distribution with the chosen number of lags.

Value

An object of class htest containing the test statistic, p-value and degrees of freedom.

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.16.0 it is built on the model context, and also refuses a multivariate lm ("needs a model with one response"), a fit whose coefficients, residuals or weights were changed after the fit, and a fit with model = FALSE whose model frame, rebuilt from the data it was fitted to, no longer holds its response, because those data changed after the fit, with or without data. A fit with model = FALSE whose data are gone is tested on data that give its response and fitted values, and refused without them. The test reads nothing from data, and since 0.16.0 matches it to the observations of the fit as every test does: data with row names of their own in any order, data that number their rows (a tibble, a data.table) in the order of the fit; other numbered data, rows in front of those of the fit and repeated row names are refused, where it gave the value of the fit. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

McLeod, A. I., & Li, W. K. (1983). Diagnostic checking ARMA time series models using squared-residual autocorrelations. Journal of Time Series Analysis, 4(4), 269–273.

Examples

 data(mtcars)
 m <- lm(mpg ~ wt + qsec, data = mtcars)
 performMcLeodLiTest(m, lags = 10)

Perform the non-constant variance (NCV) score test (deprecated)

Description

The Cook-Weisberg score test for non-constant error variance, the diagnostic returned by car::ncvTest().

Usage

performNCVTest(model, data = NULL, var_formula = NULL)

Arguments

Details

Let \hat{e}_i denote the OLS residuals and \hat{\sigma}^2 = \sum_i \hat{e}_i^2 / n the maximum-likelihood variance estimate. The test regresses the scaled squared residuals u_i = \hat{e}_i^2 / \hat{\sigma}^2 on the variance regressors Z and computes half the explained sum of squares of that auxiliary regression. Under homoscedasticity and normal errors the statistic is asymptotically chi-square with degrees of freedom equal to the number of variance regressors.

With var_formula = NULL the fitted values are the sole variance regressor, matching the default of car::ncvTest() and Stata's estat hettest.

Because the divisor 2 is the null variance of \hat{e}^2 / \sigma^2 under normality, the test is sensitive to non-normal errors; performKoenkerTest provides the studentized variant, which is robust to non-normal kurtosis.

Value

An object of class htest containing the chi-square score statistic, its degrees of freedom and the p-value.

Deprecated

performNCVTest() is deprecated since 0.14.0: use performCookWeisbergTest, which gained var_formula and is the one entry point for the Cook-Weisberg score test. With a variance formula, performNCVTest(model, var_formula) returns what performCookWeisbergTest(model, var_formula) returns; without one the two return the same statistic, degrees of freedom and p-value, and performCookWeisbergTest() labels the result as the Cook-Weisberg test. The old name still works, returns what it returned and warns.

Weighted fits

On a fit with weights the test is computed from the Pearson residuals \sqrt{w_i}\, e_i, so it tests whether the weights are adequate; see the section on weighted fits in performKoenkerTest(), which also covers weights estimated by fitWLS().

Validation

Reproduces car::ncvTest() to within 1e-8 for both the default and the var_formula forms; see ‘tests/testthat/test-pass-a-reference.R’. Releases before 0.7.0 regressed the absolute residuals on the fitted values and reported a t statistic, which is a Glejser-type test rather than the Cook-Weisberg score test.

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.16.0 it is built on the model context, and also refuses a multivariate lm ("needs a model with one response"), a fit whose coefficients, residuals or weights were changed after the fit, and a fit with model = FALSE whose model frame, rebuilt from the data it was fitted to, no longer holds its response, because those data changed after the fit, with or without data. A fit with model = FALSE whose data are gone is tested on data that give its response and fitted values, and refused without them. data is matched to the observations of the fit as by performCookWeisbergTest, with or without var_formula. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Cook, R. D., & Weisberg, S. (1983). Diagnostics for heteroscedasticity in regression. Biometrika, 70(1), 1–10. doi:10.1093/biomet/70.1.1

Fox, J., & Weisberg, S. (2019). An R Companion to Applied Regression (3rd ed.). SAGE.

See Also

performCookWeisbergTest, performKoenkerTest, performBPTest

Examples

 data(mtcars)
 m <- lm(mpg ~ wt + qsec, data = mtcars)
 # Deprecated: use performCookWeisbergTest(m) and
 # performCookWeisbergTest(m, var_formula = ~wt)
 suppressWarnings(performNCVTest(m, var_formula = ~wt))

Perform O'Brien test for equality of variances

Description

Implements O'Brien's observation-level transformed-residual test for equality of group variances.

Usage

performOBrienTest(model, data = NULL, group)

Arguments

Details

For each residual within a group, O'Brien's transformation combines its squared deviation from the group mean with the group's sample variance and sample size. A one-way ANOVA of those transformed observations yields an approximate F_{k-1,N-k} test. Every group must contain at least three observations.

Value

An object of class htest containing the F statistic, numerator and denominator degrees of freedom, and p-value.

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. It refuses a multivariate lm ("needs a model with one response"), since 0.16.0; before, it stopped with a coercion error. An mgcv::gam() fit is tested as before 0.16.0, with smooth terms too, until a later step settles which tests take such a fit; since 0.16.0 its observations are matched to the rows of data by the row names of its model frame. Before, they were matched by position, since the residuals of a gam carry no names, and data with the rows of the fit in another order gave a wrong value. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

O'Brien, R. G. (1981). A simple test for variance effects in experimental designs. Psychological Bulletin, 89(3), 570–574.

Examples

 set.seed(1702)
 n <- 25
 d <- data.frame(g = factor(rep(letters[1:3], each = n)), x = rnorm(3 * n))
 d$y <- 2 + d$x + rnorm(3 * n, sd = rep(c(1, 1, 1.8), each = n))
 m <- lm(y ~ x, data = d)
 performOBrienTest(m, d, "g")

Perform Park test for heteroscedasticity

Description

Regresses the log of squared residuals on the log of a suspected variable.

Usage

performParkTest(model, data = NULL, variable)

Arguments

Details

The slope coefficient from \log(e^2) regressed on \log x is tested for zero using a t statistic. A significant coefficient indicates scale depending on x.

Value

An object of class htest containing the t statistic, p-value and degrees of freedom.

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. It refuses a multivariate lm ("needs a model with one response"), since 0.16.0; before, it stopped with a coercion error. An mgcv::gam() fit is tested as before 0.16.0, with smooth terms too, until a later step settles which tests take such a fit; since 0.16.0 its observations are matched to the rows of data by the row names of its model frame. Before, they were matched by position, since the residuals of a gam carry no names, and data with the rows of the fit in another order gave a wrong value. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Park, R. E. (1966). Estimation with heteroscedastic error terms. Econometrica, 34(4), 888. doi:10.2307/1910108

Examples

 data(mtcars)
 m <- lm(mpg ~ wt + qsec, data = mtcars)
 performParkTest(m, mtcars, "wt")

Perform Pesaran CD test

Description

Checks for cross-sectional dependence in panel residuals.

Usage

performPesaranTest(model, data = NULL, id, time)

Arguments

Details

The statistic averages pairwise correlation coefficients of the residuals and is standardized to follow a normal distribution under the null of independence. On a balanced panel of N units observed at the same T periods it is

\mathrm{CD} = \sqrt{\frac{2T}{N(N - 1)}} \sum_{i < j} \hat{\rho}_{ij},

where \hat{\rho}_{ij} is the correlation of the residuals of units i and j. On an unbalanced panel, whose units are not all observed at the same periods, it is computed, since 0.16.0, from the periods each pair of units shares (Pesaran 2004, Section 9, equation 31), as plm::pcdtest(test = "cd") computes it:

\mathrm{CD} = \sqrt{\frac{1}{P}} \sum_{i < j} \sqrt{T_{ij}}\, \hat{\rho}_{ij},

where T_{ij} is the number of periods at which units i and j are both observed, \hat{\rho}_{ij} the correlation of their residuals over those periods, each centred at its mean over them, and P the number of pairs that share at least two periods, the pairs of the sum: N(N - 1)/2 when every pair does. On a balanced panel, where T_{ij} = T, the two are the same statistic. A pair that shares fewer than two periods has no correlation and is left out, with a warning, as plm::pcdtest() leaves it out, and a panel in which no pair shares two periods is refused. A correlation that is not defined, the residuals of a unit constant over the periods of a pair, adds nothing to the sum and the pair stays in the count, as on a balanced panel, with a warning; plm::pcdtest() returns NA there, so the statistic is that of plm::pcdtest() except where a correlation is undefined. Pesaran (2004) states its standard normal limit under the null of independence for T_i > k + 1 and T_{ij} > 3, with N large, T_i being the number of periods of unit i and k the number of regressors of its own regression.

Value

An object of class htest.

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. The model is a pooled lm of the panel, with the units and the periods as columns of data; a plm::plm() fit is not an lm and is refused. Since 0.16.0 a pair of id and time that holds two observations of the fit is refused. An unbalanced panel, as that of a fit that dropped some rows of some units, is tested since 0.16.0, from the periods each pair of units shares (Details), and gives the value of plm::pcdtest(test = "cd") on the residuals of the fit, except where a correlation is undefined (Details): before, the test stopped on it, or, when the number of periods of a unit divided that of the panel, returned a wrong value computed on its residuals recycled over the periods it lacks. A fit with na.action = na.exclude that dropped rows is read as before, by position, until a later step: the panel of the rows of data must then be balanced, and an unbalanced one is refused, with a message that points to na.omit. A fit with na.exclude that dropped no row is tested as any other. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Pesaran, M. H. (2004). General diagnostic tests for cross section dependence in panels. CESifo Working Paper 1229 (also Cambridge Working Papers in Economics 0435). Section 9, “CD tests in the case of unbalanced panels”, equation 31.

Pesaran, M. H. (2015). Testing weak cross-sectional dependence in large panels. Econometric Reviews 34(6-10), 1089–1117. doi:10.1080/07474938.2014.956623.

Croissant, Y. and Millo, G. (2008). Panel data econometrics in R: the plm package. Journal of Statistical Software 27(2), 1–43. doi:10.18637/jss.v027.i02. The statistic of an unbalanced panel is that of plm::pcdtest(test = "cd") (plm 2.6-6), whose help page says its tests are “suitable also for unbalanced panels”.

Examples

 df <- data.frame(id = rep(1:3, each = 5), time = rep(1:5, 3), x = runif(15), y = rnorm(15))
 m <- lm(y ~ x, data = df)
 performPesaranTest(m, df, "id", "time")

Quantile regression heteroscedasticity test

Description

Tests equality of regression slopes across two or more conditional quantiles using the joint Wald-type test implemented by quantreg::anova.rqs(). The procedure accounts for dependence among quantile-specific estimates from the same sample rather than treating their covariance matrices as independent.

Usage

performQuantileRegressionTest(
  model,
  data = NULL,
  taus = c(0.25, 0.75),
  se_type = c("nid", "ker"),
  iid = TRUE
)

Arguments

Details

Under a pure location-shift model with homoskedastic errors, regression slopes are equal across quantiles. Rejection therefore provides evidence against that location-shift/homoskedastic specification. The result should not be interpreted as a universal test for every possible form of heteroscedasticity.

Weighted fits are refused, although the resampling tests accept them. On a fit with weights w_i the null hypothesis is that the Pearson residuals \sqrt{w_i}\, e_i have constant variance, which is a statement about the regression of \sqrt{w_i}\, y_i on \sqrt{w_i}\, x_i. The quantile regressions this test fits are of y on x, where the errors have variance \sigma^2 / w_i by assumption, so their slopes differ across quantiles even when the weights are right. Fitting them on the weighted scale does not by itself repair that: there the conditional quantiles are the regression plus a constant, and the constant is not a column of the weighted regressors. Adding a column of ones would restore the location-shift form, but that version of the test needs its own derivation and its own size study, and has neither.

Value

An htest object containing the F-like joint statistic, numerator and denominator degrees of freedom, p-value, fitted quantiles, and quantile-specific slope estimates.

Supported models

It accepts an unweighted lm and an unweighted glm with the Gaussian family and the identity link and, since 0.16.0, a rank-deficient one: quantreg::rq() refuses a singular design, so the quantile regressions are fitted to the columns that were not aliased, which gives the value of the same model fitted without them. Since 0.16.0 it also accepts a model with an offset: the quantile regressions are fitted to the response of the fit minus its offset, on the design of the fit, which is the model with the offset, since the offset shifts every conditional quantile as it shifts the mean; before, they were refitted from the formula, which left the offset out. Since 0.16.0 it also accepts a model with a single regressor, such as lm(y ~ x); before, it stopped with "subscript out of bounds" on one, and runHeteroTests() returned its fallback, the rank permutation test, under the name quantile_regression. It refuses a weighted fit, a fitWLS() fit among them, as it refits the model without the weights, a glm of another family or link, a survey::svyglm() fit, a coefficient that is not finite for another reason, and data that lack rows or variables of the fit. It refuses a multivariate lm ("needs a model with one response"), and, since 0.16.0, an mgcv::gam() fit with smooth terms: a smooth term is usually penalised, and a penalised fit is not the least-squares fit of the design the test fits again. A later step settles which tests take such a fit; a gam without smooth terms is tested as the Gaussian glm it is. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Koenker, R., & Bassett, G. (1982). Robust tests for heteroscedasticity based on regression quantiles. Econometrica, 50(1), 43–61.

Koenker, R. (2005). Quantile Regression. Cambridge University Press.

Examples

if (requireNamespace("quantreg", quietly = TRUE)) {
  # The test needs at least 40 observations, so mtcars (32) is too small.
  model <- lm(stations ~ mag + depth, data = quakes)
  performQuantileRegressionTest(model, quakes)
}

Ramsey's RESET test for nonlinearity

Description

Adds powers of the fitted values and performs an F-test to check for neglected nonlinearity.

Usage

performRESETTest(model, data = NULL, power = 2:3)

Arguments

Details

On a fit with weights the augmented model is fitted by weighted least squares and both residual sums of squares are weighted. Releases before 0.12.0 refitted the augmented model without the weights, so the two models were not nested.

On a fit with an offset the fitted values are the design times the coefficients plus the offset, and the test is the nested F test of the model with its offset against the model augmented with the powers of those fitted values, fitted with the same offset and weights (since 0.16.0, decided 2026-10-18). It is not RESET of the model fitted to the response minus the offset, which augments that model with the powers of the design times the coefficients alone. Before 0.16.0 the augmented model was fitted without the offset, so the two models were not nested.

Value

An object of class htest.

Supported models

It accepts an lm and a glm with the Gaussian family and the identity link, weighted or not, 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. A fit with a zero weight is accepted too: the test compares two weighted least-squares fits, which both leave the observation out, and gives the value of the fit to the observations with positive weight (kept so, decided on 2026-10-17). It accepts a model with an offset, given as an argument or as an offset() term, weighted or not, and tests the model with its offset (since 0.16.0; see Details). It refuses a glm of another family or link, a survey::svyglm() fit, a coefficient that is not finite for another reason, data that lack rows or variables of the fit, and, since 0.16.0, a fit with model = FALSE whose model frame, rebuilt from the data it was fitted to, no longer holds its response, because those data changed after the fit: it read that frame and returned a value. A fit with model = FALSE whose data are gone is tested on data that give its response and fitted values (since 0.16.0; it refused such a fit until the tests computed from the model alone were built on the model context), and refused without them. Since then it also refuses numbered data in another order than the fit, rows in front of those of the fit and repeated row names, as every test does, where it gave the value of the fit. It refuses a multivariate lm ("needs a model with one response"), which it tested without data before 0.16.0, and, since 0.16.0, an mgcv::gam() fit with smooth terms: a smooth term is usually penalised, and a penalised fit is not the least-squares fit of the design the test fits again. A later step settles which tests take such a fit; a gam without smooth terms is tested as the Gaussian glm it is. Which model each test accepts, and why, is set out in test_results, under Supported models.

Examples

data(mtcars)
m <- lm(mpg ~ wt + qsec, data = mtcars)
performRESETTest(m)

Rank-based permutation heteroscedasticity test

Description

Performs a non-parametric permutation test based on the rank correlation between absolute residuals and a chosen ordering variable. Significant correlations imply systematic changes in residual spread.

Usage

performRankPermutationTest(
  model,
  data = NULL,
  order_by = NULL,
  B = 999,
  progress = interactive()
)

Arguments

Value

An htest object containing the observed Spearman correlation and a permutation-based p-value.

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. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Hollander, M., Wolfe, D. A., & Chicken, E. (2013). Nonparametric Statistical Methods (3rd ed.). Wiley.

Examples

data(mtcars)
model <- lm(mpg ~ wt + hp, data = mtcars)
set.seed(42)
performRankPermutationTest(model, mtcars, B = 199, progress = FALSE)

Scatter-plot diagnostics for heteroscedasticity

Description

Computes correlations between absolute residuals and selected variables.

Usage

performScatterDiagnostic(model, data, vars)

Arguments

Details

High correlations suggest increasing spread with the explanatory variables and motivate variance-stabilizing transformations.

Value

A named numeric vector of correlations.

Supported models

It accepts the models the tests accept: an lm and a glm with the Gaussian family and the identity link, unweighted or with positive weights, a fitWLS() fit and a rank-deficient fit. It refuses a glm of another family or link (since 0.16.0), a survey::svyglm() fit, a fit with a zero weight, and data that lack rows or (since 0.16.0) variables of the fit. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Cleveland, W. S. (1979). Robust locally weighted regression and smoothing scatterplots. Journal of the American Statistical Association, 74(368), 829–836.

Examples

 data(mtcars)
 m <- lm(mpg ~ wt + qsec, data = mtcars)
 performScatterDiagnostic(m, mtcars, c("wt", "qsec"))

Spatial heteroscedasticity test

Description

Evaluates spatial clustering in squared residuals using Moran's I with a permutation reference distribution. Significant positive autocorrelation in the squared residuals is evidence of spatially varying variance.

Usage

performSpatialHeteroTest(
  model,
  data = NULL,
  listw,
  permutations = 999,
  zero.policy = NULL
)

Arguments

Value

An htest result with Moran's I statistic applied to squared residuals.

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. It refuses a multivariate lm ("needs a model with one response"), since 0.16.0; before, it stopped with a coercion error. An mgcv::gam() fit is tested as before 0.16.0, with smooth terms too, until a later step settles which tests take such a fit. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Anselin, L. (1988). Spatial Econometrics: Methods and Models. Kluwer.

Bivand, R. S., Pebesma, E., & Gómez-Rubio, V. (2013). Applied Spatial Data Analysis with R (2nd ed.). Springer.

Examples

if (requireNamespace("spdep", quietly = TRUE)) {
  data(mtcars)
  coords <- cbind(runif(nrow(mtcars)), runif(nrow(mtcars)))
  nb <- spdep::knn2nb(spdep::knearneigh(coords, k = 4))
  lw <- spdep::nb2listw(nb)
  model <- lm(mpg ~ wt + hp, data = mtcars)
  performSpatialHeteroTest(model, mtcars, listw = lw, permutations = 199)
}

Perform Spearman rank correlation test

Description

Computes Spearman's rho between absolute residuals and fitted values.

Usage

performSpearmanTest(model, data = NULL)

Arguments

Details

The squared rank correlation coefficient times the sample size minus 1 approximates a chi-square distribution with one degree of freedom.

Value

An object of class htest containing the test statistic, p-value, degrees of freedom and rho.

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.16.0 it is built on the model context, and also refuses a multivariate lm ("needs a model with one response"), a fit whose coefficients, residuals or weights were changed after the fit, and a fit with model = FALSE whose model frame, rebuilt from the data it was fitted to, no longer holds its response, because those data changed after the fit, with or without data. A fit with model = FALSE whose data are gone is tested on data that give its response and fitted values, and refused without them. The test reads nothing from data, and since 0.16.0 matches it to the observations of the fit as every test does: data with row names of their own in any order, data that number their rows (a tibble, a data.table) in the order of the fit; other numbered data, rows in front of those of the fit and repeated row names are refused, where it gave the value of the fit. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Spearman, C. (1904). The proof and measurement of association between two things. American Journal of Psychology, 15(1), 72–101.

Examples

 data(mtcars)
 m <- lm(mpg ~ wt + qsec, data = mtcars)
 performSpearmanTest(m)

Perform Spread-Level test

Description

Tests the slope of log(|residuals|) on log(fitted).

Usage

performSpreadLevelTest(model, data = NULL)

Arguments

Details

A slope near 0 indicates constant variance. A significant slope suggests a power transformation may stabilize the spread.

Value

An object of class htest.

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.16.0 it is built on the model context, and also refuses a multivariate lm ("needs a model with one response"), a fit whose coefficients, residuals or weights were changed after the fit, and a fit with model = FALSE whose model frame, rebuilt from the data it was fitted to, no longer holds its response, because those data changed after the fit, with or without data. A fit with model = FALSE whose data are gone is tested on data that give its response and fitted values, and refused without them. The test reads nothing from data, and since 0.16.0 matches it to the observations of the fit as every test does: data with row names of their own in any order, data that number their rows (a tibble, a data.table) in the order of the fit; other numbered data, rows in front of those of the fit and repeated row names are refused, where it gave the value of the fit. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Cleveland, W. S. (1979). Robust locally weighted regression and smoothing scatterplots. Journal of the American Statistical Association, 74(368), 829–836.

Examples

 data(mtcars)
 m <- lm(mpg ~ wt + qsec, data = mtcars)
 performSpreadLevelTest(m)

Studentized Breusch–Pagan test (deprecated)

Description

performStudentizedBPTest() is deprecated since 0.14.0: use performKoenkerTest(), which computes the same statistic. It still works and warns, and returns what it returned before except on a fit returned by fitWLS(). The model frame of that fit holds the weights of the model fitWLS() refitted (none, for an unweighted model), not the estimated weights, and before 0.15.0 performStudentizedBPTest() read them from there, so it computed the statistic of the unweighted fit. Since 0.15.0 it reads the weights of the fit, as performKoenkerTest() does, and the two agree. On a model with an offset they agree since 0.16.0: before, the refit of the design left the offset out, and performStudentizedBPTest() returned the statistic of another model.

Usage

performStudentizedBPTest(model, data = NULL)

Arguments

Details

Computes Koenker's (1981) studentized form of the Breusch–Pagan Lagrange multiplier statistic, n R^2 from the regression of the squared residuals on the regressors. The classical Breusch–Pagan statistic assumes normal errors and this one does not: in the simulations of ⁠inst/validation/⁠ it holds its level under t_5 errors, where the classical statistic rejects far too often.

Following Koenker (1981) and the implementation in lmtest::bptest(), the procedure fits an auxiliary regression of e_i^2 - \hat{\sigma}^2 on the regressors from the original model (including the intercept), where e_i denotes the residuals and \hat{\sigma}^2 = \sum e_i^2 / n the mean of their squares. With \hat{g}_i the fitted values of that regression, the statistic T = n \sum \hat{g}_i^2 / \sum (e_i^2 - \hat{\sigma}^2)^2 is its n R^2. Under homoskedasticity it is asymptotically chi-squared with degrees of freedom equal to the number of regressors beyond the intercept, whether or not the errors are normal. On a weighted fit e_i is the Pearson residual \sqrt{w_i} times the raw residual and the auxiliary regression is unweighted; see the section on weighted fits in performKoenkerTest(). This departs from lmtest::bptest(), which keeps the raw residuals and weights the auxiliary regression. The residuals are those of the refit of the design of the fit to its response, with its weights and its offset, so on a model with an offset the statistic is that of the model fitted to the response minus the offset, Koenker's statistic of the model (since 0.16.0). The implementation shares the validation helpers used across the package to ensure that: (i) the model satisfies minimum sample-size thresholds via rvalidateModelInputs(), (ii) the response, design, weights and offset are those of the fit itself, with the rows of data matched to its observations (the model context, since 0.16.0), and (iii) the requirements registered for this test in rvalidateTestRequirements() are met.

Value

An object of class htest reporting the chi-squared statistic and p-value for the null hypothesis of constant error variance.

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. A model with an offset is accepted since 0.16.0: the refit of the design carries the offset, as performKoenkerTest() accounts for it. It refuses a multivariate lm ("needs a model with one response"), and, since 0.16.0, an mgcv::gam() fit with smooth terms: a smooth term is usually penalised, and a penalised fit is not the least-squares fit of the design the test fits again. A later step settles which tests take such a fit; a gam without smooth terms is tested as the Gaussian glm it is. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

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

Davidson, R., & MacKinnon, J. G. (2004). Econometric Theory and Methods. Oxford University Press. Section 7.5 derives the test as n R^2 from the regression of the squared residuals on the variance regressors.

See Also

performBPTest() for the classical LM statistic, which assumes normal errors, and performKoenkerTest(), which computes the same studentized n R^2 statistic on the squared residuals. The robust workflow in performBPTestRobust() augments the studentized statistic with bootstrap diagnostics.

Examples

data(mtcars)
mod <- lm(mpg ~ wt + qsec, data = mtcars)
# Deprecated: warns, and returns the statistic of performKoenkerTest()
suppressWarnings(performStudentizedBPTest(mod, mtcars))
performKoenkerTest(mod, mtcars)

Szroeter test for ordered alternatives

Description

Detects monotonic changes in variance when observations can be meaningfully ordered—typically by time or another covariate—using the cumulative weighting scheme proposed by Szroeter (1978).

Usage

performSzroeterTest(model, data = NULL, order_by,
  alternative = c("greater", "two.sided", "less"))

Arguments

Details

After ordering the residuals \hat{e}_{(i)} by order_by, the test forms the rank-weighted average of the squared residuals

h = \frac{\sum_{i = 1}^n i \, \hat{e}_{(i)}^2}{\sum_{i = 1}^n \hat{e}_{(i)}^2},

which is Szroeter's (1978) class of statistics evaluated at the canonical weights h_i = i. Variance that grows with order_by shifts weight onto the high ranks and drives h upwards.

Since 0.16.0 the null mean and variance of h are those of the design. Let M be the residual-maker of the design in the order of the test (the rows \sqrt{w_i}\, x_i on a weighted fit), A = \mathrm{diag}(1, \ldots, n), k the rank of the design and \nu = n - k. Under normal errors of constant variance (given the weights) the residuals make h a ratio of quadratic forms, h = \sum_j \lambda_j z_j^2 / \sum_j z_j^2, with z_j independent standard normal and \lambda_j the \nu eigenvalues of MAM on the residual space. The ratio is independent of its denominator, so its moments are ratios of moments:

E\,h = \frac{\mathrm{tr}(MA)}{\nu} = \frac{n (n + 1) / 2 - \sum_i i\, l_{(i)}}{\nu}, \qquad \mathrm{Var}\,h = \frac{2\, (\nu\, \mathrm{tr}(MAMA) - \mathrm{tr}(MA)^2)}{\nu^2 (\nu + 2)},

with l_{(i)} the leverages in the order of the test and \mathrm{tr}(MAMA) = \sum_i i^2 - 2 \sum_i i^2 l_{(i)} + \|Q_1' A Q_1\|_F^2, Q_1 an orthonormal basis of the design, so that no eigenvalue is computed and the cost is O(n k^2). The test refers

Q = \frac{h - E\,h}{\sqrt{\mathrm{Var}\,h}}

to the standard normal distribution. The moments are derived here, by the argument that gives those of the Durbin-Watson statistic, and are not taken from Szroeter (1978). Without regressors they are (n + 1) / 2 and (n^2 - 1) / (6 (n + 2)). Before 0.16.0 the test used (n + 1) / 2 and (n^2 - 1) / (6n), Szroeter's asymptotic moments, whatever the design: on a weighted fit whose leverage falls along the ordering the mean of h is above (n + 1) / 2 (by about 0.44 in the design of the validation below, at every n), and the test rejected too often. The statistic h is unchanged; Q and the p-value change for every fit, on lm(mpg ~ wt + qsec, mtcars) ordered by wt from -0.6494 and 0.7420 to -0.6626 and 0.7462.

The test needs at least two residual degrees of freedom, \nu \ge 2. With one the residuals lie on a single line, h takes the same value whatever the errors and its null variance is 0: since 0.16.0 the test refuses such a fit, for which it returned a value before.

The residuals and weights are those of the fit, and order_by is read from the rows of data that hold its observations (the model context, since 0.16.0); the package validation helpers check that the ordering variable is present, and ensure sufficient sample size and variability in the squared residuals.

Ties in order_by keep the order of the rows the model was fitted to. Unless the rows were sorted by the response, that order does not depend on the residuals, so the test keeps its level under the null hypothesis; but with a discrete ordering variable the statistic depends on the order of the rows.

runHeteroTests(tests = "szroeter"), and anything else that runs the registered diagnostic, orders the observations by the fitted values of model, so its alternative is variance that rises with the mean. With a single regressor that is the order of the regressor when its slope is positive and the reverse order when it is negative. The fitted values are computed from the coefficients, so observations with the same regressors tie exactly, as they do under a factor or a binary regressor, and the ties are broken as above. Since 0.16.0 they are read from the fit, in its order, and data is matched to the observations of the fit as this function matches it: without data the registered test reads the model frame, and data in the order of the fit whose row numbers are not those of the fit are taken. Before, it matched the fitted values to the row names of data, and refused both ("Missing rows").

Value

An object of class htest reporting Szroeter's standardised statistic Q, the sample size, the underlying rank-weighted statistic h in estimate, and the p-value from the normal reference distribution, with the null mean and variance of h given the design.

Validation

The statistic and its null moments are checked against an independent reconstruction from the full residual-maker in ‘tests/testthat/test-pass-a-reference.R’, and the moments against the eigenvalues of MAM to 1e-13. Releases before 0.7.0 divided the centred statistic by a further \sqrt{n}, shrinking Q by a factor of roughly 2 / \sqrt{n} and leaving the test with essentially no power against any alternative; see ‘NEWS.md’.

inst/validation/szroeter-design-null-size.R tested the samples of inst/validation/harrison-mccabe-size-power.R, lm(y ~ x1 + x2) ordered by x1, at 5000 replications per cell. With Gaussian errors the test rejected a true null hypothesis at the 5% level 0.048, 0.050 and 0.049 of the time at 40, 100 and 400 observations, and 0.049, 0.047 and 0.050 on a correctly weighted fit (weights 1 / x_1^2), all inside the release gate's interval, 0.042 to 0.058; with the moments of 0.15.0 it rejected 0.044, 0.048 and 0.049, and on the weighted fit 0.067, 0.056 and 0.055. Its moments agreed with those of 100,000 simulated samples of one design of each size, variances within 1%. The exact null distribution of h given the design, from the eigenvalues by Imhof's (1961) method, gave the same rates within 0.0006, at a cost of O(n^3) per call. Under t_5 errors no version holds the level (0.10 to 0.15): the null assumes normal errors, and performGlejserTest() with robust = TRUE is the test to use there.

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. It refuses a multivariate lm ("needs a model with one response"), since 0.16.0; before, it stopped with a coercion error. An mgcv::gam() fit is tested as before 0.16.0, with smooth terms too, until a later step settles which tests take such a fit; since 0.16.0 its observations are matched to the rows of data by the row names of its model frame. Before, they were matched by position, since the residuals of a gam carry no names, and data with the rows of the fit in another order gave a wrong value. The null moments of a gam with smooth terms are those of its design as an unpenalised least-squares fit, with the rank of that design, as performHarrisonMcCabeTest() takes its null distribution. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Szroeter, J. (1978). A class of parametric tests for heteroscedasticity in linear econometric models. Econometrica, 46(6), 1311–1327. doi:10.2307/1913833

Godfrey, L. G. (1988). Misspecification Tests in Econometrics. Cambridge University Press. Section 5.4 outlines the Szroeter test.

Imhof, J. P. (1961). Computing the distribution of quadratic forms in normal variables. Biometrika, 48(3-4), 419–426. doi:10.1093/biomet/48.3-4.419

See Also

performKoenkerTest() for the regressor-based alternative and performGQTest() for split-sample diagnostics on ordered data.

Examples

data(mtcars)
mod <- lm(mpg ~ wt + qsec, data = mtcars)
performSzroeterTest(mod, mtcars, order_by = "wt")

# Detect ordered heteroscedasticity in simulated data with variance increasing
# in an index variable
set.seed(404)
n <- 150
x <- sort(runif(n))
y <- 2 + 0.5 * x + rnorm(n, sd = 0.4 + 0.8 * x)
df <- data.frame(y, x)
performSzroeterTest(lm(y ~ x, data = df), df, order_by = "x")

# Two-sided version when the direction of the change is not known in advance
performSzroeterTest(mod, mtcars, order_by = "wt", alternative = "two.sided")

Variance inflation factors

Description

Compute variance inflation factors for assessing multicollinearity. Each VIF is the corresponding diagonal element of the inverse correlation matrix of the design, so the computation works for models containing factors (each contrast column receives its own VIF).

Usage

performVIFDiagnostic(model)

Arguments

Value

A named numeric vector of VIF values, one per design-matrix column (excluding the intercept). Perfectly collinear columns are reported as Inf.

Examples

data(mtcars)
m <- lm(mpg ~ wt + qsec, data = mtcars)
performVIFDiagnostic(m)

Test the functional form of the error variance

Description

Tests whether the error variance follows a stated parametric function of stated variables, \sigma_i^2 = h(z_i; \gamma). Every other test in the package has constant variance as its null hypothesis, so it can say that the variance moves with some variables but not whether an exponential or a power function of them describes it. Here the variance function is the null hypothesis, and rejecting it means that the chosen shape is inadequate.

Usage

performVarianceFormTest(
  model,
  data = NULL,
  var_formula = NULL,
  form = c("exponential", "power"),
  against = "squares"
)

Arguments

Details

The procedure has three steps.

1. The variance function is fitted to the residuals of model and the model is refitted by weighted least squares, exactly as fitWLS() does.

2. The squared standardized residuals r_i = \tilde{e}_i^2 / \hat\sigma_i^2 of the weighted fit are formed. If the variance function is right their mean does not depend on anything.

3. r_i is regressed on a constant, on g_i = \hat\sigma_i^{-2}\, \partial \sigma_i^2 / \partial \gamma, and on the added terms a_i. The statistic is the F test that the coefficients on a_i are zero.

For the exponential form g_i is z_i, and for the power form it is \log z_i. The regression is the score direction of the variance function extended by the added terms, so the test is a studentized Lagrange multiplier test of the stated form against that extension.

Step 3 keeps g_i in the regression for a reason. Estimating \gamma fits the variance along g_i, so those directions say little about the form: a Breusch–Pagan or Koenker test of the weighted fit on the same variance regressors rejects a wrong form about as often as the right one. The information is in the terms outside the fitted variance function. Conditioning on g_i is also what makes the null distribution free of the error in \hat\gamma to first order, which is the regression-based construction of Wooldridge (1990, 1991).

Value

An object of class htest with the F statistic, its numerator and denominator degrees of freedom, and the p-value for the null hypothesis that the variance function is correctly specified. estimate holds the fitted coefficients of the variance function: slopes of the log-variance for the exponential form and exponents for the power form. The element variance_function records the form, the variance regressors and the added terms that entered the test.

Interpretation and caveats

A rejection says that the added terms move the variance beyond what the stated function allows. A non-rejection says only that those terms do not; it does not establish the form. Two functions that are close over the observed range of z_i cannot be told apart, and the test will then accept both: a variance that is linear in a regressor is accepted as exponential or as a power when the regressor varies over a short range. Running the test for several forms is a comparison of fits, and choosing the form with the largest p-value is not a test.

The test judges the shape, not the fitted coefficients. Because it conditions on g_i, it does not respond to weights that are off along the variance regressors themselves. A test of constant variance on the weighted fit is sensitive to that, with the reservations about its level described under weighted fits in performKoenkerTest().

The F reference distribution assumes that the standardized errors have a constant fourth moment, the assumption behind the Koenker test. It does not assume normality.

With q variance regressors the default adds up to q (q + 1) / 2 terms. When that is large relative to the sample, name the variance regressors with var_formula or the added terms with against.

Validation

No reference implementation exists, so the statistic is checked against a reconstruction from its definition to 1e-8 (tests/testthat/test-variance-form.R) and by simulation. Over 24 designs in which the form under test was right (both forms, constant and heteroscedastic variances, 50 to 400 observations, Gaussian and t_5 errors, 5000 replications each) the rejection rate at the 5% level ran from 4.2% to 5.7%. Power against a log-variance that is quadratic in a regressor was 98% at 150 observations. Against the other log-linear form it was 22% at 150 observations and 74% at 400. See inst/validation/README.md.

Supported models

It accepts an unweighted lm, a fitWLS() fit and, since 0.16.0, an unweighted glm with the Gaussian family and the identity link and a rank-deficient fit. The Gaussian glm is the least-squares fit lm(), with the same design and residuals, from which the variance function is estimated and the model refitted, so the test gives the value of the lm, and fitWLS() refits it as that lm; it refused every glm before 0.16.0. With the exponential form, the default, the variance regressors of a rank-deficient fit span the same space with or without its aliased columns, and the test gives the value of the same model fitted without them. With the power form they are the logarithms of every column of the design, the aliased ones included, and the logarithm of an aliased column is not a combination of the logarithms of the others, so the value is not that of the model without them; it does not depend on the order of the terms, as the value of the model without them does. It refuses weights that did not come from fitWLS(), the prior weights of a Gaussian glm among them, a glm of another family or link, a survey::svyglm() fit, a coefficient that is not finite for another reason, and data that lack rows or variables of the fit. Since 0.16.0 it is built on the model context, and also refuses a multivariate lm, a fit whose coefficients, residuals or weights were changed after the fit, and a fit with model = FALSE whose model frame, rebuilt from the data it was fitted to, no longer holds its response, because those data changed after the fit: it refitted that frame and returned a value. A fit with model = FALSE whose data are gone is tested on data that give its response and fitted values. Without var_formula or against, the test reads nothing from data, which it matches as every test does, and refuses numbered data in another order, rows in front of those of the fit and repeated row names, where it gave the value of the fit. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

Wooldridge, J. M. (1990). A unified approach to robust, regression-based specification tests. Econometric Theory, 6(1), 17–43. doi:10.1017/S0266466600004898

Wooldridge, J. M. (1991). On the application of robust, regression-based diagnostics to models of conditional means and conditional variances. Journal of Econometrics, 47(1), 5–46. doi:10.1016/0304-4076(91)90076-P

Harvey, A. C. (1976). Estimating regression models with multiplicative heteroscedasticity. Econometrica, 44(3), 461–465. doi:10.2307/1913974

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

Carroll, R. J., & Ruppert, D. (1988). Transformation and Weighting in Regression. Chapman & Hall. Chapter 3 treats the estimation and checking of variance functions.

See Also

fitWLS() fits the variance functions tested here. performHarveyTest() and performParkTest() test constant variance against the exponential and the power form respectively.

Examples

# Standard deviation proportional to x, so the variance is a power of x
set.seed(2026)
n <- 300
sim <- data.frame(x = runif(n, 1, 5))
sim$y <- 1 + 2 * sim$x + sim$x * rnorm(n)
mod <- lm(y ~ x, data = sim)

# The power form is the true one; estimate is the fitted exponent
performVarianceFormTest(mod, form = "power")

# The exponential form in x is not
performVarianceFormTest(mod, form = "exponential")

# Test the variance function of a weighted fit directly
wls <- fitWLS(mod, form = "power")
performVarianceFormTest(wls)

# Test a variance function against one particular added term
performVarianceFormTest(mod, form = "exponential", against = ~ log(x))

Perform White's test for heteroscedasticity

Description

Implements White's test on a fitted linear model.

Usage

performWhiteTest(model, data = NULL, cross_products = TRUE, max_interactions = 10)

Arguments

Details

An auxiliary regression of e^2 on all regressors, their squares and cross-products produces R^2. The statistic n R^2 follows a chi-square distribution whose degrees of freedom are the rank of the auxiliary regression less its intercept (see Degrees of freedom).

The test is computed in memory whatever the size of the data. Up to 0.15.0 it switched to performWhiteTestStreaming() on data above 100 MB, which stops on any model with a factor or a poly() term; call that function to stream.

Value

An object of class htest containing the test statistic, p-value and degrees of freedom.

Degrees of freedom

Let Z hold the regressors, their squares and their cross-products. The auxiliary regression of e^2 on Z has an intercept, and the degrees of freedom are the rank of [1, Z] minus one, as lmtest::bptest() counts them given the same terms as varformula. A column of Z that is a linear combination of the intercept and the columns before it is dropped: the square of a dummy, which equals the dummy, and, when the columns of Z span the constant, one column more. They span it with the squares of the dummies of contr.sum, contr.helmert or contr.poly contrasts, the last the default for an ordered factor, the dummies of every level of a factor in a model without an intercept, an interaction with an empty cell such as that of the baseline levels, or the squares of an orthogonal poly() basis. In 0.15.0 the degrees of freedom were the rank of Z alone, one more on such a design, with the same statistic and a larger p-value: on lm(y ~ x1 + f, contrasts = list(f = "contr.sum")), with a factor of two levels, 5 where lmtest::bptest() gives 4, and on lm(y ~ x1 + o), with an ordered factor of four levels, 9 for 8.

The rank is decided, at the default tolerance of qr(), on the auxiliary design of the centred regressors, which spans the same columns with the intercept, so that it does not depend on where the regressors are located; the regression is fitted on Z, or on that design when the fit of Z would lose one of the columns kept. Before, a regressor far from zero, such as x = 10^4 + u or a date written as YYYYMMDD, lost its square, which the rank of [1, Z] took for a combination of the intercept and the regressor: on lm(y ~ x1 + x2) with 10^4 added to x1, 10.34 on 4 degrees of freedom, where the model of x1 as it was gives 11.14 on 5, as both now do. That holds for regressors that enter the model linearly, which are centred one by one. A product or a power of a variable located far from zero inside the model stays ill-conditioned once its columns are centred, and can still lose a column; how far is far depends on the spread of the variable. On 200 observations of two standard normal regressors, lm(y ~ I(x1 + 2e3) * I(x2 + 2e2)), with x1 about 2,000 standard deviations from zero, keeps every coefficient, and White's test gives 46.83 on 7 degrees of freedom where lm(y ~ x1 * x2), the same model, gives 51.90 on 8; lm(y ~ I(x1 + 2e3) + I((x1 + 2e3)^2)) gives 8.05 on 3 for 11.29 on 4. In ten such data sets, of normal and of uniform regressors, the degrees of freedom held, and the statistic moved by at most 3e-5 relative, while the shift was at most 700 standard deviations of the regressor; the column was lost from 800 to 2,000 standard deviations on, depending on the term and the data. From about 3,000 to 11,000 standard deviations on, lm() itself aliases the interaction or the square, before White's test sees it, and the test is that of the model without the term: with I(x1 + 1e4) * I(x2 + 1e3), 35.02 on 5, the value of lm(y ~ x1 + x2). A year from 1990 to 2020, about 220 standard deviations from zero, and its square, year * z or income * age are far from that, and give the value of the centred model to seven digits or more.

On a design whose auxiliary columns do not span the constant and whose regressors are not located far from zero, such as continuous regressors or factors with the default treatment contrasts, nothing changed. The same holds with cross_products = FALSE or more regressors than max_interactions, and for the tests computed through this function: performWhiteTestRobust(), whose interval, effect size and power follow the degrees of freedom, and the registered "white" of runHeteroTests(). The p-value of performWhiteTestBootstrap() depends on the statistic alone and does not change on a design that spans the constant. performWhiteTestStreaming() stops on such a design, whose auxiliary columns are collinear with its intercept, and where it gives a value it gives these degrees of freedom, but on an mgcv::gam() fit with smooth terms, whose smooth terms it does not test. On an auxiliary design whose columns are nearly collinear, as the squares and products of the basis of a smooth term can be, the columns kept can differ from those kept in 0.15.0, so the statistic can change too: there the auxiliary regression of 0.15.0 could silently lose a column to its intercept while its degrees of freedom still counted it, as on a mgcv::gam(y ~ s(x1) + x2) fit whose statistic moves by 0.28%, to the value of lmtest::bptest() with the same terms.

Too few observations

An auxiliary regression that fits the squared residuals exactly has an R^2 of 1, and the statistic is the number of observations, whatever the data. It fits them exactly with as many independent terms as observations, and with fewer when the squared residuals are a combination of its terms: on a paired design, such as lm(y ~ treat + subject) with two observations of each subject, the two residuals of a subject sum to zero, so their squares are equal and a combination of the subject dummies. On 30 subjects the test gave 60 on 30 degrees of freedom (p = 0.0009) on homoscedastic data. Since 0.16.0 it refuses such a regression, where it returned that statistic, and says what to use, as the call allows: cross_products = FALSE when the products were in and the regressors alone do not fit the squared residuals exactly, Koenker's test when they do not, and more data when the regression has as many terms as observations. Six continuous regressors have 27 squares and products, 28 terms with the intercept, so they need at least 29 observations. With fewer terms than observations the design forces the exact fit, and more data of the same design do not help: 500 subjects of two observations are refused as 30 are, and the refusal asks for more observations per level.

A fit counts as exact when its residual sum of squares is at most sqrt(.Machine$double.eps), about 1.5 \times 10^{-8}, of the total sum of squares of the squared residuals. An exact fit leaves rounding error, about 10^{-30} of the total on the paired design; squared residuals that depart from the fit by 2 \times 10^{-4} of their size leave about 2 \times 10^{-7}, and are tested.

The robust form, the registered "white", performWhiteTestBootstrap() and performWhiteTestStreaming() refuse both kinds of exact fit too. runHeteroTests() falls back on White's test without cross products, as on any refusal of White's test, unless the regressors alone fit the squared residuals exactly, as on the paired design, or that is not known, as in the streaming test: the Breusch-Pagan and Koenker tests, its other fallbacks, regress the squared residuals on the regressors and would return the number of observations too, so the runner stops with the refusal.

The refusal covers the exact fit only, and one step away from it the test still rejects almost always: in 2,000 homoscedastic draws of lm(y ~ x + subject), with normal x and errors, White's test rejected at the 5% level in 99% of them with 30 subjects of two observations and one of them observed a third time, and in all but one with 30 pairs whose weights vary within pairs. Those are fixed effects for small groups, on which the test warns (see Fixed effects for small groups).

Fixed effects for small groups

White's, Koenker's and the Breusch-Pagan tests are not valid with fixed effects for small groups: a factor with many levels and few observations of each, such as subject in lm(y ~ x + subject) with each subject observed a few times. On homoscedastic data they reject too often, and the more often the more levels the factor has and the fewer observations each level has. The excess comes from the dummies of the factor in the auxiliary regression. The residuals of a level are taken from its own mean, so under normal errors the squared residuals of a level of T observations vary between levels T / (T - 1) times as much as the auxiliary regression assumes, and with G levels the statistic exceeds its degrees of freedom by about (G - 1) / (T - 1). In the study below Koenker's statistic had that mean, and the squared residuals regressed on x alone, without the dummies, were rejected in 4.6% to 5.5% of the draws by Koenker's statistic and in 4.3% to 5.6% by the Breusch-Pagan statistic (inst/validation/fixed-effects-dummies.R). The excess grows with G at a fixed T, while the spread of the chi-squared distribution grows as \sqrt{2G}, so with a fixed number of observations per level the size goes to 1 as the number of levels grows.

Since 0.16.0 the test warns when a term of the model that holds a factor or a character variable has more levels than observations per level: L > n / L, that is L^2 > n, with L the number of distinct combinations of the factors of the term among the n observations of the fit. The term is a factor, such as subject; an interaction of factors, counted by its cells, such as a:b of a crossed a * b or school:class of a nested school / class; or a slope for each level, x:subject. A variable left out of the formula, as in y ~ . - id, does not count. The warning names the term and gives L, n and the mean number of observations per level; the statistic, the degrees of freedom and the p-value are those the test computes without it. It has class heteroTests_fixed_effects_warning, which suppressWarnings(expr, classes = "heteroTests_fixed_effects_warning") muffles, and no other warning. It is given on every call, also when cachedTest() or runHeteroTests() answers from its cache, and runHeteroTests() and runDiagnostics() give it once for each test that gives it. performWhiteTestRobust() and the registered "white" give it too, and so does performWhiteTestStreaming() where it computes: it refuses a factor that is a main effect, whose dummies and their squares are collinear, and computes White's statistic on a slope for each level. performWhiteTestBootstrap(), whose p-value is a bootstrap's, held its level on these designs and does not give it.

On a balanced design of G subjects of T observations the rule is G \le T. In the study it kept the three tests within 2.5% to 7.5% at the 5% level where it holds, and put at least one of them outside where it does not, in 68 of the 69 designs where they computed; the exception is the Breusch-Pagan test on 20 subjects of 20, at 8.0%. Unbalanced groups, a second regressor and a regressor correlated with the subject effects followed the same pattern, with the mean number of observations per level, n / G, in place of T. The rule is a warning, not a correction, and many panels fail it: on 100 firms observed over 10 years the three tests reject in 18% to 21% of the draws. Rejection at the 5% level, in percent, in 2,000 homoscedastic draws of lm(y ~ x + subject), with x, the subject effects and the errors standard normal:

A slope for each level counts too, but there the error goes the other way: with a slope for each subject, lm(y ~ x + x:subject), on 30 subjects of 3 and of 4 and on 100 of 8, the three tests rejected too rarely, White's in 0.7% to 2.2% of the draws, Koenker's in 0.1% to 0.8% and the Breusch-Pagan test in 0.4% to 2.8% (inst/validation/fixed-effects-variants-size.R). The warning therefore says that the p-value is unreliable: in the simulations usually far too small, and too large with a slope per level.

With at most max_interactions regressors White's auxiliary regression holds the products of x with the dummies, a slope for each subject, and on 5 subjects of 4 it rejected too rarely, in 1.0%. The study is inst/validation/fixed-effects-size.R, with its results in fixed-effects-size.csv (system.file("validation", package = "heteroTests")). Juhl and Sosa-Escudero (2014) derive a test of heteroskedasticity for this setting, after fixed effects estimation of a panel of many units observed a fixed number of times; it needs to be told which factor is the unit, and the package does not implement it.

Weighted fits

On a fit with weights the test is computed from the Pearson residuals \sqrt{w_i}\, e_i, so it tests whether the weights are adequate; see the section on weighted fits in performKoenkerTest(), which also covers weights estimated by fitWLS().

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. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

White, H. (1980). A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity. Econometrica, 48(4), 817–838. doi:10.2307/1912934

Greene, W. H. (2018). Econometric Analysis (8th ed.). Pearson.

Juhl, T., & Sosa-Escudero, W. (2014). Testing for heteroskedasticity in fixed effects models. Journal of Econometrics, 178(Part 3), 484–494. doi:10.1016/j.jeconom.2013.07.005

Examples

 data(mtcars)
 m <- lm(mpg ~ wt + qsec, data = mtcars)
 performWhiteTest(m, mtcars)

Bootstrap White test

Description

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

Details

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

1. Compute White's LM statistic n R^2 from the auxiliary regression on squares and cross-products of the regressors.

2. Generate B bootstrap samples by resampling the centred residuals with replacement, forming y^{*(b)} = \hat{y} + \hat{e}^{*(b)}.

3. Refit the model to each bootstrap sample and recompute White's statistic T^{*(b)} using the same auxiliary specification.

4. 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.

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.

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 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. A transformed response, such as log(y) in log(y) ~ x, is regenerated on the scale the model uses, and the refit reads it in place of log(y).

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.

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. 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

See Also

performWhiteTest() for the classical statistic and performWhiteTestRobust() for enriched diagnostics.

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)

Robust White test with bootstrap and effect sizes

Description

Extends performWhiteTest() with optional bootstrap resampling, confidence intervals, effect size reporting, and power analysis.

Usage

performWhiteTestRobust(
  model,
  data = NULL,
  method = c("standard", "reduced"),
  bootstrap = FALSE,
  B = 999,
  ci_level = 0.95,
  parallel = FALSE
)

Arguments

Details

Provides enriched inference around the White test by combining bootstrap resampling (Efron & Tibshirani, 1993) with asymptotic approximations. The robust_details element summarises interval estimates, effect sizes, and power calculations to aid decision-making.

Value

An object of class htest augmented with a robust_details list containing bootstrap, effect size, and power information.

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. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

White, H. (1980). A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroscedasticity. Econometrica, 48(4), 817–838.

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

See Also

performWhiteTest() for the base statistic and performWhiteTestBootstrap() for a lighter-weight resampling option.

Examples

data(mtcars)
mod <- lm(mpg ~ wt + qsec, data = mtcars)
performWhiteTestRobust(mod, mtcars, bootstrap = TRUE, B = 200)

Memory-efficient White test for large datasets

Description

Computes White's statistic via chunked cross-products rather than fitting the full auxiliary regression in memory. This streaming approach allows the test to scale to datasets that would otherwise exhaust available RAM.

Usage

performWhiteTestStreaming(
  model,
  data = NULL,
  chunk_size = 10000,
  cross_products = TRUE,
  max_interactions = 10,
  progress = interactive()
)

Arguments

Details

The streaming implementation iteratively builds the cross-product matrices required for the auxiliary regression without materialising the full design matrix. Each chunk contributes to X'X, X'y, and summary statistics for the response. The final n R^2 statistic is then computed exactly as in the standard White test, ensuring numerical equivalence while dramatically reducing peak memory usage. When Matrix is installed, sparse cross-products are leveraged automatically for large chunks.

Since 0.16.0 the chunks are rows of the design of the fit, read from the fit itself as performWhiteTest() reads it, so the value does not depend on the size of the chunks but by rounding; the auxiliary design, several times wider, is built chunk by chunk. Up to 0.16.0 each chunk built its own design from the formula, on its rows of data, so a term computed from the rows it is given took other values in each chunk than in the fit: I(x > median(x)) the median of the chunk, and cut(x, quantile(x)) its breaks, which stopped the test ("factor has new levels") in chunks of any size on a fit that dropped the row at the minimum of x. An mgcv::gam() fit with smooth terms is tested as before, on the design of the variables of its formula, which its model frame holds, rather than on the basis of its smooths, on which performWhiteTest() tests it.

Value

A htest object containing the chi-squared statistic, p-value, and metadata about the chunked computation.

Supported models

It accepts an lm and, since 0.16.0, 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. The Gaussian glm with prior weights w is the least-squares fit lm(..., weights = w), with the same design and residuals, and the test gives the value of the lm. It refused every glm before 0.16.0, because the checks it runs read the R-squared of summary(), which a glm does not report; they read that of the least-squares fit now. 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, data that lack rows or variables of the fit, and, since 0.16.0, a multivariate lm, a fit changed after it was fitted and a fit with model = FALSE whose data are gone, unless data gives its response and fitted values (see test_results). It stops when terms of its auxiliary regression are collinear, as the squares of dummy variables are, which performWhiteTest() drops, and, since 0.16.0, with the refusal of performWhiteTest() (see its section Too few observations) when its independent terms are at least as many as the observations, where it stopped as on collinear terms, or fit the squared residuals exactly, collinear or not, as on a paired design. Which model each test accepts, and why, is set out in test_results, under Supported models.

References

White, H. (1980). A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroscedasticity. Econometrica, 48(4), 817–838.

See Also

performWhiteTest() for the exact computation and performWhiteTestRobust() for enhanced reporting.

Examples

data(mtcars)
performWhiteTestStreaming(lm(mpg ~ wt + qsec, data = mtcars), mtcars, chunk_size = 16)

Wild bootstrap heteroscedasticity test

Description

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

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.

Value

An object of class htest containing the observed Breusch-Pagan statistic, bootstrap p-value, and additional details under the bootstrap element.

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.

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

Plot diagnostics

Description

Produce basic residual diagnostic plots.

Usage

## S3 method for class 'HeteroDiagnostic'
plot(x, plots = c("residuals_fitted", "spread_level",
  "density", "qq", "bubble_variance"), ...)

Arguments

Value

List of ggplot objects.

Plot method for power analysis results

Description

Plot method for power analysis results

Usage

## S3 method for class 'power_analysis'
plot(x, ...)

Arguments

Compare residuals before and after remediation

Description

Draws the residuals of two models against their fitted values in two panels, side by side, to show what a remediation method (a weighted fit, a transformation) did to the spread.

Usage

plotBeforeAfter(original, remedied)

Arguments

Details

A weighted fit is shown through its Pearson residuals \sqrt{w_i}\, e_i, so the plot shows whether the weighting flattened the spread. Its raw residuals would look as heteroscedastic as the original ones however good the weights were.

Each panel has its own axes. Weighted residuals and the residuals of a transformed response are not in the units of the original ones, and on a common axis one of the two panels could be squeezed flat. The panels are there to compare the shape of the two clouds, not their size.

Value

A ggplot object with one panel per model. Its data have the columns fitted, resid, model ("original" or "remedied") and panel.

Earlier versions

Up to 0.12.0 the two sets of residuals were overlaid in one panel and told apart by colour.

Examples

data(mtcars)
m1 <- lm(mpg ~ wt, data = mtcars)
m2 <- fitWLS(m1)
plotBeforeAfter(m1, m2)

Bubble plot of residual variance by covariate

Description

Displays residual magnitude against a predictor with bubble size proportional to ⁠|residual|⁠.

Usage

plotBubbleVariance(model, variable = NULL)

Arguments

Value

A ggplot object.

Examples

data(mtcars)
m <- lm(mpg ~ wt + qsec, data = mtcars)
plotBubbleVariance(m, "wt")

Generate a suite of diagnostic plots

Description

This convenience wrapper returns residual-vs-fitted and spread-level plots to help visually assess heteroscedastic patterns.

Usage

plotDiagnosticSuite(model, enhanced = FALSE)

Arguments

Value

A list with elements residuals_fitted and spread_level, each a ggplot object.

Examples

data(mtcars)
m <- lm(mpg ~ wt + qsec, data = mtcars)
plots <- plotDiagnosticSuite(m)
plots$residuals_fitted

Enhanced diagnostic plot suite (deprecated)

Description

Deprecated since 0.14.0: use plotDiagnosticSuite(model, enhanced = TRUE), which returns the same list. This name still works and warns.

Usage

plotDiagnosticSuiteEnhanced(model)

Arguments

Value

A list of ggplot objects.

Density plot of residuals

Description

Shows the distribution of residuals with a kernel density estimate.

Usage

plotResidualDensity(model)

Arguments

Value

A ggplot object.

Examples

data(mtcars)
m <- lm(mpg ~ wt + qsec, data = mtcars)
plotResidualDensity(m)

QQ plot of residuals

Description

Visualises departure from normality using a QQ plot.

Usage

plotResidualQQ(model)

Arguments

Value

A ggplot object.

Examples

data(mtcars)
m <- lm(mpg ~ wt + qsec, data = mtcars)
plotResidualQQ(m)

Plot residuals vs fitted values

Description

Generates a simple scatter plot of residuals against fitted values from a linear model. A horizontal reference line at zero is added.

Usage

plotResidualsFitted(model, enhanced = FALSE)

Arguments

Details

For a weighted fit the Pearson residuals \sqrt{w_i}\, e_i are plotted, as in plot.lm(): those are the residuals that have constant variance when the weights are right. The other residual plots in the package do the same.

With enhanced = TRUE the plot marks influential observations (Cook's distance above 4/n) in red, sizes the points by the absolute residual and adds the standard-error band of the LOESS smooth. This is what the deprecated plotResidualsFittedEnhanced() drew.

Value

A ggplot object.

Examples

data(mtcars)
m <- lm(mpg ~ wt + qsec, data = mtcars)
plotResidualsFitted(m)
plotResidualsFitted(m, enhanced = TRUE)

Enhanced residuals vs fitted plot (deprecated)

Description

Deprecated since 0.14.0: use plotResidualsFitted(model, enhanced = TRUE), which draws the same plot. This name still works and warns.

Usage

plotResidualsFittedEnhanced(model)

Arguments

Value

A ggplot object.

Spread-Level plot for variance diagnostics

Description

Plots the square root of the absolute residuals against fitted values. A lowess smooth is added to highlight trends.

Usage

plotSpreadLevel(model)

Arguments

Value

A ggplot object.

Examples

data(mtcars)
m <- lm(mpg ~ wt + qsec, data = mtcars)
plotSpreadLevel(m)

Print method for enhanced White test results

Description

Print method for enhanced White test results

Usage

## S3 method for class 'htest_enhanced'
print(x, ...)

Arguments

Test-specific sample size requirements

Description

Central catalogue of the minimum sample sizes required by heteroscedasticity diagnostics supported in the package. Each entry either provides a requirements list with human-readable rationale or a function that computes the minimum number of observations from contextual arguments (e.g. lag lengths for dynamic tests).

Usage

rTEST_REQUIREMENTS

Format

A named list keyed by lowercase test identifiers.

Details

Static entries are named lists with some combination of the following components:

min_obs

Overall minimum observation count.

min_obs_per_group

Minimum number of observations that each group should contribute.

reason

Short explanation describing why the requirement exists.

Dynamic entries are functions that receive additional arguments supplied to rvalidateSampleSize() and should return either a numeric minimum or a list containing min_obs and, optionally, reason.

Effect size calculations for heteroscedasticity diagnostics

Description

Converts chi-squared statistics produced by heteroscedasticity tests into interpretable effect sizes such as Cramer's V, the phi coefficient, or an eta-squared analogue. The helper also provides qualitative magnitude descriptors and a brief interpretation string that can be surfaced to users.

Usage

rcalculateEffectSize(
  test_result,
  model,
  data,
  type = c("cramers_v", "phi", "eta_squared")
)

Arguments

Value

A named list with elements effect_size, magnitude, practical_significance, interpretation, and type.

Deprecated

Deprecated since 0.14.0: an internal helper of the package, which leaves the exports in 1.0.0. It still works and warns.

Examples

data(mtcars)
model <- lm(mpg ~ wt + cyl, data = mtcars)
result <- performWhiteTest(model, mtcars)
rcalculateEffectSize(result, model, mtcars)

Register a custom diagnostic

Description

Adds a function to the diagnostics registry so it can be called by runHeteroTests.

Usage

registerDiagnostic(name, fun)

Arguments

Details

A diagnostic registered this way is never cached by cachedTest or runHeteroTests, nor is a test of the package registered again under its name: whether it draws random numbers, or reads anything besides the model and the data, cannot be known.

Value

Invisibly returns NULL.

See Also

listDiagnostics for the names already registered.

Examples

custom <- function(model, data) list(statistic = 0)
registerDiagnostic("custom", custom)

Register a diagnostic plot

Description

Register a diagnostic plot

Usage

registerPlot(name, fun)

Arguments

Value

Invisibly returns NULL.

Examples

registerPlot("custom_plot", function(model) ggplot2::ggplot())

Estimate test power from observed effect size

Description

Uses the observed chi-squared statistic and effect size to approximate the achieved power of a heteroscedasticity test and to suggest an increased sample size required to attain the desired power threshold.

Usage

restimate_test_power(
  test_result,
  model,
  data,
  alpha = 0.05,
  target_power = 0.8,
  max_multiplier = 5
)

Arguments

Value

A list with elements power, recommended_n, and details containing auxiliary information (degrees of freedom, non-centrality parameter, and a power curve over a grid of effect sizes).

Deprecated

Deprecated since 0.14.0: an internal helper of the package, which leaves the exports in 1.0.0. It still works and warns.

Handle missing values according to a specified strategy

Description

Provides a centralized mechanism for dealing with missing values in heteroscedasticity tests. The function can either drop incomplete cases, emit warnings, or fail fast when missingness is not permitted.

Usage

rhandleMissingValues(data, variables, strategy = "complete_cases")

Arguments

Value

A list with components data (the processed data frame), removed_cases (row indices removed), removed_count (number of removed observations), removed_fraction (proportion removed relative to the original data), removed_variables (variables with observed missingness), and loss_message (the formatted warning text).

Robust heteroscedasticity diagnostics

Description

Provides enhanced implementations of core heteroscedasticity diagnostics with bootstrap support, effect size calculations, and statistical power summaries. The module also exposes orchestration helpers for running the enhanced diagnostics in batch and utilities for validating results against reference implementations from other packages.

Process validation helper outputs

Description

Standardises how validation helper results are handled across the diagnostic implementations.

Usage

rprocessValidationResult(result)

Arguments

Value

Invisibly returns result after issuing warnings or errors.

Run enhanced heteroscedasticity diagnostics

Description

Executes one or more robust diagnostics, with an automatic bootstrap for small samples. The Breusch-Pagan entry is Koenker's studentized statistic, performBPTestRobust(studentize = TRUE).

Usage

runAdvancedDiagnostics(
  model,
  data,
  tests = "all",
  auto_enhance = TRUE,
  bootstrap_B = 999,
  parallel = FALSE,
  ci_level = 0.95
)

rrunAdvancedDiagnostics(
  model,
  data,
  tests = "all",
  auto_enhance = TRUE,
  bootstrap_B = 999,
  parallel = FALSE,
  ci_level = 0.95
)

Arguments

Value

A list with the executed results and metadata describing the enhancements that were applied.

Deprecated name

rrunAdvancedDiagnostics() is deprecated since 0.15.0: use runAdvancedDiagnostics(), the same function under the package's naming scheme. The old name still works, returns what it returned and warns; it leaves the exports in 1.0.0.

Run registered diagnostic plots

Description

Run registered diagnostic plots

Usage

runDiagnosticPlots(
  model,
  plots = c("residuals_fitted", "spread_level", "density", "qq", "bubble_variance")
)

Arguments

Value

Named list of ggplot objects

Examples

runDiagnosticPlots(lm(mpg ~ wt, mtcars))

Run a suite of model diagnostics

Description

Combines heteroscedasticity checks with tests for collinearity, nonlinearity and influential observations.

Usage

runDiagnostics(
  model,
  data = NULL,
  tests = c("white", "koenker"),
  power = 2:3,
  use_cache = TRUE,
  chunk_threshold_mb = Inf,
  chunk_size = 10000,
  progress = interactive()
)

Arguments

Value

A list with the results of all diagnostics.

Examples

data(mtcars)
runDiagnostics(mpg ~ wt + qsec, mtcars)

Run multiple heteroscedasticity diagnostics

Description

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

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.

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 the streaming White test stops 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.

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".

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

Parallel test execution helper

Description

Executes a set of registered diagnostics in parallel when multiple CPU cores are available, falling back to sequential execution otherwise.

Usage

runHeteroTestsParallel(model, data, tests, n_cores = NULL)

Arguments

Details

When n_cores exceeds one the diagnostics run on a socket cluster of that many workers (parallel::makeCluster()), started for the call, each of which loads the installed heteroTests. The model and the data travel to the workers with each task. With n_cores = 1, or when the default comes to one worker, the diagnostics run in the session.

Test i draws its random numbers from stream i of a L'Ecuyer-CMRG generator seeded by one draw from the session's generator, on the workers and in the session alike. A seeded call that runs a resampling test, such as "wild_bootstrap", therefore returns the same value each time and whatever n_cores, n_cores = 1 included. It is not the value that test returns when called directly from the same seed, which draws from the session's generator; both are bootstraps of the same data.

Value

Named list of test results in the order provided by tests.

See Also

cachedTest() for memoised execution of individual diagnostics.

Run multivariate heteroscedasticity tests

Description

Convenience wrapper for multivariate tests such as Box's M.

Usage

runMultivariateTests(data, group, tests = c("box_m"))

Arguments

Value

A named list of htest objects.

Examples

runMultivariateTests(iris[,1:4], iris$Species)

Run panel-data heteroscedasticity tests

Description

Breusch-Pagan LM and Pesaran CD tests for panel models.

Usage

runPanelTests(model, data, id, time = NULL,
          tests = c("bp_random", "pesaran"))

Arguments

Value

A named list of htest objects.

See Also

runTimeSeriesTests

Examples

 df <- data.frame(id = rep(1:3, each = 4), time = rep(1:4, 3),
                  x = runif(12), y = rnorm(12))
 m <- lm(y ~ x, data = df)
 runPanelTests(m, df, id = "id", time = "time")

Run heteroscedasticity diagnostics on the data of a survey design

Description

A convenience wrapper that takes the data out of a survey::svydesign() object, fits the model by ordinary least squares and runs runHeteroTests. It does not use the sampling weights, strata or clusters of the design, and warns when the design has unequal weights or clusters.

Usage

runSurveyHeteroTests(formula, design, tests = c("white", "koenker"), ...)

Arguments

Details

The tests are the ones runHeteroTests runs on any lm: they treat the observations as an independent sample with equal weights. They hold their level on survey data when the design can be ignored given the regressors, that is when the chance of being sampled depends only on variables in the model, and the sample has no clusters. When selection depends on the response, or observations are clustered, they reject a model with constant error variance far more often than their nominal level.

The package has no design-based test of constant variance. Sampling weights are not precision weights: they say how many population units an observation stands for and nothing about its error variance. Passing them to lm(weights = ), or testing a survey::svyglm() fit, makes the tests examine residuals multiplied by the square roots of the weights, and a homoscedastic model is then rejected whenever the weights vary with the regressors. svyglm fits are refused for that reason.

inst/validation/survey-designs-size.R measures all of this, together with two design-based tests that were tried and not adopted because they reject too often at the sample sizes examined.

Value

An object of class hetero_test_suite. Since 0.16.0 it carries attr(, "model_origin") = "survey": the model is fitted here, so summary() shows its formula rather than its call (see hetero_test_suite).

Earlier versions

Up to 0.12.0 this page described a survey-weighted fit. The weights did not reach the fit, which was by ordinary least squares, as it is now, so the reported values are unchanged.

See Also

runHeteroTests

Examples

if (requireNamespace("survey", quietly = TRUE)) {
  data(api, package = "survey")
  design <- survey::svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat)
  # Warns: the weights of this design are unequal and are not used.
  res <- runSurveyHeteroTests(api00 ~ api99 + ell, design)
  generics::tidy(res)
}

Run time-series heteroscedasticity tests

Description

Convenience wrapper for Engle's ARCH LM and McLeod-Li tests.

Usage

runTimeSeriesTests(model, lags = 1, tests = c("arch_lm", "mcleod_li"))

Arguments

Value

A named list of htest objects.

See Also

runPanelTests

Examples

 data(mtcars)
 m <- lm(mpg ~ wt + qsec, mtcars)
 runTimeSeriesTests(m, lags = 2)

Benchmark heteroTests against reference implementations

Description

run_benchmark_suite() evaluates a selection of heteroscedasticity diagnostics across varying sample sizes while measuring runtime, memory allocation, and statistical agreement with the lmtest and car packages. The companion generate_benchmark_report() helper aggregates the raw output into tidy summary tables and recommendations.

Usage

run_benchmark_suite(
  sample_sizes = c(100L, 500L, 1000L, 10000L, 100000L, 1000000L),
  tests = NULL,
  replicates = 3L,
  hetero_patterns = c("none", "linear", "group"),
  hetero_strength = 1,
  baseline_packages = c("lmtest", "car"),
  seed = 123L,
  n_predictors = 4L,
  profile_memory = TRUE,
  progress = interactive()
)

generate_benchmark_report(benchmark_results, accuracy_tolerance = 1e-4)

Arguments

Details

Datasets are generated from a linear regression model with a configurable number of predictors and variance patterns (none, linear, group, or exponential). Each diagnostic from heteroTests is executed alongside equivalent tests from lmtest and car when available. Memory profiling is enabled via bench when installed; otherwise only runtimes are collected. The output is suitable for regression testing (to detect performance regressions) and for building artefacts such as CSV benchmark reports.

Value

run_benchmark_suite() returns a list with four elements:

performance

Data frame containing runtime, memory, and inference metrics for every test, sample size, and replicate.

accuracy

Data frame recording absolute differences in p-values and test statistics relative to each baseline implementation.

summary

List of aggregated medians for quick inspection of speed, memory, and accuracy.

metadata

Metadata describing the benchmark configuration (sample sizes, baselines, seed, and availability of optional dependencies).

generate_benchmark_report() converts the raw results into aggregated speed, memory, accuracy, recommendation, and scalability tables alongside the original metadata.

Examples

if (requireNamespace("lmtest", quietly = TRUE) &&
    requireNamespace("car", quietly = TRUE)) {
  results <- run_benchmark_suite(
    sample_sizes = c(100L, 250L),
    replicates = 1L,
    profile_memory = FALSE,
    progress = FALSE
  )
  report <- generate_benchmark_report(results)
  head(report$recommendations)
}

Validate data inputs for heteroscedasticity tests

Description

Performs standard checks on input data.frames prior to running heteroscedasticity diagnostics. The routine verifies that required variables are available and that the sample size meets minimum criteria.

Usage

rvalidateDataInputs(data, required_vars = NULL, min_obs = 10)

Arguments

Value

Invisibly returns the validated data object.

Validate distributional assumptions

Description

Assesses frequently used distributional assumptions before running heteroscedasticity diagnostics. The checks cover approximate normality, positivity for log-transformed quantities, sufficient variation, and the presence of extreme outliers.

Usage

rvalidateDistributionalAssumptions(data, assumptions = list())

Arguments

Value

A list containing passed (logical flag), messages (character vector of violations), warnings (character vector of recoverable issues), and details (named list with diagnostic information for each assumption).

Validate grouping variable structure

Description

Ensures that grouping variables supplied to heteroscedasticity tests satisfy minimum size and level requirements.

Usage

rvalidateGroupingVariable(data, group_var, min_group_size = 3, min_groups = 2)

Arguments

Value

A list mirroring the structure of rvalidateDistributionalAssumptions() with information about the evaluated grouping variable.

Validate regression model inputs

Description

Ensures that a supplied regression model satisfies the basic assumptions required by the heteroscedasticity testing infrastructure. The checks include class validation, successful fitting, availability of finite residuals, and minimum sample size requirements.

Usage

rvalidateModelInputs(model, test_name, min_obs = 10, allow_aliased = FALSE)

Arguments

Value

Invisibly returns model when validation passes.

Validate test-specific sample size requirements

Description

Evaluates whether the supplied data meet the minimum observation counts required by individual heteroscedasticity diagnostics.

Usage

rvalidateSampleSize(test_name, model = NULL, data = NULL, groups = NULL, ...)

Arguments

Value

A validation result describing whether the sample-size checks passed.

Validate test-specific requirements

Description

Aggregates validation logic tailored to the supplied test_name and returns a summary of any problems detected.

Usage

rvalidateTestRequirements(test_name, model, data, ...)

Arguments

Value

Validation result combining all relevant requirements.

Validate diagnostics against reference implementations

Description

Compares the package's test results with established implementations from other packages (currently lmtest and car). Differences exceeding the supplied tolerance are flagged.

Usage

rvalidate_against_reference(test_name, model, data, tolerance = 1e-06)

Arguments

Value

A list describing the comparison, including a status field with values "ok", "mismatch", or "skipped" when the reference package is unavailable.

Deprecated

Deprecated since 0.14.0: an internal helper of the package, which leaves the exports in 1.0.0. It still works and warns.

Safely fit a linear model

Description

Wraps stats::lm() in a tryCatch block that logs the error via ht_log() and optionally attempts recovery (e.g., dropping rows with missing values) before rethrowing. This helper is used internally by diagnostic tests so failures are easier to debug.

Usage

safe_lm(formula, data, ..., .recover = TRUE)

Arguments

Value

A fitted model object.

Safe variance calculation

Description

Computes the variance of x while guarding against near-zero values that could lead to division by zero in subsequent calculations.

Usage

safe_var(x)

Arguments

Value

A variance value with a minimum of .Machine$double.eps.

Colour scale for heteroscedasticity diagnostics

Description

Provides a discrete palette used across diagnostic comparisons.

Usage

scale_colour_hetero_diagnostic(...)

scale_fill_hetero_diagnostic(...)

Arguments

Value

A ggplot2 scale.

Examples

ggplot2::ggplot(mtcars, ggplot2::aes(wt, mpg, colour = factor(cyl))) +
  ggplot2::geom_point() +
  scale_colour_hetero_diagnostic()

Simulate heteroscedastic data

Description

Utilities to generate data with non-constant variance patterns. simulate_hetero simulates a linear regression model with heteroscedastic errors using a user-supplied variance function. The helper functions prefixed by sigma_ implement common variance patterns and can be passed to sigma_func. simulate_arch1 creates a simple ARCH(1) time series.

Usage

simulate_hetero(n, beta0, beta1, sigma_func, seed = NULL)

simulate_arch1(n, mu = 0, alpha0 = 0.5, alpha1 = 0.3, seed = NULL)

sigma_linear(x)

sigma_exponential(x)

sigma_group(x)

sigma_piecewise(x)

sigma_poly(x, a = 0.5, b = 0.1, c = 0.02)

sigma_sin(x, A = 1, B = 0.5, omega = 2 * pi/10, phi = 0)

sigma_multiplicative(x, mu_func, p = 1)

sigma_spatial(coords, gamma0 = 0.5, gamma1 = 0.1)

sigma_logistic(x, L = 2, k = 1, x0 = 5)

sigma_inverse(x, a = 1, b = 1)

sigma_power(x, a = 0.5, p = 0.5)

sigma_step(x, thr = 5, low = 0.5, high = 2)

sigma_u_shape(x, a = 0.1, b = 0.05, center = 5)

sigma_exp_decay(x, a = 2, b = 0.2)

sigma_gaussian_peak(x, base = 0.5, height = 1, mu = 5, sd = 1)

sigma_piecewise_linear(x, thr = 5, slope1 = 0.1, slope2 = 0.3)

Arguments