OMNI ADRD: Multilevel Modeling

Modeling Continuous Outcome using the Basic Multilevel Model

Overview

This tutorial will cover multilevel modeling for cognitive aging research. Multilevel modeling is a statistical technique that accounts for non-independence of observations attributable to their “clustering.” For example, observations are clustered when a participant is observed at multiple time points. It is also sometimes referred to as a mixed-effects model or hierarchical linear model.

Outline

1. Background
2. Running a multilevel model in R
3. Conclusion

Background

In cognitive aging research, multilevel models are especially useful in contexts where participants have repeated measurements or the data have a natural “clustering” to the observations (e.g., cluster randomized controlled trials). In this context, the goal of the multilevel model is to partition the variance between-clusters and within-clusters. In repeated measures designs, this is called within-person and between-person variability.

The two can be defined as follows:

• Within-person variability explains how people differ from themselves across timepoints.
• Between-person variability explains how people differ from each other.

While the difference between the two is easy to define, they are easily confused in practice, and this confusion can lead to misleading substantive conclusions. Take the relationship between aging and processing speed. Suppose we measure each person’s processing speed once along with age. What we find is that being older correlates with having a lower processing speed. However, it is important to keep in mind that this is a claim about how people differ from each other rather than how people change over time. Because the model compares people against different people—as opposed to tracking one individual’s change over time—it is unable to say anything about the developmental mechanisms at play. A quantity measured once cannot vary within person, so it cannot be used to explain variation across time. Evaluating the claim requires measuring both speed and age repeatedly in the same people, so that a person’s own slowing can be related to their own aging.

The between and within-person relationship between processing speed and aging hsa been thoroughly investigated in the literature. One study (Sliwinski and Buschke, 1999) obtained repeated assessments from 302 older adults on processing speed using the Digit Symbol Substitution Test, and used multilevel modeling to partition the within-person and between-person variance. When looking at the relationship between processing speed and age, they found differing conclusions depending on the level of analysis. Between-persons, the relationship between age and processing speed is significant and negative; within-persons, the relationship was not significant.

The figure below shows the general shape of the problem in simulated data scaled to resemble the measures used in that study. Each cluster of points is one person measured repeatedly. The x-axis is age in years; the y-axis is the Digit Symbol Substitution test score (number of correct symbols).

When the clustering is ignored, the fitted line is dominated by the between-person differences: it is steep because older cohorts tend to be a lower-scoring demographic. When the clustering structure is accounted for, the within-person relationship becomes null. In extreme cases, the between-person and within-person associations can even differ in sign, a pattern often referred to as Simpson’s paradox.

If we ran a linear regression on these data without accounting for the clustering structure, we would find the following.

## 
## Call:
## lm(formula = dsst ~ age, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.9354 -2.4479 -0.1784  2.3084  6.7791 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  78.5718     8.6259   9.109 4.28e-11 ***
## age          -0.5173     0.1133  -4.567 5.09e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.438 on 38 degrees of freedom
## Multiple R-squared:  0.3543, Adjusted R-squared:  0.3373 
## F-statistic: 20.85 on 1 and 38 DF,  p-value: 5.088e-05

The slope of age is significant and negative: -0.52 DSST points per year of age. But the average within-person slope across the eight simulated people is only 0.03 and not statistically significant Reading the pooled coefficient as though it described within-person change would overstate how much of a person’s own cognitive decline is attributable to that person’s own aging. Multilevel modeling is one way to keep the two sources of variation distinct, and thus to avoid drawing spurious conclusions from their conflation.

How to specify the multilevel model

Our tutorial will be modeled after a real example from the literature, Hyun et al. (2019). They used ecological momentary assessment data to assess the effects of stress anticipation on working memory. However, the data for this tutorial are simulated, and thus should only be used for pedagogical purposes. Because Hyun et al. administered similar cognitive tasks in a short time period, they must account for each participant’s starting proficiency and their person-specific practice effects. That is, participants may significantly vary in their beginning proficiency with the task and the rate at which they improve. One way to handle these person-specific effects is through multilevel modeling.

