Longitudinal data consist of repeated measures collected on the same subjects over time. One common method for modeling such data is with a repeated measures analysis of variance (RM ANOVA). However, using a RM ANOVA treats discrete time points as categorical, with only the mean at each time point formally considered. This method works in certain controlled situations, but it does not offer enough flexibility to be used practically outside of these specific situations. Linear mixed-effect regression (LMER) is another method for analyzing longitudinal data that offers greater flexibility than the RMANOVA.

Linear Mixed Effect Regression

LMER is a method for analyzing longitudinal data, and it can be viewed as an extension of traditional multiple regression. The extension is made by introducing individual-level terms to the regression model that index variability among the subjects. These terms are known as random effects. The random effects are in contrast to the fixed effects, which are similar to the traditional regression coefficients. Fixed effects are constant among individuals and index group-level change. It is the mix of fixed and random effects that is the genesis of the term linear mixed effects regression. Random effects vary among individuals and index deviations from the group. In the context of longitudinal models, the random effects reflect variation among individual change curves (Long, p. 9).

\[y_{ij}=\left(\beta_0+U_{0j}\right)+\left(\beta_1+U_{1j}\right)\left(Time_{ij}\right)+\epsilon_{ij}\]



Advantages

  • LMER can accommodate missing data. Subjects who have response data for at least one time point can be included in the analysis. However, valid inference under these conditions is predicated on assumptions regarding the mechanism that caused the missingness (more on missing data below).

  • LMER provides control over the number and nature of terms used to model change over time. Regardless of the number of time points, the researcher can fit a low order polynomial, or fit other types of curves that directly address research questions.

  • LMER is very flexible regarding the structure of the data. Timing of observations can vary among subjects, and the distance between adjacent time points need not be equal.

  • LMER allows for various types of predictors. Predictors can be dynamic and/or static, quantitative and/or categorical.

Types of Variability

  • Within-subjects variability: changes in the response variable over time (i.e variability among the repeated measures). Any predictor that changes over time accounts for within-subjects variability.

  • Between-subject variability: variation due to individual differences. Any predictor that is constant over time-but not constant among subjects-accounts for between-subjects variability.

Random Effects

  • In LMER, each subject has their own intercept and slope.

  • An individual’s intercept and slope are expressed in relation to the group intercept and slope.

  • Random effects are represented by the deviation of the individual regression lines from the group regression line.

  • By using random effects, you can collapse the time dimension and create one-row-per-subject data, where the rows are independent of one another (thereby removing statistical dependence).

Assumptions

  • Subjects have been randomly sampled from a population (hence, random effects).

  • The random effects have a joint normal distribution.

  • The random effects are correlated, but they are independent of the random error.

  • Within-subject error, which is the difference between a subject’s observed data points and the subject’s random slope, is normally distributed, with a mean equal to 0 and constant variance.

LMER vs RM ANOVA

Let’s contrast two methods for analyzing longitudinal data: LMER and RM ANOVA. Doing so will help us understand when and when not to use them. A more detailed contrast of these two methods can be found in an article by Lohse, Shen, and Kozlowski (2020).

Model concept

  • RM ANOVA
    • Compares means of a continuous outcome stratified by one or more categorical variable(s) to the grand mean.
    • Individuals are treated as a factor with error aggregated from each individual’s mean of repeated measures and partitioned as intra-individual variance from the error term.
  • Mixed-Effects Regression
    • Accounts for correlations in data with clustered or nested structure.
    • In longitudinal models, change over time is considered a within-person factor accounting for within-person correlations across time points and estimating error as residuals from each individual’s trajectory, and between-person error is accounted for as random effects in a correlation matrix, which can be explained by fixed covariate associations with trajectory parameters.

Modeling Outcomes Over Time

  • RM ANOVA
    • Addresses questions about mean-differences, after between-subject variance is removed.
    • Discrete time points are treated as categorical, with only the mean at each time point formally considered.
  • Mixed-Effect Regression
    • Time is modeled explicitly as a trajectory for each individual.
    • The shape of the trajectory is determined by fitting progressively more complex mathematical functions.

