<jupyter_start><jupyter_text>Model 1b<jupyter_code>import itertools
import pyfir
from nipy.modalities.fmri import design_matrix, experimental_paradigm
import statsmodels.api as sm

level2_dict = []
adj_r = []
# pars = ['neutral_cue', 'left_cue', 'right_cue', 'neutral_easy_correct', 'valid_easy_correct', 'invalid_easy_correct', 'neutral_hard_correct', 'valid_hard_correct', 'invalid_hard_correct', 'easy_error', 'hard_error',]
# pars = 

mask = 'STh_L'

for subject_id in np.arange(1, 20):
    for run in [1,2,3]:

        try:
            onsets_dict = get_onsets(subject_id, run, 'model1b', shift=-3)
            pars = onsets_dict.keys()

            conditions = []
            onsets = []
            for key in onsets_dict.keys():
                conditions += [key] * len(onsets_dict[key])
                onsets += onsets_dict[key]

            paradigm =  experimental_paradigm.BlockParadigm(con_id=conditions, 
                                                            onset=onsets,
                                                            duration=[[1.]] * len(conditions))
            tr = 3.0
            frametimes = np.linspace(0, 387 * tr, 387)
            hrf_model = 'Canonical'

            X, names= design_matrix.dmtx_light(frametimes, paradigm, hrf_model=hrf_model, drift_model='polynomial', fir_delays=np.arange(11), drift_order=3)

            X /= X.max(0)
            X = pandas.DataFrame(X, columns=names)

            for i, mask in enumerate(['STh_L', 'STh_R']):
                signal = get_signal(subject_id,  run, mask)
                r = sm.GLSAR(signal, X).iterative_fit()

                adj_r.append({'subject_id':subject_id, 'mask':mask, 'r2_adj':r.rsquared_adj, 'r2':r.rsquared, 'run':run})

                for par in pars:
                    level2_dict.append({'value':r.params[par],
                                        'par':par,
                                        'subject_id':subject_id,
                                        'mask':mask,
                                        'run':run})
        except Exception as e:
            print 'Error: %s/%s/%s' %  (subject_id, mask, e)
df_glm = pandas.DataFrame(level2_dict)
df_glm = df_glm.query('par != "constant"')
df_glm['dummy'] = 1

df_glm['% signal change'] = df_glm['value']
sns.set_palette('Set1')
sns.set_context('poster')
sns.set_style('white')

# sns.factorplot(x='dummy', y='value', hue='par', data=df_glm, col='mask', col_wrap=3, kind='bar', ci=95, hue_order=['payoff_cue (left)', 'payoff_cue (right)', 'neutral_cue','difficult_rdm', 'easy_rdm'])
sns.factorplot(x='dummy', y='value', hue='par', data=df_glm, col='mask', col_wrap=3, kind='bar', ci=95)
shifts = pandas.read_pickle('/home/gdholla1/projects/bias/data/derivatives/behavior/hddm/bias_shifts.pkl')
shifts = shifts.set_index('subj_idx')
tmp = df_glm.pivot_table(index='subject_id', values='value', columns=['par', 'mask'])
tmp = tmp['left'] - tmp['neutral']
shifts = shifts.merge(tmp, left_index=True, right_index=True)

shifts = pandas.read_pickle('/home/gdholla1/projects/bias/data/derivatives/behavior/hddm/bias_shifts.pkl')
shifts = shifts.set_index('subj_idx')
tmp = df_glm.pivot_table(index='subject_id', values='value', columns=['par', 'mask'])
tmp = tmp['left'] - tmp['neutral']
shifts = shifts.merge(tmp, left_index=True, right_index=True)

results = []

for mask in masks:
    sns.lmplot('z_cue_coding', mask, shifts)
    
    r, p = sp.stats.pearsonr(shifts['z_cue_coding'], shifts[mask])
    
    results.append({'mask':mask, 'r':r, 'p':p})
    
results = pandas.DataFrame(results)
_, p_fdr, _, _ = statsmodels.sandbox.stats.multicomp.multipletests(results['p'], method='fdr_bh')
results['p.fdr'] = p_fdr
results
%%R -i shifts
library(BayesFactor)

bf = regressionBF(z_cue_coding ~ 1 + STh_L, data = shifts, whichModels='all')
print(bf)

bf = regressionBF(z_cue_coding ~ 1 + STh_R, data = shifts, whichModels='all')
print(bf)<jupyter_output><empty_output><jupyter_text>Difficulty  Include errors<jupyter_code>import itertools
import pyfir
from nipy.modalities.fmri import design_matrix, experimental_paradigm
import statsmodels.api as sm

level2_dict = []
adj_r = []
# pars = ['neutral_cue', 'left_cue', 'right_cue', 'neutral_easy_correct', 'valid_easy_correct', 'invalid_easy_correct', 'neutral_hard_correct', 'valid_hard_correct', 'invalid_hard_correct', 'easy_error', 'hard_error',]
# pars = 

mask = 'STh_L'

