rcorrelationlme4nlmelongitudinal

How does one compute the "normalized" model residuals based via lme4/merMod in R?


The nlme package gives me a way of compiting the normalized residuals using resid(fitted object, type="normalized") but lme4 has no option to do so. I cannot diagnose autocorrelation without this feature in lme4.

I do not think the R stats package resid(lme4object,type="normalized") works, and lme4-object$residuals is not correct syntax.

lmer/glmer are merMod objects.

Class "merMod" of Fitted Mixed-Effect Models Description A mixed-effects model is represented as a merPredD object and a response module of a class that inherits from class lmResp. A model with a lmerResp response has class lmerMod; a glmResp response has class glmerMod; and a nlsResp response has class nlmerMod.

Definition ‘"normalized"’

the normalized residuals (standardized residuals pre-multiplied by the inverse square-root factor of the estimated error correlation matrix) are used.

Aside, what's a "error correlation matrix"? You mean fisher-info/variance-covariance or yule-walker?


How does one compute the normalized model residuals in lme4 or by-hand?


lme4 only gives me the "residual standard deviation", "Pearson", and "deviance residuals", and the ones listed. Documentation Begin:

Description
residuals of merMod objects
Usage
## S3 method for class 'merMod'
residuals(object,
type = if (isGLMM(object)) "deviance" else "response",
scaled = FALSE, ...)
## S3 method for class 'lmResp'
residuals(object,
type = c("working", "response", "deviance", "pearson", "partial"),
...)
## S3 method for class 'glmResp'
residuals(object,
type = c("deviance", "pearson", "working", "response", "partial"),
...)
Arguments
object a fitted [g]lmer (merMod) object
type type of residuals
scaled scale residuals by residual standard deviation (=scale parameter)?
... additional arguments (ignored: for method compatibility)
Details
• The default residual type varies between lmerMod and glmerMod objects: they try to mimic
residuals.lm and residuals.glm respectively. In particular, the default type is "response",
i.e. (observed-fitted) for lmerMod objects vs. "deviance" for glmerMod objects. type="partial"
is not yet implemented for either type.
• Note that the meaning of "pearson" residuals differs between residuals.lm and residuals.lme.
The former returns values scaled by the square root of user-specified weights (if any), but
not by the residual standard deviation, while the latter returns values scaled by the estimated
standard deviation (which will include the effects of any variance structure specified in the
weights argument). To replicate lme behaviour, use type="pearson", scaled=TRUE.

https://cran.r-project.org/web/packages/lme4/lme4.pdf

https://stats.stackexchange.com/questions/80823/do-autocorrelated-residual-patterns-remain-even-in-models-with-appropriate-corre

https://cran.r-project.org/web/packages/nlme/nlme.pdf

https://search.r-project.org/CRAN/refmans/lme4/html/merMod-class.html


Solution

  • nlme::lme allows for modeling correlation and heteroscedasticity in the residual error term, while lme4::lmer doesn't (these are called "R-side" structures in SAS terminology). Specifying type = "normalized" provides residuals that account for/correct for any modeled structure in the residuals; since lme4::lmer doesn't have those structures, normalizing the residuals wouldn't do anything.

    If a linear mixed model (with R-side structures) is expressed as

    y ~ MVN(mu, Sigma_r(phi))
    mu = X beta + Z b
    b ~ MVN(0, Sigma_g(theta))
    

    where X is the fixed-effect model matrix, beta is the FE parameter vector, Z is the random-effect model matrix, b is the vector of BLUPs/conditional modes, Sigma_r and Sigma_g are the covariance matrices of the residual error and the random effects, and phi and theta are the parameter vectors defining these matrices ...

    ... then Sigma_r is the error correlation (covariance) matrix. In lmer, Sigma_r is always a homogeneous diagonal matrix (phi^2*I).