Regression models for 2-timepoint non-experimental data

\(\mathcal M_1 \colon\) The unconditional post model.

We can write the unconditional post model as

\[ \begin{align*} \text{post}_i & \sim \operatorname{Normal}(\mu, \sigma) \\ \mu & = \beta_0, \end{align*} \]

where \(\beta_0\) is both the model intercept and the estimate for the mean value of post. The focus this model places on post comes at the cost of any contextual information on what earlier values we might compare post to. Also, since the only variable in the model is post, this technically is not a 2-timepoint model. But given its connection to the models to follow, it’s worth working through.

Here’s how to fit the model with brms.

m1 <-
  brm(data = small_data_wide,
      post ~ 1,
      seed = 1,
      cores = 4, chains = 4, iter = 3500, warmup = 1000)

print(m1)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: post ~ 1 
##    Data: small_data_wide (Number of observations: 100) 
## Samples: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup samples = 10000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     1.02      0.12     0.79     1.25 1.00     6972     5683
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.17      0.08     1.02     1.36 1.00     8412     6599
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

We might compare those parameters with their sample values.

small_data_wide %>% 
  summarise(mean = mean(post),
            sd = sd(post))

## # A tibble: 1 x 2
##    mean    sd
##   <dbl> <dbl>
## 1  1.02  1.16

If you only care about computing the population estimates for \(\mu_\text{post}\) and \(\sigma_\text{post}\), this model does a great job. With no other variables in the model, this approach does a poor job telling us about growth processes.

\(\mathcal M_2 \colon\) The simple autoregressive model.

The simple model with the pre scores predicting post is a substantial improvement from the previous one. It follows the form

\[ \begin{align*} \text{post}_i & \sim \operatorname{Normal}(\mu_i, \sigma) \\ \mu_i & = \beta_0 + \beta_1 \text{pre}_i, \end{align*} \]

where the intercept \(\beta_0\) is the expected value of post when pre is at zero. As with many other regression contexts, centering the predictor pre at the mean or some other meaningful value can help make \(\beta_0\) more interpretable. Of greater interest is the \(\beta_1\) coefficient, which is the expected deviation from \(\beta_0\) for a one-unit increase in pre. But since pre and post are really the same variable y measured at two timepoints, it might be helpful if we express this model in another way. In perhaps more technical form, the simple model with pre predicting post is really an autoregressive model following the form

\[ \begin{align*} y_{ti} & \sim \operatorname{Normal}(\mu_i, \sigma) \\ \mu_i & = \beta_0 + \phi y_{t - 1,i}, \end{align*} \]

where the criterion \(y\) varies across persons \(i\) and timepoints \(t\). Here we only have two timepoints, \(\text{post} = t\) and \(\text{pre} = t - 1\). The strength of association between \(y_t\) and \(y_{t - 1}\) captured by the autoregressive parameter \(\phi\), which is often expressed in a correlation metric.

Here’s how to fit the model with brm().

m2 <-
  brm(data = small_data_wide,
      post ~ 1 + pre,
      seed = 1,
      cores = 4, chains = 4, iter = 3500, warmup = 1000)

print(m2)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: post ~ 1 + pre 
##    Data: small_data_wide (Number of observations: 100) 
## Samples: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup samples = 10000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     1.05      0.12     0.83     1.28 1.00     8995     6957
## pre           0.24      0.10     0.04     0.43 1.00     9548     6404
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.15      0.08     1.00     1.32 1.00     8407     7109
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Let’s compare the pre coefficient with the Pearson’s correlation between pre and post.

small_data_wide %>% 
  summarise(correlation = cor(pre, post))

## # A tibble: 1 x 1
##   correlation
##         <dbl>
## 1       0.232

Well look at that. Recall that the \(\beta_0\) parameter is the expected value in post when the predictor pre is at zero. Though the sample mean for pre is very close to zero, it’s not exactly so.

small_data_wide %>% 
  summarise(pre_mean = mean(pre))

## # A tibble: 1 x 1
##   pre_mean
##      <dbl>
## 1   -0.154

Here’s how to use that information to predict the mean value for post.

fixef(m2)[1, 1] + fixef(m2)[2, 1] * mean(small_data_wide$pre)

## [1] 1.015851

We can get a full posterior summary with aid from fitted().

nd <- tibble(pre = mean(small_data_wide$pre))

fitted(m2, newdata = nd)

##      Estimate Est.Error      Q2.5    Q97.5
## [1,] 1.015851 0.1142147 0.7939083 1.241721

\(\mathcal M_3 \colon\) The bivariate autoregressive model.

Though the simple autoregressive model gives us a sense of the strength of association between pre and post–and thus a sense of the stability in \(y\) over time–, it still lacks an explicit parameter for mean value of \(y\) at \(t - 1\). Enter the bivariate autoregressive model,

\[ \begin{align*} \text{post}_i & \sim \operatorname{Normal}(\mu_i, \sigma) \\ \text{pre}_i & \sim \operatorname{Normal}(\nu, \tau) \\ \mu_i & = \beta_0 + \beta_1 \text{pre}_i \\ \nu & = \gamma_0, \end{align*} \]

where \(\text{post}_i\) is still modeled as a simple linear function of \(\text{pre}_i\), but now we also include an unconditional model for \(\text{pre}_i\). This will give us an explicit comparison for where we started at the outset (\(\nu\)) and where we ended up (\(\mu_i\)). We can fit this model using the brms multivariate syntax where the two submodels are encased in bf() statements (Bürkner, 2021b). Also, be careful to use set_rescor(rescor = FALSE) to omit a residual correlation between the two. Their association is already handled with the \(\beta_1\) parameter.

m3 <-
  brm(data = small_data_wide,
      bf(post ~ 1 + pre) +
        bf(pre ~ 1) +
        set_rescor(rescor = FALSE),
      seed = 1,
      cores = 4, chains = 4, iter = 3500, warmup = 1000)

print(m3)

##  Family: MV(gaussian, gaussian) 
##   Links: mu = identity; sigma = identity
##          mu = identity; sigma = identity 
## Formula: post ~ 1 + pre 
##          pre ~ 1 
##    Data: small_data_wide (Number of observations: 100) 
## Samples: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup samples = 10000
## 
## Population-Level Effects: 
##                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## post_Intercept     1.05      0.12     0.82     1.28 1.00    10717     7470
## pre_Intercept     -0.15      0.11    -0.38     0.07 1.00    12822     7706
## post_pre           0.24      0.10     0.04     0.43 1.00    11030     7641
## 
## Family Specific Parameters: 
##            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_post     1.15      0.08     1.00     1.33 1.00    10955     7694
## sigma_pre      1.16      0.08     1.01     1.33 1.00    13273     7703
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Now that we’ve broken out the multivariate syntax, we might consider a second bivariate model.

\(\mathcal M_4 \colon\) The bivariate correlational model.

The bivariate correlational model follows the form

\[ \begin{align*} \begin{bmatrix} \text{post}_i \\ \text{pre}_i \end{bmatrix} & \sim \operatorname{MVNormal} \left (\begin{bmatrix} \mu \\ \nu \end{bmatrix}, \mathbf \Sigma \right) \\ \mu & = \beta_0 \\ \nu & = \gamma_0 \\ \mathbf \Sigma & = \mathbf{SRS} \\ \mathbf S & = \begin{bmatrix} \sigma & 0 \\ 0 & \tau \end{bmatrix} \\ \mathbf R & = \begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix}, \end{align*} \]

where means of both pre and post are modeled in intercept-only models. However, the association between the two timepoints is captured in the residual correlation \(\rho\). Yet because we have no predictors in for either variable, the “residual” correlation is really just a correlation. We might also gain some insights if we re-express this model in terms of \(y_{ti}\) and \(y_{t - 1,i}\):

