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.
# S3 method for default hlm_resid(object, ...) # S3 method for lmerMod hlm_resid( object, level = 1, standardize = FALSE, include.ls = TRUE, data = NULL, ... ) # S3 method for lme hlm_resid( object, level = 1, standardize = FALSE, include.ls = TRUE, data = NULL, ... )
| object | an object of class |
|---|---|
| ... | do not use |
| level | 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 |
| standardize | for any level, if |
| include.ls | a logical indicating if LS residuals be included in the
return tibble. |
| data | if |
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
Adam Loy loyad01@gmail.com, Jack Moran, Jaylin Lowe
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. #> # … with 170 more rows, and 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. #> # … with 170 more rows, and 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.831fm.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