Consider the following repeated measures model:

for $i = 1, \ldots, n$, $j = 1, 2$ where $n$ is the sample size, $j$ represents the index of the repeated measure, i.e., each subject has two measurements, $\mu_i$ is a normally distributed random effect, $\varepsilon_{ij}$ is a normally distributed error term, $y_{ij}$ is the continuous response, and $a_{ij}, b_{ij}$ are covariates. This is a multilevel model because of the nested structure of the data, and also non-linear in the $\beta_1$ parameter. In this post I simulate some data under this model, and try to leverage Bayesian computation techniques to estimate the parameters using the brms which is an interface to fit Bayesian generalized (non-)linear multilevel models using Stan.

We first load the required packages:

## Simulate Data (n = 500, two timepoints)

Next I simulate some multilevel data:

## Fit Non-linear Multilevel Bayesian Model

We first need to specify priors for $\beta_1$ and the random effect $\mu_i$. I have found the easiest way to do this, is to first get information for which priors may be specified using the brms::get_prior function:

The nl = TRUE argument indicates that this should be treated as a non-linear model. We need to specify priors for the population level effects as well as the standard deviation of the random effect (what is referred to as group-level effects in the output):

We can check the corresponding STAN code to verify that the prior information has been updated accordingly using the brms::make_stancode function:

Next we fit the model with 6 Markov chains:

We can check the summary to see if the model was able to accurately estimate $\beta_1$:

The estimate of interest in the output above is beta_Intercept under Population-Level Effects. This corresponds to the estimate of $\beta_1$ in the equation above. We see that the estimate is close to the true value of 3. We also see that the estimate of the standard deviation of the random effect is 2.13 [95% CI: 0.28, 3.42], indicating a strong subject specific effect (which is what we would expect since we generated the data this way).

We can also plot some model diagnostics to ensure proper mixing of the Markov chains: