Harrison-McCabe test for heteroscedasticity¶

Source: R/performHarrisonMcCabeTest.R · Help: ?performHarrisonMcCabeTest

Orders the observations by a variable and tests whether the share of the residual sum of squares that falls in the first part of the ordered sample is too small, or too large, for errors of constant variance. The p-value is exact under normal errors, given the regressors.

Usage¶

performHarrisonMcCabeTest(
  model,
  data = NULL,
  order_by,
  point = 0.5,
  alternative = c("less", "two.sided", "greater")
)

Arguments¶

model: A fitted stats::lm object, with or without weights, or a stats::glm with the Gaussian family and the identity link.
data: The data frame model was fitted to, or NULL, the default. The residuals, design and weights are read from the fit itself, and so, since 0.16.0, is its rank, decided at the tolerance of the fit (see the section on the rank), where the test decided the rank of the design again at 1e-7; order_by is read 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 it is ordered by needs data, as does an ordering 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.
order_by: Character scalar naming the column of data that orders the observations. Tied values keep the order the observations have in the fitted model.
point: The size of the first part of the ordered sample. A number strictly between 0 and 1 is a fraction of the sample, and the first part then holds floor(point * n) observations, as in lmtest::hmctest(). A whole number of 1 or more is the number of observations. Each part must hold at least as many observations as the model has coefficients.
alternative: Character scalar. "less" (the default) is the alternative that the variance increases with order_by, which makes the share small; "greater" that it decreases; "two.sided" either. The two-sided p-value is twice the smaller of the two one-sided ones, capped at 1, as in performGQTest() and in the default of skedastic::harrison_mccabe().

Value¶

An object of class htest with the share HMC as the statistic, the number of observations in the first part m and the sample size n as parameters, the p-value, and the alternative in words.

Details¶

With \(\hat{e}_{(1)}, \dots, \hat{e}_{(n)}\) the least-squares residuals in the order of order_by, the statistic of Harrison and McCabe (1979) is

\[ b = \frac{\sum_{i=1}^{m} \hat{e}_{(i)}^2}{\sum_{i=1}^{n} \hat{e}_{(i)}^2}, \]

the share of the residual sum of squares that falls in the first \(m\) observations. A variance that increases along the ordering makes \(b\) small, and one that decreases makes it large. It is the statistic of lmtest::hmctest().

The 1979 article was not consulted for this implementation. The definition follows Harrison (1981), by its first author, which restates the statistic with \(m\) the integer part of \(n / 2\) and rejects for small \(b\), and it agrees with lmtest::hmctest().

Under normal errors of constant variance, \(b\) is a ratio of quadratic forms in normal variables, and its distribution depends on the regressors. With \(M\) the residual-maker of the design in the ordered sample and \(D\) the diagonal matrix that selects the first \(m\) observations,

\[ P(b \le b_0) = P\left(\sum_j (\lambda_j - b_0) z_j^2 \le 0\right), \]

where the \(z_j\) are independent standard normal and the \(\lambda_j\) are the \(n - k\) eigenvalues of \(M D M\) on the residual space, \(k\) being the number of coefficients. They are found from the \(k\) singular values of the first \(m\) rows of an orthonormal basis of the design, and the probability is evaluated by Imhof's (1961) formula with stats::integrate(), as Harrison (1981) does for the exact variant of the test. For a constant and a linear trend, the 5% critical values this gives at 10, 20, 30 and 40 observations are those of Harrison's (1981) Table 3 to the three decimals printed there. Under the null hypothesis the mean of \(b\) is \((m - \sum_{i \le m} h_{(i)}) / (n - k)\), with \(h_{(i)}\) the leverages in the ordered sample.

The probability is computed to an absolute accuracy of about 2e-11. Smaller p-values are not resolved, and a statistic far in a tail can return a p-value of 0.

The p-value depends only on the space spanned by the regressors, so it does not change with the parametrisation of the model.

The rank¶

\(k\) is the rank of the fit, the number of its coefficients that are not aliased, decided at the tolerance the fit used: 1e-7 for an lm, or the tol it was given, and 1e-11 for a glm (since 0.16.0). Before, the rank was decided again at 1e-7, so a Gaussian glm, or an lm fitted with a smaller tol, whose design holds a column that is a combination of the others to within 1e-7 but not to within its own tolerance, was given a null distribution of a smaller rank than its residuals. With x3 = x1 + 1e-9 * rnorm(120) the glm has rank 4, and its p-value is 0.002273 where it was 0.002288, computed with rank 3. On every other fit the two ranks agree and the value is unchanged.

Direction of the alternative¶

In performGQTest() and performSzroeterTest() the statistic grows when the variance increases along the ordering, and in this test it shrinks. The alternative names the tail of each statistic, so the same hypothesis, a variance that increases with order_by, is alternative = "less" here and alternative = "greater" in those two functions. In all three it is the default, and the alternative element of the result states it in words.

Weighted fits¶

On a fit with weights the statistic is computed from the Pearson residuals \(\sqrt{w_i}\, \hat{e}_i\), and the null distribution from the weighted design \(\sqrt{w_i}\, x_i\). Both are those of the equivalent unweighted regression, so the test asks whether the weights are adequate. The weights of a Gaussian glm are its prior weights.

Heavy tails¶

The null distribution assumes normal errors. In the designs of inst/validation/harrison-mccabe-size-power.R the test rejects a true null hypothesis 11.7% to 16.3% of the time at the 5% level when the errors follow a \(t_5\) distribution, and more often the larger the sample; performGQTest() and performSzroeterTest() do the same. When the tails may be heavy, use performGlejserTest() with the ordering variable as variable and robust = TRUE, which in the same study rejects 4.7% to 4.9% of the time under those errors. Its default transformation takes the absolute value of the variable, which preserves the ordering only when the variable is positive.

Reference implementations¶

The statistic agrees with lmtest::hmctest() to 1e-8. Its p-value does not: with simulate.p = TRUE it is simulated from residuals of a regression on a constant alone, whatever the model, and with simulate.p = FALSE it is NA. skedastic::harrison_mccabe() reports a different p-value as well; inst/validation/README.md says why.

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): it has the same design, residuals and weights, and it states the normal errors of constant variance, given the weights, under which the null distribution above is exact, so the test gives the value of the lm when the two decide the rank alike (see the section on the rank). It refused every glm before 0.16.0. 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. Since 0.17.0 it refuses an mgcv::gam() fit with penalised smooth terms: the null distribution of the test is derived for the least-squares residuals of the design of the fit, the basis of its smooths, which the residuals of a penalised fit are not; before, it compared them with that distribution. A gam that is the least-squares fit of its design, as one whose smooth terms are all fixed (s(x, fx = TRUE)), is tested on that design, as before, and gives the value of the lm() of it (decisions 37 and 38 of the design of the model context). Since 0.16.0 the observations of a gam 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¶

Harrison, M. J. (1981). A comparison of the bounds, beta-approximate, and exact variants of two tests for heteroscedasticity based on ordinary least squares residuals. The Economic and Social Review, 12(4), 235–252.

Harrison, M. J., & McCabe, B. P. M. (1979). A test for heteroscedasticity based on ordinary least squares residuals. Journal of the American Statistical Association, 74(366a), 494–499. doi:10.1080/01621459.1979.10482544

Imhof, J. P. (1961). Computing the distribution of quadratic forms in normal variables. Biometrika, 48(3-4), 419–426. doi:10.1093/biomet/48.3-4.419

Examples¶

data(mtcars)
mod <- lm(mpg ~ wt + qsec, data = mtcars)
performHarrisonMcCabeTest(mod, mtcars, order_by = "wt")

# Variance that increases with x
set.seed(7)
x <- runif(80, 1, 5)
d <- data.frame(x = x, y = 1 + 2 * x + rnorm(80, sd = exp(0.3 * x)))
performHarrisonMcCabeTest(lm(y ~ x, data = d), d, order_by = "x")

# Either direction, with a third of the sample in the first part
performHarrisonMcCabeTest(lm(y ~ x, data = d), d,
  order_by = "x",
  point = 1 / 3, alternative = "two.sided"
)