<jupyter_start><jupyter_text>load data  load RTs data<jupyter_code>data = hddm.load_csv('data/study1.csv')

data = hddm.utils.flip_errors(data)
data.head(10)
# drop subject 10 due to data convergence issue
data = data.loc[~(data.subj_idx == 10),]<jupyter_output><empty_output><jupyter_text>Confirmatory Model fitting - HDDM model<jupyter_code># fit the HDDM model with drift rate parameter
# sample 10000 with 2000 burn-in samples
# save the model
m = hddm.HDDM(data, depends_on={'v': 'stim'}, p_outlier=.05)
m.find_starting_values()
m.sample(10000, burn=2000, dbname='traces_hddm.db', db='pickle')
m.save('traces_hddm')<jupyter_output><empty_output><jupyter_text>Test correlation with measured mood<jupyter_code># load data for measured mood for each participants
mood = pd.read_csv("data/study1_measured_mood.csv")
# data wrangling
subj_index_0V1A = [index for index in list(m.nodes_db.index) if "v_subj(0V1A)" in index]
subj_index_1V0A = [index for index in list(m.nodes_db.index) if "v_subj(1V0A)" in index]
node_0V1A = m.nodes_db.loc[subj_index_0V1A, 'node']
node_1V0A = m.nodes_db.loc[subj_index_1V0A, 'node']
v_0V1A_subj = []
v_1V0A_subj = []
for i in range(len(node_0V1A)):
    v_0V1A_subj.append(node_0V1A[i].trace().mean())
    v_1V0A_subj.append(node_1V0A[i].trace().mean())
mood["v_0V1A"] = v_0V1A_subj
mood["v_1V0A"] = v_1V0A_subj
mood
import pymc3 as pm
# test the correlation between drift rate and measured mood 
# using pymc3 package
# test valence

with pm.Model() as model_valence:
    alpha = pm.Normal("alpha", mu=0, sigma=10)
    beta_valence = pm.Normal("beta", mu=0, sigma=10)
    sigma = pm.HalfNormal("sigma", sigma=1)
    
    mu_valence = alpha + beta_valence * mood["movie_valence_measured"].values
    
    v_1V0A = pm.Normal('v_1V0A', mu=mu_valence, sd=sigma, observed=mood['v_1V0A'].values)
    
    trace_v = pm.sample(10000, tune=2000)
# test the correlation between drift rate and measured mood 
# using pymc3 package
# test arousal
with pm.Model() as model_arousal:
    alpha = pm.Normal("alpha", mu=0, sigma=10)
    beta_arousal = pm.Normal("beta", mu=0, sigma=10)
    sigma = pm.HalfNormal("sigma", sigma=1)
    
    mu_arousal = alpha + beta_arousal * mood["movie_arousal_measured"].values
    
    v_0V1A = pm.Normal('v_0V1A', mu=mu_arousal, sd=sigma, observed=mood['v_0V1A'].values)
    
    trace_a = pm.sample(10000, tune=2000)
# save the traces of the correlation to csv for inference testing and plotting
pd.DataFrame({"cor_valence": trace_v['beta'],
              "cor_arousal": trace_a['beta']}).to_csv("plot/study1_correlation_traces.csv")<jupyter_output><empty_output><jupyter_text>Exploratory Full Regression Analysis<jupyter_code>m_full_model = hddm.HDDMRegressor(data, ["v ~ movie_valence + movie_arousal + movie_valence*movie_arousal",
                                   "a ~ movie_valence + movie_arousal + movie_valence*movie_arousal",
                                   "t ~ movie_valence + movie_arousal + movie_valence*movie_arousal",
                                   "z ~ movie_valence + movie_arousal + movie_valence*movie_arousal"],
                            include=('v', 'a', 't', 'z'),
                            p_outlier=0.05)
m_full_model.find_starting_values()
m_full_model.sample(10000, burn=2000, dbname='traces_regression.db', db='pickle')
m_full_model.save('traces_regression')<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup