# 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

# 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')

# 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')

# 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')

# 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))
        }
    }

# 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]
        }
    }

"""
Explanation: A Diffusion Decision Model for Language-based Decisions
<img src="images/ddm_desc.jpeg" width="600">
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)
End of explanation
"""


# 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']

"""
Explanation: Performance Comparison
End of explanation
"""


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)))

rt_comps = dict()

# 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('--------------------')

"""
Explanation: Reaction Time
End of explanation
"""


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

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'])

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'])

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'])

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'])

"""
Explanation: Bootstrap Median Differences Between Subject Groups
Pilot vs Controls
End of explanation
"""


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

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'])

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'])

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'])

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'])

"""
Explanation: Controls vs Patients
End of explanation
"""


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

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'])

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'])

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'])

"""
Explanation: Bootstrap Median Differences Between Stimuli Groups (Incogruency Effect Only)
Pilot
End of explanation
"""


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

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'])

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'])

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'])

"""
Explanation: Controls
End of explanation
"""


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

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'])

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'])

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'])

"""
Explanation: Patients
End of explanation
"""


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

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)

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)))

acc_comps = dict()

# 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('--------------------')

"""
Explanation: Accuracy
End of explanation
"""


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

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'])

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'])

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'])

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'])

"""
Explanation: Bootstrap Median Differences Between Subject Groups
Pilot vs Controls
End of explanation
"""


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

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'])

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'])

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'])

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'])

"""
Explanation: Controls vs Patients
End of explanation
"""


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

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'])

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'])

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'])

"""
Explanation: Bootstrap Median Differences Between Stimuli Categories (Within Subject Groups, Incogruency Only)
Pilot
End of explanation
"""


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

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'])

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'])

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'])

"""
Explanation: Controls
End of explanation
"""


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

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'])

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'])

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'])

"""
Explanation: Patients
End of explanation
"""


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

## 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')
}

## 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'])))

"""
Explanation: 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
End of explanation
"""


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

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'])

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'])

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'])

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

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'])

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'])

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'])

"""
Explanation: Bootstrap Median Differences
http://faculty.psy.ohio-state.edu/myung/personal/course/826/bootstrap_hypo.pdf
Between groups
End of explanation
"""


# 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('----------------')

# 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('----------------')

# 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('----------------')

"""
Explanation: Within Groups
End of explanation
"""


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()))

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()))

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()))

"""
Explanation: Diffusion Decision Model Comparisons
Drift Rate Between Stimuli Groups (Within Subjects)
End of explanation
"""


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()))

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()))

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()))

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()))

"""
Explanation: Drift Rate Between Subjects & Stimuli Groups
End of explanation
"""


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()))

"""
Explanation: Non-Decision Time Between Subjects
End of explanation
"""


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()))

"""
Explanation: Bias Between Subjects
End of explanation
"""


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()))

"""
Explanation: Threshold
Between Subjects
Tested from model where stimulus type does not affect threshold
End of explanation
"""


# 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))

# 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))

# 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))

"""
Explanation: Model Deviance Information Criterion (DIC)
Measure of model fitness for hierarchical models, lower is better
End of explanation
"""


