hlm_resid takes a hierarchical linear model fit as a lmerMod or
lme object and extracts residuals and predicted values from the model,
using Least Squares (LS) and Empirical Bayes (EB) methods. It then appends
them to the model data frame in the form of a tibble inspired by the augment
function in broom. This unified framework enables the analyst to more
easily conduct an upward residual analysis during model exploration/checking.
# Default S3 method
hlm_resid(object, ...)
# S3 method for class 'lmerMod'
hlm_resid(
object,
level = 1,
standardize = FALSE,
include.ls = TRUE,
data = NULL,
...
)
# S3 method for class 'lme'
hlm_resid(
object,
level = 1,
standardize = FALSE,
include.ls = TRUE,
data = NULL,
...
)an object of class lmerMod or lme.
do not use
which residuals should be extracted: 1 for within-group
(case-level) residuals, the name of a grouping factor for between-group
residuals (as defined in flist in lmerMod objects or in
groups in lme objects)
for any level, if standardize = TRUE the
standardized residuals will be returned for any group; for level-1 only, if
standardize = "semi" then the semi-standardized level-1 residuals
will be returned
a logical indicating if LS residuals be included in the
return tibble. include.ls = FALSE decreases runtime substantially.
if na.action = na.exclude, the user must provide the data
set used to fit the model, otherwise NULL.
The hlm_resid function provides a wrapper that will extract
residuals and predicted values from a fitted lmerMod or lme
object. The function provides access to residual quantities already made available by
the functions resid, predict, and ranef, but adds
additional functionality. Below is a list of types of residuals and predicted
values that are extracted and appended to the model data.
.resid and .fittedResiduals calculated using
the EB method (using maximum likelihood). Level-1 EB residuals are interrelated
with higher level residuals. Equivalent to the residuals extracted by
resid(object) and predict(object) respectively. When
standardize = TRUE, residuals are standardized by sigma components of
the model object.
.ls.resid and .ls.fittedResiduals calculated calculated by fitting separate LS regression models for each group. Level-1 LS residuals are unconfounded by higher level residuals, but unreliable for small within-group sample sizes.
.mar.resid and .mar.fittedMarginal residuals only
consider the fixed effect portion of the estimates. Equivalent to
resid(object, level=0) in nlme, not currently implemented
within the lme4::resid function. When standardize = TRUE,
Cholesky marginal residuals are returned.
.ranef.var_nameThe group level random effects using the EB method of
estimating parameters. Equivalent to ranef(object) on the specified
level. EB residuals are preferred at higher levels as LS residuals are dependent
on a large sample size.
.ls.var_nameThe group level random effects using the LS method of
estimating parameters. Calculated using ranef on a lmList
object to compare the random effects of individual models to the global
model.
Note that standardize = "semi" is only implemented for level-1 LS residuals.
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.
David Robinson and Alex Hayes (2020). broom: Convert Statistical Analysis Objects into Tidy Tibbles. R package version 0.5.6. https://CRAN.R-project.org/package=broom
data(sleepstudy, package = "lme4")
fm.lmer <- lme4::lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
fm.lme <- nlme::lme(Reaction ~ Days, random = ~Days|Subject, sleepstudy)
# level-1 and marginal residuals
fm.lmer.res1 <- hlm_resid(fm.lmer) ## raw level-1 and mar resids
fm.lmer.res1
#> # A tibble: 180 x 10
#> .id Reaction Days Subject .resid .fitted .ls.resid .ls.fitted .mar.resid
#> <dbl> <dbl> <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 250. 0 308 -4.10 254. 5.37 244. -1.85
#> 2 2 259. 1 308 -14.6 273. -7.25 266. -3.17
#> 3 3 251. 2 308 -42.2 293. -36.9 288. -21.5
#> 4 4 321. 3 308 8.78 313. 12.0 309. 38.6
#> 5 5 357. 4 308 24.5 332. 25.6 331. 63.6
#> 6 6 415. 5 308 62.7 352. 61.7 353. 111.
#> 7 7 382. 6 308 10.5 372. 7.42 375. 68.0
#> 8 8 290. 7 308 -101. 391. -106. 397. -34.5
#> 9 9 431. 8 308 19.6 411. 12.3 418. 95.4
#> 10 10 466. 9 308 35.7 431. 26.3 440. 121.
#> # i 170 more rows
#> # i 1 more variable: .mar.fitted <dbl>
fm.lme.std1 <- hlm_resid(fm.lme, standardize = TRUE) ## std level-1 and mar resids
fm.lme.std1
#> # A tibble: 180 x 10
#> .id Reaction Days Subject .std.resid .fitted .std.ls.resid .ls.fitted
#> <dbl> <dbl> <dbl> <fct> <dbl> <dbl> <dbl> <dbl>
#> 1 1 250. 0 308 -0.160 254. 0.139 244.
#> 2 2 259. 1 308 -0.571 273. -0.175 266.
#> 3 3 251. 2 308 -1.65 293. -0.851 288.
#> 4 4 321. 3 308 0.343 313. 0.268 309.
#> 5 5 357. 4 308 0.958 332. 0.566 331.
#> 6 6 415. 5 308 2.45 352. 1.36 353.
#> 7 7 382. 6 308 0.412 372. 0.166 375.
#> 8 8 290. 7 308 -3.95 391. -2.45 397.
#> 9 9 431. 8 308 0.766 411. 0.296 418.
#> 10 10 466. 9 308 1.39 431. 0.680 440.
#> # i 170 more rows
#> # i 2 more variables: .chol.mar.resid <dbl>, .mar.fitted <dbl>
# level-2 residuals
fm.lmer.res2 <- hlm_resid(fm.lmer, level = "Subject") ## level-2 ranefs
fm.lmer.res2
#> # A tibble: 18 x 5
#> Subject .ranef.intercept .ranef.days .ls.intercept .ls.days
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 308 2.26 9.20 -7.21 11.3
#> 2 309 -40.4 -8.62 -46.4 -8.21
#> 3 310 -39.0 -5.45 -47.9 -4.35
#> 4 330 23.7 -4.81 38.3 -7.46
#> 5 331 22.3 -3.07 34.3 -5.20
#> 6 332 9.04 -0.272 12.8 -0.901
#> 7 333 16.8 -0.224 23.6 -1.33
#> 8 334 -7.23 1.07 -11.2 1.79
#> 9 335 -0.334 -10.8 11.6 -13.3
#> 10 337 34.9 8.63 38.7 8.56
#> 11 349 -25.2 1.17 -36.3 3.03
#> 12 350 -13.1 6.61 -25.6 9.04
#> 13 351 4.58 -3.02 9.74 -4.03
#> 14 352 20.9 3.54 25.0 3.10
#> 15 369 3.28 0.872 3.56 0.881
#> 16 370 -25.6 4.82 -41.0 7.59
#> 17 371 0.807 -0.988 2.23 -1.28
#> 18 372 12.3 1.28 15.6 0.831
fm.lme.res2 <- hlm_resid(fm.lme, level = "Subject", include.ls = FALSE) ##level-2 ranef, no LS
fm.lme.res2
#> # A tibble: 18 x 3
#> Subject .ranef.intercept .ranef.days
#> <chr> <dbl> <dbl>
#> 1 308 2.26 9.20
#> 2 309 -40.4 -8.62
#> 3 310 -39.0 -5.45
#> 4 330 23.7 -4.81
#> 5 331 22.3 -3.07
#> 6 332 9.04 -0.272
#> 7 333 16.8 -0.224
#> 8 334 -7.23 1.07
#> 9 335 -0.334 -10.8
#> 10 337 34.9 8.63
#> 11 349 -25.2 1.17
#> 12 350 -13.1 6.61
#> 13 351 4.58 -3.02
#> 14 352 20.9 3.54
#> 15 369 3.28 0.872
#> 16 370 -25.6 4.82
#> 17 371 0.807 -0.988
#> 18 372 12.3 1.28