\[ \begin{align*} \begin{bmatrix} y_{ti} \\ y_{t - 1,i} \end{bmatrix} & \sim \operatorname{MVNormal} \left (\begin{bmatrix} \mu_t \\ \mu_{t - 1} \end{bmatrix}, \mathbf \Sigma \right) \\ \mu_t & = \beta_t \\ \mu_{t - 1} & = \beta_{t - 1} \\ \mathbf \Sigma & = \mathbf{SRS} \\ \mathbf S & = \begin{bmatrix} \sigma_t & 0 \\ 0 & \sigma_{t - 1} \end{bmatrix} \\ \mathbf R & = \begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix}, \end{align*} \]

where what we formerly called an autoregressive coefficient in \(\mathcal M_2\), we’re now calling a correlation. Note also that this model freely estimates \(\sigma_t\) and \(\sigma_{t - 1}\). In some contexts, these are presumed to be equal. Though we won’t be imposing that constraint, here, I believe it is possible with the brms non-linear syntax (Bürkner, 2021c). Anyway, here’s how to fit the model with the brm() function .

m4 <-
  brm(data = small_data_wide,
      bf(post ~ 1) +
        bf(pre ~ 1) +
        set_rescor(rescor = TRUE),
      seed = 1,
      cores = 4, chains = 4, iter = 3500, warmup = 1000)

Note out use of the set_rescor(rescor = TRUE) syntax in the model formula. This explicitly told brm() to include the residual correlation. Here’s the summary.

print(m4)

##  Family: MV(gaussian, gaussian) 
##   Links: mu = identity; sigma = identity
##          mu = identity; sigma = identity 
## Formula: post ~ 1 
##          pre ~ 1 
##    Data: small_data_wide (Number of observations: 100) 
## Samples: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup samples = 10000
## 
## Population-Level Effects: 
##                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## post_Intercept     1.02      0.12     0.79     1.24 1.00    10992     8033
## pre_Intercept     -0.15      0.12    -0.39     0.08 1.00    10679     8012
## 
## Family Specific Parameters: 
##            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_post     1.18      0.08     1.03     1.36 1.00    11342     7819
## sigma_pre      1.16      0.08     1.02     1.34 1.00    11150     7777
## 
## Residual Correlations: 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## rescor(post,pre)     0.23      0.09     0.03     0.40 1.00    10074     7292
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Now the intercept and sigma parameters do a good job capturing the sample statistics.

small_data_wide %>% 
  pivot_longer(pre:post) %>% 
  group_by(name) %>% 
  summarise(mean = mean(value),
            sd = sd(value)) %>% 
  mutate_if(is.double, round, digits = 2)

## # A tibble: 2 x 3
##   name   mean    sd
##   <chr> <dbl> <dbl>
## 1 post   1.02  1.16
## 2 pre   -0.15  1.14

The new ‘rescor’ line at the bottom of the print() summary approximates the Pearson’s correlation of the two variables, much like the autoregressive parameter did two models up.

small_data_wide %>% 
  summarise(correlation = cor(pre, post))

## # A tibble: 1 x 1
##   correlation
##         <dbl>
## 1       0.232

Another nice quality of this model is if you subtract \(\gamma_0\) from \(\beta_0\) (i.e., \(\beta_t - \beta_{t - 1}\)), you’d end up with the posterior mean of the change score.

posterior_samples(m4) %>% 
  mutate(change = b_post_Intercept - b_pre_Intercept) %>% 
  summarise(mu_change = mean(change))

##   mu_change
## 1  1.168898

Keeping that in mind, let’s switch gears to the first of the change-score models.

\(\mathcal M_5 \colon\) The unconditional change-score model.

We’ve already saved that in our small_data_wide data as change. Here’s what the unconditional change-score model² looks like:

\[ \begin{align*} \text{change}_i & \sim \operatorname{Normal}(\mu, \sigma) \\ \mu & = \beta_0, \end{align*} \]

where \(\beta_0\) is the expected value for \(\text{change}_i\). In the terms of the last model, \(\text{change}_i = \text{post}_i - \text{pre}_i\) or, in the terms of the simple autoregressive model, \(y_{\Delta i} = y_{ti} - y_{t - 1,i}\).

m5 <-
  brm(data = small_data_wide,
      change ~ 1,
      seed = 1,
      cores = 4, chains = 4, iter = 3500, warmup = 1000)

print(m5)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: change ~ 1 
##    Data: small_data_wide (Number of observations: 100) 
## Samples: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup samples = 10000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     1.17      0.14     0.89     1.45 1.00     8075     6240
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.44      0.10     1.26     1.66 1.00     8673     6829
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Thus, this model suggests the average change from pre to post was about 1.2 units. We can compute the sample statistics for that in two ways.

small_data_wide %>% 
  summarise(change = mean(change),
            `post - pre` = mean(post - pre))

## # A tibble: 1 x 2
##   change `post - pre`
##    <dbl>        <dbl>
## 1   1.17         1.17

The major deficit in the unconditional change model is that change is disconnected from any reference points. We have no explicit way of knowing what number we changed from or what number we changed to. The next model offers a solution.

\(\mathcal M_6 \colon\) The conditional change-score model.

Instead of fitting an unconditional model of change, why not condition on the initial pre value? We might express this as

\[ \begin{align*} \text{change}_i & \sim \operatorname{Normal}(\mu_i, \sigma) \\ \mu_i & = \beta_0 + \beta_1 \text{pre}_i, \end{align*} \]

where \(\beta_0\) is now the expected value for change when pre is at zero. As with the simple autoregressive model, centering the predictor pre at the mean or some other meaningful value can help make \(\beta_0\) more interpretable. Perhaps of greater interest, the \(\beta_1\) coefficient allows us to predict different levels of change, conditional in the initial values at pre.

m6 <-
  brm(data = small_data_wide,
      change ~ 1 + pre,
      seed = 1,
      cores = 4, chains = 4, iter = 3500, warmup = 1000)

print(m6)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: change ~ 1 + pre 
##    Data: small_data_wide (Number of observations: 100) 
## Samples: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup samples = 10000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     1.05      0.12     0.82     1.28 1.00     9315     7108
## pre          -0.76      0.10    -0.96    -0.57 1.00    10559     7804
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.15      0.08     1.00     1.33 1.00     9745     7268
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Since the intercept in this model is the expected change value based on when pre == 0, it might be easiest to interpret that value when using the mean of pre.

fixef(m6)[1, 1] + fixef(m6)[2, 1] * mean(small_data_wide$pre)

## [1] 1.167458

That is the point prediction for the mean of change. Let’s compare that to the sample value.

small_data_wide %>% 
  summarise(mean = mean(change))

## # A tibble: 1 x 1
##    mean
##   <dbl>
## 1  1.17

Of course, we can get a fuller summary using the fitted() method.

nd <- tibble(pre = mean(small_data_wide$pre))

fitted(m6, newdata = nd)

##      Estimate Est.Error      Q2.5   Q97.5
## [1,] 1.167458 0.1151598 0.9428394 1.39075

Note how the coefficient for pre is about -0.76. This is a rough analogue of the negative correlation among the intercepts and slopes in the original data-generating model.

rho

## [1] -0.5

It tells us something very similar; participants with higher values at pre tended to have lower change values. Much like with the simple autoregressive model, a deficit of this model is there is no explicit parameter for the expected value of pre, which we can amend by fitting a bivariate model.

\(\mathcal M_7 \colon\) The bivariate conditional change-score model.

The bivariate conditional change-score model follows the form

\[ \begin{align*} \text{change}_i & \sim \operatorname{Normal}(\mu_i, \sigma) \\ \text{pre}_i & \sim \operatorname{Normal}(\nu, \tau) \\ \mu_i & = \beta_0 + \beta_1 \text{pre}_i \\ \nu & = \gamma_0, \end{align*} \]

where the simple linear model of \(\text{change}_i\) conditional on \(\text{pre}_i\) is coupled with an unconditional intercept-only model for \(\text{pre}_i\). We can fit this model with brms by way of the multivariate syntax, where the two submodels are encased in bf() statements and we set set_rescor(rescor = FALSE) to omit a residual correlation between the two.

m7 <-
  brm(data = small_data_wide,
      bf(change ~ 1 + pre) +
        bf(pre ~ 1) +
        set_rescor(rescor = FALSE),
      seed = 1,
      cores = 4, chains = 4, iter = 3500, warmup = 1000)

print(m7)

##  Family: MV(gaussian, gaussian) 
##   Links: mu = identity; sigma = identity
##          mu = identity; sigma = identity 
## Formula: change ~ 1 + pre 
##          pre ~ 1 
##    Data: small_data_wide (Number of observations: 100) 
## Samples: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup samples = 10000
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## change_Intercept     1.05      0.12     0.82     1.28 1.00    11908     7575
## pre_Intercept       -0.16      0.12    -0.39     0.07 1.00    13113     8030
## change_pre          -0.76      0.10    -0.96    -0.57 1.00    12023     7190
## 
## Family Specific Parameters: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_change     1.15      0.08     1.00     1.33 1.00    13064     7696
## sigma_pre        1.16      0.08     1.01     1.33 1.00    11419     7142
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Now we have an intercept for pre, we can use the model parameters to compute the expected values for pre, change, and post.

posterior_samples(m7) %>% 
  mutate(pre    = b_pre_Intercept,
         change = b_change_Intercept + b_change_pre * b_pre_Intercept) %>% 
  mutate(post = pre + change) %>% 
  pivot_longer(pre:post) %>% 
  group_by(name) %>% 
  mean_qi(value) %>% 
  mutate_if(is.double, round, digits = 2)

## # A tibble: 3 x 7
##   name   value .lower .upper .width .point .interval
##   <chr>  <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 change  1.17   0.88   1.45   0.95 mean   qi       
## 2 post    1.01   0.78   1.25   0.95 mean   qi       
## 3 pre    -0.16  -0.39   0.07   0.95 mean   qi

The posterior means are in the value column and the lower- and upper-levels of the percentile-based 95% intervals are in the .lower and .upper columns. Now compare those mean estimates with the sample means.

small_data_wide %>% 
  pivot_longer(-id) %>% 
  group_by(name) %>% 
  summarise(mean = mean(value)) %>% 
  mutate_if(is.double, round, digits = 2)

## # A tibble: 3 x 2
##   name    mean
##   <chr>  <dbl>
## 1 change  1.17
## 2 post    1.02
## 3 pre    -0.15

As handy as this model is, the \(\beta_1\) coefficient might not be in the most intuitive metric. Let’s reparameterize.

\(\mathcal M_8 \colon\) The bivariate correlational pre/change model.

The bivariate correlational pre/change model follows the form

\[ \begin{align*} \begin{bmatrix} \text{change}_i \\ \text{pre}_i \end{bmatrix} & \sim \operatorname{MVNormal} \left (\begin{bmatrix} \mu \\ \nu \end{bmatrix}, \mathbf \Sigma \right) \\ \mu & = \beta_0 \\ \nu & = \gamma_0 \\ \mathbf \Sigma & = \mathbf{SRS} \\ \mathbf S & = \begin{bmatrix} \sigma & 0 \\ 0 & \tau \end{bmatrix} \\ \mathbf R & = \begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix}, \end{align*} \]

where means of both pre and change are modeled in intercept-only models and the association between the two is captured by the correlation \(\rho\). And again, because we have no predictors in for either variable, the “residual” correlation is really just a correlation. Here’s how to fit the model with brms.

m8 <-
  brm(data = small_data_wide,
      bf(change ~ 1) +
        bf(pre ~ 1) +
        set_rescor(rescor = TRUE),
      seed = 1,
      cores = 4, chains = 4, iter = 3500, warmup = 1000)

print(m8)

##  Family: MV(gaussian, gaussian) 
##   Links: mu = identity; sigma = identity
##          mu = identity; sigma = identity 
## Formula: change ~ 1 
##          pre ~ 1 
##    Data: small_data_wide (Number of observations: 100) 
## Samples: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup samples = 10000
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## change_Intercept     1.17      0.15     0.88     1.45 1.00     8200     7551
## pre_Intercept       -0.15      0.12    -0.38     0.08 1.00     8160     7194
## 
## Family Specific Parameters: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_change     1.44      0.10     1.26     1.67 1.00     8811     7147
## sigma_pre        1.16      0.08     1.01     1.33 1.00     8996     7487
## 
## Residual Correlations: 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## rescor(change,pre)    -0.60      0.06    -0.71    -0.46 1.00     8299     7233
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Now the parameter in the ‘Residual Correlations’ section of the summary output is a close analogue to our original data-generating rho parameter.

rho

## [1] -0.5

Those with higher pre values tended to have lower change values. We can look at that with a plot of the data.

p2 <-
  small_data_wide %>% 
  ggplot(aes(x = pre, y = change)) +
  geom_point() +
  stat_ellipse(color = "grey50")

(p1 + p2) &
  coord_cartesian(xlim = range(small_data_wide$pre),
                  ylim = range(small_data_wide$change))

Here we’ve placed the scatter plot of the data-generating id-level intercepts and slopes next to the scatter plot of the pre and change scores. They are not exactly the same, but the latter are a partial consequence of the former. This is why the correlation parameter in our m8 model closely, but not exactly, resembled the data-generating rho parameter.

Okay, let’s switch gears again.

\(\mathcal M_9 \colon\) The grand-mean model.

The simplest model we might fit using the long version of the 2-timepoint data, small_data_long, is what we might call the grand-mean model, or what Hoffman (2015) called the between-person empty model. It follows the form

\[ \begin{align*} y_{ti} & \sim \operatorname{Normal}(\mu, \sigma) \\ \mu & = \beta_0, \end{align*} \]

where the intercept \(\beta_0\) is the expected value value for \(y\) across all \(i\) participants and \(t\) timepoints. We fit this with brm() much like we fit the unconditional post model and the unconditional change-score model.

m9 <-
  brm(data = small_data_long,
      y ~ 1 ,
      seed = 1,
      cores = 4, chains = 4, iter = 3500, warmup = 1000)

print(m9)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: y ~ 1 
##    Data: small_data_long (Number of observations: 200) 
## Samples: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup samples = 10000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.43      0.09     0.25     0.61 1.00     8486     6236
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.30      0.07     1.18     1.43 1.00     9455     6868
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

We might check these with the sample statistics for y.

small_data_long %>% 
  summarise(mean = mean(y),
            sd = sd(y))

## # A tibble: 1 x 2
##    mean    sd
##   <dbl> <dbl>
## 1 0.431  1.29

Though this model did to a good job describing the population values for the mean and standard deviation for y, it did a terrible job telling us about change in y, about individual differences in that change, or about anything else of interest we might have as longitudinal researchers.

\(\mathcal M_{10} \colon\) The random-intercept model.

Now that we have a grand mean, we might want to ask what kinds of variables would help explain the variation around the grand mean. From a multilevel perspective, the first source of variation of interest will be across participants, which we can express by allowing the mean to vary by participant. This is what Hoffman called the within-person empty model and the empty means, random intercept model. It follows the form

\[ \begin{align*} y_{ti} & \sim \operatorname{Normal}(\mu_i, \sigma) \\ \mu_i & = \beta_0 + u_{\text{id},i} \\ u_\text{id} & \sim \operatorname{Normal}(0, \sigma_\text{id}) \\ \sigma & = \sigma_\epsilon, \end{align*} \]

where \(\beta_0\) is now the grand mean among the participant-level means in \(y_{ti}\). The participant-level deviations from the grand mean are expressed as \(u_{\text{id},i}\), which is normally distributed with a mean at zero (these are deviations, after all) and a standard deviation of \(\sigma_\text{id}\). The \(\sigma_\epsilon\) parameter is a mixture of the variation within participants and over time. With brms, we can fit the random-intercept model model like so.

