Pilot Subjects Analysis

Dataset gathered from pilot subjects. All subjects were healthy adults. TODO - ADD DESCRIPTION FOR PILOT SUBJECTS


In [1]:
"""
Environment setup
"""
%matplotlib inline
%cd /lang_dec
import warnings; warnings.filterwarnings('ignore')
import hddm
import numpy as np
import matplotlib.pyplot as plt

from utils import model_tools


/Users/celefthe/Programming/projects/language_decision

In [2]:
# Import pilot subject data (as pandas dataframe)
pilot_data = hddm.load_csv('/lang_dec/data/pilot_clean.csv')

Reaction Time & Accuracy

Here we include the reaction time and accuracy metrics from the original dataset


In [19]:
us = pilot_data.loc[pilot_data['stim'] == 'US']
ss = pilot_data.loc[pilot_data['stim'] == 'SS']
cp = pilot_data.loc[pilot_data['stim'] == 'CP']
cs = pilot_data.loc[pilot_data['stim'] == 'CS']

plt.boxplot([ss.rt.values, cp.rt.values, cs.rt.values, us.rt.values], 
            labels=('SS', 'CP', 'CS', 'US'),)
plt.title('Comparison of Reaction Time Differences Between Stimuli Groups')
plt.show()



In [83]:
ss_accuracy = (len([x for x in ss.response.values if x >= 1]) / len(ss.response.values)) * 100 
cp_accuracy = (len([x for x in cp.response.values if x >= 1]) / len(cp.response.values)) * 100
cs_accuracy = (len([x for x in cs.response.values if x >= 1]) / len(cs.response.values)) * 100
us_accuracy = (len([x for x in us.response.values if x >= 1]) / len(us.response.values)) * 100

print("SS Accuracy: " + str(ss_accuracy) + "%")
print("CP Accuracy: " + str(cp_accuracy) + "%")
print("CS Accuracy: " + str(cs_accuracy) + "%")
print("US Accuracy: " + str(us_accuracy) + "%")


SS Accuracy: 98.12108559498957%
CP Accuracy: 92.08333333333333%
CS Accuracy: 98.52941176470588%
US Accuracy: 99.37369519832986%

In [117]:
plt.bar([1,2,3,4], 
        [ss_accuracy, cp_accuracy, cs_accuracy, us_accuracy])


Out[117]:
<Container object of 4 artists>

Does the drift rate depend on stimulus type?


In [23]:
"""
Construct a model where the drift rate depends on the stimulus type
"""
pilot_model = hddm.HDDM(pilot_data, depends_on={'v': 'stim'}, bias=True)
pilot_model.find_starting_values()
pilot_model.sample(9000, burn=200, dbname='language_decision/models/pilot', db='txt')


 [-----------------100%-----------------] 9000 of 9000 complete in 500.8 sec
Out[23]:
<pymc.MCMC.MCMC at 0x10fd1e240>

Convergence Checks

Before carrying on with analysing the output of the model, we need to check that the markov chains have properly converged. There's a number of ways to do this, which the authors of the hddm library recommend$^1$. We'll begin by visually inspecting the MCMC posterior plots.


In [136]:
pilot_model.plot_posteriors()


Plotting a
Plotting a_std
Plotting v(CP)
Plotting v(CS)
Plotting v(SS)
Plotting v(US)
Plotting v_std
Plotting t
Plotting t_std
Plotting z
Plotting z_std

PASS - No problematic patterns, such as drifts or large jumps, can be in any of the traces above. Autocorrelation also drops to zero quite quickly when considering past samples - which is what we want.

We can also formally test for model convergence using the Gelman-Rubin R statistic$^2$, which compares the within- and between-chain variance of different runs of the same model; models converge if variables are between $0.98$ and $1.02$. A simple algorithm to check this is outlined below:


In [153]:
models = []
for i in range(100):
    m = hddm.HDDM(pilot_data, depends_on={'v': 'stim'}, bias=True)
    m.find_starting_values()
    m.sample(6000, burn=20)
    models.append(m)

model_tools.check_convergence(models)


 [-----------------100%-----------------] 6000 of 6000 complete in 341.5 sec
No convergence problems detected!
Out[153]:
True

PASS - Formal testing reveals no convergence problems; Gelman-Rubin R statistic values for all model variablesfall within the desired range ($0.98$ to $1.02$)

Drift Rate Analysis

Here, we examine whether the type of stimulus significantly affects the drift rate of the decision-making process.


In [137]:
pilot_stats = pilot_model.gen_stats()
print("Threshold (a) Mean: " + str(pilot_stats['mean']['a']) + " (std: " + str(pilot_stats['std']['a']) + ")")
print("Non-Decision (t) Mean: " + str(pilot_stats['mean']['t']) + " (std: " + str(pilot_stats['std']['t']) + ")")
print("Bias (z) Mean: " + str(pilot_stats['mean']['z']) + " (std: " + str(pilot_stats['std']['z']) + ")")
print("SS Mean Drift Rate: " + str(pilot_stats['mean']['v(SS)']) + " (std: " + str(pilot_stats['std']['v(SS)']) + ")")
print("CP Mean Drift Rate: " + str(pilot_stats['mean']['v(CP)']) + " (std: " + str(pilot_stats['std']['v(CP)']) + ")")
print("CS Mean Drift Rate: " + str(pilot_stats['mean']['v(CS)']) + " (std: " + str(pilot_stats['std']['v(CS)']) + ")")
print("US Mean Drift Rate: " + str(pilot_stats['mean']['v(US)']) + " (std: " + str(pilot_stats['std']['v(US)']) + ")")


