r2mlm_comp: Compute R-squared differences between two multilevel models, automatically inputting parameter estimates.

Description

r2mlm_comp reads in two multilevel models (MLMs) (generated using lmer or nlme) under comparison (designated Model A and Model B), and outputs all R-squared measures in the Rights and Sterba (2019) framework for both models, as well as R-squared differences between the two models. Definitions of these R-squared difference measures are provided in Rights & Sterba (2020) Table 1; importantly, to detect the impact of a specific kind of term (e.g., the kind of term added to Model A to form Model B), a particular target single-source R-squared difference measure from this framework is used. For instructions on how to identify which target single-source R-squared difference measure to interpret to detect the impact of which kind of term that distinguishes Model A from B, see Rights and Sterba (2020) Table 2. Additionally, this function produces side-by-side graphical comparisons of the R-squared measures for Model A vs. Model B that can be used to visualize changes in each measure across models. This function assumes all level-1 predictors are cluster-mean-centered, for reasons described in Rights & Sterba (2020). Any number of level-1 and/or level-2 predictors is supported and any of the level-1 predictors can have random slopes. This function can be used with either the hierarchical or the simultaneous model-building approach described in Rights and Sterba (2020). This function can be used with either nested or non-nested model comparisons (in which R-squared estimates for Model A are subtracted from those for Model B).

Usage

r2mlm_comp(modelA, modelB, data = NULL, bargraph = TRUE)

Value

If the inputs are valid models, then the output will be a list and associated graphical representation of R-squared decompositions. If the models are not valid, the function will return an error prompting the user to input valid models.

Arguments

modelA, modelB: Models generated using lmer or nlme. Note that models using lmer must specify random effects at the end of the model, like so: outcome ~ 1 + fixed_effects + (random_effects | cluster_variable). Anything else (e.g., outcome ~ 1 + (random_effects | cluster_variable) + fixed_effects) will not work.
data: Optional argument, only needed if models are not hierarchical. Dataset with rows denoting observations and columns denoting variables.
bargraph: Optional bar graph output, default is TRUE.

Details

Assumes that both models are fit with lmer or both models are fit with nlme.

Examples

Run this code

# Using lme4 for your model
# The "bobyqa" optimizer is required for these particular models to converge

if (FALSE)  # Model A, no "salary" components included

modelA_lme4 <- lmer(satisfaction ~ 1 + control_c + control_m + s_t_ratio + (1
+ control_c | schoolID), data = teachsat, REML = TRUE, control =
lmerControl(optimizer = "bobyqa"))

# Model B, full model with "salary" components included

modelB_lme4 <- lmer(satisfaction ~ 1 + salary_c + control_c + salary_m +
control_m + s_t_ratio + (1 + salary_c + control_c | schoolID), data =
teachsat, REML = TRUE, control = lmerControl(optimizer = "bobyqa"))

# Compare models
r2mlm_comp(modelA_lme4, modelB_lme4)

# Compare models, optional data argument specified
r2mlm_comp(modelA_lme4, modelB_lme4, teachsat)

# Using nlme for your model

# Model A, no "salary" components included

modelA_nlme <- lme(satisfaction ~ 1 + control_c + control_m + s_t_ratio,
                  random = ~ 1 + control_c | schoolID,
                  data = teachsat,
                  method = "REML",
                  control = lmeControl(opt = "optim"))

# Model B, full model with "salary" components included

modelB_nlme <- lme(satisfaction ~ 1 + salary_c + control_c + salary_m +
                  control_m + s_t_ratio,
                  random = ~ 1 + salary_c + control_c | schoolID,
                  data = teachsat,
                  method = "REML",
                  control = lmeControl(opt = "optim"))

# Compare models
r2mlm_comp(modelA_nlme, modelB_nlme)

# Compare models, optional data argument specified
r2mlm_comp(modelA_nlme, modelB_nlme, teachsat)