m10 <-
  brm(data = small_data_long,
      y ~ 1 + (1 | id),
      seed = 1,
      cores = 4, chains = 4, iter = 3500, warmup = 1000,
      control = list(adapt_delta = .99))

print(m10)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: y ~ 1 + (1 | id) 
##    Data: small_data_long (Number of observations: 200) 
## Samples: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup samples = 10000
## 
## Group-Level Effects: 
## ~id (Number of levels: 100) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.23      0.15     0.01     0.56 1.00     2079     4188
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.43      0.09     0.24     0.61 1.00    13669     6976
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.27      0.07     1.14     1.41 1.00     6211     5523
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Our \(\beta_0\) intercept returned a very similar estimate for the grand mean, but we now interpret it as the grand mean for the within-person means for y. The variation in the id-level deviations around the grand mean, which we called \(u_{\text{id},i}\), is summarized by \(\sigma_\text{id}\) in the ‘Group-Level Effects’ section of the summary output.

Authors of many longitudinal text books (e.g., Hoffman, 2015; Singer & Willett, 2003) typically present this model as a way to directly compare the between- and within-person variation in the data by way of the intraclass correlation coefficient (ICC),

\[\text{ICC} = \frac{\text{between-person variance}}{\text{total variance}} = \frac{\sigma_\text{id}^2}{\sigma_\text{id}^2 + \sigma_\epsilon^2}.\]

When using frequentist methods, the ICC is typically expressed with a point estimate. When working with all our posterior draws, we can get full posterior distributions within the Bayesian framework.

posterior_samples(m10) %>% 
  mutate(icc = sd_id__Intercept^2 / (sd_id__Intercept^2 + sigma^2)) %>% 
  ggplot(aes(x = icc, y = 0)) +
  stat_halfeye(.width = .95) +
  scale_x_continuous("Intraclass correlation coefficient (ICC)", 
                     breaks = 0:5 / 5, expand = c(0, 0), limits = 0:1) +
  scale_y_continuous(NULL, breaks = NULL)

The ICC is a proportion, which limits it to the range of zero to one. Here it suggests that 0–20% of the variation in our data is due to differences between participants; the remaining variation occurs within them. Given that our data were collected across time, it might make sense to fit a model that explicitly accounts for time.

\(\mathcal M_{11} \colon\) The cross-classified model.

A direct extension of the random-intercept model is the cross-classified multilevel model (McElreath, 2015, Chapter 12), which we might express as

\[ \begin{align*} y_{ti} & \sim \operatorname{Normal}(\mu_{ti}, \sigma) \\ \mu_{ti} & = \beta_0 + u_{\text{id},i} + u_{\text{time},i} \\ u_\text{id} & \sim \operatorname{Normal}(0, \sigma_\text{id}) \\ u_\text{time} & \sim \operatorname{Normal}(0, \sigma_\text{time}) \\ \sigma & = \sigma_\epsilon, \end{align*} \]

where \(\sigma_\text{id}\) captures systemic differences between participants, \(\sigma_\text{time}\) captures systemic variation across the two timepoints, and \(\sigma_\epsilon\) captures the variation within participants over time. Another way to think of this model is as a Bayesian multilevel version of the repeated-measures ANOVA³, where variance is partitioned into a between level (\(\sigma_\text{id} \approx \text{SS}_\text{between}\)), a model level (\(\sigma_\text{time} \approx \text{SS}_\text{model}\)), and error (\(\sigma_\epsilon \approx \text{SS}_\text{error}\)).

m11 <-
  brm(data = small_data_long,
      y ~ 1 + (1 | time) + (1 | id),
      seed = 1,
      cores = 4, chains = 4, iter = 3500, warmup = 1000,
      control = list(adapt_delta = .9999,
                     max_treedepth = 11))

print(m11)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: y ~ 1 + (1 | time) + (1 | id) 
##    Data: small_data_long (Number of observations: 200) 
## Samples: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup samples = 10000
## 
## Group-Level Effects: 
## ~id (Number of levels: 100) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.51      0.16     0.14     0.78 1.00     1519     1545
## 
## ~time (Number of levels: 2) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.67      1.30     0.41     5.00 1.00     5240     7037
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.46      1.07    -1.78     2.81 1.00     4541     4935
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.04      0.08     0.90     1.20 1.00     2313     3750
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

The intercept is the grand mean across all measures of y. The first row in the ‘Group-Level Effects’ section is our summary for \(\sigma_\text{id}\), which gives us a sense of the variation between participants in their overall tendencies in the criterion y. If we use the posterior_samples(), we can even look at the posteriors for the \(u_{\text{id},i}\) parameters, themselves.

posterior_samples(m11) %>% 
  pivot_longer(starts_with("r_id")) %>% 
  
  ggplot(aes(x = value, y = reorder(name, value))) +
  stat_pointinterval(point_interval = mean_qi, .width = .95, size = 1/6) +
  scale_y_discrete(expression(italic(i)), breaks = NULL) +
  labs(subtitle = expression(sigma[id]~is~the~summary~of~the~variation~across~these),
       x = expression(italic(u)[id][','*italic(i)]))

Given that each of the \(u_{\text{id},i}\) parameters is based primarily on two data points (the two data points per participant), it should be no surprise they are fairly wide. Even a few more measurement occasions within participants will narrow them substantially. If you fit the same model using the original 6-timepoint data, you’ll see the 95% intervals are almost half as wide.

Perhaps of greater interest are the \(u_{\text{time},i}\) parameters. If you combine them with the intercept, you’ll get the model-based expected values at both timepoints.

posterior_samples(m11) %>% 
  transmute(pre  = b_Intercept + `r_time[0,Intercept]`,
            post = b_Intercept + `r_time[1,Intercept]`) %>% 
  pivot_longer(everything()) %>% 
  group_by(name) %>% 
  mean_qi(value) %>% 
  mutate_if(is.double, round, digits = 2)

## # A tibble: 2 x 7
##   name  value .lower .upper .width .point .interval
##   <chr> <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 post   1.01   0.78   1.24   0.95 mean   qi       
## 2 pre   -0.15  -0.38   0.09   0.95 mean   qi

The final variance parameter, \(\sigma_\text{time}\), captures the within-participant variation over time. With a model like this, it seems natural to directly compare the magnitudes of the three variance parameters, which answers the question: Where’s the variance at? Here we’ll do so with a plot.

posterior_samples(m11) %>% 
  select(sd_id__Intercept:sigma) %>% 
  set_names("sigma[id]", "sigma[time]", "sigma[epsilon]") %>% 
  pivot_longer(everything()) %>% 
  mutate(name = factor(name, levels = c("sigma[epsilon]", "sigma[time]", "sigma[id]"))) %>%
  
  ggplot(aes(x = value, y = name)) +
  tidybayes::stat_halfeye(.width = .95, size = 1, normalize = "xy") +
  scale_x_continuous("marginal posterior", expand = expansion(mult = c(0, 0.05)), breaks = c(0, 1, 2, 5)) +
  scale_y_discrete(NULL, labels = ggplot2:::parse_safe) +
  coord_cartesian(xlim = c(0, 5.25),
                  ylim = c(1.5, 3.5)) +
  theme(axis.text.y = element_text(hjust = 0))

At first glance, it might be surprising how wide the \(\sigma_\text{time}\) posterior is compared to the other two. Yet recall this parameter is summarizing the standard deviation of only two levels. If you have experience with multilevel models, you’ll know that it can be difficult to estimate a variance parameter with few levels–two levels is the extreme lower limit. This is why we had to fiddle with the adapt_delta and max_treedepth parameters within the brm() function to get the model to sample properly. Though we pulled this off using default priors, don’t be surprised if you have to use tighter priors when fitting a model like this.

\(\mathcal M_{12} \colon\) The simple liner model.

The last three models focused on the grand mean and sources of variance around that grand mean. A more familiar looking approach might be to fit a simple linear model with y conditional on time,

