Motivation

Regression explains variation in an outcome (dependent variable) using one or more predictors (independent variables).

A regression coefficient quantifies the strength of the relationship between outcome and predictor (effect size), but in units specific to the predictor, so comparison to other predictors’ coefficients is not straightforward.

Standardized effect sizes, like Cohen’s \(d\) and \(R^2\), are unit-less and thus permit direct comparison between predictors measured in different units, even across different studies.

Multilevel regression models are commonly used to model data that are hierarchically structured, in which observations collected within aggregate units are not independent.

However, currently no set of standardized effect size measures have been widely adopted by multilevel modelers.

Purpose

This workshop discusses a new set of \(R^2\) measures introduced by Rights and Sterba in their 2019 article ” Quantifying Explained Variance in Multilevel Models: An Integrative Framework for Defining R-Squared Measures”.

The authors have themselves said that using their new \(R^2\) measures requires a solid understanding of multilevel models (Shaw et al., 2023).

To this end, we will discuss:

• Linear regression review
• Multilevel Models
• New \(R^2\) measures for multilevel models
• Using the r2mlm package

Understanding variation

How is variation measured?

Imagine we have an outcome, \(y_i\), measured for units \(i=1,...,N\).

One common way to quantify variation is to measure deviations from the mean, \(y_i - \bar{y}\), where \(\bar{y}\) is the mean of \(y_i\).

3 measures of variation:

• sums of squares: sum of the squared deviations from the mean
• variance is the average squared deviation from the mean
• standard deviation, the square root of variance, is the average deviation from the mean

Linear regression

Linear regression uses a set of predictors to explain the variation in the outcome.

Below we model outcome \(y_i\) to be predicted by a single predictor, \(x_i\):

\[y_i = \beta_0 + \beta_1x_i + \epsilon_i\] where

• \(\beta_0\): intercept
• \(\beta_1\): coefficient quantifying relationship between \(x_i\) and \(y_i\)
• \(\epsilon_i\): error term, representing unexplained variation

The errors \(\epsilon_i\) are assumed to be normally distributed with a zero mean and variance \(\sigma_{\epsilon}^2\):

\[\epsilon_i \sim N(0, \sigma_{\epsilon}^2)\]

\(\sigma_{\epsilon}^2\), the error variance, quantifies how much \(y_i\) variation is unexplained by the model.

As relevant predictors are added to the model, \(\sigma_{\epsilon}^2\) should decrease.

Prediction

After fitting a model, we can estimate predicted outcome for each observations, \(\hat{y}_i\), based on the estimated coefficients.

In the intercept-only model, all predictions equal the intercept, which equals \(\bar{y}\) in this model:

\[\begin{align} \hat{y}_i &= \hat{\beta}_0 \\ &= \bar{y} \end{align} \]

For a model with predictor \(x_i\), the predictions lie on a line determined by \(\hat{\beta_0}\) and \(\hat{\beta}_1\):

\[\hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1x_i\]

Explained variation

The intercept-only model explains no variation in \(y_i\).

As predictors are added, the predictions \(\hat{y}_i\) on average move away from the sample mean \(\bar{y}\) and closer to their observed values \(y_i\).

The distance that predictions move away from the sample mean towards the observed values measures how much the model is explaining variation.

A steeper slope will move the predicted values away more from the sample mean, explaining more variation.

Residuals

Residuals, \(\hat{\epsilon}_i\), are the difference between the observed and predicted values and estimate the model errors, \(\epsilon_i\):

\[\hat{\epsilon}_i = y_i - \hat{y}_i\]

Residuals represent the variation in the outcome that cannot be explained by the model.

Residuals shrink as explained variation increases.

Below we see how the residual variance (expressed as a standard deviation) shrinks when we add a predictor to the model.

# intercept-only model (no predictors)
#   residual standard error is estimate of square root of error variance
m0 <- lm(y ~ 1, data=small)
summary(m0)

## 
## Call:
## lm(formula = y ~ 1, data = small)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.556  -6.556   1.444   5.694  16.444 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   50.556      1.979   25.55 5.31e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.396 on 17 degrees of freedom

# Adding predictor x to the model reduces residual standard error
m <- lm(y ~ x, data=small)
summary(m)

## 
## Call:
## lm(formula = y ~ x, data = small)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.113  -7.009   1.564   5.941  11.244 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  33.1008    11.6468   2.842   0.0118 *
## x             0.3574     0.2353   1.519   0.1482  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.091 on 16 degrees of freedom
## Multiple R-squared:  0.1261, Adjusted R-squared:  0.07144 
## F-statistic: 2.308 on 1 and 16 DF,  p-value: 0.1482

Partitioning the total variation

In regression, we can partition total variation into explained and residual variation:

\[\text{Total Variation = Explained Variation + Residual Variation}\]

As we add relevant predictors to the model, explained variation grows while residual variation shrinks.

\(R^2\)

\(R^2\) expresses the proportion of \(\text{Total Variation}\) that is \(\text{Explained Variation}\):

Most \(R^2\) measures can be expressed generally as:

\[R^2 = \frac{\text{Explained Variation}}{\text{Total Variation}}\] \(R^2\) is thus interpreted as the proportion of outcome variation explained by the model predictors.

Sometimes, residual variation is easier to estimate than explained variation, so some \(R^2\) are defined as :

\[R^2 = 1 - \frac{\text{Residual Variation}}{\text{Total Variation}}\]

Different \(R^2\) measures use different statistics (e.g. variances, sums of squares) to represent variation in the numerator and denominator.^†

Because the intercept-only model explains no variation, \(R^2=0\).

# R^2 for intercept-only model
summary(lm(y ~ 1, data=small))$r.squared

## [1] 0

The predictor \(x_i\) explains about 12.6% of the variation in \(y_i\).

# R^2 for x
summary(lm(y ~ x, data=small))$r.squared

## [1] 0.1260599

Linear regression \(R^2\) properties

• allows comparison between predictors and between models as a standardized effect size
• varies between 0 and 1, with values closer to 1 implying more explained variation
• only increases (or stays same) as more predictors are added to the model

^†\(R^2\) for linear regression uses sums of squares for explained and total variation

Variance explained explained

What determines how much variance is explained by a predictor?

Holding everything else constant, variance explained increases as:

• the magnitude of the regression coefficient increases
  • sex explains more variation in height than whether someone is born in the first or second half of the year (which presumably has no relationship to height)
• the variance of the predictor increases
  • sex (50% male=1) and a rare disease (4% disease=1) may have effects of similar size on height, but sex explains more variation because it varies more

# height data
height[4:10,]

##      height male before_july disease
## 4  70.33877    0           1       1
## 5  62.55333    0           0       0
## 6  64.84237    0           1       0
## 7  61.73652    0           0       0
## 8  64.03135    0           1       0
## 9  64.42657    0           0       0
## 10 66.27407    0           1       0

# male and disease have similar coefficients
# before_july has nearly zero effect
summary(lm(height ~ male + before_july + disease, data=height))

## 
## Call:
## lm(formula = height ~ male + before_july + disease, data = height)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6396 -1.7393  0.0938  1.4779  6.9888 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  64.2630     0.2760 232.816  < 2e-16 ***
## male          4.8323     0.3165  15.266  < 2e-16 ***
## before_july   0.1727     0.3165   0.546    0.586    
## disease       4.9490     0.8077   6.128 4.81e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.238 on 196 degrees of freedom
## Multiple R-squared:  0.5802, Adjusted R-squared:  0.5738 
## F-statistic:  90.3 on 3 and 196 DF,  p-value: < 2.2e-16

Below we graph the individual contributions to \(R^2\) by predictor:

Multilevel data

Multilevel data are hierarchically structured where sampled observations (level-1) can be aggregated into higher level units (level-2).

Multilevel data arise when:

• individuals are sampled in clusters (e.g. students within schools, patients within hospitals)
  • individuals are level-1 and clusters are level-2
• subjects are repeatedly measured (i.e., longitudinal data)
  • timepoints are level-1 and subjects are level-2

Level-2 units are assumed to be independent of each other but level-1 observations within the same level-2 unit are not assumed to be independent, violating the independence assumption of single-level models.

Note: Henceforth, we will refer to multilevel data arising from clustered sampling but all concepts discussed apply equally to repeated measures data as well.

Example multilevel data

We will discuss multilevel model with a simulated dataset of worker stress, with 15 companies sampled and 10 employess within each company sampled.

Variables

• company_id
• employee_id
• stress: the outcome, ranges 1-100
• salary: employee’s salary (level-1) in thousands of dollars
• pct_open: percent of company’s position that are unfilled (level-2)

The dataset is called empstress.

# empstress data
head(empstress)

##   company_id employee_id salary pct_open stress
## 1          1           1 76.594   11.082     48
## 2          1           2 80.259   11.082     41
## 3          1           3 67.790   11.082     43
## 4          1           4 82.514   11.082     47
## 5          1           5 76.245   11.082     54
## 6          1           6 78.169   11.082     41

Multilevel model notation

The multilevel model explains variation in outcome \(y_{ij}\), measured for individual \(i=1,...,n_j\) within level-2 unit \(j=1,...,J\).

Each of the \(J\) level-2 unit has \(n_j\) individuals sampled from it, and the total sample size is \(N =\sum_{j=1}^J{n_j}\).

Predictors in multilevel models can be measured at level-1 or level-2.

• Level-1 predictors, subscripted with \(_{ij}\), vary across individuals
• Level-2 predictors, subscripted with \(_j\), only vary across level-2 units and are constant across level-1 units within the same level-2 unit

For the empstress data

• We will model variation in \(stress_{ij}\), measured for employee \(i=1,...,n_j\) within company \(j=1,...,J\)
• \(J=15\) companies were sampled, and \(n_j=10\) employees were sampled for each company, resulting in \(N=150\)
• level-1 predictor \(salary_{ij}\) varies across employees
• level-2 predictor \(pct\_open_{j}\) only varies across companies

# empstress
head(empstress)

##   company_id employee_id salary pct_open stress
## 1          1           1 76.594   11.082     48
## 2          1           2 80.259   11.082     41
## 3          1           3 67.790   11.082     43
## 4          1           4 82.514   11.082     47
## 5          1           5 76.245   11.082     54
## 6          1           6 78.169   11.082     41

Means in multilevel data

For variables measured at level-1 like \(y_{ij}\), we can calculate 2 kinds of means:

• the overall (sample) mean \(y_{..}\)

\[\bar{y}_{..} = \frac{1}{N}\sum_{j=1}^J{\sum_{i=1}^{n_j}{y_{ij}}}\]

• the cluster means \(y_{.j}\) for each cluster \(j\):

\[\bar{y}_{.j} = \frac{1}{n_j}\sum_{i=1}^{n_j}{y_{ij}}\]

Below we see a graph of stress scores (blue points), its overall mean (line) and company means (larger red points).

Total variation decomposition in multilevel models

We can decompose a level-1 outcome \(y_{ij}\) into the sum of a cluster mean \(y_{.j}\) and the deviation from that cluster mean:

\[\begin{align} y_{ij} &= \bar{y}_{.j} + (y_{ij} - \bar{y}_{.j}) \\ y_{ij} &= \text{cluster mean + deviation from cluster mean} \end{align}\]

We can estimate the variation of each of the terms in the equation above:

• \(var(y_{ij})\) represents \(\text{Total Variation}\)
• \(var(\text{cluster mean})\) represents \(\text{Between-Cluster Variation}\)
• \(var(\text{deviation from cluster mean})\) represents \(\text{Within-Cluster Variation}\)

The total variation of \(y_{ij}\) in a multilevel model can thus be decomposed as:

\[\text{Total Variation = Between-Cluster Variation + Within-Cluster Variation}\]

Importantly, both \(\text{Between-Cluster Variation}\) and \(\text{Within-Cluster Variation}\) can decomposed into explained and residual variation:

\[\text{Between-Cluster Variation = Explained Between-Cluster Variation + Residual Between-Cluster Variation}\]

\[\text{Within-Cluster Variation = Explained Within-Cluster Variation + Residual Within-Cluster Variation}\]

Within-cluster variation

Describes variation of individuals within the same cluster (e.g. employees’ stress scores varying within a company).

Level-1 predictors explain within-cluster variation in the outcome (e.g. differences in employees’ salaries within a company explain differences in their stress scores).

Between-cluster variation

Describes variation of cluster means about the overall mean (e.g. company stress score means varying around the overall mean).

Level-2 predictors explain between-cluster variation in the outcome (e.g. differences in percent unfilled positions can explain differences in mean stress across companies)

Random intercept model

The random intercept model has no predictors and thus decomposes the total variation into residual within-cluster and between-cluster variation.

\[\text{Total Variation}_{ri} = \text{Residual Between-Cluster Variation + Residual Within-Cluster Variation}\]

\[\begin{align}y_{ij} &= \beta_{0j} + \epsilon_{ij} \\ \beta_{0j} &= \gamma_{00} + u_{0j} \end{align}\]

In the random intercept model:

• \(\beta_{0j}\) is the intercept for cluster \(j\) and will be an estimate of the cluster mean^†
• \(\epsilon_{ij}\) is the level-1 error, individual \(i\)’s deviation from the predicted score, here \(\hat{y}_{ij} = \hat{\beta}_{0j}\) (unexplained within-cluster variation)
• \(\gamma_{00}\) is the fixed intercept, the mean of all clusters’ intercepts
• \(u_{0j}\) is the random intercept for cluster \(j\), the deviation of cluster \(j\)’s intercept, \(\beta_{0j}\), from the mean intercept, \(\gamma_{00}\) (unexplained between-cluster variation)

# random intercept model 
m_ri <- lmer(stress ~ 1 + (1|company_id), data=empstress)
summary(m_ri)

## Linear mixed model fit by REML ['lmerMod']
## Formula: stress ~ 1 + (1 | company_id)
##    Data: empstress
## 
## REML criterion at convergence: 1107.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2379 -0.6096 -0.0342  0.5618  3.3270 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  company_id (Intercept) 116.08   10.774  
##  Residual                73.41    8.568  
## Number of obs: 150, groups:  company_id, 15
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   45.000      2.868   15.69

In the output above, we see:

• random intercept variance estimate \(\hat{\tau}_{00}^2\)=116.08, representing the total amount of \(\text{Residual Between-Cluster Variation}\)
• residual variance estimate \(\hat{\sigma}_{\epsilon}^2\)=73.41, representing the total amount of \(\text{Residual Within-Cluster Variation}\)

Overall, \(stress_{ij}\) varies more between-companies than within-companies.

^†Random effect estimates are partially pooled and shrunken towards the grand mean and will generally not equal the observed cluster means \(\bar{y}_{.j} = \sum_{i=1}^{n_j}{y_{ij}}\).

Intraclass correlation coefficient

The random intercept model is often used to calculate the intraclass correlation coefficient (ICC), the proportion of total variation in \(y_{ij}\) that is residual between-cluster variation.

The ICC is calculated by dividing \(\hat{\tau}_{00}^2\), the random intercept variance, by \(\hat{\tau}_{00}^2 + \hat{\sigma}_{\epsilon}^2\), the sum of the random intercept variance and residual variance, which estimates the total variance.

\[ICC = \frac{\hat{\tau}_{00}^2}{\hat{\tau}_{00}^2 + \hat{\sigma}_{\epsilon}^2}\]

