HDDM in HSSM - Model Fit, Link Function

[lauraglobig]

Hi,
I am currently exploring HSSM but would like to keep some of the existing HDDM features.

1. Is there a way to get model fit parameters from HDDM, such as BPIC, WAIC, DIC

2. I am trying to implement a link function for v but haven’t found an example of this. This is the link function for stimcoding I am currently using in HDDM - is there an equivalent for HSSM?

def v_link_func(x, data=data):
stim = (np.asarray(patsy.dmatrix('0 + C(s, [[1], [-1]])',
{'s': data.stim.iloc[x.index]})))
return x * stim
I noticed this came up as an issue before on this forum but couldn’t find an example.

3. Is there a way to keep the stats output as it was in HDDM, i.e., to also print CI?

Thank you!

[frankmj]

Hi Laura, HSSM uses WAIC and LOO via Arviz - both of which are better than DIC. You should be able to see an example of this in the main tutorial. But unfortunately there is currently no simple way to get DIC from HSSM because the Arviz developers decided it was not worth including. Re your first question on the hddm forum (not included here?) you can get individual subject posteriors if you use a hierarchical model (use the | subject notation - see the tutorial there too- but this might now be default all hierarchical) I will let other weigh in on issue 2, my memory is that there is link function examples in the tutorial under regression but not specifically for stim coding so this should probably be filed as an issue (for new upgrades/ features to hssm ) The plot_posteriors function in Arviz should give you a highest density interval (HDI), which is a type of CI. Other Arviz functions should as well Michael

On Tue, Aug 1, 2023 at 5:59 PM Laura K. Globig ***@***.***> wrote: Hi, I am currently exploring HSSM but would like to keep some of the existing HDDM features. 1. Is there a way to get model fit parameters from HDDM, such as BPIC, WAIC, DIC 2. I am trying to implement a link function for v but haven’t found an example of this. This is the link function for stimcoding I am currently using in HDDM - is there an equivalent for HSSM? def v_link_func(x, data=data): stim = (np.asarray(patsy.dmatrix('0 + C(s, [[1], [-1]])', {'s': data.stim.iloc[x.index]}))) return x * stim I noticed this came up as an issue before on this forum but couldn’t find an example. 1. Is there a way to keep the stats output as it was in HDDM, i.e., to also print CI? Thank you! — Reply to this email directly, view it on GitHub <#249>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAG7TFGQSCLKFFOHMGJZY3TXTF345ANCNFSM6AAAAAA3AMBLE4> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

-- Michael J Frank, PhD | Edgar L. Marston Professor Director, Carney Center for Computational Brain Science <https://www.brown.edu/carney/ccbs> Laboratory of Neural Computation and Cognition <https://www.lnccbrown.com/> Brown University website <http://ski.clps.brown.edu>

[lauraglobig]

Hi Michael
Thanks so much for the quick reply. That's super helpful!
I managed to solve the individual subject issue, so didn't include it here. An example of a stim-coding link function would be great. I was only able to find identity link functions, but may have missed it. Thanks again!
Best,
Laura

[lauraglobig]

Hi Michael,

Just following up on the link function issue for stimulus coding. I had a go at trying to adapt it but would love to get your thoughts on whether this is the correct approach. Apologies in advance for the long post - I tried to provide all the relevant information.

For reference: My model is as follows: 1) v varies as a function of intercept and trial type (0/1), 2) responses are coded as what trial type participants say they are in. 3) As such a positive z indicates a bias towards TT1, and similarly a positive vC_trialtype1 indicates a bias towards choosing trialtype1.

In HDDM I used to specify this as follows:

data = hddm.load_csv('./DDM.csv')
# create link function
def v_link_func(x, data=data):
    stim = (np.asarray(patsy.dmatrix('0 + C(s, [[1], [-1]])',
                                {'s': data.stim.iloc[x.index]})))
    return x * stim
# where stim is a way to indicate trial type where trial type=1 is stim=1 and trial type=0, is stim=2

# regression equation
v_reg = {'model': 'v ~ 1 + C(trialtype)', 'link_func': v_link_func}

# set up model
m_reg2 = hddm.HDDMRegressor(data, v_reg, informative=False, group_only_regressors=False, p_outlier=0.05, include=['z','sz'],  group_only_nodes=(['sz']))

