2

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

BornPerson
  • 123
  • 7

1 Answers1

1

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

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
  • If I specify type="normalized", that would account for the $\phi^2$ factor, but it's a dilation, and therefore doesn't change anything? In the case of autoregressive, there's a autocorrelation factor and a $\phi^2$ factor, and the n-th power of the of autocorrelation coefficient indicating time-decay affects my normalized residuals more so than a compound-symmetry structure? @BenBolker – BornPerson Aug 23 '22 at 14:40
  • Exactly right ... you can use `scaled=TRUE` to scale by the residual standard deviation. – Ben Bolker Aug 23 '22 at 14:42
  • Should this question be expanded to ask about how to compute "normalized" model residuals in geepack/gee? @BenBolker – BornPerson Aug 23 '22 at 15:24
  • If you like, but I have no idea what the answer would be. This might be verging on a [CrossValidated](https://stats.stackexchange.com) question, as I'm not even sure what the definition of normalized model residuals would be for a GEE ... I think the point of a GEE is that the parameter estimates are supposed to be *robust* to a wide range of heteroscedasticity/correlation structures, so you don't need to worry about them ... – Ben Bolker Aug 23 '22 at 18:09