misty (version 0.6.7)

multilevel.r2: R-Squared Measures for Multilevel and Linear Mixed Effects Models

Description

This function computes R-squared measures by Raudenbush and Bryk (2002), Snijders and Bosker (1994), Nakagawa and Schielzeth (2013) as extended by Johnson (2014), and Rights and Sterba (2019) for multilevel and linear mixed effects models estimated by using the lmer() function in the package lme4 or lme() function in the package nlme.

Usage

multilevel.r2(model, print = c("all", "RB", "SB", "NS", "RS"), digits = 3,
              plot = FALSE, gray = FALSE, start = 0.15, end = 0.85,
              color = c("#D55E00", "#0072B2", "#CC79A7", "#009E73", "#E69F00"),
              write = NULL, append = TRUE, check = TRUE, output = TRUE)

Value

Returns an object of class misty.object, which is a list with following entries:

call

function call

type

type of analysis

data

matrix or data frame specified in data

plot

ggplot2 object for plotting the results

args

specification of function arguments

result

list with result tables, i.e., rb for the R2 measures by Raudenbush and Bryk (2002), sb for the R2 measures by Snijders and Bosker (1994), ns for the R2 measures by Nakagawa and Schielzeth (2013), and rs for the R2 measures by Rights and Sterba (2019)

Arguments

model

a fitted model of class "lmerMod" from the lme4 package or "lme" from the nlme package.

print

a character vector indicating which R-squared measures to be printed on the console, i.e., RB for measures from Raudenbush and Bryk (2002), SB for measures from Snijders and Bosker (1994), NS for measures from Nakagawa and Schielzeth (2013) as extended by Johnson (2014), and RS for measures from Rights and Sterba (2019). The default setting is print = "RS".

digits

an integer value indicating the number of decimal places to be used.

plot

logical: if TRUE, bar chart showing the decomposition of scaled total, within-cluster, and between-cluster outcome variance into five (total), three (within-cluster), and two (between-cluster) proportions is drawn. Note that the ggplot2 package is required to draw the bar chart.

gray

logical: if TRUE, graphical parameter to draw the bar chart in gray scale.

start

a numeric value between 0 and 1, graphical parameter to specify the gray value at the low end of the palette.

end

a numeric value between 0 and 1, graphical parameter to specify the gray value at the high end of the palette.

color

a character vector, graphical parameter indicating the color of bars in the bar chart in the following order: Fixed slopes (Within), Fixed slopes (Between), Slope variation (Within), Intercept variation (Between), and Residual (Within). By default, colors from the colorblind-friendly palettes are used

write

a character string naming a text file with file extension ".txt" (e.g., "Output.txt") for writing the output into a text file.

append

logical: if TRUE (default), output will be appended to an existing text file with extension .txt specified in write, if FALSE existing text file will be overwritten.

check

logical: if TRUE (default), argument specification is checked.

output

logical: if TRUE (default), output is shown on the console.

Author

Simon Grund, Alexander Robitzsch, Oliver Luedtk, Mairead Shaw, Jason D. Rights, Sonya K. Sterba, Jessica K. Flake, and Takuya Yanagida

Details

A number of R-squared measures for multilevel and linear mixed effects models have been developed in the methodological literature (see Rights & Sterba, 2018). Based on these measures, following measures were implemented in the current function:

Raudenbush and Bryk (2002)

R-squared measures by Raudenbush and Bryk (2002) are based on the proportional reduction of unexplained variance when predictors are added. More specifically, variance estimates from the baseline/null model (i.e., \(\sigma^2_{e|b}\) and \(\sigma^2_{u0|b}\)) and variance estimates from the model including predictors (i.e., \(\sigma^2_{e|m}\) and \(\sigma^2_{u0|m}\)) are used to compute the proportional reduction in variance between baseline/null model and the complete model by:

$$R^2_1(RB) = \frac{\sigma^2_{e|b} - \sigma^2_{e|m}}{\sigma^2_{e|b}}$$

for the proportional reduction at level-1 (within-cluster) and by:

$$R^2_2(RB) = \frac{\sigma^2_{u0|b} - \sigma^2_{u0|m}}{\sigma^2_{u0|b}}$$

for the proportional reduction at level-2 (between-cluster), where \(|b\) and \(|m\) represent the baseline and full models, respectively (Hox et al., 2018; Roberts et al., 2010).

A major disadvantage of these measures is that adding predictors can increases rather than decreases some of the variance components and it is even possible to obtain negative values for \(R^2\) with these formulas (Snijders & Bosker, 2012). According to Snijders and Bosker (1994) this can occur because the between-group variance is a function of both level-1 and level-2 variance:

$$var(\bar{Y}_j) = \sigma^2_{u0} + \frac{\sigma^2_e}{n_j}$$

Hence, adding a predictor (e.g., cluster-mean centered predictor) that explains proportion of the within-group variance will decrease the estimate of \(\sigma^2_e\) and increase the estimate \(\sigma^2_{u0}\) if this predictor does not explain a proportion of the between-group variance to balance out the decrease in \(\sigma^2_e\) (LaHuis et al., 2014). Negative estimates for \(R^2\) can also simply occur due to chance fluctuation in sample estimates from the two models.

Another disadvantage of these measures is that \(R^2_2(RB)\) for the explained variance at level-2 has been shown to perform poorly in simulation studies even with \(j = 200\) clusters with group cluster size of \(n_j = 50\) (LaHuis et al., 2014; Rights & Sterba, 2019).

Moreover, when there is missing data in the level-1 predictors, it is possible that sample sizes for the baseline and complete models differ.

Finally, it should be noted that R-squared measures by Raudenbush and Bryk (2002) are appropriate for random intercept models, but not for random intercept and slope models. For random slope models, Snijders and Bosker (2012) suggested to re-estimate the model as random intercept models with the same predictors while omitting the random slopes to compute the R-squared measures. However, the simulation study by LaHuis (2014) suggested that the R-squared measures showed an acceptable performance when there was little slope variance, but did not perform well in the presence of higher levels of slope variance.

Snijders and Bosker (1994)

R-squared measures by Snijders and Bosker (1994) are based on the proportional reduction of mean squared prediction error and is computed using the formula:

$$R^2_1(SB) = \frac{\hat{\sigma}^2_{e|m} + \hat{\sigma}^2_{u0|m}}{\hat{\sigma}^2_{e|b} + \hat{\sigma}^2_{u0|b}}$$

for computing the proportional reduction of error at level-1 representing the total amount of explained variance and using the formula:

$$R^2_2(SB) = \frac{\hat{\sigma}^2_{e|m} / n_j + \hat{\sigma}^2_{u0|m}}{\hat{\sigma}^2_{e|b} / n_j + \hat{\sigma}^2_{u0|b}}$$

for computing the proportional reduction of error at level-2 by dividing the \(\hat{\sigma}^2_e\) by the group cluster size \(n_j\) or by the average cluster size for unbalanced data (Roberts et al., 2010). Note that the function uses the harmonic mean of the group sizes as recommended by Snijders and Bosker (1994). The population values of \(R^2\) based on these measures cannot be negative because the interplay of level-1 and level-2 variance components is considered. However, sample estimates of \(R^2\) can be negative either due to chance fluctuation when sample sizes are small or due to model misspecification (Snijders and Bosker, 2012).

When there is missing data in the level-1 predictors, it is possible that sample sizes for the baseline and complete models differ.

Similar to the R-squared measures by Raudenbush and Bryk (2002), the measures by Snijders and Bosker (1994) are appropriate for random intercept models, but not for random intercept and slope models. Accordingly, for random slope models, Snijders and Bosker (2012) suggested to re-estimate the model as random intercept models with the same predictors while omitting the random slopes to compute the R-squared measures. The simulation study by LaHuis et al. (2014) revealed that the R-squared measures showed an acceptable performance, but it should be noted that \(R^2_2(SB)\) the explained variance at level-2 was not investigated in their study.

Nakagawa and Schielzeth (2013)

R-squared measures by Nakagawa and Schielzeth (2013) are based on partitioning model-implied variance from a single fitted model and uses the variance of predicted values of \(var(\hat{Y}_{ij})\) to form both the outcome variance in the denominator and the explained variance in the numerator of the formulas:

$$R^2_m(NS) = \frac{var(\hat{Y}_{ij})}{var(\hat{Y}_{ij}) + \sigma^2_{u0} + \sigma^2_{e}}$$

for marginal total \(R^2_m(NS)\) and:

$$R^2_c(NS) = \frac{var(\hat{Y}_{ij}) + \sigma^2_{u0}}{var(\hat{Y}_{ij}) + \sigma^2_{u0} + \sigma^2_{e}}$$

for conditional total \(R^2_c(NS)\). In the former formula \(R^2\) predicted scores are marginalized across random effects to indicate the variance explained by fixed effects and in the latter formula \(R^2\) predicted scores are conditioned on random effects to indicate the variance explained by fixed and random effects (Rights and Sterba, 2019).

