Package 'heteroTests'

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.

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.

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)

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)

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)

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)

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.

Description

Avoids recomputing expensive diagnostics by hashing their inputs and reusing stored results.

Usage

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

Arguments

Details

The cache key is created by hashing the test name, model coefficients, residuals, data, and additional arguments. When the key already exists in the cache environment the stored result is returned immediately; otherwise the diagnostic is executed and its output stored for future calls. Set use_cache = FALSE to force recomputation when the underlying model has changed.

Value

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

See Also

clearTestCache() resets the cache between analysis sessions.

Description

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

Usage

checkData(data)

Arguments

Value

Invisible data if valid, otherwise an error is thrown.

Description

Checks that model is of class lm or glm.

Usage

checkModel(model)

Arguments

Value

Invisible model if valid, otherwise an error is thrown.

Description

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.

Usage

checkModelEnhanced(model, data = NULL)

Arguments

Details

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)
checkModelEnhanced(fit)

Description

Removes memoised outputs stored by runHeteroTests when called with use_cache = TRUE. Useful for deterministic test runs or after mutating the underlying data/model so cached summaries should be invalidated.

Usage

clearAnalysisCache()

See Also

runHeteroTests, clearTestCache

Description

Removes all objects from the internal cache environment so that subsequent diagnostic calls recompute their statistics.

Usage

clearTestCache()

See Also

cachedTest() for invoking diagnostics with memoisation.

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

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

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

Arguments

Value

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

Examples

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

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)

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)

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)

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.

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

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.

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

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

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)

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)

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 generics::tidy(), generics::glance(), generics::augment() and ggplot2::autoplot().

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)
 test(d)

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.

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

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

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.

Description

Regresses squared residuals on their lags to detect ARCH effects.

Usage

performArchLMTest(model, 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.

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)

Description

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

Usage

performBartlettTest(model, data, 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.

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

Description

Test whether multiple groups have equal covariance matrices.

Usage

performBoxMTest(data, group)

Arguments

Value

An object of class htest.

Examples

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

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

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

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)
performBreuschPaganTest(model, data)

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

Value

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

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

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

Examples

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

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,
  studentized = TRUE,
  bootstrap = FALSE,
  B = 1000,
  ci_level = 0.95,
  parallel = FALSE
)

Arguments

Details

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

Value

An augmented htest object containing a robust_details element describing the additional diagnostics.

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 performStudentizedBPTest() for the standalone studentised variant.

Examples

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

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

Value

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

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)

Description

Brown-Forsythe test using medians instead of means.

Usage

performBrownForsytheTest(model, data, 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.

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

Description

The Cook-Weisberg (1983) score test with the fitted values as the sole variance regressor.

Usage

performCookWeisbergTest(model)

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.

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)

Description

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

Usage

performDavidianCarrollTest(model, 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.

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)

Description

Non-parametric test based on ranks.

Usage

performFlignerKilleenTest(model, data, group)

Arguments

Details

Absolute residuals are ranked after adjusting for group medians. The statistic 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.

References

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

Examples

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

Description

Regresses absolute residuals on a transformation of a suspected variable. With robust = TRUE the test is the version of Im (2000) and of Machado and Santos Silva (2000), which keeps its level when the errors are asymmetric.

Usage

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

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.

The default statistic has its nominal level when positive and negative errors are equally likely. Use robust = TRUE when the residuals are 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 are recorded in inst/validation/pass-a-size-power.csv for the default and in inst/validation/glejser-skewed-errors.csv for both.

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 the default statistic 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 the default. The two p-values agree as the sample grows.

In the same study the corrected statistic is within simulation error of the 5% level from 150 observations on. At 50 observations it rejects 5.7% to 6.6% of the time under skewed errors. Under symmetric errors $\hat{m}$ is close to zero and the two statistics differ little: the corrected one rejects about half a percentage point 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.

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)
 performGlejserTest(m, mtcars, "wt")

 # The version that keeps its level under asymmetric errors
 performGlejserTest(m, mtcars, "wt", robust = TRUE)

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

Description

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

Usage

