Skip to contents

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_mt() computes its 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) 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_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"), ...)

Arguments

object

a fit of model(), model_mt() or gibbs()

type

the kind of residual, see the description

...

unused

level

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

Value

a vector with one value per row of the data of the fit (NA for a record that did not enter), and a matrix with one column per trait for model_mt().

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_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.

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.

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.