Perform Fligner-Killeen test for homogeneity of variances¶

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

The median-centred Fligner-Killeen test of homogeneity of variances, on the Pearson residuals of the fit, computed with fligner.test(): it ranks the absolute deviations of the residuals from their group medians and scores the ranks by quantiles of the normal distribution, which makes it robust to outliers and to errors that are not normal.

Usage¶

performFlignerKilleenTest(model, data = NULL, group)

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 group from the rows of data that hold the observations of the fit; the observations whose group is missing are left out, with a warning. 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 grouping variable the model does not use needs data. 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.
group: name of the grouping variable.

Value¶

An object of class htest containing the chi-squared statistic, p-value and degrees of freedom.

Details¶

The test is the median-centred version of the test of Fligner and Killeen (1976), the "F-K:med $X^2$" of Conover, Johnson and Johnson (1981), whose simulation study found it among the tests of homogeneity of variances most robust to departures from normality; it is the test fligner.test() computes. The residuals of each group are centred at the median of the group, the absolute values of the centred residuals are ranked over all $n$ observations, and the rank $r_i$ is given the score $a_i = \Phi^{-1}\{(1 + r_i / (n + 1)) / 2\}$. With $\bar a_j$ the mean score of group $j$, of size $n_j$, $\bar a$ the mean of all the scores and $s_a^2 = \sum_i (a_i - \bar a)^2 / (n - 1)$, the statistic

\[ X^2 = \sum_j n_j (\bar a_j - \bar a)^2 / s_a^2 \]

is referred to a chi-squared distribution with $k-1$ degrees of freedom, for $k$ groups. No logarithm is taken. The observations whose group is missing are left out before the statistic is computed.

Tied values are given the mean of their ranks. Since 0.16.0 the deviations are those of the residuals $r_i = \sqrt{w_i}\,(y_i - x_i^\top \hat\beta - o_i)$ computed observation by observation, as the Gaussian glm computes them, and not the residuals an lm stores, which its QR decomposition rounds differently from one observation to the next. For an mgcv::gam() fit $x_i$ is the row of its linear predictor matrix. Observations with the same values of the model's variables, weight and offset, such as repeated values of an integer response or of a response recorded with few digits, are then given the same residual, as in exact arithmetic, where their residuals agree to within rounding, even when a basis such as poly() gives them design rows that differ in their last bits; a term of the row position, such as seq_along(x), keeps them apart. This needs the model's variables: it is done when data is given, or when the model frame holds them, as it holds a variable the model uses under its own name (x, not poly(x, 3)). And absolute deviations closer than $8 \epsilon s_j$, with $\epsilon$ = .Machine$double.eps and $s_j$ the largest $|r_i|$ of group $j$ (the smaller of the two for deviations of different groups), are ranked as tied: in a group of even size the two residuals in the middle lie at the same distance from the median of the group, half their difference, a tie that the centring breaks by up to $\epsilon s_j$. On data with repeated observations the statistic can differ much from before: on lm(extra ~ group, data = sleep) it was 0.154 and is 0.213, the value of exact ties. On continuous data it moves by much less (see NEWS). Where neither changes a rank, the value is that of fligner.test() on the residuals of the fit, as before. Two limits remain. A tie between different observations that exact decimal arithmetic would give, as when two pairs of observations of a regressor recorded with few digits differ by the same amounts, is restored only when rounding leaves it within the tolerance, and two forms of one model can then differ by several percent: on 70 observations that repeat 30 values of $x$ recorded with one decimal, y ~ poly(x, 3) + g gives 1.21547, and y ~ x + I(x^2) + I(x^3) + g and y ~ poly(x, 3, raw = TRUE) + g give 1.30280, since poly() breaks in its last bits the ties between repeated values of $x$ in different groups. And in a group whose residuals are far from zero, $10^8$ times their spread, the tolerance may also tie two deviations that differ by a few units in the last place of the residuals.

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¶

Conover, W. J., Johnson, M. E., & Johnson, M. M. (1981). A comparative study of tests for homogeneity of variances, with applications to the outer continental shelf bidding data. Technometrics, 23(4), 351–361. doi:10.1080/00401706.1981.10487680

Fligner, M. A., & Killeen, T. J. (1976). Distribution-free two-sample tests for scale. Journal of the American Statistical Association, 71(353), 210–213. doi:10.1080/01621459.1976.10481517

Examples¶

 data(mtcars)
 mtcars$cyl <- factor(mtcars$cyl)
 m <- lm(mpg ~ wt, data = mtcars)
 performFlignerKilleenTest(m, mtcars, "cyl")