performGQTest(model, data, 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 original model is re-estimated on the lower and upper segments. 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).

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.

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

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

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.

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.

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

Description

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

Usage

performHartleyFmaxTest(model, data, 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.

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

Description

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

Usage

performHarveyTest(model, 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’.

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

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

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)

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)

Description

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

Usage

performKoenkerTest(model, data)

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.

Value

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

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.

References

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

Examples

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

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

Value

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

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)

Description

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

Usage

performLeveneTest(model, data, 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.

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

Description

Applies the Ljung-Box test to squared residuals.

Usage

performMcLeodLiTest(model, 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.

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)

Description

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

Usage

performNCVTest(model, 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.

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.

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)
 performNCVTest(m)

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

Description

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

Usage

performOBrienTest(model, data, 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.

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

Description

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

Usage

performParkTest(model, data, 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.

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

Description

Checks for cross-sectional dependence in panel residuals.

Usage

performPesaranTest(model, data, 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.

Value

An object of class htest.

References

Pesaran, M. H. (2004). General diagnostic tests for cross section dependence in panels. Cambridge Working Papers in Economics.

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

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

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

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,
  order_by = NULL,
  B = 999,
  progress = interactive()
)

Arguments

Value

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

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)

Description

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

Usage

performRESETTest(model, 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.

Value

An object of class htest.

Examples

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

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.

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

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,
  listw,
  permutations = 499,
  zero.policy = NULL
)

Arguments

Value

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

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

Description

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

Usage

performSpearmanTest(model)

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.

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)

Description

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

Usage

performSpreadLevelTest(model)

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.

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)

Description

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.

Usage

performStudentizedBPTest(model, data)

Arguments

Details

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 implementation shares the validation helpers used across the package to ensure that: (i) the model and data satisfy minimum sample-size thresholds via rvalidateModelInputs() and rvalidateDataInputs(), (ii) missing values are handled by rhandleMissingValues(), 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.

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)
performStudentizedBPTest(mod, mtcars)

# Detect heteroscedasticity driven by a single regressor
set.seed(321)
x <- runif(180)
y <- 5 - 1.5 * x + rnorm(180, sd = 0.4 + 0.6 * x)
df <- data.frame(y, x)
performStudentizedBPTest(lm(y ~ x, data = df), df)

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, 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$. Under homoskedasticity $h$ is centred on the mean rank $(n + 1) / 2$ with variance $(n^2 - 1) / (6n)$, so

$Q = \frac{h - (n + 1) / 2}{\sqrt{(n^2 - 1) / (6n)}}$

is asymptotically standard normal. Variance that grows with order_by shifts weight onto the high ranks and drives $Q$ upwards.

The implementation relies on the package validation helpers to align the supplied data with the model residuals, 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.

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 asymptotic normal reference distribution.

Validation

The statistic and its null variance are checked against an independent reconstruction of Szroeter (1978) in ‘tests/testthat/test-pass-a-reference.R’. 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’.

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.

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

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.

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

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)

Description

Implements White's test on a fitted linear model.

Usage

performWhiteTest(model, data, 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 with degrees of freedom equal to the number of regressors in the auxiliary model.

Value

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

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

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.

Examples

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

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, B = 1000, parallel = FALSE)

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.

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.

The weights are held fixed in every replicate. When they were estimated from the same data, as by fitWLS(), the replicates leave that estimation out. In inst/validation/weighted-fits-size.R, on data whose standard deviation is $\exp(0.3 x)$, with the variance function correctly specified, normal errors, 150 observations and B = 499, the rejection rate at the 5% level was 0.064. In the same design with the true weights it was 0.054, and with the weights estimated again in each replicate, which this function does not do, 0.042. The section on weighted fits in performWildBootstrapTest() has the larger error of Koenker's statistic.

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 forks parallel::detectCores() - 1 workers, which Windows
# does not support. The call runs only where that is at most two workers,
# the most R CMD check --as-cran allows.
if (.Platform$OS.type != "windows" && isTRUE(parallel::detectCores() <= 3L)) {
  performWhiteTestBootstrap(mod, sim, B = 199, parallel = TRUE)
}

Description

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

Usage

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