Residuals and fitted values of a mixed-model fit

The conditional residual e = y - X b - Z u is computed by the engine at the final components, together with the fitted values and, in model(), the leverage h_ii of each record in the mixed model equations (Hodges 1998): h_ii = w_i' C^-1 w_i with w_i the row of [X Z] scaled by sqrt(weight) and C the coefficient matrix, read from the selective inverse the PEV already uses, so it costs no extra factorization. model_ar1() and model_mt() compute their residuals and resid_var with one more assembly, factorization and selective inverse at the final components, at the end of every fit: the cost of about one more REML iteration, and a second pass through the peak memory of the fit (if that pass fails, the fit is kept and the residuals are NA, with the reason in message). The variance of the residual is NOT the residual variance: Var(e_i) = s2e (1 - h_ii) / w_i, which is smallest in small contemporary groups. On a simulated animal model with one record per animal (h2 = 0.4, true components, validation/residual_diagnostics_calibration.R alavanca) it averaged about 0.35 s2e in groups of two records and 0.69 s2e in groups of forty; the figures move with h2 and the layout. The raw conditional residual is heteroscedastic by construction, and only the standardized one has unit variance under the model.

Usage

# S3 method for class 'breeding_fit'
residuals(
  object,
  type = c("conditional", "standardized", "marginal", "pearson"),
  ...
)

# S3 method for class 'breeding_fit'
fitted(object, level = c("conditional", "marginal"), ...)

# S3 method for class 'breeding_fit_ar1'
residuals(
  object,
  type = c("conditional", "standardized", "normalized", "marginal", "pearson"),
  ...
)

# S3 method for class 'breeding_fit_ar1'
fitted(object, level = c("conditional", "marginal"), ...)

# S3 method for class 'breeding_fit_mt'
residuals(
  object,
  type = c("conditional", "standardized", "marginal", "pearson"),
  ...
)

# S3 method for class 'breeding_fit_mt'
fitted(object, level = c("conditional", "marginal"), ...)

# S3 method for class 'breeding_gibbs'
residuals(object, type = c("conditional", "marginal"), ...)

# S3 method for class 'breeding_gibbs'
fitted(object, level = c("conditional", "marginal"), ...)

# S3 method for class 'breeding_fit_thr'
residuals(
  object,
  data,
  type = c("surrogate", "pearson"),
  missing_code = NULL,
  seed = NULL,
  ...
)

# S3 method for class 'breeding_fit_surv'
residuals(
  object,
  data,
  censor,
  type = c("coxsnell", "martingale", "deviance", "quantile"),
  entry = NULL,
  subject = NULL,
  seed = NULL,
  ...
)

Arguments

object

a fit of model(), model_ar1(), model_mt(), gibbs(), model_threshold() or model_survival()

type

the kind of residual, see the description

...

unused

level

"conditional" for X b + Z u (default), "marginal" for the fixed part X b

data

threshold and survival only: the data of the fit

missing_code

threshold only: the missing code of the fit, if one was given

seed

threshold and survival only: seed of the random draws of the surrogate and quantile residuals (NULL leaves R's generator as it is)

censor, entry, subject

survival only: as given to model_survival()

Value

a vector with one value per row of the data of the fit (NA for a record that did not enter), a matrix with one column per trait for model_mt() and a multi-trait model_ar1(), a matrix with one column per category for the threshold "pearson", and for survival with subject = a vector named by subject.

Details

Types:

• "conditional": y - X b - Z u, in the units of the data.

• "standardized": e_i sqrt(w_i) / sqrt(s2e (1 - h_ii)), unit variance under the model given the components. A record with 1 - h_ii < 1e-8 (a contemporary group or any fixed level with a single record) has a residual of zero by construction and no standardization: it is NA. In model_ar1() and model_mt() the divisor is Var(e) = R0[t, t] - w' C^-1 w, from the absolute-unit equations.

• "marginal": y - X b, which keeps the random effects in.

• "pearson": e_i sqrt(w_i) / sqrt(s2e), the convention of nlme, NOT standardized by the leverage.

• "normalized" (model_ar1() only): the innovations of the AR(1)/CAR(1) process, (e_t - g e_(t-1)) / sqrt(s2e (1 - g^2)) with g = rho^dt along the subject's records in time order, dt the gap in steps of time_step on a grid (g negative for rho < 0 and an odd dt) and in the unit of the time on a continuous axis (time_step NA; see time in model_ar1()); in a POSIXct series of daily records dt is 1, not 86400 (the first record of a subject is e_1 / sqrt(s2e)), the "normalized" residual of nlme: the residual premultiplied by the inverse Cholesky factor of its within-subject covariance R. With several traits the vector of a record is also premultiplied by the inverse Cholesky factor of R0, and a record with a missing trait is handled by the same factorization of the covariance of the OBSERVED cells (a Kalman filter along the subject's records): the next record is conditioned on what was observed, and only the missing cells are NA. As a conditional residual it is computed at the estimates and carries the same confounding as "conditional".