The advantage of these measures is that they can never become negative and that they can also be extended to generalized linear mixed effects models (GLMM) when outcome variables are not continuous (e.g., binary outcome variables). Note that currently the function does not provide \(R^2\) measures for GLMMs, but these measures can be obtained using the r.squaredGLMM() function in the MuMIn package.

A disadvantage is that these measures do not allow random slopes and are restricted to the simplest random effect structure (i.e., random intercept model). In other words, these measures do not fully reflect the structure of the fitted model when using random intercept and slope models. However, Johnson (2014) extended these measures to allow random slope by taking into account the contribution of random slopes, intercept-slope covariances, and the covariance matrix of random slope to the variance in \(Y_{ij}\). As a result, R-squared measures by Nakagawa and Schielzeth (2013) as extended by Johnson (2014) can be used for both random intercept, and random intercept and slope models.

The major criticism of the R-squared measures by Nakagawa and Schielzeth (2013) as extended by Johnson (2014) is that these measures do not decompose outcome variance into each of total, within-cluster, and between-cluster variance which precludes from computing level-specific \(R^2\) measures. In addition, these measures do not distinguish variance attributable to level-1 versus level-2 predictors via fixed effects, and they also do not distinguish between random intercept and random slope variation (Rights and Sterba, 2019).

Rights and Sterba (2019)

R-squared measures by Rights and Sterba (2019) provide an integrative framework of R-squared measures for multilevel and linear mixed effects models with random intercepts and/or slopes. Their measures are also based on partitioning model implied variance from a single fitted model, but they provide a full partitioning of the total outcome variance to one of five specific sources:

• variance attributable to level-1 predictors via fixed slopes (shorthand: variance attributable to f1)

• variance attributable to level-2 predictors via fixed slopes (shorthand: variance attributable to f2)

• variance attributable to level-1 predictors via random slope variation/ covariation (shorthand: variance attributable to v)

• variance attributable to cluster-specific outcome means via random intercept variation (shorthand: variance attributable to m)

• variance attributable to level-1 residuals

\(R^2\) measures are based on the outcome variance of interest (total, within-cluster, or between-cluster) in the denominator, and the source contributing to explained variance in the numerator:

Total \(R^2\) measures

incorporate both within-cluster and between cluster variance in the denominator and quantify variance explained in an omnibus sense:

• \(R^{2(f_1)}_t\): Proportion of total outcome variance explained by level-1 predictors via fixed slopes.

• \(R^{2(f_2)}_t\): Proportion of total outcome variance explained by level-2 predictors via fixed slopes.

• \(R^{2(f)}_t\): Proportion of total outcome variance explained by all predictors via fixed slopes.

• \(R^{2(v)}_t\): Proportion of total outcome variance explained by level-1 predictors via random slope variation/covariation.

• \(R^{2(m)}_t\): Proportion of total outcome variance explained by cluster-specific outcome means via random intercept variation.

• \(R^{2(fv)}_t\): Proportion of total outcome variance explained by predictors via fixed slopes and random slope variation/covariation.

• \(R^{2(fvm)}_t\): Proportion of total outcome variance explained by predictors via fixed slopes and random slope variation/covariation and by cluster-specific outcome means via random intercept variation.

Within-Cluster \(R^2\) measures

incorporate only within-cluster variance in the denominator and indicate the degree to which within-cluster variance can be explained by a given model:

• \(R^{2(f_1)}_w\): Proportion of within-cluster outcome variance explained by level-1 predictors via fixed slopes.

• \(R^{2(v)}_w\): Proportion of within-cluster outcome variance explained by level-1 predictors via random slope variation/covariation.

• \(R^{2(f_1v)}_w\): Proportion of within-cluster outcome variance explained by level-1 predictors via fixed slopes and random slope variation/covariation.

Between-Cluster \(R^2\) measures

incorporate only between-cluster variance in the denominator and indicate the degree to which between-cluster variance can be explained by a given model:

• \(R^{2(f_2)}_b\): Proportion of between-cluster outcome variance explained by level-2 predictors via fixed slopes.

• \(R^{2(m)}_b\): Proportion of between-cluster outcome variance explained by cluster-specific outcome means via random intercept variation.

