Skip to contents

Residual diagnostics of a mixed-model fit

Source: R/residuals.R
residual_diagnostics.Rd

The standard checks of the residuals of a mixed (animal) model: homoscedasticity by class (by default every fixed class term of the formula, the contemporary group), against the fitted fixed part and against covariates, normality, and within-subject serial correlation for repeated records, each with an effect size and a p-value from a parametric bootstrap under the fitted model (from pivotal discrepancies for a gibbs() chain). Next to them, without a bootstrap p: the largest standardized residuals (not for a chain), with a Bonferroni p, the leverage and a generalized Cook distance, and a summary of the standardized breeding values (mean z^2, skewness and kurtosis, no p-value).

Usage

residual_diagnostics(
  fit,
  data,
  pedigree = NULL,
  genotypes = NULL,
  by = NULL,
  covariates = NULL,
  subject = NULL,
  time = NULL,
  nsim = 199L,
  seed = NULL,
  n_outliers = 10L,
  min_reliability = 0.01
)

Arguments

fit

a fit of model(), a single-trait model_ar1(), or a gibbs() chain run with residual_draws > 0

data

the data.frame of the fit, with its rows in the same order: the bootstrap rebuilds the design from it and stops if the data residuals are not the fit's; a chain's rebuilt design must give its stored fitted values

pedigree, genotypes

as in the fit. The bootstrap checks them with the data (see THE DATA GUARD); for a gibbs() chain they are not checked and enter only the standardized breeding values

by

column names of data whose classes are tested for a common residual variance. By default every fixed class term of the formula

covariates

numeric columns of data to test the residual variance against. By default the cov() terms of the formula

subject, time

columns of data for the within-subject lag-1 autocorrelation of a repeated-records fit; a model_ar1() fit uses its own

nsim

number of bootstrap replicates (or reference draws for a chain). The bootstrap solves nsim normal replicates and, when the error shows excess kurtosis, a pilot of 49 and nsim more with the matched tail (see COST)

seed

seed of R's generator for the replicates; NULL uses its current state

n_outliers

how many of the largest standardized residuals to list

min_reliability

a level whose reliability is below this has no standardized breeding value

Value

