<jupyter_start><jupyter_text>These analyses and simulations use data from Urai et al. 2017, and aim to further pinpoint the relationship between several proxies for confidence on the previous trial.If you use this work in any way, or see something interesting, please contact me first (<EMAIL>)! I'm very happy to chat with people about collaborating to better understand these patterns. Thanks.<jupyter_code>import numpy as np
import pandas as pd
import matplotlib
#matplotlib.use('TkAgg')
import matplotlib.pylab as plt
import seaborn as sns
sns.set(style="darkgrid", context="talk")
import scipy as sp
from IPython import embed as shell
import psychofit as psy
# import tqdm # progress bar
#### Urai et al. (2018): pupil dilation and reaction times scale with decision uncertainty ####
# try:
#     data_urai2017 = pd.read_csv('https://ndownloader.figshare.com/files/12753785')
#     data_urai2017['choice']        = data_urai2017.resp
#     data_urai2017['evidence']      = data_urai2017.coherence * data_urai2017.stimulus # signed evidence
#     data_urai2017['subj_idx']      = data_urai2017.subjnr

# except:
data_urai2017 = pd.read_csv('/Users/urai/Data/HDDM/NatComm/2ifc_data_hddm.csv')
data_urai2017['choice']        = np.sign(data_urai2017.response - 0.1)
data_urai2017['evidence']      = data_urai2017.stimulus
data_urai2017['sessionnr']     = data_urai2017.session
data_urai2017['blocknr']       = data_urai2017.block
data_urai2017['correct']       = 1 * (np.sign(data_urai2017.stimulus) == data_urai2017.choice)
data_urai2017['decision_pupil'] = data_urai2017.prevpupil.shift(-1)

data_urai2017['coherence'].values[data_urai2017['coherence'] <= 0.025] =  0.025
data_urai2017['coherence'] = data_urai2017['coherence'] * 100
data_urai2017['evidence'] = data_urai2017['coherence'] * np.sign(data_urai2017.stimulus)
# normalize RT for each participant
zscore = lambda x: (x - x.mean()) / x.std()
data_urai2017['rt_norm'] = data_urai2017.groupby(['subj_idx', 'sessionnr', 'blocknr'])['rt'].transform(zscore)
data_urai2017['evidence_norm'] = data_urai2017.groupby(['subj_idx', 'sessionnr', 'blocknr'])['evidence'].transform(zscore)
data_urai2017['coherence_norm'] = data_urai2017.groupby(['subj_idx', 'sessionnr', 'blocknr'])['coherence'].transform(zscore)
data_urai2017['decision_pupil'] = data_urai2017.groupby(['subj_idx', 'sessionnr', 'blocknr'])['decision_pupil'].transform(zscore)

data_urai2017.sample(n=10)
## FIRST, MAKE A QUICK OVERVIEW OF BASIC BEHAVIORAL DATA
fig, ax = plt.subplots(1, 3, figsize=(13,4))
sns.lineplot(data=data_urai2017, x='evidence', y='choice', color='k', marker='o', ci=None,  ax=ax[0])
ax[0].set(title='Psychometric function')
sns.regplot(data=data_urai2017, x='coherence', y='correct', x_bins=5, color='k', marker='o', ci=None, ax=ax[1])
ax[1].set(ylim=[0.5, 1])
sns.lineplot(data=data_urai2017, x='evidence', y='rt', color='k', marker='o', ci=None,ax=ax[2])
ax[2].set(title='Chronometric function')
fig.suptitle('Figure 1. Basic behavior in the 2IFC task', verticalalignment='baseline')
plt.tight_layout()
## RT reflects uncertainty
data_urai2017.loc[:,'rt_bin'] = pd.qcut(data_urai2017['rt_norm'], 2, labels=False) # median split

fig, ax = plt.subplots(1, 3, figsize=(13,4))
sns.lineplot(data=data_urai2017, x='evidence', y='rt_norm', hue='correct', marker='o', ax=ax[0], ci=None, palette={0:"firebrick", 1:"forestgreen"}, legend=False)
ax[0].set(title='Vevaiometric curve')
sns.regplot(data=data_urai2017, x='rt_norm', y='correct', x_bins=5, color='k', marker='o', ax=ax[1], ci=None)
ax[1].set(title='Confidence calibration')
sns.lineplot(data=data_urai2017, x='coherence', y='correct', hue='rt_bin', marker='o',ci=None, ax=ax[2], palette={0:"dimgrey", 1:"black"}, legend=False)
ax[2].set(title='Conditional psychometric', ylim=[0.5, 1])
fig.suptitle('Figure 2. Reaction time reflects uncertainty', verticalalignment='baseline')
plt.tight_layout()
# fig.savefig('rt.pdf')
## Pupil reflects uncertainty
data_urai2017.loc[:,'pupil_bin'] = pd.qcut(data_urai2017['decision_pupil'], 2, labels=False) # median split

fig, ax = plt.subplots(1, 3, figsize=(13,4))
sns.lineplot(data=data_urai2017, x='evidence', y='decision_pupil', hue='correct', marker='o', ci=None, ax=ax[0], palette={0:"firebrick", 1:"forestgreen"}, legend=False)
ax[0].set(title='Vevaiometric curve')
sns.regplot(data=data_urai2017, x='decision_pupil', y='correct', x_bins=5, color='k', marker='o', ax=ax[1], ci=None)
ax[1].set(title='Confidence calibration')
sns.lineplot(data=data_urai2017, x='coherence', y='correct', hue='pupil_bin', marker='o', ci=None, ax=ax[2], palette={0:"dimgrey", 1:"black"}, legend=False)
ax[2].set(title='Conditional psychometric', ylim=[0.5, 1])
fig.suptitle('Figure 3. Post-decision pupil dilation reflects uncertainty', verticalalignment='baseline')
plt.tight_layout()
#fig.savefig('pupil.pdf')
# RATHER THAN FITTING PSYCHOMETRIC FUNCTION SHIFTS, DO A QUICK HEURISTIC: REPETITION PROBABILITY
data_urai2017['repeat']              = 1 * (data_urai2017['choice'] == data_urai2017.choice.shift(1))
data_urai2017['previous_coherence']  = data_urai2017.coherence.shift(1)
data_urai2017['previous_choice']     = data_urai2017.choice.shift(1)
data_urai2017['previous_correct']    = data_urai2017.correct.shift(1)
data_urai2017['previous_rt_bins']    = pd.qcut(data_urai2017.rt_norm.shift(1), 3, labels=False) 
data_urai2017['previous_rt']         = data_urai2017.rt_norm.shift(1) 
data_urai2017['previous_pupil_bins'] = pd.qcut(data_urai2017.decision_pupil.shift(1), 3, labels=False) 
data_urai2017['previous_pupil']      = data_urai2017.decision_pupil.shift(1) 

# data_urai2017.loc[:,'current_evidence'] = pd.qcut(data_urai2017['coherence'], 2, labels=False) # median split
## DEFINE UPDATING BASED ON THE PSYCHOMETRIC FUNCTION SHIFT
# from ibl_pipeline.analyses import behavior as behavioral_analyses

def fit_psychfunc(df):
    
    choicedat = df.dropna(subset=['choice']) # ignore missed trials
    choicedat = choicedat.groupby('evidence').agg({'rt':'count', 'response':'mean'}).reset_index()  
    assert(all(choicedat['response'] >= 0)) # make sure we're not using signed choices
    
    if sum(choicedat.rt) > 20: # only fit if each datapoint has at least 10 trials
        pars, L = psy.mle_fit_psycho(choicedat.values.transpose(), P_model='erf_psycho', 
            parstart=np.array([choicedat['evidence'].mean(), 20., 0.05]), 
            parmin=np.array([choicedat['evidence'].min(), 0., 0.]), 
            parmax=np.array([choicedat['evidence'].max(), 100., 1.]))
        df2 = {'bias':pars[0],'threshold':pars[1], 'lapse':pars[2]}
    else:
        df2 = {'bias':np.nan, 'threshold':np.nan, 'lapse':np.nan}
    return pd.DataFrame(df2, index=[0])

def updating(df, colname):
    
    # fit psychometric functions
    psychfuncfits = df.groupby(['subj_idx', 'previous_choice', 'previous_correct', colname]).apply(fit_psychfunc).reset_index()    
    
    # make a new dataframe to compute the bias shift per group
    df2 = pd.DataFrame()
    for name, grouped in psychfuncfits.groupby(['previous_correct', 'subj_idx', colname]):
        if grouped.shape[0] == 2: # skip if there aren't 2 previous responses in this bin
            biasshift = grouped.loc[grouped.previous_choice == -1, 'bias'].item() - grouped.loc[grouped.previous_choice == 1, 'bias'].item()
            df3 = grouped.groupby(['previous_correct', 'subj_idx', colname]).mean().reset_index() # average the rest
            df3['update'] = biasshift # now insert the computed biasshift
            df2 = df2.append(df3)
            
    print(colname)
    df2 = df2.dropna(subset=['update'])
    # print(df2.columns)
    return df2

