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.
data: A base::data.frame containing the variables used to fit model and the ordering variable named by order_by. When NULL, the default, the model frame of model is used, so order_by must then be a variable of the model. The model frame holds the terms of the model rather than its variables (log(x) rather than x), so a model that transforms a variable needs data, as does a model with an offset() term.
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.

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.

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.

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"
)