Szroeter test for ordered alternatives¶

Source: R/additional_tests.R · Help: ?performSzroeterTest

Detects monotonic changes in variance when observations can be meaningfully ordered—typically by time or another covariate—using the cumulative weighting scheme proposed by Szroeter (1978).

Usage¶

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

Arguments¶

model: A fitted stats::lm object representing the mean equation under study.
data: The data frame model was fitted to, or NULL, the default. The residuals and weights are read from the fit itself, and order_by 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 that defines the ordering of observations prior to computing the statistic. The column must be numeric or coercible to an orderable vector.
alternative: Character scalar specifying the alternative hypothesis: "greater" (the default) tests for variance increasing with order_by, "less" for variance decreasing, and "two.sided" for a change in either direction. Szroeter's test targets monotone alternatives, so the one-sided form is the usual choice.

Value¶

An object of class htest reporting Szroeter's standardised statistic Q, the sample size, the underlying rank-weighted statistic h in estimate, and the p-value from the normal reference distribution, with the null mean and variance of h given the design.

Details¶

After ordering the residuals \(\hat{e}_{(i)}\) by order_by, the test forms the rank-weighted average of the squared residuals

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

which is Szroeter's (1978) class of statistics evaluated at the canonical weights \(h_i = i\). Variance that grows with order_by shifts weight onto the high ranks and drives \(h\) upwards.

Since 0.16.0 the null mean and variance of \(h\) are those of the design. Let \(M\) be the residual-maker of the design in the order of the test (the rows \(\sqrt{w_i}\, x_i\) on a weighted fit), \(A = \mathrm{diag}(1, \ldots, n)\), \(k\) the rank of the design and \(\nu = n - k\). Under normal errors of constant variance (given the weights) the residuals make \(h\) a ratio of quadratic forms, \(h = \sum_j \lambda_j z_j^2 / \sum_j z_j^2\), with \(z_j\) independent standard normal and \(\lambda_j\) the \(\nu\) eigenvalues of \(MAM\) on the residual space. The ratio is independent of its denominator, so its moments are ratios of moments:

\[ E\,h = \frac{\mathrm{tr}(MA)}{\nu} = \frac{n (n + 1) / 2 - \sum_i i\, l_{(i)}}{\nu}, \qquad \mathrm{Var}\,h = \frac{2\, (\nu\, \mathrm{tr}(MAMA) - \mathrm{tr}(MA)^2)}{\nu^2 (\nu + 2)}, \]

with \(l_{(i)}\) the leverages in the order of the test and \(\mathrm{tr}(MAMA) = \sum_i i^2 - 2 \sum_i i^2 l_{(i)} + \|Q_1' A Q_1\|_F^2\), \(Q_1\) an orthonormal basis of the design, so that no eigenvalue is computed and the cost is \(O(n k^2)\). The test refers

\[ Q = \frac{h - E\,h}{\sqrt{\mathrm{Var}\,h}} \]

to the standard normal distribution. The moments are derived here, by the argument that gives those of the Durbin-Watson statistic, and are not taken from Szroeter (1978). Without regressors they are \((n + 1) / 2\) and \((n^2 - 1) / (6 (n + 2))\). Before 0.16.0 the test used \((n + 1) / 2\) and \((n^2 - 1) / (6n)\), Szroeter's asymptotic moments, whatever the design: on a weighted fit whose leverage falls along the ordering the mean of \(h\) is above \((n + 1) / 2\) (by about 0.44 in the design of the validation below, at every \(n\)), and the test rejected too often. The statistic \(h\) is unchanged; \(Q\) and the p-value change for every fit, on lm(mpg ~ wt + qsec, mtcars) ordered by wt from -0.6494 and 0.7420 to -0.6626 and 0.7462.

The test needs at least two residual degrees of freedom, \(\nu \ge 2\). With one the residuals lie on a single line, \(h\) takes the same value whatever the errors and its null variance is 0: since 0.16.0 the test refuses such a fit, for which it returned a value before.

The residuals and weights are those of the fit, and order_by is read from the rows of data that hold its observations (the model context, since 0.16.0); the package validation helpers check that the ordering variable is present, and ensure sufficient sample size and variability in the squared residuals.

Ties in order_by keep the order of the rows the model was fitted to. Unless the rows were sorted by the response, that order does not depend on the residuals, so the test keeps its level under the null hypothesis; but with a discrete ordering variable the statistic depends on the order of the rows.

runHeteroTests(tests = "szroeter"), and anything else that runs the registered diagnostic, orders the observations by the fitted values of model, so its alternative is variance that rises with the mean. With a single regressor that is the order of the regressor when its slope is positive and the reverse order when it is negative. The fitted values are computed from the coefficients, so observations with the same regressors tie exactly, as they do under a factor or a binary regressor, and the ties are broken as above. Since 0.16.0 they are read from the fit, in its order, and data is matched to the observations of the fit as this function matches it: without data the registered test reads the model frame, and data in the order of the fit whose row numbers are not those of the fit are taken. Before, it matched the fitted values to the row names of data, and refused both ("Missing rows").

Validation¶

The statistic and its null moments are checked against an independent reconstruction from the full residual-maker in tests/testthat/test-pass-a-reference.R, and the moments against the eigenvalues of \(MAM\) to 1e-13. Releases before 0.7.0 divided the centred statistic by a further \(\sqrt{n}\), shrinking \(Q\) by a factor of roughly \(2 / \sqrt{n}\) and leaving the test with essentially no power against any alternative; see NEWS.md.

inst/validation/szroeter-design-null-size.R tested the samples of inst/validation/harrison-mccabe-size-power.R, lm(y ~ x1 + x2) ordered by x1, at 5000 replications per cell. With Gaussian errors the test rejected a true null hypothesis at the 5% level 0.048, 0.050 and 0.049 of the time at 40, 100 and 400 observations, and 0.049, 0.047 and 0.050 on a correctly weighted fit (weights \(1 / x_1^2\)), all inside the release gate's interval, 0.042 to 0.058; with the moments of 0.15.0 it rejected 0.044, 0.048 and 0.049, and on the weighted fit 0.067, 0.056 and 0.055. Its moments agreed with those of 100,000 simulated samples of one design of each size, variances within 1%. The exact null distribution of \(h\) given the design, from the eigenvalues by Imhof's (1961) method, gave the same rates within 0.0006, at a cost of \(O(n^3)\) per call. Under \(t_5\) errors no version holds the level (0.10 to 0.15): the null assumes normal errors, and performGlejserTest() with robust = TRUE is the test to use there.

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. 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. The null moments of a gam are those of its design, with the rank of that design, as performHarrisonMcCabeTest() takes its null distribution. Which model each test accepts, and why, is set out in test_results, under Supported models.

References¶

Szroeter, J. (1978). A class of parametric tests for heteroscedasticity in linear econometric models. Econometrica, 46(6), 1311–1327. doi:10.2307/1913831

Godfrey, L. G. (1988). Misspecification Tests in Econometrics. Cambridge University Press. Section 5.4 outlines the Szroeter test.

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)
performSzroeterTest(mod, mtcars, order_by = "wt")

# Detect ordered heteroscedasticity in simulated data with variance increasing
# in an index variable
set.seed(404)
n <- 150
x <- sort(runif(n))
y <- 2 + 0.5 * x + rnorm(n, sd = 0.4 + 0.8 * x)
df <- data.frame(y, x)
performSzroeterTest(lm(y ~ x, data = df), df, order_by = "x")

# Two-sided version when the direction of the change is not known in advance
performSzroeterTest(mod, mtcars, order_by = "wt", alternative = "two.sided")