I have data where subjects have been measured by multiple cognitive tests at multiple timepoints. The response variables are thus binomially distributed, and the hierarchical structure of the data is of the form subject-timepoint within subject. The hypothesis is that the tests at each timepoint model the same underlying latent trait, and I want to investigate how other parameters interact with the age trajectory of this latent trait.

The model I want to fit is part of the GLLAMM framework defined by Rabe-Hesketh and Skrondal. A simplified version, that captures the essential elements, has response distribution

where \mu_{ijk} = \text{logit}^{-1}(\nu_{ijk}). The index i captures the subject and the index j captures the timepoint (within-subject) and index k captures the $k$th measurement within subject and timepoint. The linear predictor \nu_{ijk} is defined according to the measurement model

where x_{ijk} are covariates affect the measurement (e.g., retest effects, item difficulty). z_{ijk} is a dummy variable indicating which item $(i,j,k)$th measurement is measuring. Assuming K items, we have \lambda \in \mathbb{R}^{K}. The \eta_{ij} is the person-specific ability. Hence, the term (\lambda'z_{ij}) \eta_{ij} will simplify to the familiar \lambda_{k} \eta_{ijk} for each item k. Finally, the structural model is

where \zeta_{ij}^{(2)} \sim N(0, \psi^{(2)}) is a random intercept term for each measurement occasion (across items) and \zeta_{i}^{(3)} \sim N(0, \psi^{(3)}) is a random intercept for each participant (across subjects). Finally, and this is the critical part, s(a_{ij}) is a smooth term describing how the latent ability varies with age. To final goal is to see how s(a_{ij}) interacts with other variables, but that is not relevant for this question.

The reduced form model, and droping some indices for readability, is that the linear predictor is defined according to

In the `brms`

paper on item response theory it is well described how one can use nonlinear functions, but here I would like to use GAM-style regression splines, and I wonder if it’s possible to set this up using `brms`

, or if I have to code it directly in Stan.

If I assume a parallel measurement model, i.e., \lambda = 1 for all items, I can set it up using the following

```
mod <- brm(
formula = bf(
recall | trials(10) ~ test + mo(times_tested) + s(age) + (1 | id / timepoint)
),
data = dat,
family = binomial(),
chains = 4,
iter = 2000,
cores = 4,
save_pars = save_pars(group = FALSE)
)
```

Here `test`

is a factor variable describing which item is being measured, `mo(times_tested)`

is a monotonic term aimed to capture retest effects, and `s(age)`

models how the latent variable varies with age. However, I know that the items have different discriminatory ability, so I would like to estimate a factor loading which describes how `s(age) + (1 | id / timepoint)`

impacts the outcome differentially, depending on the value of `test`

.

Would it be possible to set up the desired model using `brms`

, or am I better off starting with the `brms`

generated code and then modifying it in Stan to get what I want?