notebook.community

Edit and run

A Diffusion Decision Model for Language-based Decisions

Here, we use the hddm (Python) library$^1$ to model language-based decision making based on response-time & accuracy data from a binary forced-choice decision task. During each trial, subjects were presented with a digital image of an item (visual stimulus) and a verbal description (auditory stimulus). At the end of each trial, subjects were asked to indicate whether the two types of stimuli agreed (i.e. visual and auditory stimulus pointed to the same object) or were different (i.e. non-matching sound and image). Auditory and visual stimuli pairs fell under four distinct categories, as summarised below:

  • Condition 1 (SS): Same stimuli (see cat / hear cat)
  • Condition 2 (CP): Phonologically congruent (see cat / hear hat)
  • Condition 3 (CS): Semantically congruent (see cat / hear dog)
  • Condition 4 (US): Unrelated stimuli (see house / hear cat)




In [ ]:



    


# Environment setup
%matplotlib inline
%cd /lang_dec
# Imports
import warnings; warnings.filterwarnings('ignore')
import hddm
import scipy
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import bayesian_bootstrap.bootstrap as bootstrap
from utils import model_tools, signal_detection



    


















    






/Users/celefthe/Programming/projects/language_decision

















In [2]:



    


# Import pilot models
pilot_data = hddm.load_csv('/lang_dec/data/pilot_clean.csv')
pilot_model = hddm.HDDM(pilot_data, depends_on={'v': 'stim'}, bias=True)
pilot_model.load_db(dbname='language_decision/models/pilot', db='txt')

pilot_model_lumped = hddm.HDDM(pilot_data)
pilot_model_lumped.load_db(dbname='language_decision/models/pilot_lumped', db='txt')

#pilot_model_threshold = hddm.HDDM(pilot_data, depends_on={'v': 'stim', 'a': 'stim'})
#pilot_model_threshold.load_db(dbname='language_decision/models/pilot_threshold', db='txt')



    


















    
Out[2]:








<hddm.models.hddm_info.HDDM at 0x116568cf8>

















In [3]:



    


# Import control models
controls_data = hddm.load_csv('/lang_dec/data/controls_clean.csv')
controls_model = hddm.HDDM(controls_data, depends_on={'v': 'stim'}, bias=True)
controls_model.load_db(dbname='language_decision/models/controls', db='txt')

controls_model_lumped = hddm.HDDM(controls_data)
controls_model_lumped.load_db(dbname='language_decision/models/controls_lumped', db='txt')

#controls_model_threshold = hddm.HDDM(controls_data, depends_on={'v': 'stim', 'a': 'stim'})
#controls_model_threshold.load_db(dbname='language_decision/models/controls_threshold', db='txt')



    


















    
Out[3]:








<hddm.models.hddm_info.HDDM at 0x116f89cf8>

















In [4]:



    


# Import patient models
patients_data = hddm.load_csv('/lang_dec/data/patients_clean.csv')
patients_model = hddm.HDDM(patients_data, depends_on={'v': 'stim'}, bias=True)
patients_model.load_db(dbname='language_decision/models/patients', db='txt')

patients_model_lumped = hddm.HDDM(patients_data)
patients_model_lumped.load_db(dbname='language_decision/models/patients_lumped', db='txt')

#patients_model_threshold = hddm.HDDM(patients_data, depends_on={'v': 'stim', 'a': 'stim'})
#patients_model_threshold.load_db(dbname='language_decision/models/patients_threshold', db='txt')



    


















    
Out[4]:








<hddm.models.hddm_info.HDDM at 0x1170ba518>

















In [5]:



    


# Bootstrap function
def bayes_bootstrap(observed, stat_function, replications, resample_size):
    data_bs = bootstrap.bayesian_bootstrap(
        list(filter(None, observed)), stat_function, replications, resample_size)
    return {
        'bootstrap_data': data_bs,
        'central_tendency': {
            str(stat_function.__name__): stat_function(data_bs)
        },
        'CIs': {
            'high':  scipy.stats.mstats.mquantiles(data_bs, prob=[0.975])[0],
            'low': scipy.stats.mstats.mquantiles(data_bs, prob=[0.025])[0]
        } if stat_function.__name__ == 'median' else {
            'high':  stat_function(data_bs) + (1.96 * scipy.std(data_bs)),
            'low': stat_function(data_bs) - (1.96 * scipy.std(data_bs))
        }
    }



    











In [6]:



    


# Bootstrap comparison function
def compare_bootstrap_medians(sample1, sample2):
    median_diffs = [x - y for x in sample1 for y in sample2]
    return {
        'median_diffs': median_diffs,
        'central_tendency': scipy.median(median_diffs),
        'CIs': {
            'high':  scipy.stats.mstats.mquantiles(median_diffs, prob=[0.975])[0],
            'low': scipy.stats.mstats.mquantiles(median_diffs, prob=[0.025])[0]
        }
    }

Performance Comparison





In [7]:



    


# Gather data

pilot_us = pilot_data.loc[pilot_data['stim'] == 'US']
pilot_ss = pilot_data.loc[pilot_data['stim'] == 'SS']
pilot_cp = pilot_data.loc[pilot_data['stim'] == 'CP']
pilot_cs = pilot_data.loc[pilot_data['stim'] == 'CS']

controls_us = controls_data.loc[controls_data['stim'] == 'US']
controls_ss = controls_data.loc[controls_data['stim'] == 'SS']
controls_cp = controls_data.loc[controls_data['stim'] == 'CP']
controls_cs = controls_data.loc[controls_data['stim'] == 'CS']

patients_us = patients_data.loc[patients_data['stim'] == 'US']
patients_ss = patients_data.loc[patients_data['stim'] == 'SS']
patients_cp = patients_data.loc[patients_data['stim'] == 'CP']
patients_cs = patients_data.loc[patients_data['stim'] == 'CS']

Reaction Time





In [124]:



    


plt.figure(1)

plt.subplot(221)
vp_rt_ctrl = sns.violinplot(
    data=[controls_ss.rt.values, controls_cp.rt.values, 
          controls_cs.rt.values, controls_us.rt.values], 
    inner='box')
vp_rt_ctrl.set_xticklabels(['SS','CP','CS', 'US'])
vp_rt_ctrl.set_xlabel("Stimulus Category")
vp_rt_ctrl.set_ylabel("Reaction Time (s)")
vp_rt_ctrl.set_title("Elderly Controls")
# statistical annotation
x1, x2, x3, x4 = 0, 1, 2, 3   # columns (first column: 0, see plt.xticks())
y, h, col = 4.05, 0.05, 'k'
#Significant
plt.plot([x1, x1, x2, x2], [y, y+h, y+h, y], lw=1.5, c=col)
plt.plot([x1, x1, x3, x3], [y, y+h, y+h, y], lw=1.5, c=col)
plt.plot([x1, x1, x4, x4], [y, y+h, y+h, y], lw=1.5, c=col)
plt.text((x1+x2)*.5, y+h, "*", ha='center', va='bottom', color=col)
plt.ylim(0, 4.60)

plt.subplot(222)
vp_rt_ptns = sns.violinplot(
    data=[patients_ss.rt.values, patients_cp.rt.values, 
          patients_cs.rt.values, patients_us.rt.values], 
    inner='box')
vp_rt_ptns.set_xticklabels(['SS','CP','CS', 'US'])
vp_rt_ptns.set_xlabel("Stimulus Category")
vp_rt_ptns.set_ylabel("Reaction Time (s)")
vp_rt_ptns.set_title("Elderly Patients")
# statistical annotation
x1, x2, x3, x4 = 0, 1, 2, 3   # columns (first column: 0, see plt.xticks())
y, h, col = 4.25, 0.05, 'k'
#Significant
plt.plot([x1, x1, x2, x2], [y, y+h, y+h, y], lw=1.5, c=col)
plt.plot([x1, x1, x3, x3], [y, y+h, y+h, y], lw=1.5, c=col)
plt.plot([x1, x1, x4, x4], [y, y+h, y+h, y], lw=1.5, c=col)
plt.text((x1+x2)*.5, y+h, "*", ha='center', va='bottom', color=col)
plt.ylim(0, 4.60)

plt.subplot(223)
vp_rt_plt = sns.violinplot(
    data=[pilot_ss.rt.values, pilot_cp.rt.values, pilot_cs.rt.values, pilot_us.rt.values], 
    inner='box')
vp_rt_plt.set_xticklabels(['SS','CP','CS', 'US'])
vp_rt_plt.set_xlabel("Stimulus Category")
vp_rt_plt.set_ylabel("Reaction Time (s)")
vp_rt_plt.set_title("Young")
plt.ylim(0, 4.60)
# statistical annotation
x1, x2, x3, x4 = 0, 1, 2, 3   # columns (first column: 0, see plt.xticks())
y, h, col = 4.05, 0.05, 'k'
#Significant
plt.plot([x1, x1, x2, x2], [y, y+h, y+h, y], lw=1.5, c=col)
plt.plot([x1, x1, x3, x3], [y, y+h, y+h, y], lw=1.5, c=col)
plt.plot([x1, x1, x4, x4], [y, y+h, y+h, y], lw=1.5, c=col)
plt.text((x1+x2)*.5, y+h, "*", ha='center', va='bottom', color=col)

plt.subplots_adjust(top=2.5, bottom=0.001, left=0.10, right=2.5, hspace=0.30,
                    wspace=0.15)


print("N_pilot = " + str(len(set(pilot_data.subj_idx))))
print("N_controls = " + str(len(set(controls_data.subj_idx))))

print("Pilot medians:")
print("\t SS:" + str(scipy.median(pilot_ss.rt.values)))
print("\t CP:" + str(scipy.median(pilot_cp.rt.values)))
print("\t CS:" + str(scipy.median(pilot_cs.rt.values)))
print("\t US:" + str(scipy.median(pilot_us.rt.values)))

print("Controls medians:")
print("\t SS:" + str(scipy.median(controls_ss.rt.values)))
print("\t CP:" + str(scipy.median(controls_cp.rt.values)))
print("\t CS:" + str(scipy.median(controls_cs.rt.values)))
print("\t US:" + str(scipy.median(controls_us.rt.values)))

print("Patient medians:")
print("\t SS:" + str(scipy.median(patients_ss.rt.values)))
print("\t CP:" + str(scipy.median(patients_cp.rt.values)))
print("\t CS:" + str(scipy.median(patients_cs.rt.values)))
print("\t US:" + str(scipy.median(patients_us.rt.values)))



    


















    






N_pilot = 16
N_controls = 6
Pilot medians:
	 SS:1.0177292088
	 CP:1.20372684853
	 CS:1.10076620218
	 US:1.04213379743
Controls medians:
	 SS:1.08752833395
	 CP:1.20626424568
	 CS:1.19874492421
	 US:1.13699611987
Patient medians:
	 SS:1.43588068862
	 CP:1.69737381073
	 CS:1.80817536085
	 US:1.68778825953










    
























In [10]:



    


rt_comps = dict()



    











In [129]:



    


# Collate analysis data (actual analysis follows below)
for subj in rt_comps.keys():
    print(subj)
    for stim in rt_comps[subj].keys():
        print(stim)
        print({x: rt_comps[subj][stim][x] for x in rt_comps[subj][stim].keys() if x != 'bootstrap_data'})
    print('--------------------')



    


















    






controls
ss_cs
{'CIs': {'low': 0.099510060484133148, 'high': 0.11712019733414313}, 'central_tendency': {'median': 0.10775729840315762}}
ss_cp
{'CIs': {'low': 0.11252954218161926, 'high': 0.13030634585760703}, 'central_tendency': {'median': 0.1219901352542363}}
ss_us
{'CIs': {'low': 0.049197064934541984, 'high': 0.063114956223817117}, 'central_tendency': {'median': 0.05588531372086436}}
--------------------
pilot_controls
us
{'CIs': {'low': 0.067217717473522537, 'high': 0.075719666857572615}, 'central_tendency': {'median': 0.071644066755652602}}
cp
{'CIs': {'low': 0.067213465211475656, 'high': 0.075732170024561896}, 'central_tendency': {'median': 0.071594170937942236}}
ss
{'CIs': {'low': 0.067263889550740613, 'high': 0.075872322605174528}, 'central_tendency': {'median': 0.07181934910340726}}
cs
{'CIs': {'low': 0.067290146482134855, 'high': 0.075591179492014368}, 'central_tendency': {'median': 0.071743240810974274}}
--------------------
pilot
ss_cs
{'CIs': {'low': 0.073783205822983339, 'high': 0.079659070755213279}, 'central_tendency': {'median': 0.076602441513273362}}
ss_cp
{'CIs': {'low': 0.16038673141884374, 'high': 0.16698965363138646}, 'central_tendency': {'median': 0.16381164673077817}}
ss_us
{'CIs': {'low': 0.022407778585182996, 'high': 0.027373195997997687}, 'central_tendency': {'median': 0.024965266840354272}}
--------------------
patients
ss_cs
{'CIs': {'low': 0.29599476976612082, 'high': 0.31917081819254445}, 'central_tendency': {'median': 0.30766554123056267}}
ss_cp
{'CIs': {'low': 0.2350414122118572, 'high': 0.25700254479819007}, 'central_tendency': {'median': 0.24646836028887248}}
ss_us
{'CIs': {'low': 0.19675814170354564, 'high': 0.21911037398669286}, 'central_tendency': {'median': 0.2080033543101083}}
--------------------
controls_patients
us
{'CIs': {'low': 0.49753372910130694, 'high': 0.52314437979941997}, 'central_tendency': {'median': 0.51096195127047395}}
cp
{'CIs': {'low': 0.42348480835153168, 'high': 0.44968687309863531}, 'central_tendency': {'median': 0.43764282144593381}}
ss
{'CIs': {'low': 0.32040392983589755, 'high': 0.33719697209113753}, 'central_tendency': {'median': 0.32881700199277469}}
cs
{'CIs': {'low': 0.52794412770992949, 'high': 0.55523939797497168}, 'central_tendency': {'median': 0.54114360916037196}}
--------------------

Bootstrap Median Differences Between Subject Groups

Pilot vs Controls




In [12]:



    


subj_group = 'pilot_controls'
rt_comps[subj_group] = dict()



    











In [13]:



    


subj_group = 'pilot_controls'
stim_group = 'ss'
diffs = [y - x for x in pilot_ss.rt.values for y in controls_ss.rt.values]

rt_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))


#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(rt_comps[subj_group][stim_group]['CIs'])
print(rt_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.067263889550740613, 'high': 0.075872322605174528}
{'median': 0.07181934910340726}

















In [14]:



    


subj_group = 'pilot_controls'
stim_group = 'cs'
diffs = [y - x for x in pilot_ss.rt.values for y in controls_ss.rt.values]

rt_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))


#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(rt_comps[subj_group][stim_group]['CIs'])
print(rt_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.067290146482134855, 'high': 0.075591179492014368}
{'median': 0.071743240810974274}

















In [15]:



    


subj_group = 'pilot_controls'
stim_group = 'cp'
diffs = [y - x for x in pilot_ss.rt.values for y in controls_ss.rt.values]

rt_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))


#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(rt_comps[subj_group][stim_group]['CIs'])
print(rt_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.067213465211475656, 'high': 0.075732170024561896}
{'median': 0.071594170937942236}

















In [16]:



    


subj_group = 'pilot_controls'
stim_group = 'us'
diffs = [y - x for x in pilot_ss.rt.values for y in controls_ss.rt.values]

rt_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))


#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(rt_comps[subj_group][stim_group]['CIs'])
print(rt_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.067217717473522537, 'high': 0.075719666857572615}
{'median': 0.071644066755652602}
Controls vs Patients




In [17]:



    


subj_group = 'controls_patients'
rt_comps[subj_group] = dict()



    











In [18]:



    


subj_group = 'controls_patients'
stim_group = 'ss'
diffs = [y - x for x in controls_ss.rt.values for y in patients_ss.rt.values]

