Modeling two-level longitudinal data

Prof. Maria Tackett

Mar 04, 2024

Announcements

  • Quiz 03 - Tue, Mar 04 at 9am - Thu, Mar 06 at noon

    • Covers readings and lectures in Feb 19 - 28 lectures

    • Correlated data, multilevel models

Topics

  • Exploratory data analysis for longitudinal data
  • Fit and compare longitudinal models
  • Inference for fixed and random effects

Modeling two-level longitudinal data

Data: Charter schools in MN

Today’s data set contains standardized test scores and demographic information for schools in Minneapolis, MN from 2008 to 2010. The data were collected by the Minnesota Department of Education. Understanding the effectiveness of charter schools is of particular interest, since they often incorporate unique methods of instruction and learning that differ from public schools.

  • MathAvgScore: Average MCA-II score for all 6th grade students in a school (response variable)
  • urban: urban (1) or rural (0) location school location
  • charter: charter school (1) or a non-charter public school (0)
  • schPctfree: proportion of students who receive free or reduced lunches in a school (based on 2010 figures).
  • year08: Years since 2008

Data

schoolName year08 urban charter schPctfree MathAvgScore
RIPPLESIDE ELEMENTARY 0 0 0 0.363 652.8
RIPPLESIDE ELEMENTARY 1 0 0 0.363 656.6
RIPPLESIDE ELEMENTARY 2 0 0 0.363 652.6
RICHARD ALLEN MATH&SCIENCE ACADEMY 0 1 1 0.545 NA
RICHARD ALLEN MATH&SCIENCE ACADEMY 1 1 1 0.545 NA
RICHARD ALLEN MATH&SCIENCE ACADEMY 2 1 1 0.545 631.2

Assess missingness

Missing data is common in longitudinal data. Before starting the analysis, it is important to understand the missing data patterns. Use the skim function from the skimr R package to get a quick view of the missingness.

library(skimr)
charter |> skim() |> select(skim_variable, n_missing, complete_rate)
# A tibble: 9 × 3
  skim_variable n_missing complete_rate
  <chr>             <int>         <dbl>
1 schoolid              0         1    
2 schoolName            0         1    
3 urban                 0         1    
4 charter               0         1    
5 schPctnonw            0         1    
6 schPctsped            0         1    
7 schPctfree            0         1    
8 year08                0         1    
9 MathAvgScore        121         0.935

Closer look at missing data pattern

MathAvgScore0_miss MathAvgScore1_miss MathAvgScore2_miss n
0 0 0 540
0 0 1 6
0 1 0 4
0 1 1 7
1 0 0 25
1 0 1 1
1 1 0 35

Code

charter |>
  select(schoolid, schoolName, year08, MathAvgScore) |>
  pivot_wider(id_cols = c(schoolid, schoolName), names_from = year08,
              names_prefix = "MathAvgScore.", values_from = MathAvgScore) |> 
  mutate(MathAvgScore0_miss = if_else(is.na(MathAvgScore.0), 1, 0),
         MathAvgScore1_miss = if_else(is.na(MathAvgScore.1), 1, 0),
         MathAvgScore2_miss = if_else(is.na(MathAvgScore.2), 1, 0)) |> 
  count(MathAvgScore0_miss, MathAvgScore1_miss, MathAvgScore2_miss) |> 
  kable()

Dealing with missing data

  • Complete case analysis: Only include schools with complete data for all three years.
    • Would remove 12.6% of observations in this data.
  • Last observation carried forward: Keep the last observation from each group (school) and conduct analysis for independent observations.
  • Impute missing observations: “Fill in” values of missing observations using the typical observed trends from groups with similar covariates.
  • Apply multilevel methods: Estimate patterns using available data recognizing that trends for groups with complete data are more precise than for those with fewer measurements. This is under the condition that the probability of missingness does not depend on unobserved predictors or the response.

Dealing with missing data

For your group’s assigned approach?

  • What is an advantage?

  • What is a disadvantage?

One person post your group’s response in the lectures channel on Slack.

03:00

Strategy for building multilevel models

  • Conduct exploratory data analysis for Level One and Level Two variables.

  • Fit model with no covariates to assess variability at each level.

  • Create Level One models. Start with a single term, then add terms as needed.

  • Create Level Two models. Start with a single term, then add terms as needed. Start with equation for intercept term.

  • Begin with the full set of variance components, then remove variance terms as needed.

Alternate model building strategies in BMLR Section 8.6

Exploratory data analysis

Given the longitudinal structure of the data, we are able to answer questions at two levels

  • Level One (within school): How did average math scores for a given school change over time?

  • Level Two (between schools): What is the effect of school-specific covariates on the average math scores in 2008 and the rate of change from 2008 to 2010?

We can conduct exploratory data analysis at both levels, e.g.,

  • Univariate and bivariate EDA
  • Lattice plots
  • Spaghetti plots

Application exercise

📋 sta310-sp24.netlify.app/ae/lec-14-longitudinal


See BMLR Section 9.3 for full exploratory data analysis.

Unconditional means model

Start with the unconditional means model, a model with no covariates at any level. This is also called the random intercepts model.

Level One : \(Y_{ij} = a_{i} + \epsilon_{ij}, \hspace{5mm} \epsilon_{ij} \sim N(0, \sigma^2)\)

Level Two: \(a_i = \alpha_0 + u_i, \hspace{5mm} u_{i} \sim N(0, \sigma^2_u)\)

Write the composite model.

Intraclass correlation

The intraclass correlation is the relative variability between groups

\[\hat{\rho} = \frac{\text{Between variability}}{\text{Total variability}} = \frac{\hat{\sigma}_u^2}{\hat{\sigma}^2_u + \hat{\sigma}^2}\]

  • Fit the unconditional means model and calculate the intraclass correlation in the AE.

  • What does this value tell you?

Unconditional growth model

A next step in the model building is the unconditional growth model, a model with Level One predictors but no Level Two predictors.

Level One

\[ Y_{ij} = a_i + b_iYear08_{ij} + \epsilon_{ij}, \hspace{8mm} \epsilon_{ij} \sim N(0, \sigma^2) \]

Level Two:

\[ \begin{aligned} &a_i = \alpha_0 + u_i \hspace{20mm} b_i = \beta_0 + v_i \\ &\left[ \begin{array}{c} u_{i} \\ v_{i} \end{array} \right] \sim N \left( \left[ \begin{array}{c} 0 \\ 0 \end{array} \right], \left[ \begin{array}{cc} \sigma_{u}^{2} & \sigma_{uv} \\ \sigma_{uv} & \sigma_{v}^{2} \end{array} \right] \right) \end{aligned} \]

Unconditional growth model

  • Write the composite model.
  • Fit the unconditional growth model in the AE.
  • What can we learn from this model?

Pseudo \(R^2\)

We can use Pseudo R2 to explain changes in variance components between two models

Note: This should only be used when the definition of the variance component is the same between the two models

\[\text{Pseudo }R^2 = \frac{\hat{\sigma}^2(\text{Model 1}) - \hat{\sigma}^2(\text{Model 2})}{\hat{\sigma}^2(\text{Model 1})}\]

Calculate the \(\text{Pseudo }R^2\) in the AE to estimate the change of within school variance between the unconditional means and unconditional growth models.

Model with school-level covariates

Fit a model with school-level covariates that takes the following form:

Level One

\[\begin{equation*} Y_{ij}= a_{i} + b_{i}Year08_{ij} + \epsilon_{ij} \hspace{10mm} \epsilon_{ij} \sim N(0, \sigma^2) \end{equation*}\]

Level Two

\[\begin{align*} a_{i} & = \alpha_{0} + \alpha_{1}charter_i + \alpha_{2}urban_i + \alpha_{3}schpctfree_i + u_{i} \\ b_{i} & = \beta_{0} + \beta_{1}charter_i + \beta_{2}urban_i + v_{i} \\[10pt] &\left[ \begin{array}{c} u_{i} \\ v_{i} \end{array} \right] \sim N \left( \left[ \begin{array}{c} 0 \\ 0 \end{array} \right], \left[ \begin{array}{cc} \sigma_{u}^{2} & \sigma_{uv} \\ \sigma_{uv} & \sigma_{v}^{2} \end{array} \right] \right) \end{align*}\]

Model with school-level covariates

  • Write out the composite model.
  • Fit the model in R.
  • Use the model to describe how the average math scores differed between charter and non-charter schools.

Consider a simpler model

Would a model without the effects \(v_i\) and \(\rho_{uv}\) be preferable?

Level One

\[\begin{equation*} Y_{ij}= a_{i} + b_{i}Year08_{ij} + \epsilon_{ij} \hspace{10mm} \epsilon_{ij} \sim N(0, \sigma^2) \end{equation*}\]

Level Two

\[\begin{align*} a_{i} & = \alpha_{0} + \alpha_{1}charter_i + \alpha_{2}urban_i + \alpha_{3}schpctfree_i + u_{i} \\ & u_i \sim N(0, \sigma_u^2) \\[5pt] b_{i} & = \beta_{0} + \beta_{1}charter_i + \beta_{2}urban_i \end{align*}\]

In this model, the effect of year is the same for all schools with a given combination of charter and urban

Full and reduced models

full_model <- lmer(MathAvgScore ~ charter + urban + schPctfree + 
                     charter:year08  + urban:year08 + year08 + 
                     (year08|schoolid), REML = F, data = charter) 


reduced_model <-  lmer(MathAvgScore ~ charter + urban + schPctfree + 
                     charter:year08  + urban:year08 + year08 +
                       (1 | schoolid), REML = F, data = charter)

Notes on drop-in-deviance test

anova(full_model, reduced_model, test = "LRT") |>
  kable(digits = 3)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
reduced_model 9 9952.992 10002.11 -4967.496 9934.992 NA NA NA
full_model 11 9953.793 10013.83 -4965.897 9931.793 3.199 2 0.202

  • \(\chi^2\) test is conservative, i.e., the p-values are larger than they should be, when testing random effects at the boundary ( e.g., \(\sigma_v^2 = 0\)) or those with bounded ranges (e.g., \(\rho_{uv}\) )

  • If you observe small p-values, you can feel relatively certain the tested effects are statistically significant

  • Use bootstrap methods to obtain more accurate p-values

References

🔗 STA 310 - Spring 2024

Boedeker, Peter. 2019. “Hierarchical Linear Modeling with Maximum Likelihood, Restricted Maximum Likelihood, and Fully Bayesian Estimation.” Practical Assessment, Research, and Evaluation 22 (1): 2.
Faraway, Julian J. 2016. Extending the Linear Model with r: Generalized Linear, Mixed Effects and Nonparametric Regression Models. CRC press.
Roback, Paul, and Julie Legler. 2021. Beyond multiple linear regression: applied generalized linear models and multilevel models in R. CRC Press.
Singer, Judith D, and John B Willett. 2003. Applied Longitudinal Data Analysis: Modeling Change and Event Occurrence. Oxford university press.