Perform Koenker studentized Breusch-Pagan test¶
Source: man/performKoenkerTest.Rd · Help: ?performKoenkerTest
Implementation of the Koenker version of the Breusch-Pagan test.
Usage¶
Arguments¶
model- an object of class
lm. data- The data frame
modelwas fitted to, orNULL, the default. The response, design, weights and offset are read from the fit itself, so a model that transforms its variables, has anoffset()term, or takes itsweightsoroffsetfrom a column of the data needs nodata, and the model frame ofmodelis accepted asdata(since 0.16.0; before, such models neededdata). 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; seetest_results.
Value¶
An object of class htest containing the test statistic, p-value and degrees of freedom.
Details¶
Squared residuals are regressed on the regressors but the statistic \(n R^2\) uses a studentized form that is robust to non-normality. It is compared against a chi-squared distribution whose degrees of freedom are the rank of the auxiliary regression less its intercept, the number of regressors on a model with an 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 performKoenkerTestStreaming() on data above 100 MB; call that function to stream.
Degrees of freedom¶
The auxiliary regression of the squared residuals on the regressors \(Z\) has an intercept, and the degrees of freedom are the rank of \([1, Z]\) minus one, the rank of the auxiliary regression less its intercept, as White's test counts them (see Degrees of freedom in performWhiteTest()) and as lmtest::bptest() does. A column of \(Z\) that is a linear combination of the intercept and the columns before it is dropped. So a model with an intercept counts its regressors, and a model without an intercept whose columns span the constant, such as the dummies of every level of a factor in lm(y ~ 0 + group + x1), counts one fewer: the auxiliary regression adds an intercept that those columns already span. The rank is decided as in White's test, on the regressors centred one by one, at the default tolerance of qr(). A model whose only regressors are constant has none beyond the intercept, and is refused, as a model with an intercept alone is.
Up to 0.17.0 the degrees of freedom were the number of columns of \(Z\), one more on a model without an intercept whose columns span the constant, with the same statistic and a larger p-value: on lm(y ~ 0 + group + x1), with a factor of four levels and 120 observations, 10.52 on 5 degrees of freedom (p = 0.062), where lm(y ~ x1 + group), the same column space, and lmtest::bptest() give 10.52 on 4 (p = 0.032), as the test now does. performKoenkerTestStreaming() gives these degrees of freedom too, on such a model, which it refused up to 0.17.0 ("Koenker auxiliary regression became rank deficient"), and so does the registered "koenker" of runHeteroTests(). performBPTestRobust() with studentize = TRUE, the deprecated performStudentizedBPTest() and the registered "student_bp" counted the rank already. performWildBootstrapTest() reports these degrees of freedom since 0.18.0, where it reported the number of regressors; its p-value, a bootstrap's, does not use them. A model with an intercept, or without one whose columns do not span the constant, such as lm(y ~ 0 + x1 + x2), keeps its degrees of freedom. On the latter lmtest::bptest() given the model regresses the squared residuals on its columns without an intercept, and gives another statistic, as before; given the regressors as varformula = ~ x1 + x2 it gives the value of this test.
Up to 0.17.0 the auxiliary regression was also fitted by lm() from a formula, on the columns of \(Z\) as they are, and on three kinds of design of full rank its statistic was wrong. It is now the value of performKoenkerTestStreaming(), which had it already, and of lmtest::bptest() but on the Gaussian glm, whose design lmtest::bptest() refits by lm.fit(), dropping the regressor too; the registered "koenker" and performWildBootstrapTest(), whose statistic this is, change with it. On a model without an intercept with a regressor whose spread is less than about \(10^{-7}\) of its distance from zero, that fit dropped the regressor as a combination of its intercept, and the degrees of freedom still counted it: on the same 120 observations, lm(y ~ 0 + I(x1 + 1e8) + x2) gave 2.372 on 2 degrees of freedom (p = 0.305), where lmtest::bptest() given varformula = ~ x1 + x2 and the test now give 13.92 (p = 0.001). A Gaussian glm, fitted at a tolerance of \(10^{-11}\), keeps a regressor that lm() would alias: glm(y ~ I(1e6 + x1 / 10) + x2) gave 1.457 (p = 0.483) for 9.044 (p = 0.011), the value of lm(y ~ x1 + x2), the same column space. And a regressor named e_squared was taken for the response of the formula: lm(y ~ e_squared + x2) gave 0.014 (p = 0.993) for 9.044, the value of the same model under another name. On other designs the statistic is the one it was.
Exact fit¶
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. 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, the value of lmtest::bptest(), and on lm(y ~ x1 + f) with a factor of 10 levels of two observations, 20 on 10 (p = 0.029). Since 0.18.0 it refuses such a regression when the design forces it, where it returned that statistic, as White's test refuses it (see Too few observations in performWhiteTest()). The refusal is an error of class heteroTests_exact_fit_error, whose message says why and asks for more observations per level: the design forces the exact fit, and more levels of two observations do not help. A fit counts as exact by the rule of White's test, 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. The design forces it when the regression would fit the squares of the residuals of any response: the test pushes a fixed pseudo-response, standard normal draws from a fixed seed that leave the random numbers of the session as they were (but for the second value of a pair that the normal kind Box-Muller keeps, which R resets), through the QR decomposition of the fit, with its weights, and refuses only when the same regressors fit the squares of its residuals exactly too. A fit the data make exact is tested, as up to 0.17.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, the dummies fit the squared residuals to within \(1.1 \times 10^{-9}\) of their total sum of squares, an exact fit by the rule, and the test gives 92 on 9 degrees of freedom. A fit without a QR decomposition, such as an mgcv::gam() fit tested on its design, is refused, with a message that says only that the regression fits the squared residuals exactly. Up to 0.17.0 the test gave the warning of summary.lm() there, "essentially perfect fit: summary may be unreliable", with the value.
performKoenkerTestStreaming(), the registered "koenker", performBPTestRobust(), performStudentizedBPTest() and the registered "student_bp", which compute this statistic, refuse it too, and so do performWildBootstrapTest(), which bootstraps it, and bootstrapTestStatistic() of this test. runHeteroTests() takes no fallback, and stops with the refusal: its fallback, the Cook-Weisberg test, tests the fitted values, another hypothesis, and would return its value under this name. The other tests whose auxiliary regression can fit the squared residuals, or a function of them, exactly refuse such a fit by the same rule: the Breusch-Pagan test, the Cook-Weisberg test with a var_formula of the regressors, Harvey's test, the high-dimensional test when it keeps every principal component, and the variance form test.
The refusal covers the exact fit only. One step away from it, with one of the 30 subjects observed a third time, the test gives 59.75 on 30 (p = 0.001), and on any fit with fixed effects for small groups it rejects too often, and warns (see Fixed effects for small groups).
Fixed effects for small groups¶
Koenker's test, like White's and the Breusch-Pagan tests, is 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 it rejects 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 has mean its degrees of freedom plus \((G - 1) / (T - 1)\), as it had in the study below, within 2.3% for \(T\) of 3 or more. Regressed on x alone, without the dummies, the squared residuals were rejected in 4.6% to 5.5% of the draws (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. performKoenkerTestStreaming(), performBPTestRobust() with studentize = TRUE, the deprecated performStudentizedBPTest() and the registered "koenker" and "student_bp" give it too. performWildBootstrapTest(), which bootstraps the same statistic, 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.
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¶
lm(..., weights = w) states that the error variance is \(\sigma^2 / w_i\). On such a fit the test is computed from the Pearson residuals \(\sqrt{w_i}\, e_i\), the residuals of the equivalent unweighted regression, which have constant variance when the weights are right. The null hypothesis is then that the weights are adequate; the variance regressors stay on their original scale. Releases before 0.12.0 used the raw residuals, whose variance differs across observations by assumption, so a correctly weighted fit was rejected as often as an unweighted one. A fit with a zero weight is refused.
The reference distribution takes the weights as known. When they were estimated from the same data, as by fitWLS(), it is only approximate, and not more accurate in larger samples. With a correctly specified variance function and normal errors, the Koenker test rejected 8% of the time at the 5% level with 150 observations and 11% with 600, the White test 6% and 8%, and the Harvey test never, because it repeats the regression the weights were estimated from (inst/validation/weighted-fits-size.csv). performVarianceFormTest() accounts for the estimation, held its level in the same designs, and is the test to use on a fitWLS() fit.
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¶
Koenker, R. (1981). A note on studentizing a test for heteroscedasticity. Journal of Econometrics, 17(1), 107–112.
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