Perform Glejser test for heteroscedasticity¶

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

Regresses absolute residuals on a transformation of a suspected variable. By default, robust = TRUE since 0.15.0, the absolute residuals are first corrected for the asymmetry of the errors: this is the version of Im (2000) and of Machado and Santos Silva (2000), which keeps its level when the errors are asymmetric, from about 150 observations on (see the section on asymmetric errors for smaller samples). robust = FALSE gives Glejser's (1969) statistic, the default up to 0.14.0.

Usage¶

performGlejserTest(model, data = NULL, variable,
                   transformation = c("abs", "sqrt", "inverse", "inverse_sqrt"),
                   robust = TRUE)

Arguments¶

model: an object of class lm.
data: The data frame model was fitted to, or NULL, the default. The residuals and weights are read from the fit itself, and variable from the rows of data that hold the observations of the fit. Without data it is read from the model frame of model, which holds it only when it is a variable of the model under its own name (x, not log(x)); a model that transforms the variable needs data, as does a variable the model does not use. Since 0.16.0, 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.
variable: name of the suspected variable.
transformation: transformation applied to variable.
robust: logical. When TRUE, the default since 0.15.0, the absolute residuals are corrected for the asymmetry of the errors before the auxiliary regression. When FALSE the statistic is Glejser's (1969), which has its nominal level when positive and negative errors are equally likely, as they are for symmetric errors; it is what the function returned by default up to 0.14.0. See the sections on asymmetric errors and on the default.

Value¶

An object of class htest containing the t statistic, p-value and degrees of freedom. The method element says which of the two statistics was computed.

Details¶

A variety of transformations of the explanatory variable (e.g. absolute value, square root) can reveal a relationship between scale and the covariate. Significance of the slope parameter is assessed with a t test. For a weighted fit the residuals are the Pearson residuals.

Glejser's statistic, robust = FALSE, has its nominal level only when positive and negative errors are equally likely. The corrected statistic of the default does not need that, and is the one to use when the residuals may be skewed or their shape is unknown.

Either form conditions on one variable and one transformation. Scanning several transformations and reporting only the smallest p-value invalidates the nominal level. Simulated size and power of both statistics are recorded in inst/validation/glejser-skewed-errors.csv; those of Glejser's statistic in inst/validation/pass-a-size-power.csv, and those of the default in inst/validation/full-sweep-size-power.csv and inst/validation/weighted-fits-size.csv.

Asymmetric errors¶

Glejser's statistic takes the absolute residuals for the absolute errors. The two differ by the estimation error of the mean equation,

\[ |\hat{e}_i| \approx |e_i| - \mathrm{sign}(e_i)\, x_i^\top (\hat{\beta} - \beta), \]

and the second term averages out only when positive and negative errors are equally likely (Godfrey 1996). When they are not, and the variable tested is correlated with the regressors, the test does not in general have its nominal level, and a larger sample does not help. In the designs of inst/validation/glejser-skewed-errors.R Glejser's statistic, robust = FALSE, rejects a true null hypothesis about 13% of the time at the 5% level under exponential errors and about 17% under lognormal errors, at every sample size from 50 to 1000. For other error distributions the rate can also fall below the nominal level.

robust = TRUE replaces the regressand \(|\hat{e}_i|\) by

\[ |\hat{e}_i| - \hat{m}\, \hat{e}_i, \qquad \hat{m} = \frac{1}{n} \sum_{i=1}^{n} \mathrm{sign}(\hat{e}_i), \]

where \(\hat{m}\) is the proportion of positive residuals less the proportion of negative ones. This is the correction of Im (2000). The regressand of Machado and Santos Silva (2000), \(\hat{e}_i\,[I(\hat{e}_i \ge 0) - \hat{\eta}]\) with \(\hat{\eta}\) the proportion of non-negative residuals, is half of it when no residual is exactly zero, and the two then give the same test. A residual that is zero up to rounding error is counted as zero, so that the statistic does not depend on the sign of a rounding error.

The auxiliary regression is otherwise the same. Both references state the statistic as \(n R^2\) of that regression, referred to a chi-squared distribution with one degree of freedom. The t statistic reported here is the signed form of it, \(n R^2 = n t^2 / (t^2 + n - 2)\), and its p-value is taken from the t distribution with \(n - 2\) degrees of freedom, as for Glejser's statistic. The two p-values agree as the sample grows.

In the same study, with 5000 replications in each cell, the corrected statistic is within simulation error of the 5% level from 150 observations on: all 27 of its cells with a true null hypothesis and 150 observations or more are inside the 99% interval around 5%, from 4.21% to 5.79%, at 4.4% to 5.7%. At 50 observations it rejects 5.7% to 6.6% of the time under skewed errors, and four of those six cells are above that interval, at 5.80% to 6.58%. Under symmetric errors \(\hat{m}\) is close to zero and the two statistics differ little: the corrected one rejects 0.5 to 0.7 percentage points less often at 50 observations and gives up one or two points of power there. For a weighted fit the correction is applied to the Pearson residuals.

The default¶

Up to 0.14.0 the default was robust = FALSE. Since 0.15.0 it is robust = TRUE, so a call that does not give robust returns a different statistic and p-value from those it returned before, and its method ends in "(robust to asymmetric errors)". robust = FALSE returns what such a call returned up to 0.14.0.

Glejser's statistic needs positive and negative errors to be equally likely, which is seldom known: in the study above 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, where the corrected statistic rejects 4.6% to 5.7% of the time from 150 observations on. The change costs some power when the errors are symmetric, 1.1 and 1.9 percentage points at 50 observations under Gaussian errors in that study, and the corrected statistic is above the 5% level at 50 observations under skewed errors, as above.

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. It refuses a multivariate lm ("needs a model with one response"), since 0.16.0; before, it stopped with a coercion error. An mgcv::gam() fit is tested as before 0.16.0, with penalised smooth terms too: the test reads its residuals and does not fit it again, so it does not take them for least-squares residuals (decision 37 of the design of the model context); since 0.16.0 its observations are matched to the rows of data by the row names of its model frame. Before, they were matched by position, since the residuals of a gam carry no names, and data with the rows of the fit in another order gave a wrong value. Which model each test accepts, and why, is set out in test_results, under Supported models.

References¶

Glejser, H. (1969). A new test for heteroskedasticity. Journal of the American Statistical Association, 64(325), 316–323.

Godfrey, L. G. (1996). Some results on the Glejser and Koenker tests for heteroskedasticity. Journal of Econometrics, 72(1–2), 275–299. doi:10.1016/0304-4076(94)01722-0

Im, K. S. (2000). Robustifying Glejser test of heteroskedasticity. Journal of Econometrics, 97(1), 179–188. doi:10.1016/S0304-4076(99)00061-5

Machado, J. A. F. and Santos Silva, J. M. C. (2000). Glejser's test revisited. Journal of Econometrics, 97(1), 189–202. doi:10.1016/S0304-4076(00)00016-6

Examples¶

 data(mtcars)
 m <- lm(mpg ~ wt + qsec, data = mtcars)
 # The default, which keeps its level under asymmetric errors from about
 # 150 observations
 performGlejserTest(m, mtcars, "wt")

 # Glejser's (1969) statistic, the default up to 0.14.0
 performGlejserTest(m, mtcars, "wt", robust = FALSE)

 # Examine an alternative transformation
 performGlejserTest(m, mtcars, "wt", transformation = "inverse")