Calculating residuals from HLMs

Description

HLMresid is a function that extracts residuals from a hierarchical linear model fit using lmer. That is, it is a unified framework that extracts/calculates residuals from mer or lmerMod objects.

Usage

## Default S3 method:
HLMresid(object, ...)

## S3 method for class 'mer'
HLMresid(object, level, type = "EB", sim = NULL,
  standardize = FALSE, ...)

## S3 method for class 'lmerMod'
HLMresid(object, level, type = "EB", sim = NULL,
  standardize = FALSE, ...)

Arguments

Details

This function extracts residuals from the model, and can extract residuals estimated using least squares (LS) or Empirical Bayes (EB). This unified framework enables the analyst to more easily conduct an upward residual analysis during model exploration/checking.

The HLMresid function provides a wrapper that will extract residuals from a fitted mer or lmerMod object. The function provides access to residual quantities already made available by the functions resid and ranef, but adds additional functionality. Below is a list of types of residuals that can be extracted.

raw level-1 residuals

These are equivalent to the residuals extracted by resid if level = 1, type = "EB", and standardize = FALSE is specified. You can also specify type = "LS" for LS residuals that are not equivalent to those from resid.

standardized level-1 residuals

Specify level = 1, and standardize = TRUE. This works with both type = "EB" or "LS".

semi-standardized level-1 residuals

Specify level = 1, type = "LS" and standardize = "semi".

raw group level residuals

These are equivalent to extracting the predicted random effects for a given group using ranef. Set level to a grouping factor name and type = "EB". type = "LS" can also be specified, though this is less common.

standardized group level residuals

Set level to a grouping factor name, type = "EB", and standardized = TRUE. This will not produce standardized residuals for type = "LS".

marginal residuals

The marginal residuals can be obtained by setting level = "marginal". Only type = "EB" is implemented.

cholesky residuals

These are essentially standardized marginal residuals. To obtain cholesky residuals set level = "marginal", type = "EB", and standardize = TRUE.

Note that standardize = "semi" is only implemented for level-1 LS residuals.

Author(s)

Adam Loy loyad01@gmail.com

References

Hilden-Minton, J. (1995) Multilevel diagnostics for mixed and hierarchical linear models. University of California Los Angeles.

Houseman, E. A., Ryan, L. M., & Coull, B. A. (2004) Cholesky Residuals for Assessing Normal Errors in a Linear Model With Correlated Outcomes. Journal of the American Statistical Association, 99(466), 383–394.

See Also

LSresids, resid, ranef

Examples

library(lme4)
fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)

# level-1 residuals
all.equal(HLMresid(object = fm1, level = 1, type = "EB"), resid(fm1)) ## EB
r1LS <- HLMresid(object = fm1, level = 1, type = "LS") ## raw LS resids
head(r1LS)
r1LS.std <- HLMresid(object = fm1, level = 1, type = "LS", standardize = TRUE) ## std. LS resids
head(r1LS.std)

# level-2 residuals
all.equal(r2EB <- HLMresid(object = fm1, level = "Subject", type = "EB"), ranef(fm1)[["Subject"]])
r2EB.std <- HLMresid(object = fm1, level = "Subject", type = "EB", standardize = TRUE)
head(r2EB)
head(r2EB.std)

# marginal residuals
mr <- HLMresid(object = fm1, level = "marginal")
cholr <- HLMresid(object = fm1, level = "marginal", standardize = TRUE) # Cholesky residuals

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(HLMdiag)

Attaching package: 'HLMdiag'

The following object is masked from 'package:stats':

    covratio

> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/HLMdiag/HLMresid.mer.Rd_%03d_medium.png", width=480, height=480)
> ### Name: HLMresid.default
> ### Title: Calculating residuals from HLMs
> ### Aliases: HLMresid HLMresid.default HLMresid.lmerMod HLMresid.mer
> ### Keywords: models regression
> 
> ### ** Examples
> 
> library(lme4)
Loading required package: Matrix
> fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
> 
> # level-1 residuals
> all.equal(HLMresid(object = fm1, level = 1, type = "EB"), resid(fm1)) ## EB
[1] TRUE
> r1LS <- HLMresid(object = fm1, level = 1, type = "LS") ## raw LS resids
> head(r1LS)
  Reaction Days   LS.resid   fitted
1 249.5600    0   5.367331 244.1927
2 258.7047    1  -7.252672 265.9574
3 250.8006    2 -36.921474 287.7221
4 321.4398    3  11.953024 309.4868
5 356.8519    4  25.600421 331.2515
6 414.6901    5  61.673919 353.0162
> r1LS.std <- HLMresid(object = fm1, level = 1, type = "LS", standardize = TRUE) ## std. LS resids
> head(r1LS.std)
  Reaction Days   LS.resid   fitted  std.resid
1 249.5600    0   5.367331 244.1927  0.1388498
2 258.7047    1  -7.252672 265.9574 -0.1750999
3 250.8006    2 -36.921474 287.7221 -0.8511542
4 321.4398    3  11.953024 309.4868  0.2677905
5 356.8519    4  25.600421 331.2515  0.5657375
6 414.6901    5  61.673919 353.0162  1.3629169
> 
> # level-2 residuals
> all.equal(r2EB <- HLMresid(object = fm1, level = "Subject", type = "EB"), ranef(fm1)[["Subject"]])
[1] TRUE
> r2EB.std <- HLMresid(object = fm1, level = "Subject", type = "EB", standardize = TRUE)
Warning message:
In ranef.merMod(object, postVar = TRUE) :
  'postVar' is deprecated: please use 'condVar' instead
> head(r2EB)
    (Intercept)       Days
308    2.258565  9.1989719
309  -40.398577 -8.6197032
310  -38.960246 -5.4488799
330   23.690498 -4.8143313
331   22.260203 -3.0698946
332    9.039526 -0.2721707
> head(r2EB.std)
    (Intercept)       Days
308   0.1871093  3.9911595
309  -3.3467918 -3.7398320
310  -3.2276343 -2.3641064
330   1.9626228 -2.0887947
331   1.8441309 -1.3319357
332   0.7488732 -0.1180867
> 
> # marginal residuals
> mr <- HLMresid(object = fm1, level = "marginal")
> cholr <- HLMresid(object = fm1, level = "marginal", standardize = TRUE) # Cholesky residuals
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>