rt_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(rt_comps[subj_group][stim_group]['CIs'])
print(rt_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.32040392983589755, 'high': 0.33719697209113753}
{'median': 0.32881700199277469}

















In [19]:



    


subj_group = 'controls_patients'
stim_group = 'cs'
diffs = [y - x for x in controls_cs.rt.values for y in patients_cs.rt.values]

rt_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(rt_comps[subj_group][stim_group]['CIs'])
print(rt_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.52794412770992949, 'high': 0.55523939797497168}
{'median': 0.54114360916037196}

















In [20]:



    


subj_group = 'controls_patients'
stim_group = 'cp'
diffs = [y - x for x in controls_cp.rt.values for y in patients_cp.rt.values]

rt_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(rt_comps[subj_group][stim_group]['CIs'])
print(rt_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.42348480835153168, 'high': 0.44968687309863531}
{'median': 0.43764282144593381}

















In [21]:



    


subj_group = 'controls_patients'
stim_group = 'us'
diffs = [y - x for x in controls_us.rt.values for y in patients_us.rt.values]

rt_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(rt_comps[subj_group][stim_group]['CIs'])
print(rt_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.49753372910130694, 'high': 0.52314437979941997}
{'median': 0.51096195127047395}

Bootstrap Median Differences Between Stimuli Groups (Incogruency Effect Only)

Pilot




In [22]:



    


subj_group = 'pilot'
rt_comps[subj_group] = dict()



    











In [23]:



    


subj_group = 'pilot'
stim_group = 'ss_us'
diffs = [y - x for x in pilot_ss.rt.values for y in pilot_us.rt.values]

rt_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(rt_comps[subj_group][stim_group]['CIs'])
print(rt_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.022407778585182996, 'high': 0.027373195997997687}
{'median': 0.024965266840354272}

















In [24]:



    


subj_group = 'pilot'
stim_group = 'ss_cs'
diffs = [y - x for x in pilot_ss.rt.values for y in pilot_cs.rt.values]

rt_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(rt_comps[subj_group][stim_group]['CIs'])
print(rt_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.073783205822983339, 'high': 0.079659070755213279}
{'median': 0.076602441513273362}

















In [25]:



    


subj_group = 'pilot'
stim_group = 'ss_cp'
diffs = [y - x for x in pilot_ss.rt.values for y in pilot_cp.rt.values]

rt_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(rt_comps[subj_group][stim_group]['CIs'])
print(rt_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.16038673141884374, 'high': 0.16698965363138646}
{'median': 0.16381164673077817}
Controls




In [26]:



    


subj_group = 'controls'
rt_comps[subj_group] = dict()



    











In [27]:



    


subj_group = 'controls'
stim_group = 'ss_us'
diffs = [y - x for x in controls_ss.rt.values for y in controls_us.rt.values]

rt_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(rt_comps[subj_group][stim_group]['CIs'])
print(rt_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.049197064934541984, 'high': 0.063114956223817117}
{'median': 0.05588531372086436}

















In [28]:



    


subj_group = 'controls'
stim_group = 'ss_cs'
diffs = [y - x for x in controls_ss.rt.values for y in controls_cs.rt.values]

rt_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(rt_comps[subj_group][stim_group]['CIs'])
print(rt_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.099510060484133148, 'high': 0.11712019733414313}
{'median': 0.10775729840315762}

















In [29]:



    


subj_group = 'controls'
stim_group = 'ss_cp'
diffs = [y - x for x in controls_ss.rt.values for y in controls_cp.rt.values]

rt_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(rt_comps[subj_group][stim_group]['CIs'])
print(rt_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.11252954218161926, 'high': 0.13030634585760703}
{'median': 0.1219901352542363}
Patients




In [30]:



    


subj_group = 'patients'
rt_comps[subj_group] = dict()



    











In [31]:



    


subj_group = 'patients'
stim_group = 'ss_us'
diffs = [y - x for x in patients_ss.rt.values for y in patients_us.rt.values]

rt_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(rt_comps[subj_group][stim_group]['CIs'])
print(rt_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.19675814170354564, 'high': 0.21911037398669286}
{'median': 0.2080033543101083}

















In [32]:



    


subj_group = 'patients'
stim_group = 'ss_cs'
diffs = [y - x for x in patients_ss.rt.values for y in patients_cs.rt.values]

rt_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(rt_comps[subj_group][stim_group]['CIs'])
print(rt_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.29599476976612082, 'high': 0.31917081819254445}
{'median': 0.30766554123056267}

















In [33]:



    


subj_group = 'patients'
stim_group = 'ss_cp'
diffs = [y - x for x in patients_ss.rt.values for y in patients_cp.rt.values]

rt_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(rt_comps[subj_group][stim_group]['CIs'])
print(rt_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.2350414122118572, 'high': 0.25700254479819007}
{'median': 0.24646836028887248}

Accuracy





In [34]:



    


def get_accuracies(subjects):
    accuracies = []
    grouped_subjects = subjects.groupby('subj_idx')
    for subject in [grouped_subjects.get_group(x) for x in grouped_subjects.groups]:
        accuracies.append(len([x for x in subject.response.values if x >= 1]) / len(subject.response.values))
    return accuracies



    











In [35]:



    


pilot_ss_accuracies = get_accuracies(pilot_ss)
pilot_cp_accuracies = get_accuracies(pilot_cp)
pilot_cs_accuracies = get_accuracies(pilot_cs)
pilot_us_accuracies = get_accuracies(pilot_us)

controls_ss_accuracies = get_accuracies(controls_ss)
controls_cp_accuracies = get_accuracies(controls_cp)
controls_cs_accuracies = get_accuracies(controls_cs)
controls_us_accuracies = get_accuracies(controls_us)

patients_ss_accuracies = get_accuracies(patients_ss) 
patients_cp_accuracies = get_accuracies(patients_cp)
patients_cs_accuracies = get_accuracies(patients_cs)
patients_us_accuracies = get_accuracies(patients_us)



    











In [128]:



    


plt.figure(1)

plt.subplot(221)
box_acc_ctrl = sns.boxplot(
    data=[controls_ss_accuracies, controls_cp_accuracies, 
          controls_cs_accuracies, controls_us_accuracies])
box_acc_ctrl.set_xticklabels(['SS','CP','CS', 'US'])
box_acc_ctrl.set_xlabel("Stimulus Category")
box_acc_ctrl.set_ylabel("Accuracy")
box_acc_ctrl.set_title("Elderly Controls")
plt.ylim(0, 1.15)
# statistical annotation
x1, x2, x3, x4 = 0, 1, 2, 3   # columns (first column: 0, see plt.xticks())
y, h, col = 1.03, 0.02, 'k'
#Significant
plt.plot([x1, x1, x2, x2], [y, y+h, y+h, y], lw=1.5, c=col)
plt.plot([x1, x1, x4, x4], [y, y+h, y+h, y], lw=1.5, c=col)
plt.text((x1+x2)*.5, y+h, "*", ha='center', va='bottom', color=col)
# NS
y, h, col = 0.5, 0.02, 'k'
plt.plot([x1, x1, x3, x3], [y, y-h, y-h, y], lw=1.5, c=col)
plt.text((x1+x3)*.5, y-h, "ns", ha='center', va='top', color=col)


plt.subplot(222)
box_acc_ptns = sns.boxplot(
    data=[patients_ss_accuracies, patients_cp_accuracies, 
          patients_cs_accuracies, patients_us_accuracies])
