In a previous post, we discussed the multilevel model for change, which investigates average change in time while accounting for individual-level variation. In this blog post, we will discuss how to extend this model to include predictors.

Time-varying versus time-constant variables

When working with longitudinal data, we can consider variables to be time-varying or time-constant. For example, a survey measure of general satisfaction is a time-varying variable as it is collected in each wave. So, for each individual, this could change (although it does not have to change). The year of birth is a time-constant variable. Collecting this information once from each individual is enough, as we do not expect it to change.

There are also some situations where we can treat a variable as time constant or time-varying. For example, we could recode birth year as age, a time-varying variable. Similarly, we could take each individual’s average of general satisfaction and make that variable time constant. In such situations, treating a variable as time-constant or variable should be based on theoretical considerations, our data and the modelling strategy used. You can read more about this in the post about data structures.

Explaining Change

Here, we will expand the multilevel model to include time-constant predictors. In a previous post, we explored how best to describe the change in time of the outcome by including non-linear effects. This is an essential first step in analyzing longitudinal data. Nevertheless, we often want to understand the causes of that change or at least how different groups may have different rates of change. We will turn to this aspect, discussing the inclusion of time-constant predictors and the same Understanding Society data.

Including Time Constant Predictors

We will start by adding a time-constant predictor to our model. For example, let’s explore how the change in mental health is different for men and women (where we treat these as time constant). If we add the “gndr” variable to our model, the new formula for our regression will be:

Y_ij = γ₀₀ + γ₀₁GNDR_i + γ₁₀TIME_ij + ξ_0i + ξ_1iTIME_ij + ϵ_ij

The interpretation of the coefficients is mostly the same as before. The interpretation of the intercept (γ00) changes and now refers to the expected value of the outcome when “gndr” and “time” are 0. In our case, this would refer to males at the beginning of the study. The effect of “gndr” (γ10) refers to how different females are from males in their mental health. The interpretation of the random effects is the same as before.

We will use the linear change in time model to keep things simple, but this can be easily expanded to include non-linear change.

Let’s run this model:

library(tidyverse)
library(lme4)

load("./data/us_clean_syn.RData")

usl <- usl |>
  mutate(
    wave0 = wave - 1,
    gndr = gndr.fct,
    degree = degree.fct
  )

m3 <- lmer(
  data = usl,
  sf12mcs ~ 1 + wave0 + gndr + (1 + wave0 | pidp)
)

summary(m3)

## Linear mixed model fit by REML ['lmerMod']
## Formula: sf12mcs ~ 1 + wave0 + gndr + (1 + wave0 | pidp)
##    Data: usl
## 
## REML criterion at convergence: 935141.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.3978 -0.4295  0.1569  0.5623  3.1354 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  pidp     (Intercept) 39.181   6.259         
##           wave0        3.017   1.737    -0.37
##  Residual             59.118   7.689         
## Number of obs: 127989, groups:  pidp, 48630
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 51.05364    0.05852  872.47
## wave0       -0.41879    0.02172  -19.28
## gndrFemale  -1.37299    0.07123  -19.27
## 
## Correlation of Fixed Effects:
##            (Intr) wave0 
## wave0      -0.413       
## gndrFemale -0.681  0.000
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00462369 (tol = 0.002, component 1)

We will concentrate on the fixed part of the model. The intercept, 51.054, is the expected mental-health score for men at the first wave. Mental health decreases by 0.419 points per wave. After accounting for time, women score 1.373 points lower than men at the first wave.

This model can include multiple time constant predictors of different nature, categorical or continuous. The interpretation will be the same as in a regular regression.

Allowing for different rates of change

One crucial assumption we made so far is that the rate of change is the same for men and females. That means that the difference in mental state remains constant in time. This is often an important question that we want to explore. From a substantive point of view, having different rates of change, which would lead to convergence or divergence, is very important. To explicitly examine this in our model, we can add an interaction between gender and time. This would allow for different rates of change for men and women. We would write this model as follows:

Y_ij = γ₀₀ + γ₀₁GNDR_i + γ₁₀TIME_ij + γ₁₁GNDR_i * TIME_ij + ξ_0i + ξ_1iTIME_ij + ϵ_ij

We can easily include interactions in lmer() by adding the two variables separated by :. So our new model is:

m4 <- lmer(
  data = usl,
  sf12mcs ~ 1 + wave0 + gndr + gndr:wave0 + (1 + wave0 | pidp)
)

summary(m4)