In multilevel modeling, the model is split into two levels. The level-1 model explains within-person variability while the level-2 model explains between-person variability. In this example, level-1 corresponds to the repeated working memory tasks while level-2 would explain the characteristics of the person measured; however, this generalizes to any multilevel data structure, where level-2 corresponds to cluster level traits (between-person variables) and level-1 corresponds to the observations within each cluster (within-person variables).

In order to specify a multilevel model, you need a cluster-ID variable. In our case, this will be the participant ID. Once you have the indicator, you can specify the within-person and between-person models.

The within-person model explains the measure-to-measure variability when only looking at one person’s measurements, or more generally, the variability of the observations within the cluster. In Hyun et al. (2019), these within-person observations are the repeated EMA sessions completed by each participant over the course of the study. In notation, the within-person model is \[\text{Task Score}_{it} = \beta_{0i} + \beta_{1}\text{Stress Anticipation}_{it} + \beta_{2i}\text{Study Day}_{it} + e_{it}\] where the subscript \(i\) denotes the \(i\)th participant and the subscript \(t\) denotes the \(t\)th EMA. The outcome, \(\text{Task Score}_{it}\), is participant \(i\)’s working memory task score at time \(t\). In a more generalized form, we use \(i\) to denote the level-2 units and \(t\) for the level-1 observations.

Notice that some of the coefficients have a \(i\) subscript. This means that they are allowed to vary across people. In other words, each effect is allowed to be “person-specific.” Accordingly, \(\beta_{0i}\) represents a person-specific intercept and \(\beta_{2i}\) represents a person-specific practice effect (operationalized by the number of days in the study). Notice that \(\beta_1\) does not have the \(i\) subscript. This is because the effect of stress anticipation on task score is not modeled to vary across people, in other words, it is not person-specific.

There is an additional error term, \(e_{it}\), which is assumed to be normally distributed \[e_{it} \sim N(0,\mathbf{R}),\] \[\mathbf{R} = \mathbf{I}\sigma^{2}_{e},\] where \(\mathbf{I}\) is the identity matrix (diagonal matrix of 1s) and \(\sigma^{2}_{e}\) is the residual (within-person) variance.

The level-2 model explaining the level-1 model’s coefficients is given by \[\beta_{0i} = \gamma_{00} + \gamma_{01}\text{Age}_{i} + u_{0i},\] \[\beta_{1} = \gamma_{10},\] \[\beta_{2i} = \gamma_{20} + u_{2i}.\]

Notice that in the first equation for \(\beta_{0i}\) there is a level-2 predictor, \(\text{Age}_i\). One of the advantages of multilevel modeling is its ability to incorporate predictors at both level-1 and level-2. By modeling \(\text{Age}_i\) as a level-2 predictor, one can capture the effects of between-person variables on the outcome.

The \(\gamma\) terms represent fixed effects. These are the parameters governing each of the coefficients in the level-1 within-person model. Specifically, \(\gamma_{00}\) is the intercept for the intercept term \(\beta_{0i}\), \(\gamma_{10}\) is the intercept for \(\beta_1\), and \(\gamma_{20}\) is the intercept for \(\beta_{2i}\). \(\gamma_{01}\) represents the fixed effects of Age on the intercept. Notice that only \(\beta_{0i}\) and \(\beta_{2i}\) have the additional \(u\) term. This is because they are allowed to vary across people, while \(\beta_1\) is fixed. \(u_{0i}\) and \(u_{2i}\) represent the random effects for the intercept and slope, respectively. The random effects are the person-specific deviation from their expected intercept or slope. In multilevel models the random effects are assumed to be multivariate normal, \[\mathbf{U}_{i} \sim N(0,\mathbf{G})\] \[\mathbf{G} = \left[\begin{array} {rr} \sigma^{2}_{u0} & \sigma_{u0u2} \\ \sigma_{u2u0} & \sigma^{2}_{u2} \end{array}\right]\]