\[ \begin{align*} y_{ti} & \sim \operatorname{Normal}(\mu_{ti}, \sigma) \\ \mu_{ti} & = \beta_0 + \beta_1 \text{time}_{ti}, \end{align*} \]

where the intercept \(\beta_0\) is the expected value at the first timepoint and \(\beta_1\) captures the change in \(y_{ti}\) for the final timepoint.

m12 <-
  brm(data = small_data_long,
      y ~ 1 + time,
      seed = 1,
      cores = 4, chains = 4, iter = 3500, warmup = 1000)

print(m12)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: y ~ 1 + time 
##    Data: small_data_long (Number of observations: 200) 
## Samples: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup samples = 10000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.15      0.12    -0.38     0.08 1.00     9459     7344
## time          1.17      0.16     0.85     1.49 1.00     9582     7294
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.16      0.06     1.05     1.28 1.00     9982     7362
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Now our \(\beta_0\) and \(\beta_1\) parameters are the direct pre/post single-level analogues to the population-level \(\beta_0\) and \(\beta_1\) parameters from our original multilevel model based on the full 6-timepoint data set.

fixef(m0) %>% round(digits = 2)

##           Estimate Est.Error  Q2.5 Q97.5
## Intercept    -0.13      0.11 -0.33  0.08
## time          1.16      0.13  0.91  1.42

The major deficit in this model, which is the reason you’ll see it criticized in the methodological literature, is it ignores how the y values are nested within levels of id. The only variance parameter, \(\sigma\), was estimated under the typical assumption that the residuals are all independent of one another. Sure, the model formula accounted for the overall trend in time, but it ignored the insights revealed from many of the other models the capture between-participant correlations in intercepts and slopes. This means that if you know something about the value of one’s residual for when time == 0, you’ll also know something about where to expect their residual for when time == 1. The two are not independent. As long as we’re working with the data in the long format, we’ll want to account for this, somehow.

\(\mathcal M_{13} \colon\) The liner model with a random intercept.

A natural first step to accounting for how the y values are nested within levels of id is to fit a random-intercept model, or what Hoffman (2015) called the fixed linear time, random intercept model. It follows the form

\[ \begin{align*} y_{ti} & \sim \operatorname{Normal}(\mu_{ti}, \sigma) \\ \mu_{ti} & = \beta_{0i} + \beta_1 \text{time}_{ti} \\ \beta_{0i} & = \gamma_0 + u_{0i} \\ u_{0i} & \sim \operatorname{Normal}(0, \sigma_0) \\ \sigma & = \sigma_\epsilon, \end{align*} \]

where \(\beta_{0i}\) is the intercept and \(\beta_1\) is the time slope. Although this parameterization holds the variation in slopes constant across participants, between-participant variation is at least captured in \(\beta_{0i}\), which is decomposed into a grand mean, \(\gamma_0\), and participant-level deviations around that grand mean, \(u_{0i}\). Those participant-level deviations are summarized by the \(\sigma_0\) parameter. In this model, \(\sigma_\epsilon\) is a mixture of within-participant variation and between-participant variation in slopes.

m13 <-
  brm(data = small_data_long,
      y ~ 1 + time + (1 | id),
      seed = 1,
      cores = 4, chains = 4, iter = 3500, warmup = 1000)

print(m13)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: y ~ 1 + time + (1 | id) 
##    Data: small_data_long (Number of observations: 200) 
## Samples: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup samples = 10000
## 
## Group-Level Effects: 
## ~id (Number of levels: 100) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.52      0.15     0.14     0.78 1.00     1645     1960
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.15      0.12    -0.38     0.08 1.00     9441     7715
## time          1.17      0.15     0.88     1.46 1.00    14821     6993
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.03      0.08     0.90     1.19 1.00     2403     4360
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

The \(\beta\) parameters are very similar to those from the simple linear model, above.

fixef(m12) %>% round(digits = 2)

##           Estimate Est.Error  Q2.5 Q97.5
## Intercept    -0.15      0.12 -0.38  0.08
## time          1.17      0.16  0.85  1.49

But now look at the size of \(\sigma_0\), which suggests substantial differences in staring points. To get a sense of what this means, we’ll plot all 100 participant-level trajectories with a little help from fitted().

nd <- distinct(small_data_long, id, time)

fitted(m13, 
       newdata = nd) %>% 
  data.frame() %>% 
  bind_cols(nd) %>% 
  
  ggplot(aes(x = time, y = Estimate, group = id)) +
  geom_abline(intercept = fixef(m13)[1, 1],
              slope = fixef(m13)[2, 1],
              size = 3, color = "blue") +
  geom_line(size = 1/4, alpha = 2/3) +
  scale_x_continuous(breaks = 0:1) +
  labs(subtitle = "Random intercepts, fixed slope",
       y = "y") +
  coord_cartesian(ylim = c(-1, 2))

The bold blue line in the middle is based on the population-level intercept and slope, whereas the thinner black lines are the participant-level trajectories. To keep from overplotting, we’re only showing the posterior means, here. Because we only allowed the intercept to vary across participants, all the slopes are identical. And indeed, look at all the variation we see in the intercepts–an insight lacking in the simple linear model.

\(\mathcal M_{14} \colon\) The liner model with a random slope.

The counterpoint to the last model is to allow the time slopes, but not the intercepts, vary across participants:

\[ \begin{align*} y_{ti} & \sim \operatorname{Normal}(\mu_{ti}, \sigma) \\ \mu_{ti} & = \beta_0 + \beta_{1i} \text{time}_{ti} \\ \beta_{1i} & = \gamma_1 + u_{1i} \\ u_{1i} & \sim \operatorname{Normal}(0, \sigma_1) \\ \sigma & = \sigma_\epsilon, \end{align*} \]

where \(\beta_0\) is the intercept for all participants. Now \(\beta_{1i}\) is the population mean for the distribution of slopes, which vary across participants, the standard deviation for which is measured by \(\sigma_1\).

m14 <-
  brm(data = small_data_long,
      y ~ 1 + time + (0 + time | id),
      seed = 1,
      cores = 4, chains = 4, iter = 3500, warmup = 1000,
      control = list(adapt_delta = .95))

print(m14)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: y ~ 1 + time + (0 + time | id) 
##    Data: small_data_long (Number of observations: 200) 
## Samples: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup samples = 10000
## 
## Group-Level Effects: 
## ~id (Number of levels: 100) 
##          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(time)     0.33      0.21     0.02     0.75 1.00     2240     4641
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.15      0.11    -0.38     0.08 1.00    20513     7195
## time          1.17      0.16     0.85     1.49 1.00    18079     7133
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.13      0.06     1.01     1.26 1.00     6825     6149
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Here’s random-slopes alternative to the random-intercepts plot from the last section.

fitted(m14, 
       newdata = nd) %>% 
  data.frame() %>% 
  bind_cols(nd) %>% 
  
  ggplot(aes(x = time, y = Estimate, group = id)) +
  geom_abline(intercept = fixef(m14)[1, 1],
              slope = fixef(m14)[2, 1],
              size = 3, color = "blue") +
  geom_line(size = 1/4, alpha = 2/3) +
  scale_x_continuous(breaks = 0:1) +
  labs(subtitle = "Fixed intercept, random slopes",
       y = "y") +
  coord_cartesian(ylim = c(-1, 2))

But why choose between random intercepts or random slopes?

\(\mathcal M_{15} \colon\) The multilevel growth model with regularizing priors.

If you were modeling 2-timepoint data with conventional frequentist estimators (e.g., maximum likelihood), you can have random intercepts or random slopes, but you can’t have both; that would require data from three timepoints or more. But because Bayesian models bring in extra information by way of the priors, you can actually fit a full multilevel growth model with both random intercepts and slopes:

\[ \begin{align*} y_{ti} & \sim \operatorname{Normal}(\mu_{ti}, \sigma_\epsilon ) \\ \mu_{ti} & = \beta_0 + \beta_1 \text{time}_{ti} + u_{0i} + u_{1i} \text{time}_{ti} \\ \begin{bmatrix} u_{0i} \\ u_{1i} \end{bmatrix} & \sim \operatorname{MVNormal} \left (\begin{bmatrix} 0 \\ 0 \end{bmatrix}, \mathbf \Sigma \right) \\ \mathbf \Sigma & = \mathbf{SRS} \\ \mathbf S & = \begin{bmatrix} \sigma_0 & 0 \\ 0 & \sigma_1 \end{bmatrix} \\ \mathbf R & = \begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix}. \end{align*} \]

The trick is you have to go beyond the diffuse brms default settings for the priors for \(\sigma_0\) and \(\sigma_1\). If you have high-quality information from theory or previous studies, you can base the priors on those. Another approach is to use regularizing priors. Given standardized data, members of the Stan team like either \(\operatorname{Normal}^+(0, 1)\) or \(\operatorname{Student-t}^+(3, 0, 1)\) for variance parameters (see the Generic prior for anything section from the Prior choice recommendations wiki). In the second edition of his text, McElreath generally favored the \(\operatorname{Exponential}(1)\) prior (McElreath, 2020), which is the approach we’ll experiment with, here. It’ll also help if we use a regularizing prior on the correlation among the intercepts and slopes, \(\rho\).

m15 <-
  brm(data = small_data_long,
      y ~ 1 + time + (1 + time | id),
      seed = 1,
      prior = prior(exponential(1), class = sd) +
        prior(lkj(4), class = cor),
      cores = 4, chains = 4, iter = 3500, warmup = 1000,
      control = list(adapt_delta = .9995))

Even with our tighter priors, we still had to adjust the adapt_delta parameter to improve the quality of the MCMC sampling. Take a look at the model summary.

print(m15)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: y ~ 1 + time + (1 + time | id) 
##    Data: small_data_long (Number of observations: 200) 
## Samples: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup samples = 10000
## 
## Group-Level Effects: 
## ~id (Number of levels: 100) 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)           0.52      0.18     0.12     0.84 1.01     1012     1473
## sd(time)                0.31      0.23     0.01     0.88 1.00     1021      876
## cor(Intercept,time)    -0.05      0.33    -0.63     0.59 1.00     4039     6272
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.15      0.11    -0.38     0.07 1.00    10022     8348
## time          1.17      0.15     0.88     1.46 1.00    13870     7175
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.01      0.10     0.80     1.18 1.00     1060      946
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Though the posteriors, particularly for the \(\sigma\) and \(\rho\) parameters, are not as precise as with the 6-timepoint data, we now have a model with a summary mirroring the structure of the data-generating model. Yet compared to the data-generating values, the estimates for \(\sigma_1\) and \(\rho\) are particularly biased toward zero. Here’s a look at the trajectories.

fitted(m15,
       newdata = nd) %>% 
  data.frame() %>% 
  bind_cols(nd) %>% 
  
  ggplot(aes(x = time, y = Estimate, group = id)) +
  geom_abline(intercept = fixef(m15)[1, 1],
              slope = fixef(m15)[2, 1],
              size = 3, color = "blue") +
  geom_line(size = 1/4, alpha = 2/3) +
  scale_x_continuous(breaks = 0:1) +
  labs(subtitle = "Random intercepts AND random slopes\n(2-timepoint data)",
       y = "y") +
  coord_cartesian(ylim = c(-1, 2))

For comparision, here’s the plot for the original 6-timepoint model.

fitted(m0,
       newdata = nd) %>% 
  data.frame() %>% 
  bind_cols(nd) %>% 
  
  ggplot(aes(x = time, y = Estimate, group = id)) +
  geom_abline(intercept = fixef(m0)[1, 1],
              slope = fixef(m0)[2, 1],
              size = 3, color = "blue") +
  geom_line(size = 1/4, alpha = 2/3) +
  scale_x_continuous(breaks = 0:1) +
  labs(subtitle = "Random intercepts AND random slopes\n(6-timepoint data)",
       y = "y") +
  coord_cartesian(ylim = c(-1, 2))

There wasn’t enough information in the 2-timepoint data set to capture the complexity in the full 6-timepoint data set.

\(\mathcal M_{16} \colon\) The fixed effects with correlated error model.

Though our Bayesian 2-timepoint version of the full multilevel growth model was exciting, it’s not generally used in the wild. Even with our tighter regularizing priors, there just wasn’t enough information in the data to do the model justice. A very different and humbler approach is to combine the simple linear model with the autoregressive model,

\[ \begin{align*} y_{ti} & \sim \operatorname{Normal}(\mu_{ti}, \mathbf \Sigma) \\ \mu_{ti} & = \beta_0 + \beta_1 \text{time}_{ti} \\ \mathbf \Sigma & = \mathbf{SRS} \\ \mathbf S & = \begin{bmatrix} \sigma & 0 \\ 0 & \sigma \end{bmatrix} \\ \mathbf R & = \begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix}, \end{align*} \]

where \(\rho\) captures the correlation between the responses in the two timepoints, \(t\) and \(t - 1\), which is an alternative to the way the mixed model from above handles the dependencies (a.k.a heteroskedasticity) inherent in longitudinal data. Note how \(\mathbf \Sigma\) in this model is defined very differently from the full multilevel growth model from above. To fit this model with brms, we use the ar() syntax.

m16 <-
  brm(data = small_data_long,
      y ~ time + ar(time = time, p = 1, gr = id),
      seed = 1,
      cores = 4, chains = 4, iter = 3500, warmup = 1000)

print(m16)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: y ~ time + ar(time = time, p = 1, gr = id) 
##    Data: small_data_long (Number of observations: 200) 
## Samples: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup samples = 10000
## 
## Correlation Structures:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## ar[1]     0.23      0.10     0.04     0.43 1.00     9766     7458
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.15      0.12    -0.38     0.07 1.00    10593     7136
## time          1.17      0.15     0.89     1.46 1.00     9607     6546
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.15      0.06     1.04     1.27 1.00    10604     7665
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

This summary suggests that, after you account for the linear trend, the correlation between \(y_{ti}\) and \(y_{t - 1,i}\) is about 0.23. Though we don’t get id-specific variance parameters, this model does account for the nonindependence of the data over time. If you scroll back up, notice how similar this is to the correlation from the bivariate correlational model, m4.

\(\mathcal M_{17} \colon\) The cross-classified model with robust variances for discrete time.

One of the criticisms of the conventional repeated measures ANOVA approach is how it presumes the variances in \(y\) are constant over time. However, we can relax that constraint with a model like

\[ \begin{align*} y_{ti} & \sim \operatorname{Normal}(\mu_{ti}, \sigma_{ti}) \\ \mu_{ti} & = \beta_0 + u_{\text{id},i} + u_{\text{time},t} \\ \log (\sigma_{ti}) & = \eta_{\text{time},t} \\ u_\text{id} & \sim \operatorname{Normal}(0, \sigma_\text{id}) \\ u_\text{time} & \sim \operatorname{Normal}(0, \sigma_\text{time}), \end{align*} \]

where we are now modeling both parameters in the likelihood, \(\mu\) AND \(\sigma\). The second line shows a typical-looking model for \(\mu_{ti}\). All the excitement lies in the third line, which contains the linear model for \(\log (\sigma_{ti})\). The reason we are modeling \(\log (\sigma_{ti})\) rather than directly modeling \(\sigma\) is to avoid solutions that predict negative values for \(\sigma\). For this model, it’s unlikely we’d run into that problem. But since the brms default is to use the log link anytime we model \(\sigma\) within the distributional modeling syntax, we’ll just get used to the log link right from the start. If you are unfamiliar with link functions, they’re widely used within the generalized linear modeling framework. Logistic regression with the logit link and Poisson regression with the log link are two widely-used examples. For more on link functions and the generalized linear model, check out the texts by Agresti (2015); Gelman, Hill, and Vehtari (2020); and McElreath (2020, 2015).

