Perform White's test for heteroscedasticity¶

Source: man/performWhiteTest.Rd · Help: ?performWhiteTest

Implements White's test on a fitted linear model.

Usage¶

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

Arguments¶

model: an object of class lm.
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.
cross_products: Logical. Include cross-product terms in the auxiliary regression?
max_interactions: Maximum number of cross-product terms admitted to the auxiliary regression. Guards the auxiliary design against growing quadratically with the number of regressors.

Value¶

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

Details¶

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

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

Degrees of freedom¶

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

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

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

Too few observations¶

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

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

Since 0.18.0 a fit with fewer terms than observations is refused only when the design forces it, as on the paired design: the test pushes a fixed pseudo-response through the QR decomposition of the fit and refuses only when the same terms fit the squares of its residuals exactly too, as Koenker's and the Breusch-Pagan tests decide it (see Exact fit in performKoenkerTest()); whether the regressors alone fit the squared residuals, which sets the advice and the fallback of runHeteroTests(), is asked of those squares too. From 0.16.0 to 0.17.0 the test refused a fit the data make exact as well, and it gives its value since 0.18.0: in a one-way layout of nine groups of ten observations and one of two, with a value of the group of two typed 100 times too large, 92 on 9 degrees of freedom; and on pairs with the same x whose residuals are 1 + x and -(1 + x), whose squares x and its square fit exactly, 80, the number of observations, on 2. A fit without a QR decomposition, such as an mgcv::gam() fit tested on its design, is refused as before, with a message that says only that the regression fits the squared residuals exactly.

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

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

Fixed effects for small groups¶

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

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

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

rrrrrl Subjects Observations of each White Koenker Breusch-Pagan Warns

10 10 6.4 6.2 7.3 no

20 20 6.0 6.4 8.0 no

5 4 1.0 5.4 7.0 yes

30 3 50.7 57.3 39.6 yes

30 4 22.2 23.9 28.7 yes

100 8 23.8 23.8 25.2 yes

100 10 18.1 18.4 21.3 yes

200 4 97.8 97.6 82.5 yes

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

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

Weighted fits¶

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

Supported models¶

It accepts an lm and a glm with the Gaussian family and the identity link, unweighted or with positive weights, a fitWLS() fit and, since 0.16.0, a rank-deficient fit, which it tests on the estimable part of its design: it gives the value of the same model fitted without the aliased columns. It refuses a glm of another family or link, a survey::svyglm() fit, a fit with a zero weight, a coefficient that is not finite for another reason, and data that lack rows or variables of the fit. Which model each test accepts, and why, is set out in test_results, under Supported models.

References¶

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

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

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

Examples¶

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