\(ICC>.5\) suggests that cluster means vary more than individual scores within cluster, or \(\hat{\tau}_{00}^2 > \hat{\sigma}_{\epsilon}^2\). Here, level-2 predictors will be important to explain variation in \(y_ij\).

\(ICC<.5\) suggests that cluster means vary less than individual scores within cluster, or \(\hat{\tau}_{00}^2 < \hat{\sigma}_{\epsilon}^2\). Here, level-1 predictors will be important to explain variation in \(y_{ij}\).

Within-cluster centering

Like \(y_{ij}\), a predictor measured at level-1 can be decomposed into a cluster mean and a deviation from the cluster mean^†:

\[\begin{align} x_{ij} &= \bar{x}_{.j} + (x_{ij} - \bar{x}_{.j}) \\ x_{ij} &= \text{cluster mean + deviation from cluster mean} \end{align}\]

A variable superscripted with \(^w\), e.g., \(x_{ij}^w\), signifies its deviation from the cluster mean, and is often called group-mean-centered or cluster-mean-centered:

\[x_{ij}^w = (x_{ij} - \bar{x}_{.j})\]

and thus

\[x_{ij} = \bar{x}_{.j} + x_{ij}^w\]

Below we use mutate() to create both the company means of salary, salary_m, and the company-mean-centered salary, salary_c.

# to create group-mean centered salary, salary_c:
# 1. add cluster means of predictor salary (salary_m) to data
# 2. subtract salary_m from original salary to create salary_c
empstress <- empstress %>% 
  mutate(.by=company_id, # process data by company_id
         salary_m = mean(salary)) %>%  
  mutate(salary_c = salary-salary_m) # separate mutate so we don't process by company

As a predictor, \(x_{ij}^w\) can only explain level-1 variation because it does not vary at level-2, as all between-cluster differences have been removed (all cluster means are 0):

The cluster means \(x_{.j}\) can themselves be added as level-2 predictors to explain level-2 variation in the outcome.

^†If the cluster means of \(x_{ij}\) are all equal, which can happen with experimentally manipulated variables (e.g. using time as a predictor when everyone is measured at the same times), then \(x_{ij}\) will only explain level-1 variation so within-cluster centering is not necessary

Within-cluster and between-cluster effects

Recall that

\[x_{ij} = x_{ij}^w + \bar{x}_{.j}\]

In the model in which \(x_{ij}^w\) and \(\bar{x}_{.j}\) are predictors:

\[y_{ij} = \beta_{0j} + \beta_1x_{ij}^w + \beta_2\bar{x}_{.j} + \epsilon_{ij}\]

• \(\beta_1\) is the within-cluster effect of \(x_{ij}\), (e.g., differences in salary between employees within the same company explains differences in their stress scores)
• \(\beta_2\) is the between-cluster effect of \(x_{ij}\), (e.g. differences in mean salary across companies explains differences in mean stress across companies)

Importantly, the above model is generally not equivalent to the model in which the original \(x_{ij}\) is the predictor:

\[y_{ij} = \beta_{0j}^* + \beta_1^*x_{ij} + \epsilon_{ij}\]

Implying that:

\[\begin{align}\beta_1^* &\ne \beta_1 \\ \beta_1^* &\ne \beta_2\end{align}\]

Instead, the coefficient for the original \(x_{ij}\), \(\beta_1^*\), will be some uninterpretable combination of the within-cluster and between-cluster effects, \(\beta_1\) and \(\beta_2\), respectively.

• If \(\beta_1\) and \(\beta_2\) have opposite signs, \(\beta_1^*\) may be close to \(0\) and thus modeling \(x_{ij}\) may lead to a misleading inference of little or no effect
• Thus, it is often recommended that level-1 predictors be decomposed into their within (\(x_{ij}^w\)) and between (\(\bar{x}_{.j}\)) components so that coefficients are more cleanly interpretable.

Compare the coefficients of the model where group_centered salary_c and company means salary_m are entered as predictor to the model where the original salary is entered:

# salary_c and salary_m coefficients have clear within-company and between-company interpretations
m_salary_c_m <- lmer(stress ~ salary_c + salary_m +
                       (1 | company_id),
                     data=empstress)

summary(m_salary_c_m)

## Linear mixed model fit by REML ['lmerMod']
## Formula: stress ~ salary_c + salary_m + (1 | company_id)
##    Data: empstress
## 
## REML criterion at convergence: 1072
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.14261 -0.58580  0.02166  0.46986  3.14261 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  company_id (Intercept) 22.71    4.766   
##  Residual               64.73    8.046   
## Number of obs: 150, groups:  company_id, 15
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 125.7257    11.9571  10.515
## salary_c      0.7414     0.1696   4.371
## salary_m     -1.0053     0.1479  -6.798
## 
## Correlation of Fixed Effects:
##          (Intr) slry_c
## salary_c  0.000       
## salary_m -0.993  0.000

# salary coefficient is uninterpretable average of coefficients above
m_salary <- lmer(stress ~ salary +
                   (1 | company_id),
                 data=empstress)

summary(m_salary)

## Linear mixed model fit by REML ['lmerMod']
## Formula: stress ~ salary + (1 | company_id)
##    Data: empstress
## 
## REML criterion at convergence: 1102.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.23736 -0.62578  0.00259  0.53223  3.10754 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  company_id (Intercept) 239.80   15.485  
##  Residual                65.62    8.101  
## Number of obs: 150, groups:  company_id, 15
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   4.6013    13.3741   0.344
## salary        0.5031     0.1587   3.170
## 
## Correlation of Fixed Effects:
##        (Intr)
## salary -0.953

Fixed effects

Multilevel models can include fixed and random effects.

The intercept and any regression coefficient for a level-1 predictor can be either

• constrained to be the same across all clusters (fixed effect only)
  • e.g., \(\beta_1 = \gamma_{10}\)
• allowed to vary across clusters (fixed + random effect)
  • e.g., \(\beta_{1j} = \gamma_{10} + u_{1j}\)

In either case, the fixed effect, e.g. \(\gamma_{10}\), is interpreted as the mean coefficient across clusters.

# fixed slope of salary_c to explain variation in stress
m_salary_c_fixed <- lmer(stress ~ salary_c + 
                           (1 | company_id),
                         data=empstress)

summary(m_salary_c_fixed)

## Linear mixed model fit by REML ['lmerMod']
## Formula: stress ~ salary_c + (1 | company_id)
##    Data: empstress
## 
## REML criterion at convergence: 1091.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.22676 -0.61828 -0.01386  0.49847  3.00557 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  company_id (Intercept) 116.95   10.814  
##  Residual                64.73    8.046  
## Number of obs: 150, groups:  company_id, 15
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  45.0000     2.8685  15.688
## salary_c      0.7414     0.1696   4.371
## 
## Correlation of Fixed Effects:
##          (Intr)
## salary_c 0.000

The coefficient \(\beta_1 = \gamma_{10} \approx 0.741\) is interpreted as the fixed slope of group-mean centered salary, salary_c, which suggests that employees with higher salaries than other employees within the same company have more stress (perhaps because they have more responsibilities).

Random effects

Random effects in multilevel models are typically used to model heterogeneity (variation) around the fixed effect

• Random intercepts model heterogeneity in the cluster means of the outcome
• Random slopes model heterogeneity in the coefficients for level-1 predictors

In the model below with level-1 predictors \(x_{ij}^w\) and \(z_{ij}^w\):

\[\begin{align}y_{ij} &= \beta_{0j} + \beta_{1j}x_{ij}^w + \beta_2z_{ij}^w + \epsilon_{ij} \\ \beta_{0j} &= \gamma_{00} + u_{0j} \\ \beta_{1j} &= \gamma_{10} + u_{1j} \\ \beta_2 &= \gamma_{20} \end{align}\]