Threshold (a) Mean: 2.33606146724 (std: 0.133133292101)
Non-Decision (t) Mean: 0.55566081719 (std: 0.0741704228536)
Bias (z) Mean: 0.398816016617 (std: 0.0167336998574)
SS Mean Drift Rate: 2.66498783699 (std: 0.119452551423)
CP Mean Drift Rate: 1.72591437758 (std: 0.109362267982)
CS Mean Drift Rate: 2.21468336544 (std: 0.114123188241)
US Mean Drift Rate: 2.56013623275 (std: 0.116672600386)

In [138]:
v_SS, v_CP, v_CS, v_US = pilot_model.nodes_db.node[['v(SS)', 'v(CP)', 'v(CS)', 'v(US)']]

hddm.analyze.plot_posterior_nodes([v_SS, v_CP, v_CS, v_US])


Estimating the model in a Bayesian framework allows us to test for significance directly on the posterior, rather than having to rely on conventional frequency statistical tests:


In [139]:
print('P(SS > US) = ' + str((v_SS.trace() > v_US.trace()).mean()))
print('P(CP > SS) = ' + str((v_CP.trace() > v_SS.trace()).mean()))
print('P(CS > SS) = ' + str((v_CS.trace() > v_SS.trace()).mean()))
print('P(CP > CS) = ' + str((v_CP.trace() > v_CS.trace()).mean()))


P(SS > US) = 0.757461024499
P(CP > SS) = 0.0
P(CS > SS) = 0.000890868596882
P(CP > CS) = 0.000668151447661
  • The drift rate for CP is significantly lower than all other conditions
  • The drift rate for CS is significantly lower than SS and US, but significantly higher than CP
  • The drift rates for SS and US are not significantly different

In [140]:
"""
Distribution for the non-decision time t
"""
time_nondec = pilot_model.nodes_db.node[['t']]

hddm.analyze.plot_posterior_nodes(time_nondec)



In [141]:
"""
Distribution of bias z
"""
z = pilot_model.nodes_db.node[['z']]

hddm.analyze.plot_posterior_nodes(z)


Does the stimulus type affect the distance between the two boundaries (threshold)?

Threshold (or a) describes the relative difference in the distance between the upper and lower response boundaries of the DDM.

We explore whether stimulus type affects the threshold / distance between the two boundaries


In [3]:
pilot_model_threshold = hddm.HDDM(pilot_data, depends_on={'v': 'stim', 'a': 'stim'}, bias=True)
pilot_model_threshold.find_starting_values()
pilot_model_threshold.sample(10000, burn=200, dbname='language_decision/models/pilot_threshold', db='txt')


 [-----------------100%-----------------] 10001 of 10000 complete in 633.5 sec
Out[3]:
<pymc.MCMC.MCMC at 0x111cf1860>

Convergence checks


In [61]:
models_threshold = []
for i in range(5):
    m = hddm.HDDM(pilot_subjects, depends_on={'v': 'stim', 'a': 'stim'})
    m.find_starting_values()
    m.sample(6000, burn=20)
    models_threshold.append(m)

model_tools.check_convergence(models_threshold)


 [-----------------100%-----------------] 6000 of 6000 complete in 294.9 sec
No convergence problems detected!
Out[61]:
True

Threshold analysis

Since models converge, we can check the posteriors for significant differences in threshold between stimuli groups as we did for drift rates.


In [150]:
a_SS, a_CP, a_CS, a_US = pilot_model_threshold.nodes_db.node[['a(SS)', 'a(CP)', 'a(CS)', 'a(US)']]

hddm.analyze.plot_posterior_nodes([a_SS, a_CP, a_CS, a_US])



In [151]:
print('P(SS > US) = ' + str((a_SS.trace() > a_US.trace()).mean()))
print('P(CP > SS) = ' + str((a_CP.trace() > a_SS.trace()).mean()))
print('P(CS > SS) = ' + str((a_CS.trace() > a_SS.trace()).mean()))
print('P(CP > CS) = ' + str((a_CP.trace() > a_CS.trace()).mean()))
print('P(CS > US) = ' + str((a_CS.trace() > a_US.trace()).mean()))


P(SS > US) = 0.0581291759465
P(CP > SS) = 0.524276169265
P(CS > SS) = 0.880846325167
P(CP > CS) = 0.124053452116
P(CS > US) = 0.344432071269

No significant differences in threshold found across the different groups.


In [152]:
print("a(US) mean: " + str(a_US.trace().mean()))
print("a(SS) mean: " + str(a_SS.trace().mean()))
print("a(CS) mean: " + str(a_CS.trace().mean()))
print("a(CP) mean: " + str(a_CP.trace().mean()))


a(US) mean: 2.61989955936
a(SS) mean: 2.30079956919
a(CS) mean: 2.54098385298
a(CP) mean: 2.31454205216

Lumped Model


In [4]:
pilot_model_lumped = hddm.HDDM(pilot_data)
pilot_model_lumped.find_starting_values()
pilot_model_lumped.sample(10000, burn=200, dbname='language_decision/models/pilot_lumped', db='txt')


 [-----------------100%-----------------] 10000 of 10000 complete in 284.6 sec
Out[4]:
<pymc.MCMC.MCMC at 0x111d92978>

In [5]:
pilot_model_lumped.plot_posteriors()


Plotting a
Plotting a_std
Plotting v
Plotting v_std
Plotting t
Plotting t_std