box_acc_ptns.set_xticklabels(['SS','CP','CS', 'US'])
box_acc_ptns.set_xlabel("Stimulus Category")
box_acc_ptns.set_ylabel("Accuracy")
box_acc_ptns.set_title("Elderly Patients")
plt.ylim(0, 1.15)
# statistical annotation
x1, x2, x3, x4 = 0, 1, 2, 3   # columns (first column: 0, see plt.xticks())
y, h, col = 1.03, 0.02, 'k'
#Significant
plt.plot([x1, x1, x2, x2], [y, y+h, y+h, y], lw=1.5, c=col)
plt.plot([x1, x1, x3, x3], [y, y+h, y+h, y], lw=1.5, c=col)
plt.text((x1+x2)*.5, y+h, "*", ha='center', va='bottom', color=col)
# NS
y, h, col = 1.05, 0.05, 'k'
plt.plot([x1, x1, x4, x4], [y, y+h, y+h, y], lw=1.5, c=col)
plt.text((x1+x4)*.5, y+h, "ns", ha='center', va='bottom', color=col)

plt.subplot(223)
box_acc_plt = sns.boxplot(
    data=[pilot_ss_accuracies, pilot_cp_accuracies, pilot_cs_accuracies, pilot_us_accuracies])
box_acc_plt.set_xticklabels(['SS','CP','CS', 'US'])
box_acc_plt.set_xlabel("Stimulus Category")
box_acc_plt.set_ylabel("Accuracy")
box_acc_plt.set_title("Young")
plt.ylim(0, 1.15)
# statistical annotation
x1, x2, x3, x4 = 0, 1, 2, 3   # columns (first column: 0, see plt.xticks())
y, h, col = 1.03, 0.02, 'k'
#Significant
plt.plot([x1, x1, x2, x2], [y, y+h, y+h, y], lw=1.5, c=col)
plt.plot([x1, x1, x4, x4], [y, y+h, y+h, y], lw=1.5, c=col)
plt.text((x1+x2)*.5, y+h, "*", ha='center', va='bottom', color=col)
# NS
y, h, col = 0.8, 0.02, 'k'
plt.plot([x1, x1, x3, x3], [y, y-h, y-h, y], lw=1.5, c=col)
plt.text((x1+x3)*.5, y-h, "ns", ha='center', va='top', color=col)

plt.subplots_adjust(top=2.5, bottom=0.001, left=0.10, right=2.5, hspace=0.30,
                    wspace=0.15)


plt.show()
print("N_pilot = " + str(len(set(pilot_data.subj_idx))))
print("N_controls = " + str(len(set(controls_data.subj_idx))))
print("N_patients = " + str(len(set(patients_data.subj_idx))))

print("Pilot medians:")
print("\t SS:" + str(scipy.median(pilot_ss_accuracies)))
print("\t CP:" + str(scipy.median(pilot_cp_accuracies)))
print("\t CS:" + str(scipy.median(pilot_cs_accuracies)))
print("\t US:" + str(scipy.median(pilot_us_accuracies)))

print("Controls medians:")
print("\t SS:" + str(scipy.median(controls_ss_accuracies)))
print("\t CP:" + str(scipy.median(controls_cp_accuracies)))
print("\t CS:" + str(scipy.median(controls_cs_accuracies)))
print("\t US:" + str(scipy.median(controls_us_accuracies)))

print("Patient medians:")
print("\t SS:" + str(scipy.median(patients_ss_accuracies)))
print("\t CP:" + str(scipy.median(patients_cp_accuracies)))
print("\t CS:" + str(scipy.median(patients_cs_accuracies)))
print("\t US:" + str(scipy.median(patients_us_accuracies)))



    


















    

















    






N_pilot = 16
N_controls = 6
N_patients = 10
Pilot medians:
	 SS:0.983333333333
	 CP:0.933333333333
	 CS:1.0
	 US:1.0
Controls medians:
	 SS:0.95
	 CP:0.830459770115
	 CS:0.933333333333
	 US:1.0
Patient medians:
	 SS:0.923850574713
	 CP:0.625
	 CS:0.644252873563
	 US:0.888095238095

















In [38]:



    


acc_comps = dict()



    











In [127]:



    


# Collate analysis data (actual analysis follows below)
for subj in acc_comps.keys():
    print(subj)
    for stim in acc_comps[subj].keys():
        print(stim)
        print({x: acc_comps[subj][stim][x] for x in acc_comps[subj][stim].keys() if x != 'bootstrap_data'})
    print('--------------------')



    


















    






controls
ss_cs
{'CIs': {'low': -0.083333333333333315, 'high': 0.033333333333333326}, 'central_tendency': {'median': 0.0011494252873562871}}
ss_cp
{'CIs': {'low': 0.085632183908046, 'high': 0.40000000000000002}, 'central_tendency': {'median': 0.1333333333333333}}
ss_us
{'CIs': {'low': -0.11666666666666664, 'high': -0.033333333333333326}, 'central_tendency': {'median': -0.033333333333333326}}
--------------------
pilot_controls
us
{'CIs': {'low': -0.033333333333333326, 'high': 0.19999999999999996}, 'central_tendency': {'median': 0.16666666666666663}}
cp
{'CIs': {'low': 0.099999999999999978, 'high': 0.2195402298850575}, 'central_tendency': {'median': 0.1333333333333333}}
ss
{'CIs': {'low': 0.033333333333333326, 'high': 0.1333333333333333}, 'central_tendency': {'median': 0.066666666666666652}}
cs
{'CIs': {'low': 0.034482758620689613, 'high': 0.066666666666666652}, 'central_tendency': {'median': 0.066666666666666652}}
--------------------
pilot
ss_cs
{'CIs': {'low': -0.033333333333333326, 'high': 0.0011494252873562871}, 'central_tendency': {'median': -0.033333333333333326}}
ss_cp
{'CIs': {'low': 0.066666666666666652, 'high': 0.099999999999999978}, 'central_tendency': {'median': 0.066666666666666652}}
ss_us
{'CIs': {'low': -0.033333333333333326, 'high': -0.033333333333333326}, 'central_tendency': {'median': -0.033333333333333326}}
--------------------
patients
ss_cs
{'CIs': {'low': 0.11343660968660971, 'high': 0.3666666666666667}, 'central_tendency': {'median': 0.2614942528735632}}
ss_cp
{'CIs': {'low': 0.21149425287356327, 'high': 0.34482758620689657}, 'central_tendency': {'median': 0.28045977011494261}}
ss_us
{'CIs': {'low': -0.066666666666666652, 'high': 0.21314285714285819}, 'central_tendency': {'median': -0.00512820512820511}}
--------------------
controls_patients
us
{'CIs': {'low': 0.19047619047619047, 'high': 0.83333333333333337}, 'central_tendency': {'median': 0.46666666666666667}}
cp
{'CIs': {'low': 0.064491379310344793, 'high': 0.31149425287356325}, 'central_tendency': {'median': 0.19999999999999996}}
ss
{'CIs': {'low': -0.033333333333333326, 'high': 0.07407407407407407}, 'central_tendency': {'median': 0.033333333333333326}}
cs
{'CIs': {'low': 0.15740740740740738, 'high': 0.65769230769230769}, 'central_tendency': {'median': 0.30000000000000004}}
--------------------

Bootstrap Median Differences Between Subject Groups

Pilot vs Controls




In [40]:



    


subj_group = 'pilot_controls'
acc_comps[subj_group] = dict()



    











In [41]:



    


subj_group = 'pilot_controls'
stim_group = 'ss'
diffs = [x - y for x in pilot_ss_accuracies for y in controls_ss_accuracies]

acc_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(acc_comps[subj_group][stim_group]['CIs'])
print(acc_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.033333333333333326, 'high': 0.1333333333333333}
{'median': 0.066666666666666652}

















In [42]:



    


subj_group = 'pilot_controls'
stim_group = 'cs'
diffs = [x - y for x in pilot_cs_accuracies for y in controls_cs_accuracies]

acc_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(acc_comps[subj_group][stim_group]['bootstrap_data'])
print(acc_comps[subj_group][stim_group]['CIs'])
print(acc_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.034482758620689613, 'high': 0.066666666666666652}
{'median': 0.066666666666666652}

















In [43]:



    


subj_group = 'pilot_controls'
stim_group = 'cp'
diffs = [x - y for x in pilot_cp_accuracies for y in controls_cp_accuracies]

acc_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(acc_comps[subj_group][stim_group]['bootstrap_data'])
print(acc_comps[subj_group][stim_group]['CIs'])
print(acc_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.099999999999999978, 'high': 0.2195402298850575}
{'median': 0.1333333333333333}

