The conditional residual is confounded with the random effects (it carries R P Z u; Hilden-Minton 1995, Loy and Hofmann 2015), so a plot of it against the conditional fitted value X b + Z u shows a trend even under the model: compare against the fixed part fitted(fit, level = "marginal") instead.

gibbs() returns the posterior mean of the realized residual y - W loc_s, which by linearity is y - X E[b] - Z E[u] and shares the same confounding; type = "marginal" adds back Z E[u], and fitted() gives X E[b] + Z E[u] or X E[b]. For family = "probit" the residual is that of the liability (the drawn liability minus W loc_s) and the fitted values are on the liability scale. Only these two types exist for a chain: it has no leverage, and its residual draws are what residual_diagnostics() uses instead.

The threshold and survival methods need the data again and return plug-in residuals, computed with the shrunken random effects as if they were known: they are as confounded as the Gaussian conditional residual and their distribution under the model is not calibrated (see Details). model_threshold() (ordinal mode): "surrogate", the surrogate residual of Liu and Zhang (2018), a draw of the liability from the normal with mean x'b + z'u truncated to the interval of the observed category, minus that mean (N(0, 1) under the model with the true parameters); "pearson", the matrix of (1{y = c} - p_c) / sqrt(p_c (1 - p_c)) by category. model_survival(): "coxsnell", the cumulative hazard of each subject summed over its pieces; "martingale", event - coxsnell (Therneau et al. 1990); "deviance", sign(M) sqrt(-2 (M + event log(event - M))); "quantile", the normalized randomized survival probability of Li, Wu and Feng (2021), qnorm(u) with u = S(t) at an event and u ~ U(0, S(c)) at a censoring.

Measured on simulated data (plug-in, with the PEV of the latent effect as the only departure from the true parameters; 1000 data sets per cell, validation/residual_diagnostics_calibration.R plugin): the Kolmogorov-Smirnov test of the quantile residual against N(0, 1) rejected 0.054, 0.062 and 0.243 of the data sets at a PEV of 0, 0.05 and 0.15 (Weibull, 26% censoring, 3000 subjects), and the same test of the surrogate residual 0.043, 0.067 and 0.175 (probit, 5000 records); a PEV of 0.15 or more is ordinary with one record per animal. The Shapiro-Wilk test of the surrogate residual stayed at 0.043 with a PEV of 0.15, because it is scale-free and blind to that mixture. No test of these residuals is offered: read them as descriptive.

References

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

Loy, A. & Hofmann, H. (2015). Are you normal? The problem of confounded residual structures in hierarchical linear models. Journal of Computational and Graphical Statistics 24:1191-1209.

Pinheiro, J.C. & Bates, D.M. (2000). Mixed-Effects Models in S and S-PLUS. Springer.

Liu, D. & Zhang, H. (2018). Residuals and diagnostics for ordinal regression models: a surrogate approach. Journal of the American Statistical Association 113:845-854.

Therneau, T.M., Grambsch, P.M. & Fleming, T.R. (1990). Martingale-based residuals for survival models. Biometrika 77:147-160.

Li, L., Wu, T. & Feng, C. (2021). Model diagnostics for censored regression via randomized survival probabilities. Statistics in Medicine 40:1482-1497.