## NOW FIT THE PSYCHOMETRIC FUNCTION SHIFT, DEPENDING ON PREVIOUS COHERENCE
coherence_update = updating(data_urai2017, 'previous_coherence')
pupil_update     = updating(data_urai2017, 'previous_pupil_bins')
rt_update        = updating(data_urai2017, 'previous_rt_bins')
## HOW DO PUPIL, RT AND COHERENCE AFFECT NEXT-TRIAL UPDATING?
fig, ax = plt.subplots(1, 3, figsize=(13,4), sharey=True)
kwargs = {'y':'update', 'ci':68, 'marker':'o', 'estimator':np.nanmean, 'err_style':"bars"}
sns.lineplot(data=coherence_update[coherence_update.previous_correct == 1], x='previous_coherence',  ax=ax[0], **kwargs)
sns.lineplot(data=rt_update[rt_update.previous_correct == 1], x='previous_rt_bins', ax=ax[1], **kwargs)
sns.lineplot(data=pupil_update[pupil_update.previous_correct == 1], x='previous_pupil_bins', ax=ax[2], **kwargs)
fig.suptitle('Figure 4. Decision uncertainty reduces repetition', verticalalignment='baseline')
plt.tight_layout()
fig.savefig('updating.pdf')
## DO THE RT-BINNING PER LEVEL OF PREVIOUS COHERENCE
def updating_2split(df, colname1, colname2):
    
    # fit psychometric functions
    psychfuncfits = df.groupby(['subj_idx', 'previous_choice', 'previous_correct', colname1, colname2]).apply(fit_psychfunc).reset_index()    
    
    # make a new dataframe to compute the bias shift per group
    df2 = pd.DataFrame()
    for name, grouped in psychfuncfits.groupby(['previous_correct', 'subj_idx', colname1, colname2]):
        if grouped.shape[0] == 2: # skip if there aren't 2 previous responses in this bin
            biasshift = grouped.loc[grouped.previous_choice == -1, 'bias'].item() - grouped.loc[grouped.previous_choice == 1, 'bias'].item()
            df3 = grouped.groupby(['previous_correct', 'subj_idx', colname1, colname2]).mean().reset_index() # average the rest
            df3['update'] = biasshift # now insert the computed biasshift
            df2 = df2.append(df3)
            
    df2 = df2.dropna(subset=['update'])
    print(colname1)
    print(colname2)
    # print(df2.columns)
    return df2
# NOW DO THE SPLIT INTO PUPIL OR RT BINS, BUT PER LEVEL OF PREVIOUS COHERENCE
def shiftrt_percoh(x):
    return pd.qcut(x, 4, labels=False)
    
data_urai2017['previous_rt_bins_percoh'] = data_urai2017.groupby(['previous_coherence'])['previous_rt'].apply(shiftrt_percoh)
rt_update_percoh = updating_2split(data_urai2017, 'previous_rt_bins_percoh', 'previous_coherence')

data_urai2017['previous_pupil_bins_percoh'] = data_urai2017.groupby(['previous_coherence'])['previous_pupil'].apply(shiftrt_percoh)
pupil_update_percoh = updating_2split(data_urai2017, 'previous_pupil_bins_percoh', 'previous_coherence')
rt_update_percoh.head(n=10)
fig, ax = plt.subplots(1, 2, figsize=(13,4), sharey=True)
sns.lineplot(data=rt_update_percoh[rt_update_percoh.previous_correct == 1], x='previous_rt_bins_percoh', 
             hue='previous_coherence', legend=False, ax=ax[0], palette='viridis', **kwargs)
sns.lineplot(data=pupil_update_percoh[pupil_update_percoh.previous_correct == 1], x='previous_pupil_bins_percoh', 
             hue='previous_coherence', legend='brief', palette='viridis',  ax=ax[1], **kwargs)
lgd = ax[1].legend(loc='center right', bbox_to_anchor=(1.75, 0.7), ncol=1) # move box outside
fig.suptitle('Figure 5. More updating with lower previous contrast, and also with lower previous decision uncertainty')
fig.savefig('updating_percoh.pdf')
sns.lmplot(data=data_urai2017, x='rt', y='correct', hue='coherence', 
           x_bins=10, palette='viridis', truncate=True, ci=None)<jupyter_output><empty_output><jupyter_text>* Stronger coherence -> more correct / more confidence* Within each coherence level, fast RT -> more correct / more confidence*BUT** Stronger coherence -> less updating on next trial* Within each coherence level, fast RT -> more updating on next trial<jupyter_code>## test anova