The decomposition of the total outcome variance can be visualized in a bar chart by specifying plot = TRUE. The first column of the bar chart decomposes scaled total variance into five distinct proportions (i.e., \(R^{2(f_1)}_t\), \(R^{2(f_2)}_t\), \(R^{2(f)}_t\), \(R^{2(v)}_t\), \(R^{2(m)}_t\), \(R^{2(fv)}_t\), and \(R^{2(fvm)}_t\)), the second column decomposes scaled within-cluster variance into three distinct proportions (i.e., \(R^{2(f_1)}_w\), \(R^{2(v)}_w\), and \(R^{2(f_1v)}_w\)), and the third column decomposes scaled between-cluster variance into two distinct proportions (i.e., \(R^{2(f_2)}_b\), \(R^{2(m)}_b\)).

Note that the function assumes that all level-1 predictors are centered within cluster (i.e., group-mean or cluster-mean centering) as has been widely recommended (e.g., Enders & Tofighi, D., 2007; Rights et al., 2019). In fact, it does not matter whether a lower-level predictor is merely a control variable, or is quantitative or categorical (Yaremych et al., 2021), cluster-mean centering should always be used for lower-level predictors to obtain an orthogonal between-within partitioning of a lower-level predictor's variance that directly parallels what happens to a level-1 outcome (Hoffman & Walters, 2022). In the absence of cluster-mean-centering, however, the function provides total \(R^2\) measures, but does not provide any within-cluster or between-cluster \(R^2\) measures.

By default, the function only computes R-squared measures by Rights and Sterba (2019) because the other R-squared measures reflect the same population quantity provided by Rights and Sterba (2019). That is, R-squared measures \(R^2_1(RB)\) and \(R^2_2(RB)\) by Raudenbush and Bryk (2002) are equivalent to \(R^{2(f_1v)}_w\) and \(R^{2(f_2)}_b\), R-squared measures \(R^2_1(SB)\) and \(R^2_2(SB)\) are equivalent to \(R^{2(f)}_t\) and \(R^{2(f_2)}_b\), and R-squared measures \(R^2_m(NS)\) and \(R^2_c(NS)\) by Nakagawa and Schielzeth (2013) as extended by Johnson (2014) are equivalent to \(R^{2(f)}_t\) and \(R^{2(fvm)}_t\) (see Rights and Sterba, Table 3).

Note that none of these measures provide an \(R^2\) for the random slope variance explained by cross-level interactions, a quantity that is frequently of interest (Hoffman & Walters, 2022).

References

Enders, C. K., & Tofighi, D. (2007). Centering predictor variables in cross-sectional multilevel models: A new look at an old issue. Psychological Methods, 12, 121-138. https://doi.org/10.1037/1082-989X.12.2.121

Hoffmann, L., & Walter, W. R. (2022). Catching up on multilevel modeling. Annual Review of Psychology, 73, 629-658. https://doi.org/10.1146/annurev-psych-020821-103525

Hox, J., Moerbeek, M., & van de Schoot, R. (2018). Multilevel Analysis: Techniques and Applications (3rd ed.) Routledge.

Johnson, P. C. D. (2014). Extension of Nakagawa & Schielzeth’s R2 GLMM to random slopes models. Methods in Ecology and Evolution, 5(9), 944-946. https://doi.org/10.1111/2041-210X.12225

LaHuis, D. M., Hartman, M. J., Hakoyama, S., & Clark, P. C. (2014). Explained variance measures for multilevel models. Organizational Research Methods, 17, 433-451. https://doi.org/10.1177/1094428114541701

Nakagawa, S., & Schielzeth, H. (2013). A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods in Ecology and Evolution, 4(2), 133-142. https://doi.org/10.1111/j.2041-210x.2012.00261.x

Raudenbush, S. W., & Bryk, A. S., (2002). Hierarchical linear models: Applications and data analysis methods. Sage.

Rights, J. D., Preacher, K. J., & Cole, D. A. (2020). The danger of conflating level-specific effects of control variables when primary interest lies in level-2 effects. British Journal of Mathematical and Statistical Psychology, 73(Suppl 1), 194-211. https://doi.org/10.1111/bmsp.12194

Rights, J. D., & Sterba, S. K. (2019). Quantifying explained variance in multilevel models: An integrative framework for defining R-squared measures. Psychological Methods, 24, 309-338. https://doi.org/10.1037/met0000184

Roberts, K. J., Monaco, J. P., Stovall, H., & Foster, V. (2011). Explained variance in multilevel models (pp. 219-230). In J. J. Hox & J. K. Roberts (Eds.), Handbook of Advanced Multilevel Analysis. Routledge.

Snijders, T. A. B., & Bosker, R. (1994). Modeled variance in two-level models. Sociological methods and research, 22, 342-363. https://doi.org/10.1177/0049124194022003004

Snijders, T. A. B., & Bosker, R. J. (2012). Multilevel analysis: An introduction to basic and advanced multilevel modeling (2nd ed.). Sage.

Yaremych, H. E., Preacher, K. J., & Hedeker, D. (2021). Centering categorical predictors in multilevel models: Best practices and interpretation. Psychological Methods. Advance online publication. https://doi.org/10.1037/met0000434

See Also

multilevel.cor, multilevel.descript, multilevel.icc, multilevel.indirect

Examples

if (FALSE) {
# Load misty, lme4, nlme, and ggplot2 package
library(misty)
library(lme4)
library(nlme)
library(ggplot2)

# Load data set "Demo.twolevel" in the lavaan package
data("Demo.twolevel", package = "lavaan")

#----------------------------------------------------------------------------
#'
# Cluster mean centering, center() from the misty package
Demo.twolevel$x2.c <- center(Demo.twolevel$x2, type = "CWC",
                             cluster = Demo.twolevel$cluster)

# Compute group means, cluster.scores() from the misty package
Demo.twolevel$x2.b <- cluster.scores(Demo.twolevel$x2,
                                     cluster = Demo.twolevel$cluster)

# Estimate multilevel model using the lme4 package
mod1a <- lmer(y1 ~ x2.c + x2.b + w1 + (1 + x2.c | cluster), data = Demo.twolevel,
              REML = FALSE, control = lmerControl(optimizer = "bobyqa"))

# Estimate multilevel model using the nlme package
mod1b <- lme(y1 ~ x2.c + x2.b + w1, random = ~ 1 + x2.c | cluster, data = Demo.twolevel,
             method = "ML")

#----------------------------------------------------------------------------
#'
# Example 1a: R-squared measures according to Rights and Sterba (2019)
multilevel.r2(mod1a)
#'
# Example 1b: R-squared measures according to Rights and Sterba (2019)
multilevel.r2(mod1b)
#'
# Example 1a: Write Results into a text file
multilevel.r2(mod1a, write = "ML-R2.txt")

#-------------------------------------------------------------------------------

# Example 2: Bar chart showing the decomposition of scaled total, within-cluster,
# and between-cluster outcome variance
multilevel.r2(mod1a, plot = TRUE)

# Bar chart in gray scale
multilevel.r2(mod1a, plot = TRUE, gray = TRUE)

# Save bar chart, ggsave() from the ggplot2 package
ggsave("Proportion_of_Variance.png", dpi = 600, width = 5.5, height = 5.5)

#-------------------------------------------------------------------------------

# Example 3: Estimate multilevel model without random slopes
# Note. R-squared measures by Raudenbush and Bryk (2002), and  Snijders and
# Bosker (2012) should be computed based on the random intercept model
mod2 <- lmer(y1 ~ x2.c + x2.b + w1 + (1 | cluster), data = Demo.twolevel,
             REML = FALSE, control = lmerControl(optimizer = "bobyqa"))

# Print all available R-squared measures
multilevel.r2(mod2, print = "all")

#-------------------------------------------------------------------------------

# Example 4: Draw bar chart manually
mod1a.r2 <- multilevel.r2(mod1a, output = FALSE)

# Prepare data frame for ggplot()
df <- data.frame(var = factor(rep(c("Total", "Within", "Between"), each = 5),
                              level = c("Total", "Within", "Between")),
                 part = factor(c("Fixed Slopes (Within)", "Fixed Slopes (Between)",
                                 "Slope Variation (Within)", "Intercept Variation (Between)",
                                 "Residual (Within)"),
                 level = c("Residual (Within)", "Intercept Variation (Between)",
                           "Slope Variation (Within)", "Fixed Slopes (Between)",
                           "Fixed Slopes (Within)")),
                 y = as.vector(mod1a.r2$result$rs$decomp))

# Draw bar chart in line with the default setting of multilevel.r2()
ggplot(df, aes(x = var, y = y, fill = part)) +
  theme_bw() +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c("#E69F00", "#009E73", "#CC79A7", "#0072B2", "#D55E00")) +
  scale_y_continuous(name = "Proportion of Variance", breaks = seq(0, 1, by = 0.1)) +
  theme(axis.title.x = element_blank(),
        axis.ticks.x = element_blank(),
        legend.title = element_blank(),
        legend.position = "bottom",
        legend.box.margin = margin(-10, 6, 6, 6)) +
  guides(fill = guide_legend(nrow = 2, reverse = TRUE))
}