Prudent choice of time metric is warranted in cognitive aging research. In our motivating example, \(\text{Study Day}_{it}\) is time in study (days since enrollment), which measures task familiarity in order to capture retest or practice effects. Chronological age is another option to model time at level-1 and is most useful when indexing to age-related developmental change. Baseline age can also be incorporated as a level-2 predictor of intercept to model cross-sectional differences between age groups (as we do above), though these differences also include cohort-effects. Because the three metrics answer different questions, the choice should follow from the research question. One caveat is that because \(\text{Chronological Age}_{it} = \text{Baseline Age}_i + \text{Study Day}_{it}\), the model can only include two of the three predictors—including all three is not statistically identified

The level-1 and level-2 models can be combined, forming the “combined model.” The combined model is the specification that most software packages use (including the lme4 package in R, which we will later introduce). The combined model is obtained by substituting the level-2 equations for their corresponding coefficient. For a simple example, the level-1 and level-2 equations \[Y_{it} = \beta_{0i} + e_{it}\] \[\beta_{0i} = \gamma_{00} + u_{0i} \] become \[Y_{it} = \gamma_{00} + u_{0i} + e_{it}\] in the combined model. This is achieved by substituting \(\gamma_{00} + u_{0i}\) for \(\beta_{0i}\) in the level-1 equation.

For a more complex example, the combined equation for our motivating example would be \[\text{Task Score}_{it} = (\gamma_{00} + \gamma_{01}\text{Age}_{i} + u_{0i}) + (\gamma_{10})\text{Stress Anticipation}_{it} + (\gamma_{20} + u_{2i})\text{Study Day}_{it} + e_{it}\] \[\text{Task Score}_{it} = \gamma_{00} + \gamma_{01}\text{Age}_{i} + \gamma_{10}\text{Stress Anticipation}_{it} + \gamma_{20}\text{Study Day}_{it} + u_{0i} + u_{2i}\text{Study Day}_{it} + e_{it}.\] It is convention to order the combined equation such that the fixed components come first (the \(\gamma\) terms), followed by the random components (the \(u\) and \(e\) terms).

B. Running a multilevel model in R

To specify the multilevel model in R, we will use the lme4 and lmerTest packages.

Preliminaries

Loading Libraries

Loading libraries used in this script.

library(dplyr)      # for data manipulation
library(ggplot2)    # for data visualization
library(lme4)       # for multilevel models
library(lmerTest)   # for p-values

Loading Data

First, we load in the simulated dataset

dataset <- readRDS("OMNI_MLM_data.rds")

head(dataset)

##   participant_id study_day stress_anticipation      age test_error
## 1              1         1          -0.1510197 50.55294   3.278502
## 2              1         2          -0.1053900 50.55294   2.613113
## 3              1         3           0.7594594 50.55294   2.687086
## 4              1         4          -1.2426694 50.55294   2.593084
## 5              1         5          -0.5898462 50.55294   3.297520
## 6              1         6          -2.0391118 50.55294   2.587156

Set-up for basic multilevel model with a continuous outcome

Before running any formal model, it is best to inspect the data. First check out the distribution of the outcome variable (test_error) to confirm that it is continuous:

#pdf("TestErrorHistogram.pdf",height=3, width=3.5)
ggplot(data=dataset, aes(x=test_error)) +
  geom_histogram(aes(y=after_stat(density)),
                   binwidth=.25,
                   colour="black", fill="white",
                 na.rm = TRUE) +
  labs(x = "Test Error", y="Density") +
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=16))

#dev.off()

Now look at a subset of persons to see what each participant’s trajectory of test_error across study_day looks like.

#faceted plot
ggplot(data=dataset[which(dataset$participant_id <= 25),], aes(x=study_day, y=test_error)) +
  geom_point(na.rm=TRUE) +
  stat_smooth(method="lm", fullrange=TRUE, na.rm=TRUE) +
  xlab("Study Day") + ylab("Test Error") +
  facet_wrap( ~ participant_id) +
  theme(axis.title=element_text(size=16),
        axis.text=element_text(size=14),
        strip.text=element_text(size=14))

## `geom_smooth()` using formula = 'y ~ x'

From the plot, we can see that there appears to be substantial variation in the intercepts of each participant’s trajectory, and that the slopes across study day also differ from person to person. Some participants improve steadily over the course of the study while others change very little, which suggests person-specific practice effects. This is the pattern that motivates giving study_day both a fixed effect and a random effect in the models below.

Running the models

Centering

Before fitting the model, it is worth deciding how each predictor should be centered. Two conventions are common, within-person centering and grand-mean centering.

Within-person centering subtracts each person’s person-level mean from their observations so that the predictor expresses how far a given predictor departs from that person’s typical level. It removes the between-person information from the predictor entirely, leaving only within-person variation. To get the between-person effect back, you must include the person-level means as a level-2 predictor of the intercept.

Grand-mean centering subtracts the mean of all observations from the level-1 predictors. This shifts the variable so that zero is the total sample average. The coefficient of the grand-mean centered predictor does not cleanly separate between- and within-person variation, and is often called the “total effect.” To get a clean estimate of the within effects, you must include the person-level mean as a level-2 predictor of the intercept. In contrast to within-person centering, the person level mean at level-2 represents the contextual effect, or the difference between the between and within effects.

Here, we choose to within-person center stress_anticipation, which is entered as a person-level mean (suffix _pm) and a within-person-centered version (suffix _c).

dataset <- dataset %>%
  group_by(participant_id) %>%
  mutate(stress_anticipation_pm = mean(stress_anticipation, na.rm = TRUE),
         stress_anticipation_c  = stress_anticipation - stress_anticipation_pm) %>%
  ungroup()

To fit the multilevel model in R, we use the lmer() function in the lme4 package. The lmer() function uses the combined model, so the first step is to obtain the combined model as shown in the previous section. Separate the terms in the combined model into fixed components and random components. In the lmer() function, the fixed effects are specified by adding the variable name (or a 1 for intercepts), the same as the lm() function in base R. Random effects are specified in parentheses. In the parentheses, begin by listing the variables with random effects variances, then insert a |, and then insert your level-2 id. For example, (random effects variables | id).

To start, let’s fit the unconditional means model,

#unconditional means model
model0_fit <- lmer(formula = test_error ~ 1 + (1|participant_id),
              data=dataset,
              na.action=na.exclude)
summary(model0_fit)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: test_error ~ 1 + (1 | participant_id)
##    Data: dataset
## 
## REML criterion at convergence: 2066.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.83682 -0.62848 -0.00221  0.64149  3.02126 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  participant_id (Intercept) 0.5166   0.7187  
##  Residual                   0.1228   0.3504  
## Number of obs: 1680, groups:  participant_id, 240
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   2.25044    0.04718 238.99999    47.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Let’s extract the random effects variances with the VarCorr() function:

VarCorr(model0_fit)

##  Groups         Name        Std.Dev.
##  participant_id (Intercept) 0.71874 
##  Residual                   0.35036

We can then compute the intra-class correlation (ICC) as the ratio of the random intercept variance (between-person) to the total variance (between + within). The ICC gives the proportion of total variance attributable to the clustering structure. Large ICCs are generally indicative of the appropriateness of a multilevel model, while smaller ICCs could indicate that a multilevel model is unnecessary (although, a low ICC does not necessarily mean a multilevel model is theoretically or statistically inappropriate).

Store the random effect variances, which will be the first column of the VarCorr object (see above).

RandomEffects <- as.data.frame(VarCorr(model0_fit))
RandomEffects

##              grp        var1 var2      vcov     sdcor
## 1 participant_id (Intercept) <NA> 0.5165941 0.7187448
## 2       Residual        <NA> <NA> 0.1227518 0.3503595

