Researchers studying longitudinal data routinely center their predictors to isolate between- and within-cluster contrasts (Enders and Tofighi 2007). This within-cluster centering is usually an easy data-manipulation step. However, centering variables on the observed means can bias the resulting estimates, a problem that is avoided with latent mean centering. In this entry, we study how to latent-mean-center variables in multilevel models using brms.

# Packageslibrary(knitr)library(brms)library(ggthemes)library(scales)library(kableExtra)library(posterior)library(tidyverse)# Options for samplingoptions(brms.backend ="rstan",mc.cores = parallel::detectCores(logical =FALSE) )# Function for tablesk2 <-function(x, escape =TRUE) { x %>%kbl(digits =2, escape = escape) %>%kable_classic_2("striped", full_width =FALSE, html_font ="Arial")}# Plotting themetheme_set(theme_few() +theme(axis.title.y =element_blank(),legend.title =element_blank(), panel.grid.major =element_line(linetype ="dotted", linewidth = .1),legend.position ="bottom", legend.justification ="left" ))# Download and uncompress McNeish and Hamaker materials if not yet donepath <-"materials/materials.zip"if (!file.exists(path)) {dir.create("materials", showWarnings =FALSE)download.file("https://files.osf.io/v1/resources/wuprx/providers/osfstorage/5bfc839601593f0016774697/?zip=",destfile = path )unzip(path, exdir ="materials")}

Warning

I am currently studying this model implementation in more detail. It seems that there are some pieces missing, and results regarding some parameters seem unreliable. So if you want to use this in your own work, please ensure appropriate model performance with a parameter recovery simulation. Take everything below with a grain of salt. I’ll update here once I learn more about this model’s performance.

The earlier drafts and mistakes I made in coding the brms model up can be found in the Git history of this file 😄

Introduction

Within cluster centering, or person-mean centering (psychologists’ clusters are typically persons), is an easy but essential data processing step. For example consider the example data of 100 people’s ratings of urge to smoke and depression, collected over 50 days with one response per day (McNeish and Hamaker 2020)^{1}, shown in Table 1 and Figure 1.

^{1} Grab a free copy at https://osf.io/j56bm/download. I couldn’t figure if this example data is real or simulated, or what the measurement instruments were.

Code

d <-read_csv("materials/Data/Two-Level Data.csv", col_names =c("urge", "dep", "js", "hs", "person", "time")) %>%select(-hs, -js) %>%relocate(person, time, 1)d <- d %>%# Grand mean center both variablesmutate(across(c(urge, dep), list(c =~ . -mean(.)))) %>%group_by(person) %>%mutate(across(c(urge_c, dep_c), list(# Between-person center (= person's mean)b =~mean(.), # Within-person center (= deviation from person's mean)w =~ . -mean(.) ),.names ="{.col}{.fn}" ) ) %>%ungroup()

Table 1: Example longitudinal data (McNeish & Hamaker, 2020)

person

time

urge

dep

urge_c

dep_c

urge_cb

urge_cw

dep_cb

dep_cw

1

1

0.34

0.43

0.33

0.42

-0.31

0.64

-0.19

0.61

1

2

-0.48

-0.68

-0.49

-0.69

-0.31

-0.18

-0.19

-0.51

1

3

-4.44

-1.49

-4.45

-1.50

-0.31

-4.14

-0.19

-1.31

1

4

-4.19

-0.74

-4.20

-0.75

-0.31

-3.89

-0.19

-0.56

1

5

-0.91

-0.52

-0.92

-0.53

-0.31

-0.61

-0.19

-0.34

2

1

1.65

0.68

1.64

0.66

-0.48

2.12

-0.02

0.68

2

2

0.31

1.49

0.30

1.48

-0.48

0.78

-0.02

1.49

2

3

0.46

0.03

0.45

0.02

-0.48

0.94

-0.02

0.03

2

4

-1.09

-1.02

-1.10

-1.03

-0.48

-0.62

-0.02

-1.02

2

5

1.67

1.07

1.66

1.06

-0.48

2.14

-0.02

1.07

Note: ^{} Only displaying two participants' first five observations. _c: grand mean centered; _cb: between-person centered (ie. person mean); _cw: within-person centered.

The first four variables in Table 1 are the original data values, indicating the person and timepoint of measurement, urge to smoke, and depression. I’ve then created the grand-mean, between-person, and within-person variables by simple data transformations. Between-person centered variables are person-specific means, and within-person centered variables are deviations around that person’s mean.

However, the person-mean is an unknown quantity, and centering on the observed value rather than an estimate of the true “latent” quantity can be problematic. Specifically, observed mean centering leads to Nickell’s (negative bias in autoregressive effects) and Lüdtke’s (bias in other time-varying effects) biases (McNeish and Hamaker 2020, 617–18). Essentially these problems arise from not considering that the person means are unobserved, latent quantities, but instead treating them as values known without uncertainty.

So, what to do? McNeish and Hamaker (2020) and others discuss latent mean centering, which accounts for uncertainty in the person-means appropriately, and thus debiases the estimated coefficients. Latent mean centering is done inside the model, and means treating the means as estimated parameters. However, I have only been able to find examples that do this latent mean centering in MPlus, such as (McNeish and Hamaker 2020). Therefore my goal here is to reproduce their model with the free and open source software Stan front-end brms.

Single-level AR(1) model

To begin with, we replicate the authors’ basic N=1 model predicting the urge to smoke from the urge to smoke on a previous measurement occasion, and the current level of depression. Because we are modelling one person’s data only, there is no need for centering, but this model serves as a useful starting point for our quest.

Following (McNeish and Hamaker 2020), we assume that Urge at time \(t\) is normally distributed around a mean \(\mu_t\) with standard deviation \(\sigma\). We then model the mean on an intercept, on Urge at the previous measurement occasion, and on the current level of depression^{2}:

^{2} I’ve used the more common \(\phi\) (phi) throughout than the \(\varphi\) used by M&P

This is straightforward. We first create a lagged urge variable, and then fit the model. Notice though that this will lead to one missing data point because the first value doesn’t have a lagged value. We confirm that our estimates are in line with those reported in the paper

Code

# Estimate modelfit_p5 <-brm( urge ~ urge1 + dep,family =gaussian(),data = d %>%# Pick one individual (same as used in M&H2020)filter(person ==5) %>%# Create lagged urgemutate(urge1 =lag(urge)),file ="fit_p5")

Code

# Show table of coefficients' posterior summariesas_draws_df(fit_p5) %>%select(1:4) %>%mutate(sigma_sq = sigma^2, .keep ="unused") %>%summarise_draws(median, ~quantile2(., probs =c(.025, .975))) %>%mutate(variable =str_c( variable, c(" ($\\alpha$)", " ($\\phi$)", " ($\\beta$)", " ($\\sigma^2$)") ),across(c(median, q2.5, q97.5), ~number(., .01)),`Result (brms)`=str_glue("{median} [{q2.5}, {q97.5}]"),Authors =c("0.07 [-0.24, 0.39]", "0.35 [0.25, 0.46]", "2.43 [2.12, 2.75]", "1.36 [0.92, 2.20]" ) ) %>%select(-c(median:q97.5)) %>%k2(escape =FALSE) %>%footnote("I've no idea how to render the math here, SRY", footnote_as_chunk =TRUE)

variable

Result (brms)

Authors

b_Intercept ($\alpha$)

0.08 [-0.25, 0.40]

0.07 [-0.24, 0.39]

b_urge1 ($\phi$)

0.35 [0.24, 0.47]

0.35 [0.25, 0.46]

b_dep ($\beta$)

2.43 [2.12, 2.76]

2.43 [2.12, 2.75]

sigma_sq ($\sigma^2$)

1.34 [0.92, 2.05]

1.36 [0.92, 2.20]

Note: ^{} I've no idea how to render the math here, SRY

Multilevel AR(1) model

Above, we modelled a single person’s urge to smoke on their previous urge to smoke and current depression. Here, we attempt to model 100 individuals’ data in a single multilevel model. Before worrying about latent mean centering, we can estimate this model using the observed mean centered values shown in Table Table 1. The authors’ model of these data is

where we now have a subscript \(i\) for participants, parameters with bars (population-level) and without (person-specific, with subscripts \(i\)), and the variance-covariance matrix for the latter, where all covariances are set to zero as in (McNeish and Hamaker 2020). Notice that the the subscripted parameters are deviations with a mean of zero, so we can talk about e.g. \(\bar{\alpha} + \alpha_2\) as person 2’s intercept. \(\text{Urge}^c_{(t-1)i}\) and \(\text{Dep}^c_{ti}\) are the within-person centered values of urge to smoke on the previous timepoint and depression, respectively (_cw values in Table Table 1.) We can fit this model with brms with a small modification of the previous model

We now estimated the model using observed person mean centering. These estimates are very close to the ones reported in authors’ Table 4, because there are so many observations per person. (Note this is probably not the best example because of this, and also because I don’t really know what the data are from the paper). Importantly, we have not estimated the latent depression variable, so that is missing

The problem then boils down to figuring out how to get the quantities \(\text{Urge}^c_{(t-1)i}\) and \(\text{Dep}^c_{ti}\). Usually, we calculate them from data as deviations from the person’s observed mean, like we did above in Table 1. However, here’ we want to use latent-mean centering:

It turns out that specifying this model with latent mean centering is fairly straightforward with brms. First, we will need to specify a non-linear formula where we name all parameters, and then another one that specifies that one of the predictors is a parameter too. Thanks to Mauricio Garnier-Villarreal, Ethan McCormick, Simon Brauer, Joran Jongerling, and others who helped out with my Stan discourse question to figure out the syntax!

Here goes. We specify a formula of urge on the named parameters and predictors, as you do with brms’ nonlinear formulas. Then in the subsequent lines, each parameter is specified their own model. The trick is to predict the latent means inside another nlf(), and then the predictor there in another model formula. That’s it! And because this is a nonlinear formula, we need to assign some priors.

Looking carefully, there is some small difference in the intercepts of urge and mean depression. Their variances, however, are identical to M&H. I think this might have something to do with how MPlus / brms works and how the priors are specified, so I am not worried about that. Or maybe with how the lagged variable is treated.

Conclusion

Because it is this easy to specify latent means in brms, I think I will be using them much more often from now on, especially if my sample size per person is small. I don’t think this will make much of a difference after that sample size is greater than, say, the magic number 30.

Enders, Craig K., and Davood Tofighi. 2007. “Centering Predictor Variables in Cross-Sectional Multilevel Models: A New Look at an Old Issue.”Psychological Methods 12 (2): 121–38. https://doi.org/10.1037/1082-989X.12.2.121.

McNeish, Daniel, and Ellen L. Hamaker. 2020. “A Primer on Two-Level Dynamic Structural Equation Models for Intensive Longitudinal Data in Mplus.”Psychological Methods 25 (5): 610–35. https://doi.org/10.1037/met0000250.