From the BAMBI documentation I discovered that I need to specify my own link family and create v_link, v_linkinv and v_linkinv_backend (https://bambinos.github.io/bambi/_modules/bambi/families/link.html#Link), where the backend function needs to be something that works with PyMC backend. While v_link and v_linkinv worked fine with the old specification - the backend was giving me some issues, so I tried to convert the function as follows.

def v_link(x, data=data):
    stim = np.asarray("0 + C(s, [[1], [-1]])", {"s": data.stim.loc[x.index]})
    return x * stim

def v_linkinv(x, data=data):
    stim = np.asarray("0 + C(s, [[1], [-1]])", {"s": data.stim.loc[x.index]})
    return x / stim

# Define a custom function to compute stim based on data
def compute_stim(data):
    s_values = data.stim.values
    print (s_values)
    stim = np.zeros_like(s_values)
    stim[s_values == 1] = 1
    stim[s_values == 2] = -1
    print(stim)
    return stim

def v_linkinv_backend(x, data=data):
    stim = compute_stim(data)
    return x / stim

# specify model

stimcoded_ddm_model = hssm.HSSM(
    data=data,
    model="ddm",
    include=[
        {
            "name": "v",
            "formula": "v ~ 1 + C(trialtype) + (1+ C(trialtype)|participant_id)",
             "link": bm.families.Link(
                 "our_link",
                 link=v_link,
                 linkinv=v_linkinv,
                 linkinv_backend=v_linkinv_backend,
            ),
           #  "link": "identity",
        },
        {
            "name": "a",
            "formula": "a ~ 1 + (1|participant_id)",
            "link": "identity",
        },
        {
            "name": "t",
            "formula": "t ~ 1 + (1|participant_id)",
            "link": "identity",
        },
        {
            "name": "z",
            "formula": "z ~ 1 + (1|participant_id)",
            "link": "identity",
        },
    ],
)

My two questions are:

1. Is this the correct way to translate the link function ?
2. Is the implementation of the backend/linkinv appropriate here?

I'd appreciate any feedback or thoughts on this. Happy to provide sample simulated data too, if that's helpful?
Thanks for your help!

[frankmj]

Hi Laura, I looked into this a bit and don't think that link function will quite work (I believe it will return x for the entire data frame rather than only for each entry, my guess is that you will not recover parameters using that as it will end up averaging over the stims). The link function and inverse link functions are generally used to allow parameter transformations (eg a parameter that ranges from 0-1 might be sampled along a more continuous space and then converted back). For now, you can actually just use an identity link function and instead just recode the "stim" column to be +1 and -1 for each trial type the subject is in. (or you could code "trialtype" that way and replace stim with that below) Then in your regression just do (for a single participant first) : "formula": "v ~ 1 + stim", "link": "identity", now the fitted stim coefficient will give the drift rate that is positive on stim=1 trials and negative on stim=-1 trials. The intercept will tell you if there is an overall drift bias but one can also just do 0 + stim if you want to force the v's to be symmetrical as does the original stim coding model. Including z in the model will also allow the starting point bias to be estimated (I confirmed that both the drift and starting point biases are recovered in this set up). To make this hierarchical just do "formula": "v ~ 1 + 1|participant_id + stim|participant_id", Also more generally (for others) one can also include different conditions (e.g. if the difficulty in reporting trialtype is manipulated, like in the HDDMStimRegressor tutorial example <https://hddm.readthedocs.io/en/latest/tutorial_regression_stimcoding.html>). In that case the dataset would have an additional column (e.g. "cond"), and the model to fit would look like "formula": "v ~ 1 + stim:C(cond)", This will fit separate stim-coded drift rates for each condition and I also confirmed it recovers them and starting points from simulated data. Michael Michael J Frank, PhD | Edgar L. Marston Professor Director, Carney Center for Computational Brain Science <https://www.brown.edu/carney/ccbs> Laboratory of Neural Computation and Cognition <https://www.lnccbrown.com/> Brown University website <http://ski.clps.brown.edu>

…

On Sat, Aug 5, 2023 at 8:20 PM Laura K. Globig ***@***.***> wrote: Hi Michael, Just following up on the link function issue for stimulus coding. I had a go at trying to adapt it but would love to get your thoughts on whether this is the correct approach. Apologies in advance for the long post - I tried to provide all the relevant information. For reference: My model is as follows: 1) v varies as a function of intercept and trial type (0/1), 2) responses are coded as what trial type participants say they are in. 3) As such a positive z indicates a bias towards TT1, and similarly a positive vC_trialtype1 indicates a bias towards choosing trialtype1. In HDDM I used to specify this as follows: `# # load the data data = hddm.load_csv('./DDM.csv') create link function def v_link_func(x, data=data): stim = (np.asarray(patsy.dmatrix('0 + C(s, [[1], [-1]])', {'s': data.stim.iloc[x.index]}))) return x * stim where stim is a way to indicate trial type where trial type=1 is stim=1 and trial type=0, is stim=2 regression equation v_reg = {'model': 'v ~ 1 + C(trialtype)', 'link_func': v_link_func} set up model m_reg2 = hddm.HDDMRegressor(data, v_reg, informative=False, group_only_regressors=False, p_outlier=0.05, include=['z','sz'], group_only_nodes=(['sz']))` From the BAMBI documentation I discovered that I need to specify my own link family and create v_link, v_linkinv and v_linkinv_backend ( https://bambinos.github.io/bambi/_modules/bambi/families/link.html#Link), where the backend function needs to be something that works with PyMC backend. While v_link and v_linkinv worked fine with the old specification - the backend was giving me some issues, so I tried to convert the function as follows. ` import bambi as bm # needed for custom link def v_link(x, data=data): if data.stim.loc[x.index] == 1: return x * 1 else: return x * -1 def v_linkinv(x, data=data): if data.stim.loc[x.index] == 1: return 1 / (x * 1) else: return 1 / (x * -1) def v_linkinv_backend(x, data=data): return v_linkinv(x, data=data) specify model stimcoded_ddm_model = hssm.HSSM( data=data, model="ddm", include=[ { "name": "v", "formula": "v ~ 1 + C(trialtype) + (1+ C(trialtype)|participant_id)", "link": bm.families.Link( "our_link", link=v_link, linkinv=v_linkinv, linkinv_backend=v_linkinv_backend, ), # "link": "identity", }, { "name": "a", "formula": "a ~ 1 + (1|participant_id)", "link": "identity", }, { "name": "t", "formula": "t ~ 1 + (1|participant_id)", "link": "identity", }, { "name": "z", "formula": "z ~ 1 + (1|participant_id)", "link": "identity", }, ], ) ` My two questions are: 1. Is this the correct way to translate the link function ? 2. Is the implementation of the backend/linkinv appropriate here? I'd appreciate any feedback or thoughts on this. Happy to provide sample simulated data too, if that's helpful? Thanks for your help! — Reply to this email directly, view it on GitHub <#249 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAG7TFAQ7GFVAH6X5JKDDA3XT3PLHANCNFSM6AAAAAA3AMBLE4> . You are receiving this because you commented.Message ID: ***@***.***>

[lauraglobig]

Thanks so much! That is very helpful!