Next, compute the ICC. It is the ratio of the random intercept variance (between-person var) over the total variance (between + within var):

ICC_between <- RandomEffects[1,4]/(RandomEffects[1,4]+RandomEffects[2,4]) 
ICC_between

## [1] 0.8080041

From the unconditional means model, the ICC indicated that among the total variance in test error, 80.8% is attributable to between-person sources of variation whereas 19.2% is attributable to within-person sources of variation.

OK - let’s add the predictors stress_anticipation and study_day to examine their effects on test error. Per the centering discussion above, stress_anticipation is entered as two separate terms—the within-person centered version (stress_anticipation_c) and the person-level mean (stress_anticipation_pm)—which lets the model estimate a within-person effect and a between-person effect separately. We do not want stress_anticipation_c to have a random effects variance, so we do not put it in the parentheses; however, we do want study_day to have a random effects variance, so we include it in the parentheses along with the intercept term, 1.

# fit model
model1_fit <- lmer(formula = test_error ~ 1 + study_day +
                      stress_anticipation_c + stress_anticipation_pm +
                      (1 + study_day|participant_id),
                    data=dataset,
                    na.action=na.exclude)
summary(model1_fit)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## test_error ~ 1 + study_day + stress_anticipation_c + stress_anticipation_pm +  
##     (1 + study_day | participant_id)
##    Data: dataset
## 
## REML criterion at convergence: 1765.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1676 -0.5746 -0.0119  0.6027  3.2743 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev. Corr  
##  participant_id (Intercept) 0.433741 0.65859        
##                 study_day   0.007912 0.08895  -0.08 
##  Residual                   0.080838 0.28432        
## Number of obs: 1680, groups:  participant_id, 240
## 
## Fixed effects:
##                          Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)             2.318e+00  4.526e-02  2.381e+02  51.216   <2e-16 ***
## study_day              -1.689e-02  6.709e-03  2.391e+02  -2.517   0.0125 *  
## stress_anticipation_c   6.699e-02  7.980e-03  1.331e+03   8.395   <2e-16 ***
## stress_anticipation_pm -1.013e-01  1.160e-01  2.380e+02  -0.874   0.3830    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                (Intr) stdy_d strss_ntcptn_c
## study_day      -0.223                      
## strss_ntcptn_c  0.011 -0.019               
## strss_ntcptn_p  0.004  0.000 -0.007

# save predicted scores
dataset$pred_m1 <- predict(model1_fit)

Interpreting the Results

Because stress_anticipation is entered as a within-person centered (_c) and person-mean (_pm) pair, it yields two coefficients: a within-person effect and a between-person effect.

• Fixed Effects:
  • (Intercept): The average value for test_error at day 0 and a person’s own average stress_anticipation
  • study_day: test_error changed over the course of the study. For every unit increase in study_day, test_error changed by -0.017 (p = 0.012)
  • stress_anticipation_c (within-person): The within-person effect was 0.067, which was statistically significantly different from zero (p < .001). On sessions when individuals anticipated higher stress than their own usual level, their test_error tended to be higher than usual.
  • stress_anticipation_pm (between-person): The between-person effect was -0.101, which was statistically not significantly different from zero (p = 0.383).
• Random Effects:
  • sd((Intercept)): There was a substantial extent of between-person differences in the average value of test_error
  • sd(study_day): There was a substantial extent of between-person differences in the average within-person practice effect (study_day slope)
  • cov((Intercept),study_day): The correlation between the random intercept and random slope indicates whether those who started with higher test_error also tended to improve more or less quickly over the study.

We can also get confidence intervals for the fixed and random effects. Depending on model complexity, the confint() function can sometimes take a while to run

# Get confidence intervals for both fixed and random effects
confint(model1_fit, signames = FALSE) #use signames = FALSE for interpretable variable names

## Computing profile confidence intervals ...