from statsmodels.formula.api import ols
formula = 'update ~ previous_coherence*previous_rt_bins_percoh'
lm = ols(formula, rt_update_percoh).fit()
print(lm.summary())

## test anova
formula = 'update ~ previous_coherence*previous_pupil_bins_percoh'
lm = ols(formula, pupil_update_percoh).fit()
print(lm.summary())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:                 update   R-squared:                       0.042
Model:                            OLS   Adj. R-squared:                  0.037
Method:                 Least Squares   F-statistic:                     8.103
Date:                Thu, 06 Jun 2019   Prob (F-statistic):           2.73e-05
Time:                        15:01:24   Log-Likelihood:                -1731.7
No. Observations:                 561   AIC:                             3471.
Df Residuals:                     557   BIC:                             3489.
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
=================================================================================================================
                 [...]<jupyter_text>**Conclusion**: the weaker the previous coherence, the larger the choice update (i.e. more information about the bound placement). However, _within_ each level of coherence, faster RTs (so a more confident choice) also lead to stronger updating. Even when ignoring previous coherence, and just grouping by RT, this captures the<jupyter_code>### CAN ALSO USE STATSMODELS!
import statsmodels.api as sm
import statsmodels.formula.api as smf

data_urai2017 = data_urai2017.dropna(subset=['previous_coherence', 'previous_rt', 'previous_pupil'])
designM = data_urai2017[['evidence', 'previous_coherence', 'previous_choice']]
designM['prev_coh_resp'] = designM.previous_choice * designM.previous_coherence
designM = sm.add_constant(designM, prepend=True)

logit_mod = sm.Logit(data_urai2017.response, designM)
logit_res = logit_mod.fit(disp=0)
print(logit_res.summary())

# md = smf.mixedlm("response ~ evidence", data_urai2017, groups=data_urai2017["subj_idx"], family=sm.families.Binomial())
# mdf = data_urai2017.fit()
# print(mdf.summary())
### CAN ALSO USE STATSMODELS!
import statsmodels.api as sm
import statsmodels.formula.api as smf

designM = data_urai2017[['evidence', 'previous_rt', 'previous_choice']]
designM['prev_rt_resp'] = designM.previous_choice * designM.previous_rt
designM = sm.add_constant(designM, prepend=True)

logit_mod = sm.Logit(data_urai2017.response, designM)
logit_res = logit_mod.fit(disp=0)
print(logit_res.summary())
### CAN ALSO USE STATSMODELS!
import statsmodels.api as sm
import statsmodels.formula.api as smf

designM = data_urai2017[['evidence', 'previous_pupil', 'previous_choice']]
designM['prev_pupil_resp'] = designM.previous_choice * designM.previous_pupil
designM = sm.add_constant(designM, prepend=True)

logit_mod = sm.Logit(data_urai2017.response, designM)
logit_res = logit_mod.fit(disp=0)
print(logit_res.summary())
## NOW COMBINE ALL 3 IN ONE MODEL..

designM = data_urai2017[['evidence', 'previous_coherence', 'previous_rt', 'previous_pupil', 'previous_choice']]
designM['prev_coh_resp'] = designM.previous_choice * designM.previous_coherence
designM['prev_rt_resp'] = designM.previous_choice * designM.previous_rt
designM['prev_pupil_resp'] = designM.previous_choice * designM.previous_pupil
designM = sm.add_constant(designM, prepend=True)

logit_mod = sm.Logit(data_urai2017.response, designM)
logit_res = logit_mod.fit(disp=0)
print(logit_res.summary())<jupyter_output>Logit Regression Results                           
==============================================================================
Dep. Variable:               response   No. Observations:                64343
Model:                          Logit   Df Residuals:                    64334
Method:                           MLE   Df Model:                            8
Date:                Sun, 02 Jun 2019   Pseudo R-squ.:                  0.2950
Time:                        16:34:59   Log-Likelihood:                -31440.
converged:                       True   LL-Null:                       -44595.
                                        LLR p-value:                     0.000
======================================================================================
                         coef    std err          z      P>|z|      [0.025      0.975]
--------------------------------------------------------------------------------------
const                 -0.014[...]
{'fraction_non_alphanumeric': 0.1607131840335616, 'fraction_numerical': 0.042383860818064036, 'mean_word_length': 4.273259596616786, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 19, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '240318', 'license': 'MIT', 'n_tokens_mistral': 6006, 'n_tokens_neox': 5193, 'n_words': 1091, 'path': 'newstuff/.ipynb_checkpoints/confidence_updating_4Armin-checkpoint.ipynb', 'repo_name': 'anne-urai/pupilUncertainty', 'size': '16221'}