for subject_id in np.arange(1, 20):
    for run in [1,2,3]:

        try:
            onsets_dict = get_onsets(subject_id, run, 'model1a', shift=-3)
            pars = onsets_dict.keys()

            conditions = []
            onsets = []
            for key in onsets_dict.keys():
                conditions += [key] * len(onsets_dict[key])
                onsets += onsets_dict[key]

            paradigm =  experimental_paradigm.BlockParadigm(con_id=conditions, 
                                                            onset=onsets,
                                                            duration=[[1.]] * len(conditions))
            tr = 3.0
            frametimes = np.linspace(0, 387 * tr, 387)
            hrf_model = 'Canonical with Derivative'

            X, names= design_matrix.dmtx_light(frametimes, paradigm, hrf_model=hrf_model, drift_model='polynomial', fir_delays=np.arange(11), drift_order=3)

            X /= X.max(0)
            X = pandas.DataFrame(X, columns=names)

            for i, mask in enumerate(['STh_L', 'STh_R']):
                signal = get_signal(subject_id,  run, mask)
                r = sm.GLSAR(signal, X).iterative_fit()

                adj_r.append({'subject_id':subject_id, 'mask':mask, 'r2_adj':r.rsquared_adj, 'r2':r.rsquared, 'run':run})

                for par in pars:
                    level2_dict.append({'value':r.params[par],
                                        'par':par,
                                        'subject_id':subject_id,
                                        'mask':mask,
                                        'run':run})
        except Exception as e:
            print 'Error: %s/%s/%s' %  (subject_id, mask, e)
df_glm = pandas.DataFrame(level2_dict)
df_glm = df_glm.query('par != "constant"')
df_glm['dummy'] = 1

df_glm['% signal change'] = df_glm['value']
sns.set_palette('Set1')
sns.set_context('poster')
sns.set_style('white')

# sns.factorplot(x='dummy', y='value', hue='par', data=df_glm, col='mask', col_wrap=3, kind='bar', ci=95, hue_order=['payoff_cue (left)', 'payoff_cue (right)', 'neutral_cue','difficult_rdm', 'easy_rdm'])
sns.factorplot(x='dummy', y='value', hue='par', data=df_glm, col='mask', col_wrap=3, kind='bar', ci=95)

df_glm_include_errors = df_glm.copy()
df_glm_include_errors['include errors'] = True
shifts = pandas.read_pickle('/home/gdholla1/projects/bias/data/derivatives/behavior/hddm/bias_shifts.pkl')
shifts = shifts.set_index('subj_idx')
tmp = df_glm.pivot_table(index='subject_id', values='value', columns=['par', 'mask'])

tmp = tmp['difficult_rdm'] - tmp['easy_rdm']
shifts = shifts.merge(tmp, left_index=True, right_index=True)
tmp = pandas.melt(shifts, id_vars=['v_cue_coding', 'v_difficulty', 'z_cue_coding'], var_name='mask')

tmp['mask'] = tmp['mask'].apply(lambda x: x.replace('STh', 'STN'))

# fac = sns.FacetGrid(tmp, col='mask', col_wrap=3)

# # def lmplot(x, y, color, **kwargs):
# #     sns.lmplot(x, y, **kwargs)

def make_r_title(x, y, **kwargs):
    r, p = sp.stats.pearsonr(x, y)
    
    title = plt.gca().get_title()
    
    ax = plt.gca()

    y = 0.85 - len(ax.texts) * 0.05
    ax.annotate('r = %.2f' % (r), xy=(0.05,y), xycoords='axes fraction')


sns.set_palette(sns.color_palette('Set1')[1:])

tmp['v_difficulty'] *= -1

fac = sns.lmplot('v_difficulty', 'value', data=tmp, col='mask', col_wrap=3)
fac.map(make_r_title, 'v_difficulty', 'value', )
# fac.set_titles('{col_name}')
fac.set_xlabels('Drift rate reduction')
fac.set_ylabels('hard RDM - easy RDM\n(% signal change)')

results = []

for mask in masks:
    r, p = sp.stats.pearsonr(shifts['v_difficulty'], shifts[mask])
    
    results.append({'mask':mask, 'r':r, 'p':p})
    
results = pandas.DataFrame(results)
_, p_fdr, _, _ = statsmodels.sandbox.stats.multicomp.multipletests(results['p'], method='fdr_bh')
results['p.fdr'] = p_fdr
results
%%R -i shifts
library(BayesFactor)

bf = regressionBF(v_difficulty ~ 1 + STh_L, data = shifts, whichModels='all')
print(bf)

bf = regressionBF(v_difficulty ~ 1 + STh_R, data = shifts, whichModels='all')
print(bf)<jupyter_output><empty_output><jupyter_text>Without errors<jupyter_code>import itertools
import pyfir
from nipy.modalities.fmri import design_matrix, experimental_paradigm
import statsmodels.api as sm

level2_dict = []
adj_r = []
# pars = ['neutral_cue', 'left_cue', 'right_cue', 'neutral_easy_correct', 'valid_easy_correct', 'invalid_easy_correct', 'neutral_hard_correct', 'valid_hard_correct', 'invalid_hard_correct', 'easy_error', 'hard_error',]
# pars = 

mask = 'STh_L'