Anyway, the linear model for \(\mu_{ti}\) is exactly the same as with the original cross-classified model, m11; it includes a grand mean (\(\beta_0\)) and two kinds of deviations around that grand mean (\(u_{\text{id},i}\) and \(u_{\text{time},t}\)). The model for \(\log (\sigma_{ti})\) contains an intercept, which varies across the two levels of time, \(\eta_{\text{time},t}\). Here’s how to fit the model with brms::brm().

m17 <-
  brm(data = small_data_long,
      bf(y ~ 1 + (1 | time) + (1 | id),
         sigma ~ 0 + factor(time)),
      seed = 1,
      cores = 4, chains = 4, iter = 3500, warmup = 1000,
      control = list(adapt_delta = .999,
                     max_treedepth = 12))

print(m17)

##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: y ~ 1 + (1 | time) + (1 | id) 
##          sigma ~ 0 + factor(time)
##    Data: small_data_long (Number of observations: 200) 
## Samples: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup samples = 10000
## 
## Group-Level Effects: 
## ~id (Number of levels: 100) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.53      0.15     0.15     0.79 1.00     1376     1194
## 
## ~time (Number of levels: 2) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.66      1.26     0.41     5.00 1.00     5133     6854
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept             0.45      1.07    -1.78     2.66 1.00     4653     5056
## sigma_factortime0     0.01      0.10    -0.18     0.20 1.00     3068     5577
## sigma_factortime1     0.03      0.10    -0.17     0.22 1.00     2678     4488
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Even though we didn’t explicitly ask to use the log link in our brm() syntax, you can look at the second line in the print() output to see that it was automatically used. Though I won’t explore how to do so, here, one can fit this model without the log link. Anyway, the primary focus in this model is the sigma_factortime0 and sigma_factortime1 lines in the ‘Population-Level Effects’ section of the output. Those are the summaries for the \(\sigma\) parameters, conditional on whether time == 0 or time == 1. Though is might be difficult to evaluate parameters on the log scale, we can simply exponentiate them to convert them back to their natural metric.

fixef(m17)[2:3, c(1, 3:4)] %>% exp()

##                   Estimate      Q2.5    Q97.5
## sigma_factortime0 1.013928 0.8352872 1.223336
## sigma_factortime1 1.032290 0.8446723 1.247435

In this case, it looks like the two parameters are largely overlapping. If we work directly with the posterior draws, we can compute a formal difference score and plot the results.

posterior_samples(m17) %>% 
  mutate(`sigma[time==0]` = exp(b_sigma_factortime0),
         `sigma[time==1]` = exp(b_sigma_factortime1),
         `sigma[time==1]-sigma[time==0]` = exp(b_sigma_factortime1) - exp(b_sigma_factortime0)) %>% 
  pivot_longer(`sigma[time==0]`:`sigma[time==1]-sigma[time==0]`) %>% 
  
  ggplot(aes(x = value, y = name)) + 
  stat_halfeye(.width = .95) +
  scale_y_discrete(NULL, labels = ggplot2:::parse_safe) +
  coord_cartesian(ylim = c(1.5, 3.1)) +
  xlab("marginal posterior") +
  theme(axis.text.y = element_text(hjust = 0))

In this case, it looks like there was little difference between \(\sigma_{\text{time} = 0}\) and \(\sigma_{\text{time} = 1}\). This shouldn’t be a surprise; we simulated that data that way. However, it won’t always be like this in real-world data.

\(\mathcal M_{18} \colon\) The simple liner model with robust variance for linear time.

Our cross-classified approach treated time as a factor. Here we’ll treat it as a continuous variable in the models of both \(\mu\) and \(\sigma\). This will be a straight extension of the simple linear model, m12, following the form

\[ \begin{align*} y_{ti} & \sim \operatorname{Normal}(\mu_{ti}, \sigma_{ti}) \\ \mu_{ti} & = \beta_0 + \beta_1 \text{time}_{ti} \\ \log(\sigma_{ti}) & = \eta_0 + \eta_1 \text{time}_{ti}, \end{align*} \]

where the intercept \(\beta_0\) is the expected value at the first timepoint and \(\beta_1\) captures the change in \(y_i\) for the final timepoint. Now \(\log(\sigma_{ti})\) has a similar linear model, where \(\eta_0\) is the expected log of the standard deviation at the first timepoint and \(\eta_1\) captures the change standard deviation for the final timepoint.

m18 <-
  brm(data = small_data_long,
      bf(y ~ 1 + time, 
         sigma ~ 1 + time),
      seed = 1,
      cores = 4, chains = 4, iter = 3500, warmup = 1000)

print(m18)

##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: y ~ 1 + time 
##          sigma ~ 1 + time
##    Data: small_data_long (Number of observations: 200) 
## Samples: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup samples = 10000
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept          -0.15      0.11    -0.38     0.06 1.00    11620     6673
## sigma_Intercept     0.14      0.07     0.00     0.28 1.00    11090     7442
## time                1.17      0.16     0.85     1.49 1.00    11137     7818
## sigma_time          0.01      0.10    -0.18     0.21 1.00    11157     7620
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Even though this model looks very different from the last one, we can wrangle the posterior draws a little to make a similar plot comparing \(\sigma\) at the two timepoints.

posterior_samples(m18) %>% 
  mutate(`sigma[time==0]` = exp(b_sigma_Intercept),
         `sigma[time==1]` = exp(b_sigma_Intercept + b_sigma_time * 1),
         `sigma[time==1]-sigma[time==0]` = exp(b_sigma_Intercept) - exp(b_sigma_Intercept + b_sigma_time)) %>% 
  pivot_longer(`sigma[time==0]`:`sigma[time==1]-sigma[time==0]`) %>% 
  
  ggplot(aes(x = value, y = name)) + 
  stat_halfeye(.width = .95) +
  scale_y_discrete(NULL, labels = ggplot2:::parse_safe) +
  coord_cartesian(ylim = c(1.5, 3.1)) +
  xlab("marginal posterior") +
  theme(axis.text.y = element_text(hjust = 0))

Though this model is robust to differences in \(\sigma\) based on timepoint, it still ignores systemic differences across participants. The next model tackles that that limitation in spades.

\(\mathcal M_{19} \colon\) The liner model with correlated random intercepts for \(\mu\) and \(\sigma.\)

Here we return to the multilevel model framework to accommodate participant-level differences for both \(\mu\) and \(\sigma\):

\[ \begin{align*} y_{ti} & \sim \operatorname{Normal}(\mu_{ti}, \sigma_{ti}) \\ \mu_{ti} & = \beta_0 + \beta_1 \text{time}_{ti} + u_{0i} \\ \log(\sigma_{ti}) & = \eta_0 + u_{2i} \\ \begin{bmatrix} u_{0i} \\ u_{2i} \end{bmatrix} & \sim \operatorname{MVNormal} \left (\begin{bmatrix} 0 \\ 0 \end{bmatrix}, \mathbf \Sigma \right) \\ \mathbf \Sigma & = \mathbf{SRS} \\ \mathbf S & = \begin{bmatrix} \sigma_0 & 0 \\ 0 & \sigma_2 \end{bmatrix} \\ \mathbf R & = \begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix}, \end{align*} \]

where \(\beta_0\) is the grand mean for the intercepts and \(u_{0i}\) are the participant-level deviations around that grand mean. \(\beta_1\) is the time slope, which is invariant across all participants in this model. The \(\eta_0\) parameter is the grand mean for the log standard deviations and \(u_{2i}\) captures the participant-level deviations around that grand mean. In the fourth line, we learn that \(u_{0i}\) and \(u_{2i}\) are multivariate normal, with a mean vector of two zeros and a variance/covariance matrix \(\mathbf \Sigma\). As is typical within the brms framework, we decompose \(\mathbf \Sigma\) into a variance matrix \(\mathbf S\) and correlation matrix \(\mathbf R\). Of particular interest is the \(\rho\) parameter, which captures the correlation in the participant-level intercepts and participant-level standard deviations. Here’s how to fit the model with brms.