• \(\beta_{0j}\) is the intercept for cluster \(j\), and is the sum of the fixed (mean) intercept \(\gamma_{00}\) and each cluster’s random intercept \(u_{0j}\)
• \(\beta_{1j}\) is the slope coefficient for \(x_{ij}^w\) for cluster \(j\), and is the sum of the fixed mean slope \(\gamma_{10}\) and each cluster’s random slope \(u_{1j}\)
• \(\beta_{2}\) is the slope coefficient for \(z_{ij}^w\) for all clusters as there is only a fixed component, \(\gamma_{20}\)

Random effects like \(u_{0j}\) and \(u_{1j}\) are typically assumed to have a multivariate normal distribution, each with zero means and a covariance matrix quantifying their variances and covariances:

\[\Big(\begin{matrix} u_{0j} \\ u_{1j} \end{matrix}\Big) \sim N\Bigg( \Big(\begin{matrix} 0 \\ 0 \end{matrix}\Big), \Big(\begin{matrix} \tau_{00}^2 & \tau_{01} \\ \tau_{01} & \tau_{11}^2 \end{matrix}\Big) \Bigg)\]

The elements of the covariance matrix \(\tau_{00}^2\), \(\tau_{11}^2\), \(\tau_{01}=\tau_{10}\) are parameters estimated in the multilevel model.

# fixed and random slope of salary_c to explain variation in stress
m_salary_c_random <- lmer(stress ~ salary_c + 
                           (1 + salary_c | company_id),
                         data=empstress)

summary(m_salary_c_random)

## Linear mixed model fit by REML ['lmerMod']
## Formula: stress ~ salary_c + (1 + salary_c | company_id)
##    Data: empstress
## 
## REML criterion at convergence: 1087.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.31751 -0.58010 -0.05927  0.47469  2.48187 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev. Corr
##  company_id (Intercept) 117.4425 10.837       
##             salary_c      0.3329  0.577   0.19
##  Residual                59.3291  7.703       
## Number of obs: 150, groups:  company_id, 15
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  45.0000     2.8679  15.691
## salary_c      0.8178     0.2245   3.643
## 
## Correlation of Fixed Effects:
##          (Intr)
## salary_c 0.126

• fixed slope coefficient \(\gamma_{10} \approx 0.818\) is interpreted as the mean slope of salary_c across companies.
• The variance of the slopes is estimated at \(tau_{11}^2 \approx 0.333\), suggesting that a company’s slope \(\beta_{1j}\) deviates from the fixed slope \(\gamma_{10}\) by about \(\sqrt{0.333} \approx 0.577\).

# use coef() to see each cluster's intercepts and slopes
coef(m_salary_c_random)

## $company_id
##    (Intercept)    salary_c
## 1     46.17442  0.48449509
## 2     27.79865  0.79443093
## 3     41.50957  0.95495935
## 4     31.01874  0.73785965
## 5     53.72239  0.71527322
## 6     51.56692  0.87882989
## 7     44.06388  0.89459952
## 8     52.49329  1.26991889
## 9     62.81527  1.60314719
## 10    39.82638  0.58183403
## 11    39.10364  0.78680564
## 12    27.18478  0.55404316
## 13    55.85613 -0.08883465
## 14    50.70981  0.86780080
## 15    51.15612  1.23142129
## 
## attr(,"class")
## [1] "coef.mer"

Random slopes may provide a better fit of level-1 predictors to the data than the fixed slope.

Excercise 1 dataset

To review multilevel modeling and practice using the r2mlm package, we will use a dataset that simulates the following experiment:

Subjects are college students who were asked to study a different list of 100 words every night for 2 weeks for between 0 and 2 hours per night, at their discretion, and then asked to recall as many words as possible the next morning.

30 subjects, 14 measurements each

Variables

• subj_id: subject ID
• words: the outcome, number of words recalled each morning
• sleep: the number of hours the subject slept the night before
• study: the number of hours the subject studied the night before
• gpa: subject’s grade point average (level 2)

Exercise 1

1. Run the random intercept model with words as the outcome and random intercepts by subj_id.
2. Calculate the ICC for words.
3. Use dplyr::mutate() to add the following to the data set:
   • cluster means for sleep and study (call them sleep_m and study_m)
   • subject-mean-centered sleep and subject-mean-centered study (call them sleep_c and study_c)
4. Run a multilevel model with fixed and random slopes for sleep_c
   • interpret the coefficient
   • take note of the changes in residual and random intercept variances from the random intercept model.
5. Add sleep_m to the model in (4)
   • interpret its coefficient,
   • and take note of the changes in variances again.

Total variance decomposition in multilevel models

Rights and Sterba (2019) decompose the model-implied total variance of the outcome as follows:

\[\widehat{var}(y_{ij}) = f_1 + f_2 + v + m + \hat{\sigma}_\epsilon^2\] Explained variance components

• \(f_1\) is variance explained by level-1 fixed effects^†
• \(f_2\) is variance explained by level-2 fixed effects
• \(v\) is variance explained by random (varying) slopes
• \(m\) is random intercept (mean) variance (\(\hat{\tau}_{00}^2\))^††

Residual variance components

• \(\hat{\sigma}_{\epsilon}^2\) is level-1 residual variance

This complete variance decomposition allowed Rights and Sterba (2019) to develop a general framework for defining multilevel \(R^2\) measures that are clearly interpretable and that subsumes most multilevel \(R^2\) measures previously defined in the literature.

^†In this decomposition, all level-1 predictors are assumed to be cluster-mean-centered

^†† Up until now we have interpreted random intercept variation to be residual or unexplained between-cluster variation, but the r2mlm pacakge interprets it as variation explained by cluster IDs.

Within-cluster variance decomposition

Within this framework, we can also fully decompose the model-implied within-cluster variance decomposition:

\[\widehat{var}(y_{ij}-\bar{y}_{.j}) = f_1 + v + \hat{\sigma}_\epsilon^2\]

Explained variance components

• \(f_1\): (group-centered) level-1 predictor fixed (mean) effects explain within-cluster differences
• \(v\) random (varying) slopes may provide a closer fit to outcomes than the fixed slope and thus explain more variance

Residual variance component

• \(\hat{\sigma}_{\epsilon}^2\) quantifies unexplained within-cluster variation

Below we see graphs that show how the fixed slope and random slope explain increasingly more within-company variation, reflected by the decreasing residuals.

Notice how the residual variances decrease of the 3 models depicted above decreases as we add salary_c fixed slope and random slope to the model (while random intercept variance does not decrease).

# VarCorr() returns random effect and residual standard deviations
#  wrapping VarCorr() in print() with comp="Variance" returns variances
# random intercept model
print(VarCorr(m_ri), comp="Variance") 

##  Groups     Name        Variance
##  company_id (Intercept) 116.083 
##  Residual                73.412

# fixed slope of salary_c
print(VarCorr(m_salary_c_fixed), comp="Variance")  

##  Groups     Name        Variance
##  company_id (Intercept) 116.951 
##  Residual                64.732

# fixed and random slope of salary_c
print(VarCorr(m_salary_c_random), comp="Variance")  

##  Groups     Name        Variance Cov    
##  company_id (Intercept) 117.4425        
##             salary_c      0.3329   1.216
##  Residual                59.3291

Between-cluster variance decomposition

\[\widehat{var}(\bar{y}_{.j}-\bar{y}_{..}) = f_2 + m \] Explained variance component

• \(f_2\): level-2 predictor fixed effects explain between-cluster differences

Residual variance component

• \(m\) (\(\tau_{00}^2\)) quantifies unexplained between-cluster variation

To demonstrate, we will regress stress on level-2 pct_open:

# level-2 pct_open explains between-company variation in mean stress
m_pct_open <- lmer(stress ~ pct_open + 
                     (1 | company_id),
                   data=empstress)

summary(m_pct_open)

## Linear mixed model fit by REML ['lmerMod']
## Formula: stress ~ pct_open + (1 | company_id)
##    Data: empstress
## 
## REML criterion at convergence: 1103.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2748 -0.6105 -0.0242  0.5497  3.2955 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  company_id (Intercept) 105.77   10.284  
##  Residual                73.41    8.568  
## Number of obs: 150, groups:  company_id, 15
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  30.3667    10.0800   3.013
## pct_open      1.4508     0.9616   1.509
## 
## Correlation of Fixed Effects:
##          (Intr)
## pct_open -0.962

Compare the random intercept variances (and residual variances) for the model with just random intercepts to the model where we add pct_open as a predictor.

# random intercept model
print(VarCorr(m_ri), comp="Variance")

##  Groups     Name        Variance
##  company_id (Intercept) 116.083 
##  Residual                73.412

# level-2 pct_open as predictor
print(VarCorr(m_pct_open), comp="Variance")

##  Groups     Name        Variance
##  company_id (Intercept) 105.770 
##  Residual                73.412

The residuals below represent unexplained between-company variation in mean stress (i.e., random intercept variation), and we see how pct_open explains away some of this unexplained variation.

Rights and Sterba (2019) \(R^2\) Notation

The new multilevel \(R^2\) measures can include different terms in the numerator and denominator, so notation is required to distinguish among them:

• Superscripts denote sources of explained variance in the numerator (\(f_1\), \(f_2\), \(v\), or \(m\))
• Subscripts denote whether the variance being explained in the denominator is total, within-cluster, or between-cluster (\(t\), \(w\), or \(b\))

\(R^{2(f_1)}_w\) is thus the proportion of within-cluster variance explained by level-1 fixed effects.

Proportion of total variance explained by fixed effects

Recall that the model-implied total variance of \(y_{ij}\) is:

\[\widehat{var}(y_{ij}) = f_1 + f_2 + v + m + \hat{\sigma}_\epsilon^2\]

Proportion of total variance explained by level-1 (within-cluster) fixed effects

\[R^{2(f_1)}_t = \frac{f_1}{f_1 + f_2 + v + m + \hat{\sigma}_\epsilon^2}\]

Proportion of total variance explained by level-2 (between-cluster) fixed effects

\[R^{2(f_2)}_t = \frac{f_2}{f_1 + f_2 + v + m + \hat{\sigma}_\epsilon^2}\]

Proportion of total variance explained by fixed effects

\[\begin{align}R^{2(f)}_t &= R^{2(f_1)}_t + R^{2(f_2)}_t \\ &= \frac{f_1+f_2}{f_1 + f_2 + v + m + \hat{\sigma}_\epsilon^2}\end{align}\]

Superscript \(f\) represents level-1 and level-2 fixed effects together.

Proportion of total variance explained by random variation

Proportion of total variance explained by random slope variation

\[R^{2(v)}_t = \frac{v}{f_1 + f_2 + v + m + \hat{\sigma}_\epsilon^2}\]

Proportion of total variance explained by random intercept variation

\[R^{2(m)}_t = \frac{m}{f_1 + f_2 + v + m + \hat{\sigma}_\epsilon^2}\]

Proportion of total variance explained by predictors and model

Proportion of total variance explained by predictors (fixed effects and random slopes)

\[\begin{align}R^{2(fv)}_t &= R^{2(f_1)}_t + R^{2(f_2)}_t + R^{2(v)}_t \\ &= \frac{f_1+f_2+v}{f_1 + f_2 + v + m + \hat{\sigma}_\epsilon^2}\end{align}\]

Proportion of total variance explained by model (all fixed and random effects)

\[\begin{align}R^{2(fvm)}_t &= R^{2(f_1)}_t + R^{2(f_2)}_t + R^{2(v)}_t + R^{2(m)}_t \\ &= \frac{f_1+f_2+v+m}{f_1 + f_2 + v + m + \hat{\sigma}_\epsilon^2}\end{align}\]

Proportion of within-cluster variance explained

Recall that the model-implied within cluster variance is:

\[\widehat{var}(y_{ij}-\bar{y}_{.j}) = f_1 + v + \hat{\sigma}_\epsilon^2\]

Proportion of within-cluster variance explained by level-1 fixed effects

\[R^{2(f_1)}_w = \frac{f_1}{f_1 + v + \hat{\sigma}_\epsilon^2}\]

Proportion of within-cluster variance explained by random slopes

\[R^{2(v)}_w = \frac{v}{f_1 + v + \hat{\sigma}_\epsilon^2}\]

Proportion of within-cluster variance explained by level-1 predictors (fixed effects and random slopes)

\[\begin{align}R^{2(f_1v)}_w &= R^{2(f_1)}_w + R^{2(v)}_w \\ &= \frac{f_1 + v}{f_1 + v + \hat{\sigma}_\epsilon^2}\end{align}\]

Proportion of between-cluster variance explained

Recall that the model-implied between-cluster variance is given by:

\[\widehat{var}(\bar{y}_{.j}-\bar{y}_{..}) = f_2 + m \]

Proportion of between-cluster variance explained by level-2 fixed effects

\[R^{2(f_2)}_b = \frac{f_2}{f_2 + m}\]

Proportion of between-cluster variance explained by random intercept (cluster mean) variation

\[R^{2(m)}_b = \frac{m}{f_2 + m}\]

We do not estimate \(R^{2(f_2m)}_b\) because it always equals 1.

The r2mlm package

r2mlm calculates these new \(R^2\) measures in R.

In this section we discuss the following r2mlm functions:

• r2mlm(): calculates and visualizes all possible \(R^2\) measures for a model
• r2mlm_ci(): calculates bootstrapped confidence intervals for \(R^2\) measures
• r2mlm_comp()compares 2 models and calculates differences in all \(R^2\) measures

The r2mlm package calculates \(R^2\) measures for multilevel models run using lmer() from the lme4 package^†.

# load packages
library(lme4)
library(r2mlm)

^†or using nlme() from the accompanying package nlme

Note: Much of the material in this section is adapted from Shaw et al. (2023), an article discussing the features and usage of the r2mlm package

Example data

For this section we use a simulated dataset related to teacher job satisfaction that loads with r2mlm called teachsat.

Sample is 300 schools (level-2), each with 30 teachers (level-1), for a total \(N=9000\).

Variables

• schoolID
• teacherID
• satisfaction: teacher job satisfaction (outcome), 1-10 scale
• control_c: school-mean-centered teacher’s rating of control over curriculum
• salary_c: school-mean-centered teacher’s salary (in thousands of dollars)
• control_m: school mean rating of teacher’s rating of control over curriculum
• salary_m: school mean teacher’s salary
• s_t_ratio: student-to-teacher ratio for the school

# teachsat multilevel data
head(teachsat)

##   schoolID teacherID satisfaction  control_c  salary_c control_m salary_m
## 1        1         1     8.582175  2.7239726  7.274769  5.474631 77.69538
## 2        1         2     7.428820 -0.6604034 13.321313  5.474631 77.69538
## 3        1         3     7.990401 -1.3002809 20.131688  5.474631 77.69538
## 4        1         4     5.100618 -0.6278236 -9.810052  5.474631 77.69538
## 5        1         5     7.256876  1.5190521  1.205659  5.474631 77.69538
## 6        1         6     5.154663 -1.9307584  1.778745  5.474631 77.69538
##   s_t_ratio
## 1        30
## 2        30
## 3        30
## 4        30
## 5        30
## 6        30

Random intercept model

We start with the random intercept model to quantify the amount of variation within and between clusters.

\[\begin{align} y_{ij} &= \beta_{0j} + \epsilon_{ij} \\ \beta_{0j} &= \gamma_{00} + u_{0j} \end{align}\]

# random intercept model
m_m <- lmer(satisfaction ~ 1 + (1|schoolID),
             data=teachsat)

summary(m_m)

## Linear mixed model fit by REML ['lmerMod']
## Formula: satisfaction ~ 1 + (1 | schoolID)
##    Data: teachsat
## 
## REML criterion at convergence: 29426.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8004 -0.6500 -0.0083  0.6370  3.4825 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  schoolID (Intercept) 0.8376   0.9152  
##  Residual             1.3955   1.1813  
## Number of obs: 9000, groups:  schoolID, 300
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  5.99852    0.05429   110.5

We calculate the ICC by dividing the random intercept variance by the sum of the random intercept variance and residual variance:

#ICC
.8376/(.8376+1.3955)

## [1] 0.375084

• 37.5% of the total variation in satisfaction is between-cluster (differences in mean satisfaction across schools)
• 62.5% is within-cluster (differences in satisfaction between teachers in the same school).

To calculate all relevant \(R^2\) using r2mlm(), simply supply the lmer object as an argument:

# R^2 for random intercept model
r2mlm(m_m)

## $Decompositions
##                     total within between
## fixed, within   0.0000000      0      NA
## fixed, between  0.0000000     NA       0
## slope variation 0.0000000      0      NA
## mean variation  0.3750827     NA       1
## sigma2          0.6249173      1      NA
## 
## $R2s
##         total within between
## f1  0.0000000      0      NA
## f2  0.0000000     NA       0
## v   0.0000000      0      NA
## m   0.3750827     NA       1
## f   0.0000000     NA      NA
## fv  0.0000000      0      NA
## fvm 0.3750827     NA      NA

The output is divided into 3 sections:

Graphs

The three bars, left to right, represent total, within-cluster, and between-cluster variation of outcome.

• The only sources of variation in satisfaction are random intercept variation (blue stripes) and residual variation (white)
• random intercept variation (\(m\)) explains a bit less than 40% (~38%) of total variation, and all between-cluster variation (since there are no level 2 fixed effects)

$Decompositions

Proportion of variance explained by \(f_1\), \(f_2\), \(v\), \(m\), and \(\sigma_{\epsilon}^2\), respectively, for each of total, within-cluster, and between-cluster variance.

• mean variation (random intercept variation) accounts for 37.5% of the total variance of satisfaction
• sigma2 (residual variation) accounts for the remaining 62.5% of total variance
• All of the within-school variation is residual variation
• All of the between-school variation is random intercept variation
• NA means this source (row) cannot explain the variation captured by the column (e.g. \(f2\) cannot explain within-school variance)

$R2s

Displays all \(R^2\) measures calculable for a model.

• \(R^{2(m)}_t \approx 0.375\), 37.5% of the total variance is explained by random intercept variation
• \(R^{2(fvm)}_t \approx 0.375\), 37.5% of the total variance is explained by the model
• \(R^{2(m)}_b = 1\), 100% of the between-school variance is explained by random intercept variation

Adding fixed effects

Now we add fixed slopes of level-1 predictor salary_c and level-2 predictor salary_m to the random intercept model:

\[\begin{align} y_{ij} &= \beta_{0j} + \beta_1salary\_c + \epsilon_{ij} \\ \beta_{0j} &= \gamma_{00} +\gamma_{01}salary\_m + u_{0j} \\ \beta_1 &= \gamma_{10} \end{align}\]

# fixed effects and random intercept model
m_f.m <- lmer(satisfaction ~ salary_c + salary_m + (1 | schoolID),
             data=teachsat)
summary(m_f.m)

## Linear mixed model fit by REML ['lmerMod']
## Formula: satisfaction ~ salary_c + salary_m + (1 | schoolID)
##    Data: teachsat
## 
## REML criterion at convergence: 26428.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.4422 -0.6710 -0.0137  0.6522  3.7630 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  schoolID (Intercept) 0.5984   0.7736  
##  Residual             0.9980   0.9990  
## Number of obs: 9000, groups:  schoolID, 300
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 3.104077   0.267652   11.60
## salary_c    0.071763   0.001219   58.88
## salary_m    0.039240   0.003575   10.98
## 
## Correlation of Fixed Effects:
##          (Intr) slry_c
## salary_c  0.000       
## salary_m -0.985  0.000

The salary_c coefficient suggests that teachers that have higher salaries than other teachers within the same school are more satisfied.

The salary_m coefficient suggests that schools with higher mean teacher salaries have higher mean satisfaction.

Both the \(\hat{\tau}_{00}^2=0.598\) and \(\hat{\sigma}_{\epsilon}^2=0.998\) have been reduced from the random intercept model (\(0.838\) and \(1.396\), respectively), suggesting a significant proportion of each of the variances have been explained by the fixed effects.

Let’s see the \(R^2\) measures for this model:

# R^2 for fixed effects and random intercept model
r2mlm(m_f.m)

## $Decompositions
##                     total    within   between
## fixed, within   0.1720308 0.2780736        NA
## fixed, between  0.1135477        NA 0.2977533
## slope variation 0.0000000 0.0000000        NA
## mean variation  0.2678004        NA 0.7022467
## sigma2          0.4466212 0.7219264        NA
## 
## $R2s
##         total    within   between
## f1  0.1720308 0.2780736        NA
## f2  0.1135477        NA 0.2977533
## v   0.0000000 0.0000000        NA
## m   0.2678004        NA 0.7022467
## f   0.2855784        NA        NA
## fv  0.2855784 0.2780736        NA
## fvm 0.5533788        NA        NA

Graphs

• The contributions of the level-1 fixed effect (solid red) and level-2 fixed effect (solid blue) are now visualized
• More than 50% of the total variance has been explained, and the 2 fixed effects combine to explain a bit less than 30%
• The fixed effect of salary_c (level-1) explains a little less than 30% of the within-school variance
• The fixed effect of salary_m (level-2) explains approximately 30% of the between-school variance

$Decompositions

• The level-2 fixed effect explains more of the total variance than the level-1 fixed effect
• mean variation and sigma2, level-2 and level-1 residual variances, respectively, make up less of the total variance compared to the random intercept model because the fixed effects are explaining some of the unexplained variance away

$R2s

• The model explains 55.3% of the total variance, reflected by \(R^{2(fvm)}_t \approx .553\)
• We can verify that the \(R^2\) measures with multiple sources of explained variance are sums of \(R^2\) with single terms, for example:

\[R^{2(f)}_t = R^{2(f_1)}_t + R^{2(f_2)}_t\]

and

\[R^{2(fvm)}_t = R^{2(f_1)}_t + R^{2(f_2)}_t + R^{2(v)}_t + R^{2(m)}_t\]

Adding a random slope

The final source of explained variation to examine with r2mlm() is \(v\), random slopes.

We allow the effect of salary_c to vary by school via random slopes, which may explain more within-school variation in teacher satisfaction.

\[\begin{align} y_{ij} &= \beta_{0j} + \beta_{1j}salary\_c + \epsilon_{ij} \\ \beta_{0j} &= \gamma_{00} +\gamma_{01}salary\_m + u_{0j} \\ \beta_{1j} &= \gamma_{10} + u_{1j} \end{align}\]

# random slope, fixed effects and random intercept model
m_v.f.m <- lmer(satisfaction ~ salary_c + salary_m + (1 + salary_c | schoolID),
             data=teachsat)
summary(m_v.f.m)