Variability Over Time

  • RM ANOVA
    • Assumes common, identically-time data collections. This can lead to increased variability in the discrete time points that is really due to variation in when data were collected.
  • Mixed-Effect Regression
    • Can accommodate variability in spacing of time points and in the actual timing of individual data collection.
    • Can account for increased heterogeneity of the data over time, but residuals still need to be homogeneous.

Data Missing on the Outcome

  • RM ANOVA
    • Missing data generally cannot be accommodated.
    • If data are missing at random (MAR), multiple imputation can be used for estimation.
    • If data are missing not at random (MNAR), listwise deletion will reduce statistical power and potentially introduce bias.
  • Mixed-Effect Regression
    • Data that are MAR can be accommodated without exclusion or imputation.
    • Data that are MNAR can be fit, but factors associated with missingness need to identified and included in the model.
    • Data missing at key moments may lead to poor model fit.

Data Missing on the Explanatory Variables

  • RM ANOVA
    • Missing between-person data in explanatory variables/covariates cannot be accommodated.
    • Cases need to either dropped from the model or imputed.
  • Mixed-Effect Regression
    • Missing between-person data in explanatory variables/covariates cannot be accommodated.
    • Cases need to either dropped from the model or imputed.

Including Covariates that Change Over Time

  • RM ANOVA
    • Time-varying covariates cannot be included in an RM ANOVA model.
  • Mixed-Effect Regression
    • Time varying covariates can be included, but you need to careful about collinearity and variance at both the between- and within-subject levels.

Pre-Modeling: Missing Data

Prior to statistical modeling it’s critical to develop an understanding of the data you’re working with. Specifically, it’s important to investigate which mechanisms of missing data are potentially at play.

Visualizations can provide insight about which mechanism is likely at play. It is important to note that there is almost never a clear answer, and it is up to the researcher to use their best judgment and decide which mechanism best explains the observed data.

library(ggplot2)
library(dplyr)
library(lme4)
library(lmerTest)
data <- read.csv("../ACRM_2018/data/data_session1.csv", stringsAsFactors = TRUE)
head(data)
str(data)
## 'data.frame':    720 obs. of  19 variables:
##  $ subID         : Factor w/ 40 levels "s01","s02","s03",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ month         : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ sex           : Factor w/ 2 levels "f","m": 1 1 1 1 1 1 1 1 1 1 ...
##  $ female        : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ bmi_z         : num  -0.44 -0.44 -0.44 -0.44 -0.44 -0.44 -0.44 -0.44 -0.44 -0.44 ...
##  $ rasch_FIM     : num  12.6 24.9 40.3 43.9 47.7 ...
##  $ rasch_FIM_MAR : num  12.6 24.9 NA NA NA ...
##  $ rasch_FIM_MNAR: num  12.6 24.9 40.3 43.9 47.7 ...
##  $ rasch_FIM_LOCF: num  12.6 24.9 40.3 43.9 47.7 ...
##  $ inpatient_time: Factor w/ 7 levels "1wk","2wks","3wks",..: 5 5 5 5 5 5 5 5 5 5 ...
##  $ los           : int  32 32 32 32 32 32 32 32 32 32 ...
##  $ AIS_grade     : Factor w/ 3 levels "C1-4","C5-8",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ therapy_hours : num  14.6 14.7 14.6 14.5 8.34 ...
##  $ PT_hours      : num  7.3 7.3 7.3 7.2 4.17 ...
##  $ OT_hours      : num  7.4 7.4 7.3 7.3 4.23 ...
##  $ SUM_PT_hours  : num  33.3 33.3 33.3 33.3 33.3 ...
##  $ SUM_OT_hours  : num  33.6 33.6 33.6 33.6 33.6 ...
##  $ AVE_PT_hours  : num  6.65 6.65 6.65 6.65 6.65 ...
##  $ AVE_OT_hours  : num  6.73 6.73 6.73 6.73 6.73 ...
set.seed(123)
data$rasch_FIM_MCAR <- data$rasch_FIM
MCAR_sample <- as.numeric(sample(row.names(data), 100))
data$rasch_FIM_MCAR[MCAR_sample] <- NA
data[,c(1, 2, 6:8, 20)]
ggplot(data, aes(month, rasch_FIM, group = subID, col = subID)) +
  geom_line() +
  theme_classic() +
  theme(legend.position = "none")