## Linear mixed model fit by REML ['lmerMod']
## Formula: sf12mcs ~ 1 + wave0 + gndr + gndr:wave0 + (1 + wave0 | pidp)
##    Data: usl
## 
## REML criterion at convergence: 935117.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.3894 -0.4295  0.1573  0.5625  3.1512 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  pidp     (Intercept) 39.170   6.259         
##           wave0        3.006   1.734    -0.37
##  Residual             59.116   7.689         
## Number of obs: 127989, groups:  pidp, 48630
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)      51.20087    0.06459 792.713
## wave0            -0.55125    0.03281 -16.802
## gndrFemale       -1.63514    0.08628 -18.952
## wave0:gndrFemale  0.23567    0.04376   5.385
## 
## Correlation of Fixed Effects:
##             (Intr) wave0  gndrFm
## wave0       -0.565              
## gndrFemale  -0.749  0.423       
## wv0:gndrFml  0.423 -0.750 -0.564
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00399381 (tol = 0.002, component 1)

With the interaction included, the intercept remains the expected first-wave score for men (51.201). Men decline by 0.551 points per wave, and women start 1.635 points lower. The positive interaction of 0.236 means that women decline more slowly. Their estimated slope is about -0.316 (-0.551 + 0.236), so the gender gap narrows across the four waves.

We can also see this visually by using the predicted scores:

gender_predictions <- expand_grid(
  wave0 = 0:3,
  gndr = factor(levels(usl$gndr), levels = levels(usl$gndr))
) |>
  mutate(predicted = predict(m4, newdata = pick(everything()), re.form = NA))

ggplot(gender_predictions, aes(wave0 + 1, predicted, colour = gndr)) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 2) +
  theme_bw() +
  labs(x = "Wave", y = "Predicted mental health", colour = "Gender")

Men begin with higher predicted mental-health scores. The lines are not parallel: the positive interaction means that women decline more slowly, so the difference becomes smaller across the four waves.

We can expand the model to see an important example of interaction effects. Below, we include the effect of degree, another variable we treat as a time constant, as well as the interaction with time. Given the small effect we observed earlier, we keep the main effect of sex but exclude the interaction with time.

m5 <- lmer(
  data = usl,
  sf12mcs ~ 1 + wave0 + gndr + degree + degree:wave0 +
    (1 + wave0 | pidp)
)

summary(m5)

## Linear mixed model fit by REML ['lmerMod']
## Formula: sf12mcs ~ 1 + wave0 + gndr + degree + degree:wave0 + (1 + wave0 |  
##     pidp)
##    Data: usl
## 
## REML criterion at convergence: 934129.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.3917 -0.4289  0.1551  0.5610  3.1577 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  pidp     (Intercept) 38.931   6.239         
##           wave0        3.002   1.733    -0.37
##  Residual             59.121   7.689         
## Number of obs: 127868, groups:  pidp, 48573
## 
## Fixed effects:
##                       Estimate Std. Error t value
## (Intercept)           51.72881    0.08395 616.162
## wave0                 -0.55139    0.03650 -15.107
## gndrFemale            -1.37171    0.07120 -19.267
## degreeNo degree       -1.01131    0.09068 -11.153
## wave0:degreeNo degree  0.19794    0.04543   4.357
## 
## Correlation of Fixed Effects:
##             (Intr) wave0  gndrFm dgrNdg
## wave0       -0.503                     
## gndrFemale  -0.473 -0.004              
## degreeNdegr -0.718  0.468 -0.003       
## wv0:dgrNdgr  0.404 -0.804  0.005 -0.568
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.0105688 (tol = 0.002, component 1)

In the final model, the intercept of 51.729 is the expected first-wave score for men with a degree. Respondents without a degree score 1.011 points lower at the first wave. The degree group declines by 0.551 points per wave, while the positive interaction of 0.198 indicates a slower decline among respondents without a degree. The education gap therefore narrows across the four waves.

education_predictions <- expand_grid(
  wave0 = 0:3,
  gndr = factor("Male", levels = levels(usl$gndr)),
  degree = factor(levels(usl$degree), levels = levels(usl$degree))
) |>
  mutate(predicted = predict(m5, newdata = pick(everything()), re.form = NA))

ggplot(education_predictions, aes(wave0 + 1, predicted, colour = degree)) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 2) +
  theme_bw() +
  labs(x = "Wave", y = "Predicted mental health", colour = "Education")

Conclusions

The multilevel model for change can be easily expanded to include time-constant predictors. For the most part, the interpretation is similar to that from multiple regression. Care is needed with the interpretation when including interactions with time to allow for different rates of change for subgroups. You have seen above some example code for using predicted scores and visualizations to better understand your model’s results. You can also check how to include time varying predictors in this blogpost as well as investigating nonlinear change in time here.