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()orgibbs()- type
the kind of residual, see the description
- ...
unused
- level
"conditional"forX b + Z u(default),"marginal"for the fixed partX 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 with1 - 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 isNA. Inmodel_mt()the divisor isVar(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 ofnlme, 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.