## Linear mixed model fit by REML ['lmerMod']
## Formula: satisfaction ~ salary_c + salary_m + (1 + salary_c | schoolID)
##    Data: teachsat
## 
## REML criterion at convergence: 25376.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3085 -0.6585 -0.0081  0.6432  3.6531 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  schoolID (Intercept) 0.604140 0.77726      
##           salary_c    0.002291 0.04786  0.04
##  Residual             0.826908 0.90934      
## Number of obs: 9000, groups:  schoolID, 300
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 3.086462   0.267511   11.54
## salary_c    0.071689   0.002987   24.00
## salary_m    0.039479   0.003573   11.05
## 
## Correlation of Fixed Effects:
##          (Intr) slry_c
## salary_c  0.006       
## salary_m -0.985  0.000

Although the salary_c slope variance appears small at 0.0023, the standard deviation is 0.0479 (\(\sqrt{0.0023} \approx 0.0479\)), more than 50% the size of the coefficient itself (0.0717).

\(\hat{\sigma}_{\epsilon}^2=0.827\) has again decreased compared to the model with just fixed effects (\(0.998\)), suggesting that \(v\) explains some of the total and within-school variance.

Now let’s see all possible \(R^2\) measures calculated for this model:

# model has all possible sources of explained variance
r2mlm(m_v.f.m)

## $Decompositions
##                      total    within   between
## fixed, within   0.17107144 0.2776837        NA
## fixed, between  0.11452685        NA 0.2982981
## slope variation 0.07624695 0.1237643        NA
## mean variation  0.26940742        NA 0.7017019
## sigma2          0.36874735 0.5985520        NA
## 
## $R2s
##          total    within   between
## f1  0.17107144 0.2776837        NA
## f2  0.11452685        NA 0.2982981
## v   0.07624695 0.1237643        NA
## m   0.26940742        NA 0.7017019
## f   0.28559828        NA        NA
## fv  0.36184523 0.4014480        NA
## fvm 0.63125265        NA        NA

Graphs

• The model now explains a bit more than 60% of the total variance
• salary_c slope variation (red stripes) explains a little less than 10% of the total variance
• salary_c slope variation explains a bit more than 10% of the within-cluster variance, which is now about 40% explained by the model

$Decompositions

• salary_c slope variation is the smallest source of explained variance
• No values are 0, meaning all sources of variation are incorporated into this model

$R2s

• \(R^{2(fvm)}_t \approx 0.631\), the model explains about 63.1% of the total variance
• \(R^{2(fv)}_t \approx 0.362\), the predictors together explain about 36.2% of the total variance (via fixed and random slopes)
• \(R^{2(fv)}_w \approx 0.401\), salary_c explains about 40.1% of the within-school variance (via its fixed and random slope)

Confidence intervals for multilevel \(R^2\) measures

The function r2mlm_ci() calculates bootstrapped confidence intervals for each \(R^2\)^†.

Using r2mlm_ci() requires the bootmlm package, which resamples multilevel data using the bootstrap.

The bootmlm package is still in development, so is not yet available for download from CRAN, but can be downloaded from github using the following code:

if (!require("remotes")) {
  install.packages("remotes")
}
remotes::install_github("marklhc/bootmlm")

The following arguments are required when calling r2mlm_ci()

• model: an lmer() or nlme() object
• nsim: number of bootstrap samples (e.g. 500 or 1000)
• boottype: "parametric", which assumes normally distributed residuals, or "residual", which does not
• confinttype: type of confidence interval, one of "norm" (normal), "basic", or "perc" (percentile)

Below, we request 95% percentile confidence intervals using 1000 residual bootstrap samples:

# set random seed to reproduce results
set.seed(1932)
# this takes a long time to run!
r2mlm_ci(m_v.f.m, nsim=1000, boottype="residual", confinttype = "perc",
         progress=F) # suppresses a progress bar from appearing in output

##                 2.5 %     97.5 %
## f1_total   0.14921313 0.19447585
## f2_total   0.08144853 0.15445686
## v_total    0.06288597 0.09266709
## m_total    0.23435219 0.30509080
## f_total    0.25205535 0.32365693
## fv_total   0.32817014 0.40108991
## fvm_total  0.60408394 0.65694501
## f1_within  0.24738847 0.30914729
## v_within   0.10441981 0.14750445
## fv_within  0.37114323 0.43443146
## f2_between 0.22192875 0.38347671
## m_between  0.61652329 0.77807125

^† As of writing this workshop, the r2mlm_ci() function appears to be still somewhat in development.

Comparing models

Although r2mlm() combines individual sources of explained variance of the same type (e.g., \(f1\) represents all level-1 fixed effects), we may be interested in the unique contribution of a single variable or smaller subset of variables.

Or, we may wish to compare \(R^2\) measures between multilevel models with completely different sets of predictors.

For both purposes we can use r2mlm_comp() to calculate the difference in all \(R^2\) measures between 2 models.

Example: calculating unique contributions to \(R^2\)

Assume we are interested in calculating how much various \(R^2\) measures increase when we add the following to our previous model with salary_c and salary_m as predictors:

• control_c fixed and random slope
• control_m fixed slope

First, we fit the model:

# model with additional predictors
m_all <-  lmer(satisfaction ~ salary_c + salary_m + control_c + control_m +
                 (1 + salary_c + control_c | schoolID),
             data=teachsat)
summary(m_all)

## Linear mixed model fit by REML ['lmerMod']
## Formula: satisfaction ~ salary_c + salary_m + control_c + control_m +  
##     (1 + salary_c + control_c | schoolID)
##    Data: teachsat
## 
## REML criterion at convergence: 22463.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3888 -0.6513 -0.0009  0.6350  3.9643 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr       
##  schoolID (Intercept) 0.609125 0.7805              
##           salary_c    0.002362 0.0486    0.07      
##           control_c   0.028816 0.1698    0.04 -0.13
##  Residual             0.553253 0.7438              
## Number of obs: 9000, groups:  schoolID, 300
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 2.692568   0.345176   7.801
## salary_c    0.072820   0.002960  24.604
## salary_m    0.040016   0.003557  11.251
## control_c   0.278388   0.011053  25.186
## control_m   0.072244   0.042491   1.700
## 
## Correlation of Fixed Effects:
##           (Intr) slry_c slry_m cntrl_c
## salary_c   0.008                      
## salary_m  -0.787  0.000               
## control_c  0.004 -0.107  0.000        
## control_m -0.637  0.000  0.044  0.000 
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00825855 (tol = 0.002, component 1)

• the large t value value for control_c suggests strong evidence that it explains a significant portion of within-school variance in satisfaction
• the small t value for control_m suggests a lack of strong evidence that it explains much between-school variance (even though its coefficient is larger than that for salary_m).

Nest we use r2mlm_comp to calculate the difference in \(R^2\) measures between the two models:

# calculate R^2 difference for individual contributions
r2mlm_comp(m_v.f.m, m_all)

## $`Model A R2s`
##          total    within   between
## f1  0.17107144 0.2776837        NA
## f2  0.11452685        NA 0.2982981
## v   0.07624695 0.1237643        NA
## m   0.26940742        NA 0.7017019
## f   0.28559828        NA        NA
## fv  0.36184523 0.4014480        NA
## fvm 0.63125265        NA        NA
## 
## $`Model B R2s`
##         total    within   between
## f1  0.2567172 0.4194249        NA
## f2  0.1180306        NA 0.3042572
## v   0.1102090 0.1800596        NA
## m   0.2698998        NA 0.6957428
## f   0.3747479        NA        NA
## fv  0.4849568 0.5994845        NA
## fvm 0.7548567        NA        NA
## 
## $`R2 differences, Model B - Model A`
##            total     within      between
## f1  0.0856457841 0.14174120           NA
## f2  0.0035037875         NA  0.005959122
## v   0.0339620354 0.05629528           NA
## m   0.0004924254         NA -0.005959122
## f   0.0891495715         NA           NA
## fv  0.1231116070 0.19803648           NA
## fvm 0.1236040324         NA           NA

The output contains:

• Visualizations of the variance decompositions for both models
• Visual comparisons between the models of the total, within, and between variance decompositions
  • control_c appears to have increased the contributions of both fixed slopes and slope variation to within-school variance
  • control_m appears to very slightly incresae the contribution of fixed slopes to between-school variance
• \(R^2\) measures for each model
  • the fixed slope of salary_c explains about 8.6% more of the total variance (than the reduced model), and about 14.2% more of the within-school variance
  • the random slope of salary_c explains about 3.4% more of the total variance and about 5.6% more of the within-school variance
  • the fixed slope of salary_m explains about 0.35% more of total variance, and about 0.6% of the between-school variance
  • together, salary_c and salary_m explain about 12.4% more of the total variance

Example: comparing different sets of predictors

Suppose we are interested in whether salary or curriculum control explains more of the variation in teacher satisfaction.

To assess this, we will compare 2 models:

• salary_c with fixed and random slopes, salary_m with fixed slope
• control_c with fixed and random slopes, control_m with fixed slope

If you compare 2 non-nested models with r2mlm_comp(), you must specify the dataset as well.

# salary model
m_salary <- lmer(satisfaction ~ salary_c + salary_m +
                   (1 + salary_c | schoolID),
                 data=teachsat)

# curriculum control model
m_control <- lmer(satisfaction ~ control_c + control_m +
                   (1 + control_c | schoolID),
                 data=teachsat)

# compare models
r2mlm_comp(m_salary, m_control, data=teachsat)

## $`Model A R2s`
##          total    within   between
## f1  0.17107144 0.2776837        NA
## f2  0.11452685        NA 0.2982981
## v   0.07624695 0.1237643        NA
## m   0.26940742        NA 0.7017019
## f   0.28559828        NA        NA
## fv  0.36184523 0.4014480        NA
## fvm 0.63125265        NA        NA
## 
## $`Model B R2s`
##          total     within     between
## f1  0.07973695 0.12847381          NA
## f2  0.00128136         NA 0.003377757
## v   0.03112197 0.05014435          NA
## m   0.37807113         NA 0.996622243
## f   0.08101831         NA          NA
## fv  0.11214028 0.17861816          NA
## fvm 0.49021141         NA          NA
## 
## $`R2 differences, Model B - Model A`
##           total      within    between
## f1  -0.09133449 -0.14920993         NA
## f2  -0.11324548          NA -0.2949203
## v   -0.04512498 -0.07361995         NA
## m    0.10866371          NA  0.2949203
## f   -0.20457997          NA         NA
## fv  -0.24970495 -0.22282988         NA
## fvm -0.14104124          NA         NA

• Model A is the salary model and Model B is the curriculum control model
• salary_c explains more within-school variation than control_c
• salary_m explains more between-school variation than salary_m
• Overall, salary explains about 14.1% more total variance than curriculum control

Using uncentered level-1 predictors

Despite general recommendations to group-mean center level-1 predictors, r2mlm() can be used with uncentered level-1 predictors.

Uncentered level-1 predictors potentially explain variation at both levels, thus:

• \(f_1\) and \(f_2\) cannot be calculated, so \(f\), the total variance explained by fixed effects is used instead
• only total variance is decomposed, because within-cluster variance and between-cluster variance are defined using \(f_1\) and \(f_2\), respectively
• random slope variation \(v\) can also explain variation at both levels
• the model implied total variance is now

\[\widehat{var}(y_{ij})_{x_{uncentered}} = f + v + m + \hat{\sigma}_\epsilon^2\] * Only \(R^{2(f)}_t\), \(R^{2(v)}_t\), \(R^{2(m)}_t\), \(R^{2(fv)}_t\), and \(R^{2(fvm)}_t\) will be reported (see Table 5 of Rights and Sterba(2019) for definitions)

Below, we sum control_c and control_m to create uncentered control, which we then model as a predictor:

# create uncentered level-1 control
teachsat$control <- teachsat$control_c + teachsat$control_m


# fixed and random slopes for control
m_control_un <- lmer(satisfaction ~ control + 
                              (1 + control | schoolID),
                            data=teachsat)

summary(m_control_un)

## Linear mixed model fit by REML ['lmerMod']
## Formula: satisfaction ~ control + (1 + control | schoolID)
##    Data: teachsat
## 
## REML criterion at convergence: 28017.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2656 -0.6153 -0.0006  0.6128  3.9539 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  schoolID (Intercept) 1.50889  1.2284        
##           control     0.02637  0.1624   -0.63
##  Residual             1.14230  1.0688        
## Number of obs: 9000, groups:  schoolID, 300
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  4.69931    0.08008   58.68
## control      0.26333    0.01178   22.34
## 
## Correlation of Fixed Effects:
##         (Intr)
## control -0.705

We see that we only get decompositions and \(R^2\) measures related to total variance for a model with uncentered level-1 predictors:

# r2mlm with uncentered level-1 predictor
r2mlm(m_control_un)

## $Decompositions
##                      total
## fixed           0.10497843
## slope variation 0.03992823
## mean variation  0.38130827
## sigma2          0.47378506
## 
## $R2s
##          total
## f   0.10497843
## v   0.03992823
## m   0.38130827
## fv  0.14490667
## fvm 0.52621494

• About 10.5% of the total variance of satisfaction is explained by the fixed slope of control
• About 4% of the total variance is explained by variation of the random slope of control
• After accounting for the fixed and random slopes of control, the random intercept variation explains about 38.1% of the total variance
• The model explains about 52.6% of the total variance of satisfaction

Which \(R^2\) to report?

Rights and Sterba (2019) explicitly caution against reporting \(R^2\) measures that combine several sources (e.g. \(R^{2(fvm)}_t\)) without reporting the \(R^2\) measures that capture a single source.

They provide an example where \(R^{2(fvm)}_t=0.8\) is reported alone, and thus it is unclear how much each of, \(f_1\), \(f_2\), \(v\) and \(m\) contribute.

The following figure illustrates the issue, where a model includes the fixed and random slope of a level-1 predictor and random intercepts:

Although \(R^{2(f,v,m)}_t=0.8\) for all models, the interpretations are very different:

• the left-most image depicts the level-1 fixed slope being the sole source of explained variance, suggesting a strong mean effect of the predictor but no random variation around this mean effect
• the second image from the left depicts level-1 random slopes being the sole source of explained variance, suggesting a zero mean (fixed) effect of the predictor
• the third image from the left depicts random intercept variation being the sole source of explained variance, with neither fixed nor random effects of the predictor
• the right-most image depicts equal contributions of level-1 fixed slope, random slope variation and random intercept variation

The interpretation of \(R^{2(f,v,m)}_t=0.8\) could be thus clarified by also reporting \(R^{2(f)}_t\), \(R^{2(v)}_t\), and \(R^{2(m)}_t\).

Excerise 2

1. Use r2mlm() on the model in (5) in Exercise 1 (fixed and random slopes of sleep_c, fixed slope for sleep_m) and interpret all output
2. Use r2mlm_comp() to determine whether (sleep_c + sleep_m) or (study_c + study_m) explain more variation in words. Model fixed and random slopes for the subject-mean-centered predictors.

References

Original theoretical article introducing new \(R^2\) measures

Rights, J. D., & Sterba, S. K. (2019). Quantifying explained variance in multilevel models: An integrative framework for defining R-squared measures. Psychological methods, 24(3), 309.

Article that reviews multilevel models and discusses usage of r2mlm in detail

Shaw, M., Rights, J. D., Sterba, S. S., & Flake, J. K. (2023). r2mlm: An R package calculating R-squared measures for multilevel models. Behavior Research Methods, 55(4), 1942-1964.

Extension: \(R^2\) for multilevel models with more than 2 levels

Rights, J. D., & Sterba, S. K. (2023). R-squared measures for multilevel models with three or more levels. Multivariate behavioral research, 58(2), 340-367.