In [44]:



    


subj_group = 'pilot_controls'
stim_group = 'us'
diffs = [x - y for x in pilot_us_accuracies for y in controls_us_accuracies]

acc_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(acc_comps[subj_group][stim_group]['CIs'])
print(acc_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': -0.033333333333333326, 'high': 0.19999999999999996}
{'median': 0.16666666666666663}
Controls vs Patients




In [45]:



    


subj_group = 'controls_patients'
acc_comps[subj_group] = dict()



    











In [46]:



    


subj_group = 'controls_patients'
stim_group = 'ss'
diffs = [y - x for x in patients_ss_accuracies for y in controls_ss_accuracies]

acc_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(acc_comps[subj_group][stim_group]['CIs'])
print(acc_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': -0.033333333333333326, 'high': 0.07407407407407407}
{'median': 0.033333333333333326}

















In [47]:



    


subj_group = 'controls_patients'
stim_group = 'cs'
diffs = [y - x for x in patients_cs_accuracies for y in controls_cs_accuracies]

acc_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(acc_comps[subj_group][stim_group]['CIs'])
print(acc_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.15740740740740738, 'high': 0.65769230769230769}
{'median': 0.30000000000000004}

















In [48]:



    


subj_group = 'controls_patients'
stim_group = 'cp'
diffs = [y - x for x in patients_cp_accuracies for y in controls_cp_accuracies]

acc_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(acc_comps[subj_group][stim_group]['CIs'])
print(acc_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.064491379310344793, 'high': 0.31149425287356325}
{'median': 0.19999999999999996}

















In [49]:



    


subj_group = 'controls_patients'
stim_group = 'us'
diffs = [y - x for x in patients_us_accuracies for y in controls_us_accuracies]

acc_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(acc_comps[subj_group][stim_group]['CIs'])
print(acc_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.19047619047619047, 'high': 0.83333333333333337}
{'median': 0.46666666666666667}

Bootstrap Median Differences Between Stimuli Categories (Within Subject Groups, Incogruency Only)

Pilot




In [50]:



    


subj_group = 'pilot'
acc_comps[subj_group] = dict()



    











In [51]:



    


subj_group = 'pilot'
stim_group = 'ss_us'
diffs = [x - y for x in pilot_ss_accuracies for y in pilot_us_accuracies]

acc_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(acc_comps[subj_group][stim_group]['CIs'])
print(acc_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': -0.033333333333333326, 'high': -0.033333333333333326}
{'median': -0.033333333333333326}

















In [52]:



    


subj_group = 'pilot'
stim_group = 'ss_cs'
diffs = [x - y for x in pilot_ss_accuracies for y in pilot_cs_accuracies]

acc_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(acc_comps[subj_group][stim_group]['CIs'])
print(acc_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': -0.033333333333333326, 'high': 0.0011494252873562871}
{'median': -0.033333333333333326}

















In [53]:



    


subj_group = 'pilot'
stim_group = 'ss_cp'
diffs = [x - y for x in pilot_ss_accuracies for y in pilot_cp_accuracies]

acc_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(acc_comps[subj_group][stim_group]['CIs'])
print(acc_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.066666666666666652, 'high': 0.099999999999999978}
{'median': 0.066666666666666652}
Controls




In [54]:



    


subj_group = 'controls'
acc_comps[subj_group] = dict()



    











In [55]:



    


subj_group = 'controls'
stim_group = 'ss_us'
diffs = [x - y for x in controls_ss_accuracies for y in controls_us_accuracies]

acc_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(acc_comps[subj_group][stim_group]['CIs'])
print(acc_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': -0.11666666666666664, 'high': -0.033333333333333326}
{'median': -0.033333333333333326}

















In [56]:



    


subj_group = 'controls'
stim_group = 'ss_cs'
diffs = [x - y for x in controls_ss_accuracies for y in controls_cs_accuracies]

acc_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(acc_comps[subj_group][stim_group]['CIs'])
print(acc_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': -0.083333333333333315, 'high': 0.033333333333333326}
{'median': 0.0011494252873562871}

















In [57]:



    


subj_group = 'controls'
stim_group = 'ss_cp'
diffs = [x - y for x in controls_ss_accuracies for y in controls_cp_accuracies]

acc_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(acc_comps[subj_group][stim_group]['CIs'])
print(acc_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.085632183908046, 'high': 0.40000000000000002}
{'median': 0.1333333333333333}
Patients




In [58]:



    


subj_group = 'patients'
acc_comps[subj_group] = dict()



    











In [59]:



    


subj_group = 'patients'
stim_group = 'ss_us'
diffs = [x - y for x in patients_ss_accuracies for y in patients_us_accuracies]

acc_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(acc_comps[subj_group][stim_group]['CIs'])
print(acc_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': -0.066666666666666652, 'high': 0.21314285714285819}
{'median': -0.00512820512820511}

















In [60]:



    


stim_group = 'ss_cs'
diffs = [x - y for x in patients_ss_accuracies for y in patients_cs_accuracies]

acc_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(acc_comps[subj_group][stim_group]['CIs'])
print(acc_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.11343660968660971, 'high': 0.3666666666666667}
{'median': 0.2614942528735632}

















In [61]:



    


stim_group = 'ss_cp'
diffs = [x - y for x in patients_ss_accuracies for y in patients_cp_accuracies]

acc_comps[subj_group][stim_group] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(rt_comps[subj_group][stim_group]['bootstrap_data'])
print(acc_comps[subj_group][stim_group]['CIs'])
print(acc_comps[subj_group][stim_group]['central_tendency'])



    


















    






{'low': 0.21149425287356327, 'high': 0.34482758620689657}
{'median': 0.28045977011494261}

Signal Detection Theory & d'

Hit: SS + "same" response
False Alarm: US/CS/CP + "same" response

Miss: SS + "diff" response
Correct Rejection: US/CS/CP + "diff" response

Stim Group 1: SS
Stim Group 2: US / CS / CP





In [62]:



    


def get_d_primes(dataset, stim1, stim2, include_id=False):
    d_primes = dict()
    subject_ids = set(dataset.subj_idx)
    for subject_id in subject_ids:
        stim1_data = dataset.loc[
            dataset['subj_idx'] == subject_id].loc[
            dataset['stim'] == str(stim1)]
        stim1_trials = len(stim1_data)
        hits = len(stim1_data.loc[
            dataset['response'] == 1.0])
        
        stim2_data = dataset.loc[
            dataset['subj_idx'] == subject_id].loc[
            dataset['stim'] == str(stim2)]
        stim2_trials = len(stim2_data)
        
        fas = len(stim2_data.loc[
            dataset['response'] == 0.0])
        
        if not stim1_trials or not stim2_trials:
            d_primes[subject_id] = None # N/A placeholder value
            continue
        
        d_prime = signal_detection.signal_detection(
            n_stim1=stim1_trials, 
            n_stim2=stim2_trials, 
            hits=hits, 
            false_alarms=fas)['d_prime']
        d_primes[subject_id] = d_prime
       
    if not include_id:
        return list(d_primes.values())
    return d_primes



    











In [63]:



    


## Pilot d'
pilot_dprimes = {
    'ss_us': get_d_primes(pilot_data, 'SS', 'US'),
    'ss_cs': get_d_primes(pilot_data, 'SS', 'CS'),
    'ss_cp': get_d_primes(pilot_data, 'SS', 'CP')
}

## Controls d'
controls_dprimes = {
    'ss_us': get_d_primes(controls_data, 'SS', 'US'),
    'ss_cs': get_d_primes(controls_data, 'SS', 'CS'),
    'ss_cp': get_d_primes(controls_data, 'SS', 'CP')
}

## Patients d'
patients_dprimes = {
    'ss_us': get_d_primes(patients_data, 'SS', 'US'),
    'ss_cs': get_d_primes(patients_data, 'SS', 'CS'),
    'ss_cp': get_d_primes(patients_data, 'SS', 'CP')
}



    











In [189]:



    


## Plotting d's
plt.figure(1)

plt.subplot(221)
box_dprime_ctrl = sns.boxplot(
    data=[controls_dprimes['ss_cp'],
          controls_dprimes['ss_cs'],
          controls_dprimes['ss_us']])
box_dprime_ctrl.set_xticklabels(['SS vs CP','SS vs CS', 'SS vs US'])
box_dprime_ctrl.set_xlabel("Stimulus Pair")
box_dprime_ctrl.set_ylabel("d prime")
box_dprime_ctrl.set_title("Elderly Controls")
box_dprime_ctrl.set_ylim(-2.5, 4.5)
# statistical annotation
x1, x2, x3 = 0, 1, 2   # columns (first column: 0, see plt.xticks())
y, h, col = 4.1, 0.05, 'k'
#Significant
plt.plot([x1, x1, x3, x3], [y, y+h, y+h, y], lw=1.5, c=col)
plt.text((x1+x3)*.5, y+h, "*", ha='center', va='bottom', color=col)


plt.subplot(222)
box_dprime_ptns = sns.boxplot(
    data=[patients_dprimes['ss_cp'],
          patients_dprimes['ss_cs'],
          controls_dprimes['ss_us']])
box_dprime_ptns.set_xticklabels(['SS vs CP','SS vs CS', 'SS vs US'])
box_dprime_ptns.set_xlabel("Stimulus Pair")
box_dprime_ptns.set_ylabel("d prime")
box_dprime_ptns.set_title("Elderly Patients")
box_dprime_ptns.set_ylim(-2.5, 4.5)
# statistical annotation
x1, x2, x3 = 0, 1, 2   # columns (first column: 0, see plt.xticks())
y, h, col = 4.1, 0.05, 'k'
#Significant
plt.plot([x1, x1, x3, x3], [y, y+h, y+h, y], lw=1.5, c=col)
plt.text((x1+x3)*.5, y+h, "*", ha='center', va='bottom', color=col)


plt.subplot(223)
box_dprime_plt = sns.boxplot(
    data=[pilot_dprimes['ss_cp'],
          pilot_dprimes['ss_cs'],
          pilot_dprimes['ss_us']])
box_dprime_plt.set_xticklabels(['SS vs CP','SS vs CS', 'SS vs US'])
box_dprime_plt.set_xlabel("Stimulus Pair")
box_dprime_plt.set_ylabel("d prime")
box_dprime_plt.set_title("Young")
box_dprime_plt.set_ylim(-2.5, 4.5)
# statistical annotation
x1, x2, x3 = 0, 1, 2   # columns (first column: 0, see plt.xticks())
y, h, col = 2.5, 0.05, 'k'
#Significant
plt.plot([x1, x1, x2, x2], [y, y-h, y-h, y], lw=1.5, c=col)
plt.text((x1+x2)*.5, y-h-0.05, "*", ha='center', va='top', color=col)
plt.plot([x1, x1, x3, x3], [y-0.5, y-0.5-h, y-.5-h, y-.5], lw=1.5, c=col)
plt.text((x1+x3)*.5, y-.5-h-0.05, "*", ha='center', va='top', color=col)
plt.plot([x2, x2, x3, x3], [0.5+y, 0.5+y-h, 0.5+y-h, 0.5+y], lw=1.5, c=col)
plt.text((x2+x3)*.5, y+.5-h-0.05, "*", ha='center', va='top', color=col)

plt.subplots_adjust(top=2.5, bottom=0.001, left=0.10, right=2.5, hspace=0.30,
                    wspace=0.25)

print("Pilot medians:")
print("\t SS vs CS:" + str(scipy.median(pilot_dprimes['ss_cs'])))
print("\t SS vs CP:" + str(scipy.median(pilot_dprimes['ss_cp'])))
print("\t SS vs US:" + str(scipy.median(pilot_dprimes['ss_us'])))

print("Controls medians:")
print("\t SS vs CS:" + str(scipy.median(controls_dprimes['ss_cs'])))
print("\t SS vs CP:" + str(scipy.median(controls_dprimes['ss_cp'])))
print("\t SS vs US:" + str(scipy.median(controls_dprimes['ss_us'])))

print("Patient medians:")
print("\t SS vs CS:" + str(scipy.median(patients_dprimes['ss_cs'])))
print("\t SS vs CP:" + str(scipy.median(list(filter(None, patients_dprimes['ss_cp'])))))
print("\t SS vs US:" + str(scipy.median(patients_dprimes['ss_us'])))



    


















    






Pilot medians:
	 SS vs CS:3.96195987
	 SS vs CP:3.33500058186
	 SS vs US:3.96195987
Controls medians:
	 SS vs CS:3.15981319365
	 SS vs CP:2.62354586968
	 SS vs US:3.79554552511
Patient medians:
	 SS vs CS:1.67593772856
	 SS vs CP:1.51009468399
	 SS vs US:2.01398100743

Bootstrap Median Differences

http://faculty.psy.ohio-state.edu/myung/personal/course/826/bootstrap_hypo.pdf

Between groups




In [130]:



    


## Pilot vs Controls d'
dp_pilot_controls = dict()



    











In [131]:



    


stim_key = 'ss_us'
diffs = [x - y for x in pilot_dprimes[stim_key] for y in controls_dprimes[stim_key]]
dp_pilot_controls[stim_key] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(dp_pilot_controls[stim_key]['bootstrap_data'])
print(dp_pilot_controls[stim_key]['CIs'])
print(dp_pilot_controls[stim_key]['central_tendency'])



    


















    






{'low': 0.29413059836906852, 'high': 0.8664930697142128}
{'median': 0.33282868977188995}

















In [132]:



    


stim_key = 'ss_cs'
diffs = [x - y for x in pilot_dprimes[stim_key] for y in controls_dprimes[stim_key]]
dp_pilot_controls[stim_key] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(dp_pilot_controls[stim_key]['bootstrap_data'])
print(dp_pilot_controls[stim_key]['CIs'])
print(dp_pilot_controls[stim_key]['central_tendency'])



    


















    






{'low': 0.60353023970399189, 'high': 1.1993217594861032}
{'median': 0.92108988651002743}

















In [133]:



    


stim_key = 'ss_cp'
diffs = [x - y for x in pilot_dprimes[stim_key] for y in controls_dprimes[stim_key]]
dp_pilot_controls[stim_key] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(dp_pilot_controls[stim_key]['bootstrap_data'])
print(dp_pilot_controls[stim_key]['CIs'])
print(dp_pilot_controls[stim_key]['central_tendency'])



    


















    






{'low': 0.58964897972426977, 'high': 1.4934523578551717}
{'median': 0.94108928758385746}

















In [69]:



    


## Controls vs Patients d'
dp_controls_patients = dict()



    











In [70]:



    


stim_key = 'ss_us'
diffs = [x - y for x in controls_dprimes[stim_key] for y in patients_dprimes[stim_key]]
dp_controls_patients[stim_key] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(dp_controls_patients[stim_key]['bootstrap_data'])
print(dp_controls_patients[stim_key]['CIs'])
print(dp_controls_patients[stim_key]['central_tendency'])



    


















    






{'low': 0.62695928814095891, 'high': 2.1931348314821921}
{'median': 1.3266417617679784}

















In [71]:



    


stim_key = 'ss_cs'
diffs = [x - y for x in controls_dprimes[stim_key] for y in patients_dprimes[stim_key]]
dp_controls_patients[stim_key] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(dp_controls_patients[stim_key]['bootstrap_data'])
print(dp_controls_patients[stim_key]['CIs'])
print(dp_controls_patients[stim_key]['central_tendency'])



    


















    






{'low': 0.95841282815938955, 'high': 2.4374510688381914}
{'median': 1.5315445228373428}

















In [72]:



    


stim_key = 'ss_cp'
diffs = [x - y for x in controls_dprimes[stim_key] 
         for y in list(filter(None, patients_dprimes[stim_key]))]
dp_controls_patients[stim_key] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))

#plt.hist(dp_controls_patients[stim_key]['bootstrap_data'])
print(dp_controls_patients[stim_key]['CIs'])
print(dp_controls_patients[stim_key]['central_tendency'])



    


















    






{'low': 0.55065982450242856, 'high': 1.7758748076109185}
{'median': 1.0744674729460513}
Within Groups




In [144]:



    


# pilot
dp_pilot = dict()

diffs = [x - y for x in pilot_dprimes['ss_cp'] 
         for y in list(filter(None, pilot_dprimes['ss_cs']))]
dp_pilot['cp_cs'] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))
print('cp_cs')
print(dp_pilot['cp_cs']['CIs'])
print(dp_pilot['cp_cs']['central_tendency'])
print('----------------')

diffs = [x - y for x in pilot_dprimes['ss_cp'] 
         for y in list(filter(None, pilot_dprimes['ss_us']))]
dp_pilot['cp_us'] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))
print('cp_us')
print(dp_pilot['cp_us']['CIs'])
print(dp_pilot['cp_us']['central_tendency'])
print('----------------')

diffs = [x - y for x in pilot_dprimes['ss_cs'] 
         for y in list(filter(None, pilot_dprimes['ss_us']))]
dp_pilot['cs_us'] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))
print('cs_us')
print(dp_pilot['cs_us']['CIs'])
print(dp_pilot['cs_us']['central_tendency'])
print('----------------')



    


















    






cp_cs
{'low': -0.62695928814095803, 'high': -0.47541870257483354}
{'median': -0.57236247134514517}
----------------
cp_us
{'low': -0.72314301917913015, 'high': -0.58826119673813793}
{'median': -0.62695928814095891}
----------------
cs_us
{'low': -0.29413059836906941, 'high': -0.015269042965854851}
{'median': -0.29413059836906852}
----------------

















In [145]:



    


# controls
dp_controls = dict()

diffs = [x - y for x in controls_dprimes['ss_cp'] 
         for y in list(filter(None, controls_dprimes['ss_cs']))]
dp_controls['cp_cs'] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))
print('cp_cs')
print(dp_controls['cp_cs']['CIs'])
print(dp_controls['cp_cs']['central_tendency'])
print('----------------')

diffs = [x - y for x in controls_dprimes['ss_cp'] 
         for y in list(filter(None, controls_dprimes['ss_us']))]
dp_controls['cp_us'] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))
print('cp_us')
print(dp_controls['cp_us']['CIs'])
print(dp_controls['cp_us']['central_tendency'])
print('----------------')

diffs = [x - y for x in controls_dprimes['ss_cs'] 
         for y in list(filter(None, controls_dprimes['ss_us']))]
dp_controls['cs_us'] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))
print('cs_us')
print(dp_controls['cs_us']['CIs'])
print(dp_controls['cs_us']['central_tendency'])
print('----------------')



    


















    






