Custom likelihood functions with carry-over information in TuringGLM, usecase: behavioral modeling

[PTWaade]

I am a behavioral modeler, in the sense that I specify probabilistic input-output rules that are used to describe behavioral data. These models then have parameters, which are estimated separately for participants and compared across groups. This is usually done in STAN, JAGS, JAX, or PyMC, and I do it in Turing.
An example of a model would be a simple Rescorla-Wagner, which learns about its environment with some learning rate.
Often, parameters are estimated with Turing, and then used in a GLM to find group differences or correlations with psychiatric symptoms etc. It is of course optimal to include this GLM in the same model as the behavioural model. To do this, however, the GLM has to be specified in the same Turing model as the behavioral model. This both takes time, is inaccessible to researchers not used to writing their own GLMs in Turing or similar, and does not take advantage of the modularity inherent in this kind of work. Often, the same GLM can be used with different behavioral models, and behavioral models can also be used outside of GLMs to for example do agent-based simulation (and because of Julia's solution the 2-language problem, this can be done with the same code function implementing the behavioral model).

To do this, it would be very powerful if a behavioral model can be specified directly in a GLM formula, for example using TuringGLM. A syntax could, tentatively, look something like this:

actions ~ behavioural_model(observations, behavioural_parameter)
behavioural_parameter ~ 1 + demographic_parameter + (1|ID)
etc.

Note that most behavioral models carry over information from observation to observation, for example as they learn about their environment, and that this needs to be accounted for.
This would allow researchers to specify state-of-the art hierarchical behavioural models and their linear models in one model, without needing expertise in programming in Turing, and without spending time on implementing the hierarchical structure anew every time; it makes this approach accessible to researchers used to working with classic linear regressions; and it allows for modular use of behavioural models either as part of GLMs or on their own.

How hard would this be to implement? Would it fit well within TuringGLM, or should it be a separate package? Does something similar already exist, that I have not encountered?

[torfjelde]

What we could quite easily make happen is something like:

using Turing, TuringGLM

@model function behavioural_model(observations, behavioural_parameter)
    # ...
end

# Define the model for `behavioural_parameter` using TuringGLM.jl.
# `const` to make it efficient
const behavioural_parameter_model = turing_model(
    @formula(behavioural_parameter ~ 1 + demographic_parameter),
    ...
)

@model function full_model(observations)
    # ...
    @submodel behavioural_parameter = behavioural_parameter_model(demographic_parameter)
    actions ~ behavioural_model(observations, behavioural_parameter)
end

where we're just making use of the @submodel macro to embed one Turing model into another (this is already possible).

Would this be similar to what you're looking for?

I think the syntax that you want of being able to do

behavioural_parameter ~ 1 + demographic_parameter + (1|ID)

directly inside a @model is somewhat non-trivial unfortunately, at least with the current TuringGLM.jl implementation.

It's definitively feasible to implement something like

@formula behavioural_parameter ~ 1 + demographic_parameter + (1|ID)

that would achieve this, but my instinct tells me it wouldn't be trivial.