m19 <-
  brm(data = small_data_long,
      bf(y ~ 1 + time + (1 |x| id),
         sigma ~ 1 + (1 |x| id)),
      seed = 1,
      cores = 4, chains = 4, iter = 3500, warmup = 1000,
      control = list(adapt_delta = .9))

You may have noticed the |x| parts in the formula lines for y and sigma. What that did was tell brms we wanted those parameters to be correlated. That is, that’s how we estimated the \(\rho\) parameter. There was nothing special about including x between the vertical lines. We could have used any other character. The important thing is that we used the same character in both. Anyway, here’s the model summary.

print(m19)

##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: y ~ 1 + time + (1 | x | id) 
##          sigma ~ 1 + (1 | x | id)
##    Data: small_data_long (Number of observations: 200) 
## Samples: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup samples = 10000
## 
## Group-Level Effects: 
## ~id (Number of levels: 100) 
##                                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)                      0.52      0.16     0.14     0.78 1.00     1496     1443
## sd(sigma_Intercept)                0.11      0.08     0.01     0.29 1.00     3032     4589
## cor(Intercept,sigma_Intercept)    -0.34      0.49    -0.97     0.80 1.00     5949     6097
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept          -0.15      0.12    -0.38     0.08 1.00    10102     7923
## sigma_Intercept     0.02      0.08    -0.13     0.17 1.00     2344     3931
## time                1.17      0.15     0.88     1.45 1.00    14868     8009
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Once again, brms used the log link for \(\sigma\). If you want to see \(\eta_0\) in its natural \(\sigma\) metric, exponentiate.

fixef(m19)["sigma_Intercept", c(1, 3:4)] %>% exp()

##  Estimate      Q2.5     Q97.5 
## 1.0171275 0.8759026 1.1849572

The sd(sigma_Intercept) row in the ‘Group-Level Effects’ section shows the variation in those \(\log \sigma\)’s. It might be easier to appreciate them in a plot.

posterior_samples(m19) %>% 
  pivot_longer(starts_with("r_id__sigma")) %>% 
  mutate(sigma_i = exp(sd_id__sigma_Intercept + value)) %>% 
  
  ggplot(aes(x = sigma_i, y = reorder(name, sigma_i))) +
  stat_pointinterval(point_interval = mean_qi, .width = .95, size = 1/6) +
  scale_y_discrete(expression(italic(i)), breaks = NULL) +
  labs(subtitle = expression(sigma[2]~summarizes~the~variation~across~these),
       x = expression(sigma[italic(i)]))

In this case, there was not a lot of variation in \(\sigma\) across participants. This is because we simulated the data that way. Though it may be hard to model participant-level variances with 2-timepoint data, I have found it comes in handy in real-world data sets based on more measurement occasions.

Finally, it might be useful to consider our \(\rho\) parameter, which suggested a mild negative correlation between the \(u_{0i}\) and \(u_{2i}\) deviations. Here’s how you might visualize that in a plot.

bind_rows(
  posterior_samples(m19) %>% select(starts_with("r_id[")) %>% set_names(1:100),
  posterior_samples(m19) %>% select(starts_with("r_id__sigma") %>% set_names(1:100))
) %>% 
  mutate(iter = rep(1:c(n() / 2), times = 2),
         type = rep(c("intercept", "log_sigma"), each = n() / 2)) %>% 
  pivot_longer(-c(iter, type)) %>% 
  pivot_wider(names_from = type, values_from = value) %>% 
  
  ggplot(aes(x = intercept, y = log_sigma, group = name)) +
  stat_ellipse(geom = "polygon", level = .01, alpha = 1/4) +
  labs(x = expression(italic(u)[0][italic(i)]),
       y = expression(log(italic(u)[2][italic(i)])))

Each of the ovals is a 1% ellipse of the bivariate posterior for \(u_{0i}\) and \(\log(u_{2i})\). Notice how using ellipses helps reveal the differences in the between- and within-person patterns.

This approach where residual variance parameters vary across participants has its origins in the work of Donald Hedeker and colleagues (Hedeker et al., 2008, 2012). More recently, Philippe Rast and colleagues (particularly graduate student, Donald Williams) have adapted this approach for use within the Stan/brms ecosystem (Williams, Liu, et al., 2019; Williams, Rouder, et al., 2019; Williams, Martin, et al., 2019; Williams, Zimprich, et al., 2019).

\(\mathcal M_{20} \colon\) The liner model with a random slope for \(\mu\) and uncorrelated random intercept for \(\sigma\).

Though we can find interesting things when we allow the random components in the \(\mu\) and \(\sigma\) models, we don’t have to think of them as covarying. Here we fit an extension of the linear model with a random time slope, where we add an orthogonal random intercept for \(\sigma\):

\[ \begin{align*} y_{ti} & \sim \operatorname{Normal}(\mu_{ti}, \sigma_{ti}) \\ \mu_{ti} & = \beta_0 + \beta_1 \text{time}_{ti} + u_{1i} \\ \log(\sigma_{ti}) & = \eta_0 + u_{2i} \\ u_{1i} & \sim \operatorname{Normal}(0, \sigma_1) \\ u_{2i} & \sim \operatorname{Normal}(0, \sigma_2), \end{align*} \]

where the two random components, \(u_{1i}\) and \(u_{2i}\), are now modeled with separate normal distributions, \(\operatorname{Normal}(0, \sigma_1)\) and \(\operatorname{Normal}(0, \sigma_2)\).

m20 <-
  brm(data = small_data_long,
      bf(y ~ 1 + time + (0 + time | id),
         sigma ~ 1 + (1 | id)),
      seed = 1,
      cores = 4, chains = 4, iter = 3500, warmup = 1000)

Note that in sharp contrast with our syntax for the previous model, this time we did not employ the |x| syntax in the formula lines for y and sigma. By omitting the |x| syntax, we omitted the correlation among those two random effects. Here’s the summary.

print(m20)

##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: y ~ 1 + time + (0 + time | id) 
##          sigma ~ 1 + (1 | id)
##    Data: small_data_long (Number of observations: 200) 
## Samples: 4 chains, each with iter = 3500; warmup = 1000; thin = 1;
##          total post-warmup samples = 10000
## 
## Group-Level Effects: 
## ~id (Number of levels: 100) 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(time)                0.33      0.21     0.02     0.75 1.00     2031     4346
## sd(sigma_Intercept)     0.12      0.09     0.00     0.32 1.00     2463     3519
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept          -0.15      0.11    -0.37     0.07 1.00    13451     7317
## sigma_Intercept     0.11      0.06    -0.02     0.23 1.00     4612     4412
## time                1.18      0.16     0.86     1.49 1.00    12600     7054
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

We can see the \(\sigma_1\) and \(\sigma_2\) summaries in the ‘Group-Level Effects’ section. It’s hard to compare them directly, because one is based on parameters in the log metric. But we can at least get a sense of what these parameters are summarizing by plotting the bivariate posterior for \(u_{1i}\) and \(\log(u_{2i})\).

bind_rows(
  posterior_samples(m20) %>% select(starts_with("r_id[")) %>% set_names(1:100),
  posterior_samples(m20) %>% select(starts_with("r_id__sigma") %>% set_names(1:100))
) %>% 
  mutate(iter = rep(1:c(n() / 2), times = 2),
         type = rep(c("slope", "log_sigma"), each = n() / 2)) %>% 
  pivot_longer(-c(iter, type)) %>% 
  pivot_wider(names_from = type, values_from = value) %>% 
  
  ggplot(aes(x = slope, y = log_sigma, group = name)) +
  stat_ellipse(geom = "polygon", level = .01, alpha = 1/4) +
  labs(x = expression(italic(u)[1][italic(i)]),
       y = expression(log(italic(u)[2][italic(i)])))

Notice how, this time, the 1% ellipses suggest no clear association between these two dimensions.