cp_cs
{'low': -1.4133221495555461, 'high': 0.10578631622924256}
{'median': -0.72898361090395358}
----------------
cp_us
{'low': -1.7374010524565033, 'high': -0.31688257307629453}
{'median': -1.0995034601610798}
----------------
cs_us
{'low': -1.06897886631943, 'high': 0.064346393365945609}
{'median': -0.62695928814095891}
----------------

















In [148]:



    


# patients
dp_patients = dict()

diffs = [x - y for x in list(filter(None, patients_dprimes['ss_cp'])) 
         for y in list(filter(None, patients_dprimes['ss_cs']))]
dp_patients['cp_cs'] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))
print('cp_cs')
print(dp_patients['cp_cs']['CIs'])
print(dp_patients['cp_cs']['central_tendency'])
print('----------------')

diffs = [x - y for x in list(filter(None, patients_dprimes['ss_cp'])) 
         for y in list(filter(None, patients_dprimes['ss_us']))]
dp_patients['cp_us'] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))
print('cp_us')
print(dp_patients['cp_us']['CIs'])
print(dp_patients['cp_us']['central_tendency'])
print('----------------')

diffs = [x - y for x in list(filter(None, patients_dprimes['ss_cs'])) 
         for y in list(filter(None, patients_dprimes['ss_us']))]
dp_patients['cs_us'] = bayes_bootstrap(diffs, scipy.median, 1000, len(diffs))
print('cs_us')
print(dp_patients['cs_us']['CIs'])
print(dp_patients['cs_us']['central_tendency'])
print('----------------')



    


















    






cp_cs
{'low': -0.67977961758755123, 'high': 0.40964165090354621}
{'median': -0.15224087191336999}
----------------
cp_us
{'low': -1.5371833597567561, 'high': -0.072426190032757853}
{'median': -0.92513281576029494}
----------------
cs_us
{'low': -1.43476843955915, 'high': 0.065079298229978932}
{'median': -0.80776768103767815}
----------------

Diffusion Decision Model Comparisons

Drift Rate Between Stimuli Groups (Within Subjects)





In [73]:



    


pilot_vSS, pilot_vCP, pilot_vCS, pilot_vUS = pilot_model.nodes_db.node[['v(SS)', 'v(CP)', 'v(CS)', 'v(US)']]

hddm.analyze.plot_posterior_nodes([pilot_vSS, pilot_vCP, pilot_vCS, pilot_vUS])

print("Drift Rate Analysis: Pilot")
print('P(v[SS] > v[US]) = ' + str((pilot_vSS.trace() > pilot_vUS.trace()).mean()))
print('P(v[SS] > v[CP]) = ' + str((pilot_vSS.trace() > pilot_vCP.trace()).mean()))
print('P(v[SS] > v[CS]) = ' + str((pilot_vSS.trace() > pilot_vCS.trace()).mean()))
print('P(v[CP] > v[CS]) = ' + str((pilot_vCP.trace() > pilot_vCS.trace()).mean()))
print('P(v[CP] > v[US]) = ' + str((pilot_vCP.trace() > pilot_vUS.trace()).mean()))
print('P(v[CS] > v[US]) = ' + str((pilot_vCS.trace() > pilot_vUS.trace()).mean()))



    


















    






Drift Rate Analysis: Pilot
P(v[SS] > v[US]) = 0.760795454545
P(v[SS] > v[CP]) = 1.0
P(v[SS] > v[CS]) = 0.998181818182
P(v[CP] > v[CS]) = 0.0
P(v[CP] > v[US]) = 0.0
P(v[CS] > v[US]) = 0.0110227272727










    
























In [74]:



    


controls_vSS, controls_vCP, controls_vCS, controls_vUS = controls_model.nodes_db.node[['v(SS)', 'v(CP)', 'v(CS)', 'v(US)']]

hddm.analyze.plot_posterior_nodes([controls_vSS, controls_vCP, controls_vCS, controls_vUS])

print("Drift Rate Analysis: Controls")
print('P(v[SS] > v[US]) = ' + str((controls_vSS.trace() > controls_vUS.trace()).mean()))
print('P(v[SS] > v[CP]) = ' + str((controls_vSS.trace() > controls_vCP.trace()).mean()))
print('P(v[SS] > v[CS]) = ' + str((controls_vSS.trace() > controls_vCS.trace()).mean()))
print('P(v[CP] > v[CS]) = ' + str((controls_vCP.trace() > controls_vCS.trace()).mean()))
print('P(v[CP] > v[US]) = ' + str((controls_vCP.trace() > controls_vUS.trace()).mean()))
print('P(v[CS] > v[US]) = ' + str((controls_vCS.trace() > controls_vUS.trace()).mean()))



    


















    






Drift Rate Analysis: Controls
P(v[SS] > v[US]) = 0.315
P(v[SS] > v[CP]) = 0.989431818182
P(v[SS] > v[CS]) = 0.800909090909
P(v[CP] > v[CS]) = 0.06125
P(v[CP] > v[US]) = 0.00511363636364
P(v[CS] > v[US]) = 0.0948863636364










    
























In [75]:



    


patients_vSS, patients_vCP, patients_vCS, patients_vUS = patients_model.nodes_db.node[['v(SS)', 'v(CP)', 'v(CS)', 'v(US)']]