for subject_id in np.arange(1, 20):
    for run in [1,2,3]:

        try:
            onsets_dict = get_onsets(subject_id, run, 'model1c', shift=-3)
            pars = onsets_dict.keys()

            conditions = []
            onsets = []
            for key in onsets_dict.keys():
                conditions += [key] * len(onsets_dict[key])
                onsets += onsets_dict[key]

            paradigm =  experimental_paradigm.BlockParadigm(con_id=conditions, 
                                                            onset=onsets,
                                                            duration=[[1.]] * len(conditions))
            tr = 3.0
            frametimes = np.linspace(0, 387 * tr, 387)
            hrf_model = 'Canonical with Derivative'

            X, names= design_matrix.dmtx_light(frametimes, paradigm, hrf_model=hrf_model, drift_model='polynomial', fir_delays=np.arange(11), drift_order=3)

            X /= X.max(0)
            X = pandas.DataFrame(X, columns=names)

            for i, mask in enumerate(['STh_L', 'STh_R']):
                signal = get_signal(subject_id,  run, mask)
                r = sm.GLSAR(signal, X).iterative_fit()

                adj_r.append({'subject_id':subject_id, 'mask':mask, 'r2_adj':r.rsquared_adj, 'r2':r.rsquared, 'run':run})

                for par in pars:
                    level2_dict.append({'value':r.params[par],
                                        'par':par,
                                        'subject_id':subject_id,
                                        'mask':mask,
                                        'run':run})
        except Exception as e:
            print 'Error: %s/%s/%s' %  (subject_id, mask, e)
df_glm = pandas.DataFrame(level2_dict)
df_glm = df_glm.query('par != "constant"')
df_glm['dummy'] = 1

df_glm['% signal change'] = df_glm['value']
sns.set_palette('Set1')
sns.set_context('poster')
sns.set_style('white')

# sns.factorplot(x='dummy', y='value', hue='par', data=df_glm, col='mask', col_wrap=3, kind='bar', ci=95, hue_order=['payoff_cue (left)', 'payoff_cue (right)', 'neutral_cue','difficult_rdm', 'easy_rdm'])
sns.factorplot(x='dummy', y='value', hue='par', data=df_glm, col='mask', col_wrap=3, kind='bar', ci=95)

df_glm_no_errors = df_glm.copy()
shifts = pandas.read_pickle('/home/gdholla1/projects/bias/data/derivatives/behavior/hddm/bias_shifts.pkl')
shifts = shifts.set_index('subj_idx')
tmp = df_glm.pivot_table(index='subject_id', values='value', columns=['par', 'mask'])

tmp = tmp['hard rdm (correct)'] - tmp['easy rdm (correct)']
shifts = shifts.merge(tmp, left_index=True, right_index=True)
from statsmodels.sandbox.stats import multicomp
import statsmodels
tmp = pandas.melt(shifts, id_vars=['v_cue_coding', 'v_difficulty', 'z_cue_coding'], var_name='mask')

tmp['mask'] = tmp['mask'].apply(lambda x: x.replace('STh', 'STN'))

# fac = sns.FacetGrid(tmp, col='mask', col_wrap=3)

# # def lmplot(x, y, color, **kwargs):
# #     sns.lmplot(x, y, **kwargs)

def make_r_title(x, y, **kwargs):
    r, p = sp.stats.pearsonr(x, y)
    
    title = plt.gca().get_title()
    
    ax = plt.gca()

    y = 0.85 - len(ax.texts) * 0.05
    ax.annotate('r = %.2f' % (r), xy=(0.05,y), xycoords='axes fraction')

tmp['v_difficulty'] *= -1

sns.set_palette(sns.color_palette('Set1')[1:])
    
fac = sns.lmplot('v_difficulty', 'value', data=tmp, col='mask', col_wrap=3)
fac.map(make_r_title, 'v_difficulty', 'value', )
# fac.set_titles('{col_name}')
fac.set_xlabels('Drift rate reduction')
fac.set_ylabels('hard RDM - easy RDM\n(% signal change)')
results = []

for mask in masks:
    r, p = sp.stats.pearsonr(shifts['v_difficulty'], shifts[mask])
    
    results.append({'mask':mask, 'r':r, 'p':p})
    
results = pandas.DataFrame(results)
_, p_fdr, _, _ = statsmodels.sandbox.stats.multicomp.multipletests(results['p'], method='fdr_bh')
results['p.fdr'] = p_fdr
results
%%R -i shifts
library(BayesFactor)

bf = regressionBF(v_difficulty ~ 1 + STh_L, data = shifts, whichModels='all')
print(bf)


bf = regressionBF(v_difficulty ~ 1 + STh_R, data = shifts, whichModels='all')
print(bf)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15861858190709047, 'fraction_numerical': 0.01115525672371638, 'mean_word_length': 3.010110294117647, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '272347', 'license': 'MIT', 'n_tokens_mistral': 4866, 'n_tokens_neox': 4575, 'n_words': 818, 'path': 'notebooks/model_based_analyses/Interindividual BOLD-HDDM correlations with whole masks.ipynb', 'repo_name': 'Gilles86/bias_project', 'size': '13088'}