08_GLM-ROI-SSTMSIT-split_striatum.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run GLMS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ['OPENBLAS_NUM_THREADS'] = '1'\n",
    "os.environ['MKL_NUM_THREADS'] = '1'\n",
    "os.environ['NUMEXPR_NUM_THREADS'] = '1'\n",
    "os.environ['OMP_NUM_THREADS'] = '1'\n",
    "\n",
    "import pandas as pd\n",
    "pd.set_option('display.max_rows', 500)\n",
    "pd.set_option('display.max_columns', 500)\n",
    "import numpy as np\n",
    "import shutil\n",
    "import itertools\n",
    "import statannot\n",
    "\n",
    "import glob\n",
    "import re\n",
    "import seaborn as sns\n",
    "import statsmodels.api as sm\n",
    "import statsmodels\n",
    "import scipy\n",
    "import nideconv\n",
    "from nideconv import GroupResponseFitter\n",
    "from scipy import signal\n",
    "\n",
    "import pickle as pkl\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.pyplot import figure\n",
    "import matplotlib.backends.backend_pdf\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "def flatten(t):\n",
    "    return [item for sublist in t for item in sublist]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def to_psc(x): # calculate percent signal change\n",
    "    return x / x.mean() * 100 - 100\n",
    "\n",
    "def load_events_confounds(sub, ses, task, run, include_physio=True, include_cosines=True):\n",
    "    event_fn = f'../derivatives/event_files/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_events.tsv'\n",
    "#    regressor_fn = f'../derivatives/behavior/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-model-regressors.tsv'\n",
    "    confounds_fn = f'../derivatives/fmriprep/fmriprep/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run}_desc-confounds_timeseries.tsv'\n",
    "    \n",
    "    events = pd.read_csv(event_fn, sep='\\t', index_col=None)\n",
    "    events['duration'] = .001\n",
    "            \n",
    "    # get confounds, cosines\n",
    "    confounds = pd.read_csv(confounds_fn, sep='\\t').fillna(method='bfill')\n",
    "    include_confs = ['trans_x', 'trans_y', 'trans_z', 'rot_x', 'rot_y', 'rot_z', 'dvars', 'framewise_displacement', 'global_signal']\n",
    "    if include_cosines:\n",
    "        include_confs += [x for x in confounds.columns if 'cos' in x]\n",
    "    confounds = confounds[include_confs]\n",
    "    \n",
    "    if include_physio:\n",
    "        run_idx = run#+1\n",
    "        retroicor_fn = f'../derivatives/retroicor/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_run-{run_idx}_desc-retroicor_regressors.tsv'\n",
    "        if not os.path.exists(retroicor_fn):\n",
    "            ## take first 20 aCompCor components\n",
    "            print(\"No retroicor found, including 20 a_comp_cor components\")\n",
    "            a_comp_cor = pd.read_csv(confounds_fn, sep='\\t')[['a_comp_cor_' + str(x).zfill(2) for x in range(20)]]\n",
    "#             confounds = pd.concat([confounds, a_comp_cor], axis=1)\n",
    "        else:\n",
    "            retroicor = pd.read_csv(retroicor_fn, sep='\\t', header=None).iloc[:,:20]  ## 20 components in total\n",
    "            retroicor.columns = ['cardiac_' + str(x) for x in range(6)] + ['respiratory_' + str(x) for x in range(8)] + ['respiratoryxcardiac_' + str(x) for x in range(4)] + ['HRV', 'RVT']\n",
    "            confounds = pd.concat([confounds, retroicor], axis=1)\n",
    "\n",
    "    return events, confounds\n",
    "\n",
    "def load_timeseries(atlas_type='MASSP', ses='sstmsit', task='msit'):\n",
    "\n",
    "    signal_fns = sorted(glob.glob(f'../derivatives/extracted_signals/sub-*/ses-{ses}/func/*task-{task}*{atlas_type}-signals.tsv'))\n",
    "\n",
    "    regex = re.compile(f'.*sub-(?P<sub>\\d+)_ses-{ses}_task-{task}_run-(?P<run>\\d)_desc-{atlas_type}-signals.tsv')\n",
    "    dfs = []\n",
    "    for signal_fn in signal_fns:\n",
    "        signals = pd.read_csv(signal_fn, sep='\\t')\n",
    "        gd = regex.match(signal_fn).groupdict()\n",
    "\n",
    "        if 'time' in signals.columns:\n",
    "            signals = signals.rename(columns={'time': 'volume'})\n",
    "\n",
    "\n",
    "        signals = signals.set_index(['volume']).apply(to_psc).reset_index()  # to PSC\n",
    "        signals['time'] = signals['volume'] * 1.38\n",
    "\n",
    "        del signals['volume']\n",
    "        signals['sub'] = gd['sub']\n",
    "        signals['run'] = int(gd['run'])\n",
    "        signals = signals.set_index(['sub', 'run', 'time'])\n",
    "        dfs.append(signals)\n",
    "\n",
    "    df = pd.concat(dfs)\n",
    "    \n",
    "    return df\n",
    "\n",
    "def sort_data(df, ses, task):\n",
    "    \n",
    "    all_events = []\n",
    "    all_confounds = []\n",
    "\n",
    "    for sub, run in df.reset_index().set_index(['sub', 'run']).index.unique():\n",
    "    #     if sub == '010' and run == '2':    # no RETROICOR\n",
    "    #         continue\n",
    "        if sub == '007' and run == '1':    # no RETROICOR\n",
    "            continue\n",
    "        events, confounds = load_events_confounds(sub, ses, task, run, include_physio=True)\n",
    "        events['sub'] = sub\n",
    "        events['run'] = run\n",
    "        confounds['sub'] = sub\n",
    "        confounds['run'] = run\n",
    "\n",
    "        all_events.append(events)\n",
    "        all_confounds.append(confounds)\n",
    "\n",
    "    events = pd.concat(all_events).set_index(['sub', 'run'])\n",
    "    confounds = pd.concat(all_confounds).set_index(['sub', 'run'])\n",
    "\n",
    "    events = events.rename(columns={'trial_type': 'event_type'})\n",
    "\n",
    "    events['duration'] = 0.001  # stick function\n",
    "    if task == 'msit': \n",
    "        events.loc[events['event_type']=='timing_feedback','duration'] = 0.4\n",
    "        \n",
    "    # make psc\n",
    "    confounds['global_signal'] = (confounds['global_signal']-confounds['global_signal'].mean())/confounds['global_signal'].std()\n",
    "\n",
    "    # change subject number ordering\n",
    "    subs = df.reset_index()['sub'].unique()\n",
    "    mapping = {y:x+1 for x,y in enumerate(subs)}\n",
    "\n",
    "    events = events.reset_index()\n",
    "    events['sub'] = events['sub'].replace(mapping).astype(str)\n",
    "    events = events.set_index(['sub', 'run'])\n",
    "\n",
    "    df = df.reset_index()\n",
    "    df['sub'] = df['sub'].replace(mapping).astype(str)\n",
    "    df = df.set_index(['sub', 'run', 'time'])\n",
    "\n",
    "    confounds = confounds.reset_index()\n",
    "    confounds['sub'] = confounds['sub'].replace(mapping).astype(str)\n",
    "    confounds = confounds.set_index(['sub', 'run'])\n",
    "\n",
    "    events.index = events.index.rename('subject', level=0)\n",
    "    df.index = df.index.rename('subject', level=0)\n",
    "    confounds.index = confounds.index.rename('subject', level=0)\n",
    "    \n",
    "    return events, df, confounds\n",
    "\n",
    "# # compares timecourse extraction for each extraction method. e.g. nilearn, manual, high passed or not.\n",
    "# def make_figs(subs, ses, task, params, show_plots=False, overwrite=False):\n",
    "\n",
    "#     if ses == 'rlsat':\n",
    "#         runs = [1,2,3]\n",
    "#     else:    \n",
    "#         runs = [1,2]\n",
    "\n",
    "#     # load timeseries all subs for extraction methods\n",
    "#     df_signals = []\n",
    "#     for p1, p2, p3, p4 in params:\n",
    "#         df = load_timeseries(atlas_type=p1, manual_extraction=p2, use_hp_data=p3, ses=ses, task=task)\n",
    "#         events, df, confounds = sort_data(df, ses, task)\n",
    "#         df['extraction_method'] = p4\n",
    "# #         df = df.sort_index() # speeds up index and stops 'impact performance' warning # but it fucks up the ordering\n",
    "#         df_signals.append(df)\n",
    "\n",
    "#     df_signals = pd.concat(df_signals).set_index(['extraction_method'],append=True)\n",
    "\n",
    "#     # loop through structures and make timeseries and correlation matrix\n",
    "#     for _sub in df_signals.index.unique(level='subject'):\n",
    "        \n",
    "#         sub = subs[int(_sub)-1]\n",
    "#         fn_ext_out = f'../derivatives/extracted_signals_plots/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_signals.pdf'\n",
    "#         fn_corr_out = f'../derivatives/extracted_signals_plots/sub-{sub}/ses-{ses}/func/sub-{sub}_ses-{ses}_task-{task}_corr.pdf'\n",
    "#         if os.path.isfile(fn_ext_out) and not overwrite:\n",
    "#             print(f'sub {sub} already ran, moving on..')\n",
    "#             continue\n",
    "#         os.makedirs(os.path.dirname(fn_ext_out), exist_ok=True)\n",
    "# #         _sub = sub.lstrip('0')\n",
    "#         # timeseries\n",
    "#         pdf = matplotlib.backends.backend_pdf.PdfPages(fn_ext_out)\n",
    "#         for roi in df_signals.columns:\n",
    "#             for run in runs:\n",
    "\n",
    "#                 roi_corr = df_signals.loc[_sub,run,:,'nilearn'].corrwith(df_signals.loc[_sub,run,:,'man_hp'])[roi]\n",
    "#                 fig = plt.figure(figsize=(40,10))\n",
    "#                 ax = plt.axes()\n",
    "#                 ax.set(xlabel='time (s)', ylabel='PSC', title = f'sub {sub} run {run} {roi} timecourse. Corr = {roi_corr}')\n",
    "#                 for ext in df_signals.index.unique(level='extraction_method'):\n",
    "\n",
    "#                     x = np.arange(0,len(df_signals.loc[_sub,run,:,ext][roi])*1.38,1.38)\n",
    "#                     plt.plot(x,df_signals.loc[_sub,run,:,ext][roi], label=ext)\n",
    "\n",
    "#                 plt.legend(loc='upper right')\n",
    "#                 pdf.savefig(fig)\n",
    "#                 if not show_plots:\n",
    "#                     plt.close()\n",
    "#         pdf.close()            \n",
    "#         # correlation matrix of nilearn signals vs manual extraction\n",
    "#         pdf = matplotlib.backends.backend_pdf.PdfPages(fn_corr_out)\n",
    "#         corr = pd.concat([df_signals.loc[_sub,:,:,'nilearn'],df_signals.loc[_sub,:,:,'man_hp']], axis=1).corr()\n",
    "#         fig = plt.figure(figsize=(30,23))\n",
    "#         sns.heatmap(corr, vmax=1,vmin=-1,\n",
    "#                 xticklabels=corr.columns,\n",
    "#                 yticklabels=corr.columns)\n",
    "#         pdf.savefig(fig)\n",
    "#         pdf.close()    \n",
    "#         if not show_plots:\n",
    "#             plt.close()\n",
    "\n",
    "#     return 0\n",
    "\n",
    "# function to plot the deconvolutions \n",
    "def plot_deconvs():\n",
    "    return"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## load timeseries\n",
    "\n",
    "separate motor and event glms due to differenes in trial types that have responses\n",
    "\n",
    "Regions of interest (SST):\n",
    "\n",
    "rIFG\n",
    "ACC\n",
    "M1\n",
    "SMA\n",
    "PaCG\n",
    "Ins\n",
    "GPe\n",
    "GPi\n",
    "RN\n",
    "SN\n",
    "STN\n",
    "STR\n",
    "Tha\n",
    "VTA\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# collect together all regions of interest for each sub.\n",
    "subs = [x.split('/')[-2].split('-')[-1] for x in sorted(glob.glob('../derivatives/extracted_signals/sub-*/ses-sstmsit'))]\n",
    "\n",
    "# remove subs that did not do task correctly\n",
    "# ['004', '010', '008', '019', '013', '027']\n",
    "\n",
    "to_remove = ['004', '010', '008', '019', '013', '027']\n",
    "subs = [x for x in subs if x not in to_remove]\n",
    "print(len(subs)) # should be 31\n",
    "\n",
    "# subs = ['002']\n",
    "runs = [1,2]\n",
    "tasks = ['sst']\n",
    "for sub in subs:\n",
    "#     print(sub)\n",
    "    for run in runs:\n",
    "#         print(run)\n",
    "        for task in tasks:\n",
    "            all_sigs = []#pd.DataFrame()\n",
    "            signals = sorted(glob.glob(f'../derivatives/extracted_signals/sub-{sub}/ses-sstmsit/func/*task-{task}_run-{run}*s.tsv'))\n",
    "#             print(signals)\n",
    "            if signals == []:\n",
    "                print(f'no data for sub {sub} run {run} task {task}')\n",
    "                continue\n",
    "            signals = [x for x in signals if 'ALL' not in x]\n",
    "            for index, sig in enumerate(signals):\n",
    "                if index == 0:\n",
    "                    all_sigs = pd.read_csv(sig, sep='\\t')\n",
    "                else:\n",
    "                    all_sigs = pd.concat([all_sigs, pd.read_csv(sig, sep='\\t')],axis=1)\n",
    "            \n",
    "            all_sigs = all_sigs.loc[:,~all_sigs.columns.duplicated()]\n",
    "            save_name = f'../derivatives/extracted_signals/sub-{sub}/ses-sstmsit/func/sub-{sub}_ses-sstmsit_task-{task}_run-{run}_desc-ALL-signals.tsv'\n",
    "            all_sigs.to_csv(save_name, '\\t')\n",
    "#             all_sigs.to_csv(sig.replace('THAL','ALL'),'\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def butter_highpass(cutoff, fs, order=5):\n",
    "    nyq = 0.5 * fs\n",
    "    normal_cutoff = cutoff / nyq\n",
    "    b, a = signal.butter(order, normal_cutoff, btype='high', analog=False)\n",
    "    return b, a\n",
    "\n",
    "def butter_highpass_filter(data, cutoff, fs, order=5):\n",
    "    b, a = butter_highpass(cutoff, fs, order=order)\n",
    "    y = signal.filtfilt(b, a, data)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Load timeseries\n",
    "atlas_type = 'ALL'\n",
    "signal_fns = sorted(glob.glob(f'../derivatives/extracted_signals/sub-*/ses-sstmsit/func/*task-sst*{atlas_type}*s_hp.tsv'))\n",
    "signal_fns = [x for x in signal_fns if not any(s in x for s in to_remove)] # remove subs that didnt do task right\n",
    "signal_fns\n",
    "\n",
    "filter_out_confounds = True\n",
    "filter_hp = False # already highpassed so no need\n",
    "\n",
    "# excluded_runs = [('009', '2'), #] #,  # fell asleep\n",
    "#                  ('010', '2'),  # physio failed   # TODO: find a better solution here!\n",
    "#                  ('010', '3'),  # physio failed   # TODO: find a better solution here!\n",
    "#                  ('018', '3')]  # physio failed   # TODO: find a better solution here!\n",
    "\n",
    "excluded_runs=[('018','2'),\n",
    "               ('025','2'),\n",
    "               ('041','1')]\n",
    "\n",
    "regex = re.compile(f'.*sub-(?P<sub>\\d+)_ses-sstmsit_task-sst_run-(?P<run>\\d)_desc-{atlas_type}-signals_hp.tsv')\n",
    "dfs = []\n",
    "for signal_fn in signal_fns:\n",
    "    signals = pd.read_csv(signal_fn, sep='\\t')\n",
    "    gd = regex.match(signal_fn).groupdict()\n",
    "    if tuple(gd.values()) in excluded_runs:\n",
    "        # run was excluded\n",
    "        continue\n",
    "    \n",
    "    if 'time' in signals.columns:\n",
    "        # if there's a column named 'time', it's called 'time' but it's really volume number..\n",
    "        signals = signals.rename(columns={'time': 'volume'})\n",
    "    signals = signals.set_index('volume')\n",
    "    \n",
    "    # psc?\n",
    "    signals = signals.apply(to_psc)  # to PSC\n",
    "\n",
    "    # filter out confounds?\n",
    "    if filter_out_confounds:\n",
    "        _, confounds = load_events_confounds(sub=gd['sub'], ses='sstmsit', task='sst', run=gd['run'])\n",
    "        confounds['intercept'] = 1   # add intercept!\n",
    "        betas, residuals, rank, s, = np.linalg.lstsq(a=confounds, b=signals)\n",
    "        signals_hat = confounds@betas\n",
    "        signals_hat.index = signals.index\n",
    "        signals_hat.columns = signals.columns\n",
    "        signals -= signals_hat   # residuals\n",
    "        \n",
    "    # high pass?\n",
    "    if filter_hp:\n",
    "        signals = signals.apply(lambda x: butter_highpass_filter(x, cutoff=1/128, fs=1/1.38) + x.mean(), axis=0)\n",
    "    \n",
    "    # index to time\n",
    "    signals.index *= 1.38\n",
    "    signals.index.name = 'time'\n",
    "    signals['subject'] = gd['sub']\n",
    "    signals['run'] = int(gd['run'])\n",
    "        \n",
    "    signals = signals.reset_index().set_index(['subject', 'run', 'time'])\n",
    "    dfs.append(signals)\n",
    "df = pd.concat(dfs)\n",
    "df = df.reindex(sorted(df.columns), axis=1)\n",
    "\n",
    "if atlas_type == 'ALL':\n",
    "    df.rename(columns = {'lM1':'M1-l', 'rM1':'M1-r', 'Caudate-l':'CN-l', 'Caudate-r':'CN-r', 'Putamen-l':'PUT-l', 'Putamen-r':'PUT-r'}, inplace = True)\n",
    "#     rois_ = sorted(['IFG-l','IFG-r','ACC-l','ACC-r','M1-l','M1-r','SMA-l','SMA-r','PaCG-l','PaCG-r','Ins-l','Ins-r','GPe-l','GPe-r','GPi-l','GPi-r','RN-l','RN-r','SN-l','SN-r',\n",
    "#              'STN-l','STN-r','Str-l','Str-r','Tha-l','Tha-r','VTA-l','VTA-r'])\n",
    "#     rois_ = ['ACC-l','ACC-r','IFG-l','IFG-r','Ins-l','Ins-r','M1-l','M1-r','pSG-l','pSG-r','SMA-l','SMA-r','SPL-l','SPL-r','GPe-l','GPe-r','GPi-l','GPi-r','RN-l','RN-r','SN-l','SN-r',\n",
    "#              'STN-l','STN-r','Str-l','Str-r','Tha-l','Tha-r','VTA-l','VTA-r']\n",
    "    rois_ = ['ACC-l','ACC-r','IFG-l','IFG-r','Ins-l','Ins-r','M1-l','M1-r','pSG-l','pSG-r','SMA-l','SMA-r','SPL-l','SPL-r','CN-l','CN-r','GPe-l','GPe-r','GPi-l','GPi-r','PUT-l','PUT-r','RN-l','RN-r','SN-l','SN-r',\n",
    "             'STN-l','STN-r','Tha-l','Tha-r','VTA-l','VTA-r']\n",
    "    df = df[rois_]\n",
    "    \n",
    "print(len(df.index.get_level_values(0).unique()))\n",
    "\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## load event, confounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_n = 'all_events'\n",
    "all_events = []\n",
    "all_confounds = []\n",
    "\n",
    "for sub, run in df.reset_index().set_index(['subject', 'run']).index.unique():\n",
    "    events, confounds = load_events_confounds(sub, 'sstmsit', 'sst', run, include_physio=True)\n",
    "    events['subject'] = sub\n",
    "    events['run'] = run\n",
    "    confounds['subject'] = sub\n",
    "    confounds['run'] = run\n",
    "    \n",
    "    all_events.append(events)\n",
    "    all_confounds.append(confounds)\n",
    "    \n",
    "events = pd.concat(all_events).set_index(['subject', 'run'])\n",
    "confounds = pd.concat(all_confounds).set_index(['subject', 'run'])\n",
    "\n",
    "events = events.rename(columns={'trial_type': 'event_type'})\n",
    "\n",
    "events['duration'] = 0.001  # stick function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## glm function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit_glm(timeseries, events, confounds, include_events, include_rois, t_r=1.38, oversample_design_matrix=10, concatenate_runs=False, fit_type='ols'):\n",
    "    events_1 = events.reset_index().set_index(['subject', 'run', 'event_type']).loc[(slice(None), slice(None), include_events),:]\n",
    "    events_1 = events_1.reset_index().set_index(['subject', 'run', 'event_type'])\n",
    "    events_1.onset -= t_r/2   # stc\n",
    "\n",
    "    glm1 = GroupResponseFitter(timeseries.copy()[include_rois],\n",
    "                               events_1, \n",
    "                               confounds=confounds.copy().reset_index() if confounds is not None else None,\n",
    "                               input_sample_rate=1/t_r, \n",
    "                               oversample_design_matrix=oversample_design_matrix,\n",
    "                               concatenate_runs=concatenate_runs)\n",
    "    for event_type in include_events:\n",
    "#         if event_type in ['stimulus_value_difference', 'feedback_PE']:\n",
    "#             glm1.add_event(event_type, covariates='modulation', basis_set='canonical_hrf', interval=[0, 19.32])\n",
    "#         else:\n",
    "        glm1.add_event(event_type, basis_set='canonical_hrf_with_time_derivative',  interval=[0, 19.32])\n",
    "\n",
    "    glm1.fit(type=fit_type)\n",
    "    return glm1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# fit event glm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if atlas_type == 'ATAG':\n",
    "    gm_nuclei = ['M1', 'PreSMA', 'IFG'] #'ACC',\n",
    "    include_rois = [hemi + roi for roi in gm_nuclei for hemi in ['l','r']]\n",
    "    include_rois.remove('lIFG')\n",
    "elif atlas_type=='MASSP':\n",
    "    gm_nuclei = ['Amg', 'Cl', 'GPe', 'GPi', 'PAG', 'PPN', 'RN', 'SN', 'STN', 'Str', 'Tha', 'VTA']\n",
    "    include_rois = [roi + '-' + hemi for roi in gm_nuclei for hemi in ['l', 'r']]\n",
    "elif atlas_type == 'CORT':\n",
    "    gm_nuclei = ['ACC','aSG','IFG','Ins','PaCG','pCC','postcG','precC','precGy','pSG','SMA','SPL']\n",
    "    include_rois = [roi + '-' + hemi for roi in gm_nuclei for hemi in ['l', 'r']]\n",
    "elif atlas_type == 'ALL':\n",
    "    #     gm_nuclei = sorted(['IFG','ACC','M1','SMA','PaCG','Ins','GPe','GPi','RN','SN','STN','Str','Tha','VTA'])\n",
    "    gm_nuclei = ['ACC','IFG','Ins','M1','pSG','SMA','SPL','CN','GPe','GPi','PUT','RN','SN','STN','Tha','VTA']\n",
    "    include_rois = [roi + '-' + hemi for roi in gm_nuclei for hemi in ['l', 'r']]\n",
    "\n",
    "# include_events=['response_left','response_right','fs','ss','go']\n",
    "include_events = ['fs','ss','go']\n",
    "glm1 = fit_glm(timeseries=df, events=events, confounds=confounds, include_events=include_events, include_rois=include_rois)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## fit motor glm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if atlas_type == 'ATAG':\n",
    "    gm_nuclei = ['M1', 'PreSMA', 'IFG'] #'ACC',\n",
    "    include_rois = [hemi + roi for roi in gm_nuclei for hemi in ['l','r']]\n",
    "    include_rois.remove('lIFG')\n",
    "elif atlas_type=='MASSP':\n",
    "    gm_nuclei = ['Amg', 'Cl', 'GPe', 'GPi', 'PAG', 'PPN', 'RN', 'SN', 'STN', 'Str', 'Tha', 'VTA']\n",
    "    include_rois = [roi + '-' + hemi for roi in gm_nuclei for hemi in ['l', 'r']]\n",
    "elif atlas_type == 'CORT':\n",
    "    gm_nuclei = ['ACC','aSG','IFG','Ins','PaCG','pCC','postcG','precC','precGy','pSG','SMA','SPL']\n",
    "    include_rois = [roi + '-' + hemi for roi in gm_nuclei for hemi in ['l', 'r']]\n",
    "elif atlas_type == 'ALL':\n",
    "#     gm_nuclei = sorted(['IFG','ACC','M1','SMA','PaCG','Ins','GPe','GPi','RN','SN','STN','Str','Tha','VTA'])\n",
    "    gm_nuclei = ['ACC','IFG','Ins','M1','pSG','SMA','SPL','CN','GPe','GPi','PUT','RN','SN','STN','Tha','VTA']\n",
    "    include_rois = [roi + '-' + hemi for roi in gm_nuclei for hemi in ['l', 'r']]\n",
    "\n",
    "# include_events=['response_left','response_right','fs','ss','go']\n",
    "include_events = ['response_left','response_right']\n",
    "glm_m = fit_glm(timeseries=df, events=events, confounds=confounds, include_events=include_events, include_rois=include_rois)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy\n",
    "\n",
    "glm1_ar1 = copy.deepcopy(glm1)\n",
    "\n",
    "tc = glm1.get_subjectwise_timecourses()\n",
    "tc = tc.reset_index().melt(id_vars=['subject', 'event type', 'covariate', 'time'], var_name='roi', value_name='signal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# motor\n",
    "all_betas = glm_m.response_fitters.groupby(['subject', 'run']).apply(lambda x: x.iloc[0].betas)                   # Extract all betas to dataframe\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['response_left', 'response_right'], 'intercept', 'HRF')] # Select only betas of interest\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True).groupby(level=0).diff()                   # Take difference right - left\n",
    "contrast_betas = contrast_betas.xs('response_right',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_right-left', ignore_index=False)\n",
    "\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "\n",
    "## recode to contralateral vs ipsilateral\n",
    "contrast_betas['beta_contra-ipsi'] = contrast_betas['beta_right-left'].copy()\n",
    "contrast_betas.loc[contrast_betas.hemisphere == 'r', 'beta_contra-ipsi'] *= -1\n",
    "contrast_betas_motor = contrast_betas.copy()\n",
    "\n",
    "## event glms\n",
    "all_betas = glm1.response_fitters.groupby(['subject', 'run']).apply(lambda x: x.iloc[0].betas)                   # Extract all betas to dataframe\n",
    "\n",
    "# fs - go\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['go', 'fs'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True).groupby(level=0).diff()                   # Take difference fs - go\n",
    "contrast_betas = contrast_betas.xs('fs',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_fs-go', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "contrast_betas_fsgo = contrast_betas.copy()\n",
    "\n",
    "# ss - fs\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['ss', 'fs'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "# subjectwise_betas[['STN-l', 'STN-r', 'Str-l', 'Str-r']].melt(ignore_index=False).reset_index().to_csv('./download-for-steven/betas_SAT.csv')\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True).groupby(level=0).diff()                   # Take difference fs - ss\n",
    "contrast_betas = contrast_betas.xs('fs',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_fs-ss', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "contrast_betas_fsss = contrast_betas.copy()\n",
    "\n",
    "# go - ss\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['go', 'ss'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "# subjectwise_betas[['STN-l', 'STN-r', 'Str-l', 'Str-r']].melt(ignore_index=False).reset_index().to_csv('./download-for-steven/betas_SAT.csv')\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True).groupby(level=0).diff()                   # Take difference ss - go\n",
    "contrast_betas = contrast_betas.xs('ss',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_ss-go', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "contrast_betas_ssgo = contrast_betas.copy()\n",
    "\n",
    "##########\n",
    "###########\n",
    "###########\n",
    "# now look just at absolute activtions with trial types\n",
    "############\n",
    "\n",
    "# fs\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['fs'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True)#.groupby(level=0).diff()                   # Take against baseline fs\n",
    "# contrast_betas = contrast_betas.xs('fs',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_fs', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "betas_fs = contrast_betas.copy()\n",
    "\n",
    "# ss\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['ss'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True)#.groupby(level=0).diff()                   # Take against baseline ss\n",
    "# contrast_betas = contrast_betas.xs('fs',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_ss', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "betas_ss = contrast_betas.copy()\n",
    "\n",
    "# go\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['go'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True)#.groupby(level=0).diff()                   # Take against baseline go\n",
    "# contrast_betas = contrast_betas.xs('fs',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_go', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "betas_go = contrast_betas.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# look at significance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import display_html\n",
    "from itertools import chain,cycle\n",
    "def display_side_by_side(*args,titles=cycle([''])):\n",
    "    html_str=''\n",
    "    for df,title in zip(args, chain(titles,cycle(['</br>'])) ):\n",
    "        html_str+='<th style=\"text-align:center\"><td style=\"vertical-align:top\">'\n",
    "        html_str+=f'<h2>{title}</h2>'\n",
    "        html_str+=df.to_html().replace('table','table style=\"display:inline\"')\n",
    "        html_str+='</td></th>'\n",
    "    display_html(html_str,raw=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# absolute activation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = betas_fs.groupby('ROI')['beta_fs'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p']).reindex(rois_)\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "# stat_df.loc[stat_df['p_corrected']<=0.0001,'significance']='****'\n",
    "stat_fs = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "tmp = betas_ss.groupby('ROI')['beta_ss'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p']).reindex(rois_)\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "# stat_df.loc[stat_df['p_corrected']<=0.0001,'significance']='****'\n",
    "stat_ss = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "tmp = betas_go.groupby('ROI')['beta_go'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p']).reindex(rois_)\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "# stat_df.loc[stat_df['p_corrected']<=0.0001,'significance']='****'\n",
    "stat_go = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "stat_fs, stat_ss, stat_go = stat_fs.reindex(rois_), stat_ss.reindex(rois_), stat_go.reindex(rois_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display_side_by_side(stat_fs,stat_ss,stat_go, titles=['fs','ss','go'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# contrast significance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = contrast_betas_fsgo.groupby('ROI')['beta_fs-go'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p']).reindex(rois_)\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "# stat_df.loc[stat_df['p_corrected']<=0.0001,'significance']='****'\n",
    "stat_fsgo = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "tmp = contrast_betas_fsss.groupby('ROI')['beta_fs-ss'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p']).reindex(rois_)\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "# stat_df.loc[stat_df['p_corrected']<=0.0001,'significance']='****'\n",
    "stat_fsss = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "tmp = contrast_betas_ssgo.groupby('ROI')['beta_ss-go'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p']).reindex(rois_)\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "# stat_df.loc[stat_df['p_corrected']<=0.0001,'significance']='****'\n",
    "stat_ssgo = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "tmp = contrast_betas_motor.groupby('ROI')['beta_contra-ipsi'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p']).reindex(rois_)\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "# stat_df.loc[stat_df['p_corrected']<=0.0001,'significance']='****'\n",
    "stat_leftright = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "stat_fsgo, stat_fsss, stat_ssgo, stat_leftright = stat_fsgo.reindex(rois_), stat_fsss.reindex(rois_), stat_ssgo.reindex(rois_), stat_leftright.reindex(rois_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display_side_by_side(stat_leftright,stat_fsss,stat_fsgo,stat_ssgo, titles=['contra-ipsilateral','fs > ss','fs > go','ss > go'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## plot all 4 contrasts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats\n",
    "\n",
    "def mean_confidence_interval(data, confidence=0.95):\n",
    "    a = 1.0 * np.array(data)\n",
    "    n = len(a)\n",
    "    m, se = np.mean(a), scipy.stats.sem(a)\n",
    "    h = se * scipy.stats.t.ppf((1 + confidence) / 2., n-1)\n",
    "    #return m, m-h, m+h\n",
    "    return m+h # only need highest \n",
    "\n",
    "def mean_confidence_interval_lower(data, confidence=0.95):\n",
    "    a = 1.0 * np.array(data)\n",
    "    n = len(a)\n",
    "    m, se = np.mean(a), scipy.stats.sem(a)\n",
    "    h = se * scipy.stats.t.ppf((1 + confidence) / 2., n-1)\n",
    "    #return m, m-h, m+h\n",
    "    return m-h # only need highest \n",
    "\n",
    "def add_asteriks(ax_n, data, col_name, significance=0):\n",
    "    l_count, r_count = 0, 0\n",
    "    \n",
    "    \n",
    "#     new = contrast_betas_inccon.groupby(['ROI_nohemi','hemisphere'])['beta_inc-con'].apply(mean_confidence_interval).reset_index().set_index('ROI_nohemi')\n",
    "#     new = new.loc[sorted(new.index,key=str.casefold)].drop_duplicates().reset_index().set_index(['ROI_nohemi','hemisphere'])\n",
    "    roi_data = data.groupby(['ROI_nohemi','hemisphere'])[col_name].apply(mean_confidence_interval).reset_index().set_index('ROI_nohemi')\n",
    "#     roi_data = roi_data.loc[sorted(roi_data.index,key=str.casefold)].drop_duplicates().reset_index().set_index(['ROI_nohemi','hemisphere'])\n",
    "    roi_data = roi_data.loc[roi_data.index].drop_duplicates().reset_index().set_index(['ROI_nohemi','hemisphere']).reindex(gm_nuclei,level=0)\n",
    "    \n",
    "    roi_data_lower = data.groupby(['ROI_nohemi','hemisphere'])[col_name].apply(mean_confidence_interval_lower).reset_index().set_index('ROI_nohemi')\n",
    "#     roi_data_lower = roi_data_lower.loc[sorted(roi_data_lower.index,key=str.casefold)].drop_duplicates().reset_index().set_index(['ROI_nohemi','hemisphere'])\n",
    "    roi_data_lower = roi_data_lower.loc[roi_data_lower.index].drop_duplicates().reset_index().set_index(['ROI_nohemi','hemisphere']).reindex(gm_nuclei,level=0)\n",
    "    \n",
    "#     print(roi_data_lower)\n",
    "#     for index, (y_val, region) in enumerate(zip(data.groupby(['ROI_nohemi','hemisphere'])[col_name].apply(mean_confidence_interval), data.ROI.unique())):\n",
    "    for index, (y_val, region, y_val_lower) in enumerate(zip(roi_data[col_name], data.ROI.unique(), roi_data_lower[col_name])):\n",
    "#         print(significance)\n",
    "#         print(region)\n",
    "        if index%2 == 0:\n",
    "            x_val = -0.26 + l_count\n",
    "            l_count+=1\n",
    "            if (aster:=len(significance[index])) == 2:\n",
    "                x_val = x_val - 0.06\n",
    "            elif aster == 3:\n",
    "                x_val = x_val - 0.12\n",
    "        else:\n",
    "            x_val = 0.14 + r_count\n",
    "            r_count+=1\n",
    "            if (aster:=len(significance[index])) == 2:\n",
    "                x_val = x_val - 0.06\n",
    "            elif aster == 3:\n",
    "                x_val = x_val - 0.12\n",
    "        if y_val >= 0:\n",
    "            ax[ax_n].text(x_val, y_val+0.005, significance[index])\n",
    "        else:\n",
    "            ax[ax_n].text(x_val, y_val_lower-0.017, significance[index])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "f, axes = plt.subplots(2,2, figsize=(30,17), sharex=False,sharey=True)\n",
    "ax = axes.ravel()\n",
    "f.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.1, hspace=0.27)\n",
    "# sns.set_style('darkgrid')\n",
    "\n",
    "sns.set_theme() # to make style changable from defaults use this line of code befor using set_style\n",
    "with sns.axes_style(\"darkgrid\"):\n",
    "\n",
    "#     # 1 ## MOTOR RESPONSE\n",
    "#     sns.barplot(x='ROI_nohemi', y='beta_contra-ipsi', hue='hemisphere', \n",
    "#                 data=contrast_betas_motor, ax=ax[0],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "#     ax[0].set_ylabel('% change', fontsize=22)\n",
    "#     ax[0].set_xlabel('')\n",
    "#     ax[0].legend_.remove()\n",
    "#     ax[0].set_title('Contralateral - ipsilateral', fontsize=26)\n",
    "#     ax[0].set(ylim=(-0.17,0.25))\n",
    "#     add_asteriks(0, contrast_betas_motor, 'beta_contra-ipsi', stat_leftright['significance'])\n",
    "#     # ax[0].set_xticklabels(ax[0].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "#     ax[0].set_yticklabels(ax[0].get_yticks().round(2),fontsize=20)\n",
    "#     ax[0].set_xticklabels(ax[0].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "\n",
    "    ## 2\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_fs-go', hue='hemisphere', \n",
    "                data=contrast_betas_fsgo, ax=ax[0],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    ax[0].set_ylabel('% change', fontsize=22)\n",
    "    ax[0].set_xlabel('')\n",
    "    ax[0].legend_.remove()\n",
    "    ax[0].set_title('FS - GO', fontsize=26)\n",
    "    ax[0].set(ylim=(-0.17,0.25))\n",
    "    add_asteriks(0, contrast_betas_fsgo, 'beta_fs-go', stat_fsgo['significance'])\n",
    "    ax[0].set_xticklabels(ax[0].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    # ax[1].set_xticklabels(ax[1].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    ax[0].set_yticklabels(ax[0].get_yticks().round(2),fontsize=20)\n",
    "\n",
    "    ## 3\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_fs-ss', hue='hemisphere', \n",
    "                data=contrast_betas_fsss, ax=ax[1],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    ax[1].set_ylabel('')\n",
    "    ax[1].set_xlabel('')\n",
    "    ax[1].legend_.remove()\n",
    "    ax[1].set_title('FS - SS', fontsize=26)\n",
    "    ax[1].set(ylim=(-0.17,0.25))\n",
    "    add_asteriks(1, contrast_betas_fsss, 'beta_fs-ss', stat_fsss['significance'])\n",
    "    ax[1].set_xticklabels(ax[1].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "#     ax[1].set_yticklabels(ax[1].get_yticks().round(2),fontsize=20)\n",
    "\n",
    "    ## 4\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_ss-go', hue='hemisphere', \n",
    "                data=contrast_betas_ssgo, ax=ax[2],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    ax[2].set_ylabel('% change', fontsize=22)\n",
    "    ax[2].set_xlabel('ROI', fontsize=22)\n",
    "#     ax[3].legend_.remove()\n",
    "    ax[2].set_title('SS - GO', fontsize=26)\n",
    "    ax[2].set(ylim=(-0.17,0.25))\n",
    "    add_asteriks(2, contrast_betas_ssgo, 'beta_ss-go', stat_ssgo['significance'])\n",
    "    ax[2].set_xticklabels(ax[2].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    ax[2].set_yticklabels(ax[2].get_yticks().round(2),fontsize=20)\n",
    "    plt.setp(ax[2].get_legend().get_texts(), fontsize='22') # for legend text\n",
    "    plt.setp(ax[2].get_legend().get_title(), fontsize='18') # for legend title\n",
    "\n",
    "    for x in range(4) :[t.set_color(i) for (i,t) in zip(['orange']*7+['red']*9,ax[x].get_xticklabels())]\n",
    "        \n",
    "    ax[3].axis('off')\n",
    "\n",
    "\n",
    "    f.savefig('figure_download_scott/GLM_ROI_SST_contrasts.pdf', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, axes = plt.subplots(2,2, figsize=(30,17), sharex=False,sharey=True)\n",
    "ax = axes.ravel()\n",
    "f.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.1, hspace=0.27)\n",
    "# sns.set_style('darkgrid')\n",
    "\n",
    "sns.set_theme() # to make style changable from defaults use this line of code befor using set_style\n",
    "with sns.axes_style(\"darkgrid\"):\n",
    "\n",
    "    ## 1 ## FS\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_fs', hue='hemisphere', \n",
    "                data=betas_fs, ax=ax[0],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    ax[0].set_ylabel('% change', fontsize=22)\n",
    "    ax[0].set_xlabel('')\n",
    "    ax[0].legend_.remove()\n",
    "    ax[0].set_title('FS', fontsize=26)\n",
    "    ax[0].set(ylim=(-0.15,0.35))\n",
    "    add_asteriks(0, betas_fs, 'beta_fs', stat_fs['significance'])\n",
    "    # ax[0].set_xticklabels(ax[0].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    ax[0].set_yticklabels(ax[0].get_yticks().round(2),fontsize=20)\n",
    "    ax[0].set_xticklabels(ax[0].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "\n",
    "    ## 2 SS\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_ss', hue='hemisphere', \n",
    "                data=betas_ss, ax=ax[1],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    ax[1].set_ylabel('')\n",
    "    ax[1].set_xlabel('')\n",
    "    ax[1].legend_.remove()\n",
    "    ax[1].set_title('SS', fontsize=26)\n",
    "    ax[1].set(ylim=(-0.15,0.35))\n",
    "    add_asteriks(1, betas_ss, 'beta_ss', stat_ss['significance'])\n",
    "    ax[1].set_xticklabels(ax[1].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "\n",
    "    ## 3 GO\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_go', hue='hemisphere', \n",
    "                data=betas_go, ax=ax[2],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    ax[2].set_ylabel('% change', fontsize=22)\n",
    "    ax[2].set_xlabel('', fontsize=22)\n",
    "#     ax[2].legend_.remove()\n",
    "    ax[2].set_title('GO', fontsize=26)\n",
    "    ax[2].set(ylim=(-0.15,0.35))\n",
    "    add_asteriks(2, betas_go, 'beta_go', stat_go['significance'])\n",
    "    # ax[2].set_xticklabels(ax[2].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    ax[2].set_yticklabels(ax[2].get_yticks().round(2),fontsize=20)\n",
    "    ax[2].set_xticklabels(ax[2].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    plt.setp(ax[2].get_legend().get_texts(), fontsize='22') # for legend text\n",
    "    plt.setp(ax[2].get_legend().get_title(), fontsize='18') # for legend title\n",
    "    \n",
    "    for x in range(4) :[t.set_color(i) for (i,t) in zip(['orange']*7+['red']*9,ax[x].get_xticklabels())]\n",
    "        \n",
    "    ax[3].set_visible(False)\n",
    "    \n",
    "    f.savefig('figure_download_scott/GLM_ROI_SST_supplementary.pdf', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load and save glm info for model based analysis\n",
    "# absolute values for FS, SS and GO trials for each region\n",
    "all_betas = glm1.response_fitters.groupby(['subject', 'run']).apply(lambda x: x.iloc[0].betas)                   # Extract all betas to dataframe\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "fs_betas = subjectwise_betas.loc[(slice(None),['fs'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "ss_betas = subjectwise_betas.loc[(slice(None),['ss'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "go_betas = subjectwise_betas.loc[(slice(None),['go'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "fs_betas.to_csv('scott/fs_betas_tsv.tsv',sep='\\t')\n",
    "ss_betas.to_csv('scott/ss_betas_tsv.tsv',sep='\\t')\n",
    "go_betas.to_csv('scott/go_betas_tsv.tsv',sep='\\t')\n",
    "# fs_reg = \n",
    "\n",
    "\n",
    "\n",
    "# # fs - go\n",
    "# subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "# subjectwise_betas = subjectwise_betas.loc[(slice(None),['fs'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "# contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True).groupby(level=0).diff()                   # Take difference SPD-ACC\n",
    "# contrast_betas = contrast_betas.xs('fs',level=1)\n",
    "# contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_fs-go', ignore_index=False)\n",
    "# if atlas_type == 'ATAG':\n",
    "#     contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "#     contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "# else:\n",
    "#     contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "#     contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "# contrast_betas_fsgo = contrast_betas.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MSIT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regions of interest (MSIT):\n",
    "    \n",
    "IFG\n",
    "ACC\n",
    "SMA\n",
    "SPL\n",
    "pSG\n",
    "Ins\n",
    "GPe\n",
    "GPi\n",
    "RN\n",
    "SN\n",
    "STN\n",
    "STR\n",
    "Tha\n",
    "VTA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# collect together all regions of interest for each sub.\n",
    "subs = [x.split('/')[-2].split('-')[-1] for x in sorted(glob.glob('../derivatives/extracted_signals/sub-*/ses-sstmsit'))]\n",
    "\n",
    "# remove subs that did not do task correctly\n",
    "# ['004', '010', '008', '019', '013', '027']\n",
    "\n",
    "to_remove = ['025']\n",
    "subs = [x for x in subs if x not in to_remove]\n",
    "print(len(subs)) # should be 36\n",
    "\n",
    "# subs = ['002']\n",
    "runs = [1,2]\n",
    "tasks = ['msit']\n",
    "for sub in subs:\n",
    "#     print(sub)\n",
    "    for run in runs:\n",
    "#         print(run)\n",
    "        for task in tasks:\n",
    "            all_sigs = []#pd.DataFrame()\n",
    "            signals = sorted(glob.glob(f'../derivatives/extracted_signals/sub-{sub}/ses-sstmsit/func/*task-{task}_run-{run}*hp.tsv'))\n",
    "#             print(signals)\n",
    "            if signals == []:\n",
    "                print(f'no data for sub {sub} run {run} task {task}')\n",
    "                continue\n",
    "            signals = [x for x in signals if 'ALL' not in x]\n",
    "            for index, sig in enumerate(signals):\n",
    "                if index == 0:\n",
    "                    all_sigs = pd.read_csv(sig, sep='\\t')\n",
    "                else:\n",
    "                    all_sigs = pd.concat([all_sigs, pd.read_csv(sig, sep='\\t')],axis=1)\n",
    "            \n",
    "            all_sigs = all_sigs.loc[:,~all_sigs.columns.duplicated()]\n",
    "            save_name = f'../derivatives/extracted_signals/sub-{sub}/ses-sstmsit/func/sub-{sub}_ses-sstmsit_task-{task}_run-{run}_desc-ALL-signals_hp.tsv'\n",
    "            all_sigs.to_csv(save_name, '\\t')\n",
    "#             all_sigs.to_csv(sig.replace('THAL','ALL'),'\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Load timeseries\n",
    "atlas_type = 'ALL'\n",
    "signal_fns = sorted(glob.glob(f'../derivatives/extracted_signals/sub-*/ses-sstmsit/func/*task-msit*{atlas_type}*s_hp.tsv'))\n",
    "signal_fns\n",
    "\n",
    "filter_out_confounds = True\n",
    "filter_hp = False\n",
    "\n",
    "excluded_runs=[('007','1'),\n",
    "               ('025','2'),\n",
    "               ('025','2')]\n",
    "\n",
    "regex = re.compile(f'.*sub-(?P<sub>\\d+)_ses-sstmsit_task-msit_run-(?P<run>\\d)_desc-{atlas_type}-signals_hp.tsv')\n",
    "dfs = []\n",
    "for signal_fn in signal_fns:\n",
    "    signals = pd.read_csv(signal_fn, sep='\\t')\n",
    "    gd = regex.match(signal_fn).groupdict()\n",
    "    if tuple(gd.values()) in excluded_runs:\n",
    "        # run was excluded\n",
    "        continue\n",
    "    \n",
    "    if 'time' in signals.columns:\n",
    "        # if there's a column named 'time', it's called 'time' but it's really volume number..\n",
    "        signals = signals.rename(columns={'time': 'volume'})\n",
    "    signals = signals.set_index('volume')\n",
    "    \n",
    "    # psc?\n",
    "    signals = signals.apply(to_psc)  # to PSC\n",
    "\n",
    "    # filter out confounds?\n",
    "    if filter_out_confounds:\n",
    "        _, confounds = load_events_confounds(sub=gd['sub'], ses='sstmsit', task='msit', run=gd['run'])\n",
    "        confounds['intercept'] = 1   # add intercept!\n",
    "        betas, residuals, rank, s, = np.linalg.lstsq(a=confounds, b=signals)\n",
    "        signals_hat = confounds@betas\n",
    "        signals_hat.index = signals.index\n",
    "        signals_hat.columns = signals.columns\n",
    "        signals -= signals_hat   # residuals\n",
    "        \n",
    "    # high pass?\n",
    "    if filter_hp:\n",
    "        signals = signals.apply(lambda x: butter_highpass_filter(x, cutoff=1/128, fs=1/1.38) + x.mean(), axis=0)\n",
    "    \n",
    "    \n",
    "    # index to time\n",
    "    signals.index *= 1.38\n",
    "    signals.index.name = 'time'\n",
    "    signals['subject'] = gd['sub']\n",
    "    signals['run'] = int(gd['run'])\n",
    "        \n",
    "    signals = signals.reset_index().set_index(['subject', 'run', 'time'])\n",
    "    dfs.append(signals)\n",
    "df = pd.concat(dfs)\n",
    "\n",
    "if atlas_type == 'ALL':\n",
    "    df.rename(columns = {'lM1':'M1-l', 'rM1':'M1-r', 'Caudate-l':'CN-l', 'Caudate-r':'CN-r', 'Putamen-l':'PUT-l', 'Putamen-r':'PUT-r'}, inplace = True)\n",
    "#     rois_ = sorted(['IFG-l','IFG-r','ACC-l','ACC-r','SPL-l','SPL-r','SMA-l','SMA-r','pSG-l','pSG-r','Ins-l','Ins-r','GPe-l','GPe-r','GPi-l','GPi-r','RN-l','RN-r','SN-l','SN-r',\n",
    "#              'STN-l','STN-r','Str-l','Str-r','Tha-l','Tha-r','VTA-l','VTA-r'])\n",
    "    rois_ = ['ACC-l','ACC-r','IFG-l','IFG-r','Ins-l','Ins-r','M1-l','M1-r','pSG-l','pSG-r','SMA-l','SMA-r','SPL-l','SPL-r','CN-l','CN-r','GPe-l','GPe-r','GPi-l','GPi-r','PUT-l','PUT-r','RN-l','RN-r','SN-l','SN-r',\n",
    "             'STN-l','STN-r','Tha-l','Tha-r','VTA-l','VTA-r']\n",
    "    df = df[rois_]\n",
    "\n",
    "print(len(df.index.get_level_values(0).unique()))    \n",
    "    \n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## load event, confounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_n = 'all_events'\n",
    "all_events = []\n",
    "all_confounds = []\n",
    "\n",
    "for sub, run in df.reset_index().set_index(['subject', 'run']).index.unique():\n",
    "    events, confounds = load_events_confounds(sub, 'sstmsit', 'msit', run, include_physio=True)\n",
    "    events['subject'] = sub\n",
    "    events['run'] = run\n",
    "    confounds['subject'] = sub\n",
    "    confounds['run'] = run\n",
    "    \n",
    "    all_events.append(events)\n",
    "    all_confounds.append(confounds)\n",
    "    \n",
    "events = pd.concat(all_events).set_index(['subject', 'run'])\n",
    "confounds = pd.concat(all_confounds).set_index(['subject', 'run'])\n",
    "\n",
    "events = events.rename(columns={'trial_type': 'event_type'})\n",
    "\n",
    "events['duration'] = 0.001  # stick function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## glm function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit_glm(timeseries, events, confounds, include_events, include_rois, t_r=1.38, oversample_design_matrix=10, concatenate_runs=False, fit_type='ols'):\n",
    "    events_1 = events.reset_index().set_index(['subject', 'run', 'event_type']).loc[(slice(None), slice(None), include_events),:]\n",
    "    events_1 = events_1.reset_index().set_index(['subject', 'run', 'event_type'])\n",
    "    events_1.onset -= t_r/2   # stc\n",
    "\n",
    "    glm1 = GroupResponseFitter(timeseries.copy()[include_rois],\n",
    "                               events_1, \n",
    "                               confounds=confounds.copy().reset_index() if confounds is not None else None,\n",
    "                               input_sample_rate=1/t_r, \n",
    "                               oversample_design_matrix=oversample_design_matrix,\n",
    "                               concatenate_runs=concatenate_runs)\n",
    "    for event_type in include_events:\n",
    "#         if event_type in ['stimulus_value_difference', 'feedback_PE']:\n",
    "#             glm1.add_event(event_type, covariates='modulation', basis_set='canonical_hrf', interval=[0, 19.32])\n",
    "#         else:\n",
    "        glm1.add_event(event_type, basis_set='canonical_hrf_with_time_derivative',  interval=[0, 19.32])\n",
    "\n",
    "    glm1.fit(type=fit_type)\n",
    "    return glm1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## fit events"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if atlas_type == 'ATAG':\n",
    "    gm_nuclei = ['M1', 'PreSMA', 'IFG'] #'ACC',\n",
    "    include_rois = [hemi + roi for roi in gm_nuclei for hemi in ['l','r']]\n",
    "    include_rois.remove('lIFG')\n",
    "elif atlas_type=='MASSP':\n",
    "    gm_nuclei = ['Amg', 'Cl', 'GPe', 'GPi', 'PAG', 'PPN', 'RN', 'SN', 'STN', 'Str', 'Tha', 'VTA']\n",
    "    include_rois = [roi + '-' + hemi for roi in gm_nuclei for hemi in ['l', 'r']]\n",
    "elif atlas_type == 'CORT':\n",
    "    gm_nuclei = sorted(['ACC','aSG','IFG','Ins','PaCG','pCC','postcG','precC','precGy','pSG','SMA','SPL'], key=str.casefold)\n",
    "    include_rois = [roi + '-' + hemi for roi in gm_nuclei for hemi in ['l', 'r']]\n",
    "elif atlas_type == 'ALL':\n",
    "#     gm_nuclei = sorted(['IFG','ACC','SPL','SMA','pSG','Ins','GPe','GPi','RN','SN','STN','Str','Tha','VTA'], key=str.casefold)\n",
    "    gm_nuclei = ['ACC','IFG','Ins','M1','pSG','SMA','SPL','CN','GPe','GPi','PUT','RN','SN','STN','Tha','VTA']\n",
    "    include_rois = [roi + '-' + hemi for roi in gm_nuclei for hemi in ['l', 'r']]\n",
    "\n",
    "# include_events=['response_left','response_right','fs','ss','go']\n",
    "# include_events = ['fs','ss','go']\n",
    "include_events = ['con','simon','flanker','inc']\n",
    "glm1 = fit_glm(timeseries=df, events=events, confounds=confounds, include_events=include_events, include_rois=include_rois)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## fit motor glm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if atlas_type == 'ATAG':\n",
    "    gm_nuclei = ['M1', 'PreSMA', 'IFG'] #'ACC',\n",
    "    include_rois = [hemi + roi for roi in gm_nuclei for hemi in ['l','r']]\n",
    "    include_rois.remove('lIFG')\n",
    "elif atlas_type=='MASSP':\n",
    "    gm_nuclei = ['Amg', 'Cl', 'GPe', 'GPi', 'PAG', 'PPN', 'RN', 'SN', 'STN', 'Str', 'Tha', 'VTA']\n",
    "    include_rois = [roi + '-' + hemi for roi in gm_nuclei for hemi in ['l', 'r']]\n",
    "elif atlas_type == 'CORT':\n",
    "    gm_nuclei = sorted(['ACC','aSG','IFG','Ins','PaCG','pCC','postcG','precC','precGy','pSG','SMA','SPL'], key=str.casefold)\n",
    "    include_rois = [roi + '-' + hemi for roi in gm_nuclei for hemi in ['l', 'r']]\n",
    "elif atlas_type == 'ALL':\n",
    "#     gm_nuclei = sorted(['IFG','ACC','SPL','SMA','pSG','Ins','GPe','GPi','RN','SN','STN','Str','Tha','VTA'], key=str.casefold)\n",
    "    gm_nuclei = ['ACC','IFG','Ins','M1','pSG','SMA','SPL','CN','GPe','GPi','PUT','RN','SN','STN','Tha','VTA']\n",
    "    include_rois = [roi + '-' + hemi for roi in gm_nuclei for hemi in ['l', 'r']]\n",
    "\n",
    "# include_events=['response_left','response_right','fs','ss','go']\n",
    "include_events = ['response_index','response_middle','response_ring']\n",
    "glm_m = fit_glm(timeseries=df, events=events, confounds=confounds, include_events=include_events, include_rois=include_rois)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy\n",
    "\n",
    "glm1_ar1 = copy.deepcopy(glm1)\n",
    "\n",
    "tc = glm1.get_subjectwise_timecourses()\n",
    "tc = tc.reset_index().melt(id_vars=['subject', 'event type', 'covariate', 'time'], var_name='roi', value_name='signal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# compare contrasts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the betas\n",
    "all_betas = glm1.response_fitters.groupby(['subject', 'run']).apply(lambda x: x.iloc[0].betas)                   # Extract all betas to dataframe\n",
    "\n",
    "# inc con\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['con', 'inc'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True).groupby(level=0).diff()                   # Take difference SPD-ACC\n",
    "contrast_betas = contrast_betas.xs('inc',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_inc-con', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "contrast_betas_inccon = contrast_betas.copy()\n",
    "\n",
    "# inc simon\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['simon', 'inc'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True).groupby(level=0).diff()                   # Take difference SPD-ACC\n",
    "contrast_betas = contrast_betas.xs('inc',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_inc-simon', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "contrast_betas_incsimon = contrast_betas.copy()\n",
    "\n",
    "# inc flanker\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['flanker', 'inc'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True).groupby(level=0).diff()                   # Take difference SPD-ACC\n",
    "contrast_betas = contrast_betas.xs('inc',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_inc-flanker', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "contrast_betas_incflanker = contrast_betas.copy()\n",
    "\n",
    "# simon flanker\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['flanker', 'simon'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True).groupby(level=0).diff()                   # Take difference SPD-ACC\n",
    "contrast_betas = contrast_betas.xs('simon',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_simon-flanker', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "contrast_betas_simonflanker = contrast_betas.copy()\n",
    "\n",
    "# simon con\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['con', 'simon'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True).groupby(level=0).diff()                   # Take difference SPD-ACC\n",
    "contrast_betas = contrast_betas.xs('simon',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_simon-con', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "contrast_betas_simoncon = contrast_betas.copy()\n",
    "\n",
    "# flanker con\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['con', 'flanker'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True).groupby(level=0).diff()                   # Take difference SPD-ACC\n",
    "contrast_betas = contrast_betas.xs('flanker',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_flanker-con', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "contrast_betas_flankercon = contrast_betas.copy()\n",
    "\n",
    "##########\n",
    "###########\n",
    "###########\n",
    "# now look just at absolute activtions with trial types\n",
    "############\n",
    "\n",
    "# Inc\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['inc'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True)#.groupby(level=0).diff()                   # Take difference SPD-ACC\n",
    "# contrast_betas = contrast_betas.xs('fs',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_inc', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "betas_inc = contrast_betas.copy()\n",
    "\n",
    "# sim\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['simon'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True)#.groupby(level=0).diff()                   # Take difference SPD-ACC\n",
    "# contrast_betas = contrast_betas.xs('fs',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_simon', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "betas_simon = contrast_betas.copy()\n",
    "\n",
    "# flanker\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['flanker'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True)#.groupby(level=0).diff()                   # Take difference SPD-ACC\n",
    "# contrast_betas = contrast_betas.xs('fs',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_flanker', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "betas_flanker = contrast_betas.copy()\n",
    "\n",
    "# con\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['con'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True)#.groupby(level=0).diff()                   # Take difference SPD-ACC\n",
    "# contrast_betas = contrast_betas.xs('fs',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_con', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "betas_con = contrast_betas.copy()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "betas_con"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# look at significance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = betas_inc.groupby('ROI')['beta_inc'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p'])\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "\n",
    "stat_inc = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "tmp = betas_flanker.groupby('ROI')['beta_flanker'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p'])\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "\n",
    "stat_flanker = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "tmp = betas_simon.groupby('ROI')['beta_simon'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p'])\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "\n",
    "stat_simon = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "tmp = betas_con.groupby('ROI')['beta_con'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p'])\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "\n",
    "stat_con = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "stat_inc, stat_flanker, stat_simon, stat_con = stat_inc.reindex(rois_), stat_flanker.reindex(rois_), stat_simon.reindex(rois_), stat_con.reindex(rois_)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display_side_by_side(stat_inc,stat_simon,stat_flanker,stat_con, titles=['inc','simon','flanker','con'])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = contrast_betas_inccon.groupby('ROI')['beta_inc-con'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p'])\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "# stat_df.loc[stat_df['p_corrected']<=0.0001,'significance']='****'\n",
    "\n",
    "stat_inccon = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "tmp = contrast_betas_incsimon.groupby('ROI')['beta_inc-simon'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p'])\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "# stat_df.loc[stat_df['p_corrected']<=0.0001,'significance']='****'\n",
    "\n",
    "stat_incsimon = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "tmp = contrast_betas_incflanker.groupby('ROI')['beta_inc-flanker'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p'])\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "# stat_df.loc[stat_df['p_corrected']<=0.0001,'significance']='****'\n",
    "\n",
    "stat_incflanker = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "tmp = contrast_betas_simonflanker.groupby('ROI')['beta_simon-flanker'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p'])\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "# stat_df.loc[stat_df['p_corrected']<=0.0001,'significance']='****'\n",
    "\n",
    "stat_simonflanker = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "tmp = contrast_betas_simoncon.groupby('ROI')['beta_simon-con'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p'])\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "# stat_df.loc[stat_df['p_corrected']<=0.0001,'significance']='****'\n",
    "\n",
    "stat_simoncon = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "tmp = contrast_betas_flankercon.groupby('ROI')['beta_flanker-con'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p'])\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "# stat_df.loc[stat_df['p_corrected']<=0.0001,'significance']='****'\n",
    "\n",
    "stat_flankercon = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "stat_inccon, stat_incsimon, stat_incflanker, stat_simonflanker, stat_simoncon, stat_flankercon = stat_inccon.reindex(rois_), stat_incsimon.reindex(rois_), stat_incflanker.reindex(rois_), stat_simonflanker.reindex(rois_), stat_simoncon.reindex(rois_), stat_flankercon.reindex(rois_)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display_side_by_side(stat_inccon,stat_incsimon,stat_incflanker,stat_simoncon,stat_simonflanker,stat_flankercon, titles=['inc > con','inc > simon','inc > flanker','simon > con', 'simon > flanker', 'flanker > con'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## plot all 6 contrasts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "f, axes = plt.subplots(3,2, figsize=(30,25), sharex=False,sharey=True)\n",
    "ax = axes.ravel()\n",
    "f.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.1, hspace=0.27)\n",
    "\n",
    "sns.set_theme() # to make style changable from defaults use this line of code befor using set_style\n",
    "with sns.axes_style(\"whitegrid\"):\n",
    "    ## 1\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_inc-con', hue='hemisphere', \n",
    "                data=contrast_betas_inccon, ax=ax[0],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    sns.set(font_scale=1)\n",
    "    ax[0].set_ylabel('% change', fontsize=22)\n",
    "    ax[0].set_xlabel('')\n",
    "    ax[0].legend_.remove()\n",
    "    ax[0].set_title('INC - CON', fontsize=26)\n",
    "    # ax[0].set(ylim=(-0.095,0.175))\n",
    "    ax[4].set(ylim=(-0.1,0.2))\n",
    "    add_asteriks(0, contrast_betas_inccon, 'beta_inc-con', stat_inccon['significance'])\n",
    "    # ax[0].set_xticklabels(ax[0].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    ax[0].set_yticklabels(ax[0].get_yticks().round(2),fontsize=20)\n",
    "    ax[0].set_xticklabels(ax[0].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    # plt.legend(loc='upper left')\n",
    "    # plt.setp(ax[4].get_legend().get_texts(), fontsize='17') # for legend text\n",
    "    # plt.setp(ax[4].get_legend().get_title(), fontsize='20') # for legend title\n",
    "\n",
    "    ## 2\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_inc-simon', hue='hemisphere', \n",
    "                data=contrast_betas_incsimon, ax=ax[1],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    ax[1].set_ylabel('')\n",
    "    ax[1].set_xlabel('')\n",
    "#     ax[1].legend_.remove()\n",
    "    ax[1].set_title('INC - SIM',fontsize=26)\n",
    "    ax[1].set(ylim=(-0.1,0.2))\n",
    "    add_asteriks(1, contrast_betas_incsimon, 'beta_inc-simon', stat_incsimon['significance'])\n",
    "    # ax[1].set_xticklabels(ax[1].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    # ax[1].set_yticklabels(ax[1].get_yticks().round(2),fontsize=20)\n",
    "    ax[1].set_xticklabels(ax[1].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    plt.setp(ax[1].get_legend().get_texts(), fontsize='22') # for legend text\n",
    "    plt.setp(ax[1].get_legend().get_title(), fontsize='18') # for legend title\n",
    "\n",
    "\n",
    "    ## 3\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_inc-flanker', hue='hemisphere', \n",
    "                data=contrast_betas_incflanker, ax=ax[2],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    ax[2].set_ylabel('% change',fontsize=22)\n",
    "    ax[2].set_xlabel('')\n",
    "    ax[2].set_title('INC - FLA',fontsize=26)\n",
    "    ax[2].legend_.remove()\n",
    "    ax[2].set(ylim=(-0.1,0.2))\n",
    "    add_asteriks(2, contrast_betas_incflanker, 'beta_inc-flanker', stat_incflanker['significance'])\n",
    "    # ax[0].set_xticklabels(ax[2].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    ax[2].set_yticklabels(ax[2].get_yticks().round(2),fontsize=20)\n",
    "    ax[2].set_xticklabels(ax[2].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "\n",
    "\n",
    "    ## 4\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_simon-con', hue='hemisphere', \n",
    "                data=contrast_betas_simoncon, ax=ax[3],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    ax[3].set_ylabel('')\n",
    "    ax[3].set_xlabel('')\n",
    "    ax[3].legend_.remove()\n",
    "    ax[3].set_title('SIM - CON', fontsize=26)\n",
    "    ax[3].set(ylim=(-0.1,0.2))\n",
    "    add_asteriks(3, contrast_betas_simoncon, 'beta_simon-con', stat_simoncon['significance'])\n",
    "    ax[3].set_xticklabels(ax[3].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "\n",
    "\n",
    "    ## 5\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_simon-flanker', hue='hemisphere', \n",
    "                data=contrast_betas_simonflanker, ax=ax[4],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    ax[4].set_ylabel('% change', fontsize=22)\n",
    "    ax[4].set_xlabel('ROI', fontsize=22)\n",
    "    ax[4].legend_.remove()\n",
    "    ax[4].set_title('SIM - FLA', fontsize=26)\n",
    "    ax[4].set(ylim=(-0.1,0.2))\n",
    "    add_asteriks(4, contrast_betas_simonflanker, 'beta_simon-flanker', stat_simonflanker['significance'])\n",
    "    ax[4].set_xticklabels(ax[4].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    ax[4].set_yticklabels(ax[4].get_yticks().round(2),fontsize=20)\n",
    "\n",
    "    ## 6\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_flanker-con', hue='hemisphere', \n",
    "                data=contrast_betas_flankercon, ax=ax[5],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    ax[5].set_ylabel('')\n",
    "    ax[5].set_xlabel('ROI', fontsize=22)\n",
    "    ax[5].legend_.remove()\n",
    "    ax[5].set_title('FLA - CON', fontsize=26)\n",
    "    ax[5].set(ylim=(-0.1,0.2))\n",
    "    add_asteriks(5, contrast_betas_flankercon, 'beta_flanker-con', stat_flankercon['significance'])\n",
    "    ax[5].set_xticklabels(ax[5].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    ax[5].set_yticklabels(ax[5].get_yticks().round(2),fontsize=20)\n",
    "\n",
    "    for x in range(6) :[t.set_color(i) for (i,t) in zip(['orange']*7+['red']*9,ax[x].get_xticklabels())]\n",
    "\n",
    "    f.savefig('figure_download_scott/GLM_ROI_MSIT_contrasts.pdf', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, axes = plt.subplots(2,2, figsize=(30,17), sharex=False,sharey=True)\n",
    "ax = axes.ravel()\n",
    "f.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.1, hspace=0.27)\n",
    "\n",
    "sns.set_theme() # to make style changable from defaults use this line of code befor using set_style\n",
    "with sns.axes_style(\"whitegrid\"):\n",
    "    ## 1\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_inc', hue='hemisphere', \n",
    "                data=betas_inc, ax=ax[0],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    sns.set(font_scale=1)\n",
    "    ax[0].set_ylabel('% change', fontsize=22)\n",
    "    ax[0].set_xlabel('')\n",
    "    ax[0].legend_.remove()\n",
    "    ax[0].set_title('INC', fontsize=26)\n",
    "    # ax[0].set(ylim=(-0.095,0.175))\n",
    "    ax[0].set(ylim=(-0.18,0.35))\n",
    "    add_asteriks(0, betas_inc, 'beta_inc', stat_inc['significance'])\n",
    "    # ax[0].set_xticklabels(ax[0].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    ax[0].set_yticklabels(ax[0].get_yticks().round(2),fontsize=20)\n",
    "    # plt.legend(loc='upper left')\n",
    "    # plt.setp(ax[0].get_legend().get_texts(), fontsize='17') # for legend text\n",
    "    # plt.setp(ax[0].get_legend().get_title(), fontsize='20') # for legend title\n",
    "    ax[0].set_xticklabels(ax[0].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "\n",
    "    ## 2\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_flanker', hue='hemisphere', \n",
    "                data=betas_flanker, ax=ax[1],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    ax[1].set_ylabel('')\n",
    "    ax[1].set_xlabel('')\n",
    "    ax[1].legend_.remove()\n",
    "    ax[1].set_title('FLA',fontsize=26)\n",
    "    ax[1].set(ylim=(-0.18,0.35))\n",
    "    add_asteriks(1, betas_flanker, 'beta_flanker', stat_flanker['significance'])\n",
    "    ax[1].set_xticklabels(ax[1].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "\n",
    "    ## 3\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_simon', hue='hemisphere', \n",
    "                data=betas_simon, ax=ax[2],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    ax[2].set_ylabel('% change',fontsize=22)\n",
    "    ax[2].set_xlabel('')\n",
    "    ax[2].set_title('SIM',fontsize=26)\n",
    "    ax[2].legend_.remove()\n",
    "    ax[2].set(ylim=(-0.18,0.35))\n",
    "    add_asteriks(2, betas_simon, 'beta_simon', stat_simon['significance'])\n",
    "    # ax[0].set_xticklabels(ax[2].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    ax[2].set_yticklabels(ax[2].get_yticks().round(2),fontsize=20)\n",
    "    ax[2].set_xticklabels(ax[2].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "\n",
    "    ## 4\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_con', hue='hemisphere', \n",
    "                data=betas_con, ax=ax[3],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    ax[3].set_ylabel('')\n",
    "    ax[3].set_xlabel('')\n",
    "#     ax[3].legend_.remove()\n",
    "    ax[3].set_title('CON', fontsize=26)\n",
    "    ax[3].set(ylim=(-0.18,0.35))\n",
    "    add_asteriks(3, betas_con, 'beta_con', stat_con['significance'])\n",
    "    ax[3].set_xticklabels(ax[3].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    plt.setp(ax[3].get_legend().get_texts(), fontsize='22') # for legend text\n",
    "    plt.setp(ax[3].get_legend().get_title(), fontsize='18') # for legend title\n",
    "\n",
    "    for x in range(4) :[t.set_color(i) for (i,t) in zip(['orange']*7+['red']*9,ax[x].get_xticklabels())]\n",
    "\n",
    "    f.savefig('figure_download_scott/GLM_ROI_MSIT_supplementary.pdf', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load and save glm info for model based analysis\n",
    "# absolute values for INC, SIMON and Flanker and INCONGRUENT trials for each region\n",
    "all_betas = glm1.response_fitters.groupby(['subject', 'run']).apply(lambda x: x.iloc[0].betas)                   # Extract all betas to dataframe\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "inc_betas = subjectwise_betas.loc[(slice(None),['inc'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "con_betas = subjectwise_betas.loc[(slice(None),['con'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "inc_betas.to_csv('scott/inc_betas_tsv.tsv', sep='\\t')\n",
    "con_betas.to_csv('scott/con_betas_tsv.tsv', sep='\\t')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# check IFG activation using HCP_MMP1 atlas instead of HarvardOxford (reviewers wants this for resubmission)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def butter_highpass(cutoff, fs, order=5):\n",
    "    nyq = 0.5 * fs\n",
    "    normal_cutoff = cutoff / nyq\n",
    "    b, a = signal.butter(order, normal_cutoff, btype='high', analog=False)\n",
    "    return b, a\n",
    "\n",
    "def butter_highpass_filter(data, cutoff, fs, order=5):\n",
    "    b, a = butter_highpass(cutoff, fs, order=order)\n",
    "    y = signal.filtfilt(b, a, data)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Load timeseries\n",
    "atlas_type = 'ALL'\n",
    "signal_fns = sorted(glob.glob(f'../derivatives/extracted_signals/sub-*/ses-sstmsit/func/*task-sst*{atlas_type}*hp.tsv'))\n",
    "signal_fns = [x for x in signal_fns if not any(s in x for s in to_remove)] # remove subs that didnt do task right\n",
    "signal_fns\n",
    "\n",
    "filter_out_confounds = True\n",
    "filter_hp = False # already highpassed so no need\n",
    "\n",
    "# excluded_runs = [('009', '2'), #] #,  # fell asleep\n",
    "#                  ('010', '2'),  # physio failed   # TODO: find a better solution here!\n",
    "#                  ('010', '3'),  # physio failed   # TODO: find a better solution here!\n",
    "#                  ('018', '3')]  # physio failed   # TODO: find a better solution here!\n",
    "\n",
    "excluded_runs=[('018','2'),\n",
    "               ('025','2'),\n",
    "               ('041','1')]\n",
    "\n",
    "regex = re.compile(f'.*sub-(?P<sub>\\d+)_ses-sstmsit_task-sst_run-(?P<run>\\d)_desc-{atlas_type}-signals_hp.tsv')\n",
    "dfs = []\n",
    "for signal_fn in signal_fns:\n",
    "    signals = pd.read_csv(signal_fn, sep='\\t')\n",
    "    gd = regex.match(signal_fn).groupdict()\n",
    "    if tuple(gd.values()) in excluded_runs:\n",
    "        # run was excluded\n",
    "        continue\n",
    "    \n",
    "    if 'time' in signals.columns:\n",
    "        # if there's a column named 'time', it's called 'time' but it's really volume number..\n",
    "        signals = signals.rename(columns={'time': 'volume'})\n",
    "    signals = signals.set_index('volume')\n",
    "    \n",
    "    # psc?\n",
    "    signals = signals.apply(to_psc)  # to PSC\n",
    "\n",
    "    # filter out confounds?\n",
    "    if filter_out_confounds:\n",
    "        _, confounds = load_events_confounds(sub=gd['sub'], ses='sstmsit', task='sst', run=gd['run'])\n",
    "        confounds['intercept'] = 1   # add intercept!\n",
    "        betas, residuals, rank, s, = np.linalg.lstsq(a=confounds, b=signals)\n",
    "        signals_hat = confounds@betas\n",
    "        signals_hat.index = signals.index\n",
    "        signals_hat.columns = signals.columns\n",
    "        signals -= signals_hat   # residuals\n",
    "        \n",
    "    # high pass?\n",
    "    if filter_hp:\n",
    "        signals = signals.apply(lambda x: butter_highpass_filter(x, cutoff=1/128, fs=1/1.38) + x.mean(), axis=0)\n",
    "    \n",
    "    # index to time\n",
    "    signals.index *= 1.38\n",
    "    signals.index.name = 'time'\n",
    "    signals['subject'] = gd['sub']\n",
    "    signals['run'] = int(gd['run'])\n",
    "        \n",
    "    signals = signals.reset_index().set_index(['subject', 'run', 'time'])\n",
    "    dfs.append(signals)\n",
    "df = pd.concat(dfs)\n",
    "df = df.reindex(sorted(df.columns), axis=1)\n",
    "\n",
    "if atlas_type == 'ALL':\n",
    "    df.rename(columns = {'lM1':'M1-l', 'rM1':'M1-r'}, inplace = True)\n",
    "#     rois_ = sorted(['IFG-l','IFG-r','ACC-l','ACC-r','M1-l','M1-r','SMA-l','SMA-r','PaCG-l','PaCG-r','Ins-l','Ins-r','GPe-l','GPe-r','GPi-l','GPi-r','RN-l','RN-r','SN-l','SN-r',\n",
    "#              'STN-l','STN-r','Str-l','Str-r','Tha-l','Tha-r','VTA-l','VTA-r'])\n",
    "    rois_ = ['ACC-l','ACC-r','IFGhcp-l','IFGhcp-r','Ins-l','Ins-r','M1-l','M1-r','pSG-l','pSG-r','SMA-l','SMA-r','SPL-l','SPL-r','GPe-l','GPe-r','GPi-l','GPi-r','RN-l','RN-r','SN-l','SN-r',\n",
    "             'STN-l','STN-r','Str-l','Str-r','Tha-l','Tha-r','VTA-l','VTA-r']\n",
    "    df = df[rois_]\n",
    "\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## load event, confounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_n = 'all_events'\n",
    "all_events = []\n",
    "all_confounds = []\n",
    "\n",
    "for sub, run in df.reset_index().set_index(['subject', 'run']).index.unique():\n",
    "    events, confounds = load_events_confounds(sub, 'sstmsit', 'sst', run, include_physio=True)\n",
    "    events['subject'] = sub\n",
    "    events['run'] = run\n",
    "    confounds['subject'] = sub\n",
    "    confounds['run'] = run\n",
    "    \n",
    "    all_events.append(events)\n",
    "    all_confounds.append(confounds)\n",
    "    \n",
    "events = pd.concat(all_events).set_index(['subject', 'run'])\n",
    "confounds = pd.concat(all_confounds).set_index(['subject', 'run'])\n",
    "\n",
    "events = events.rename(columns={'trial_type': 'event_type'})\n",
    "\n",
    "events['duration'] = 0.001  # stick function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## glm function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit_glm(timeseries, events, confounds, include_events, include_rois, t_r=1.38, oversample_design_matrix=10, concatenate_runs=False, fit_type='ols'):\n",
    "    events_1 = events.reset_index().set_index(['subject', 'run', 'event_type']).loc[(slice(None), slice(None), include_events),:]\n",
    "    events_1 = events_1.reset_index().set_index(['subject', 'run', 'event_type'])\n",
    "    events_1.onset -= t_r/2   # stc\n",
    "\n",
    "    glm1 = GroupResponseFitter(timeseries.copy()[include_rois],\n",
    "                               events_1, \n",
    "                               confounds=confounds.copy().reset_index() if confounds is not None else None,\n",
    "                               input_sample_rate=1/t_r, \n",
    "                               oversample_design_matrix=oversample_design_matrix,\n",
    "                               concatenate_runs=concatenate_runs)\n",
    "    for event_type in include_events:\n",
    "#         if event_type in ['stimulus_value_difference', 'feedback_PE']:\n",
    "#             glm1.add_event(event_type, covariates='modulation', basis_set='canonical_hrf', interval=[0, 19.32])\n",
    "#         else:\n",
    "        glm1.add_event(event_type, basis_set='canonical_hrf_with_time_derivative',  interval=[0, 19.32])\n",
    "\n",
    "    glm1.fit(type=fit_type)\n",
    "    return glm1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# fit event glm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if atlas_type == 'ATAG':\n",
    "    gm_nuclei = ['M1', 'PreSMA', 'IFG'] #'ACC',\n",
    "    include_rois = [hemi + roi for roi in gm_nuclei for hemi in ['l','r']]\n",
    "    include_rois.remove('lIFG')\n",
    "elif atlas_type=='MASSP':\n",
    "    gm_nuclei = ['Amg', 'Cl', 'GPe', 'GPi', 'PAG', 'PPN', 'RN', 'SN', 'STN', 'Str', 'Tha', 'VTA']\n",
    "    include_rois = [roi + '-' + hemi for roi in gm_nuclei for hemi in ['l', 'r']]\n",
    "elif atlas_type == 'CORT':\n",
    "    gm_nuclei = ['ACC','aSG','IFG','Ins','PaCG','pCC','postcG','precC','precGy','pSG','SMA','SPL']\n",
    "    include_rois = [roi + '-' + hemi for roi in gm_nuclei for hemi in ['l', 'r']]\n",
    "elif atlas_type == 'ALL':\n",
    "    #     gm_nuclei = sorted(['IFG','ACC','M1','SMA','PaCG','Ins','GPe','GPi','RN','SN','STN','Str','Tha','VTA'])\n",
    "    gm_nuclei = ['ACC','IFGhcp','Ins','M1','pSG','SMA','SPL','GPe','GPi','RN','SN','STN','Str','Tha','VTA']\n",
    "    include_rois = [roi + '-' + hemi for roi in gm_nuclei for hemi in ['l', 'r']]\n",
    "\n",
    "# include_events=['response_left','response_right','fs','ss','go']\n",
    "include_events = ['fs','ss','go']\n",
    "glm1 = fit_glm(timeseries=df, events=events, confounds=confounds, include_events=include_events, include_rois=include_rois)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## fit motor glm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if atlas_type == 'ATAG':\n",
    "    gm_nuclei = ['M1', 'PreSMA', 'IFG'] #'ACC',\n",
    "    include_rois = [hemi + roi for roi in gm_nuclei for hemi in ['l','r']]\n",
    "    include_rois.remove('lIFG')\n",
    "elif atlas_type=='MASSP':\n",
    "    gm_nuclei = ['Amg', 'Cl', 'GPe', 'GPi', 'PAG', 'PPN', 'RN', 'SN', 'STN', 'Str', 'Tha', 'VTA']\n",
    "    include_rois = [roi + '-' + hemi for roi in gm_nuclei for hemi in ['l', 'r']]\n",
    "elif atlas_type == 'CORT':\n",
    "    gm_nuclei = ['ACC','aSG','IFG','Ins','PaCG','pCC','postcG','precC','precGy','pSG','SMA','SPL']\n",
    "    include_rois = [roi + '-' + hemi for roi in gm_nuclei for hemi in ['l', 'r']]\n",
    "elif atlas_type == 'ALL':\n",
    "#     gm_nuclei = sorted(['IFG','ACC','M1','SMA','PaCG','Ins','GPe','GPi','RN','SN','STN','Str','Tha','VTA'])\n",
    "    gm_nuclei = ['ACC','IFGhcp','Ins','M1','pSG','SMA','SPL','GPe','GPi','RN','SN','STN','Str','Tha','VTA']\n",
    "    include_rois = [roi + '-' + hemi for roi in gm_nuclei for hemi in ['l', 'r']]\n",
    "\n",
    "# include_events=['response_left','response_right','fs','ss','go']\n",
    "include_events = ['response_left','response_right']\n",
    "glm_m = fit_glm(timeseries=df, events=events, confounds=confounds, include_events=include_events, include_rois=include_rois)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy\n",
    "\n",
    "glm1_ar1 = copy.deepcopy(glm1)\n",
    "\n",
    "tc = glm1.get_subjectwise_timecourses()\n",
    "tc = tc.reset_index().melt(id_vars=['subject', 'event type', 'covariate', 'time'], var_name='roi', value_name='signal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_betas = glm1.response_fitters.groupby(['subject', 'run']).apply(lambda x: x.iloc[0].betas)                   # Extract all betas to dataframe\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['go', 'fs', 'ss'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "# contrast_betas = subjectwise_betas\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='signal_change', ignore_index=False)\n",
    "contrast_betas = contrast_betas.reset_index(level=1)\n",
    "\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "    \n",
    "f, axes = plt.subplots(3,5,figsize=(25,10))\n",
    "\n",
    "# for i, ax in zip(range(10), axes.flat):\n",
    "#     sns.barplot(data[i,0,0,0], hist=False, ax=ax)\n",
    "# plt.show()\n",
    "\n",
    "ax = axes.ravel()\n",
    "\n",
    "\n",
    "for i, r in enumerate(gm_nuclei):\n",
    "    sns.barplot(x='ROI_nohemi', y='signal_change', hue='event type', \n",
    "                data=contrast_betas.loc[(contrast_betas.hemisphere=='l') & (contrast_betas.ROI_nohemi==r)], ax=ax[i])\n",
    "    ax[i].legend_.remove()\n",
    "\n",
    "# box_pairs = []\n",
    "# for index in contrast_betas['ROI_nohemi'].unique():\n",
    "#     box_pairs = box_pairs + (list(itertools.combinations(list(itertools.combinations([index,'go', 'fs', 'ss'],2))[:3],2)))\n",
    "    \n",
    "# # box_pairs = box_pairs[:1]\n",
    "# statannot.add_stat_annotation(\n",
    "#     ax,\n",
    "#     data=contrast_betas, \n",
    "#     x='ROI_nohemi',\n",
    "#     y='signal_change',\n",
    "#     hue='event type',\n",
    "# #     hue_order=['go','fs','ss'],\n",
    "#     box_pairs=box_pairs,\n",
    "#     test=\"t-test_paired\", #options=['t-test_ind', 't-test_welch', 't-test_paired','Mann-Whitney', 'Mann-Whitney-gt', 'Mann-Whitney-ls','Levene', 'Wilcoxon', 'Kruskal']\n",
    "#     text_format=\"star\",\n",
    "#     loc=\"inside\",\n",
    "# #     text_offset=-0.5,\n",
    "# #     line_offset_to_box=-0.5,\n",
    "# #     line_offset=-0.5,\n",
    "# #     use_fixed_offset=True\n",
    "# #     line_height=1\n",
    "# )\n",
    "\n",
    "# ax.set_ylabel('% change')\n",
    "# ax.set_xlabel('ROI')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_betas = glm1.response_fitters.groupby(['subject', 'run']).apply(lambda x: x.iloc[0].betas)                   # Extract all betas to dataframe\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['go', 'fs', 'ss'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "# contrast_betas = subjectwise_betas\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='signal_change', ignore_index=False)\n",
    "contrast_betas = contrast_betas.reset_index(level=1)\n",
    "\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "    \n",
    "f, ax = plt.subplots(1,1,figsize=(25,10))\n",
    "sns.barplot(x='ROI_nohemi', y='signal_change', hue='event type', \n",
    "            data=contrast_betas, ax=ax)\n",
    "\n",
    "box_pairs = []\n",
    "for index in contrast_betas['ROI_nohemi'].unique():\n",
    "    box_pairs = box_pairs + (list(itertools.combinations(list(itertools.combinations([index,'go', 'fs', 'ss'],2))[:3],2)))\n",
    "    \n",
    "# box_pairs = box_pairs[:1]\n",
    "statannot.add_stat_annotation(\n",
    "    ax,\n",
    "    data=contrast_betas, \n",
    "    x='ROI_nohemi',\n",
    "    y='signal_change',\n",
    "    hue='event type',\n",
    "#     hue_order=['go','fs','ss'],\n",
    "    box_pairs=box_pairs,\n",
    "    test=\"t-test_paired\", #options=['t-test_ind', 't-test_welch', 't-test_paired','Mann-Whitney', 'Mann-Whitney-gt', 'Mann-Whitney-ls','Levene', 'Wilcoxon', 'Kruskal']\n",
    "    text_format=\"star\",\n",
    "    loc=\"inside\",\n",
    "#     text_offset=-0.5,\n",
    "#     line_offset_to_box=-0.5,\n",
    "#     line_offset=-0.5,\n",
    "#     use_fixed_offset=True\n",
    "#     line_height=1\n",
    ")\n",
    "\n",
    "ax.set_ylabel('% change')\n",
    "ax.set_xlabel('ROI')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_betas = glm_m.response_fitters.groupby(['subject', 'run']).apply(lambda x: x.iloc[0].betas)                   # Extract all betas to dataframe\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['response_left','response_right'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "# contrast_betas = subjectwise_betas\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='signal_change', ignore_index=False)\n",
    "contrast_betas = contrast_betas.reset_index(level=1)\n",
    "\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "    \n",
    "f, ax = plt.subplots(1,1,figsize=(25,10))\n",
    "sns.barplot(x='ROI', y='signal_change', hue='event type', \n",
    "            data=contrast_betas, ax=ax)\n",
    "\n",
    "box_pairs = []\n",
    "for index in contrast_betas['ROI'].unique():\n",
    "    box_pairs = box_pairs + (list(itertools.combinations(list(itertools.combinations([index,'response_left','response_right'],2))[:2],2)))\n",
    "    \n",
    "# box_pairs = box_pairs[:1]\n",
    "# print(box_pairs)\n",
    "statannot.add_stat_annotation(\n",
    "    ax,\n",
    "    data=contrast_betas, \n",
    "    x='ROI',\n",
    "    y='signal_change',\n",
    "    hue='event type',\n",
    "#     hue_order=['go','fs','ss'],\n",
    "    box_pairs=box_pairs,\n",
    "    test=\"t-test_paired\", #options=['t-test_ind', 't-test_welch', 't-test_paired','Mann-Whitney', 'Mann-Whitney-gt', 'Mann-Whitney-ls','Levene', 'Wilcoxon', 'Kruskal']\n",
    "    text_format=\"star\",\n",
    "    loc=\"inside\",\n",
    "#     text_offset=-0.5,\n",
    "#     line_offset_to_box=-0.5,\n",
    "#     line_offset=-0.5,\n",
    "#     use_fixed_offset=True\n",
    "#     line_height=1\n",
    ")\n",
    "\n",
    "ax.set_ylabel('% change')\n",
    "ax.set_xlabel('ROI')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# motor\n",
    "all_betas = glm_m.response_fitters.groupby(['subject', 'run']).apply(lambda x: x.iloc[0].betas)                   # Extract all betas to dataframe\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['response_left', 'response_right'], 'intercept', 'HRF')] # Select only betas of interest\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True).groupby(level=0).diff()                   # Take difference right - left\n",
    "contrast_betas = contrast_betas.xs('response_right',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_right-left', ignore_index=False)\n",
    "\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "\n",
    "## recode to contralateral vs ipsilateral\n",
    "contrast_betas['beta_contra-ipsi'] = contrast_betas['beta_right-left'].copy()\n",
    "contrast_betas.loc[contrast_betas.hemisphere == 'r', 'beta_contra-ipsi'] *= -1\n",
    "contrast_betas_motor = contrast_betas.copy()\n",
    "\n",
    "## event glms\n",
    "all_betas = glm1.response_fitters.groupby(['subject', 'run']).apply(lambda x: x.iloc[0].betas)                   # Extract all betas to dataframe\n",
    "\n",
    "# fs - go\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['go', 'fs'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True).groupby(level=0).diff()                   # Take difference fs - go\n",
    "contrast_betas = contrast_betas.xs('fs',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_fs-go', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "contrast_betas_fsgo = contrast_betas.copy()\n",
    "\n",
    "# ss - fs\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['ss', 'fs'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "# subjectwise_betas[['STN-l', 'STN-r', 'Str-l', 'Str-r']].melt(ignore_index=False).reset_index().to_csv('./download-for-steven/betas_SAT.csv')\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True).groupby(level=0).diff()                   # Take difference fs - ss\n",
    "contrast_betas = contrast_betas.xs('fs',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_fs-ss', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "contrast_betas_fsss = contrast_betas.copy()\n",
    "\n",
    "# go - ss\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['go', 'ss'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "# subjectwise_betas[['STN-l', 'STN-r', 'Str-l', 'Str-r']].melt(ignore_index=False).reset_index().to_csv('./download-for-steven/betas_SAT.csv')\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True).groupby(level=0).diff()                   # Take difference ss - go\n",
    "contrast_betas = contrast_betas.xs('ss',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_ss-go', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "contrast_betas_ssgo = contrast_betas.copy()\n",
    "\n",
    "##########\n",
    "###########\n",
    "###########\n",
    "# now look just at absolute activtions with trial types\n",
    "############\n",
    "\n",
    "# fs\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['fs'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True)#.groupby(level=0).diff()                   # Take against baseline fs\n",
    "# contrast_betas = contrast_betas.xs('fs',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_fs', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "betas_fs = contrast_betas.copy()\n",
    "\n",
    "# ss\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['ss'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True)#.groupby(level=0).diff()                   # Take against baseline ss\n",
    "# contrast_betas = contrast_betas.xs('fs',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_ss', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "betas_ss = contrast_betas.copy()\n",
    "\n",
    "# go\n",
    "subjectwise_betas = all_betas.groupby(level=[0,2,3,4]).mean()                                                    # Mean across runs (fixed-effects)\n",
    "subjectwise_betas = subjectwise_betas.loc[(slice(None),['go'], 'intercept', 'HRF')]              # Select only betas of interest\n",
    "\n",
    "contrast_betas = subjectwise_betas.reset_index(level=[2,3], drop=True)#.groupby(level=0).diff()                   # Take against baseline go\n",
    "# contrast_betas = contrast_betas.xs('fs',level=1)\n",
    "contrast_betas = contrast_betas.melt(var_name='ROI', value_name='beta_go', ignore_index=False)\n",
    "if atlas_type == 'ATAG':\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[0])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[1:])\n",
    "else:\n",
    "    contrast_betas['hemisphere'] = contrast_betas['ROI'].apply(lambda x: x[-1])\n",
    "    contrast_betas['ROI_nohemi'] = contrast_betas['ROI'].apply(lambda x: x[:-2])\n",
    "betas_go = contrast_betas.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# look at significance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import display_html\n",
    "from itertools import chain,cycle\n",
    "def display_side_by_side(*args,titles=cycle([''])):\n",
    "    html_str=''\n",
    "    for df,title in zip(args, chain(titles,cycle(['</br>'])) ):\n",
    "        html_str+='<th style=\"text-align:center\"><td style=\"vertical-align:top\">'\n",
    "        html_str+=f'<h2>{title}</h2>'\n",
    "        html_str+=df.to_html().replace('table','table style=\"display:inline\"')\n",
    "        html_str+='</td></th>'\n",
    "    display_html(html_str,raw=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# absolute activation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = betas_fs.groupby('ROI')['beta_fs'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p']).reindex(rois_)\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "# stat_df.loc[stat_df['p_corrected']<=0.0001,'significance']='****'\n",
    "stat_fs = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "tmp = betas_ss.groupby('ROI')['beta_ss'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p']).reindex(rois_)\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "# stat_df.loc[stat_df['p_corrected']<=0.0001,'significance']='****'\n",
    "stat_ss = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "tmp = betas_go.groupby('ROI')['beta_go'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p']).reindex(rois_)\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "# stat_df.loc[stat_df['p_corrected']<=0.0001,'significance']='****'\n",
    "stat_go = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "stat_fs, stat_ss, stat_go = stat_fs.reindex(rois_), stat_ss.reindex(rois_), stat_go.reindex(rois_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display_side_by_side(stat_fs,stat_ss,stat_go, titles=['fs','ss','go'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# contrast significance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = contrast_betas_fsgo.groupby('ROI')['beta_fs-go'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p']).reindex(rois_)\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "# stat_df.loc[stat_df['p_corrected']<=0.0001,'significance']='****'\n",
    "stat_fsgo = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "tmp = contrast_betas_fsss.groupby('ROI')['beta_fs-ss'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p']).reindex(rois_)\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "# stat_df.loc[stat_df['p_corrected']<=0.0001,'significance']='****'\n",
    "stat_fsss = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "tmp = contrast_betas_ssgo.groupby('ROI')['beta_ss-go'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p']).reindex(rois_)\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "# stat_df.loc[stat_df['p_corrected']<=0.0001,'significance']='****'\n",
    "stat_ssgo = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "tmp = contrast_betas_motor.groupby('ROI')['beta_contra-ipsi'].apply(lambda x: scipy.stats.ttest_1samp(x, 0))\n",
    "stat_df = pd.DataFrame.from_dict(dict(zip(tmp.index, tmp.values)), orient='index', columns=['t', 'p']).reindex(rois_)\n",
    "stat_df = pd.concat([stat_df, pd.DataFrame(statsmodels.stats.multitest.fdrcorrection(stat_df['p'], method='i'), index=['fdr_significant', 'p_corrected'], columns=stat_df.index).T], axis=1)\n",
    "\n",
    "stat_df['significance'] = ''\n",
    "stat_df.loc[stat_df['p_corrected']<=0.05,'significance']='*'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.01,'significance']='**'\n",
    "stat_df.loc[stat_df['p_corrected']<=0.001,'significance']='***'\n",
    "# stat_df.loc[stat_df['p_corrected']<=0.0001,'significance']='****'\n",
    "stat_leftright = stat_df.loc[sorted(stat_df.index, key=str.casefold)]\n",
    "\n",
    "stat_fsgo, stat_fsss, stat_ssgo, stat_leftright = stat_fsgo.reindex(rois_), stat_fsss.reindex(rois_), stat_ssgo.reindex(rois_), stat_leftright.reindex(rois_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display_side_by_side(stat_leftright,stat_fsss,stat_fsgo,stat_ssgo, titles=['contra-ipsilateral','fs > ss','fs > go','ss > go'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## plot all 4 contrasts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats\n",
    "\n",
    "def mean_confidence_interval(data, confidence=0.95):\n",
    "    a = 1.0 * np.array(data)\n",
    "    n = len(a)\n",
    "    m, se = np.mean(a), scipy.stats.sem(a)\n",
    "    h = se * scipy.stats.t.ppf((1 + confidence) / 2., n-1)\n",
    "    #return m, m-h, m+h\n",
    "    return m+h # only need highest \n",
    "\n",
    "def mean_confidence_interval_lower(data, confidence=0.95):\n",
    "    a = 1.0 * np.array(data)\n",
    "    n = len(a)\n",
    "    m, se = np.mean(a), scipy.stats.sem(a)\n",
    "    h = se * scipy.stats.t.ppf((1 + confidence) / 2., n-1)\n",
    "    #return m, m-h, m+h\n",
    "    return m-h # only need highest \n",
    "\n",
    "def add_asteriks(ax_n, data, col_name, significance=0):\n",
    "    l_count, r_count = 0, 0\n",
    "    \n",
    "    \n",
    "#     new = contrast_betas_inccon.groupby(['ROI_nohemi','hemisphere'])['beta_inc-con'].apply(mean_confidence_interval).reset_index().set_index('ROI_nohemi')\n",
    "#     new = new.loc[sorted(new.index,key=str.casefold)].drop_duplicates().reset_index().set_index(['ROI_nohemi','hemisphere'])\n",
    "    roi_data = data.groupby(['ROI_nohemi','hemisphere'])[col_name].apply(mean_confidence_interval).reset_index().set_index('ROI_nohemi')\n",
    "#     roi_data = roi_data.loc[sorted(roi_data.index,key=str.casefold)].drop_duplicates().reset_index().set_index(['ROI_nohemi','hemisphere'])\n",
    "    roi_data = roi_data.loc[roi_data.index].drop_duplicates().reset_index().set_index(['ROI_nohemi','hemisphere']).reindex(gm_nuclei,level=0)\n",
    "    \n",
    "    roi_data_lower = data.groupby(['ROI_nohemi','hemisphere'])[col_name].apply(mean_confidence_interval_lower).reset_index().set_index('ROI_nohemi')\n",
    "#     roi_data_lower = roi_data_lower.loc[sorted(roi_data_lower.index,key=str.casefold)].drop_duplicates().reset_index().set_index(['ROI_nohemi','hemisphere'])\n",
    "    roi_data_lower = roi_data_lower.loc[roi_data_lower.index].drop_duplicates().reset_index().set_index(['ROI_nohemi','hemisphere']).reindex(gm_nuclei,level=0)\n",
    "    \n",
    "#     print(roi_data_lower)\n",
    "#     for index, (y_val, region) in enumerate(zip(data.groupby(['ROI_nohemi','hemisphere'])[col_name].apply(mean_confidence_interval), data.ROI.unique())):\n",
    "    for index, (y_val, region, y_val_lower) in enumerate(zip(roi_data[col_name], data.ROI.unique(), roi_data_lower[col_name])):\n",
    "#         print(significance)\n",
    "#         print(region)\n",
    "        if index%2 == 0:\n",
    "            x_val = -0.26 + l_count\n",
    "            l_count+=1\n",
    "            if (aster:=len(significance[index])) == 2:\n",
    "                x_val = x_val - 0.06\n",
    "            elif aster == 3:\n",
    "                x_val = x_val - 0.12\n",
    "        else:\n",
    "            x_val = 0.14 + r_count\n",
    "            r_count+=1\n",
    "            if (aster:=len(significance[index])) == 2:\n",
    "                x_val = x_val - 0.06\n",
    "            elif aster == 3:\n",
    "                x_val = x_val - 0.12\n",
    "        if y_val >= 0:\n",
    "            ax[ax_n].text(x_val, y_val+0.005, significance[index])\n",
    "        else:\n",
    "            ax[ax_n].text(x_val, y_val_lower-0.017, significance[index])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "f, axes = plt.subplots(2,2, figsize=(30,17), sharex=False,sharey=True)\n",
    "ax = axes.ravel()\n",
    "f.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.1, hspace=0.27)\n",
    "# sns.set_style('darkgrid')\n",
    "\n",
    "sns.set_theme() # to make style changable from defaults use this line of code befor using set_style\n",
    "with sns.axes_style(\"darkgrid\"):\n",
    "\n",
    "#     # 1 ## MOTOR RESPONSE\n",
    "#     sns.barplot(x='ROI_nohemi', y='beta_contra-ipsi', hue='hemisphere', \n",
    "#                 data=contrast_betas_motor, ax=ax[0],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "#     ax[0].set_ylabel('% change', fontsize=22)\n",
    "#     ax[0].set_xlabel('')\n",
    "#     ax[0].legend_.remove()\n",
    "#     ax[0].set_title('Contralateral - ipsilateral', fontsize=26)\n",
    "#     ax[0].set(ylim=(-0.17,0.25))\n",
    "#     add_asteriks(0, contrast_betas_motor, 'beta_contra-ipsi', stat_leftright['significance'])\n",
    "#     # ax[0].set_xticklabels(ax[0].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "#     ax[0].set_yticklabels(ax[0].get_yticks().round(2),fontsize=20)\n",
    "#     ax[0].set_xticklabels(ax[0].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "\n",
    "    ## 2\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_fs-go', hue='hemisphere', \n",
    "                data=contrast_betas_fsgo, ax=ax[0],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    ax[0].set_ylabel('% change', fontsize=22)\n",
    "    ax[0].set_xlabel('')\n",
    "    ax[0].legend_.remove()\n",
    "    ax[0].set_title('FS - GO', fontsize=26)\n",
    "    ax[0].set(ylim=(-0.17,0.25))\n",
    "    add_asteriks(0, contrast_betas_fsgo, 'beta_fs-go', stat_fsgo['significance'])\n",
    "    ax[0].set_xticklabels(ax[0].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    # ax[1].set_xticklabels(ax[1].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    ax[0].set_yticklabels(ax[0].get_yticks().round(2),fontsize=20)\n",
    "\n",
    "    ## 3\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_fs-ss', hue='hemisphere', \n",
    "                data=contrast_betas_fsss, ax=ax[1],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    ax[1].set_ylabel('')\n",
    "    ax[1].set_xlabel('')\n",
    "    ax[1].legend_.remove()\n",
    "    ax[1].set_title('FS - SS', fontsize=26)\n",
    "    ax[1].set(ylim=(-0.17,0.25))\n",
    "    add_asteriks(1, contrast_betas_fsss, 'beta_fs-ss', stat_fsss['significance'])\n",
    "    ax[1].set_xticklabels(ax[1].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "#     ax[1].set_yticklabels(ax[1].get_yticks().round(2),fontsize=20)\n",
    "\n",
    "    ## 4\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_ss-go', hue='hemisphere', \n",
    "                data=contrast_betas_ssgo, ax=ax[2],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    ax[2].set_ylabel('% change', fontsize=22)\n",
    "    ax[2].set_xlabel('ROI', fontsize=22)\n",
    "#     ax[3].legend_.remove()\n",
    "    ax[2].set_title('SS - GO', fontsize=26)\n",
    "    ax[2].set(ylim=(-0.17,0.25))\n",
    "    add_asteriks(2, contrast_betas_ssgo, 'beta_ss-go', stat_ssgo['significance'])\n",
    "    ax[2].set_xticklabels(ax[2].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    ax[2].set_yticklabels(ax[2].get_yticks().round(2),fontsize=20)\n",
    "    plt.setp(ax[2].get_legend().get_texts(), fontsize='22') # for legend text\n",
    "    plt.setp(ax[2].get_legend().get_title(), fontsize='18') # for legend title\n",
    "\n",
    "    for x in range(4) :[t.set_color(i) for (i,t) in zip(['orange']*7+['red']*8,ax[x].get_xticklabels())]\n",
    "        \n",
    "    ax[3].axis('off')\n",
    "\n",
    "\n",
    "    f.savefig('figure_download_scott/GLM_ROI_SST_contrasts.pdf', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, axes = plt.subplots(2,2, figsize=(30,17), sharex=False,sharey=True)\n",
    "ax = axes.ravel()\n",
    "f.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.1, hspace=0.27)\n",
    "# sns.set_style('darkgrid')\n",
    "\n",
    "with sns.axes_style(\"whitegrid\"):\n",
    "\n",
    "    ## 1 ## FS\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_fs', hue='hemisphere', \n",
    "                data=betas_fs, ax=ax[0],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    ax[0].set_ylabel('% change', fontsize=22)\n",
    "    ax[0].set_xlabel('')\n",
    "    ax[0].legend_.remove()\n",
    "    ax[0].set_title('FS', fontsize=26)\n",
    "    ax[0].set(ylim=(-0.15,0.35))\n",
    "    add_asteriks(0, betas_fs, 'beta_fs', stat_fs['significance'])\n",
    "    # ax[0].set_xticklabels(ax[0].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    ax[0].set_yticklabels(ax[0].get_yticks().round(2),fontsize=20)\n",
    "    ax[0].set_xticklabels(ax[0].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "\n",
    "    ## 2 SS\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_ss', hue='hemisphere', \n",
    "                data=betas_ss, ax=ax[1],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    ax[1].set_ylabel('')\n",
    "    ax[1].set_xlabel('')\n",
    "    ax[1].legend_.remove()\n",
    "    ax[1].set_title('SS', fontsize=26)\n",
    "    ax[1].set(ylim=(-0.15,0.35))\n",
    "    add_asteriks(1, betas_ss, 'beta_ss', stat_ss['significance'])\n",
    "    ax[1].set_xticklabels(ax[1].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "\n",
    "    ## 3 GO\n",
    "    sns.barplot(x='ROI_nohemi', y='beta_go', hue='hemisphere', \n",
    "                data=betas_go, ax=ax[2],palette=(\"darkblue\",\"cornflowerblue\"))\n",
    "    ax[2].set_ylabel('% change', fontsize=22)\n",
    "    ax[2].set_xlabel('', fontsize=22)\n",
    "#     ax[2].legend_.remove()\n",
    "    ax[2].set_title('GO', fontsize=26)\n",
    "    ax[2].set(ylim=(-0.15,0.35))\n",
    "    add_asteriks(2, betas_go, 'beta_go', stat_go['significance'])\n",
    "    # ax[2].set_xticklabels(ax[2].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    ax[2].set_yticklabels(ax[2].get_yticks().round(2),fontsize=20)\n",
    "    ax[2].set_xticklabels(ax[2].get_xticklabels(),rotation = 30,fontsize=20)\n",
    "    plt.setp(ax[2].get_legend().get_texts(), fontsize='22') # for legend text\n",
    "    plt.setp(ax[2].get_legend().get_title(), fontsize='18') # for legend title\n",
    "    \n",
    "    for x in range(4) :[t.set_color(i) for (i,t) in zip(['orange']*7+['red']*8,ax[x].get_xticklabels())]\n",
    "        \n",
    "    ax[3].set_visible(False)\n",
    "    \n",
    "    f.savefig('figure_download_scott/GLM_ROI_SST_supplementary.pdf', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py38",
   "language": "python",
   "name": "py38"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}