##                                                2.5 %       97.5 %
## sd_(Intercept)|participant_id             0.59318146  0.727045049
## sd_study_day|participant_id               0.07832306  0.100277850
## cor_study_day.(Intercept)|participant_id -0.23206917  0.080095953
## sigma                                     0.27320556  0.295969761
## (Intercept)                               2.22914026  2.406522773
## study_day                                -0.03006143 -0.003714111
## stress_anticipation_c                     0.05132896  0.082624288
## stress_anticipation_pm                   -0.32859946  0.125919576

The labels for the random effect variances are given by the “sd_” and “cor_” prefixes for standard deviations and correlations of the corresponding variable.

Adding a Person-level Predictor

We now add age to the model to examine whether between-person differences in age are related to average test error. Adding a level-2 predictor of the intercept is done by adding the variable as a fixed effect, seen below.

You can center age to a substantively meaningful value to improve interpretability of the intercept; in this example, we center at 40.

dataset <- dataset %>%
  mutate(age_c40 = age - 40)

# fit model
model2_fit <- lmer(formula = test_error ~ 1 + study_day +
                      stress_anticipation_c + stress_anticipation_pm + age_c40 +
                      (1 + study_day|participant_id),
                    data=dataset,
                    na.action=na.exclude)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.00472808 (tol = 0.002, component 1)
##   See ?lme4::convergence and ?lme4::troubleshooting.

summary(model2_fit)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## test_error ~ 1 + study_day + stress_anticipation_c + stress_anticipation_pm +  
##     age_c40 + (1 + study_day | participant_id)
##    Data: dataset
## 
## REML criterion at convergence: 1737.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1699 -0.5841 -0.0147  0.6018  3.3283 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev. Corr  
##  participant_id (Intercept) 0.367347 0.60609        
##                 study_day   0.007913 0.08895  -0.09 
##  Residual                   0.080838 0.28432        
## Number of obs: 1680, groups:  participant_id, 240
## 
## Fixed effects:
##                          Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)             2.199e+00  4.607e-02  2.411e+02  47.720  < 2e-16 ***
## study_day              -1.689e-02  6.709e-03  2.391e+02  -2.517   0.0125 *  
## stress_anticipation_c   6.711e-02  7.975e-03  1.333e+03   8.415  < 2e-16 ***
## stress_anticipation_pm -1.034e-01  1.074e-01  2.371e+02  -0.963   0.3365    
## age_c40                 2.282e-02  3.586e-03  2.370e+02   6.362 1.02e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                (Intr) stdy_d strss_ntcptn_c strss_ntcptn_p
## study_day      -0.220                                     
## strss_ntcptn_c  0.010 -0.019                              
## strss_ntcptn_p  0.005  0.000 -0.007                       
## age_c40        -0.407  0.000  0.004         -0.003        
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00472808 (tol = 0.002, component 1)
##   See ?lme4::convergence and ?lme4::troubleshooting.

# save predicted scores
dataset$pred_m2 <- predict(model2_fit)

• Fixed Effects:
  • (Intercept): The average value for test_error at day 0 for a 40-year-old person at their own average stress_anticipation
  • study_day: For every unit increase in study_day, test_error is expected to change by -0.017 (p = 0.012)
  • stress_anticipation_c (within-person): The within-person effect was 0.067, which was statistically significantly different from zero (p < .001). On sessions when individuals anticipated higher stress than their own usual level, their test_error tended to be higher.
  • stress_anticipation_pm (between-person): The between-person effect was -0.103, which was statistically not significantly different from zero (p = 0.336).
  • age_c40: The effect of age was 0.023 (p < .001). For every year a person’s age was above 40, their test_error tended to be higher.
• Random Effects:
  • sd((Intercept)): There was a substantial extent of between-person differences in the average value of test_error
  • sd(study_day): There was a substantial extent of between-person differences in the average within-person practice effect (study_day slope)
  • cov((Intercept),study_day): The correlation between the random intercept and random slope indicates whether those who started with higher test_error also tended to improve more or less quickly over the study.