ggplot(data, aes(month, rasch_FIM_MCAR, group = subID, col = subID)) +
  geom_line() +
  theme_classic() +
  theme(legend.position = "none")

ggplot(data, aes(month, rasch_FIM_MAR, group = subID, col = subID)) +
  geom_line() +
  theme_classic() +
  theme(legend.position = "none")

ggplot(data, aes(month, rasch_FIM_MNAR, group = subID, col = subID)) +
  geom_line() +
  theme_classic() +
  theme(legend.position = "none")

data$missing <- ifelse(is.na(data$rasch_FIM_MNAR), 1, 0)
summary(data$missing)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  1.0000  0.6181  1.0000  1.0000
# Plotting missing values as a function of IV
ggplot(data=data, aes(x=month, y=missing))+
  geom_point(aes(col=subID), position = position_jitter(width=0.1, height=0.1))+
  stat_smooth(method="lm", col="black", se = F)+
  theme_classic()+
  scale_x_continuous(name="Month", limits=c(0,18))+
  scale_y_continuous(name="Missing")+
  theme(legend.position = "none")

# Since we are missing the DV to begin with,
# we need to shift our plotting strategy
ggplot(data=data, aes(x=rasch_FIM_MNAR))+
  geom_histogram(aes(y = ..density..), fill="grey", col="black", alpha=0.5)+
  geom_density(fill="lightblue", col="black", alpha=0.5)+
  theme_classic()+
  scale_x_continuous(name="Rash FIM", limits=c(0,100))+
  theme(legend.position = "none")

Pre-Modeling: Graphing

# Random sample 
data %>% 
  filter(subID %in% sample(subID, 10)) %>% 
  ggplot(aes(month, rasch_FIM, group = subID, col = subID)) +
  geom_line() +
  theme_classic()

# Facet wrap with linear model
subset(data, as.numeric(subID) < 17) %>% 
ggplot(aes(month, rasch_FIM, group = subID)) +
  geom_line() + 
  geom_smooth(method = "lm", col = "red", se = F) +
  facet_wrap(~subID)

# Facet wrap with polynomial
subset(data, as.numeric(subID) < 17) %>% 
ggplot(aes(month, rasch_FIM, group = subID)) +
  geom_line() + 
  stat_smooth(method = "lm", formula = y ~ poly(x, 2), size = .5) +
  facet_wrap(~subID) +
  theme_bw()

# Graph by missing values 
data %>%
  group_by(subID) %>% 
  mutate(rash_NA_flag = ifelse(is.na(rasch_FIM_MNAR), 1, 0),
         Total_NA = sum(rash_NA_flag)) %>%
  filter(Total_NA < 10) %>% 
  ggplot(aes(month, rasch_FIM, group = subID)) +
  geom_line() + 
  stat_smooth(method = "lm", formula = y ~ poly(x, 3), col = "red") +
  facet_wrap(~subID)

# Difference in outcomes between high and low missing value groups?
data %>%
  group_by(subID) %>% 
  mutate(rash_NA_flag = ifelse(is.na(rasch_FIM_MNAR), 1, 0),
         Total_NA = sum(rash_NA_flag),
         Missing_Cat = ifelse(Total_NA >= 11, "High", "Low")) %>%
  ggplot(aes(month, rasch_FIM_MNAR, group = subID, col = Missing_Cat)) +
  geom_line(alpha = .3) +
  geom_point(alpha = .3) +
  scale_color_brewer(palette = "Dark2") +
  theme_bw() +
  geom_smooth(aes(month, rasch_FIM, group = Missing_Cat), 
              method = "lm",
              formula = y ~ poly(x, 4))

# Plotting means
ggplot(data, aes(month, rasch_FIM)) +
  geom_point(size = .5, alpha = .5) +
  stat_summary(fun.y = mean, geom = "point", col = "red", size = 2) +
  stat_summary(fun.y = mean, geom = "line", linetype = 2, size = 1.5) +
  theme_classic()

ggplot(data, aes(month, rasch_FIM)) +
    stat_summary(aes(line = AIS_grade), fun.y = mean, geom = "line") +
    stat_summary(aes(shape = AIS_grade), fun.y = mean, geom = "point", size = 2) +
  theme_bw()

The Fundamental LMER Models

1. Unconditional Model: Time as a predictor of the DV

2. Intercept Effect: Time and a Categorical Variable as predictors of the DV

3. Intercept and Slope Effects: Time, Categorical Variable, and Time x Categorical Variable as predictors of DV

Unconditional Model: Random intercepts

UM_RI <- lmer(rasch_FIM ~ month + (1 | subID), data)
summary(UM_RI)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rasch_FIM ~ month + (1 | subID)
##    Data: data
## 
## REML criterion at convergence: 4915.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4066 -0.4538  0.2031  0.6423  2.3951 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subID    (Intercept) 100.88   10.044  
##  Residual              43.89    6.625  
## Number of obs: 720, groups:  subID, 40
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  24.43435    1.66957  45.39941   14.63   <2e-16 ***
## month         2.15474    0.04759 679.00000   45.28   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##       (Intr)
## month -0.271
fixef(UM_RI) # get the fixed effects
## (Intercept)       month 
##   24.434347    2.154744
ranef(UM_RI)$subID # random effects
coef(UM_RI)$subID # fixed effect + random effects
# We'll make a smaller dataset with the predictions for these 10: 
first10 <- data[c(1:180), ]
PRED <- data.frame(subID = c("s01","s02","s03","s04","s05","s06","s07","s08","s09","s10"),
                   Intercepts = c(coef(UM_RI)$subID[c(1:10), 1])) 
ggplot(first10, aes(x = month, y = rasch_FIM)) +
  geom_point(fill = "grey", pch = 21, size = 2, stroke = 1.25) +
  geom_line(aes(group = subID)) + facet_wrap( ~ subID, ncol = 3) +
  scale_x_continuous(name = "Time from Admission (Months)", breaks = c(0:18)) +
  scale_y_continuous(name = "Rasch-Scaled FIM Score (0-100)", limits = c(0, 100)) +
  geom_abline(aes(intercept = Intercepts, slope = 0), col = "red", lwd = 1.5, PRED)

Unconditional Model: Random intercepts & Random Slopes

UM_RI_RS <- lmer(rasch_FIM ~ month + (month | subID), data)
PRED_RI_RS <- data.frame(subID = c("s01","s02","s03","s04","s05","s06","s07","s08","s09", "s10"), 
                         Intercepts = c(coef(UM_RI_RS)$subID[c(1:10), 1]), 
                         Slopes = c(coef(UM_RI_RS)$subID[c(1:10), 2])) 
ggplot(first10, aes(x = month, y = rasch_FIM)) +
  geom_point(fill = "grey", pch = 21, size = 2, stroke = 1.25) +
  geom_line(aes(group = subID)) + facet_wrap( ~ subID, ncol = 3) +
  scale_x_continuous(name = "Time from Admission (Months)", breaks = c(0:18)) +
  scale_y_continuous(name = "Rasch-Scaled FIM Score (0-100)", limits = c(0, 100)) +
  geom_abline(aes(intercept = Intercepts, slope = Slopes), col = "red", lwd = 1.5, PRED_RI_RS) 

anova(UM_RI)
lmerTest::ranova(UM_RI)
# Model Comparison
anova(UM_RI, UM_RI_RS)
# REML differs from ML in that it includes a bias adjustment. The bias adjustment can be desirable when the sample size is not large. However, the bias adjustment is generally undesirable for model comparison. Use ML deviance instead of REML deviance when making model comparisons.

References

  1. Long, Jeffrey D.. Longitudinal Data Analysis for the Behavioral Sciences Using R. SAGE Publications. Kindle Edition.
  2. Lohse, K. R., Shen, J., & Kozlowski, A. J. (2020). Modeling Longitudinal Outcomes: A Contrast of Two Methods. Journal of Motor Learning and Development, 1(aop), 1-21. doi: https://doi.org/10.1123/jmld.2019-0007