Memory-efficient White test for large datasets¶

Source: R/performance.R · Help: ?performWhiteTestStreaming

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

Usage¶

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

Arguments¶

model: A fitted stats::lm object representing the mean specification to be diagnosed.
data: The data frame model was fitted to, or NULL, the default. The response, design, weights and offset are read from the fit itself, so a model that transforms its variables, has an offset() term, or takes its weights or offset from a column of the data needs no data, and the model frame of model is accepted as data (since 0.16.0; before, such models needed data). Data that are given are matched to the observations of the fit: rows with row names of their own by name, in any order; rows that only number themselves (a tibble, a data.table, a data frame whose row names were reset) by the row numbers of the fit, or in the order of the fit when they hold one row per observation. The rows the fit dropped may be there or not. Data whose rows cannot be shown to hold the observations of the fit are refused; see test_results.
chunk_size: Positive integer specifying the number of observations per chunk. Smaller values reduce memory usage at the expense of additional iteration overhead.
cross_products: Logical scalar indicating whether to include all pairwise cross-products of the regressors in the auxiliary regression. Defaults to TRUE and should remain enabled unless dimensionality makes the regression unstable.
max_interactions: Single positive integer giving the maximum number of original predictors for which cross-products are generated. When the number of regressors exceeds this threshold, cross-products are dropped to avoid explosive growth in columns. Defaults to 10.
progress: Logical flag indicating whether a progress bar should be displayed while streaming the data. Defaults to interactive().

Value¶

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

Details¶

The streaming implementation builds the auxiliary design, the regressors, their squares and their cross-products, chunk by chunk, without materialising it whole. Each chunk is reduced by a QR decomposition stacked on the factor of the chunks before it, which holds the cross-products of every row without forming them. The rank of the auxiliary regression and its $n R^2$ statistic are computed from that factor as in performWhiteTest(), ensuring numerical equivalence while reducing peak memory usage: the regressors are centred at their means, the rank is that of the auxiliary design at the default tolerance of qr(), and a term that is a combination of the intercept and the terms before it is dropped (see Degrees of freedom in performWhiteTest()). Since 0.18.0 the factor leaves out what is left of such a term when that is less than $10^{-13}$ of its norm: rounding error, for the square of a dummy variable. That changes the value of the test by rounding alone: by up to $2 \cdot 10^{-14}$ relative on 100,000 observations with a factor of four levels, less than the $2.2 \cdot 10^{-14}$ to $4.5 \cdot 10^{-14}$ by which the streaming test differs from performWhiteTest() there, before and since. What is left of the squared residuals, however little, is never left out. In 0.17.0 the decomposition made a reflection from that rounding error, the next such term left about its square, and so on down, and on a design with many such terms, as the squares and products of dummies are, it could divide by a norm below the reciprocal of the largest double and stop ("NA/NaN/Inf in foreign function call"), as it did on a one-way layout under a BLAS that rounds otherwise than the reference BLAS of R.

Since 0.17.0 the test drops such terms: the square of a dummy variable, which is the dummy, that of the first column of an orthogonal poly() basis, or that of x beside I(x^2), as performWhiteTest() drops them, and gives its degrees of freedom: 17.83 on 8 on lm(y ~ x1 + group) with a factor of four levels. Up to 0.16.0 it accumulated the cross-products of the uncentred auxiliary design, decided their rank by qr() of them, which squares the condition of the design, and stopped on such terms ("Auxiliary regression became rank deficient"), so on any model with a factor, and on terms that only looked collinear in those cross-products: the squares of a regressor 100 standard deviations from zero, or the terms of lm(y ~ x1 * x2) with x1 between 0.5 and 2. It took the squared residuals for constant when their total sum of squares about their mean was at most .Machine$double.eps, whatever their scale, and so refused residuals of about $10^{-5}$: on 120 observations of such residuals it now gives 12.787 on 5 degrees of freedom, the value of performWhiteTest(). It refuses squared residuals whose total sum of squares is at most .Machine$double.eps of their sum of squares, which agree to about eight digits. And in more than one chunk, on a regressor a beside a regressor a_sq, it put the square of a, which it names a_sq too, in place of that of a_sq in every chunk after the first: 27.76 in chunks of 7 where performWhiteTest() gives 27.64. It no longer uses the sparse cross-products of Matrix for large chunks.

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

Supported models¶

It accepts an lm and, since 0.16.0, a glm with the Gaussian family and the identity link, unweighted or with positive weights, a fitWLS() fit and, since 0.16.0, a rank-deficient fit, which it tests on the estimable part of its design: it gives the value of the same model fitted without the aliased columns. The Gaussian glm with prior weights w is the least-squares fit lm(..., weights = w), with the same design and residuals, and the test gives the value of the lm. It refused every glm before 0.16.0, because the checks it runs read the R-squared of summary(), which a glm does not report; they read that of the least-squares fit now. It refuses a glm of another family or link, a survey::svyglm() fit, a fit with a zero weight, a coefficient that is not finite for another reason, data that lack rows or variables of the fit, and, since 0.16.0, a multivariate lm, a fit changed after it was fitted and a fit with model = FALSE whose data are gone, unless data gives its response and fitted values (see test_results). Since 0.17.0 it drops the terms of its auxiliary regression that are collinear, as the squares of dummy variables are, as performWhiteTest() drops them, where it stopped. Since 0.16.0 it stops with the refusal of performWhiteTest() (see its section Too few observations) when its independent terms are at least as many as the observations, where it stopped as on collinear terms, or fit the squared residuals exactly, collinear or not, as on a paired design. Since 0.18.0 it refuses a fit with fewer terms than observations only when the design forces it, decided from the squared pseudo-residuals it carries in the same pass, as performWhiteTest() decides it; a fit the data make exact gives a value, where it was refused from 0.16.0. Which model each test accepts, and why, is set out in test_results, under Supported models.

References¶

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

Examples¶

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