Releases

The changes in each version, from the package's NEWS.md.

Development version

Packaging

• The logo shows in the README, and the R CMD check NOTE "Non-standard file/directory found at top level: 'assets'" is gone. The logo is in man/figures/logo.png, where R tooling looks for it, and assets/ is left out of the built package.

Documentation

• The documentation site, https://diogoribeiro7.github.io/heteroTests/, moves from pkgdown to Material for MkDocs and shows the package logo. Its API reference is generated from the help pages by tools/build-docs.R, in the sections the pkgdown index had, and tests/testthat/test-reference-index.R checks that every export is in one; the vignettes become its user guide. No function changes.

0.15.0

0.15.0 makes the interface match its description before 1.0.0. Every test takes (model, data, ...), with data = NULL in second position and performBoxMTest(data, group) the documented exception; performBPTestRobust() chooses the studentized statistic with studentize, as performHarveyTest() and lmtest::bptest() do; rrunAdvancedDiagnostics() is renamed runAdvancedDiagnostics(); and ?test_results documents what every test returns and the inputs every test refuses. Most of these changes break calls, not values: an old call form that 0.15.0 replaces still works in 0.15.x, warns with the call to use instead and, apart from the changes in bold, returns what it returned in 0.14.0, and it stops working in 1.0.0 (the table under Calls that changed gives each). The changes in bold change what calls return: every resampling function draws 999 replicates by default; the default battery of runHeteroTests() is White and Koenker; performGlejserTest() corrects for asymmetric errors by default; the tests refuse a glm that is not a least-squares fit, which most of them returned a value for, and so do the registry names "ncv" and "student_bp", deprecated in 0.14.0; they refuse a model with an NA coefficient and data that lack rows or variables of the fit, which some of them returned a value for, and performRESETTest() refuses a survey::svyglm() fit, as the other tests did already; the Goldfeld-Quandt, RESET, White bootstrap and wild bootstrap tests read the prior weights of a Gaussian glm, and the quantile regression test refuses a weighted one; the tests that read data refuse data whose rows they cannot match to the observations of the fit, and runHeteroTests() and runDiagnostics() stop on such data rather than fall back on another test; the robust White and Breusch-Pagan tests and the White bootstrap give the value on the data of the fit for rows they had renumbered; seeded resampling with parallel = TRUE is reproducible, and runHeteroTestsParallel() gives the same value with one worker as with several; the bootstrap tests estimate the weights of a fitWLS() fit again in each replicate; and the studentized Breusch-Pagan statistic reads the weights of a fitWLS() fit.

Resampling tests draw 999 replicates by default