hddm.analyze.plot_posterior_nodes([patients_vSS, patients_vCP, patients_vCS, patients_vUS])

print("Drift Rate Analysis: Patients")
print('P(v[SS] > v[US]) = ' + str((patients_vSS.trace() > patients_vUS.trace()).mean()))
print('P(v[SS] > v[CP]) = ' + str((patients_vSS.trace() > patients_vCP.trace()).mean()))
print('P(v[SS] > v[CS]) = ' + str((patients_vSS.trace() > patients_vCS.trace()).mean()))
print('P(v[CP] > v[CS]) = ' + str((patients_vCP.trace() > patients_vCS.trace()).mean()))
print('P(v[CP] > v[US]) = ' + str((patients_vCP.trace() > patients_vUS.trace()).mean()))
print('P(v[CS] > v[US]) = ' + str((patients_vCS.trace() > patients_vUS.trace()).mean()))



    


















    






Drift Rate Analysis: Patients
P(v[SS] > v[US]) = 0.922045454545
P(v[SS] > v[CP]) = 0.983068181818
P(v[SS] > v[CS]) = 0.984659090909
P(v[CP] > v[CS]) = 0.486136363636
P(v[CP] > v[US]) = 0.205227272727
P(v[CS] > v[US]) = 0.209318181818

Drift Rate Between Subjects & Stimuli Groups





In [76]:



    


pilot_vSS = pilot_model.nodes_db.node['v(SS)']
controls_vSS = controls_model.nodes_db.node['v(SS)']
patients_vSS = patients_model.nodes_db.node['v(SS)']

hddm.analyze.plot_posterior_nodes([pilot_vSS, controls_vSS, patients_vSS])

print("Drift Rate Analysis: SS")
print('P(patients_vSS > controls_vSS) = ' + str((patients_vSS.trace() > controls_vSS.trace()).mean()))
print('P(controls_vSS > pilot_vSS) = ' + str((controls_vSS.trace() > pilot_vSS.trace()).mean()))



    


















    






Drift Rate Analysis: SS
P(patients_vSS > controls_vSS) = 0.145227272727
P(controls_vSS > pilot_vSS) = 0.00181818181818










    
























In [77]:



    


pilot_vUS = pilot_model.nodes_db.node['v(US)']
controls_vUS = controls_model.nodes_db.node['v(US)']
patients_vUS = patients_model.nodes_db.node['v(US)']

hddm.analyze.plot_posterior_nodes([pilot_vUS, controls_vUS, patients_vUS])

print("Drift Rate Analysis: US")
print('P(patients_vUS > controls_vUS) = ' + str((patients_vUS.trace() > controls_vUS.trace()).mean()))
print('P(controls_vUS > pilot_vUS) = ' + str((controls_vUS.trace() > pilot_vUS.trace()).mean()))



    


















    






Drift Rate Analysis: US
P(patients_vUS > controls_vUS) = 0.00340909090909
P(controls_vUS > pilot_vUS) = 0.0201136363636










    
























In [78]:



    


pilot_vCS = pilot_model.nodes_db.node['v(CS)']
controls_vCS = controls_model.nodes_db.node['v(CS)']
patients_vCS = patients_model.nodes_db.node['v(CS)']

hddm.analyze.plot_posterior_nodes([pilot_vCS, controls_vCS, patients_vCS])

print("Drift Rate Analysis: CS")
print('P(patients_vCS > controls_vCS) = ' + str((patients_vCS.trace() > controls_vCS.trace()).mean()))
print('P(controls_vCS > pilot_vCS) = ' + str((controls_vCS.trace() > pilot_vCS.trace()).mean()))



    


















    






Drift Rate Analysis: CS
P(patients_vCS > controls_vCS) = 0.00670454545455
P(controls_vCS > pilot_vCS) = 0.00352272727273










    
























In [79]:



    


pilot_vCP = pilot_model.nodes_db.node['v(CP)']
controls_vCP = controls_model.nodes_db.node['v(CP)']
patients_vCP = patients_model.nodes_db.node['v(CP)']

hddm.analyze.plot_posterior_nodes([pilot_vCP, controls_vCP, patients_vCP])

print("Drift Rate Analysis: CP")
print('P(patients_vCP > controls_vCP) = ' + str((patients_vCP.trace() > controls_vCP.trace()).mean()))
print('P(controls_vCP > pilot_vCP) = ' + str((controls_vCP.trace() > pilot_vCP.trace()).mean()))



    


















    






Drift Rate Analysis: CP
P(patients_vCP > controls_vCP) = 0.147954545455
P(controls_vCP > pilot_vCP) = 0.00181818181818

Non-Decision Time Between Subjects





In [80]:



    


pilot_t = pilot_model.nodes_db.node['t']
controls_t = controls_model.nodes_db.node['t']
patients_t = patients_model.nodes_db.node['t']

hddm.analyze.plot_posterior_nodes([pilot_t, controls_t, patients_t])

print("Non-Decision Time Analysis")
print('P(patients_t > controls_t) = ' + str((patients_t.trace() > controls_t.trace()).mean()))
print('P(controls_t > pilot_t) = ' + str((controls_t.trace() > pilot_t.trace()).mean()))



    


















    






Non-Decision Time Analysis
P(patients_t > controls_t) = 0.126704545455
P(controls_t > pilot_t) = 0.909886363636

Bias Between Subjects





In [81]:



    


pilot_z = pilot_model.nodes_db.node['z']
controls_z = controls_model.nodes_db.node['z']
patients_z = patients_model.nodes_db.node['z']

hddm.analyze.plot_posterior_nodes([pilot_z, controls_z, patients_z])

print("Bias Analysis")
print('P(patients_z > controls_z) = ' + str((patients_z.trace() > controls_z.trace()).mean()))
print('P(controls_z > pilot_z) = ' + str((controls_z.trace() > pilot_z.trace()).mean()))



    


















    






Bias Analysis
P(patients_z > controls_z) = 0.000227272727273
P(controls_z > pilot_z) = 0.996590909091

Threshold

Between Subjects

Tested from model where stimulus type does not affect threshold





In [82]:



    


pilot_a = pilot_model.nodes_db.node['a']
controls_a = controls_model.nodes_db.node['a']
patients_a = patients_model.nodes_db.node['a']

hddm.analyze.plot_posterior_nodes([pilot_a, controls_a, patients_a])

print("Threshold Analysis: Between-Subjects")
print('P(patients_a > controls_a) = ' + str((patients_a.trace() > controls_a.trace()).mean()))
print('P(controls_a > pilot_a) = ' + str((controls_a.trace() > pilot_a.trace()).mean()))



    


















    






Threshold Analysis: Between-Subjects
P(patients_a > controls_a) = 0.998522727273
P(controls_a > pilot_a) = 0.0260227272727

Model Deviance Information Criterion (DIC)

Measure of model fitness for hierarchical models, lower is better





In [83]:



    


# Pilot models
print("Stimulus Model DIC: " + str(pilot_model.dic))
print("Lumped Model DIC: " + str(pilot_model_lumped.dic))
#print("Stimulus+Threshold Model DIC: " + str(pilot_model_threshold.dic))



    


















    






Stimulus Model DIC: 795.609087068
Lumped Model DIC: 1016.28785087

















In [84]:



    


# Control models
print("Stimulus Model DIC: " + str(controls_model.dic))
print("Lumped Model DIC: " + str(controls_model_lumped.dic))
#print("Stimulus+Threshold Model DIC: " + str(controls_model_threshold.dic))



    


















    






Stimulus Model DIC: 800.516022989
Lumped Model DIC: 869.320751214

















In [85]:



    


# Patient models
print("Stimulus Model DIC: " + str(patients_model.dic))
print("Lumped Model DIC: " + str(patients_model_lumped.dic))
#print("Stimulus+Threshold Model DIC: " + str(patients_model_threshold.dic))



    


















    






Stimulus Model DIC: 2619.7750154
Lumped Model DIC: 3222.23689013

Content source: CPernet/LanguageDecision

Similar notebooks:

notebook.community | gallery | about