Item response theory model with smoothing splines in brms

Modeling
1

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

y_{ijk} \sim \text{Binom}(10, \mu_{ijk})

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

\nu_{ijk} = x_{ijk}'\beta + (\lambda'z_{ij}) \eta_{ij}

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

\eta_{ij} = s(a_{ij}) + \zeta_{ij}^{(2)} + \zeta_{i}^{(3)}

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

\nu = x'\beta + \lambda (s(a) + \zeta^{(2)} + \zeta^{(3)})

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?

2

Turned out it wasn’t too hard. Following the VerbAgg example in the JSS paper referred to in my original question, this did the job:

formula_va_2pl <- bf(
  r2 ~ exp(logalpha) * (ability + difficulty),
  ability ~ 1 + s(age, k = 4, bs = "cr") + (1 | id),
  difficulty ~ 1 + (1 |i| item),
  logalpha ~ 1 + (1 |i| item),
  nl = TRUE
)
1 Like