I’ve got a model that I’m trying to figure out how to fit in `brms`

if possible. The general idea comes from a specific use case of a linear logistic test model from item response theory. To walk through this quickly, I’m going to attempt to write out the formulas (apologies if these don’t abide by appropriate standard notations, but hopefully they are at least accurate).

The linear logistic test model is usually treated as an extension of the one-parameter logistic regression, which is shown on the logit scale as:

\eta_{ij} = \theta_j + \beta_i

The log-odds, \eta, of person j answering item i correctly is a function of the person’s ability, \theta_j, and the item easiness, \beta_i. The linear logistic test model arises when the item easiness parameter is estimated from covariates such that,

\beta_i = \Sigma_{1}^{k}\omega_kX_k+\delta_i

Here, the item easiness parameter is just a function of various item covariates, X_k, that each have their own coefficients, \omega_k, and then each item can have its own residual estimated by a random intercept, \delta_i

My understanding is that the linear logistic test model would be specified in `brms`

using the following formula:

```
llt1pm <- brm(resp ~ cov1 + cov2 + cov3 + (1 | item) + (1 | id),
data = df, family = bernoulli(link = "logit"))
```

Extending this to the two-parameter logistic model is straightforward with `brms`

as well:

```
llt2pm <- brm(bf(resp ~ exp(logalpha)*theta + beta, nl = TRUE,
theta ~ 0 + (1 | id),
logalpha ~ 1 + cov1 + cov2 + cov3 + (1 | item),
beta ~ 1 + cov1 + cov2 + cov3 + (1 | item)),
data = df, family = bernoulli(link = "logit"))
```

I’m interested in something akin to differential item functioning but at the level of the item covariates. So, the equation for a one-parameter model would be something like this:

\eta_{ijg} = \theta_j + \omega_gG+ \omega_zZ_j + \delta_i +\omega_1X_1 + \omega_2X_2 + (\delta_j + \omega_gG + \omega_zZ_j)X_3

In effect, the item covariates (X_k) still operate as normal while there is still a permitted random residual associated with each item, \delta_i, but the coefficient or effect for one of the item covariates (X_3 in this case) is permitted to vary over persons. Specifically, it is of interest the extent to which differences in the effect of the item covariate can be separated into an effect due to a person-covariate, Z_j, a group-indicator, G, and then a residual person-level term, \delta_j. Since the person-level covariates are included in the item covariate effect estimation, these same parameters are included as part of the base model to account for potential between-person differences in \theta that may be due to the person-covariate and group indicator. Also, the group indicator is not a dummy variable but is instead coded as -1 and 1 since the emphasis is on interpreting the effect in terms of the difference between two groups.

I cannot seem to find a way to specify this model in `brms()`

. I’m more interested in the two-parameter extension of the formula as well. My first thought was that the specification should look something like this:

```
llt2pm <- brm(bf(resp ~ exp(logalpha)*theta + beta, nl = TRUE,
theta ~ z + g + (1 | id),
logalpha ~ 1 + cov1 + cov2 + ba*cov3 + (1 | item),
beta ~ 1 + cov1 + cov2 + bb*cov3 + (1 | item),
ba ~ z + g + (1 | id),
bb ~ z + g + (1 | id)),
data = df, family = bernoulli(link = "logit"))
```

This, of course, does not work because of the combined linear and non-linear components on the `beta`

and `logalpha`

parameters. I cannot, however, figure out how to properly specify the model so that `brms`

will run the model.

Does anyone know how to specify this model in `brms`

? If it’s not possible, then I’m also open to other thoughts about fitting the model.