an object of class breeding_residual_diagnostics: tests (a data.frame with check, against, statistic, p_value, p_method, effect, primary and p_adjusted, plus draws_below_0.05 for a chain), classes (the phi table of each factor), normality, acf (lags 1 and 2, with the replicates' mean and SD and which lag is tested), outliers, random (by term: summary, largest and z), error_tail (bootstrap only: excess_kurtosis, se, p_value and df; with normal replicates se is NA and df Inf, and excess_kurtosis is NA unless the pre-test rejected and the estimate was too small to match), notes, and the counts.

Details

WHY A BOOTSTRAP. The residuals of an animal model are not independent: with the contemporary group as a fixed effect the residuals of a group sum to zero, and their variance depends on the leverage. The classical tests assume independent observations, and their asymptotic p-values break in both directions. On simulated homoscedastic animal models (validation/residual_diagnostics_calibration.R assintotico, 1000 data sets per cell, rejection at 5%), Fligner-Killeen and Brown-Forsythe on the standardized residual rejected 1.000 with contemporary groups of two records and 0.000 with groups of three to six, and a naive test of the lag-1 autocorrelation of an iid fit with a permanent environment effect (r1 against N(0, 1 / pairs)) rejected 1.000 (lag), because the residual loses the subject's mean. Here every statistic is computed by the same code on the data and on nsim replicates y* = X b + Z u* + e*, with u* and e* drawn from the fitted model (components held at their estimates) and solved with the same factorization of the mixed model equations (the simulation approach to the residuals of a linear mixed model of Schutzenmeister and Piepho 2012, turned into p-values). With normal errors that Monte Carlo p-value is valid under the fitted model for any layout of the groups, conditional on the estimated components. With heavy-tailed errors it rests on the tail matching below and is approximate: it was measured in five layouts only, and failed in one (see MEASURED TYPE I). The draws use R's generator.

THE DATA GUARD. The bootstrap rebuilds the design from data, pedigree and genotypes, factorizes the mixed model equations once at the estimated components and solves the data with that factor. Unless all three are the ones of the fit, the residuals differ from the fit's and it stops: a shuffled column, another pedigree with the same animals, or other genotypes are refused before any replicate. (With a component at zero the matrix of that term does not enter the solution, so it cannot be told apart, and it does not change the replicates either.) A gibbs() chain has no factorization to solve the data with, and only the data and the level names are checked (see the chain paragraph below).

THE TAIL OF THE ERROR. With the contemporary group fixed, one large error spreads over its whole group (the residuals of a group sum to zero), so heavy tails alone make the groups look unequal: with t5 errors and groups of two to four records, the dispersion check rejected 0.099, and the Koenker check against the fitted fixed part 0.187, of homoscedastic data sets against normal replicates (cauda, 1000 data sets). So the error tail is tested first: the within-class kurtosis of the standardized residuals is compared with the normal replicates, and when it exceeds them at 5% the replicates of the dispersion, Koenker and lag checks are drawn again with a t error matched to the excess kurtosis of the error. That excess kurtosis is estimated from the same statistic, which is linear in the fourth cumulant of the error, with the slope measured on a pilot of 49 replicates with Laplace errors (excess kurtosis 3); the t has 4 + 6 / kurtosis df. $error_tail holds the estimate, its bootstrap standard error (the spread of the statistic among the matched replicates), the p-value of the excess and the df (Inf with normal replicates). The estimate is noisy and biased down for a tail as heavy as t5 (the sample kurtosis of such a tail is small most of the time): with t5 errors, whose excess kurtosis is 6, it averaged 2.4 to 4.7 by layout where the tail was matched, and the pre-test detected the tail in only 0.65 (groups of two to four records) to 1.00 (sire model) of the data sets; the rest kept normal replicates (tipo1). With normal errors the tail was matched in 0.041 to 0.056 of the data sets, which only makes the dispersion check conservative there. When the pre-test rejects but the pilot gives no positive slope, or the estimate is 0.05 or less, the replicates stay normal and $notes and print() say so; it did not happen in any of the 10,000 data sets of tipo1. The normality check is always against normal replicates. In classes of fewer than ten records the statistic is the kurtosis of all the standardized residuals, which also absorbs part of a heteroscedasticity and makes the dispersion check conservative there.

MEASURED TYPE I (validation/residual_diagnostics_calibration.R tipo1, 1000 simulated data sets per cell, five stored layouts: contemporary groups of 2 to 40, of 2 to 4, repeated records with a permanent environment, a sire model with groups of 40, record weights; normal and t5 errors; rejection at p <= 0.05, Monte Carlo SE 0.007): dispersion 0.045 to 0.064; koenker 0.040 to 0.067, except against fitted_fixed with t5 errors and groups of two to four records, 0.094 (outside the gate; see below); normality 0.051 to 0.059 with normal errors; lag1 0.039 (lag, iid fit with a permanent environment, eight records per subject, 1000 data sets) and lag2 of the innovations 0.060 (500 data sets, model_ar1() with rho = 0.3). Only these layouts and sizes (600 to 2400 records) were measured.

THE CHECKS (column check of $tests):

  • dispersion (primary, one per factor of by): with c_i = 1 - h_ii (the effective degrees of freedom of the record), the Bartlett-type ratio -sum_k m_k log(a_k / a), m_k = sum c over class k, a_k = sum(c s) / m_k and a the pooled mean, computed on the squared Fligner-Killeen scores s_i = Phi^-1(1/2 + R_i / (2 (n + 1)))^2 of the ranks R_i of |r_i| (Conover, Johnson and Johnson 1981) instead of on r_i^2, so a large error counts by its rank and not by its square (the statistic is the same for any monotone transform of |r|). The tail still moves it through the way one error spreads over its group, which is what the tail matching above is for. $classes lists, by level, phi_k = sum(c r^2) / sum(c) relative to the pooled value: the ratio of the residual variance of the class to the common one.

  • brown_forsythe (secondary): Levene's test with the median (Brown and Forsythe 1974), on the classes with at least three records.

  • lr_phi (secondary): the same ratio on r_i^2. Normal theory: with t5 errors it rejected 0.080 to 0.287 of homoscedastic data sets with the tail matched (tipo1), and 0.37 to 0.47 against normal replicates (cauda); use it only with near-normal residuals.

  • koenker (primary, against fitted_fixed always and against each covariate): n R^2 of the scores s_i of dispersion regressed on x and x^2, the studentized test of Koenker (1981) on rank scores (a large error, spread over its contemporary group, meets the fitted fixed part of that group at a high leverage, and the bounded scores keep it from dominating r^2). With groups of two to four records and t5 errors the check against fitted_fixed stays liberal (0.094 at 5%): the fitted fixed part of a small group carries the mean of its own errors, and under heavy tails that ties it to the size of the residuals in a way the matched tail reproduces only in part, and not at all in the data sets where the pre-test missed the tail. $notes says so when the tail is matched and most records are in classes of fewer than five. Never against the conditional fitted value X b + Z u: the residual is confounded with Z u.

  • normality (primary): Shapiro-Wilk W (lower is worse) up to 5000 records, Anderson-Darling A2 above, with the bootstrap p. The asymptotic Shapiro-Wilk p of the standardized residual rejected 0.040 to 0.061 in the stored layouts with normal errors (tipo1), close to nominal there, but nothing guarantees it for dependent residuals, and the bootstrap p costs nothing more.

  • lag1 (primary, with subject = and time =): the lag-1 autocorrelation of the standardized residual between consecutive records of a subject, two-sided around its bootstrap mean, which is NOT zero: removing the subject's mean pushes it to about -0.105 with eight records per subject. A residual from an iid fit does not estimate the AR coefficient; model_ar1() does, and the likelihood-ratio test of model_ar1() against model() is the test of the AR(1) itself.

  • lag2 (primary, model_ar1() only): the lag-2 autocorrelation of the innovations (see residuals.breeding_fit()), the check of a serial correlation the AR(1) does not capture. The lag 1 of the innovations is reported in $acf but not tested: rho is estimated from it, and replicates with rho held at its estimate do not repeat that fit (the residual autocorrelations of a fitted autoregression lose degrees of freedom at the fitted lags, Box and Pierce 1970), so its bootstrap p is conservative: across data sets that autocorrelation had a SD of 0.0069, against 0.0186 among the replicates, and with a normal approximation to the bootstrap distribution none of 500 data sets would have rejected at 5% (largest |z| 1.32, lag).

READING IT. A residual variance that differs among groups matters for selection, which then favors the more variable groups (Hill 1984). With many records a small p answers whether there is ANY departure, and a trivial one is enough (Li et al. 2024); read the effect column: for dispersion the standard deviation across levels of the log of the residual standard deviation, in excess of what the replicates give under homoscedasticity, on levels with at least five effective degrees of freedom (0.17 is half the groups with twice the residual variance; 0.35 a residual SD varying by a factor of about 1.4 across groups); for koenker the relative change of the residual SD per SD of the covariate; for normality the excess kurtosis (and $normality has the skewness and the ratio of tails beyond 3 SD to the normal's 0.27%); for the lags the autocorrelation minus its bootstrap mean. The dispersion effect is truncated at zero, so it is positive under homoscedasticity too: it averaged 0.019 to 0.036 with normal errors and 0.061 to 0.090 with t5 errors in the type I runs (tipo1), and values of that size are noise. In an animal model with one record per animal it is also attenuated, because the breeding value of an animal absorbs part of the extra residual variance of its group: with h2 = 0.4 it averaged 0.112 for a true 0.173 and 0.222 for 0.35, while a sire model with forty records per group gave 0.168 and 0.331 (poder, 500 data sets each). With classes of fewer than five effective degrees of freedom there is no effect (NA). A non-rejection with small groups does not prove homoscedasticity: with twice the residual variance in half of contemporary groups of 2 to 40 records (640 records) the dispersion check rejected 0.34 of the time, and with groups of two to four 0.046, no more than under the model (0.12 with a residual SD varying by a factor of about 1.4). With groups that small the kurtosis of all the residuals absorbs the heteroscedasticity and the tail gets matched; brown_forsythe had 0.10 and 0.35 there (poder, 500 data sets). Several checks are run at once: p_adjusted is Holm's adjustment over the primary rows.

The p-value is conditional on the estimated components: a component at a boundary is held there.

Records with 1 - h_ii < 1e-8 (a contemporary group of one record: the residual is zero by construction) and classes with a single record are left out and counted in $notes.

OUTLIERS ($outliers): the n_outliers largest |r|, with a Bonferroni p min(1, 2 n Phi(-|r|)), the leverage and a generalized Cook distance r^2 h / ((1 - h) sum(h)): the change of all the solutions, fixed AND random, when the record is deleted with the components held, in the metric of the mixed model equations (Hodges 1998), which in an animal model is dominated by the record's own breeding value. It is not the Cook distance of the fixed effects alone (Christensen, Pearson and Johnson 1992).

RANDOM EFFECTS ($random): the standardized prediction z_j = u_j / sqrt(s2 K_jj - PEV_j), unit variance under the model (Var(u^) = G - C^uu, Henderson 1975), with K_jj the diagonal of the covariance matrix the mixed model equations of the fit used, read from the selective inverse of its inverse when the design is rebuilt: 1 + F for a pedigree term, H_ii for every animal of a single step (genotyped or not; the prior variance accuracy() uses for a non-genotyped animal is 1 + F, which is not its variance under H), the diagonal of K for a kernel() term, and 1 for an iid term. A level with a reliability 1 - PEV / (s2 K_jj) under min_reliability (an animal without information) has no standardization and is left out and counted. Scalar terms only; a reaction norm is skipped.

A gibbs() chain is diagnosed by pivotal discrepancy measures (Johnson 2007; Yuan and Johnson 2012): each stored draw of the realized residual e_s / sqrt(s2e_s) is iid N(0, 1) under the model (a posterior draw of the realized errors, not an estimate), so the same statistics have an exact reference distribution, simulated here from nsim standard normal vectors, with no leverage and no confounding. One p-value per draw; the reported one is min(1, 2 median), since the median of p-values that are each uniform is at most alpha with probability at most 2 alpha (Markov), so it is a conservative test. The argument needs proper priors and is approximate under the flat prior of b. The posterior predictive p-value of Gelman, Meng and Stern (1996) is not used: it is not uniform under the model and concentrates near 0.5. Measured (validation/residual_diagnostics_calibration.R pivotal, 200 simulated data sets per cell, 20 draws, nsim = 199): under the model (contemporary groups of 2 to 40, 640 records) the mean fraction of draws with p < 0.05 (column draws_below_0.05) was 0.044 to 0.055 for every check, and the combined p rejected 0.000 to 0.010 at 5%; with a residual SD varying by a factor of about 1.4 across groups (2000 records) the combined p of dispersion rejected 1.000, of normality 0.72 (the mixture of variances is not normal) and of koenker against the covariate 0.005. For family = "probit" the draws are of the latent residual (variance 1): under the model, with the contemporary group random and by = "cg" (pivotal_probit, 640 records; with the group fixed, a binary at the median leaves 7 to 27 of the 56 groups, median 17 over 200 draws of the data, with every record in one category, which gibbs() refuses), the mean fraction of draws with p < 0.05 was 0.042 to 0.047 and the combined p rejected none of 200 data sets. The chain is checked against the data too: the design rebuilt from data times the posterior means must give the fitted values the chain stored, and each random group must have the levels of the chain. The check runs on the engine's standardized scale, so a covariate with a large offset (a timestamp) passes, and a whole covariate column changed by a x + c with a > 0 is not seen: those data give the same standardized design. Beyond that the pedigree and the genotypes are NOT checked against the chain (it keeps nothing of its K to compare with): they enter only the K_jj of the standardized breeding values, so another pedigree or other genotypes with the same animals change those z without an error, and only those. The standardized breeding values of a chain also use the posterior means of the components and the posterior variances in place of the PEV, and are approximate and not calibrated.

COST. The design is built, and the mixed model equations and each K^-1 factorized, once per call; each replicate is then one draw and one solve with that factor. In a single step the build (G, its inverse, H^-1 and the diagonal of H) is most of the cost: with 2000 genotyped animals among 3600 records it took 6.9 to 9.3 s, and the whole check with nsim = 39 7.6 to 8.7 s, against 9.8 s for the fit at given components (tempo_su, one thread, three runs on the same data). The replicates come in batches of at most 1.25e7 / n (two matrices of about 100 MB each), whose statistics are kept and the replicates dropped, so memory does not grow with nsim. When the tail is matched the work is 2 nsim + 49 replicates instead of nsim. In an animal model of 24,000 records in 777 contemporary groups, one thread, the same fit and data three times, with nsim = 199 and normal replicates (the tail was not matched): 84, 94 and 100 s (tempo; the model() fit itself took 734 s), on a machine shared with other jobs.

Not available: model_mt() (the residuals are; the bootstrap of the multi-trait fit is not), a multi-trait model_ar1(), and the threshold and survival fits, whose residuals are plug-in and not calibrated (see residuals.breeding_fit()). The validation of the genetic trend of an evaluation (Boichard et al. 1995, the basis of the Interbull trend tests) and of its Mendelian sampling variance across years (Tyriseva et al. 2018) check an evaluation over time, not the residuals of one fit, and are not here. Neither is the score-based test of heterogeneity in linear mixed models of Wang et al. (2021), which was not compared with these checks.

References

Koenker, R. (1981). A note on studentizing a test for heteroscedasticity. Journal of Econometrics 17:107-112.

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:351-361.

Box, G.E.P. & Pierce, D.A. (1970). Distribution of residual autocorrelations in autoregressive-integrated moving average time series models. Journal of the American Statistical Association 65:1509-1526.

Henderson, C.R. (1975). Best linear unbiased estimation and prediction under a selection model. Biometrics 31:423-447.

Brown, M.B. & Forsythe, A.B. (1974). Robust tests for the equality of variances. Journal of the American Statistical Association 69:364-367.

Wang, T., Merkle, E.C., Anguera, J.A. & Turner, B.M. (2021). Score-based tests for detecting heterogeneity in linear mixed models. Behavior Research Methods 53:216-231.

Schutzenmeister, A. & Piepho, H.-P. (2012). Residual analysis of linear mixed models using a simulation approach. Computational Statistics and Data Analysis 56:1405-1416.

Hill, W.G. (1984). On selection among groups with heterogeneous variance. Animal Production 39:473-477.

Johnson, V.E. (2007). Bayesian model assessment using pivotal quantities. Bayesian Analysis 2:719-734.

Yuan, Y. & Johnson, V.E. (2012). Goodness-of-fit diagnostics for Bayesian hierarchical models. Biometrics 68:156-164.

Boichard, D., Bonaiti, B., Barbat, A. & Mattalia, S. (1995). Three methods to validate the estimation of genetic trend for dairy cattle. Journal of Dairy Science 78:431-437.

Tyriseva, A.-M., Fikse, W.F., Mantysaari, E.A., Jakobsen, J., Aamand, G.P., Durr, J. & Lidauer, M.H. (2018). Validation of consistency of Mendelian sampling variance. Journal of Dairy Science 101:2187-2198.

Hodges, J.S. (1998). Some algebra and geometry for hierarchical models, applied to diagnostics. Journal of the Royal Statistical Society B 60:497-536.

Christensen, R., Pearson, L.M. & Johnson, W. (1992). Case-deletion diagnostics for mixed models. Technometrics 34:38-45.

Gelman, A., Meng, X.-L. & Stern, H. (1996). Posterior predictive assessment of model fitness via realized discrepancies. Statistica Sinica 6:733-807.

Li, W., Cook, D., Tanaka, E. & VanderPlas, S. (2024). A plot is worth a thousand tests: assessing residual diagnostics with the lineup protocol. Journal of Computational and Graphical Statistics.