Bias estimation in HSSM

[RaoXie]

Hello,

I just transit from HDDM to HSSM and I want to model the bias term 'z' in regression.

About the experiment, we have 3 difficulty level, and from the model-free analysis, we suspect there is some bias towards 1 stimulus across levels.

In HDDM, it can be achieved by link function, either via accuracy as 'response' variable and stimulus as regressor, or actual response as 'response' variable, as indicated in this tutorial

Is there any equivalent way of estimating bias in HSSM? I tried a self-defined link function, but doesn't seem to work. However, I tried this regression model (pic below and somehow it works, is this a correct way to estimate bias? Is there better way to do it in HSSM?

Any comments or suggestions would be appreciated! Thanks in advance!

Best, Rao

[AlexanderFengler]

Hey Rao,

If I read correctly, the tutorial link you shared refers to stimcoding in HDDM, which is not exactly what you would be doing in defining your regression model as shown below.

When you say it "works", could you share some info on traces and what values cond as well as stimulus can take?

Best,
Alex

[RaoXie]

Hello Alex,

Thank you very much for your reply. Basically I am wondering how to model bias term in HSSM when response is accuracy coded. In HDDM, I used link function, which is not supported in HSSM anymore.

Here is how I did it in HDDM, and I believe the z_link_func is important to make bias term stimulus-specific and sum up to 1:

n_subjects = 10
trials_per_level = 150 # and per stimulus
level1a = {'v':.3, 'a':2, 't':.3, 'sv':0, 'z':.5, 'sz':0, 'st':0}
level2a = {'v':.4, 'a':2, 't':.3, 'sv':0, 'z':.6, 'sz':0, 'st':0}
level3a = {'v':.5, 'a':2, 't':.3, 'sv':0, 'z':.7, 'sz':0, 'st':0}

data_a, params_a = hddm.generate.gen_rand_data({'level1': level1a,
                                                'level2': level2a,
                                                'level3': level3a},
                                                size=trials_per_level,
                                                subjs=n_subjects)

level1b = {'v':.5, 'a':2, 't':.3,'sv': 0, 'z':.5, 'sz': 0, 'st': 0}
level2b = {'v':.6, 'a':2, 't':.3,'sv': 0, 'z':.4, 'sz': 0, 'st': 0}
level3b = {'v':.6, 'a':2, 't':.3,'sv': 0, 'z':.3, 'sz': 0, 'st': 0}

data_b, params_b = hddm.generate.gen_rand_data({'level1': level1b,
                                                'level2': level2b,
                                                'level3': level3b},
                                                size=trials_per_level,
                                                subjs=n_subjects)

data_a['stimulus'] = Series(np.ones((len(data_a))), index=data_a.index)
data_b['stimulus'] = Series(np.ones((len(data_b)))*2, index=data_a.index)
mydata = data_a.append(data_b, ignore_index=True)

def z_link_func(x, data=mydata):
    stim = (np.asarray(dmatrix('0 + C(s, [[0], [1]])',
                              {'s': data.stimulus.loc[x.index]},return_type='dataframe'))
    )
    # Apply z = (1 - x) to flip them along 0.5
    z_flip = np.subtract(stim, x.to_frame())
    # The above inverts those values we do not want to flip,
    # so invert them back
    z_flip[stim == 0] *= -1
    return z_flip

z_reg = {'model': 'z ~ 1 + C(condition)', 'link_func': z_link_func}
v_reg = {'model': 'v ~ 1 + C(condition) + C(stimulus)', 'link_func': lambda x: x}
reg_descr = [z_reg, v_reg]

m_reg = hddm.HDDMRegressor(mydata, reg_descr, include=['v', 'a', 't', 'z'])
m_reg.sample(2000, burn=200)

In HSSM, I think response variable is by default accuracy coded (please correct me if I am wrong), here is how I am currently doing it, the parameters can be estimated relatively good (traces are attached below), but I am not sure if this is correct.

1. Now that 'z' in each condition and for each stimulus is estimated separately, how is 'z' guaranteed to sum up to 1, e.g., z_stim1_condition1 + z_stim2_condition1?
2. If this is completely nonsense, how should I build the model for 'z' in the current scenario?

n_subjects = 15  # number of subjects
n_trials = 150  # number of trials per subject

# Specify parameters
v_1a, a_1a, z_1a, t_1a = [0.3, 2, 0.5, 0.3]
v_2a, a_2a, z_2a, t_2a = [0.4, 2, 0.6, 0.3]
v_3a, a_3a, z_3a, t_3a = [0.5, 2, 0.7, 0.3]

v_1b, a_1b, z_1b, t_1b = [0.6, 2, 0.5, 0.3]
v_2b, a_2b, z_2b, t_2b = [0.8, 2, 0.4, 0.3]
v_3b, a_3b, z_3b, t_3b = [0.2, 2, 0.3, 0.3]

true_values_all= np.concatenate((np.repeat([[v_1a, a_1a, z_1a, t_1a]], axis=0, repeats=n_trials), 
                                np.repeat([[v_2a, a_2a, z_2a, t_2a]], axis=0, repeats=n_trials), 
                                np.repeat([[v_3a, a_3a, z_3a, t_3a]], axis=0, repeats=n_trials), 
                                np.repeat([[v_1b, a_1b, z_1b, t_1b]], axis=0, repeats=n_trials), 
                                np.repeat([[v_2b, a_2b, z_2b, t_2b]], axis=0, repeats=n_trials), 
                                np.repeat([[v_3b, a_3b, z_3b, t_3b]], axis=0, repeats=n_trials)), axis=0) 

cond = np.concatenate((np.ones(n_trials), np.ones(n_trials)*2, np.ones(n_trials)*3, np.ones(n_trials), np.ones(n_trials)*2, np.ones(n_trials)*3), axis=0)
stim = np.concatenate((np.ones(n_trials*3)*-1, np.ones(n_trials*3)), axis=0)

data_list = []
for i in range(n_subjects):
    # Simulate data
    sim_out = simulator(
        theta=true_values_all,  # parameter list
        model="ddm",  # specify model (many are included in ssms)
        n_samples=1,  # number of samples for each set of parameters
    )

    data_list.append(
        pd.DataFrame(
        np.column_stack([sim_out["rts"][:, 0], sim_out["choices"][:, 0], cond, stim, np.repeat(i, axis=0, repeats=len(sim_out["rts"][:, 0]))]),
        columns=["rt", "response", "condition", "stimulus", "subject"],
        )
    )
mydata = pd.concat(data_list)


m_reg = hssm.HSSM(
    data=mydata,
    model='ddm',
    include=[
        {
            "name": "v",
            "formula": "v ~ 1 + C(condition) + C(stimulus) + C(condition)*C(stimulus)",
            "link": "identity",
        },
        {
            "name": "z",
            "formula": "z ~ 1 + C(condition) + C(stimulus) + C(condition)*C(stimulus)",
            "link": "identity",
        }
    ],
)

infer_data_reg_v = m_reg.sample(
    sampler="mcmc", chains=1, cores=1, draws=2000, tune=2000, progressbar=True,
)

Thank you very much for your help and time!
Best, Rao