If we wanted to include a level-2 predictor of a slope—for instance, if we wanted to include age_c40 as a predictor of the effect of study_day—we would include a multiplicative term. The following code includes age_c40 as a predictor of the study_day slope in addition to the intercept.

model3_fit <- lmer(formula = test_error ~ 1 + study_day + age_c40 + study_day*age_c40 +
                      stress_anticipation_c + stress_anticipation_pm +
                      (1 + study_day|participant_id),
                    data=dataset,
                    na.action=na.exclude)
summary(model3_fit)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: test_error ~ 1 + study_day + age_c40 + study_day * age_c40 +  
##     stress_anticipation_c + stress_anticipation_pm + (1 + study_day |  
##     participant_id)
##    Data: dataset
## 
## REML criterion at convergence: 1750.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1701 -0.5839 -0.0142  0.6018  3.3284 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev. Corr  
##  participant_id (Intercept) 0.367520 0.60623        
##                 study_day   0.007958 0.08921  -0.09 
##  Residual                   0.080837 0.28432        
## Number of obs: 1680, groups:  participant_id, 240
## 
## Fixed effects:
##                          Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)             2.199e+00  4.631e-02  2.371e+02  47.473  < 2e-16 ***
## study_day              -1.696e-02  7.397e-03  2.380e+02  -2.292   0.0228 *  
## age_c40                 2.280e-02  3.695e-03  2.371e+02   6.170 2.93e-09 ***
## stress_anticipation_c   6.711e-02  7.976e-03  1.333e+03   8.414  < 2e-16 ***
## stress_anticipation_pm -1.034e-01  1.074e-01  2.370e+02  -0.963   0.3365    
## study_day:age_c40       1.255e-05  5.902e-04  2.380e+02   0.021   0.9831    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                (Intr) stdy_d ag_c40 strss_ntcptn_c strss_ntcptn_p
## study_day      -0.241                                            
## age_c40        -0.417  0.100                                     
## strss_ntcptn_c  0.009 -0.014  0.006                              
## strss_ntcptn_p  0.005  0.000 -0.003 -0.007                       
## stdy_dy:_40     0.100 -0.417 -0.241 -0.009          0.000

# save predicted scores
dataset$pred_m3 <- predict(model3_fit)

• Fixed Effects:
  • (Intercept): The average value for test_error at day 0 for a 40-year-old person at their own average stress_anticipation
  • study_day: When age_c40 = 0 (i.e., at age 40), every unit increase in study_day is expected to change test_error by -0.017 (p = 0.023).
  • age_c40: The effect of age was 0.023 (p < .001). At day 0, for every year a person’s age was above 40, their test_error tended to be higher, holding stress_anticipation_c and stress_anticipation_pm constant.
  • study_day:age_c40: The interaction between study_day and age_c40 indicates whether the practice effect (i.e., the within-person change in test_error across days) differs as a function of age. A coefficient of 0 (p = 0.983) suggests that for each one-year increase in age above 40, the study_day slope shifts by that amount. This effect was statistically non-significant.
  • stress_anticipation_c (within-person): The within-person effect was 0.067, which was statistically significantly different from zero (p < .001). On sessions when individuals anticipated higher stress than their own usual level, their test_error tended to be higher.
  • stress_anticipation_pm (between-person): The between-person effect was -0.103, which was statistically not significantly different from zero (p = 0.336).
• Random Effects:
  • sd((Intercept)): There was a substantial extent of between-person differences in the average value of test_error (after accounting for age_c40, study_day, their interaction, and stress_anticipation_c).
  • sd(study_day): There was a substantial extent of between-person differences in the within-person practice effect (study_day slope), even after age_c40 was included as a predictor of that slope. This remaining variance reflects individual differences in the rate of change that are not explained by age_c40 alone.
  • cov((Intercept),study_day): The correlation between the random intercept and random slope indicates whether those who started with higher test_error also tended to improve more quickly over the study, above and beyond what is explained by age_c40.

Conclusion

This tutorial illustrated application of the basic multilevel model for continuously distributed outcomes.