Every resampling function draws 999 replicates by default, so its default p-value changes. A resampling p-value is (1 + #) / (B + 1), the number of replicates at least as large as the observed statistic, plus one, over the number of replicates plus one. A test at level alpha that rejects when it is at most alpha is exact only when alpha * (B + 1) is an integer: 999 is at the 5% and 1% levels, and 1000 is at neither. The defaults were a mix:

Function	Default before	Default now
performWildBootstrapTest()	B = 499	B = 999
performWhiteTestBootstrap()	B = 1000	B = 999
bootstrapTestStatistic() and rbootstrap_test_statistic()	B = 1000	B = 999
performWhiteTestRobust() and performBPTestRobust(), with bootstrap = TRUE	B = 1000	B = 999
runAdvancedDiagnostics() and rrunAdvancedDiagnostics()	bootstrap_B = 500	bootstrap_B = 999
performSpatialHeteroTest()	permutations = 499	permutations = 999
performRankPermutationTest()	B = 999	unchanged

With set.seed(1) before each call, the default p-value was and is:

Call	Before	Now
performWildBootstrapTest(lm(mpg ~ wt + qsec, data = mtcars))	0.2400	0.2330
performWhiteTestBootstrap(lm(dist ~ speed, data = cars))	0.1668	0.1670
bootstrapTestStatistic(performKoenkerTest, lm(dist ~ speed, data = cars), cars)$p_value	0.0799	0.0800
performWhiteTestRobust(lm(dist ~ speed, data = cars), bootstrap = TRUE)$p.value_bootstrap	0.1688	0.1690
performBPTestRobust(lm(dist ~ speed, data = cars), bootstrap = TRUE)$p.value_bootstrap	0.0799	0.0800
runAdvancedDiagnostics(lm(mpg ~ wt + qsec, data = mtcars), mtcars), White	0.0259	0.0290
runAdvancedDiagnostics(lm(mpg ~ wt + qsec, data = mtcars), mtcars), Breusch-Pagan	0.2335	0.2420
performSpatialHeteroTest(lm(mpg ~ wt + hp, data = mtcars), listw = lw)	0.5100	0.5080

where lw joins each car to its four nearest neighbours in coordinates drawn by set.seed(7); cbind(runif(32), runif(32)). runAdvancedDiagnostics() bootstraps both tests because mtcars has fewer than 50 rows. The registered tests "wild_bootstrap", "white_bootstrap" and "spatial_hetero", and so runHeteroTests() and the functions built on it, change with them. The B in the parameter of the two bootstrap tests, robust_details$bootstrap$B of the robust tests, which is recorded even when bootstrap = FALSE, and metadata$bootstrap_B of runAdvancedDiagnostics() report the new default. Pass the old value to reproduce a result of 0.14.0.

Of the recorded evidence, only the full sweep calls a resampling test with its default, performSpatialHeteroTest(). Its row was regenerated from the same seeds: under the null it rejected 0.045 of the time and rejects 0.050, and against the alternative 0.362 and 0.350 (inst/validation/full-sweep-size-power.csv). The other studies give B explicitly.

White and Koenker are the default battery

runHeteroTests() and the functions that share its default run White's test and Koenker's studentized Breusch-Pagan test, c("white", "koenker"), where they ran White's and the classical Breusch-Pagan test, c("white", "breusch_pagan"). The classical statistic divides by 2 sigma^4, the variance of the squared errors only when they are normal. In the full sweep, with 150 observations and 400 replications, it rejects a true null hypothesis 28.75% of the time at the 5% level when the errors have a t5 distribution (115 of 400 samples), against 4.75% (19 of 400) for Koenker's statistic; under Gaussian errors they reject 5.75% and 4.50%, and both reject every sample against sd = x^2 (inst/validation/full-sweep-size-power.csv). Pass A, at 5000 replications and 100 observations, agrees: 21.86% against 4.24% under t5 errors. What Koenker's statistic gives up is a little power when the errors are normal: in Pass A it rejects 72.70% of the time against sigma^2 = exp(0.4 x1), where the classical test rejects 77.44% (inst/validation/pass-a-size-power.csv).

The default changes in runHeteroTests(), runDiagnostics(), summary() and the deprecated test() on a HeteroDiagnostic object, compareModelDiagnostics(), compareTestResults(), runSurveyHeteroTests() and the tests the dashboard selects at start, and so in generateDiagnosticReport() and generateHeteroRecommendations(), which call runHeteroTests() with its default. The classical test still runs by name, tests = "breusch_pagan", or c("white", "breusch_pagan") for the old battery. The second element of a default suite is now called koenker, so code that reads $breusch_pagan from it gets NULL. For lm(mpg ~ wt + qsec, data = mtcars), runHeteroTests(m, mtcars) returned White's test (11.82 on 5 degrees of freedom, p = 0.0373) and the classical Breusch-Pagan test (3.135 on 2, p = 0.2086), and returns White's test and Koenker's (3.086 on 2, p = 0.2138).

Two other functions follow the new default. suggestRemediation() suggests variance modelling when the Koenker test rejects as well as when the classical test does; it read only "breusch_pagan", so it would never have suggested it for a default suite. suggestDiagnosticsForProfile() recommends Koenker's test, where it recommended the classical one, beside White's as the baseline for every data profile. runAdvancedDiagnostics() is unchanged: with its default auto_enhance = TRUE its Breusch-Pagan entry is already the studentized statistic, but auto_enhance = FALSE turns the studentization off with the bootstrap and runs the classical test, which is listed under 0.16.0. None of the recorded evidence calls these functions with their default battery, so no recorded number moves.

The Glejser test corrects for asymmetric errors by default

performGlejserTest() computes the statistic of Im (2000) and of Machado and Santos Silva (2000) by default, robust = TRUE, where it computed Glejser's (1969), so the statistic and p-value it reports by default change. Glejser's statistic regresses the absolute residuals on the variable tested, and has its nominal level only when positive and negative errors are equally likely. In inst/validation/glejser-skewed-errors.csv, at 5000 replications, 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. The corrected statistic, added in 0.13.0, subtracts the mean sign of the residuals times each residual from its absolute value. Under the same skewed errors it rejects 4.6% to 5.7% of the time from 150 observations on, and all 27 of its null cells from 150 observations, with symmetric errors too, are inside the release-gate interval, [0.0421, 0.0579]. At 50 observations it is not: under skewed errors it rejects 5.7% to 6.6% of the time, and four of those six cells are above the interval, at 0.0580, 0.0592, 0.0626 and 0.0658, where Koenker's statistic on the same regressor is at 0.0586 to 0.0714. Under symmetric errors the two statistics differ little: the corrected one rejects 0.5 to 0.7 percentage points less often at 50 observations, and against sigma^2 = exp(gamma x1) with Gaussian errors it has 1.1 and 1.9 points less power at 50 observations, and 0.8 and 0.1 at 150.

The method of the result says which statistic was computed, and by default now reads "Glejser test for heteroscedasticity (robust to asymmetric errors)". Pass robust = FALSE for Glejser's statistic: it returns what the default returned in 0.14.0, checked with identical() against 0.14.0. The default was and is:

Call	0.14.0	Now
performGlejserTest(lm(mpg ~ wt + qsec, data = mtcars), mtcars, "wt")	t = -0.450 (p = 0.656)	t = -0.395 (p = 0.695)
performGlejserTest(lm(dist ~ speed, data = cars), cars, "speed")	t = 2.059 (p = 0.045)	t = 1.943 (p = 0.058)

performGlejserTest() is not a registered test, so runHeteroTests(), the dashboard and the recommendations do not call it, and no other function of the package changes with it. Of the recorded evidence, the full sweep and the known-weights block of the weighted fits call it with its default. Their two rows were regenerated from the same seeds, and with robust = FALSE the same runs reproduce the rows of 0.14.0 byte for byte: in the full sweep the size moves from 0.0600 to 0.0575, the size under t5 errors from 0.0575 to 0.0550, and the power stays 1.000 (inst/validation/full-sweep-size-power.csv); with known weights the rejection rate moves from 0.0475 to 0.0450 (inst/validation/weighted-fits-size.csv). Pass A and the Glejser study name Glejser's statistic, and now give robust = FALSE; their numbers do not change.

Inputs every test refuses

Calls that returned a value for a glm that is not a least-squares fit, for a model with an NA coefficient, or for data that cannot be matched to the fit now stop with an error. tests/testthat/test-input-validation.R gives every test these inputs, and each is refused with one message (?test_results). Before 0.15.0 these calls returned a value:

Input	Tests that returned a value and now stop	Message
A glm of a family other than gaussian, or with a link other than the identity, such as glm(am ~ wt, family = binomial)	performWhiteTest(), performWhiteTestRobust(), performWhiteTestBootstrap(), performBPTest() and performBreuschPaganTest(), performBPTestStreaming(), performBPTestRobust(), performKoenkerTest(), performKoenkerTestStreaming(), performStudentizedBPTest(), performGQTest(), performParkTest(), performSzroeterTest(), performWildBootstrapTest(), performQuantileRegressionTest(), performRankPermutationTest(), performHighDimensionalTest(), performSpatialHeteroTest(), performSpearmanTest(), performLeveneTest(), performGlejserTest(), performBrownForsytheTest(), performFlignerKilleenTest(), performBartlettTest(), performHarveyTest(), performArchLMTest(), performMcLeodLiTest(), performHartleyFmaxTest(), performOBrienTest(), performCookWeisbergTest(), performNCVTest(), performSpreadLevelTest(), performPesaranTest(), performBPRandomEffectsTest(), performDavidianCarrollTest(), performRESETTest(): 36 names, all but the three tests that already refused every glm	"<test>() needs a model fitted by stats::lm(): a glm with the binomial family and the logit link is not a least-squares fit."
A model with an NA coefficient, such as lm(mpg ~ wt + I(2 * wt), data = mtcars)	performSpreadLevelTest(), performDavidianCarrollTest(), performPesaranTest(), performBPRandomEffectsTest(), performRESETTest()	"Model coefficients must be finite before running ..."
data that lack rows the model was fitted to	performVarianceFormTest()	"... requires data to contain the rows used to fit the model. Missing rows: ..."
data that lack a variable of the model	performVarianceFormTest(), performPesaranTest(), performBPRandomEffectsTest()	"The following variables are missing from input data: ..."
A survey::svyglm() fit, of any family	performRESETTest(), the only test that did not refuse one	"The model is a survey::svyglm() fit, whose residuals are scaled by the sampling weights. ..."

The values these calls returned tested nothing the help pages describe. A binomial or Poisson fit has residuals whose variance depends on the mean by construction: performKoenkerTest(glm(am ~ wt + hp, data = mtcars, family = binomial), mtcars) returned 0.963 (p = 0.618), and performRESETTest() on the same fit F = 86.1 (p < 0.0001). The residuals of a survey::svyglm() fit are scaled by the sampling weights, which are not precision weights, and every other test refused one already. A model with an aliased term was computed by the five tests as if the term were not there, and refused by the others. The variance-form test returned a value for data that lack rows of the fit or a variable of the model, and the Pesaran and random-effects tests for data that lack a variable, which they did not look for. performPesaranTest() on data that lack rows stopped with "number of items to replace is not a multiple of replacement length" in the designs tried, and now names the rows, as the others do. In test-input-validation.R, on the design it uses, all 36 names returned a value for a binomial and for a Poisson glm before 0.15.0.

A glm with the Gaussian family and the identity link is a least-squares fit and is still accepted. Without prior weights, and with positive ones, every test that accepts it gives the result of the lm in that design to 1e-8, the five tests below since they read its weights. The exception is the Fligner-Killeen statistic, which moves by about 1e-6 of its value: with an even number of residuals the two in the middle tie at exactly half their difference from the median, and rounding, which differs between the two fits by about 1e-15, decides which ranks first. A Gaussian glm with a zero prior weight is not yet refused as the same lm is: the tests that read the Pearson residuals, such as the Koenker, White, Glejser, Spearman, GQ and wild bootstrap tests, return a value for it, and performWhiteTestBootstrap() stops with "Model residuals must be finite". It is recorded under 0.16.0. performWhiteTestStreaming(), performHarrisonMcCabeTest() and performVarianceFormTest() refuse every glm, as before.

The registry names deprecated in 0.14.0, "ncv" and "student_bp", and the "student_bp" entry of test_factory refuse such a glm too, with the message of performNCVTest() and performStudentizedBPTest(), and so does rvalidate_against_reference("studentized_bp", ...). They call the implementations behind those two functions rather than the functions, and returned a value for it through runHeteroTests(tests = ...), runDiagnostics(), cachedTest(), runHeteroTestsParallel() and test_factory$run_test(). runHeteroTests() falls back on "ncv" when White's, Koenker's, the classical, the quantile regression or the rank permutation test stops, so runHeteroTests(<binomial glm>, data) returned the Cook-Weisberg statistic under the names white and koenker, tagged as a fallback, where 0.14.0 returned White's and the classical test. Decided on 2026-10-16: it now stops with the refusal of performWhiteTest(), as runHeteroTests(<binomial glm>, data, tests = "cook_weisberg") stops with that of performCookWeisbergTest(), and with grouped data it returns each test as failed, with its refusal. runDiagnostics() stops at the same point, where it stopped later, at performRESETTest(). No name holds a fallback computed on such a fit. On an lm, a Gaussian glm and a weighted lm, the two names, the factory entry and the default battery return values identical() to those they returned before.

performGQTest(), performRESETTest(), performWhiteTestBootstrap() and performWildBootstrapTest() now read the prior weights of a Gaussian glm, and performQuantileRegressionTest() refuses a weighted one, as it refuses a weighted lm. These five read the weights of an lm, to refit with them or to refuse a weighted fit, and found none on a glm, so they took a weighted Gaussian glm for an unweighted fit. Decided on 2026-10-15: they read its prior weights, and give what they give for the same model fitted by lm(..., weights = w). test-input-validation.R asserts it to 1e-8, for the statistic and the p-value, on four seeded designs, the resampling tests from the same seed and the same B; the largest differences there are 6e-14 in a statistic and 1e-14 in a p-value. On glm(dist ~ speed, data = cars, weights = 1 / speed), the two bootstraps after set.seed(1):

Call	0.14.0	Now, as for the lm
performGQTest(order_by = "speed")	5.42 (p = 0.0004)	2.65 (p = 0.023)
performRESETTest()	F = -21.5 (p = 1.000)	F = 1.37 (p = 0.264)
performWhiteTestBootstrap(B = 999)	p = 0.695	p = 0.701
performWildBootstrapTest(B = 999)	p = 0.378	p = 0.373
performQuantileRegressionTest()	2.05 (p = 0.155)	Error: "performQuantileRegressionTest() does not support weighted fits: ...", the message for the weighted lm

The observed statistics of the two bootstraps, 0.819 and 0.766, were already those of the lm; their replicates were drawn as for an unweighted fit. On glm(mpg ~ wt + hp, data = mtcars, weights = 1 / wt) RESET returned F = -6.13 (p = 1.000), and returns 6.11 (p = 0.0065), the value of the lm.

No recorded evidence fits a glm, a model with an NA coefficient or data that do not match the fit: the rows of the full sweep for the six tests whose checks changed besides the glm one, the Pesaran, random-effects, variance-form, RESET, spread-level and Davidian-Carroll tests, run again from their seeds, are identical.

Data in another order than the fit

Calls that returned wrong values on data holding the rows of the fit in another order now stop with an error. The tests match the rows of data to the residuals by row name. A tibble, a data.table, or a data frame whose row names were reset, as dplyr::arrange() resets them, has row names that only number its rows, so such data were matched by position: in another order than the fit, or with rows in front of those of the fit, each residual was read with another observation's row, without an error. With

set.seed(1); n <- 120
d <- data.frame(x = runif(n, 1, 5), z = rnorm(n))
d$y <- 1 + d$x + d$z + rnorm(n, sd = d$x)
m <- lm(y ~ x + z, d)

performKoenkerTest(m, d) gives 17.56, and performKoenkerTest(m, dplyr::arrange(d, x)) gave 0.114. It now stops with "performKoenkerTest() cannot match the rows of data to the observations of the model: data has no row names of its own (a tibble, a data.table, or a data frame whose row names were reset, as by dplyr::arrange()), so its rows are matched to the observations by position, and in that order its values are not those of the fit: y, x, z differ from the model frame. Pass the rows of the fit in the order of the fit, or a data frame with its row names."

Before a test reads data, it takes the rows that the row names pick for the observations of the fit and compares them, row by row, with every column of the model frame of the fit. A column that data holds under its own name, such as y or x, is compared value by value, each number to 1e-8 of its own size, and dates and times exactly. Any other column is computed from data as the fit computed it: a term such as log(y), log(x), poly(x, 2), ns(x, 3) or offset(o) through the terms of the fit, and weights or an offset given as arguments from the expressions of the call. It is computed on all the rows of data, as model.frame() computed it before dropping rows, or else on the matched rows, and compared value by value too. So rows that agree on the response but not on a transformed regressor, a weight or an offset are told apart, as when the data are sorted on a 0/1 response. A term that depends on the rows it is computed from, such as I((x - mean(x))^2), I(x > median(x)) or rank(x), comes out otherwise on the rows a fit with subset = or na.omit kept, which are the right rows; it is not compared, and neither is a column that names a variable data does not hold or whose computation fails. A column computed from a vector of the workspace, such as weights = wts, is read by no test from data and needs nothing more. For the others, the columns that are compared must tell the rows of the fit apart.

The guard refuses data it cannot confirm, even when they are the data of the fit. It does so in two cases. One is when no column of the model frame can be compared, as with lm(d$y ~ d$x). The other is when the compared columns tie on rows of the fit that differ in a term that depends on the rows, or in a term of variables of data that cannot be computed. That happens with a discrete response and a row-dependent regressor on the rows an na.omit fit kept, as in lm(score ~ g + I(x > median(x)) + z2), and with y01 ~ myf(x) after rm(myf) or ns() after detach("package:splines"). Pass the full data the model was fitted on, with the rows it dropped, or its variables as columns of data. Data whose rows differ are refused, and so are repeated row names. Data with row names of their own, such as d[sample(n), ], are still matched by name, the tibble or data.table the model was fitted on is accepted, and so are rows after those of the fit. On eighteen fits (unweighted, weighted, with an offset given as an argument or as a term, log(y), na.omit with three rows dropped, a Gaussian glm, a fitWLS() fit, log(x) + poly(x, 2), model = FALSE, a discrete response with weights or an offset from the workspace, and seven with a term computed from the rows under na.omit or subset =), every exported test on the data of the fit, with data = NULL, on the same rows shuffled with their row names, on the rows the fit kept, and on the tibble and the data.table the model was fitted on returns what it returned, identical(), except the three tests below on the shuffled rows, which returned wrong values, and the robust and bootstrap tests that stopped with an error on the tibble of a subset = fit and now return a value (below); some errors changed their message.

The guard covers the 23 tests that returned wrong values on such data: performWhiteTest(), performWhiteTestStreaming(), performWhiteTestRobust(), performWhiteTestBootstrap(), performBPTest() and performBreuschPaganTest(), performBPTestStreaming(), performKoenkerTest(), performKoenkerTestStreaming(), performRankPermutationTest(), performHighDimensionalTest(), performParkTest(), performSzroeterTest(), performHarrisonMcCabeTest(), performGlejserTest(), performHartleyFmaxTest(), the five group tests (performLeveneTest(), performBrownForsytheTest(), performFlignerKilleenTest(), performBartlettTest(), performOBrienTest()), performPesaranTest() and performBPRandomEffectsTest(), and also performCookWeisbergTest() with var_formula, which read the variables of the variance formula from the same rows. Six tests that start the same way, performGQTest(), performStudentizedBPTest(), performBPTestRobust(), performWildBootstrapTest(), performSpatialHeteroTest() and performQuantileRegressionTest(), returned the right value on reordered numbered data, the quantile test to rounding (a relative difference of 3e-15 on the data above), because they take everything but an ordering from the fit or refit on the rows; they refuse it too, so that every test that reads data treats it alike, and so do the registry names "szroeter", "wild_bootstrap", "quantile_regression" and "student_bp" in runHeteroTests(), which returned the right value too. Exempting them would leave the cases with rows in front of those of the fit, on which they were wrong. The tests that read nothing from data, such as performSpearmanTest() or performRESETTest(), are unchanged. bootstrapTestStatistic() refused such data already (0.13.0); with one of these tests as test_function it now stops with the test's message, before its own check.

A model fitted with model = FALSE keeps no model frame, and the tests rebuild it from the data the fit refers to. When those data changed after the fit, the rebuilt frame held other rows, and the tests computed on it: after m <- lm(y ~ x + z, d, model = FALSE) and d <- dplyr::arrange(d, x), performKoenkerTest(m), performKoenkerTest(m, d) and Koenker's test in runHeteroTests(m) gave 0.114 instead of 17.56. The response of the rebuilt frame is now compared with that of the fit, and its design, times the coefficients, plus its offset, with the fitted values, and such a frame is refused ("cannot use the model frame of model").

A tibble that has more rows than the fit and incomplete rows was renumbered when a test dropped those rows, and its rows then matched others: on a fit to the first 100 rows of a tibble of 120, three of them incomplete, Koenker's statistic was 0.0165 where the rows of the fit give 7.93. It is refused, with the advice to pass as.data.frame(data).

performWhiteTestRobust(), performWhiteTestBootstrap() and performBPTestRobust() give the value on the data of the fit where they gave another one. They call White's or the Breusch-Pagan test again on the rows they matched to the fit, and a tibble they had subset numbered those rows afresh, so the inner test matched them to other observations. With a fit to set.seed(4); ds <- d[sample(n), ] and data = tibble::as_tibble(d), valid data, performWhiteTestRobust() and the observed statistic of performWhiteTestBootstrap() were 3.09 and performBPTestRobust(studentize = FALSE) 0.140, where the rows of the fit give 21.79 and 37.06. They now pass the rows under the row names of the fit, and return those values.

performPesaranTest(), performBPRandomEffectsTest() and performScatterDiagnostic() read data by position, the i-th row with the i-th residual, so they gave wrong values on data with row names of their own in another order too. With d$id <- rep(1:12, each = 10), d$time <- rep(1:10, 12) and set.seed(2); dp <- d[sample(n), ], the Pesaran statistic was -0.987 on d and -0.392 on dp, the random-effects statistic 1.19 and 0.506, and performScatterDiagnostic(m, ., "x") 0.468 and -0.080. They refuse such data ("the rows of data, matched by name, are not in the order of the fit"), data whose row names are not those of the fit, which cannot show their order, and data with another number of rows than the fit, with the message the random-effects test gave already ("requires data with 117 observations to match the fitted model, got 120"), where the Pesaran test stopped with "arguments imply differing number of rows" and the scatter diagnostic with "incompatible dimensions".

runHeteroTests() and runDiagnostics() stop with this error, and do not fall back on another test for it. On dplyr::arrange(d, x) above, runHeteroTests(m, ...) returned White's statistic 1.50 and Koenker's 0.114, both computed on misaligned rows (21.79 and 17.56 on d). When one of their tests stops, the runners run a substitute under its name, the last of them the Cook-Weisberg test, which reads nothing from data; for this refusal there is none, so no name can hold a test computed on, or regardless of, such data. The error keeps the class of the refusal, heteroTests_data_rows_error.

The guard is a stopgap. 0.16.0 replaces it with the model context of the roadmap, which takes the response, design and weights from the fit and matches the rows of data by name or by verified position, so that data in any order give the value on the data of the fit.

Seeded resampling in parallel is reproducible

bootstrapTestStatistic(), performWhiteTestBootstrap(), the robust tests with bootstrap = TRUE and runHeteroTestsParallel() return the same value for the same seed with parallel = TRUE, and a different value from the one they returned in parallel before. The replicates computed in worker processes drew from generators that set.seed() in the session did not reach, so a seeded call returned another p-value each time it ran. Replicate i, or test i in runHeteroTestsParallel(), now draws from stream i of a L'Ecuyer-CMRG generator (parallel::nextRNGStream()) seeded by one draw from the session's generator, on the forked path (parallel::mclapply()) and on the socket cluster alike, so the value depends on the seed and not on the number of workers or the path. With set.seed(1), bootstrapTestStatistic(performKoenkerTest, lm(dist ~ speed, data = cars), cars, parallel = TRUE, n_cores = 2)$p_value gave 0.0740 and then 0.0760 on two runs, and gives 0.0830 on every run.

A serial call draws every replicate from the session's generator, as before, and returns what it returned: 0.0800 for that call with parallel = FALSE. The parallel and the serial value differ because their draws differ; both are bootstraps of the same data, and tests/testthat/test-parallel-reproducibility.R checks that their replicates have the same distribution. The help pages say this under Parallel replicates.

runHeteroTestsParallel() with one worker returns a different value for a resampling test. With n_cores = 1, and when the default comes to one worker, for one test or on a machine with two cores or fewer, it runs the tests in the session, where they drew from the session's generator. Test i draws from stream i there too now, so the value is the one the workers give, whatever n_cores. With set.seed(99), runHeteroTestsParallel(lm(mpg ~ wt + qsec, data = mtcars), mtcars, tests = c("koenker", "wild_bootstrap"), n_cores = 1) gave a wild bootstrap p-value of 0.2180 in 0.14.0, that of performWildBootstrapTest() called from the same seed with the 499 replicates of then (0.2120 with the 999 of now), and gives 0.2110, as with n_cores = 2 or 3.

performWhiteTestBootstrap(parallel = TRUE) forked its workers with parallel::mclapply() on every platform, so on Windows, which cannot fork, it stopped with 'mc.cores' > 1 is not supported on Windows whenever the machine had more than two cores. It runs its replicates on a socket cluster there, as bootstrapTestStatistic() does, and gains n_cores, the number of workers, which the result does not depend on: with set.seed(1), performWhiteTestBootstrap(lm(dist ~ speed, data = cars), parallel = TRUE, n_cores = 2) gives 0.1750 on every run, against 0.1670 serially.

No recorded evidence runs a resampling test with parallel = TRUE.

Bootstrap tests estimate the weights of a fitWLS() fit again

On a fit returned by fitWLS(), performWildBootstrapTest(), performWhiteTestBootstrap() and bootstrapTestStatistic() estimate the weights again in each replicate by default, so the p-values they return on such a fit change. Since 0.13.0 they held the weights fitWLS() estimated fixed in every replicate, as if they were known. The replicates then leave out the estimation that the observed statistic went through, and the tests rejected a true null hypothesis too often, more so with more observations. In inst/validation/weighted-fits-size.R, the data have standard deviation exp(0.3 x) and fitWLS() fits its default exponential variance function, which is correctly specified. With normal errors, 2000 replications and B = 499, the rejection rates at the 5% level are:

Test	n	Weights held fixed, before	Estimated again, now	True weights
performWildBootstrapTest()	150	0.0865	0.0605	0.0540
performWildBootstrapTest()	600	0.1050	0.0515	0.0480
performWhiteTestBootstrap()	150	0.0640	0.0415	0.0540
bootstrapTestStatistic() with performKoenkerTest()	150	0.0825	0.0450	0.0455

(inst/validation/weighted-fits-size.csv, the rows whose weights are "estimated, held fixed", "estimated, re-estimated" and "true"). At 2000 replications the 99% interval around 0.05 is [0.037, 0.063], and all four rates with the weights estimated again are inside it. The release gate is set at 5000 replications, and against its interval, [0.042, 0.058], the wild bootstrap with 150 observations, 0.0605, is outside, above it, and the White bootstrap, 0.0415, just below it. Against [0.041, 0.057], the interval around 0.049, the rate of an exact test with B = 999, the wild bootstrap with 150 observations is outside and the other three are inside. No study of the new default at 5000 replications with B = 999 has been run; ROADMAP.md lists it under 0.16.0.

Those four rows were measured in 0.12.0 by repeating the resampling outside the package, and were not run again. The package reproduces the repetition on 200 of the 2000 samples of each: the same p-values, the same replicates and the same random draws (inst/validation/weighted-fits-refit-reproduction.csv).

Each replicate is refitted as fitWLS() fitted the model: by least squares, as the model fitWLS() was given, with its weights if it had any. The variance function fitWLS() recorded is then fitted again to the residuals, with the same form and on the regressors fitWLS() evaluated. Last, the replicate is refitted with the weights that gives. 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, as for a refit that fails; a warning gives the number that failed. fitWLS() records the weights of a weighted model it refits, as initial_weights in its "variance_function" attribute, for this. A fit made by an earlier version has no initial_weights, and its model frame, which still holds those weights, is read instead.

The new argument reestimate_weights chooses. NULL, the default, estimates the weights again on a fitWLS() fit and holds them fixed on any other model. reestimate_weights = FALSE holds them fixed and returns what 0.14.0 returned. TRUE on a model that fitWLS() did not return is an error. An lm, a weighted lm and a Gaussian glm carry no variance function, and they return what they returned. Two kinds of fitWLS() fit keep the values of 0.14.0 under the default. One on which fitWLS() fell back to equal weights estimated no weights. The other has a variance regressor computed from the response. Each replicate regenerates the response and holds the design fixed, so such a regressor is one that uses a variable of the response which the right-hand side of the model does not use, or one that is the response term itself, as ~ . gives. Examples are ~ y or ~ . for y ~ x + z, and ~ y or ~ x for I(y / x) ~ z. Its weights cannot be estimated again from the regressors fitWLS() recorded, so the default holds them fixed and warns. The rule is conservative for ~ x with I(y / x) ~ z: giving the response as a column, r ~ z, says that x stays. On both kinds of fit, TRUE is an error. In bootstrapTestStatistic() and rbootstrap_test_statistic() the argument comes after ..., so it has to be named in full.

With set.seed(1) before each call, for fit <- fitWLS(lm(dist ~ speed, data = cars)):

Call	Before, and with reestimate_weights = FALSE	Now
performWildBootstrapTest(fit, cars)	0.6730	0.6990
performWhiteTestBootstrap(fit, cars)	0.7640	0.7290
bootstrapTestStatistic(performKoenkerTest, fit, cars)$p_value	0.7020	0.6750

The functions that call these with their default change with them on a fitWLS() fit. These are the registered tests "wild_bootstrap" and "white_bootstrap", which runHeteroTests(), runHeteroTestsParallel() and the functions built on them run by name, and performWhiteTestRobust() and performBPTestRobust() with bootstrap = TRUE, and so runAdvancedDiagnostics() when it bootstraps. They have no argument of their own; to hold the weights fixed, call the bootstrap function directly with reestimate_weights = FALSE.

Of the recorded evidence, only the rows "estimated, held fixed" of weighted-fits-size.R call these functions on a fitWLS() fit. The script passes reestimate_weights = FALSE there, so they do not move, and it calls the package for the rows "estimated, re-estimated".

The studentized Breusch-Pagan statistic on a fitWLS() fit

On a fit returned by fitWLS(), performStudentizedBPTest(), performBPTestRobust() with its default studentize = TRUE, and the registered "student_bp" computed the statistic of the unweighted fit. They now compute that of the weighted fit, which is Koenker's. They read the weights from the model frame, and the model frame of a fitWLS() fit is that of the model fitWLS() refitted. It holds that model's weights, none for an unweighted model, and not the estimated ones. On a correctly specified fit the statistic therefore reported the heteroscedasticity that the weights remove. In the design of inst/validation/weighted-fits-size.R, with set.seed(s) for s from 1 to 5, x and z uniform on (1, 5), y = 1 + 2 x + 0.5 z + exp(0.3 x) e with normal e, 150 observations and fit <- fitWLS(lm(y ~ x + z)), performStudentizedBPTest(fit) returned p = 0.0019, 0.0015, 4.2e-07, 0.00097 and 0.013. It now returns 0.50, 0.86, 0.22, 0.56 and 0.53, the p-values of performKoenkerTest(). Since 0.15.0 the three read the weights of the fit, as the other tests do (decided on 2026-10-16).

The values change on a fitWLS() fit only:

• the statistic and p-value of performStudentizedBPTest();
• the statistic, the asymptotic p-value and the bootstrap p-value of performBPTestRobust(), and so the Breusch-Pagan entry of runAdvancedDiagnostics();
• "student_bp".

The model frame of a fitWLS() fit is unchanged. Nothing else changes. A corpus of seeded calls was run on the development version before and after the change: the corpus of the review of #129, extended to every test on the models that fitWLS() did not return.

• On an lm, a weighted lm, an lm with an offset given either way or with weights and an offset, a Gaussian glm with and without weights, and a model of log(y), all 472 calls return what they returned, with the same warnings and random draws. 464 are identical(). In the other 8 the result holds the glm, whose family closures compare equal only up to their environments.
• On the six fitWLS() fits, 228 of the 234 calls of the 39 other tests and runners return what they returned: 226 identical(), and 2 up to the initial_weights that a fit of a weighted model now records. The six that change are performStudentizedBPTest().

The statistic now equals Koenker's, so on a fitWLS() fit it has Koenker's level, which was already documented. The weights are taken as known, and with the variance function correctly specified and normal errors it rejects 0.077 of the time at the 5% level with 150 observations and 0.110 with 600 (inst/validation/weighted-fits-size.csv, performKoenkerTest(), weights "estimated"). The bootstrap of performBPTestRobust() now compares the statistic with replicates computed the same way. In an exploratory run by the reviewer of #129, not recorded in inst/validation/ (20 samples of a correctly specified design, B = 99), the 0.14.0 path rejected all 20 and the new default 1 of 20, with a mean p-value of 0.498. The statistic still drops an offset, as before; that is listed for 0.16.0 in ROADMAP.md. No recorded evidence calls these on a fitWLS() fit: weighted-fits-size.R runs performStudentizedBPTest() and performBPTestRobust() on fits with known weights only.

Calls that changed

Old form	New form	Old form still works
performArchLMTest(model, lags)	performArchLMTest(model, data = NULL, lags = 1)	In 0.15.x, with a warning, when lags is a number given second by position; not in 1.0.0
performMcLeodLiTest(model, lags)	performMcLeodLiTest(model, data = NULL, lags = 10)	As above
performDavidianCarrollTest(model, degree)	performDavidianCarrollTest(model, data = NULL, degree = 2)	In 0.15.x, with a warning, when degree is a number given second by position; not in 1.0.0
performCookWeisbergTest(model, var_formula)	performCookWeisbergTest(model, data = NULL, var_formula = NULL)	In 0.15.x, with a warning, when var_formula is a formula given second by position; not in 1.0.0
performNCVTest(model, var_formula) (deprecated)	performNCVTest(model, data = NULL, var_formula = NULL)	As above, with the deprecation warning of 0.14.0 as well
performHarveyTest(model, auxiliary, studentize)	performHarveyTest(model, data = NULL, auxiliary, studentize)	In 0.15.x, with a warning, when auxiliary is a string given second by position, alone or followed by studentize; not in 1.0.0
performRESETTest(model, power)	performRESETTest(model, data = NULL, power = 2:3)	In 0.15.x, with a warning, when power is a number given second by position; not in 1.0.0
performSpearmanTest(model)	performSpearmanTest(model, data = NULL)	Yes: data is optional
performSpreadLevelTest(model)	performSpreadLevelTest(model, data = NULL)	Yes: data is optional
performWhiteTest(model, data, ...) and the 28 other tests that took data second	performWhiteTest(model, data = NULL, ...)	Yes: data is now optional
performBPTestRobust(model, data, studentized = TRUE)	performBPTestRobust(model, data, studentize = TRUE)	In 0.15.x, with a warning; not in 1.0.0
rrunAdvancedDiagnostics(model, data, ...)	runAdvancedDiagnostics(model, data, ...)	In 0.15.x, with a warning; the old name goes in 1.0.0

Calls that name their arguments, such as performArchLMTest(model, lags = 3), are unchanged and do not warn.

The (model, data, ...) convention

Every test now takes (model, data, ...), with data = NULL in second position. The nine tests above that took no data, which read everything they need from the model, use the model frame when it is omitted, as performVarianceFormTest() and fitWLS() have done since 0.12.0. Giving them the data does not change the result, except that performCookWeisbergTest(model, data, var_formula) and performNCVTest() evaluate var_formula in data, matched to the residuals by row name, so the variance model can use a variable the mean model does not; with data = NULL it is evaluated in the model frame, as before. performBoxMTest(data, group) is the documented exception: it tests the covariance matrices of the data, not the residuals of a model. tests/testthat/test-model-data-convention.R checks the signatures, the warnings and that each old and new call returns the value of the old one. The registry entries "cook_weisberg" and "spread_level" pass the data on.

The other 29 tests, which read variables from data and required it, can now be called without it: data = NULL stands for the model frame of model, so performWhiteTest(model) returns what performWhiteTest(model, stats::model.frame(model)) returns, which is what runHeteroTests(model) passes them. A call that gives data returns what it returned. The model frame holds the terms of the model rather than its variables (log(x) rather than x), so a model that transforms a variable or has an offset() term still needs data, as does a test on a column that is not a variable of the model, such as the group of performLeveneTest(). On a weighted fit whose weights are a column of the data (weights = w), the eight tests that rebuild the design from data, performWhiteTest(), performWhiteTestRobust(), performWhiteTestBootstrap(), performBPTest() and performBreuschPaganTest(), performKoenkerTest(), performRankPermutationTest() and performHighDimensionalTest(), stop with "object 'w' not found" without it, as they do when given the model frame, and so do they on a fit with offset = o. runHeteroTests() and runDiagnostics() then return the Cook-Weisberg test, which they fall back on when a test stops, under the names white and koenker (in 0.14.0 under white and breusch_pagan), so pass data for such a fit. One call returns another value: performRESETTest(model, NULL), which read NULL as power and returned NaN, now computes the test with the default powers, and a data frame given second, which was read as one power per column, is now read as data. In a seeded comparison of 30 tests in 39 settings on five designs, among them a weighted fit and a fit with NA rows, 565 of the 571 calls that returned a value before return an identical() value, the other six being such performRESETTest() calls, and a call without data is identical() to one given the model frame in all 195 cases.

studentize everywhere

The argument that chooses the studentized statistic is called studentize throughout, as in lmtest::bptest() and performHarveyTest(). performBPTestRobust() called it studentized; the old name is still accepted in 0.15.x, does what studentize does and warns that it stops working in 1.0.0. The field robust_details$studentized of the result keeps its name, so the returned object does not change.

runAdvancedDiagnostics()

rrunAdvancedDiagnostics() is renamed runAdvancedDiagnostics(), the last of the r-prefixed exports 0.14.0 reviewed, so that it follows the camelCase names of the other runners. The old name still works, returns what it returned and warns with .Deprecated() that it goes in 1.0.0, with the names deprecated in 0.14.0.

The return contract

A new help page, ?test_results, documents every field a perform*Test() returns: the five every test returns (statistic, parameter, p.value, method, data.name), the optional alternative and estimate, and the fields a few tests add, such as robust_details or bootstrap. tests/testthat/test-return-contract.R runs every exported test, and the settings that add fields, and fails on a field the page does not document or one that breaks it. Six tests needed a change to comply, and no number changes in any of them:

• the parameter of performSpearmanTest(), performSpreadLevelTest(), performBartlettTest(), performFlignerKilleenTest() and performBPRandomEffectsTest() gains its name, df; it had none, so print() showed t = -0.41855, = 30, p-value = 0.6785;
• performBoxMTest() gains the data.name every other test returns.

Over 156 seeded results of every exported test on four designs, one of them a weighted fit, 132 are identical() to those of 0.14.0 and the other 24, of these six tests, are identical() once the added name or data.name is removed.

Input validation

tests/testthat/test-input-validation.R gives every exported test the same four bad inputs: a glm, a model with an NA coefficient, data that lack rows the model was fitted to, and data that lack a variable of the model. Where a test already stopped, it now stops with the message the others give:

• the eight tests that read everything from the model check the data they are given since 0.15.0, with the messages of the tests that use data; the model frame is accepted as it is;
• performBPRandomEffectsTest() names the rows data lacks, where it stopped with "arguments must have same length";
• performWhiteTestStreaming() refuses a glm, where it stopped with "missing value where TRUE/FALSE needed", and it, performHarrisonMcCabeTest() and performVarianceFormTest() refuse it with one message.

No value changes: over the 217 calls of the (model, data, ...) comparison and the 156 results of the return-contract comparison, every result is identical() to 0.14.0, apart from the df names and the data.name above. The tests that returned a value for one of the four inputs refuse it too, which changes what those calls return (Inputs every test refuses, above).

The advice for a small sample

warnDiagnosticAssumptions(), and so generateHeteroRecommendations(), told the user of a sample of fewer than 100 observations to "rely on wild bootstrap or permutation-based tests for reliable inference". In the package's validation study of the multiplier, inst/validation/wild-multiplier-size.R, with B = 999, the wild bootstrap holds its level with Gaussian and t5 errors at 30, 50 and 100 observations, but not with skewed ones: with centred exponential errors it rejects a true null hypothesis 8.2%, 7.3% and 6.1% of the time at the 5% level with the Rademacher multiplier, and 7.8%, 7.0% and 5.9% with the Mammen multiplier. The warning now says so, as the help page of performWildBootstrapTest() does since this release, and asks the user to check the residuals for skewness and to read p-values near the level with caution. Only the text changes: which warnings are given, and when, does not.

0.14.0

0.14.0 reviews the exported names before 1.0.0: it deprecates 22 exports and the registry names "student_bp" and "ncv", and removes nothing. performBreuschPaganTest(), performStudentizedBPTest() and performNCVTest() give way to performBPTest(), performKoenkerTest() and performCookWeisbergTest(), which gains var_formula; the three *Enhanced functions to enhanced = TRUE on their base names; rbootstrap_test_statistic(), test(), simulate_type_I_errors() and simulate_power_analysis() to bootstrapTestStatistic(), runDiagnostics(), simulate_size() and simulate_power(); and twelve internal helpers are to leave the exports (the table under Deprecated names gives each old name and what to use instead). No exported function returns a different value: each deprecated function returns what it returned in 0.13.0 and warns, and performStudentizedBPTest() keeps its own implementation, which disagrees with performKoenkerTest() on a fitWLS() fit until the model frame of such a fit is corrected in 0.16.0. The dashboard and benchmark functions, which move to a companion package, are not deprecated yet because that package has no name, and the deprecated names are removed in 1.0.0.

Deprecated names

0.14.0 reviews the exported names before 1.0.0. Nothing is removed: each name below still works, returns what it returned before and warns with .Deprecated(), naming its replacement where it has one. The names go in 1.0.0.

Deprecated	Use instead	Note
performBreuschPaganTest()	performBPTest()	The same function under a second name
performStudentizedBPTest()	performKoenkerTest()	Same statistic, except on a fitWLS() fit (below)
runHeteroTests(tests = "student_bp")	tests = "koenker"	Warns once per call of runHeteroTests() or runDiagnostics()
performNCVTest(model, var_formula)	performCookWeisbergTest(model, var_formula)	performCookWeisbergTest() gains var_formula; with one, it returns what performNCVTest() returned
performNCVTest(model)	performCookWeisbergTest(model)	Same statistic, degrees of freedom and p-value; labelled as the Cook-Weisberg test
runHeteroTests(tests = "ncv")	tests = "cook_weisberg"	Warns once per call of runHeteroTests() or runDiagnostics()
plotResidualsFittedEnhanced(model)	plotResidualsFitted(model, enhanced = TRUE)	enhanced = FALSE, the default, keeps the plot the base name drew
plotDiagnosticSuiteEnhanced(model)	plotDiagnosticSuite(model, enhanced = TRUE)	As above
checkModelEnhanced(model, data)	checkModel(model, data, enhanced = TRUE)	checkModel(model) still checks the class only
std_error(), std_warning()	None	Internal helpers
checkData(), validateTestInputs()	None	Internal helpers
ht_enable_log_capture(), ht_log_history(), ht_clear_log_history()	None	ht_set_log_level() stays
TestFactory, test_factory	registerDiagnostic(), runHeteroTests()	Objects cannot warn: deprecated on their help page and here only
rbootstrap_test_statistic()	bootstrapTestStatistic()	The same function; each name records its own call in the result
rcalculateEffectSize(), restimate_test_power(), rvalidate_against_reference()	None	Internal helpers
test(x, tests)	runDiagnostics(x$model, x$data, tests)	Its name masks devtools::test(); summary() of a HeteroDiagnostic object is unchanged
simulate_type_I_errors()	simulate_size()	Joins the snake_case simulation family
simulate_power_analysis()	simulate_power()	As above

performStudentizedBPTest() and performKoenkerTest() agree to 1e-8 in statistic, degrees of freedom and p-value in the designs of tests/testthat/test-koenker-studentized.R: a Gaussian null, sd = x^2 and a t5 null at 150 observations, a factor regressor and a weighted fit. Both refuse a rank-deficient fit. On a fit returned by fitWLS() they do not agree: the model frame of that fit still holds the weights of the model fitWLS() refitted, not the estimated ones, and performStudentizedBPTest() reads them from there. On the fitWLS() design of the same file it returns 42.60, where performKoenkerTest() and an lm() fitted with the estimated weights return 1.75. performKoenkerTest() is the correct one. The deprecated name and "student_bp" keep their old implementation in 0.14.0, so no value changes; the model frame is to be corrected in 0.16.0.

rrunAdvancedDiagnostics() stays exported; its name is still under review.

The simulation functions, simulate_*() and sigma_*(), keep their snake_case names on purpose, as a documented family; the rest of the package uses camelCase.

Reference index

The reference index of the pkgdown site groups the functions: the tests of the variance by family, and the model checks (performRESETTest(), performVIFDiagnostic(), performInfluenceDiagnostics(), performScatterDiagnostic() and checkModel()) apart from them, so that they are not mistaken for heteroscedasticity tests. No function changes.

The planned export list

tests/testthat/test-public-api.R lists the exports planned for 1.0.0, the names deprecated in 0.14.0 and the three that move to a companion package, from the export inventory in ROADMAP.md. An export in none of these lists fails the test, so a new export has to be added to the plan on purpose.

0.13.0

The bootstrap tests change value in three cases where their replicates were wrong: on an unweighted fit with an offset given as an argument, which all three dropped; in rbootstrap_test_statistic(), on data whose rows are in another order than the fit's or include rows the fit dropped; and in performWhiteTestBootstrap(), on a transformed response such as log(y), where the p-value was 1 whatever the data. Four more results change value: fitWLS() keeps the offset of the model it refits, and performVarianceFormTest() changes with it; the Szroeter test registered for runHeteroTests(tests = "szroeter"), and so for the functions built on it, orders the observations by the fitted values; runHeteroTests() and runDiagnostics() given a survey design no longer test a fit weighted by its sampling weights; and tidy() returns one row per test. Apart from these, no exported test returns a different value for a call that ran in 0.12.0, though rbootstrap_test_statistic() now refuses data it cannot match to the fit. The bootstrap tests also accept a weighted fit, holding the weights of a fitWLS() fit fixed, and the package gains the Harrison-McCabe test and, as performGlejserTest(robust = TRUE), a Glejser test that keeps its level under skewed errors.

Bootstrap tests on a fit with an offset

rbootstrap_test_statistic(), performWhiteTestBootstrap() and performWildBootstrapTest() keep an offset given as an argument, and change value on an unweighted fit with one. lm(..., offset = o) keeps o in the model frame and not in the formula, and each replicate of an unweighted fit was refitted from the formula alone, so the replicates were those of a model without the offset. A weighted fit already kept it, and an offset written in the formula, offset(o), was always evaluated again from the data. For lm(dist ~ speed, data = cars, offset = 0.06 * speed^2), with set.seed(1) and B = 999, rbootstrap_test_statistic() with performKoenkerTest() gave p = 0.095 and gives 0.106, and performWhiteTestBootstrap() gave p = 0.222 and gives 0.201, the values the same model gives with the offset in the formula. performWildBootstrapTest() dropped such an offset too, where it ran at all (below): for this model its replicates change and its p-value, 0.079, does not. In inst/validation/offset-refits-size.R, where the offset is 0.5 x^2, the White bootstrap rejected 0.012 of the time at the 5% level and now rejects 0.044. rbootstrap_test_statistic() with performKoenkerTest() rejected 0.047, because Koenker's statistic happens to miss what the dropped offset leaves in the residuals, and now rejects 0.040. That is inside the 99% interval around 0.045, [0.037, 0.053], the rate at which a test whose level is exact gives p < 0.05 with B = 199, which is what the cells with B = 199 are judged against; read against 0.05 it would be below the release gate, [0.042, 0.058]. With B = 999, its default, whose exact rate is 0.049, it rejects 0.046, inside the release gate and inside the interval around 0.049, [0.041, 0.057].

performWildBootstrapTest() runs on unweighted fits it stopped on. It evaluated the call of the fit again in the model frame, which holds log(wt) and not wt, and (offset) in place of the name of an offset given as an argument. So lm(mpg ~ log(wt), data = mtcars) stopped with object 'wt' not found, and so did a transformed response, an offset written in the formula and an offset given as an argument, unless the model frame happened to hold the variables they are computed from. Each replicate is now refitted on the design stored with the fit, with its offset. In inst/validation/offset-refits-size.R it rejects 0.044 of the time at the 5% level with an offset given as an argument, and 0.048 with log(x) as well. Calls that ran before give identical values, except with an offset given as an argument (above).

rbootstrap_test_statistic() on reordered data and dropped rows

rbootstrap_test_statistic() matches the rows of data to the fit by name, and changes value when they are in another order or include rows the fit does not have. It took the fitted values of an unweighted fit by position, so data holding the fitted rows in another order regenerated each row from another row's fitted value, without an error. For lm(mpg ~ wt + qsec, data = mtcars) with mtcars sorted by hp, set.seed(1) and B = 999, the p-value with performKoenkerTest() was 0.004 and is 0.230; with mtcars as given it is 0.219. In inst/validation/offset-refits-size.R, with the rows in a random order, it rejected 0.257 of the time at the 5% level and now rejects 0.049.

Rows of data the model was not fitted to, such as rows it dropped for missing values, are now left out, weighted or not, for the observed statistic and the replicates alike; the rows that stay keep the order of data, so data holding exactly the fitted rows gives the same values as before. Pairs resampling drew those rows too, and each replicate was refitted to a different number of rows. For lm(Ozone ~ Solar.R + Wind, data = airquality), whose fit drops 42 rows, pairs resampling with performKoenkerTest(), set.seed(1) and B = 999 refitted replicates of 92 to 126 rows and gave the interval 2.884 to 18.904; every replicate now has the 111 rows of the fit, and the interval is 3.134 to 19.599. Null resampling refused such data and now runs. Data that lacks a row the model was fitted to is refused.

Matching by name cannot see two things. A tibble or a data.table has no row names of its own; when it has as many rows as the fit, as the complete cases do, they are taken in the order of the fit. And a variable the formula finds outside data, such as offset(o) with o in the environment of the formula, keeps the order of the fit, so it cannot follow data in another order nor the rows that pairs resampling draws; given as an argument, the same offset does. The model is therefore refitted to the rows of data, and data whose refit does not give the residuals of the fit, to a relative 1e-8, is refused, weighted or not: rows in another order taken by position, values that differ from those of the fit, or such a variable. These gave wrong replicates without an error. A term computed from all the data, such as I(x - mean(x)), does not count as a difference. The package's own tests refuse a variable that is not a column of the data, so only a test function of the caller's own could reach that case.

performWhiteTestRobust(bootstrap = TRUE) and performBPTestRobust(bootstrap = TRUE) resample through rbootstrap_test_statistic(), so their bootstrap p-values change with it, here and with an offset (above).

The White bootstrap on a transformed response

performWhiteTestBootstrap() regenerates a transformed response, and changes value on a model with one. For a response such as log(y) in log(y) ~ x, it wrote each regenerated response into a column named log(y), which the refit did not read: the refit computed log(y) from y again, so every replicate was the data of the fit, every replicate statistic was the observed one, and the p-value was 1 whatever the data. For lm(log(Ozone) ~ Solar.R + Wind + Temp, data = airquality), with set.seed(1) and B = 999, all 999 replicate statistics were the observed one, 34.715, and the p-value was 1; it is now 0.004, the value the same model gives with log(Ozone) held in a column of its own, where performWhiteTest() gives 0.00007. The regenerated response now goes into a column of its own, which the refit reads in place of log(y). In inst/validation/white-bootstrap-transformed-size.R, with constant variance on the log scale, the test rejects 0.047 of the time at the 5% level at 5000 replications, inside the interval [0.037, 0.053] around 0.045, the rate at which a test whose level is exact gives p < 0.05 with B = 199, and 0.892 against a standard deviation that grows with a regressor. Models with a plain response give identical values.

fitWLS() on a model with an offset

fitWLS() keeps the offset of the model it refits, and changes value on a model with one. It refitted the model with the estimated weights without the offset, whether that was given as an argument, lm(..., offset = o), or written in the formula, offset(o), so the weighted fit was that of a model without it. For lm(dist ~ speed, data = cars, offset = 0.06 * speed^2) it returned the coefficients -13.996 and 3.685; it now returns -4.356 and 2.038, those lm() gives for the same call with the estimated weights. The weights themselves do not change: they come from the residuals of the model, which carry its offset. performVarianceFormTest() refits through the same code and changes with it: on that model F = 0.290 (p = 0.593) before and F = 0.773 (p = 0.384) now. In the design of the guard added to tests/testthat/test-fitwls-offset.R, an exponential variance function that is correctly specified and the offset 2 x^2, it rejected 70% of the time at the 5% level when the offset was a variable of the calling environment, which the old code could find, and rejects 4.8% now. autoCompareRemediations() and runHeteroTests(tests = "variance_form") change with them. Models without an offset give identical values.

fitWLS() and performVarianceFormTest() also run on models they stopped on. The weighted refit evaluated the call of the fit again in its model frame, which holds log(wt) and not wt, and weights and an offset given as arguments as (weights) and (offset). So fitWLS(lm(mpg ~ log(wt), data = mtcars)) stopped with object 'wt' not found, a model fitted with weights that are a column of the data with object 'w' not found, and a model whose offset is a column of the data with object 'o' not found. The refit now uses the design stored with the fit. On a model fitted with weights, fitWLS() replaces them with the estimated ones, as it did when the old refit could find them; those calls give identical values.

The Szroeter test run by name

runHeteroTests(tests = "szroeter") changes value. The registered Szroeter test ordered the observations by the second column of the data, whatever it was: the response, a regressor or a column the model does not use, so the result depended on the order of the columns. For lm(mpg ~ wt + hp, data = mtcars) with mtcars as given it ordered by cyl, which the model does not use, and gave Q = -0.709 (p = 0.761); depending on which column came second, Q ran from -1.419 to 1.333. It now orders by the fitted values, so the alternative is that the variance rises with the mean, and gives Q = 1.099 (p = 0.136) whatever the order of the columns. 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, and ties keep the order of the rows the model was fitted to.

Everything that runs the registered test changes with it: runDiagnostics(), test() and summary() on a HeteroDiagnostic, compareModelDiagnostics(), compareTestResults(), runSurveyHeteroTests(), runHeteroTestsParallel(), cachedTest() and the dashboard, when given "szroeter". performSzroeterTest() does not change, and no number recorded in inst/validation/ moves: the studies call it with an explicit ordering.

Survey designs and their sampling weights

Survey designs are no longer tested as if their sampling weights were precision weights, and runSurveyHeteroTests() says what it does.

The help page of runSurveyHeteroTests() described a survey-weighted fit. The weights never reached the fit: they were passed to lm() as a local variable, where lm() does not look, and a tryCatch() refitted without them, so the function ran the ordinary tests on an ordinary least-squares fit. It still does, and apart from the Szroeter test (above) its values are unchanged. The help page now says that the design is not used and when the tests hold their level regardless, and the function warns when the design has unequal weights or clusters.

runHeteroTests(formula, design) and runDiagnostics(formula, design) took a different route. They fitted the model with survey::svyglm() and tested that fit, whose residuals are multiplied by the square roots of the sampling weights. Sampling weights say how many population units an observation stands for, not how variable its error is. In inst/validation/survey-designs-size.R, on homoscedastic data sampled with probabilities that depend on the regressor, that statistic rejected 97% to 100% of the time; on heteroscedastic data it rejected 6% to 12% of the time where the unweighted test rejected 98% to 100%. These two calls change value: they now do what runSurveyHeteroTests() does, with the same warning. A svyglm fit passed to a test directly is refused.

The package has no design-based test of constant variance. The same study tried two, a Wald and a score test on the regression of the squared residuals on the regressors, estimated with the design. They are the only procedures near their level when selection depends on the response or the sample is clustered, where the function itself rejects 37% to 100% of the time. In the designs the ordinary tests handle correctly they rejected 6% to 25% of the time, and 3 of their 40 size cells fell inside the Monte Carlo band. They were not adopted.

tidy() of a test with two parameters

tidy() returns one row per test. A result with two parameters, such as the F statistic of performVarianceFormTest() or the wild bootstrap with its degrees of freedom and number of replications, came back as two rows, and every table and plot built from the tidy output counted the test twice. The second parameter is now in a new column, parameter2, and estimate is given when a test reports exactly one. glance() follows.

Bootstrap tests on a weighted fit

performWildBootstrapTest(), performWhiteTestBootstrap() and rbootstrap_test_statistic() accept a weighted fit. In 0.12.0 they refused one, because each replicate was refitted without the weights. A weighted fit is ordinary least squares on sqrt(w) y and sqrt(w) X, where the null hypothesis is constant variance, and the three now resample there: they draw from the Pearson residuals, divide each draw by sqrt(w) to return it to the scale of its observation, and refit with the weights. performWhiteTestRobust(bootstrap = TRUE) and performBPTestRobust(bootstrap = TRUE), which resample through rbootstrap_test_statistic(), accept a weighted fit too. Accepting a weighted fit changes nothing for an unweighted one; there the same seed gives the same values as before, except in the three cases corrected above: an offset given as an argument, data in another order or with rows the fit dropped, and a transformed response in the White bootstrap.

In inst/validation/weighted-fits-size.R, at 2000 replications with B = 499, the three reject a correctly weighted fit 0.045 to 0.054 of the time at the 5% level in two designs, and reject weights that are too weak in every replication.

On a fitWLS() fit they over-reject. They hold the estimated weights fixed in each replicate, as the asymptotic tests take them as known, and share their error. In a design where the variance function is correctly specified and the true weights give 0.054 with 150 observations and 0.048 with 600, the three reject 0.064 to 0.086 of the time at the 5% level with 150 observations, and the wild bootstrap 0.105 with 600; performKoenkerTest() rejects 0.077 and 0.110. Estimating the weights again in each replicate, which the package does not do, brings the three to 0.042 to 0.060 with 150 observations and the wild bootstrap to 0.051 with 600. That is why the question is open: ROADMAP.md proposes an argument that does it, and leaves its default to 0.14.0. A default that estimates the weights again will change the values these tests return on a fitWLS() fit.

performQuantileRegressionTest() still refuses a weighted fit, and its help page says why.

The Harrison-McCabe test

performHarrisonMcCabeTest() adds the test of Harrison and McCabe (1979) beside performGQTest() and performSzroeterTest(), the other two tests along an ordering of the observations. The statistic is the share of the residual sum of squares that falls in the first part of the sample ordered by a variable, and it agrees with lmtest::hmctest() to 1e-8. The p-value is exact under normal errors given the regressors, by Imhof's (1961) formula, and reproduces the exact critical values published by Harrison (1981). The 1979 article was not read: the definition rests on Harrison (1981), by its first author, and on lmtest. lmtest::hmctest() simulates its p-value from a regression on a constant alone, whatever the model, and skedastic::harrison_mccabe() drops repeated eigenvalues before applying Imhof's formula, which changes its p-value. Weighted fits are supported.

In inst/validation/harrison-mccabe-size-power.R, at the 5% level, the test rejects a true null hypothesis 4.3% to 5.3% of the time under Gaussian errors at 40 to 400 observations, with and without weights. Under t5 errors it rejects 11.7% to 16.3%, as the Goldfeld-Quandt and Szroeter tests do, and the help page points to performGlejserTest(robust = TRUE) there. Against the two alternatives of Pass A it has 2.3 to 4.8 percentage points less power than the Goldfeld-Quandt test and 5.4 to 7.9 less than Szroeter's. No value any other function returns changes.

A Glejser test for skewed errors

performGlejserTest() gains robust = TRUE, the Glejser test of Im (2000) and of Machado and Santos Silva (2000), which keeps its level when the errors are skewed. No other R implementation of it was found.

Glejser's statistic has its nominal level when positive and negative errors are equally likely. When they are not, and the variable tested is correlated with the regressors, it does not in general, and a larger sample does not help. In inst/validation/glejser-skewed-errors.R, at the 5% level and with 50 to 1000 observations, it rejected a true null hypothesis 12% to 14% of the time under exponential errors and 15% to 18% under lognormal errors, which is what the asymptotic theory gives for those designs. The help page said the test was not valid there and pointed to Koenker's test.

With robust = TRUE the auxiliary regression is run on |e| - m e, where m is the proportion of positive residuals less the proportion of negative ones. In the same study it rejected 4.4% to 5.7% of the time from 150 observations on, and 5.7% to 6.6% at 50 under skewed errors. Under Gaussian errors it gives up one or two percentage points of power at 50 observations. Against Koenker's statistic on the same regressor it has 1 to 4 points less power under Gaussian errors and 1 to 10 points more under exponential errors.

The default is robust = FALSE, and every value the function returned before is unchanged.

Measured and not adopted

Two of the candidate tests in ROADMAP.md were measured and are not in the package, because they do not hold their level. The studies ship in inst/validation/, and the roadmap says what either would need.

• The modified Wald test for group-wise heteroscedasticity in a fixed-effects panel, as Stata's xttest3 computes it. In inst/validation/modified-wald-size-power.R, at 5000 replications, it rejects a true null hypothesis at the 5% level 10% to 100% of the time in balanced panels of 10 to 100 units and 5 to 100 periods, under Gaussian and t5 errors.
• The consistent test of a variance function of Dette, Neumeyer and Van Keilegom (2007), as a prototype that follows the authors' manuscript. In inst/validation/variance-form-dnvk-size-power.R, at 5000 replications, its Cramér-von Mises statistic rejects a true exponential or power variance function that is not constant 1.9% to 3.6% of the time at the 5% level, at 50 to 400 observations under Gaussian and t5 errors.

Also

• listDiagnostics() returns the names runHeteroTests(tests = ) accepts, sorted: the diagnostics the package registers and any added with registerDiagnostic(). No exported function listed them before. Given a name that is not registered, runHeteroTests() now stops with Unknown tests: <name>. listDiagnostics() returns the registered names. It used to say only Unknown tests: <name>, in an error raised from the internal .run_hetero_core().

• runHeteroTests() with a fitted model and grouped data now runs. The per-group refit passed the weights to lm() as a local variable, where lm() does not look, and stopped with invalid type (closure) for variable '(weights)' for every lm, glm and parsnip fit, weighted or not. Only a formula with grouped data worked. Each group is now fitted by evaluating the model's own call on it, so a weighted model uses that group's rows of the weights, a glm keeps its family and a subset is applied within the group. A model whose weights, offset or subset are vectors in the workspace cannot be split by group and is refused with a message that says so.

• autoplot() of a test suite showed nothing when the tests rejected. The bars were the p-values on a linear axis from 0 to 1, so a test with p = 0.002 had no visible bar, and the colouring of significant tests was decided by isTRUE(p < 0.05), which is FALSE for a suite of more than one test. The plot now draws one horizontal bar per test whose length is -log10(p), on an axis labelled with p-values, with a dashed line at a new alpha argument (default 0.05) and the p-value printed beside each bar. The grouped method draws one panel per group; it used to add a panel variable for the suggestions column as well.

• plotBeforeAfter() draws the two models in two panels, each with its own axes, where it used to overlay them in one panel in two similar blues. Weighted residuals and the residuals of a transformed response are not in the units of the original ones, so a common axis could flatten one of them.

• print() on the result of suggestRemediation() now uses its print method. The method existed but was not registered in NAMESPACE, so results printed as a bare list with an attr(,"class") line. This also affects the dashboard's suggestions panel and the remediation section of generated reports, which print the same object.

• suggestDiagnosticsForProfile() had two definitions, and only the one loaded last ever ran. The other still recommended "hc_covariance", a diagnostic removed in 0.8.0, for data of 10,000 rows or more; had it been the one loaded, passing its recommendations to runHeteroTests(tests = ) would have stopped with Unknown tests: hc_covariance. It is gone, and the function returns exactly what it returned before. New tests check that every name it can return is a registered diagnostic, and that no function in the package has two different definitions.

• The cache that keeps the result of each input validation, so that the same model and data are not checked twice, now keeps the 1000 most recent results. It kept every one, and a test on a new model, as in a simulation, never reuses one, so it grew by three entries with each test: after 10^5 calls of performWhiteTest() on new models with 150 observations it held 300,001 entries and R's memory had grown by 655 MB, where it now holds 1000 and has grown by 3 MB. Simulation workers had reached 7.6 to 8.9 GB each. A test repeated on the same model still reuses the result, and no value changes.

• The help page of performStudentizedBPTest() called performKoenkerTest() "the absolute-residual variant", which it has not been since 0.6.5. Both compute Koenker's studentized n R^2 statistic from the regression of the squared residuals on the regressors, and the page now says so. It also cited a section of Davidson and MacKinnon (2004) that does not exist; the test is in Section 7.5. The vignette on statistical theory gave the classical Breusch-Pagan statistic as n R^2, which is Koenker's, and said that performKoenkerTest() uses the absolute residuals; it now gives each statistic as the package computes it. No value changes.

• The help page of performWhiteTestBootstrap() gave its p-value as the share of replicates whose statistic is at least the observed one. Since 0.7.0 the function has returned (1 + #{T_b >= T_obs}) / (B_eff + 1), where B_eff counts the replicates with a finite statistic, and the page now says so. Its example of parallel = TRUE now runs only where that starts at most two processes, the limit of R CMD check --as-cran. No value changes.

• The release workflow runs R CMD check --as-cran --run-donttest on the source tarball it is about to tag, and creates no tag or GitHub release if that check warns or errors. The check used to run only weekly and on demand, and it failed on the example above on 2026-09-28 and 2026-10-05 without holding up 0.12.0.

0.12.0

Tests on weighted fits change value, and the package gains a test of the shape of the variance. No test changes value on an unweighted fit: re-running the full sweep reproduces its 32 existing rows digit for digit. fitWLS() with its default arguments returns the weights it returned in 0.11.2 for any model with an intercept.

Tests on a weighted fit read the wrong residuals

lm(..., weights = w) states that the error variance is sigma^2 / w_i. The raw residuals of such a fit are heteroscedastic by assumption, and every test read them. A model fitted with exactly the right weights was therefore rejected as often as the same model fitted with none, and the workflow the README and the vignettes describe -- diagnose, refit with fitWLS(), test again -- could not succeed.

inst/validation/weighted-fits-size.R measures it: data with sd = x^2, 150 observations, a model fitted with weights 1 / x^4, 400 replications. Run against 0.11.2, the 27 heteroscedasticity diagnostics that ran all rejected in 400 samples out of 400, and three could not run at all. The 27 that accept a weighted fit now reject between 2.2% and 7.7% of the time, against a Monte Carlo standard error of 1.1%. A fit from fitWLS() fared no better: with the variance function correctly specified, performKoenkerTest() on 0.11.2 rejected the weighted fit 99.9% of the time.

The tutorial vignette showed the fault on real data. For quakes, compareModelDiagnostics(list(base_model, wls_model)) reported a White statistic of 125.6 for the least-squares fit and 249.2 for the weighted one, and Breusch-Pagan statistics of 191.9 and 463.7, directly above the sentence "The weighted least squares fit typically reduces the heteroscedasticity metrics." It now reports 42.2 and 40.8 for the weighted fit.

On a weighted fit the tests are now computed from the Pearson residuals sqrt(w_i) * e_i. Those are the residuals of the equivalent unweighted regression of sqrt(w) y on sqrt(w) X, and they have constant variance when the weights are right. car::ncvTest() and plot.lm() make the same choice. The null hypothesis becomes that the weights are adequate.

• Every residual-based test reads the Pearson residuals. Variance regressors, grouping variables and ordering variables stay on their original scale.
• performNCVTest() and performCookWeisbergTest() now reproduce car::ncvTest() on a weighted fit to 1e-8. They did not before.
• performGQTest() refits each segment with its weights and compares weighted residual sums of squares. It reproduces lmtest::gqtest() on the transformed regression to 1e-8 at fraction 0.1, 0.2 and 0.3.
• performRESETTest() fits the augmented model by weighted least squares. Before, it compared a weighted restricted fit with an unweighted augmented one, which are not nested, and rejected 23.0% of the time under the null in the design above. It now rejects 5.5%.
• performArchLMTest() and performMcLeodLiTest() read a variance that rises along the sample as volatility clustering, whatever the weights: 72.3% and 79.7% on a correctly weighted fit. Now 4.0% and 4.8%.
• performStudentizedBPTest() used the raw residuals and weighted the auxiliary regression instead, the treatment in lmtest::bptest(), which tests the raw residuals and is therefore not a reference for weighted fits. It now agrees with performKoenkerTest() on weighted fits as it does on unweighted ones.
• performStudentizedBPTest(), performBPTestRobust() and performWildBootstrapTest() failed with object 'w' not found whenever the weights were a column of the data. They rebuilt the design matrix from the model frame, where that column is stored as (weights).
• The residual plots (plotResidualsFitted(), plotSpreadLevel(), plotBeforeAfter() and the others) and the variance-pattern summary of the recommendation engine use the Pearson residuals of a weighted fit. plotBeforeAfter() drew a weighted remedy as heteroscedastic as the original however good the weights were.
• performWildBootstrapTest(), performWhiteTestBootstrap(), rbootstrap_test_statistic() (and so bootstrap = TRUE in the two *Robust tests) and performQuantileRegressionTest() refit the mean model without the weights. They now refuse a weighted fit, as rbootstrap_test_statistic() already refused a glm, instead of returning a statistic for a different model.
• A fit with a zero weight is refused: the observation is not part of the fit and its Pearson residual is an artificial zero.

Weights that were estimated are not weights that are known

The reference distributions take the weights as given. That is exact for weights known in advance and approximate for a fitWLS() fit, whose weights come from the same residuals, and the approximation does not improve with the sample size. With the variance function correctly specified and Gaussian errors, at the 5% level (2000 replications, standard error 0.5%):

Test on a fitWLS() fit	n = 150	n = 600
performKoenkerTest()	7.7%	11.0%
performWhiteTest()	6.3%	8.3%
performBPTest()	10.6%	11.8%
performHarveyTest()	0.0%	0.0%
performVarianceFormTest()	4.8%	5.2%

Harvey's test never rejects because it repeats the regression the weights were estimated from. The statistics remain useful for comparing a weighted fit with the original, which is what compareModelDiagnostics() does, but their p-values on a fitWLS() fit are not calibrated. The help page of performKoenkerTest() says so, and the last row is the test built for that case.

performVarianceFormTest(): is the variance function the right shape?

Every other test in the package has constant variance as its null hypothesis. Rejecting it says that the variance moves with some variables, not whether an exponential or a power function of them describes it. The new test takes a variance function as its null hypothesis.

It fits the variance function, refits the model by weighted least squares, and regresses the squared standardized residuals on the regressors of the variance function and on terms added to them: by default their squares and pairwise products, or the terms of a one-sided formula given as against. The statistic is the F test of the added terms, the regression-based construction of Wooldridge (1990, 1991). Keeping the variance regressors in that regression is what makes the test valid with estimated weights, and it is also where the information is. A Koenker test on the weighted fit asks about directions the estimation has already fitted: with variance x^3 fitted as exponential, at 600 observations, it rejected 6.9% of the time, less often than the 11.0% it rejects when the fitted form is right. The new test rejected 92.9%. The White test, which also contains squares, rejected 87.9%, but its level on a fitWLS() fit is not 5%.

inst/validation/variance-form-size-power.R runs both forms against five true variance functions at 50, 150 and 400 observations with Gaussian and t_5 errors, 5000 replications a cell:

• When the form under test is the true one, or the variance is constant, the rejection rate runs from 4.2% to 5.7% over the 24 cells, and all of them are inside the release-gate interval [0.042, 0.058].
• Against a log-variance that is quadratic in x, power is 51.7% at 50 observations and 98.0% at 150 for the exponential form.
• Against the other log-linear form, power grows more slowly: 22.4% at 150 observations and 74.1% at 400 for an exponential null and variance x^3.
• Against sigma^2 = 0.5 + x, which is close to both forms over the range of x, it rejects between 4.3% and 6.6% of the time. The test cannot separate shapes the data cannot separate, and the help page says so.

The test accepts either the unweighted model with var_formula and form, or a fitWLS() result, whose variance function it then tests. It is registered as "variance_form" for runHeteroTests().

fitWLS() fits the variance function it is asked for

fitWLS(model, data = NULL, var_formula = NULL, form = c("exponential", "power")). var_formula names the variance regressors and form the shape: exponential, Harvey's (1976) log sigma^2 = g0 + z'g, or power, Park's (1966) sigma^2 proportional to z^d. The defaults reproduce 0.11.2. The fit carries the variance function as its "variance_function" attribute, next to the "variance_model" attribute it already had.

One value changes. The variance function now always has its own constant. fitWLS() regressed log(e^2) on the model's design matrix, so for a mean model without an intercept the log-variance had no constant either, which pins the variance at 1 where the regressors are zero. Models with an intercept are unaffected.

Withdrawn before release

An additive variance function, sigma^2 = a + z'b, was implemented as a third form and removed. Nothing keeps its fitted variances positive. Under its own null, at 5000 replications, the test rejected 6.1% of the time with 150 and with 400 observations, outside the release gate, and the fitted variances were not all positive in 17% of samples at 50 observations and 6% at 150. It needs a constrained estimator before it can be offered.

Also

• The Park (1966) reference pointed at another paper. Its DOI is 10.2307/1910108 and it is a one-page note, Econometrica 34(4), 888. The help page cited 10.2307/1909774, which is Henshaw (1966), with pages 888-898, and the theory vignette gave 888-908.
• tests/testthat/test-weighted-fits.R and test-variance-form.R pin the corrections by simulated size and power as well as by reconstruction, since no stored value and no comparison with lmtest::bptest() would have caught the fault.
• The full sweep covers performVarianceFormTest(): size 5.8%, power 97.8% against a log-variance quadratic in x.
• The tutorial vignette gains a step that tests the variance function of its weighted fit.

0.11.2

Changes requested by the CRAN review of the 0.11.1 submission. No statistic, default or exported function changes.

• The Description field no longer says "in R", and now cites the methods the package implements, with DOIs: White (1980), Breusch and Pagan (1979), Koenker (1981), Goldfeld and Quandt (1965), Harvey (1976), Cook and Weisberg (1983), Engle (1982), Bartlett (1937), Brown and Forsythe (1974) and Hartley (1950).
• No example reaches an unexported object through :::. The six help pages that did -- for rhandleMissingValues(), rTEST_REQUIREMENTS, rvalidateDataInputs(), rvalidateDistributionalAssumptions(), rvalidateGroupingVariable() and rvalidateModelInputs() -- document internal objects, where :: would fail. Their examples are removed rather than the objects exported.

• No example uses \dontrun{}. The recommendations example runs in well under a second and is unwrapped. The report example takes about ten seconds and needs pandoc, so it runs under \donttest{} behind rmarkdown::pandoc_available() and writes to tempdir(). The dashboard example is guarded by interactive(), since printing the Shiny application it returns starts a blocking server.

• The \donttest{} part of the performWhiteTestBootstrap() example uses 499 bootstrap replications rather than 999. With that part included, as R CMD check --as-cran runs it, the example took 5.0 to 5.7 seconds, over the five-second limit for examples; it now takes about 1.7 seconds.

0.11.1

Two panel p-values change value, though not their conclusions.

The panel p-values no longer underflow to zero

performBPRandomEffectsTest() and performPesaranTest() computed their upper tails as 1 - pchisq() and 1 - pnorm(). Both collapse to exactly zero once the statistic is large: at LM = 182.9 the subtraction returns 0 where the upper tail returns 1.1e-41, which is the value plm::plmtest() reports. A p-value of exactly zero is wrong, and it discards the magnitude that anything working on a log scale or applying a multiplicity correction needs. Both now pass lower.tail = FALSE.

Both are now checked against plm

Neither statistic had a reference implementation among the packages the accuracy table compares against, which is why both were wrong until 0.11.0 and why simulated size was the only thing that caught them. plm has had both all along.

performBPRandomEffectsTest() reproduces plm::plmtest(type = "bp") exactly, and performPesaranTest() reproduces plm::pcdtest(model = "pooling"), both to 1e-8 across N = 20/30/40 and T = 5/6/8.

The model argument is the part worth recording. pcdtest() defaults to model = NULL, which tests the series rather than any model's residuals and gives 0.817 where this package gives 1.580. performPesaranTest() takes a fitted lm, whose residuals are identical to the pooling residuals, so "pooling" is the comparison that puts both on the same quantity; "within" and "random" give 1.220 and 1.258, which are different estimands rather than disagreements.

plm and withr join Suggests. The tests attach plm with withr::local_package() rather than calling through ::, because pcdtest() calls plm() internally without qualifying the namespace and fails when the package is only installed.

Continuous integration

• R CMD check --as-cran no longer runs on every merge. CRAN submission is not the near-term target, and that job costs about an hour: it installs a LaTeX toolchain and every check dependency, then rebuilds the vignettes. It had also stopped completing -- three consecutive runs on main were cancelled or left queued inside it -- so it was a gate producing a signal nobody received. It now runs weekly and on workflow_dispatch.
• A concurrency group is added, keyed by pull request number for pull requests and by commit for everything else. A ref-keyed group would have put pushes, the weekly schedule and manual dispatches on main into one queue, and GitHub keeps only one pending run per group, so a third arrival would cancel the pending one.

0.11.0

Two panel statistics were wrong against their published definitions. Both change, hence the minor version.

Where these came from

A size and power sweep over all 32 exported perform*Test() functions, each driven by a null and an alternative appropriate to what it actually tests. Twenty-six of the twenty-eight heteroscedasticity tests hold their nominal level: at 400 replications and n = 150, size ranges from 3.0% to 6.2% against a Monte Carlo standard error of 1.1%, with power between 95% and 100% against sd = x^2. The two below did not, and neither has a reference implementation among the packages the accuracy table compares against, so only simulation would have found them.

performBPRandomEffectsTest() rejected a third of the time under the null

The statistic was

T^2 / (2(T - 1)) * sum_i (sum_t e_it)^2 / sum_it e_it^2

against Breusch and Pagan's (1980) equation 5,

nT / (2(T - 1)) * [ sum_i (sum_t e_it)^2 / sum_it e_it^2 - 1 ]^2

The bracketed ratio is close to one under the null, and the published statistic measures its squared departure from one. Dropping the - 1 and the square leaves a quantity that sits at T^2 / (2(T - 1)) -- 3.6 for T = 6, against a chi-square(1) critical value of 3.841. Measured median statistic 3.548, and size 32.5% against a nominal 5%. It is now 5.8%, with power 100% against random intercepts.

The documentation was wrong in the same place. The help page was titled a test for random effects but described a test for "heteroscedasticity in random effects models", and the roxygen claimed it tested whether the effect's variance varies across units. It tests whether an individual effect is present at all, and does not respond to heteroscedasticity: against errors with sd = x^2 the corrected statistic rejects 3.2% of the time, which is the nominal level rather than power.

performPesaranTest() never rejected

The statistic was sqrt(N(N-1)/(2T)) times the mean pairwise residual correlation, against Pesaran's (2004, 2015)

CD = sqrt( 2T / (N(N-1)) ) * sum_{i<j} rho_ij

which is the same quantity multiplied by T. The implemented version was therefore exactly 1/T too small: at T = 8 its standard deviation was 0.133 where the published form gives 1.067, and since the reference distribution is standard normal, size was 0.0%. It is now 6.5%, with power 100% against a common time factor.

The pairwise correlations are also computed with a single cor() call rather than a double loop that grew a vector by c().

Also

• performRESETTest was checked in the same sweep and is correct: size 5.2%, power 100% against an omitted quadratic term.
• performGQTest reproduces lmtest::gqtest() exactly at matched fraction -- the difference is zero at 0.1, 0.2 and 0.3. It cannot reproduce that function's default, because lmtest omits no central observations by default and performGQTest() requires fraction strictly between 0 and 1.
• tests/testthat/test-panel-statistics.R pins both corrections by simulated size rather than by a stored value, since a stored value would not have caught either fault.

0.10.0

rbootstrap_test_statistic() now resamples under the null. Its p-value changes, and so does the p.value_bootstrap field of the two robust tests, hence the minor version.

The bootstrap p-value could not detect anything

The helper resampled rows of the data with replacement, refitted, and recomputed the statistic. Those replicates follow the statistic's distribution under whatever variance structure the data actually have, not under the null: they are centred on the observed statistic, so comparing the two returns roughly one half however strong the heteroscedasticity is.

Measured over 400 samples at n = 120, B = 199, testing Koenker's statistic:

procedure	size, homoscedastic	power, sd = x^2
pairs resampling (before)	0.0%	0.0%
null-imposed (now)	5.0%	100%
performWildBootstrapTest()	2.5%	100%
Koenker, asymptotic	5.0%	100%

Zero rejections under a strong alternative. The p-value was exported directly and also reached users as p.value_bootstrap from performWhiteTestRobust() and performBPTestRobust() when called with bootstrap = TRUE.

resample now takes "null", which is the default. The response is regenerated as fitted + e*, with e* drawn with replacement from the residuals after dividing by sqrt(1 - h_i) and centring, so the regenerated data satisfy homoscedasticity by construction. Size is 5.0% (Monte Carlo standard error 1.1%) and power is 99.2% against sd = 0.5x and 100% against sd = x^2. Size holds at 5.0% under homoscedastic t_5 errors as well, so the calibration does not depend on Gaussian tails -- which matters here, because Bartlett's and Hartley's tests reject about 24% of the time against the same heavy-tailed null (see inst/validation/README.md).

The leverage correction is not decorative. OLS residuals have variance sigma^2 (1 - h_i), so resampling them unscaled under-disperses the regenerated errors; a version without it rejected about 13% of the time under the null. This is the construction performWildBootstrapTest() already used.

resample = "pairs" is still available and still returns the replicates and the interval, which describe the statistic's variability perfectly well. It now returns NA for p_value rather than a number that cannot be interpreted.

Null-imposed resampling is restricted to a Gaussian lm whose response is a plain variable, and errors otherwise. Two cases previously ran and returned a plausible p-value from replicates that meant nothing: a transformed response such as log(y) ~ x, where fitted() is on the log scale and the data column holds y, so refitting took the logarithm twice and returned p = 1; and a glm, where the refit used safe_lm() and dropped the family and link, returning p = 0.952 for a Poisson model. resample = "pairs" resamples rows and is unaffected by the transformed-response case, so it remains available there.

A glm is now refused by both strategies. Each replicate is refitted with least squares whichever strategy is used, so on a Poisson fit of counts the coefficients move from (0.457, 0.406) on the log link to (-0.931, 2.281) on the identity scale: the replicates describe a different model from the one passed. The helper documented glm support it never had.

Percentile intervals are not confidence intervals

The ci component was documented as a "percentile confidence interval for the statistic". It is a percentile interval of the bootstrap distribution: it summarises where the resampled statistic falls, is not an interval for a parameter, and carries no coverage guarantee. The documentation says so. This closes the corresponding roadmap item.

0.9.0

fitWLS() now estimates its weights from a variance model. The numbers it returns change, hence the minor version.

fitWLS() produced unusable standard errors

The weights were the inverse squared residuals of the initial fit, w_i = 1 / e_i^2. A squared residual is a one-degree-of-freedom estimate of sigma_i^2, far too noisy to invert, and inverting it hands the weight to whichever observations the first fit happened to reproduce most closely. On quakes that meant five of a thousand points carrying 99.5% of the total weight, with a largest-to-median weight ratio of about 1e+07.

Two consequences, measured over 2000 replications of y = 1 + 2x + e, sd(e) = 0.5x, at n = 200:

estimator	SD(beta1)	coverage of nominal 95%
OLS	0.1027	93.7%
fitWLS() before	0.1022	10.4%
fitWLS() now	0.0837	96.1%
oracle weights	0.0827	95.1%

The point estimate was unbiased but no more efficient than OLS, so the function delivered none of the efficiency that is the purpose of weighting; and its intervals covered the truth about a tenth of the time, because the weighted residual sum of squares collapses towards n and the reported sigma is approximately sqrt(n / (n - p)) whatever the data.

fitWLS() now regresses log(e^2) on the model's own design matrix and weights by 1 / exp(fitted), the standard feasible-GLS recipe. That recovers the efficiency (SD 19% below OLS, close to the oracle) and calibrates the intervals. Zero residuals are floored before the logarithm by the same helper the log-variance tests use, and an unusable variance model falls back to equal weights, which degrades to the original OLS fit rather than failing.

The estimated variances are attached to the result as the "variance_model" attribute.

Knock-on corrections

• compareModelDiagnostics() can now run White's test on a WLS fit. It previously reported NA there, correctly, because the degenerate weighting made the weighted fit perfectly explained; that cell has now held three different things across three releases and the test records all three.
• autoCompareRemediations() compares AIC and RMSE across OLS, WLS and robust fits. The WLS row was computed from the degenerate fit and is now meaningful.
• The WLS section of real_world_case_studies.Rmd no longer needs its warning that the weighted variance profile was flat by arithmetic rather than by fit. On quakes the profile now falls from a 6.6-fold spread across thirds of the fitted range to 2.95-fold, which reflects the data.

Metadata

• The maintainer affiliation is now recorded identically in CITATION.cff, inst/CITATION, .zenodo.json, README.md and the paper. Only .zenodo.json carried the current name; the rest still had the former one.

0.8.1

CRAN preparation. No change to any statistic or to the public API.

The bundled Boston dataset is gone

boston_housing was a byte-identical copy of MASS::Boston -- all 506 rows and all 14 columns. Shipping it under this package's Apache License (>= 2.0) put a verbatim copy of GPL-2 | GPL-3 material under Apache terms, which is at best a question a CRAN reviewer would ask, and cran-comments.md volunteered it.

The justification recorded there -- that the copy let examples run "without attaching MASS" -- did not hold either: MASS is a hard Imports dependency, so it is installed regardless.

Examples, vignettes and the tutorial notebooks now use R's built-in quakes data, which lives in datasets and is attached by default. Nothing is redistributed, so the licensing question disappears rather than being argued. The two datasets the package still ships, diagnostic_data and hetero_data, are simulated.

The replacement was chosen on evidence, not convenience. Against stations ~ mag + depth (n = 1000, R-squared 0.74) White gives p = 2e-25, Breusch-Pagan p = 2e-42 and Koenker p = 7e-25, and the HC3 standard error for mag is 1.34 times the OLS one -- enough for the notebooks' robust-versus- classical comparison to keep its point. The response is a count, so its variance rises with its mean: the textbook mechanism, which makes it a better teaching example than the original. Candidates with stronger p-values (ChickWeight, Theoph, DNase, Loblolly, Orange) were rejected because they are repeated-measures data, where residuals are correlated within subject and these tests' independence assumption does not hold.

Dropping the dataset also drops the black column, the transformed racial-composition variable that led scikit-learn to remove this data in 1.2.

Packaging

• Three references used \url{https://doi.org/...} where CRAN asks for \doi{}, the form the rest of the package already used. Fixed in performStudentizedBPTest, performSzroeterTest and performWhiteTestBootstrap, and in the roxygen sources so a regeneration produces the same form.
• README no longer lists two limitations that were fixed in 0.8.0: the renv.lock drift and the bundled dataset.

0.8.0

The API review that the validation effort was sequenced towards. Six exports are removed. This is a breaking change, hence the minor version.

Three returned nothing but a migration error, having been withdrawn once their statistics were shown not to work:

removed	use instead
performRiceTest()	performSzroeterTest(), performGQTest()
performCurryWalshTest()	performSpatialHeteroTest()
performHCCovarianceTest()	performBPTest(), performKoenkerTest(), performWhiteTest(); sandwich::vcovHC() for robust covariance

Three computed a valid statistic under a name that promised a different method, duplicating a test that remains:

removed	identical to
performOrderedLMTest()	performKoenkerTest() -- sorting the rows and refitting returns the same residuals permuted, so its order_by argument could not affect the result
performCameronTrivediTest()	the e^2 ~ yhat + yhat^2 auxiliary F test, covered by performWhiteTest() and performNCVTest(); it was never Cameron and Trivedi's information-matrix test
performModifiedBartlettTest()	performBartlettTest() -- the correction factor its documentation called a modification is part of the standard definition

The package has not been released on CRAN, so no published contract is broken. tests/testthat/test-public-api.R records the removals and checks that every replacement named above is still exported.

Two dead fallback routes went with them: hc_covariance was registered as a recovery path for quantile_regression and high_dimensional, and hc_covariance_hc0 as a custom fallback, both calling a function that only ever raised an error.

Packaging

• The pkgdown workflow now deploys. It called build_site_github_pages(), which writes docs/, and then stopped: there was no deploy step and no gh-pages branch, so the job reported success while nothing was published. That is why the URLs in DESCRIPTION and inst/CITATION returned 404.
• renv.lock now records all 38 declared non-base dependencies, up from 14. renv::restore() previously pinned about a third of them, and the Docker build only worked because remotes::install_local() resolves the remainder from CRAN at whatever version is current. The R version moves to 4.5.1, which is where the recorded versions come from.

0.7.2

Pass C of the statistical validation matrix: the methods with the least reference coverage. Six were sound and one needed a fix. Two are withdrawn because their statistics cannot detect heteroscedasticity at all, and two more are valid statistics under names that promise something else.

Withdrawn

• performRiceTest() compared Rice's (1984) difference-based variance estimator with the mean squared residual and referred the ratio to an F distribution. For independent residuals E[(e_i - e_{i-1})^2] is s2_i + s2_{i-1}, so the numerator estimates the mean of the variances and so does the denominator: the ratio sits at one under any variance pattern. Simulated means were 1.00 under homoscedasticity, under sigma = 0.2 + 1.2x, under sigma = exp(x) (a fiftyfold spread) and under a twentyfold step change, and empirical size was 0.000. Adding an ordering argument does not repair it. Use performSzroeterTest() or performGQTest() instead.

• performCurryWalshTest() passed an unstandardised Moran's I to 2 * (1 - pnorm(abs(I))). Moran's I has null expectation -1/(n-1) and a weights-dependent variance, and is bounded near one, so it essentially never reached 1.96: it rejected in 0.0% of samples under strong heteroscedasticity and returned p = 0.875 on clustered variance. performSpatialHeteroTest() already does this correctly via spdep::moran.mc().

Both retain their exports and signal a migration error, following the withdrawn HC covariance precedent.

Corrected

• performDavidianCarrollTest() failed with NA/NaN/Inf in 'y' when a residual was exactly zero. It now uses the shared log-residual floor, as Harvey and Park do.

Documented rather than removed

• performOrderedLMTest() sorts the rows and refits, which returns the same residuals permuted, so order_by cannot change the statistic and the result is identical to performKoenkerTest(). It now warns that the argument has no effect.
• performCameronTrivediTest() is exactly the e^2 ~ yhat + yhat^2 F test, not the information-matrix test of Cameron and Trivedi (1990), which concerns overdispersion in count models. The help page now says what it computes.

Removing either redundant export is left to the API review.

Validated without change

performRankPermutationTest(), performHighDimensionalTest(), performWildBootstrapTest(), performQuantileRegressionTest(), performSpatialHeteroTest() and performPesaranTest() hold their nominal level and have power.

0.7.1

Pass B of the statistical validation matrix: the group-variance tests. This section was missing from the changelog when 0.7.1 shipped and is recorded here.

Statistical corrections

• Bug fix (critical): performOBrienTest() rejected on every input. O'Brien's transformation produces one score per observation, whose group mean is that group's variance. The implementation computed a single value per group and repeated it across that group's observations, so the scores had no within-group variability: the one-way ANOVA had a residual sum of squares of zero, giving F on the order of 1e30 and p = 0 for any data, homoscedastic or not. Empirical size is now 4.8%.

• Bug fix (critical): performHartleyFmaxTest() used the wrong null distribution. The statistic is the largest of k group variances over the smallest, but the p-value came from pf(F, df, df), the distribution of a single variance ratio; the number of groups was computed and never used. Empirical size at a nominal 5%, groups of 30, was 0.0998 at two groups, 0.2238 at three, 0.3497 at four and 0.5626 at six. It now uses the maximum-F-ratio distribution and holds 0.05 throughout.

• performOBrienTest() also called safe_lm() without a data argument, so every call errored once the transformation was corrected.

Validated without change

performLeveneTest(), performBrownForsytheTest(), performBartlettTest() and performFlignerKilleenTest() reproduce car::leveneTest() with mean and median centring, stats::bartlett.test() and stats::fligner.test() to within 1e-8.

Documentation

• performModifiedBartlettTest() is Bartlett's test under another name: the correction factor its documentation described as a modification is part of the standard definition. It is now documented as a compatibility alias.
• Bartlett's per-group minimum was relaxed from three observations to two, which is what a sample variance requires.

0.7.0

This release is the first pass of a package-wide statistical validation effort. Pass A covers the classical regression diagnostics. Every test in this group was checked against an established implementation where one exists, and otherwise against a reconstruction of the primary reference, followed by Monte Carlo size and power measurement. Reference equivalence is asserted in tests/testthat/test-pass-a-reference.R; the simulation evidence is reproducible via inst/validation/pass-a-size-power.R and recorded in inst/validation/pass-a-size-power.csv.

Several of the corrections below change the value a test returns. They are breaking changes, and the minor version has been bumped accordingly.

Statistical corrections

• Bug fix (critical): performSzroeterTest() was standardised incorrectly and had no power. Szroeter's (1978) rank-weighted statistic h = sum(i * e_(i)^2) / sum(e_(i)^2) has null mean (n + 1) / 2 and null variance (n^2 - 1) / (6n). The previous implementation rescaled h by the mean rank and then divided the centred statistic by a further sqrt(n), using an ad-hoc variance (n + 1)(2n + 1) / (12n). The net effect was to shrink the standardised statistic by a factor of roughly 2 / sqrt(n) — a factor of 5 at n = 100. In 5,000 replications the test rejected in 0.0% of samples under the null and under every alternative tried: it could not reject at any sample size or effect size. The statistic is now Szroeter's Q, which attains 4.5% size and 87–91% power in the same designs. The test also gained an alternative argument ("greater", the default, for variance increasing with the ordering variable), and reports h in estimate.

• Bug fix: performNCVTest() did not compute the NCV test. It regressed the absolute residuals on the fitted values and reported a t statistic — a Glejser-type diagnostic — while documenting itself as equivalent to car::ncvTest(). It now computes the Cook–Weisberg score statistic: half the explained sum of squares from regressing e^2 / (RSS / n) on the variance model, referred to a chi-squared distribution. It reproduces car::ncvTest() to within 1e-8. A new var_formula argument selects a variance model other than the fitted values, matching car::ncvTest(var.formula = ).

• Bug fix: performCookWeisbergTest() returned the Koenker statistic. It computed n R^2 from regressing the raw squared residuals on the fitted values — that is, Stata's estat hettest, iid fitted — while documenting itself as Stata's estat hettest, fitted. It now computes the Cook–Weisberg score statistic and delegates to performNCVTest() so the two exports cannot drift apart. Note the trade-off this exposes: the score test assumes normal errors, and in simulation it over-rejects under heavy tails where the previous (mislabelled) statistic did not. performKoenkerTest() remains the kurtosis-robust choice and is now cross-referenced from both help pages.

• Change: performHarveyTest() now follows Harvey (1976). The variance regressors default to the model's own explanatory variables rather than the fitted values and their square, and the statistic is the classical ESS / (pi^2 / 2) referred to a chi-squared distribution, where pi^2 / 2 is the null variance of log(chi^2_1). The previous behaviour is available as auxiliary = "fitted", and the studentized F form — which estimates that variance from the data instead of assuming normality, and holds size better under heavy tails — as studentize = TRUE.

• Harvey and NCV now take their degrees of freedom from the realised rank of the auxiliary design, so a rank-deficient variance model is not credited with degrees of freedom for columns that lm() aliased away.

• performHarveyTest() and performParkTest() now warn when a residual is numerically zero and has been floored before taking logarithms. A floored observation contributes log(.Machine$double.eps) ~ -36, roughly sixteen standard deviations below the mean of the null log(chi^2_1) distribution, and can dominate the auxiliary regression on its own.

• performParkTest(), performGlejserTest(), performNCVTest() and performCookWeisbergTest() now name their parameter element, so print.htest() labels the degrees of freedom instead of printing = 118.

Validated without change

performArchLMTest() and performMcLeodLiTest() were audited and required no change. ARCH LM reproduces Engle's T R^2 auxiliary regression with T = n - q (the FinTS::ArchTest() convention) and McLeod–Li reproduces the Ljung–Box statistic on squared residuals with df = m, both to within 1e-10, with empirical size of 4.4% and 5.6% and power of 99% and 97% against an ARCH(1) alternative.

Diagnostic workflow

• Bug fix: compareModelDiagnostics() could report one test's statistic under another test's name. When a requested diagnostic fails, runHeteroTests() substitutes a fallback diagnostic and tags the result with the substitute's name. compareModelDiagnostics() read the statistic without checking that tag, so a comparison table column headed white could contain, for example, an NCV statistic. It now reports NA and warns when a substitution occurred, since a numeric table has nowhere to carry the caveat.

  This masking is also why two existing tests changed expectations. Both were passing because the previously unvalidated performNCVTest() succeeded on degenerate models — including a perfectly explained fit — and so served as a fallback that always worked. With NCV validating its input, the underlying failures are now surfaced rather than papered over.

Source integrity

• Removed two shadowed duplicate definitions. performHCCovarianceTest() and performQuantileRegressionTest() each had two definitions in the source tree. The obsolete implementations in R/modern_diagnostics.R were loaded first and then silently overwritten by the corrected versions in R/zz_deprecated_hc_covariance.R and R/quantile_regression_joint.R. Runtime behaviour was correct only because of collation order; a partial sourcing, a file rename or a collation change would have resurrected the withdrawn HC pseudo-test and the superseded quantile comparison. The obsolete definitions and their roxygen blocks are gone, along with the now-unreachable .hc_adjustment() helper whose only caller was the deleted HC implementation.

Resampling inference

• Bootstrap and permutation p-values now consistently use the finite-simulation convention (1 + #{T_b >= T_obs}) / (B_eff + 1), where B_eff counts the replicates that converged. performWhiteTestBootstrap(), performWhiteTestEnhanced() and the shared rbootstrap_test_statistic() engine previously used the raw proportion #/B, which can return exactly zero — not an attainable p-value from a finite number of replicates — and is anti-conservative in the tail. performWildBootstrapTest() and performRankPermutationTest() already used the corrected form. tests/testthat/test-bootstrap-pvalue-convention.R pins all of them.

Packaging and continuous integration

• Fixed the vignette build, which was breaking R CMD build and therefore every platform's check job. Three vignettes failed to build: bibliography.Rmd called algorithms_bibliography(), which has never existed as a function (the name belongs to a documentation-only help topic); comprehensive_guide.Rmd called simulate_hetero(n = 50) without the required beta0, beta1 and sigma_func arguments; and using_heteroTests.Rmd passed mtcars as the data for a model fitted on boston_housing, referenced the non-existent performGoldfeldQuandtTest(), contained an illustrative chunk using an undefined your_data placeholder, and had a malformed chunk fence that left ten lines of code outside any chunk. R CMD build now completes and all eight vignettes render.
• inst/CITATION modernised: 25 citEntry() calls converted to bibentry() (deprecated since R 4.2), a duplicate citHeader() that silently replaced the first removed, and the package version now read from meta$Version. The previous utils::packageDescription("heteroTests")$Version raised "$ operator is invalid for atomic vectors" during CRAN incoming checks when the package was not installed.

• The CI R-version matrix was rebuilt around R-devel, R-release and R-oldrel-1, with macOS and Windows on R-release, replacing the fixed 4.0/4.1/4.2/4.3 set. The minimum version declared in Depends: now has its own Ubuntu 22.04 job that installs hard dependencies only, so it reports on the package's own R 4.0 compatibility rather than on whether modern Suggests still support R 4.0.

• Fixed a broken example that failed R CMD check. performQuantileRegressionTest() documented itself with mtcars, which has 32 observations, while the test requires at least 40. The example now uses boston_housing.

• Non-ASCII characters removed from R/messages.R and R/validation.R. The two occurrences inside string literals are now written as the escapes \u00b2 and \u2022; the en dashes in roxygen comments became ASCII hyphens.
• Resolved 18 duplicated \alias entries. man/simulate_hetero.Rd is a complete hand-written family page for simulate_hetero(), simulate_arch1() and every sigma_*() helper, but roxygen2 also generated a thinner page per helper, so each name was documented twice. The generated stubs are removed — several were incomplete, for example sigma_poly.Rd omitted a, b and c — and the corresponding roxygen blocks are marked @noRd so they do not reappear. The same applied to generate_benchmark_report, which is documented on the run_benchmark_suite page.
• Reconciled \usage sections that had drifted from their function signatures: performWhiteTest() (max_interactions), plot.HeteroDiagnostic() (plots) and summary.HeteroDiagnostic() (tests).
• Documented previously undocumented arguments in generateHeteroRecommendations.Rd, plot.power_analysis.Rd, print.htest_enhanced.Rd, test.HeteroDiagnostic.Rd and simulate_hetero.Rd, and merged the duplicated n entry on the last of these.

• .Rbuildignore now excludes .zenodo.json, CITATION.cff and CONTRIBUTING.md, which are repository metadata rather than package content. inst/validation/ is deliberately shipped so the size and power evidence travels with the package.

• Fixed a latent runtime failure in the dashboard. launchDiagnosticDashboard() used the magrittr pipe %>% in four places without importing it. DT and plotly are Suggests and are called with :: rather than attached, so the pipe was only ever resolved if one of them happened to be attached by something else. The calls are now nested.

• A second example failed R CMD check: performWhiteTestBootstrap() used mtcars (32 observations) while the bootstrap requires 50. It now uses a simulated 120-observation dataset, with the heavier replication counts moved under \donttest{}.
• Fixed 15 broken \link[stats:htest]{htest} cross-references. There is no htest topic in stats, so every one of these was a dangling link; they are now plain \code{htest}.
• Completed the NAMESPACE imports. AIC, cooks.distance, cov, lm.fit, model.frame, model.response, nobs, quantile, resid, terms and uniroot from stats, and head, object.size and tail from utils, were used without being imported. Column names referenced inside aes() are declared through utils::globalVariables().

• Removed the Author: and Maintainer: fields from DESCRIPTION, which had drifted from the values derived from Authors@R (the derived form carries the ORCID). Unescaped 15 \& sequences in .Rd files, and added _pkgdown.yml and LICENSE to .Rbuildignore.

• Fixed three defects in the benchmarking system, all of which surfaced once its example was actually reachable under R CMD check --run-donttest. run_benchmark_suite() aborted with "no rows to aggregate" whenever every measured difference was NA, because aggregate() drops NA responses and the guard only tested for zero rows. generate_benchmark_report() then failed with "argument 1 is not a vector" when called with profile_memory = FALSE, since it passed a NULL column to order(). And ht_extract_metrics() had no branch for car::ncvTest(), which returns ChiSquare/Df/p rather than the htest field names, so every accuracy row measured against car was silently NA. With that branch in place the benchmark now reports a p-value difference of about 1e-15 between performNCVTest() and car::ncvTest() -- an independent confirmation of the correction described above.

• Fixed the Breusch-Pagan benchmark baseline. performBPTest() is the classical statistic, but the benchmark compared it against lmtest::bptest() with its default studentize = TRUE, i.e. the Koenker form. The pairing is now studentize = FALSE. With this and the car::ncvTest() extractor in place, all three baselines agree with our implementations to machine precision (0, 1.7e-16 and 1.2e-15), so the benchmark now doubles as a live reference-equivalence check.

• Fixed simulate_power_analysis(), whose effect size had no effect. It built the variance function as effect_size * sigma_func(x), but heteroscedasticity tests are invariant to multiplying every standard deviation by a constant, so every effect size produced statistically identical data and the reported power differed only by Monte Carlo noise. The effect size now interpolates between a constant variance and the supplied pattern, so 0 is homoscedastic and 1 is the pattern itself. Power against sigma_linear with White's test now runs 0.04 / 0.07 / 0.33 / 0.96 across effect sizes 0 / 0.1 / 0.4 / 1.0, recovering the nominal level at zero effect.

• The Type I error test now sizes its tolerance from the Monte Carlo standard error rather than a fixed 0.03 margin, which at 200 replications was under two standard errors and failed by chance in roughly one run in seven.

• Raised the declared minimum R version to 4.1. Depends: R (>= 4.0.0) was unsatisfiable: ggplot2 and scales are hard Imports and both now require R >= 4.1, so the package could not be installed on R 4.0 whatever its own code did. The minimum-version CI job checks 4.1 accordingly.

• Fixed two CI jobs that had never run successfully. The spell check called spelling::spell_check_package(ignore = ), an argument that does not exist, so the step errored every time; inst/WORDLIST is read via use_wordlist instead. The benchmark job gated on heteroTests being no more than 50% slower than bare lmtest/car calls, which it cannot satisfy -- it runs a validation layer those do not, and at n = 100 the baseline time rounds to zero, making the ratio infinite. The timing table is now reported as an artefact and the job gates on accuracy, which is the property that must not regress.

• Fixed the pkgdown deploy, which had been failing on main since before this release. man/algorithms_bibliography.Rd carried \usage{NULL}, which is what roxygen emits for a documentation-only topic and which pkgdown cannot parse. A topic that documents no function needs no usage section. This is why the package site was never published, and so why the URLs in DESCRIPTION and inst/CITATION returned 404.

• The automated release workflow now installs dependencies before R CMD build, which renders the vignettes and therefore needs the vignette builder. That step is gated on the version changing, so the omission had never fired until the 0.7.0 bump.

• Fixed generateDiagnosticReport(), which could not generate a report. The YAML front matter it writes put two mappings on one line (output: html_document:), so rmarkdown aborted with "mapping values are not allowed in this context". html_document is now nested under output. The example is wrapped in \dontrun{}, so only a check passing --run-dontrun executed it, which is why it went unnoticed.

• The CRAN submission job now installs a LaTeX toolchain, and the pkgdown job installs rsconnect. Both were failing on missing dependencies rather than on anything in the package: the former reported pdflatex is not available, which R CMD check surfaces as a probable Rd problem, and the latter needs rsconnect to inspect inst/tutorials.

Documentation

• performGlejserTest() documents that the test is not asymptotically valid under asymmetric errors (Godfrey 1996; Im 2000), and that scanning several transformations and reporting the smallest p-value invalidates the nominal level.
• performParkTest() documents that its t statistic is valid only asymptotically, the auxiliary error being a strongly skewed log(chi^2_1) variate, and that a non-significant result speaks only to the single variable tested.

0.6.5

• Bug fix (statistical correctness): performBPTest()/performBreuschPaganTest() now compute the classical Breusch-Pagan statistic (half the explained sum of squares from the scaled squared residuals), matching lmtest::bptest(studentize = FALSE). It previously returned the studentized n R^2 statistic despite being documented as the classical, normality-based test.
• Bug fix (statistical correctness): performKoenkerTest() now implements Koenker's (1981) studentized Breusch-Pagan statistic as n R^2 from regressing the squared residuals on the regressors, matching lmtest::bptest(studentize = TRUE) and performStudentizedBPTest(). It previously regressed the absolute residuals, which is not a recognized test and produced p-values without a valid chi-squared justification.
• The streaming variants performBPTestStreaming() and performKoenkerTestStreaming() were updated to stay algebraically equivalent to their exact counterparts.
• Bug fix: resolved a duplicate %||% operator definition whose two copies had conflicting semantics. Depending on collation order this could cause runHeteroTests() to abort with a "length > 1" coercion error when given a multi-column data frame. Both copies now share a single vector-safe definition.
• Bug fix (statistical correctness): performWildBootstrapTest() now generates its bootstrap reference distribution under the homoscedastic null by resampling leverage-standardised, centred residuals i.i.d. before applying the wild multiplier. The previous implementation multiplied the original (heteroscedastic) residuals by the wild weights, so the bootstrap distribution inherited the heteroscedasticity being tested for and the test had essentially no power (returning p ~ 0.5 on strongly heteroscedastic data). It now controls size and detects heteroscedasticity reliably; regression tests cover both.
• Packaging: R6 moved from Suggests to Imports (it is used at package load time), parallel added to Imports (used via ::), and the optional digest dependency (used for cache hashing) is now declared in Suggests. These resolve R CMD check "undeclared dependency" and load-time failures.
• Bug fix: performWhiteTest() now drops perfectly collinear auxiliary columns (via QR pivoting) and reports df equal to the realised rank, instead of aborting. White's test consequently works on models containing factors, where a squared dummy column is collinear with the dummy itself.
• Bug fix: performVIFDiagnostic() now computes variance inflation factors from the inverse correlation matrix of the design, so it works for models with factor predictors (it previously built auxiliary formulas from contrast-expanded column names that do not exist in the model frame, e.g. cyl6, and errored).
• ht_data_cleaning_suggestions() now flags infinite values found in the supplied data, not only those mentioned in the error message.
• Bug fixes surfaced by R CMD check: corrected invalid namespace qualifications (stats::all.vars, stats::summary, utils::close are base functions), fixed a suggestRemediation() call in generateDiagnosticReport() that passed unsupported arguments, and declared the (guarded) optional mgcv dependency.
• Data sets boston_housing, diagnostic_data and hetero_data are now shipped as .rda files (generated by scripts under data-raw/) rather than R scripts, removing the LazyData note and making the bundled data reproducible.
• Testing: the suite is consolidated onto a single test-*.R naming convention, pins Config/testthat/edition: 3, drops the deprecated context() calls, removes three stale top-level-skip()ped files that provided no coverage, and adds behavioural regression tests for the fixes above.

0.6.4

• Added "Comprehensive Guide to Heteroscedasticity Testing" and "Troubleshooting Common Issues" vignettes.
• Implemented performWhiteTestStreaming for large datasets.
• Added wild bootstrap, HC covariance, quantile regression, permutation, high-dimensional, and spatial heteroscedasticity diagnostics.
• Integrated broom tidiers, ggplot2 theming/autoplot support, tidymodels workflows, survey designs, grouped data pipelines, and data.table/dtplyr compatibility for diagnostic workflows.
• Introduced an intelligent recommendation engine that profiles datasets, interprets diagnostics with confidence levels, and delivers remediation guidance plus decision-tree style reports for non-statisticians.
• Implemented robust error handling with adaptive fallbacks, actionable cleaning suggestions, and configurable logging for complex diagnostic workflows.

0.6.3

• Added automated report generation via generateDiagnosticReport.

0.6.2

• Added studentized Breusch-Pagan, bootstrap White and Szroeter tests.

0.6.1

• Introduced TestFactory R6 class with metadata support.
• Added exported test_factory instance registering core tests.
• DESCRIPTION now lists R6 in Suggests.