heteroTests News

0.11.2

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

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

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

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

Metadata

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

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

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

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

Corrected

Documented rather than removed

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

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

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

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

Source integrity

Resampling inference

Packaging and continuous integration

Documentation

0.6.5

0.6.4

0.6.3

0.6.2

0.6.1