Downloads: 1
[d4b612]:
/
ML_alg.ipynb
2061 lines (2060 with data), 140.6 kB
{
"cells": [
{
"cell_type": "markdown",
"id": "4b74d181-ef9b-4a68-aa78-9c3a69b59c61",
"metadata": {},
"source": [
"# ML algorithm to assess improvement after treatment"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "aa3fc6e2-913b-482f-a024-9aec6a5c72b7",
"metadata": {},
"outputs": [],
"source": [
"import warnings\n",
"warnings.filterwarnings('ignore')\n",
"%config Completer.use_jedi = False\n",
"\n",
"import glob\n",
"import numpy as np\n",
"import pandas as pd\n",
"\n",
"import scipy\n",
"import nilearn\n",
"import nilearn.plotting\n",
"import nilearn.input_data\n",
"import matplotlib.pyplot as plt\n",
"import seaborn as sns\n",
"from xgboost import XGBClassifier\n",
"import xgboost as xg\n",
"from sklearn.model_selection import cross_val_score\n",
"from sklearn.model_selection import StratifiedKFold\n",
"from sklearn.model_selection import train_test_split\n",
"from sklearn.metrics import accuracy_score\n",
"#import shap\n",
"from sklearn.model_selection import permutation_test_score\n",
"\n",
"\n",
"from sklearn.linear_model import Ridge\n",
"from sklearn.linear_model import RidgeCV"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "72526c7e-4d6e-4c73-8da4-51055748c00d",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>scr_id</th>\n",
" <th>group</th>\n",
" <th>treatment</th>\n",
" <th>age</th>\n",
" <th>sex</th>\n",
" <th>pre_pcl4_score</th>\n",
" <th>post_pcl4_score</th>\n",
" <th>delta_pcl4_score</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>CPT2001</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>30</td>\n",
" <td>1</td>\n",
" <td>70.0</td>\n",
" <td>74.0</td>\n",
" <td>4.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>CPT2002</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>39</td>\n",
" <td>0</td>\n",
" <td>71.0</td>\n",
" <td>31.0</td>\n",
" <td>-40.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>CPT2005</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>29</td>\n",
" <td>0</td>\n",
" <td>62.0</td>\n",
" <td>40.0</td>\n",
" <td>-22.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>CPT2013</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>32</td>\n",
" <td>0</td>\n",
" <td>45.0</td>\n",
" <td>68.0</td>\n",
" <td>23.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>CPT2017</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>30</td>\n",
" <td>0</td>\n",
" <td>53.0</td>\n",
" <td>55.0</td>\n",
" <td>2.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" scr_id group treatment age sex pre_pcl4_score post_pcl4_score \\\n",
"0 CPT2001 PTSDCOMB CPT 30 1 70.0 74.0 \n",
"1 CPT2002 PTSDCOMB CPT 39 0 71.0 31.0 \n",
"2 CPT2005 PTSDCOMB CPT 29 0 62.0 40.0 \n",
"3 CPT2013 PTSDCOMB CPT 32 0 45.0 68.0 \n",
"4 CPT2017 PTSDCOMB CPT 30 0 53.0 55.0 \n",
"\n",
" delta_pcl4_score \n",
"0 4.0 \n",
"1 -40.0 \n",
"2 -22.0 \n",
"3 23.0 \n",
"4 2.0 "
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# load PCL4 data\n",
"# sex: 0 = male, 1 = female\n",
"df_pcl = pd.read_csv('cpt_pclData.csv')\n",
"df_pcl.head()"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "69ea64df-a45c-45fa-a37d-b1e0cc953c4f",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"# remove CPT and add sub- to the dataframe, to align with BIDS\n",
"for i in df_pcl.iterrows():\n",
" sub = i[1].scr_id.split('CPT')\n",
" sub = 'sub-' + sub[1]\n",
" df_pcl.at[i[0], 'scr_id'] = sub"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "1b9dea21-d234-4184-a3bb-f89787525345",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"<class 'pandas.core.frame.DataFrame'>\n",
"RangeIndex: 43 entries, 0 to 42\n",
"Data columns (total 8 columns):\n",
" # Column Non-Null Count Dtype \n",
"--- ------ -------------- ----- \n",
" 0 scr_id 43 non-null object \n",
" 1 group 43 non-null object \n",
" 2 treatment 43 non-null object \n",
" 3 age 43 non-null int64 \n",
" 4 sex 43 non-null int64 \n",
" 5 pre_pcl4_score 41 non-null float64\n",
" 6 post_pcl4_score 41 non-null float64\n",
" 7 delta_pcl4_score 41 non-null float64\n",
"dtypes: float64(3), int64(2), object(3)\n",
"memory usage: 2.8+ KB\n"
]
}
],
"source": [
"df_pcl.info()"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "18d64725-3d9b-4d06-a8d9-e5b7f2610ca7",
"metadata": {},
"outputs": [],
"source": [
"# add responders / non responders to the data\n",
"responder = []\n",
"for i in np.arange(len(df_pcl)):\n",
" delta = df_pcl.at[i, 'delta_pcl4_score']\n",
" if delta<=-10:\n",
" responder.append('yes')\n",
" else:\n",
" responder.append('no')\n",
"df_pcl['responder'] = responder"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "57bffe00-d1f9-47d0-8bce-318d63f9ce37",
"metadata": {
"tags": []
},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>scr_id</th>\n",
" <th>group</th>\n",
" <th>treatment</th>\n",
" <th>age</th>\n",
" <th>sex</th>\n",
" <th>pre_pcl4_score</th>\n",
" <th>post_pcl4_score</th>\n",
" <th>delta_pcl4_score</th>\n",
" <th>responder</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>sub-2001</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>30</td>\n",
" <td>1</td>\n",
" <td>70.0</td>\n",
" <td>74.0</td>\n",
" <td>4.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>sub-2002</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>39</td>\n",
" <td>0</td>\n",
" <td>71.0</td>\n",
" <td>31.0</td>\n",
" <td>-40.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>sub-2005</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>29</td>\n",
" <td>0</td>\n",
" <td>62.0</td>\n",
" <td>40.0</td>\n",
" <td>-22.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>sub-2013</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>32</td>\n",
" <td>0</td>\n",
" <td>45.0</td>\n",
" <td>68.0</td>\n",
" <td>23.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>sub-2017</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>30</td>\n",
" <td>0</td>\n",
" <td>53.0</td>\n",
" <td>55.0</td>\n",
" <td>2.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>sub-2027</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>39</td>\n",
" <td>0</td>\n",
" <td>57.0</td>\n",
" <td>32.0</td>\n",
" <td>-25.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>sub-2028</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>40</td>\n",
" <td>0</td>\n",
" <td>48.0</td>\n",
" <td>31.0</td>\n",
" <td>-17.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td>sub-2032</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>23</td>\n",
" <td>0</td>\n",
" <td>50.0</td>\n",
" <td>17.0</td>\n",
" <td>-33.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>8</th>\n",
" <td>sub-2038</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>37</td>\n",
" <td>0</td>\n",
" <td>66.0</td>\n",
" <td>62.0</td>\n",
" <td>-4.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>9</th>\n",
" <td>sub-2042</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>30</td>\n",
" <td>0</td>\n",
" <td>70.0</td>\n",
" <td>76.0</td>\n",
" <td>6.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>10</th>\n",
" <td>sub-2044</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>27</td>\n",
" <td>0</td>\n",
" <td>72.0</td>\n",
" <td>48.0</td>\n",
" <td>-24.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>11</th>\n",
" <td>sub-2045</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>28</td>\n",
" <td>1</td>\n",
" <td>49.0</td>\n",
" <td>37.0</td>\n",
" <td>-12.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>12</th>\n",
" <td>sub-2048</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>46</td>\n",
" <td>0</td>\n",
" <td>57.0</td>\n",
" <td>55.0</td>\n",
" <td>-2.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>13</th>\n",
" <td>sub-2052</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>25</td>\n",
" <td>0</td>\n",
" <td>43.0</td>\n",
" <td>28.0</td>\n",
" <td>-15.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>14</th>\n",
" <td>sub-2053</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>26</td>\n",
" <td>0</td>\n",
" <td>44.0</td>\n",
" <td>36.0</td>\n",
" <td>-8.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>15</th>\n",
" <td>sub-2054</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>22</td>\n",
" <td>0</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>16</th>\n",
" <td>sub-2055</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>28</td>\n",
" <td>0</td>\n",
" <td>40.0</td>\n",
" <td>50.0</td>\n",
" <td>10.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>17</th>\n",
" <td>sub-2058</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>35</td>\n",
" <td>0</td>\n",
" <td>59.0</td>\n",
" <td>55.0</td>\n",
" <td>-4.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>18</th>\n",
" <td>sub-2059</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>41</td>\n",
" <td>0</td>\n",
" <td>59.0</td>\n",
" <td>46.0</td>\n",
" <td>-13.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>19</th>\n",
" <td>sub-2063</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>29</td>\n",
" <td>0</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>NaN</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>20</th>\n",
" <td>sub-2004</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>39</td>\n",
" <td>0</td>\n",
" <td>67.0</td>\n",
" <td>65.0</td>\n",
" <td>-2.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>21</th>\n",
" <td>sub-2008</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>51</td>\n",
" <td>0</td>\n",
" <td>55.0</td>\n",
" <td>21.0</td>\n",
" <td>-34.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>22</th>\n",
" <td>sub-2010</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>38</td>\n",
" <td>0</td>\n",
" <td>52.0</td>\n",
" <td>52.0</td>\n",
" <td>0.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>23</th>\n",
" <td>sub-2012</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>24</td>\n",
" <td>0</td>\n",
" <td>55.0</td>\n",
" <td>57.0</td>\n",
" <td>2.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>24</th>\n",
" <td>sub-2015</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>34</td>\n",
" <td>0</td>\n",
" <td>59.0</td>\n",
" <td>63.0</td>\n",
" <td>4.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>25</th>\n",
" <td>sub-2021</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>41</td>\n",
" <td>0</td>\n",
" <td>55.0</td>\n",
" <td>42.0</td>\n",
" <td>-13.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>26</th>\n",
" <td>sub-2022</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>27</td>\n",
" <td>0</td>\n",
" <td>49.0</td>\n",
" <td>31.0</td>\n",
" <td>-18.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>27</th>\n",
" <td>sub-2023</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>32</td>\n",
" <td>0</td>\n",
" <td>67.0</td>\n",
" <td>52.0</td>\n",
" <td>-15.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>28</th>\n",
" <td>sub-2024</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>36</td>\n",
" <td>0</td>\n",
" <td>63.0</td>\n",
" <td>67.0</td>\n",
" <td>4.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>29</th>\n",
" <td>sub-2025</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>27</td>\n",
" <td>0</td>\n",
" <td>54.0</td>\n",
" <td>18.0</td>\n",
" <td>-36.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>30</th>\n",
" <td>sub-2026</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>28</td>\n",
" <td>0</td>\n",
" <td>70.0</td>\n",
" <td>73.0</td>\n",
" <td>3.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>31</th>\n",
" <td>sub-2033</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>54</td>\n",
" <td>0</td>\n",
" <td>43.0</td>\n",
" <td>24.0</td>\n",
" <td>-19.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>32</th>\n",
" <td>sub-2034</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>25</td>\n",
" <td>0</td>\n",
" <td>55.0</td>\n",
" <td>46.0</td>\n",
" <td>-9.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>33</th>\n",
" <td>sub-2036</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>36</td>\n",
" <td>0</td>\n",
" <td>80.0</td>\n",
" <td>83.0</td>\n",
" <td>3.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>34</th>\n",
" <td>sub-2037</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>25</td>\n",
" <td>0</td>\n",
" <td>55.0</td>\n",
" <td>33.0</td>\n",
" <td>-22.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>35</th>\n",
" <td>sub-2039</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>45</td>\n",
" <td>0</td>\n",
" <td>63.0</td>\n",
" <td>24.0</td>\n",
" <td>-39.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>36</th>\n",
" <td>sub-2043</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>33</td>\n",
" <td>0</td>\n",
" <td>35.0</td>\n",
" <td>19.0</td>\n",
" <td>-16.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>37</th>\n",
" <td>sub-2047</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>20</td>\n",
" <td>0</td>\n",
" <td>70.0</td>\n",
" <td>60.0</td>\n",
" <td>-10.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>38</th>\n",
" <td>sub-2050</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>24</td>\n",
" <td>0</td>\n",
" <td>31.0</td>\n",
" <td>21.0</td>\n",
" <td>-10.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>39</th>\n",
" <td>sub-2051</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>31</td>\n",
" <td>0</td>\n",
" <td>56.0</td>\n",
" <td>44.0</td>\n",
" <td>-12.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>40</th>\n",
" <td>sub-2056</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>26</td>\n",
" <td>1</td>\n",
" <td>50.0</td>\n",
" <td>31.0</td>\n",
" <td>-19.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>41</th>\n",
" <td>sub-2062</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>21</td>\n",
" <td>0</td>\n",
" <td>55.0</td>\n",
" <td>52.0</td>\n",
" <td>-3.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>42</th>\n",
" <td>sub-2064</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>35</td>\n",
" <td>0</td>\n",
" <td>66.0</td>\n",
" <td>57.0</td>\n",
" <td>-9.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" scr_id group treatment age sex pre_pcl4_score post_pcl4_score \\\n",
"0 sub-2001 PTSDCOMB CPT 30 1 70.0 74.0 \n",
"1 sub-2002 PTSDCOMB CPT 39 0 71.0 31.0 \n",
"2 sub-2005 PTSDCOMB CPT 29 0 62.0 40.0 \n",
"3 sub-2013 PTSDCOMB CPT 32 0 45.0 68.0 \n",
"4 sub-2017 PTSDCOMB CPT 30 0 53.0 55.0 \n",
"5 sub-2027 PTSDCOMB CPT 39 0 57.0 32.0 \n",
"6 sub-2028 PTSDCOMB CPT 40 0 48.0 31.0 \n",
"7 sub-2032 PTSDCOMB CPT 23 0 50.0 17.0 \n",
"8 sub-2038 PTSDCOMB CPT 37 0 66.0 62.0 \n",
"9 sub-2042 PTSDCOMB CPT 30 0 70.0 76.0 \n",
"10 sub-2044 PTSDCOMB CPT 27 0 72.0 48.0 \n",
"11 sub-2045 PTSDCOMB CPT 28 1 49.0 37.0 \n",
"12 sub-2048 PTSDCOMB CPT 46 0 57.0 55.0 \n",
"13 sub-2052 PTSDCOMB CPT 25 0 43.0 28.0 \n",
"14 sub-2053 PTSDCOMB CPT 26 0 44.0 36.0 \n",
"15 sub-2054 PTSDCOMB CPT 22 0 NaN NaN \n",
"16 sub-2055 PTSDCOMB CPT 28 0 40.0 50.0 \n",
"17 sub-2058 PTSDCOMB CPT 35 0 59.0 55.0 \n",
"18 sub-2059 PTSDCOMB CPT 41 0 59.0 46.0 \n",
"19 sub-2063 PTSDCOMB CPT 29 0 NaN NaN \n",
"20 sub-2004 PTSDCOMB PCT 39 0 67.0 65.0 \n",
"21 sub-2008 PTSDCOMB PCT 51 0 55.0 21.0 \n",
"22 sub-2010 PTSDCOMB PCT 38 0 52.0 52.0 \n",
"23 sub-2012 PTSDCOMB PCT 24 0 55.0 57.0 \n",
"24 sub-2015 PTSDCOMB PCT 34 0 59.0 63.0 \n",
"25 sub-2021 PTSDCOMB PCT 41 0 55.0 42.0 \n",
"26 sub-2022 PTSDCOMB PCT 27 0 49.0 31.0 \n",
"27 sub-2023 PTSDCOMB PCT 32 0 67.0 52.0 \n",
"28 sub-2024 PTSDCOMB PCT 36 0 63.0 67.0 \n",
"29 sub-2025 PTSDCOMB PCT 27 0 54.0 18.0 \n",
"30 sub-2026 PTSDCOMB PCT 28 0 70.0 73.0 \n",
"31 sub-2033 PTSDCOMB PCT 54 0 43.0 24.0 \n",
"32 sub-2034 PTSDCOMB PCT 25 0 55.0 46.0 \n",
"33 sub-2036 PTSDCOMB PCT 36 0 80.0 83.0 \n",
"34 sub-2037 PTSDCOMB PCT 25 0 55.0 33.0 \n",
"35 sub-2039 PTSDCOMB PCT 45 0 63.0 24.0 \n",
"36 sub-2043 PTSDCOMB PCT 33 0 35.0 19.0 \n",
"37 sub-2047 PTSDCOMB PCT 20 0 70.0 60.0 \n",
"38 sub-2050 PTSDCOMB PCT 24 0 31.0 21.0 \n",
"39 sub-2051 PTSDCOMB PCT 31 0 56.0 44.0 \n",
"40 sub-2056 PTSDCOMB PCT 26 1 50.0 31.0 \n",
"41 sub-2062 PTSDCOMB PCT 21 0 55.0 52.0 \n",
"42 sub-2064 PTSDCOMB PCT 35 0 66.0 57.0 \n",
"\n",
" delta_pcl4_score responder \n",
"0 4.0 no \n",
"1 -40.0 yes \n",
"2 -22.0 yes \n",
"3 23.0 no \n",
"4 2.0 no \n",
"5 -25.0 yes \n",
"6 -17.0 yes \n",
"7 -33.0 yes \n",
"8 -4.0 no \n",
"9 6.0 no \n",
"10 -24.0 yes \n",
"11 -12.0 yes \n",
"12 -2.0 no \n",
"13 -15.0 yes \n",
"14 -8.0 no \n",
"15 NaN no \n",
"16 10.0 no \n",
"17 -4.0 no \n",
"18 -13.0 yes \n",
"19 NaN no \n",
"20 -2.0 no \n",
"21 -34.0 yes \n",
"22 0.0 no \n",
"23 2.0 no \n",
"24 4.0 no \n",
"25 -13.0 yes \n",
"26 -18.0 yes \n",
"27 -15.0 yes \n",
"28 4.0 no \n",
"29 -36.0 yes \n",
"30 3.0 no \n",
"31 -19.0 yes \n",
"32 -9.0 no \n",
"33 3.0 no \n",
"34 -22.0 yes \n",
"35 -39.0 yes \n",
"36 -16.0 yes \n",
"37 -10.0 yes \n",
"38 -10.0 yes \n",
"39 -12.0 yes \n",
"40 -19.0 yes \n",
"41 -3.0 no \n",
"42 -9.0 no "
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df_pcl"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "42b9370d-98df-403e-9c15-77d6c2d64aec",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<AxesSubplot:xlabel='delta_pcl4_score', ylabel='Count'>"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEHCAYAAACk6V2yAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAY4UlEQVR4nO3df5RU5Z3n8fenoaVRMQkILrFpG5kclECm1dJJxGQnaEDFhdHFCR6S1V1cJJPEH2zMIu5EdzeTYzZzHIxJZDoxIbsyOkbFTTJrQlTQmTOO0ihgK/gDI00TRWQUMCih5bt/1AWroX8UTd2uqsvndU4fqu6v5/tgnw/Xp249jyICMzPLnppyF2BmZulwwJuZZZQD3swsoxzwZmYZ5YA3M8uogeUuoNDxxx8fjY2N5S7DzKxqrFq16s2IGN7VvooK+MbGRlpaWspdhplZ1ZC0sbt9HqIxM8soB7yZWUY54M3MMqqixuDNzPpqz549tLe3895775W7lFTU1dVRX19PbW1t0ec44M0sE9rb2xkyZAiNjY1IKnc5JRURbNu2jfb2dkaPHl30eR6iMbNMeO+99xg2bFjmwh1AEsOGDTvk/ztJNeAlXSfpOUmtku6WVJdme2Z2ZMtiuO/Tl76lFvCSTgSuBnIRMR4YAMxMqz0zM+ss7SGagcBgSQOBo4HfpdyemRkAb7/9Nj/4wQ8q7lr92U5qAR8Rm4G/BtqA14DtEbHswOMkzZHUIqll69ataZVjVpEaG+qRVBU/jQ315f7rOiTdheX7779fsmuVWqnbSe0pGkkfAaYDo4G3gZ9J+kJE3FV4XEQ0A80AuVzOy0vZEWXjps3Eo98qdxlF0aQF5S7hkMyfP58NGzbQ1NREbW0txx57LCNHjmT16tU8++yzzJ8/nxUrVrB7926+/OUvc9VVV/HOO+8wffp03nrrLfbs2cM3v/lNpk+f3ulan/vc55g6dSo33XQTJ5xwAqtXr+aSSy5hwoQJ3Hbbbbz77rs8+OCDjBkzhq1btzJ37lza2toAWLhwIRMnTuTmm2+mra2NV155hba2Nq699lquvvrqg9r5zne+c1h/B2k+Jnke8NuI2Aog6QHgbOCuHs8yMyuBW265hdbWVlavXs2KFSuYOnUqra2tjB49mubmZj70oQ+xcuVKdu/ezcSJE5k8eTKjRo1i6dKlHHfccbz55pt88pOfZNq0aZ2uBbBixQrWrFnDunXrGDp0KCeffDJXXnklTz31FLfddhu33347Cxcu5JprruG6667jnHPOoa2tjSlTprBu3ToA1q9fz/Lly9m5cydjx47lS1/60kHtHK40A74N+KSko4F3gXMBzyRmZmVx1lln7X+GfNmyZaxdu5b77rsPgO3bt/PSSy9RX1/PggULePzxx6mpqWHz5s1s2bKly+udeeaZjBw5EoAxY8YwefJkACZMmMDy5csBePjhh3n++ef3n7Njxw527twJwNSpUxk0aBCDBg1ixIgR3bZzOFIL+Ih4UtJ9wNNAB/AMyVCMmVl/O+aYY/a/jghuv/12pkyZ0umYxYsXs3XrVlatWkVtbS2NjY3dPns+aNCg/a9ramr2v6+pqaGjowOAvXv38sQTTzB48OAezx8wYMD+c0op1adoIuKmiDglIsZHxBcjYnea7ZmZ7TNkyJD9d8sHmjJlCnfccQd79uwB4MUXX+T3v/8927dvZ8SIEdTW1rJ8+XI2btzY67V6MnnyZL73ve/tf9/b0Etf2+mOv8lqZpk0bNgwJk6cyPjx47n++us77bvyyisZN24cp59+OuPHj+eqq66io6ODWbNm0dLSQi6XY8mSJZxyyim9Xqsn3/3ud2lpaeETn/gE48aNY9GiRX2uuS8UUTkPruRyufCCH3YkkVRVT9FUUl4caN26dZx66qnlLiNVXfVR0qqIyHV1vO/gzcwyygFvZpZRDngzs4xywJuZZZQD3swsoxzwZmYZ5YA3s8wa1XBSSWfUHNVwUq9tvv7668ycOZMxY8Ywbtw4LrzwQl588UUGDx5MU1MT48aNY+7cuaxZs4ampiaampoYOnQoo0ePpqmpifPOO69k/fearGaWWe2b2rh12Qslu968yWN73B8RXHzxxVx++eXcc889QP7bq1u2bGHMmDGsXr2ajo4OJk2axIYNG/Z/s/WKK67goosuYsaMGSWrFXwHb2ZWMsuXL6e2tpa5c+fu39bU1MSoUaP2vx84cCBnn302L7/8cur1OODNzEqktbWVM844o8djdu3axSOPPMKECRNSr8dDNGZm/WDfQh6SmD59OhdccEHqbTrgzcxK5OMf//j+OeYPtG8Mvj95iMbMrEQmTZrE7t27+eEPf7h/28qVK/dPO9zffAdvZplVP6qh1ydfDvV6PZHE0qVLufbaa7nllluoq6ujsbGRhQsXlqyGQ5Hmottjgb8v2HQy8I2IWJhWm2ZmhTa19f+d80c/+lHuvffeg7a3trZ2e87ixYtTqSXNJfteAJoAJA0ANgNL02rPzMw6668x+HOBDRFRnoEoM7MjUH8F/Ezg7n5qy8zM6IeAl3QUMA34WTf750hqkdSydevWtMsxMzti9Mcd/AXA0xGxpaudEdEcEbmIyA0fPrwfyjEzOzL0R8BfhodnzMz6XaoBL+lo4HPAA2m2Y2bWlcaG+pJOF9zYUN9rmwMGDKCpqYnx48dz6aWXsmvXLqDraYTTnjI41S86RcQuYFiabZiZdWfjps3Eo98q2fU0aUGvxwwePHj/lASzZs1i0aJFXHfddV1OI7xjx45Upwz2N1nNzFLy6U9/mrVr13Y7jXDaPBeNmVkKOjo6eOihh5gwYUJR0winwQFvZlZC7777Lk1NTeRyORoaGpg9e3bZavEQjZlZCRWOwe/T0zTCafIdvJlZyrqbRvixxx5LtV3fwZtZZp006sSinnw5lOv1RbmmEXbAm1lmvdrW3u9tvvPOO11u724a4X3SmDLYQzRmZhnlgDczyygHvJllRkSUu4TU9KVvDngzy4S6ujq2bduWyZCPCLZt20ZdXd0hnecPWc0sE+rr62lvbyer60rU1dVRX9/7ZGeFHPBmlgm1tbWMHj263GVUFA/RmJlllAPezCyjHPBmZhnlgDczy6i0l+z7sKT7JK2XtE7Sp9Jsz8zMPpD2UzS3Ab+KiBmSjgKOTrk9MzNLpBbwko4DPgNcARARfwD+kFZ7ZmbWWZpDNCcDW4GfSHpG0o8kHZNie2ZmViDNgB8InA7cERGnAb8H5h94kKQ5kloktWT1G2jVrrGhHklV8VN31MCy13AoP2ZpSnMMvh1oj4gnk/f30UXAR0Qz0AyQy+WyN4lEBmzctJl49FvlLqMomrSgamoFSroYhdmBUruDj4jXgU2SxiabzgWeT6s9MzPrLO2naL4KLEmeoHkF+I8pt2dmZolUAz4iVgO5NNswM7Ou+ZusZmYZ5YA3M8soB7yZWUY54M3MMsoBb2aWUQ54M7OMcsCbmWWUA97MLKMc8GZmGeWANzPLKAe8mVlGOeDNzDLKAW9mllEOeDOzjHLAm5lllAPezCyjHPBmZhmV6opOkl4FdgLvAx0R4dWdzMz6SdprsgJ8NiLe7Id2zMysgIdozMwyKu07+ACWSQrgbyOi+cADJM0B5gA0NDSkXI5ZhVENmrSg3FUUR74frDZpB/zEiPidpBHAbyStj4jHCw9IQr8ZIJfLRcr1mFWW2MutzXeWu4qizJszu9wl2CFK9Z/kiPhd8ucbwFLgrDTbMzOzD6QW8JKOkTRk32tgMtCaVntmZtZZmkM0JwBLJe1r5+8i4lcptmdmZgVSC/iIeAX447Sub2ZmPfPH4mZmGeWANzPLKAe8mVlGFRXwkiYWs83MzCpHsXfwtxe5zczMKkSPT9FI+hRwNjBc0ryCXccBA9IszMzMDk9vj0keBRybHDekYPsOYEZaRZmZ2eHrMeAj4jHgMUmLI2JjP9VkZmYlUOwXnQZJagYaC8+JiElpFGVmZoev2ID/GbAI+BH51ZnMzKzCFRvwHRFxR6qVmJlZSRX7mOQvJP2FpJGShu77SbUyMzM7LMXewV+e/Hl9wbYATi5tOWZmVipFBXxEjE67EDMzK62iAl7Sf+hqe0T879KWY2ZmpVLsEM2ZBa/rgHOBpwEHvJlZhSp2iOarhe8lfQj4P6lUZGZmJdHX6YJ3AR8r5kBJAyQ9I+mXfWzLzMz6oNgx+F+Qf2oG8pOMnQrcW2Qb1wDryE9QZmZm/aTYMfi/LnjdAWyMiPbeTpJUD0wF/gqY18vhZmZWQkUN0SSTjq0nP6PkR4A/FHn9hcDXgb3dHSBpjqQWSS1bt24t8rJmZtabYld0+nPgKeBS4M+BJyX1OF2wpIuANyJiVU/HRURzROQiIjd8+PAiyzYzs94UO0RzI3BmRLwBIGk48DBwXw/nTASmSbqQ/KOVx0m6KyK+cDgFm5lZcYp9iqZmX7gntvV2bkTcEBH1EdEIzAQedbibmfWfYu/gfyXp18DdyfvPA/8vnZLMzKwUeluT9Y+AEyLiekmXAOcAAp4AlhTbSESsAFb0vUwzMztUvQ3RLAR2AkTEAxExLyKuI3/3vjDd0szM7HD0FvCNEbH2wI0R0UJ++T4zM6tQvQV8XQ/7BpeyEDMzK63eAn6lpP984EZJs4Een283M7Py6u0pmmuBpZJm8UGg54CjgItTrMvMzA5TjwEfEVuAsyV9FhifbP6HiHg09crMzOywFDsf/HJgecq1mJlZCfV1PngzM6twDngzs4xywJuZZZQD3swsoxzwZmYZ5YA3M8soB7yZWUY54M3MMsoBb2aWUQ54M7OMSi3gJdVJekrSGknPSfrvabVlZmYHK3ZN1r7YDUyKiHck1QL/JOmhiPiXFNs0M7NEagEfEQG8k7ytTX4irfbMzKyzNO/gkTSA/DzyfwR8PyKe7OKYOcAcgIaGhj631dhQz8ZNm/t8fn86adSJvNrWXu4yiqcaNGlBuasoTjXVapayVAM+It4HmiR9mPzCIeMjovWAY5qBZoBcLtfnO/yNmzYTj37rcMrtN1UXQLGXW5vvLHcVRZk3Z3bV1Ar5es3S0i9P0UTE28AK4Pz+aM/MzNJ9imZ4cueOpMHAecD6tNozM7PO0hyiGQn8NBmHrwHujYhfptiemZkVSPMpmrXAaWld38zMeuZvspqZZZQD3swsoxzwZmYZ5YA3M8soB7yZWUY54M3MMsoBb2aWUQ54M7OMcsCbmWWUA97MLKMc8GZmGeWANzPLKAe8mVlGOeDNzDLKAW9mllEOeDOzjHLAm5llVJprso6StFzSOknPSbomrbbMzOxgaa7J2gH8l4h4WtIQYJWk30TE8ym2aWZmidTu4CPitYh4Onm9E1gHnJhWe2Zm1lmad/D7SWokvwD3k13smwPMAWhoaDiMRmrQpAV9P78/qQZJ5a7C7NBU2e9t/agGNrVtLHcZZZV6wEs6FrgfuDYidhy4PyKagWaAXC4XfW4o9nJr8519Pr0/zZszm1uXvVDuMoo2b/LYcpdglSD2+ve2yqT6FI2kWvLhviQiHkizLTMz6yzNp2gE3Amsi4hb02rHzMy6luYd/ETgi8AkSauTnwtTbM/MzAqkNgYfEf8EVM8nMmZmGeNvspqZZZQD3swsoxzwZmYZ5YA3M8soB7yZWUY54M3MMsoBb2aWUQ54M7OMcsCbmWWUA97MLKMc8GZmGeWANzPLKAe8mVlGOeDNzDLKAW9mllEOeDOzjHLAm5llVJprsv5Y0huSWtNqw8zMupfmHfxi4PwUr29mZj1ILeAj4nHgX9O6vpmZ9Sy1RbeLJWkOMAegoaGhzNWYWWaoBknlrqIo9aMa2NS2seTXLXvAR0Qz0AyQy+WizOWYWVbEXm5d9kK5qyjKvMljU7mun6IxM8soB7yZWUal+Zjk3cATwFhJ7ZJmp9WWmZkdLLUx+Ii4LK1rm5lZ7zxEY2aWUQ54M7OMcsCbmWWUA97MLKMc8GZmGeWANzPLKAe8mVlGOeDNzDLKAW9mllEOeDOzjHLAm5lllAPezCyjHPBmZhnlgDczyygHvJlZRjngzcwyygFvZpZRqQa8pPMlvSDpZUnz02zLzMw6S3NN1gHA94ELgHHAZZLGpdWemZl1luYd/FnAyxHxSkT8AbgHmJ5ie2ZmVkARkc6FpRnA+RFxZfL+i8CfRMRXDjhuDjAneTsWeKGPTR4PvNnHcyuB6y+/au+D6y+/cvThpIgY3tWOgSk2qi62HfSvSUQ0A82H3ZjUEhG5w71Oubj+8qv2Prj+8qu0PqQ5RNMOjCp4Xw/8LsX2zMysQJoBvxL4mKTRko4CZgI/T7E9MzMrkNoQTUR0SPoK8GtgAPDjiHgurfYowTBPmbn+8qv2Prj+8quoPqT2IauZmZWXv8lqZpZRDngzs4zKRMBL+pqkkHR8wbYbkikSXpA0pZz19UTS/5S0VtJqScskfbRgX8X3QdJ3JK1P+rBU0ocL9lVD/ZdKek7SXkm5A/ZVfP1QnVOCSPqxpDcktRZsGyrpN5JeSv78SDlr7ImkUZKWS1qX/P5ck2yvrD5ERFX/kH8U89fARuD4ZNs4YA0wCBgNbAAGlLvWbuo/ruD11cCiauoDMBkYmLz+NvDtKqv/VPJfsFsB5Aq2V0v9A5LaTgaOSmoeV+66iqj7M8DpQGvBtv8FzE9ez9/3u1SJP8BI4PTk9RDgxeR3pqL6kIU7+L8Bvk7nL1FNB+6JiN0R8VvgZfJTJ1SciNhR8PYYPuhHVfQhIpZFREfy9l/If98Bqqf+dRHR1benq6J+qnRKkIh4HPjXAzZPB36avP4p8Gf9WdOhiIjXIuLp5PVOYB1wIhXWh6oOeEnTgM0RseaAXScCmwretyfbKpKkv5K0CZgFfCPZXFV9SPwn4KHkdTXWX6ha6q+WOotxQkS8BvkABUaUuZ6iSGoETgOepML6kOZUBSUh6WHg33Sx60ZgAfkhgoNO62Jb2Z4H7akPEfF/I+JG4EZJNwBfAW6igvrQW/3JMTcCHcCSfad1cXzF1t/VaV1sq8RniqulzkySdCxwP3BtROyQuvrPUT4VH/ARcV5X2yVNID82uib5S60HnpZ0FhU2TUJ3fejC3wH/QD7gK6YPvdUv6XLgIuDcSAYfqaL6u1Ex9feiWuosxhZJIyPiNUkjgTfKXVBPJNWSD/clEfFAsrmi+lC1QzQR8WxEjIiIxohoJP+LfnpEvE5+SoSZkgZJGg18DHiqjOV2S9LHCt5OA9Ynr6uiD5LOB/4rMC0idhXsqor6e1At9WdpSpCfA5cnry8Huvu/q7JT/q7yTmBdRNxasKuy+lDuT6NL+Kn2qyRP0STvbyT/dMELwAXlrq+Huu8HWoG1wC+AE6upD+Q/fNwErE5+FlVZ/ReTvznYDWwBfl1N9Sd1Xkj+KY4N5Iedyl5TETXfDbwG7En+/mcDw4BHgJeSP4eWu84e6j+H/FDY2oLf/QsrrQ+eqsDMLKOqdojGzMx65oA3M8soB7yZWUY54M3MMsoBb2aWUQ54M7OMcsBbxZF0s6SvFbNf0hWFUyynXNerB0xJPUDSM5J+2R/tmx0qB7xVuyuAfgn4LlxDfhbBfiVpQH+3adXJAW8VQdKNyaIVD5Ofnx1JYyT9StIqSf8o6ZQDzpkB5IAlyYIpgyV9Q9JKSa2SmtXD7E+SVkhaKOmfk+PPSrYfK+knkp5NFjL5912cWw9MBX5URN+ulvR8cq17empD0mXJtlZJ3y64xjuS/oekJ4FPSfqCpKeSfv+tQ9+64oC3spN0Bvk5VE4DLgHOTHY1A1+NiDOArwE/KDwvIu4DWoBZEdEUEe8C34uIMyNiPDCY/CRoPTkmIs4G/gL4cbLtL4HtETEhIj4BPNrFeQvJr0Owt4guzgdOS641t7s2kqGmbwOTgCbgTEl/tq9O8otj/AmwDfg8MDEimoD3yU81bdZJxc8maUeETwNLI5msTNLPgTrgbOBnBTfhg4q41mclfR04GhgKPEd+jp/u3A35BSgkHaf8koPnkf8Hh2TfW4UnSLoIeCMiVkn60yJqWkv+/zIeBB5Mth3UhqTPACsiYmvSzhLyKx89SD7E708OPxc4A1iZ/N0MpsJnXrTycMBbpThwUqQa4O3kDrUokurI3+XnImKTpJvJ/0NxKO0G+TnWe5qkaSIwTdKFyfWPk3RXRHyhm+Onkg/qacBfSvp4N230NJn4exHxfsFxP42IG3o43sxDNFYRHgcuTsbQhwD/DtgF/FbSpZCfnlXSH3dx7k7ya2LCB2H+ZrIQw4wi2v58cv1zyA+ZbAeWkV94hWRfp4WTI+KGiKiP/DTVM4FHuwt3STXAqIhYTn5I58PAsd208STwbyUdn4ypXwY81sVlHwFmSBqRnDtU0klF9NWOMA54K7vIr2359+SnXL0f+Mdk1yxgtqQ15IdaulprdDGwSNJq8lP+/hB4lvywxsoimn9L0j8Di8hPWQvwTeAjyQeda4DPHnKnPjAAuEvSs8AzwN9ExNtdtRH5Jd5uAJaTXzz76ehixamIeB74b8AySWuB35BfBNqsE08XbEcsSSuAr0VES7lrMUuD7+DNzDLKH7Ja5kn6PvkPRgvdFhF/2k/t/KSU7ZgVy0M0ZmYZ5SEaM7OMcsCbmWWUA97MLKMc8GZmGfX/ATdK5T5GkSy1AAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"sns.histplot(x='delta_pcl4_score', hue='treatment', data=df_pcl)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "3b3313d7-999a-40cd-ad84-9ced147f234a",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>scr_id</th>\n",
" <th>group</th>\n",
" <th>treatment</th>\n",
" <th>age</th>\n",
" <th>sex</th>\n",
" <th>pre_pcl4_score</th>\n",
" <th>post_pcl4_score</th>\n",
" <th>delta_pcl4_score</th>\n",
" <th>responder</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>sub-2001</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>30</td>\n",
" <td>1</td>\n",
" <td>70.0</td>\n",
" <td>74.0</td>\n",
" <td>4.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>sub-2002</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>39</td>\n",
" <td>0</td>\n",
" <td>71.0</td>\n",
" <td>31.0</td>\n",
" <td>-40.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>sub-2005</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>29</td>\n",
" <td>0</td>\n",
" <td>62.0</td>\n",
" <td>40.0</td>\n",
" <td>-22.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>sub-2013</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>32</td>\n",
" <td>0</td>\n",
" <td>45.0</td>\n",
" <td>68.0</td>\n",
" <td>23.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>sub-2017</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>30</td>\n",
" <td>0</td>\n",
" <td>53.0</td>\n",
" <td>55.0</td>\n",
" <td>2.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>5</th>\n",
" <td>sub-2027</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>39</td>\n",
" <td>0</td>\n",
" <td>57.0</td>\n",
" <td>32.0</td>\n",
" <td>-25.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>sub-2028</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>40</td>\n",
" <td>0</td>\n",
" <td>48.0</td>\n",
" <td>31.0</td>\n",
" <td>-17.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>7</th>\n",
" <td>sub-2032</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>23</td>\n",
" <td>0</td>\n",
" <td>50.0</td>\n",
" <td>17.0</td>\n",
" <td>-33.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>8</th>\n",
" <td>sub-2038</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>37</td>\n",
" <td>0</td>\n",
" <td>66.0</td>\n",
" <td>62.0</td>\n",
" <td>-4.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>9</th>\n",
" <td>sub-2042</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>30</td>\n",
" <td>0</td>\n",
" <td>70.0</td>\n",
" <td>76.0</td>\n",
" <td>6.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>10</th>\n",
" <td>sub-2044</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>27</td>\n",
" <td>0</td>\n",
" <td>72.0</td>\n",
" <td>48.0</td>\n",
" <td>-24.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>11</th>\n",
" <td>sub-2045</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>28</td>\n",
" <td>1</td>\n",
" <td>49.0</td>\n",
" <td>37.0</td>\n",
" <td>-12.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>12</th>\n",
" <td>sub-2048</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>46</td>\n",
" <td>0</td>\n",
" <td>57.0</td>\n",
" <td>55.0</td>\n",
" <td>-2.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>13</th>\n",
" <td>sub-2052</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>25</td>\n",
" <td>0</td>\n",
" <td>43.0</td>\n",
" <td>28.0</td>\n",
" <td>-15.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>14</th>\n",
" <td>sub-2053</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>26</td>\n",
" <td>0</td>\n",
" <td>44.0</td>\n",
" <td>36.0</td>\n",
" <td>-8.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>16</th>\n",
" <td>sub-2055</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>28</td>\n",
" <td>0</td>\n",
" <td>40.0</td>\n",
" <td>50.0</td>\n",
" <td>10.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>17</th>\n",
" <td>sub-2058</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>35</td>\n",
" <td>0</td>\n",
" <td>59.0</td>\n",
" <td>55.0</td>\n",
" <td>-4.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>18</th>\n",
" <td>sub-2059</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>CPT</td>\n",
" <td>41</td>\n",
" <td>0</td>\n",
" <td>59.0</td>\n",
" <td>46.0</td>\n",
" <td>-13.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>20</th>\n",
" <td>sub-2004</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>39</td>\n",
" <td>0</td>\n",
" <td>67.0</td>\n",
" <td>65.0</td>\n",
" <td>-2.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>21</th>\n",
" <td>sub-2008</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>51</td>\n",
" <td>0</td>\n",
" <td>55.0</td>\n",
" <td>21.0</td>\n",
" <td>-34.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>22</th>\n",
" <td>sub-2010</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>38</td>\n",
" <td>0</td>\n",
" <td>52.0</td>\n",
" <td>52.0</td>\n",
" <td>0.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>23</th>\n",
" <td>sub-2012</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>24</td>\n",
" <td>0</td>\n",
" <td>55.0</td>\n",
" <td>57.0</td>\n",
" <td>2.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>24</th>\n",
" <td>sub-2015</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>34</td>\n",
" <td>0</td>\n",
" <td>59.0</td>\n",
" <td>63.0</td>\n",
" <td>4.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>25</th>\n",
" <td>sub-2021</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>41</td>\n",
" <td>0</td>\n",
" <td>55.0</td>\n",
" <td>42.0</td>\n",
" <td>-13.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>26</th>\n",
" <td>sub-2022</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>27</td>\n",
" <td>0</td>\n",
" <td>49.0</td>\n",
" <td>31.0</td>\n",
" <td>-18.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>27</th>\n",
" <td>sub-2023</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>32</td>\n",
" <td>0</td>\n",
" <td>67.0</td>\n",
" <td>52.0</td>\n",
" <td>-15.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>28</th>\n",
" <td>sub-2024</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>36</td>\n",
" <td>0</td>\n",
" <td>63.0</td>\n",
" <td>67.0</td>\n",
" <td>4.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>29</th>\n",
" <td>sub-2025</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>27</td>\n",
" <td>0</td>\n",
" <td>54.0</td>\n",
" <td>18.0</td>\n",
" <td>-36.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>30</th>\n",
" <td>sub-2026</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>28</td>\n",
" <td>0</td>\n",
" <td>70.0</td>\n",
" <td>73.0</td>\n",
" <td>3.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>31</th>\n",
" <td>sub-2033</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>54</td>\n",
" <td>0</td>\n",
" <td>43.0</td>\n",
" <td>24.0</td>\n",
" <td>-19.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>32</th>\n",
" <td>sub-2034</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>25</td>\n",
" <td>0</td>\n",
" <td>55.0</td>\n",
" <td>46.0</td>\n",
" <td>-9.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>33</th>\n",
" <td>sub-2036</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>36</td>\n",
" <td>0</td>\n",
" <td>80.0</td>\n",
" <td>83.0</td>\n",
" <td>3.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>34</th>\n",
" <td>sub-2037</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>25</td>\n",
" <td>0</td>\n",
" <td>55.0</td>\n",
" <td>33.0</td>\n",
" <td>-22.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>35</th>\n",
" <td>sub-2039</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>45</td>\n",
" <td>0</td>\n",
" <td>63.0</td>\n",
" <td>24.0</td>\n",
" <td>-39.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>36</th>\n",
" <td>sub-2043</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>33</td>\n",
" <td>0</td>\n",
" <td>35.0</td>\n",
" <td>19.0</td>\n",
" <td>-16.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>37</th>\n",
" <td>sub-2047</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>20</td>\n",
" <td>0</td>\n",
" <td>70.0</td>\n",
" <td>60.0</td>\n",
" <td>-10.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>38</th>\n",
" <td>sub-2050</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>24</td>\n",
" <td>0</td>\n",
" <td>31.0</td>\n",
" <td>21.0</td>\n",
" <td>-10.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>39</th>\n",
" <td>sub-2051</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>31</td>\n",
" <td>0</td>\n",
" <td>56.0</td>\n",
" <td>44.0</td>\n",
" <td>-12.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>40</th>\n",
" <td>sub-2056</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>26</td>\n",
" <td>1</td>\n",
" <td>50.0</td>\n",
" <td>31.0</td>\n",
" <td>-19.0</td>\n",
" <td>yes</td>\n",
" </tr>\n",
" <tr>\n",
" <th>41</th>\n",
" <td>sub-2062</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>21</td>\n",
" <td>0</td>\n",
" <td>55.0</td>\n",
" <td>52.0</td>\n",
" <td>-3.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" <tr>\n",
" <th>42</th>\n",
" <td>sub-2064</td>\n",
" <td>PTSDCOMB</td>\n",
" <td>PCT</td>\n",
" <td>35</td>\n",
" <td>0</td>\n",
" <td>66.0</td>\n",
" <td>57.0</td>\n",
" <td>-9.0</td>\n",
" <td>no</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" scr_id group treatment age sex pre_pcl4_score post_pcl4_score \\\n",
"0 sub-2001 PTSDCOMB CPT 30 1 70.0 74.0 \n",
"1 sub-2002 PTSDCOMB CPT 39 0 71.0 31.0 \n",
"2 sub-2005 PTSDCOMB CPT 29 0 62.0 40.0 \n",
"3 sub-2013 PTSDCOMB CPT 32 0 45.0 68.0 \n",
"4 sub-2017 PTSDCOMB CPT 30 0 53.0 55.0 \n",
"5 sub-2027 PTSDCOMB CPT 39 0 57.0 32.0 \n",
"6 sub-2028 PTSDCOMB CPT 40 0 48.0 31.0 \n",
"7 sub-2032 PTSDCOMB CPT 23 0 50.0 17.0 \n",
"8 sub-2038 PTSDCOMB CPT 37 0 66.0 62.0 \n",
"9 sub-2042 PTSDCOMB CPT 30 0 70.0 76.0 \n",
"10 sub-2044 PTSDCOMB CPT 27 0 72.0 48.0 \n",
"11 sub-2045 PTSDCOMB CPT 28 1 49.0 37.0 \n",
"12 sub-2048 PTSDCOMB CPT 46 0 57.0 55.0 \n",
"13 sub-2052 PTSDCOMB CPT 25 0 43.0 28.0 \n",
"14 sub-2053 PTSDCOMB CPT 26 0 44.0 36.0 \n",
"16 sub-2055 PTSDCOMB CPT 28 0 40.0 50.0 \n",
"17 sub-2058 PTSDCOMB CPT 35 0 59.0 55.0 \n",
"18 sub-2059 PTSDCOMB CPT 41 0 59.0 46.0 \n",
"20 sub-2004 PTSDCOMB PCT 39 0 67.0 65.0 \n",
"21 sub-2008 PTSDCOMB PCT 51 0 55.0 21.0 \n",
"22 sub-2010 PTSDCOMB PCT 38 0 52.0 52.0 \n",
"23 sub-2012 PTSDCOMB PCT 24 0 55.0 57.0 \n",
"24 sub-2015 PTSDCOMB PCT 34 0 59.0 63.0 \n",
"25 sub-2021 PTSDCOMB PCT 41 0 55.0 42.0 \n",
"26 sub-2022 PTSDCOMB PCT 27 0 49.0 31.0 \n",
"27 sub-2023 PTSDCOMB PCT 32 0 67.0 52.0 \n",
"28 sub-2024 PTSDCOMB PCT 36 0 63.0 67.0 \n",
"29 sub-2025 PTSDCOMB PCT 27 0 54.0 18.0 \n",
"30 sub-2026 PTSDCOMB PCT 28 0 70.0 73.0 \n",
"31 sub-2033 PTSDCOMB PCT 54 0 43.0 24.0 \n",
"32 sub-2034 PTSDCOMB PCT 25 0 55.0 46.0 \n",
"33 sub-2036 PTSDCOMB PCT 36 0 80.0 83.0 \n",
"34 sub-2037 PTSDCOMB PCT 25 0 55.0 33.0 \n",
"35 sub-2039 PTSDCOMB PCT 45 0 63.0 24.0 \n",
"36 sub-2043 PTSDCOMB PCT 33 0 35.0 19.0 \n",
"37 sub-2047 PTSDCOMB PCT 20 0 70.0 60.0 \n",
"38 sub-2050 PTSDCOMB PCT 24 0 31.0 21.0 \n",
"39 sub-2051 PTSDCOMB PCT 31 0 56.0 44.0 \n",
"40 sub-2056 PTSDCOMB PCT 26 1 50.0 31.0 \n",
"41 sub-2062 PTSDCOMB PCT 21 0 55.0 52.0 \n",
"42 sub-2064 PTSDCOMB PCT 35 0 66.0 57.0 \n",
"\n",
" delta_pcl4_score responder \n",
"0 4.0 no \n",
"1 -40.0 yes \n",
"2 -22.0 yes \n",
"3 23.0 no \n",
"4 2.0 no \n",
"5 -25.0 yes \n",
"6 -17.0 yes \n",
"7 -33.0 yes \n",
"8 -4.0 no \n",
"9 6.0 no \n",
"10 -24.0 yes \n",
"11 -12.0 yes \n",
"12 -2.0 no \n",
"13 -15.0 yes \n",
"14 -8.0 no \n",
"16 10.0 no \n",
"17 -4.0 no \n",
"18 -13.0 yes \n",
"20 -2.0 no \n",
"21 -34.0 yes \n",
"22 0.0 no \n",
"23 2.0 no \n",
"24 4.0 no \n",
"25 -13.0 yes \n",
"26 -18.0 yes \n",
"27 -15.0 yes \n",
"28 4.0 no \n",
"29 -36.0 yes \n",
"30 3.0 no \n",
"31 -19.0 yes \n",
"32 -9.0 no \n",
"33 3.0 no \n",
"34 -22.0 yes \n",
"35 -39.0 yes \n",
"36 -16.0 yes \n",
"37 -10.0 yes \n",
"38 -10.0 yes \n",
"39 -12.0 yes \n",
"40 -19.0 yes \n",
"41 -3.0 no \n",
"42 -9.0 no "
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# remove nans\n",
"df_pcl = df_pcl.dropna()\n",
"df_pcl"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "6b78f789-bb69-48d5-8536-38c05e82ac3f",
"metadata": {},
"outputs": [],
"source": [
"respond_list = list(df_pcl['scr_id'][df_pcl['responder']=='yes'])\n",
"resp_list = []\n",
"for subject in respond_list:\n",
" \n",
" sub = subject.split('sub-')[1]\n",
" resp_list.append(sub)\n",
"\n",
"\n",
"nonrespond_list = list(df_pcl['scr_id'][df_pcl['responder']=='no'])\n",
"non_list = []\n",
"for subject in nonrespond_list:\n",
" \n",
" sub = subject.split('sub-')[1]\n",
" non_list.append(sub)\n"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "e46926c3-1b83-48a6-8203-0d2344d935b0",
"metadata": {
"tags": []
},
"outputs": [],
"source": [
"resp_func = ['/gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_%s/modelestimate/results/zstat1.nii.gz'% (sub) for sub in resp_list]\n",
"nonResp_func = ['/gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_%s/modelestimate/results/zstat1.nii.gz'% (sub) for sub in non_list]"
]
},
{
"cell_type": "markdown",
"id": "82b831a8-d7ca-484b-87e0-4af610be750f",
"metadata": {},
"source": [
"# now lets grab trauma vs. neutral"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "784038dd-dfb2-4251-8a25-3b2bceb4c66b",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAekAAADJCAYAAAAHFcoVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABgkklEQVR4nO29eZhdVZX//b3zvTUPqVQqISGEQY0ICAotP2zRVwVtBd4WBWWmZRQEaUVAojIItmJESCcMIcxCABsivsJP7RZbRbuRwSAzEchEpkpV3bpD3fn947BWffepUxkgSd2qrM/z5Mmte8+wzzn77LXX2msI1Wq1GgzDMAzDqDvCY90AwzAMwzCCMSFtGIZhGHWKCWnDMAzDqFNMSBuGYRhGnWJC2jAMwzDqFBPShmEYhlGnmJA2DMMwjDrFhLRhGIZh1CkmpA3DMAyjTjEhbRiGYRh1iglpwzAMw6hTTEgbhmEYRp1iQtowDMMw6hQT0oZhGIZRp5iQNgzDMIw6xYS0YRiGYdQpJqQNwzAMo04xIW0YhmEYdYoJacMwDMOoU0xIG4ZhGNuc888/H+eff/5YN2PcEx3rBhiGYRgTj2eeeWasmzAhME3aMAxjBzJz5kz85je/GetmGG+DmTNnIpVKoampCVOmTMHJJ5+MTCazXc9pQtowDMMwtpCHH34YmUwGzzzzDJ5++mlcffXV2/V8JqQNwzAMYyuZMmUKDjvssO1u1jchbRiGYRhbycqVK/HII49gjz322K7nMSFtGIZhGFvIUUcdhebmZkyfPh2TJ0/GZZddtl3PZ0LaMAzDMLaQhx56CIODg3jsscfw4osvYsOGDdv1fCakDcMwDGMr+chHPoKTTz4ZX//617freSxO2jAMYwdTKpUwNDSkf0ejUUSjNhyPN84//3zMnDkTzzzzDPbbb7/tcg7TpA3DMHYwn/70p5FKpfTfd7/73bFukvE26Orqwoknnogrrrhiu50jVKvVatvt6IZhGMZOyaGHHgoAeOyxx8a0HeMds68YhmEYm+X2228HAHR2diKVSun3oudls1kceeSRmzzGkiVL0NjYCAAIhUL6fT6fR29vLwDgpJNO2qbtHu+YJm0YhmGMYPHixZgyZQoAIJFIIBz2VkfD4bAK2Gq1qttXq1X9+5lnnsHixYsBAMccc4yu14bDYT2O/A14gl72rVarKBQKAIA1a9bgmGOO2V6XOC6wNWnDMAzDqFNMk94GSDm2a6+9dkzbYRjjFXuH6oPrr78enZ2dAIDddtsN8Xhcf4tEIgCAxsZGxGIxAMOmbsAzd5955u8AAKVSGc8++2MAwAMPzEdbWxsAoKWlRc3dwLDJu1QqIZvNAgAqlYr+XiwW8dprrwEAent7ce655267ix0n2Jr0NsBKshnGO8PeoR3Hz372MwDA5MmTAQCxWExNzT09PbpdKBRSIVqpVHSbcrmMpqYmAF7omISSsRk8Ehn+nEgk9HM4HNZQs2QyiXK5DMBbky6VSnpumRCEQiGdKPT09ODxxx/X48j269atAwB87nOfe8f3ph7ZrkJ65syZWLhwIT7+8Y9vz9MY24GZM2di7dq1iEQiaGpqwuGHH4558+bpy2kYhgEAp576awBQgetnzpznVehyPPhPf/rZHdPAcY5p0saoPPzww/j4xz+ONWvW4LDDDsPVV1+N733ve2PdLMMwtpIHHngAANDa2ooZM2YAGHbaikQiamKeMmUKVq9eDcDTbgU2a5fLZdWMTzzxEd2XnchCoXDg52q1qtufeOIjus811xzgCHnZJp/P6zZTp05FIpHQNovJXa7n17/+NQYGBgAARx999FbcnfrGhLSxWXZUSTbDMLY9v/vd79RLOxwOq/lYhFwsFkMulwPgabryez6f122i0agKy3g8rsKSvbJZkFN0lfPZv71MFOLxuO4fDocdM7i0h7XwhoYGNXfL/42NjRoa9rvf/Q4f+chH3tb9qjfMu9vYLDuqJJthGIbhYpq0MSpHHXUUQqEQMpkMPvaxj233kmyGYWw7xMkqmUyqhlkqlRyzNOBprqIZF4tFtLa2AoB6WwOuht3e3q5rzF5stKcBb0mgkKdJ13Rf0Yy/+tU/4Yc/3F/bIJp0tVrV9kQiEW0nx1rL9SQSCdX88/m8Xv/BBx+82XbVM6ZJG6Oyo0uyGYaxbViyZIkmDmlsbERraytaW1vV0zoUCjnJR2KxGGKxGBobG5FMJpFMJtHc3Owco6enBz09PUgmk7qvt3/trX/DcDYx/uzhbV+tVlGr1VCr1RCJRPT4jY2Net7m5mZtj4R+iSCW88v1JBIJvU4+xpIlS7bfjd4BmJA2NsuOKslmGMb44Itf/Dm++MWfb9NjnnDCL3HCCb/cpsecCGx3c7eVZJsY7IiSbIZhbBuamprUJFwul9UU3dnZqU5iQiwW021DoZCOz0NDQ0gmkwC8mGoxcbPzV6lUCjRz83e1Ws3RpuW3Uqmkx4xEIrj//n8G4Hl2y/fcBja5+032gOdMJulEy+WyXtN4Dxvd7tLy05/+tPP3t771LVx55ZXb+7TGNoZLskkyBMMw6otf/OIXADyBxSFNmUwGAJBKpTT710gz9OiJR8LhsApXzq3NEwBOfuLP6c05uoOOwwU7otGoeqOXy2XHG1yENyPHK5fLep28rh0KhfS+fOYznxnt1tUt21VIv/7669vz8MZ2JOjZLViwYMc3xDCMuuO66z4EADjrrP/eZsc85piHAAD33ff/brNjTgTM7mwYhjEB+MMf/oBJkyYB8LRKTrPJCUfE8Uq03iCNelMMDQ3h7LN//9Z5KqhWK3ocTiMqFItFx1Q+HA8dcTRv2eYLX3hQBTXnDvcjx5H/i8WinpevPx6Po6GhAYB3jw455JCtut6xxoS0YRjGOObRRx8FALS1takwKhaLKBaLANwQq3g8HiicOaSJvxehx2vMg4ODZM6uqUncFfaht/4BkUgUkcjw8dnczfnAg0zZ/B23Qby6GU6ywolY4vG4CvtcLqf36/DDDx9xvnrEvLsNwzCMLUYqXRk7BtOkDcMwxjHivRyLxdQBq7W11XGcam5uBuA5ZQWZt0WjHU1zDfpbjucSeuv7YdN3NBrR74Ga45gWdMzFi49yzsmaP5vR9Yx6nqiW2eQ0o9FoVK0KpVJp3Hl7myZtGIYxTpk3b56atqvVKgYHBzE4OIh8Po9EIqEJPuQzm7U3Bycb8ZvGI5EIIhFvTVmShkSjUf0e8HJ2y27yfTQa1e1DoZB+7ze3y3mDQq1GQ8z6/mvO5/N6X6rVqt6vefPmbfGxxxLTpA3DMIwxZdGiT2j6T8PFhLRhGMZWcOihh451E5T169erFsqm7FAo5FSwCkomsrVe3c8/76UF9nJvi2c1b1Gjv2vI5bySl3/7208g5m7vlN7nvfYaFspf+MKdbzvJVdD11Go1NfdzwhX+vlqtagnPseaxxx4b9TcT0oZhGOOM9evXA3C9tcU8DUBN0QILqaACG1sTjhUKAXKIWi3YHM3Cu1YDQiEW6sMZx4I9w0dnS6/Df0zeViYvtVpN72NXV9cWnX8sMCFtGIaxFWxK69lR3HfffQCAqVOnqqArFAoqjBobGzVrV0tLiwqsUqk0Ii1oKpXSEKWgMChBjn300T9zajkHabKVSuUtDRp473vPDXRI+7d/20+zn3V2dqKxsREANrluLo5jxWIR+Xze+a2hocERwOl0GgCwZs0arejF4WjlchmrV3va/he+8IVRzznWmJA2DMMYZ7S0tABwPbpF+ACe4OVUmyLc0uk0BgYGALglHsUrOpVKkeOXq92yV3alUtXjitbMsdCs1XNaUAD49rf3BODFMksbY7HYqNq0HKdSqahg7u3t1ZSinKNchH44HNZjt7S0qPD22hnR88t9rGdMSBuGYYwzRDBJGUnAracsHtSAKzDlbwCqDVerVRXc/mIbvM+wSTz8VlgVUC6H1IObt+d1cP+aOGu7knyFhahfqHNubmlnPp/XiQfvGzQx6OzsdK5ZhHuhUNgqb/exwoS0YdQh9957LwBv0JHBiGNEeeCVMBgAOOmkk3ZwSw1jy/nBD96Pr3/9Sf07mfw7AODBB+vX3DzWmJA2DMMYR1x//fXYfffdAXiaYX9/v37m8oxiwmbtOJVKaWIP+R8YTkoSpHHL98ceu0T/Fu2V16RZK/UmlCPTjl500V8dTV2Oc+yxP8d99x2l5wqKj67VhlOQsmlf1tNTqZRTbUu8uHt7e7U6FucuL5VKus3111+Pc889d8Q56wET0oZRB9xzzz1ad521ZC5IALjadFC2pgULFjhrh6effvr2bLYxBrDlpFqt6jptoVBwBLAIusmTJzvCTTyZWRCyaZq/Cyo9GQoNC91kMulkNpO16lqtrGZwLxmKt30kMpxlTBKbyDGDLEPc/2OxmK45c0w1C302jYvndn9/PwYHBwG4+b39+9YrJqQNY4y46667dFAtFAqOAOY6vjKAsOlb/gbcCkJMrVbDddddp8fhgZ2rEsnxzzrrrG16fcbE4t57jwQAR6M2tj/1v2puGIZhOMRiMcRiMU2pGQ6HUSqVkM/nkc/nEQ6HUSgUVLsWjZStL5ySU/7x77xPpVLRqlUirAG8lQ40hmg09lZs9nAqUOHKK9+r33upQb3t7777M852cnzPY3xke2u1mtNWTkHqb2+xWNTrD4fDel9kSUD2l/tYz5gmbRg7kDvvvNPJhCQmbhmAADcRhV+TFuQYAEZ4zrLJjzVv2ScSiTi1huW8119/vWNyFG0+EomY2byO4LVkLjpRKBS0kIbkyRbk2VcqFfVuFgGXSCQC16QjkYiaudmK4zeJy1J0tVrzrWND9+Xvw+Hhvj287fAacjgc1u35+vgY5XJZr0N+TyQSjtCWdpbLZb0vuVxOj8PXwfe03jAhbRjbiTvvvBOAJ4zZQ5sHzKBBsFKpOLGqMtMfbcCSQgWAN/AEZWXi7Tk8ZbRY2Gq16rRTihFEo1GceeaZ7+S2GBOAn/70swCAL37x55vdds6c5wNTfh5//P+n/e2ee45wJp7GMCakDcMwxhEXXHABfvnLXwLwtEe2ikjMdK1W08ldpVLR7FsDAwM60ZOJWiKRQHd3t34W/JWpRNuMRCL40pceHtEuMVfL+Yfzew9PFhcs+Ed85St/0H1E2IuJGoCThMU/0RTtee3atfqZy2yKQ1ljY6Nefz6f1/uSy+Wc+yXHuOCCC0ZcT71gQtowtiF33nmnmrA5daLAntv+2riiVbC2zeuEfscxHpxGyxLlP7Zsz05kQW3zH4dN6Fzi75xzzgk8n7F9EaHLQpqfaygUUgGUzWb1c7FY1BSZkkikWq2it7cXgOcJLn3J74wof8uxNsVoObRHq1FdLBZHrWXNFiZpZz6f13ZImtPGxkbdNpvNOnnBeSmJhbvcx3rGhLRhbANuv/12AN7gIRoBawAyUJXL5VHDq0ZDBh7ePigjk2zLgpzPyx7jPJjzhED25fAbXhcXpx45709+4uVnlgHWhLaxKc4883eqPZ900qNj3JrxgQlpwzCMcYYk54hGo44GKhOuaDSq8dPVahUbNnhlJvv7+5FMJgG4sdGijRaLRU0O4vd74DSiQXCcsred65kt28h57777MzpZLBaL+jkejzvn4KIa0k7O4y2Wgd7eXo2jnjRpkrY3HA47zpryfSaT0ftYz5iQNox3yC233KImtHK5PMJZi+OS2fmLBzWOV+Z1uEgk4jiLyXH8sdFBqUO5DawNV6tV1ar923P8tMDetgx7pMv28+bN02Obg9n248tf/jIA4KabblKhyn2rtbVV+2RfX5+uyXZ3d6ulR4QbRxP44/Dl+Q4NDWHt2rX6vRwjHo+PWs9Zylhy3+Y+HwqFNOHIwMCAbjN16lSdSPiFtQjYXC6n7Zc17NbWVhXcr732Gtrb2wF4SwKS8zubzTom9vEQtWBC2jDeBnfccYeuPbNgDoVCPueZkekVWTAzPEDK4Mie3rwPn8cfajPaoCl/j1aOcFNZl/zXIMdhEzrgDXyitcyfP1/3Ofvss0c9tjH+uOGGjwAAvvrVP41xSyY+JqQNwzDGKcViUbXacDis5uBVq1Zh8uTJAIAZM2Y4joKiVbNmLLDzWT6f17zgqVTKcR6UfX/0ow/ga1/7Xz2O3wIDAJdd9m71tP7e9/6uv69ZswYbN24E4E00ZcKXy+XUQaytrU339U9W5fiiMddqNTQ1NQHwKl9Jrehly5aptWG0vOD1jAlpw9gKbrvtNgCugxjHMrM2zA5fHPYigw2vwzFsmo5Go46WzrBjFw+OQd/7z8Me5n7Pb/l+NKc21sj98a9+87wwf/58Sz9qGG8DE9JbwV133QXAW4f5whc2XVrtvvvu00H8+OOP3+5tM7Y/t912mzqabEnhiyBzN3tWx2KxwFk9fyepH+UcXEdYvvdSMw6/ytwWDqfh9WZ/ukU/vA7t/z3I69t/jfIdr8XLcRYsWGCC+h2yYMEC/SzjDC91TJo0CdOnTwcwMmmOIE5WsVgMHR0dANwsYytXrnQmdxKylclk9Llu2LAB119/MADg7LN/j+99b28A3rM/44zhsCfhpps+quvDGzZscOpfy/Hz+bxaBDKZDHbddVc9jrSTq1kFJe8BoNefy+XUcY4z6cXjcb2P9dwfTUgbhmGMM6ZMmaKfxfScTqdVqE6bNs1x7pLvOVWsmIA7OzvVTByLxbBy5UoArsDO5XJoaWkB4Fp3wuEw1qxZA8AVktFoVBOjRKNR9VOoVqu6vd9iJG3MZDKO34OYrXfZZRdtJwA1ibOTpQjgarWq1z9t2jQ1q9dqNZ0MSAKXeseE9GaQxA3Tp0/HrFmzAHgd6vHHHwcA/OUvf9Ftr7vuOnzgAx8AAMycOVM75pIlS7BixQoAFkc6HrnlllsAeANVkCkXcLVL1nzlO9ZkeGDaXFgLp/Dk9UI5LjA82AqsvcqAFIvF9Lzs3MUDKMNt9nuAB6UUDTK38+dYLOYc54YbbtBjjQcPW2Pz3HjjoTjxxEcAAFdfvc8Yt2biYELaMAxjHHD//fcDAHp6ejBp0iQA3iRIQpBeeOEFZzmEP8tELpFI6PcygWMNtL+/f8SkTJBohubmZjUfM4ODg7jqqvfptqLJDg0NqVYtdZ2Fb3zjKQDAokWf0JAqhv05+vv71ektmUyqNUFM45wxj/09+Pqr1SpmzpwJAGhqatLt//CHP+DNN98EAHz+858PvP6xwoT0W0gxBM4PWy6XtSNMmTLFqQokGsK0adP04U6bNk3XSXgtccqUKaoxzZ8/3ynALh3phBNO2O7XaGw9Cxcu1MGJ17P82idrnewwJr9xnChrotJf/KkWWWNmc57Aa8w8OPlNjoLfOYyd0XitOuj7TaUvDbIosKmSr5XXqtlx7aabbsKbb76Jnp4eGMEsWrQIe+2114jvK5WKCttUKqXCKB6P63MqlUr6OZlMqlCXbTOZjMYl9/f3Bz53Nh+3tLQEWksGBgb0OMVi0enz8v3AwEBgCGAkEnEy9XHOANmehXQ+n1fTt/SlfD7vpOSVY8Tjca2CVSwW9X6xjwcA7X+LFi3CqaeeOqKNY4UJaQAHHXStmme4ilBvb6/OWKPRqLr6y2+ANxuV7eXhA9D1G8BLJiDhEOFwGJ2dnQDclIwPP/wwPvvZz26X6zO2HknzKS89MNJZSwRUuVx24pr9cJYlFpxs7vYXM+ABThgt73GpVAo0TXM5PxaSvC4ZiURUy+FygTzI8rWyxzib4rntvObJAzInd+FYbvEkX716NRYtWgQAdTVIGtuPH/5w/7FuQt1jQvotLr30OQDA97+/7xi3xDAMw2X69Ok6sWlubtbJV6FQ0O+TyaQz+WHrimTiikQiqlWyhUhMzaOF7fHEyp9KkyeIfmuP/zOH/PHv7CzGE0ppv7RH2hkKhbT9Yo0qlUp6nf6JLlsS+Dy8DCCmePEKrxd2aiG9ePFiAG7+26amJkdrkA6QSCS0c6fTadUWuMIR1w2u1WqqTXN5tK6uLjU3cUGDcDis7TnmmGO20xUbW4qsczFcio/X/DjXMGcCE/i7arXqpHEUeAklGo1qv+Nzynm5PX7YYzYSieiANDQ0FFj5CsCoVZT4PNJm/zG5/f5j8+Dovw42Z/K5RJO/9dZbccopp4y4vp2NG2+8EQDwnve8R60wzc3N+pnHnMbGRnVQzefzTngTO/DJ99LHE4mECjfu17ykw0K/WCw6XtT8Tkhb2CIZjUadcY4tNHIcMXXLuRg+l7QzlUrpRCXIORKAM0b39fUBcAVwOBxW83kqldJjJxIJve9nnHEGxpqdSkj/4he/0I7Z39+PXXbZBQDwox/tqg/rvPP+jLlzPwjAm3VJB2lpadHOy+sdXH81kUg460CyfUtLizMQyjFTqZR2zkqlopOAm266SWMYGxoa8JnPfGZb3wpjE9xyyy06QPFLzxMyfyUpFp7+AgWcHckvuBjpS4lEQidy1WpV28ICjZON+KtecRwoa1NBa898XTzg8sSCNR5e6+RyiEEpTTlud7R1bf93/I7deuutAGDCug654IInAJjlcUewUwlpwzCM8YQ4M9VqNSdPujhNDQ0NqeLBDl2JRMLxYRDtWSZ/ABytV5QFfwY5nlDKZC2ZTOLmmz8GwDO3ByX4aW9v18lfe3u7Wl1Ye25qatKJaalU0m3Yocsf2iftbG5uHmHtSqVSag73K1Ly+Uc/ehM/+EG3Xgc7n7FWX09OjBNaSP/iF78AMNzRp02bplrJa6+9pporx5kuWvQJ1SCKxaJ6JTY1NTnONvLQGxsbAz+zBy2b0IeGhvR84XDYSTgg7Wlra8Nuu+0GwHtpnnzySQBQL3LTrLcv+XzeGcDk5WUNtFarOU5fsk0sFhthVi6VSo6pkM2QrIFzdjAhFos555ffOKEEVwfizE3+En1BZmV2epN95PugjGqlUsmJYPBX3/J7zLJmzl7fXAmM3ye+v2b6hmbY4j5Wq9X0GTQ0NKiZlvsqVzpjC1AsFnOeg2zLsfUCf+7q6lJrY7FYVIEajUZV8LO5uVKp6PeVSkXby23kBCatra3o6uoC4FWqkuQj/vYEtVmuJxaL6bXyOd1a6sPHbGhocAQzv8Ny3+uBCSukf/WrX6lwZvOcPNhUKjVi1gh4HVdmV+wIwR6uqVTKifdjV382ZfPaixwzlUrpMeXlEnhffmnkOHI9v/rVr/DJT37y7dwWYxPcdNNN+jlIcCUSCUcYsf+BDA5B1aZ4kGTvZ/aFYLM2z/w5HIWPk0wmnf4lxxwaGnIEuUwI2cEokUjoufL5vCPgeSIaJKQLhYKTCMXvtc6TlEqlErhswKZvGbT5WmQb9p43jG3BFVfMHusmbDUTVkgbhmGMZ26//XbsvvvuALzJj0z0eXLn115lArR06VLVBru7ux1lQxBNs7e317H6yaSoublZtVtgOBFJNpsNzLwXjUYdxUO2Z6cwNrez01upVFJNva2tTX1y1q9fryZsPn5vby9aW1ud+8XXVi6XsWrVKgDAxo0bnXvETnQ8oZX7y6GLt99+O0466SSMJRNOSP/yl78E4HVMeWgyyw+FQqo1sPmNTUmsDfMM3++FGuT8w1oPbwPAMfcIfC52BPLH1XKJOMDLtSvX+elPf3or75ARxG233eYMJv64ZWBkDmTZPp/Paz/I5/OOVin/y+DCyXJ4Da9cLjtpO9mbmk3Asi9bW2Q7+Z89adlhTfZlE3oqlXJM0uwIyfWyeY2QTeVB18pmVl5nZIKqf8n94f/lcz152+4o2tranOcXFIPP/YP7YTKZ1MIUo5l+ue/x2rAIq66uLu3j6XTacSSU48l5ATeKxb+EwX2Sk68I7LldKBQ0MoYnCZlMxkn+4y9kUy6XnagIuf5sNqv9MJ/P6xKmP6kKL8HIfZfJwlgyoYT0fffdhxkzZgBwByh+mPIgEomEY+aTB8edj6nValq9JZ1O677FYjEwFd2GDRu0o7W3twd6s8ZiMT1vqVTSY7JZlTsSD2wyS77vvvs2W5HLGB2pbOaP/RRYKHHf4PC8arWqgwAPSLz2LLDpltdv+di1Ws1x5GGtSb4vFArOxI/X5zhWlL/nNsqgzANSPp/X9qVSKR1Mh4aG9Hse8EulkraHKxKNls2MB0EeWDkGlyckLKhl3xtvvHGnEtSGMaGEtGEYxkSBfRD8ubAlEyJbdzi8s6OjQyc5iUTCmUSJ+ZgrXAktLS2qXGzYsEEndo2NjaOWXmXriygP4vgFeJNIybjIebz5WP5CMuvWrQPgWQSkPYVCAel0Wtss38v1xONxZ7Ir96ujo0OPx7kA4vG4Tv42btzo5ADnUMCxZkIJaXa44ZSbYubjtI7JZFJTe4pXNRDsNQh46yuSVL5cLuvDW7VqlW63cuVKTJs2TbdhbYjXT4JM5eVyWQPuOzs79Zi5XM552QDXlMQd3th65KX3L2cIbGZks2+pVNL+xbHD7BUd9J3fE5y1Ybb8sNNUkEOh31tbYOcvf8Utzl/M24tpvVgsOkUOxMrjX7YJcnyT94HzA7CzpbRfYJM5L/UIQdcGeO+MOPhN9OpZ6XRa32+OHGHBMWXKFDVPcz5rYHj8qlQqjolbnn9QUpFoNIrnnvOyL7a2tqog9CctCcrmVavVdAzLZrN6/mw262T8CtqX+wpbegYGBrB8+XIAwNSpU50wNH4vAddz3d+n5L4kEglnbV/KZg4MDDhWHHkPZHwYSyaEkH700UcBeA9Cbm4mk9GHKIMNxxFyYvpcLqdOCzwockIHnmkVi0U9Ng94xWJRt4lEIjr4bdy4UQdITrvHgjyXyzmZc2TfWq2mHUVMsiykm5qa9PoPP/zwt30Pd0YWLVqkA14ymXSEJD9XgdcE/Q470k+Clkv8a7fS19gTmp1YgGHtxj95kOPzIMSFWiKRiA5a/kpIPJjyGrZ8bm1tdUz4nGiH83EHTVS47SxEpI3xeFzfQw495PKfvL1fsHDbR6vSZBgTkQkhpA3DMCYKCxYsAODF8cpkbWhoSCcqbW1t6jnd3NysBXvS6bRmUVy+fLmGbJbLZSd3NvvlyHnk98HBQZ1MsfWvXC6r4sAWFGB4UlYoFEaE0sln/p6dDAUO1ysWi87EUSyhg4OD+j23mbPeBTk1ZjIZ9VVKJBJ6v7LZrN7HUqmE/v5+AK4lKxKJ6PM466yzMBaMeyF9//33a2dqbGzUB5fL5XRdhJM78Gw+qOPKb4DrfVssFnX2HxQLK595xi/bs3MZH1/ODbhVtpLJpL5AnHNZHNcmT57slFsTTfv++++vu1qo9ch//Md/APD6CGuaQYlFWKMDXNOafF8oFJylDU7pKf/LgJRMJrUPcmy0P3d3UAy//AZ4pmOOzQ6qgjU0NISrr34dX//6NOe3SqWi7fVX9pJ+x4NcoVBwju8fXAE3FwF7JAdZJPznZPO8P8WpbCPwPbvhhhtw5plnBt6n8Yz0j7a2NqdCmTgEcgTBxo0bNR91R0cH1q9fD8CzsMlnTkvMy3kidPl5xmIxHYfYksgx9/487EEe+f5j8rZB23N7/KZ1Wc9m6yhbG6U/VKtVJ5mJmN6bmpq0DR0dHfp5zZo1erx8Pq/9jHNcJJPJwBz5O5JxK6QXLlwIwLvp0pG5niqHT8lMq1qtqnBramrS74vFonOMoCQObK7jgURMl/KZ14j5ZQhKEsEetGweTyQSjlOHtLO720tn19zc7NRtld8jkYjely9/+ctbeCd3PuTlLZfLeh+BYY2AByReT2bP7Wg06gwO/Hw5y5fAHtrsbxBUHYjDsXjdjr2iObyJC12ImfrGG4dDZlatWuWY2f1hVhyWyAO3DFrJZNJpZ5D39miTVvbvEK2Qc4rzvqMJb36Xg7KmGcZEZtwKacMwjInGwoULHf+VIM9tntDl83mdpHV2dqrJtqurCy+//DIAOM5kwLBiwb437OPCvggyiWXnSM4HwVZFzqzIGjvHe3Pb+TP7ScTjcb0mnsg2NTWpeZpTo8r18KS4Wq2qJWGvvfbSbTs7O52YaU5swlbWoLjuhQsXjonyM+6FdDKZdMIL5EazuXDt2rUAXEet5uZmvfnA6MkV5CF2dnaqlsF1S3mb9vZ2NXWmUiltV3Nzc6AmAAyHLPirEklnZE9y2bapqckxRbIp1dg0Dz74YKADnr8CFMcRc91d1qo5WY7f7At4g9DNN2dw6qnDXs2cUIdh06Lf7CvfZ7NZ7b/+sqjcB6ZMmYLvfGcKdtllF7W+sAnPf145Zl9fH15//XUAXqSCHL+5udnJzMRx0pzwBPAGTLZCsGAJCvmRCmFiJmdP5SDLAlOtVjFv3jwAwDnnnDPi9/FIW1ub3jO24iQSiRFJSAA3C1cymdTnzWvS7L0t5wCgns0siKLRKF544QUAwG677eYsXXD/FPi94YnEwMCA9oN8Pu+sbwcl2/HH1/N5ly1bBgB43/vepwI5l8vpezxlyhQAXnYyTqTDRTLkvoRCISfXOefK4ERFvM7NUTVjwbgU0gsWLAgUWLz2l0gk9CGyRyzHCEp4AYcF+DsgrxNKp+MBdLQi5xzawuY6HtwB19QnA146ndaO5F8DBLyBkNd1uKqMzMIXLFgwZo4O9czg4KCzhMFZvkYr/RgUeuJf22MzLYfKnHlm64gEJjx752IcQWbyyZMnq6MLC0AO2YrFYjqQ8GDMHuD+er1BQq+9vV0Hs4MOOkhrEy9dulQnjclkUgd5Trojvw8NDel7xV7cPPByghi5fwKvwfrjeQX/O2oYE5lxKaQNwzAmIi0tLTrxyeVyjt8Dh81xSJ9YJhobG9USweksa7Wa44sjk0C2frAVTtJp8uSpVCoFFv3h47C/DTtQcmpj9tngSbJcI+BOkuPxuLaHHX8BOGUp5dpEAy4UCk5VQbkvHCbL5/Jrz0FV5djyuiMZl0K6u7tbtQvuSJxxhsMEODCezSHSWXp6egI9ujlXcjqddtz1gyr0cLKLwcFBxxQo2gU7pnFnr1ar6o2ey+V0m4GBAce5TK6HOzp/FtPPaOlNd1Yk13kmkwmMHQbcJDOc+J8JyrXOWu1oplmB8wuz5YWXNVpbW/H+978fALSfAyNTZfL3Av/u115HW3Lh+8HmfAnn2W233fDaa68BAJ544gktmdrd3e2ktZXzjFY6MSimmk37HBHBqUL918jvp7T9+uuvx7nnnjvi3owXfvaznwEAJk2a5HjeB1GpVJzlMLl/q1ev1lCj7u5ux4IoQruxsVHXakWItre3O86Mks2M075yv/WnuWULoxyHBRp/ZudEf054fue48hwLWM77LX1dolu6urqcpEHSzzo6OvT71atXa1/N5/PaV/33mt8Jrusgz+lzn/scdhTjSkg/9NBDALzMM0FhI62trc5ahnRMFrTiWMFrdH5HCCGbzWolFa4JzPF40uEBb+1bEsLzOmWhUNDwqWnTpukMjwddTjA/NDSknzn5vVwPD7681lMoFBxBLvfrqKOO2sydnfjwvQ1abohEIs7gyBEBIojYTMsaQa1Wc9YR/X2JBQ6vgzU2NjrCTZ7lwQcfrGtvPHj5wwZZWHHtam77aNnHuJ8E1bdmQV4sFjXMZ8aMGXj66acBAM8995zT3+U8nOxEjs2Taf9EgrVF0XJ4cOQQMPnN30ZLcGJMVMaVkDYMw5iIiFbLzonsJ+EvMSl/cwhmS0uLWuMmTZqkGjEAJxTOX4WKa4TzhJLPUygUnBBF+cyZ99jczb43fkcw9sORyXMsFnMmwOysxceUSVw0GlUrkxyDr43DDMPhsKZ05vzfoVDIyWnAE+agaom8/Y5kXAlpCWqPRCKqabDjWKlU0gfQ0NCg20hwPj8ITkbhdw5i85F0gL6+PifNqOzrdxyTjjY0NKTm7kgkom1YvXo19txzT/2ePRp5nYnru8q+cj2FQsHprOJEVy6X1Vkum83q/TKG16/YEW80sy9XpwLcalbs0BdU/Yz7ksBmPTHZzZvXi0gkg3PO8QaaaDSKgw8+GIBryvMnjmDk/P39/c46m2itXAqVLTVc7YrvQzKZ1D4myzPyO1udPvCBDwAAdt99d/z617922pJKpfQYrMnzNmzy9N9rrqoVlKLUn3SGLQXj2dObPeZF6A4NDTnx90EJbjinAiep8QtDgcugcqIS9tyW83MOiIaGBqckJQvmIK9zfypXhvsSO0sGlbCMRCL67rKw5zZz7nleSuKxlcPKZHte3mHK5bKT3pnzUIzFMuK4EdJ33HEH3vOe9wDwHiKbfuWG5nI5fRgspDiTF6+XyUDEszU2LxaLRcerlTOXSUfjGVtXV5cTDiZm7Vwup8dhM7jfvDh16lQAnlmeSwZKp5UJxvr16/VzPB53HBvkvvDAfscdd+DEE0/cwjs98fjlL3/pTJjYE5q9u4XGxkZn7VngfTmBCA+OgPeS33TToD6Dk0/2vKBvv72IWi2Hr341rr4P0u/23nvvzTqm8MD35ptvaiGEdevWaTtTqZT2+56eHtWmGhoatD9Wq1Vd9slms7pvJpPRSV5vb6+a3+X9AdyCI83NzVom9fe//z0AYNmyZU7qRq4nLftVq1XHSz0oKUkoFHLKx4ZCIcyf34dqtYozzmjR+zFa2JphTBTGjZA2DMOYqLDntjhKNTY26uSeFRC2uvmztLFmGJTNjkNJBZ78DQ4OOrXDWZuWCSqbvtnJCwiuB+73q2D/CU7zGVTRL51OO5nnZCLr9wwH3OUBTo7Ck2G+X2yl4hwa5XJZTemc44KtrDuScSOk+aGw2RGA42TFaytippCH39zc7IQmCH5HFu7oQQlC/E44/Jm3Efy5mNnsx99JO8UZB3CTQHB4AZvn+V4IrO37X8qdDQ678KeyZO2OvTgFdoTyOxiy49i///tw/dzTTmt663vv70WL8o7j2IIF/TjrLC8sRI6xyy67BPYLXk6pVCr405/+BAB49dVXAzVQXnNbuXKlJnTYc889deBJJBKOhi1pUrlfFYtFLeU6MDCg2jR7m/Og9ZGPfAQAMH36dPz2t78F4HnVyvvGiVgKhYL2dak+FvQOBSW3qNVquOGGAd2GE8UEmS7HCyyk2cwv96+lpUWtHENDQ4HhWGwR9I8t7MzHFdaAkYlueDwUiwsAZ8mG7zVXXeNjBOUUGC3vd7FY1OvzOwkGLfcUCgUd93ic43sn8JjOEQTsGc9JSzo6OpxMb/z+jwXjRkgXCgVnDZBvmD+JCOCuafFMkxmtQ7Pnq8wY+cXwF1fghygPmkO8OEd0IpFw2us/t/87zorFafCkXZwfmYU1ZysS4b6zIWFXuVzOWR7gUCR+YYM8tNmr368FsJD8l38ZDpOrVqv48pcbkcvlcMcdJZx0Uhy77LIL5s59U/edP78PkUgal1yym+4ncB/kNdvHHnsMr7zyirad+7qYshsaGpyJhDz7l19+WeNNp06d6uSJl2WfbDar5mnWcsrlskYx5PN5nUTyueRezJw5E//8z/8MAFiyZIne987OTkdrkXbFYjGkUikV+OwXIG0Jh8P6HIJyfstnS2xiTETGjZA2DMOYqHBqTTF3x2IxtYpwJje/aXa0DIhBWQnXr1/vaMeA55wqEy+e9L/yyiuaOvSDH/ygk6mPJ0vSltG0Tv/nIGc49i2KRCJ44oknAHgpP9nRViZ369evV81baGtrw7Rp0/Q8bG30W6TknJybnK048j1XEEun02MyERw3Qjqfz2vMcmtrq3NDg9ZbONXnppJLAK75GoAz4xf4QbNJimuoAnDMLdxJhHK5PGpCitHa6b+OSCTiJIlnsyufR2Kzg8zhOwPyMrKTF5u4/SVK2VmMHcGC1rb4WfGSBwDHG/SUUyKYPXs2IpEIvve9Trz66qsAhq0qQVmW+G82Ff797393zHBcopVr8fK1CkNDQ5qQpFKpaNILdoYrFApOTDiX/ZPjptNppyqbtIHPJd8dc8wxeOCBB3DTTYOIRPowZ84eek2lUgnz5vUiGu3HN785Y0QWK342Yp04/fRmFAoFLFo0HHKzaFEeX/5y4wgte7whHsyNjY0qfFpaWjQ50eDgoApsYKTlBXALXPByDGf82rhxoz5vEfoNDQ0jagcAniOhLHk88cQT2H333QF448nMmTN1WzZry77cD7PZrI7LbJLnOPrXX39dHWSXLVum5+UMafF4XNvZ3t6u1yH3ZePGjep8yw6JsVhM70s6ndbPfkuMtKWlpUUtU4ODg/o8SqWSPqcdybgR0pxZK5vNOo4TbPYT3omrvH89GXBLSXL6OX+GMoHXkPwD0DsdTLi2NTtx+B0hOBPZzggnkwnKJsbPIR6PO+t2nMZwtOQnQQKbMzPtt99+ALyBVxKA5HI5x2tZ+lF/f7866fijDSTDF2d3amhocAZC6QOpVEpD8tjpZ2hoSGNoX3rpJT2+DLaAO2ngMDD2i+B7s3r1at1X1k55splIJPC5z30ODz/8f7Fy5SpcddVruo2/hrTsx/c66L0Kh8P48pcb9f7xgG/e3cZEZNwIacMwjIkKe1SzeZWdmUS7KxQKOkFjJ65UKuXEt3NMsWiAxWJRtxHzcrFY1AnOwMCAkxzlve99LwBvErl06VIAnpLiWYSAI444wjE7s1VPNFx/pj5pV1NTE37+858DcC0InZ2det5wOIw33ngDgDfpEysN56S45JJn9ToWL95Dj8fWAblmrmDI8frxeFzvdSwW0/uVTqcdi9JY1DAfN0J6+vTp2kkB14TJWnVQkpHNmbv9wfyy7yWXPIu5cz8IwHtwX/3qn/R48pJceeV79Tjf+tbfdHs2yXNe5tEcxyKRyKiagL/9vJbD5lDuQOyVLmEcOxv8TDfnlMdOVvwcUqmUYwZnT1belxMrSLIP2e+3v/1tYI5gNpM/++yz+OAHPxjYrscee0yPLc+aS5QCw7HMkyZN0mfP1xqPx9XM19/fj+effx6AF/EgmncikQgcwDiZCCegCIVC6si2zz77APD6HZtfU6kUIpEIOjs7dPmFza/VahXXXLMKfvgdSSaT+mx4ucjvnSvMnTsXF1xwwYhj1jOcHYyvSd7dVCqlzzWZTAYumfF99SPHqVQqI/p5Pp/X4zU0NGiYU61W0z4wdepUFbrPPvusHmPp0qV417veBWBkHQN+ThxRIe/F0qVL9TihUAj77rsvALfISLFYdPqxCHK2oPJ7LtfJViTev6urK1DQ8nvJpYK5XyWTyTEpB1z3Qvree+8F4GU34sGHPwetz/DaYxD+3zi8Sx7YDTd8xMnN/f3ve53ooouW6ncXX7wUK1Z4Ha2treYMPtKpi8UiLrrorwCAW275uJO8ndszmse6v82hUMgJYeHfOdSHTbZyH4899thR78lEg9d1BV62AIbXzvjehUIhfcljsZizzsbmVfawlwnk/vvvr445Ei7FnvZ8TqZcLuO//uu/AHiDlGyzfPnywH3ZNDwwMKDfNzU1BRZC4HSPqVRKS1K+/vrrKuBZMOdyucAJL9/LcrmsQl3Wu2fPnh24tNPY2Kj3xS9HarWavkPXXbceZ53V5pwzGo06z4kdooLOJb8fccQ9AICf//yLMIzxSt0LacMwjInMnDlz8NnPfhbA6PXm/d8HbcP+AP4KUzJBymQyTvEW2Z9TK8skr7293dFIOTveCy+8AMDzbxDnNk6XzKUiAdcpU0zJL730kh7/Pe95j27PsfjZbFbbXiqVHAuifxJeLpd12+7ubiecjx2NOStkEDz5Y293VuDmzJmDK664InD/bU3dC2kxxdVqNSfWmJ1aOH46yPzFSRHYU5o/ByWs4ByyF174tOOgFmRC/9739saFF3oOQjfc8BFccIEXRnDFFbNVC/cnDuDz8kslBHmt8rWFw2HVftiqwJWc8vm83sedic3l5eVtKpWK3jsebAA4mrT/3gPeYCam6ldeeQV//atnNeH+xVqemNsGBgbUzMnVzvr6+lRDTaVSTlUrdqYKyvHsv35/X5J2i2ba2tqqWvAee+yhfWlwcNBZo+PiA0HOXayN8/3hdy+opGGQUJJkJV/5Sgf+/d83IhwOayrQarWKW28dGnEu73hyv0v4+tfHjwbtjzCQ+8RjHmcT85dYDTJ9s2Cq1WpObWeBrSbS91taWpzKZlyClfuDVEVbtmwZ/vjHPwLwEtpIuxobG/WYvD5cq9V0ez4OZxzz18uWZZr169erdbJSqYxIoxsKhfU6ud/7hTTfR/6Oo3GCUkBzfo4d6aRY90KaiyHIQ+eEFBwewmbEaDQ6Ym0uSDgKQR2d17hDobBPMI885qWXPodIJDg07OKLPRP5PffMHPW8Qe3kz3I9XIVG/gbcNUt21uDwgp2JoBA4FpIsxNgMzsKF12M53CQcDqtH84c+9CE8/vjjAIAVK1aodiACmEObVq5cqS94NpvVASiTyThJdILyeJdKJfXF6O3tdcztItRXr16tpuxUKuW8P3wfxNNbqgMBnt8HD2YyCeDjFAoF3aa9vV3vAe8XDofxxS96DkF33/0ZPT8niDnttCYsXJgd4T0v1Go1zJvnheFUq1XN6FatVnHKKUncemvBlx3QnTQfccQ940ZIG9uGiy6aMdZN2C7UvZA2DMOYyPgzKLI/RFB4KRBsBmdLC+Cas0U77urqchKn8LH8sMmcHc7EIRDwFAcxMb/wwgsap9zb26tOXq+88oqasFevXq3HaWtr0+Ow5c9vBZUJYE9Pj6a55etraHgJAHDggQc692e0e8EWyyDrZSwW02eQSqWcrJOcxnRHMW6ENJsjhoaGnNhg7siisTQ3NztmYDkGO1YFreX4zZligvnRjw5QUza/UGxm4wd91ln/rZ7enO7QrzWw5huUF5o1ELmezs5Ox/tQ7gWnKOVyl+M5ycM7ISj/LuA+P3/qWPl95DN216Q6OjpwyCGHAAB+9atfqQm9ra0Nu+yyC4DhpB4bN25UjTOfz6v2Wq1WdSCLxWIaD+15Q3fq8biqlAxYiUQCy5cvH9GutWvXai7u7u5ubQv3By7pygNxX1+favaJREIHorVr1zoJcd797ncDcD2Sg/MAhNDS0oJoNIpcLocTToji5pszan6NRCTVbuitfyMH1aC/b7utgJFyhfcdP6ZuALjmmmvUEsMJQRKJRKCZlr2lR8ttzXnVuf/n83knvh3wxhW2InHlMum3sVgs0MTLntIrVqzQhFMc0fLyyy87a8lB+8r5AE8Ayr79/f3OMhQLSemf0jcLhYJj7uaiHgxb2Hh51J+gyH8d3N5rrrlmxL3YXtS9kOZOKgMhe9aykG5padEZUCKRGDVXt//v0VK9cUH1f/3Xv2DYxD28zcj8wcMdWTr7RRf9FT/60Qf0mOyMEbTG6W8jm1jl2rgTy4Drn7Dw7DuoFu1ER8oo3nTTTSOS9gNueBWH87BnfKlUcvaVtdcPf/jDmhu8qalJBdeaNWt0oJDBo7Gx0XHoEQHMJmu/l/66desAeOvW0sbOzk7VVDj/dX9/v7P8Id+vWrVKJ3NtbW16HevWrdO2cRKX/v5+JyRGhH1fX58zsD355JMAPM1m1qxZAIaLwjQ1NeGII+7Rdezf/e53OP/8bjz//POYP38AJ5wQRXNz21shioMjErcEaTz+CVPQpNufsc/M3TsPV131vrFuwnZl5xu5DcMw6gyO75VJd6lUcjyk2ZtZJuDsiOrXBDm/NztIcsy5/C8KBYeu+kNExaqYy+V036GhIZ24Shy8f1//Z7EwcVlMVl5YcSmVSk71P1as5Hv2HGdrA98LhreRySqX3GSr2uDgoG5TrVYtLWgQcoMikYjeIH/WGOk8XG0qCDYdcxIQwPWylm1OO+23KBSGk6aIluxP8xlkpqvVamoe//7399WO4s/zHHReNo+PZqqW62SNijsxO/vUarWdNn834N1zzqnOjnZcfo9zswdldJo0aRI+/OEPAwAefvhh9Zjv6enR7ePxuA54Mugkk0ksW7YMADBr1iw9f39/vzPYBWmOqVQKBxxwAADPfC45jQFoJibANclz/xJtuL+/39FS5X6wM9zMmTNVI3/ttde0uAKXWg2FQtrmvr4+NdG/9JK3Lrhy5Up88pNDuPFGT9Cce+4bOPHEmCbGWLwYiMVKmDMnpVq0Nyi7yWb85vOg5YeRSWpkn9q406Jl6WKvvfZyhCgnKmJLonw/depUJ/mIP4EN4PbhWCw2IvlRJpPR3+PxuI4tnNedQ57K5bI6Hk6aNEnzwGezWU2SwyZ2yb0OeHH08l4sX75cjzN58mTHPM/LUNI2f01t/1IVp4mtVqsjljvlHgls+q9UKno8XlYaHBx0JkfynHYkdS+kZZDhUoIsjJqamrRTdXV1BSYBEfwOEuzZK3CoQygUnLmM1yw5lzFngWJvdM7iw9mkeEDlzripNgty7K6uLhXAlUpF1ziLxaIzedlZs44Brue2f2kjqFIOl0kEvAEEAA4++GA88MAD+p2Ej6xdu1bXflmbkH6Zy+X0RW9qatLBlrUm/8RP1ns/8IEPqOB86qmndPvp06fj05/+NICRpVNFcL744otOsQ/WUFjTYG1CnIAGBwd1QGxoaNC+mc/nnfrsTz31FADg5pszOOecTgwNDWH+/D69p7FYDPfcA0SjMcRi4beuFbjyymV6r4cdnkK0dlnTSfGmKg+NJtQNY6JQ90LaMAxjoiOx6jNnznTipDkWPsj7OZ1OO6Gmor2ydzNrslySkc3QMomv1Wo6CWMlQpz/BJnMTZ8+XSe3TU1N2GMPL3e2TBQBzxokXtlNTU163s7OTqxYsQKAN+mV47e3tzuWIVaCZOLY3NzseIbLtfE1B90L9oBPp9OOhZGLvrClkvNwyHPakYRqm5h2HnrooTuwKcFIWk6OjR4tzaM/Pah/Bs6dblP5vJ9/fgMdQ0xrwbP5Wg3IZlcCABobdwnwOgVmzWpygvM3pe1z+7fkOgDXKzLIEa1Wq2nnFgefnYm1a9cGmkiB4H7AL3IsFtM1tPXr1zsWHO6PvC7GAyLgxhZzYQAe9Pwe5+Jlzfmbw+GwU0BB8IeS8LFkEPL3jaDrr9Vqugaaz+cdqxAPYHxda9YMx5N724YRDgflE+D2AmLeFhNtLueZHRsapuo2w89p6zTjWq2GAw+cuVX7CPvttx+uvfbat7XvtuCBBx7QZ5xMJvV5ZLNZx1mWk8hIXDxnx/InR2KTrXhgS/8rFArOmjQ738p6cywW0+WPfD6vyxt77723827xBODcc88FAFx//fVOVjAeg//2t78BAN71rnep0J0yZYpea19fn05OstmsniuZTOpkQqyU06ZNc0zsQWFqXAKVIyF4gsMhWIODg3pv0uk0jj766FGf3fai7jXpoGxhfjNxUAIPngFxNRjeVuBjP//8Bsye7Zkan3tuPQ0ubtzd8L7Bn5nXXstir73a9FxB+IWHwEUd+F7wGjevpfLxNpeVaGehu7tbBxj//Q+a7LAgam5u1nWzeDzuOOOIEOMCLRzWESTE2YeAM5tlMhk9RjKZdNaq5TjNzc2B/cdvKucBSQawRCLhVD0KgoUx40/6IoMWh/8May3hAIHsz9ftnyQFvRM12i6kx2QBH8SUKeExqVRkGNuLTQppqb4zltx1110AvJSFMijybKylpUXXDKvVqs6MBgYG9GWVGVdbW5tTFF1gE9AXvvAgFi8+CgDw+c//R2D1H78wXbr0WgDA3nuf53zPpqIbbzwUgGfikfbwwFoul53BU37LZDJqWmJvRtHu2tvbdYBct26dUx6O8/FKabnjjz8eOyO33347AHciBwzHSVcqFRVilUpF8xHvv//++M1vfgPAE5Jr164F4M3axSoRi8VUcMViMRXq4lDz0ksvaaWg/v5+NZlNmjRJnc+eeuopPcZBBx2kwm/58uXqZ3DQQQc5pkjuh7yeLoI2nU6rdpJOp/Hyyy8DgF4DMDK5w+67767bsNVF2tDQ0KDvWDwe10xgwyX/Eo4mHTTxZE1QJtwvvjgfALDnnmcAcKso8SSELQU8ERWT5imnJLFx40ZcddVVI847Hujv79cQPXZs9Fd9kufNqYt5glUsFrWfNzQ0OIlL/CGC1WrVSfnJZSWDqrexQ1sqlVJrZ2Njo7aBna+SyaRag0KhkDoAsw8Rv5OVSkXPWy6X9bfGxkbHmuBP38wm7mq1qu8zp1r1V7gLCl/lex2Px3UfGYd3NHWvSUu86KxZs/RhlUolvekNDQ16E3t7e/VG5vP5Eflqw+Gw/s6xn2wyWrz4KOfF50E8HBYtvEom5uH41Eqlinh82OMwFJKHHnLSJ/KAE2R6CYfDambiCYQIaQ5dqFQq+lI3NDQ44RtCtVrV+7izwlYGHrBEWHD+62nTpmH27NkAvDKTe+21FwBvssfpPaXPcFgJD0JyvLa2Nud5sLe4bLPLLrvoRIrDTpYtW4bDDjsMwMgkFkFVqsLhsD7rp556Sid+hxxyCD70oQ8BAB555BFHG+bJovQfOZ+0V65b7tNttxUA5EdYcUKh4LS2vAzgav4h0rJdaxKvKXLO5KCYd3mWixblUSxGME5lNFatWoU999wTgHffOVUsX7dYSLg0aKVScSrvcWpbMaGHQqERETCVSkWVl2g0qgI7n887isFoS0MyzsbjcfWW5uWV3t5ebdfUqVOddjFy/EKhoOctlUpOXgm5H1zDOqg0ayaTcZaJWF5wchS5bs5fzuWEeSlJlgl2NCOzZxiGYRiGURfUvSYtxduXLFmiMypOYxeLxXRWNTg4qJ/ZjCczrUKhoFoGp9DkSkOpVEpnposXH6XaxDHHPIRqdThV57AG4TqxueklJe4vhpNP/r96TJ5Bsvcmm59EA5o8ebK2n3+X6+RYxmQyqfdlaGhIZ5LZbFbv486KmPnvvfdeRyvj8CPJnrXbbrvhP//zPwF4yyxcCUi051WrVunMulQq6ffTpk3T2bloz01NTfqcM5mMmtI5X3FHR4ejJTz77LMA3NSb/qQMQfWyN27cqHWsOfdzJpNRJ6QpU6bg9ddf1335fnCMrvQxrmjkJnOoIR4frlwkxxgtHaNYosJh13zN2vNo5vFKpYKTT06gWKzirRUwXc8+5RRP85H+/p3v7JgSgtuD73znO/jtb38LwDX58/JGPB7XsZC9rnt7ewOTeHD99FAoNCJ+mE29kUhENWNexuHxksNFBwcH1YqUzWZ1iWfjxo2OLw3HUsv7MTAwEOiNXSgUdN+WlhZHs2ct3J8umTVpdobz3xOxPKZSKb2PbFZnXwt+Bt/5zncwFtS9kBZWrFihZkc2X/BAywNEUDnJjRs3OknixZSdTCb1GMcd9ws88MDndJthJ5yqb8042Agh5m/2Bg+HI5Dx9dhjlzgpJxcu/H+0DZxgQtq5ceNG7VRshmdzt7Q9FosFVnKSMAcDOPbYY7F48WIA3n2WAWDPPfdUk+Bjjz2myw0NDQ2OmVY+d3d3a3/YZ599NMnB8uXLseuuuwKAJh5JJpPaL0ulkrOWLYKc+2s2m9Vj77rrrtpGTlfKA2soFNJzPf74444JWvrJ+vXrdXLQ3NysA2tfX58T8sMCW9qWz+f1MyfScN8vN3WtfjtKNIUrmIf/9zuunXZakzrbyW8nnui9J4lEArfeOoR4PI5qtYq7767h2GPHv4OkTP6mT5+u4wAAx0wr5u5CoaBrwoVCQZ8Te2nL34DrFMlVBXkdmkOqhMHBQd2PnSbj8bg6ZUYiEd2nq6tL+y0X9cjlcjqZam1tdfyM5Pnm83mdUDY1NTmCWdrJFf/8kxBB9svn8060htyv7u5u5z4GrYOzN/xYYeZuwzAMw6hTxo0mfc4552DJkiUAvBl0UBUoNp1x9RLRMth5J5lM6kwzl8vpbHDOnD1w9NE/AwD8x398Hl/60sNvnTOOcnm4GsqwU02VNOOiatiRyLBjTyIRd2Z40o45c/bAccf9AoCnJV1//cEAPE2KUwBK+9kEw7GJXLWGPZRlZnjOOeds8X3eGTjmmGNGfPeb3/wGf/3rXwF4Ti9iZq5UKk7aQbmnHR0d6ui3fv16TJs2DYD37NhaA3jPVhwWE4mEYzHhkC4xM65atcopLRiUlMEfGSDpGNmrn7XtN998U5165syZo/s+8MADo8Z+/uAHP9DjyzVVq1Xceefwss8ww/2R2yuMpkmLKdv7u+Zo4tJ+dpJjotEovvIVzzs+m83in/95CDuwguB2Y+lSr/b8tGnTnKgBztctpNNpp5KTP9824PUzXoLg8FU5tixtcK36zs5Ofe6JREKdv3K5nJMXQCx9a9euVefH5uZm7dtvvvmmjteJREKXhjh/QC6XcyJ25Fqbmpp0/OOMirwEExRGydfMjor+6Ac2fbN1kr3I5XmMFeNGSAPAkUceCcDLmywPcXBw0AkhkfCSwcFBTdsoZpdsNqtmvkgkog+c8+LG43GtWFWtVvUlufHGQ/HVr3prfZ75b1hIywDleUB6bbnhhn/UVJEbNmzQjvG1r/2vdoaLL17qrCfyiyKD2Lp163QyIS/MihUr1BzU0tKiLwOnskyn03q/jM3z8Y9/HHPnzgXg5j3nNVsOialUKtq/RLgDcEyMYkorl8vaBxsbGx2zIfsQSATA+vXr9dgcm81Vu7jAwP/8z/9o2Bd7AVerVVx66aWbvO5NJWe48MILAXhl+WSQi8ViOOOMBBYuzDr+F2ym9ucn8OCwrKojpN/6FgBw3nleOCWbtyORiGNuFzjferFYxPe///1NXqthjEfGlZA2DMOY6IgF484779TSpBwaFY/HnUlLUPEM+Q3wJoBco1omTzKh5HAtztdQqVTUT4MdtUqlkk5iOa6bM4UVi0WntK4oFeyH0d/f76xtiwKVTCb1vJVKxbE2STs5KxhPBvmaRRnyl4FlnyDZPh6PO8qcnHP16tX6PMaKcSmkP/vZz+qsub293TGrbC6Jwbe//W0Ante0dC5gWBOYPn26k7XsO9/xnNWWLVum2nYqlUI47J2zUChiOPVhGLfc4jmCTZo0ST1o165diyuuGDYDSarRUCiEc87x2hCJRLSqUW9vr3qhX3755Zu8nksuuUTz4sbjcTXBXnTRRZvczxiJPPdoNKovPldcY0cadrKaNWuW3vempib9LI5abNbjAYUdshoaGnTJpVarqeOa3xFG+unQ0JCeX+KfAc/RUEzfbFl5J3z961/HD3/4Q72+hoYGfOMb7Y5zIldiu/nmzIhrZVM9w/HSQA0//vEanH56s5OwY1NZ4oJKE04Unn32WTUN8/WxA2MkEnGyD/JSCptvObGTPBcRhJlMRj/n83k9Xjab1e856qWjo0M/b9iwQc3vnCilublZn19zc7P2j3Q6rfv29/frslJHR4cKSbZkpdNpHd/z+by2JxaL6WdOYMLXzN7fQfeFc5rzcgo7sUmUxVgyLoU08PaF0GhCT8yCU6dO1TCTfD6vA3M+n8e3v73nW8d4FdGo1wGvumpvnHba8GfpmCtXrlQTJABcd92H3mr3X1GpeJ3kiitm47/+678AAFdeeeXbup7xmlmpHmGBwxWgpA8UCgUnHITLN0r4FmfqYqHPXqmclEcGmhdffNHxs5DBzu/tLANSR0eHagi8TUtLixY5kFCebcE3vvENAMCCBQuc7FQyyWCt6V//1bumcrmMefOGS2uyhzEAnHvuJB1Yf/xjL6znjDNa9DhB+Za5BCGfXybfhjHRGLdC2jAMYyLzgx/8AHfccQcAL8Us53KQSQtnAuPwqUKhMGpmLy5SAXgTIpnscC0AnvwlEgmnnrSYrzds2ODUhJZz9fX1OXXHuUY1Zz9kM7hM/uR/aYO0h+uac/a/0WLr+fpFk/Zn6uNr4ipjYtUaa1M3YEJaEU32kksuwcc+9jEAwOuvv64OPFyn97LL3qX1cC++eKmTFEA6VH9/v3r09vT0qGn929/eExde+DQAz6PYNOH6gwtfDA4OOjmI2TFMXup4PK4mcfaIDcqdzLS1talp+o033nAqGLEHL8eD8voaD068Fica6Pve9763dwM2wVlnnaXWqMbGRjXLy//hcNhZLw2FNuKCC3pQqVQwd+6bzrEkxacImnA47CwpyP0bGBjQOt1cSvH888/f5tdXbzzxxBMAvLSx0sfYMsNRLOFwWO8NF3Xxr1VzZAjgplrliACG661znoqOjg7NwzA4OKi5LBoaGvQ4nLp52bJl2t7p06c7eR34+P5zS9s4jt9/HYDrsc7XxMsxvK7N1i5OFiX3/cQTTxxxL3Y0JqR9sNB84IEHVAAnk0n1vs1kMpg794MAPGF8xhnDmXjEQzuTyej28XgcK1d65SzXrl2LP/3pqzviUoythAuiBDmlxONx9UsIhUKqBYTDYWfQ9CeL4LJ5a9eu1UHtzTff1AGLZ/WTJk1ycsZz0h2hVCo5gxcnseBJxfYgyLQs5R07OjqcnPmXXuoV7CiXy7jool1xzTWrUKvVcMEFPc7+/rCtQqGgIT9r167V60ulUrj44ou3/UUZRp1iQtowDKNOuf766wF4lr4Pf/jDADyHLk73KlpiS0uL46AoGuZolhzWRMUHI5PJOFWl2MQusBbb3t6uk6n169drOOKUKVN0m9WrV6v5uFqtqs9Ge3u7TmhTqZSTjpTroHNKXbmOhoaGQDO3/N7a2upMqOV+5XI5bW8ikXD8TeR4jz/+uN73esCE9CbgGNL58+fr55aWFtWSOYFKrVbTDtDf36+pIgHg7LPP3gEtNrYF+XzeCd+QgSSXyzkmbK7EI1owby9ks1ldh/vzn/+smjEnWWDP1NmzZzs5jUWrZ49V1t45iQV7+3L7tjdsepYcx62trarNJ5NJhMNhnHvuJCetrT+HtzhtDgwMaBKLa665ZkddRt1y6aWXahz/7rvv7pSEZA9sic3ntWVg2FLBpVo5FEssFclkUp9BOBzW8YwjFDj9cLlc1iiGcrmsz2zFihXa/1asWKHPubOzU7f312AQYcwJgbimQXd3t14TLzvxejqb2OVesHc7x9nn83knpemyZcv0XtcTJqS3EL+QvfXWWwF4nV463bJly7QznnLKKTu2gcY74tprr9UXOJvNOqEZXJeWi5bIzN+foKToS3vFhTk+8YlP4H//938BuOXxgOEcxJ2dnU5msaBsdYCbgU6IRqParrGqf3vZZZdt8vcf//jHOlEBvIx48j7ZZNYwXExIG4ZhjAOkkt28efPUQaulpUUnZQ0NDY6pmi0tQUViRDMWqyDgaZcyceQQUi70wU5e8XhcfTN6enr0mP39/YH+PF1dXbo9hyP29/c7XuhiKRgYGEBHh5f6NZFI6PW1tbU51ekAN3eB3APZj9Oishe5HOPll1+u20qBJqTfJqwpS5jEqaeeOlbNMd4h/gQmYp6eOnWqDjZ9fX36snd1dengwKa6SqWiGrFo1LVaTfdLJBL4h3/4BwDAH//4RydHsCSlYVN2rVZztGdpi197Zo1cPovpsd742te+NtZNGNecc845ml8BGE4/y46CXBGL4/sTiYT2VbHWcEawWq2mvzc2Nuqx/evWnFtbzsulHydPnqxrz3vttZea0FOplNNOMYmXSiU9VzweD+znyWRS2xmNRrWdbLliL26O5w9ypuRcB/Vc38CEtGHADSsZGhrScLtkMukkIpFQIy6IUSwWnTzTQXBJSpkMTJo0SXPNl0olvPvd73baI7DAZriAgBCLxXTQkrAlwzDGLyakDcMwxhkS0rnHHnto+lfWQDmdbS6X0224II9ot6LByn4y+Wtra3Nyd8v3XDGKU8PWajWdpK5Zs0YniWvWrHEcbYVwOOw4urGpXiap4XDY0bZZU2ZfETm2OLpxnm85jmwj96Jarep9rGdMSBsGPIclCbvI5XI6MPT19TlFDuRlZ42VzdmyHTCsVfNaGYevcL7g5uZmzWPMgwtr3qVSyVlb9Oe0BjzTpXip1rMJz3hnPPXUUwC8cpbSP7hIBCfz4OUSKf0JQIUcOzAODQ05JTE59EqOw9m5crmck9teSrZyUY9kMqnfl0olrF+/Xo/P+bo5BIvPKecaHBzUY/L7JhMALsPJGcz4XhQKBSfts9zHekhaMhompA3DRyqV0oEhFArpy84OM7lcTr/3py7kBCXyPxeCkO9LpZJ+7u7u1m24KhEwPGhxpilpm8ACX2r6GoYx/jEhbRiGMc6QDG/XXnst3v/+9wPwNE2ufc6Z3/h7rg4l+3FNejFxFwoFNQ0nk0lncsnarUwuGxsbVUvv6enR8/T09KjDFlfq4wlnIpFw/DnkvIVCwamgxccRJzLWmLkyFle+YhO/HO/ZZ5/V+1jPmJA2jLc499xzAXhZnkQb7uvrw8yZMwFILmq3dCTgpu4EMCJRBycbiUajulbHg9Qee+zh5BrmBA18PtmHczNHo1EdvF555RUzc+9EnH/++bjnnnsAQOsMAJ5w41AnLvspfVHMxJxZrLGxUfuhmLFlPxGGnIs7HA6rAOSUtE1NTXqepqYm/d7vxc3hY7J9oVDQfAQcEtbZ2ekU2PALaU7kw3WoBwcH9b0ZGhrC6tWr9d6NB0xIG4YPEdaAl45R1qe5+lAoFFJBWqvVnJhNdqQB3DWxvr4+rFq1CoA3GMkA1NHRoYNiKBRSoc9mbz4PO+DEYjEd1F555ZVtcxMMw6gLTEgbhmGMY6RiEyf4CIVCOunjPN7sVyGaaENDg2qd5XJZt2VHxVKppFpqc3Ozaq+ZTEYnscVi0SlJKXAJzWw2q7HUTU1Nmjs7Ho/rBLhUKukkNZFIOGl0pQ1sHWCLFn+W6/dr1XK/vvjFL27VfR4rTEgbxia49NJLtabs/vvv79TaZS9t1qploJIBY926ddiwYQMAb5DifN2c55tjo+VzNpvVAaaxsdGJx2bP2z//+c8AgPPOO2+bXbsxPpCc3ldddZUW4QCGBWW5XFaP5nK5rIJOzM6cYKRWq6nXN5fBzGQyuk7c2tqq/balpcWpycy55QX2uo5EIirguSBIqVTSZSD2MGcv7YaGBifXu1wTt5fN4TxRkDY+/fTTer/GCyakDWMzXHjhhQCAq6++GrvssgsAL+MY1z6WwWFoaAh///vfAQzHoXIhAWB4UCkWi5qViWNDo9GoU7DjmWeeAeB5gO+55566vYSy/OUvfxk362uGYWwdJqQNwzAmAJdccgluueUWAMB73/teJ3e3mKTz+bxT5hLwNFpOWiIJQdLptBMbLZabcDisTmehUEg18lKp5JxTLE1chSqVSjnma7E2DQ4OqibN6TpjsZhOWDlnQDqddip0SdvZusRm+5dfflnv0XjDhLRhbCEXX3yxfr788sudfMgCh4a0trYCcNcHq9Wqfi6Xy5qvm2OtgWFzYSQSwYwZMwB4ISMvvvgiAG/AEvPjeBx4jO3Dv/zLvwAAfvrTn2pUAhey4EIZIiBjsZgKt1qtpiblRCLhZCNj8zL3fY5skD7JkQvsTJlMJnV7zv6Vy+VUqHZ2duq+XJ6Vzen5fF7fI74mOfbQ0JAuE73++ut6X8YjJqQN423w7W9/Wz9ffvnlOsg1NzfrQMFrdTLrz+fzOqhEo1Etc1oqlZy65AKHfFUqFVx00UXb65IMw6hDTEgbhmFMML70pS9hwYIFALwqVBwOKBNKduYSLTkSiaj3dUNDgxMKKCbodDqt/hBNTU26vd9BLOhzOBzWyWs6nVaTe6VSUctTS0uLkwpX9i0UCo7jpDi7iVafz+edJChi4j7rrLPe0b0ca0xIG8Y2QAYzruPLAw0nIRFSqZQ6lPH6HKcO5QHUtGhjaxDhNGfOHMyePRsAMGvWLGdtGfD6GIc/ifCLxWKOh7isQ8diMfT39wPw6kDLNslkUvsw15Pu7+/H5MmTtV0iSPv7+1Xwd3Z2qrCPRqPOmrSsi/N6eiKR0G3knJFIRJ02n3/+eVxxxRVv/+bVESakDeMdEgqFVAhz3m0ZvLiUJQvs7u5up9qPwINUtVpVrcUwjJ0PE9KGYRgTmCuuuAIXXHABAK+GOccSA56pWfwk8vm8aq7t7e26bS6XU0esaDSqWjWbo9mZa926dfr9unXr1JmSq3C1traqBSocDmvWvHg8rg5r5XJZc33n83m1ArCzprQ3nU5ropLxFgu9KcKb32TbsXHjRnR1deGQQw7R715++WUceeSR6OrqQkdHBw477DC89NJLO7JZxjikUCjg1FNPRUtLC6ZMmTKmL+WcOXNQLpdRLpeRy+VQKpVQKpWQy+W0lF+1WkW1WkUymcSkSZP0nwxakpVJ4qXlXzqdRn9/v5oXDePtMHfuXMydOxcPPvggXn31Vbz66quaLSyTySCbzSKbzWrxjEKhgEwmo325XC5jzZo1WLNmDTZs2KD9XdLjlkol9Pb24o033sAbb7yB448/Ht3d3eju7sbxxx+v3/f29ur2YoEql8vYsGGDHl+SqJRKJWQyGW1PqVTSdnLb5XoefPBBvc6JxA7VpL/5zW/iPe95j7Mu19/fjyOOOAK33normpubcfnll+PII4/UUBPDCOK73/0uXnnlFbzxxhtYs2YNPvrRj2L27Nk4/PDDx6Q9sl589dVXOyEhgDfTF+1h5syZmDVrFgA3dGRoaEjfCzZ9L1++3NaiDWMnZrNC+oc//CH+/Oc/42c/+5l+d+655yISiWxVma8//elP+Nvf/obTTz9dA+4B4MADD8SBBx6of3/ta1/DlVdeid7eXidezpg4LFu2DB/84Afxm9/8Bvvvvz9Wr16NffbZBw888AAOPfTQLTrGHXfcgVtvvRXt7e1ob2/Haaedhttuu23MhLRhjAfmz5+vn+fMmQMA2GeffTSZSU9Pjzoz+hOJcEy/eIun02lNeZvP53H88ccHnvfYY48FANx1113YuHEjAM/0LudtaWlBe3v7iPM2NjZqe4rFIt58800AwNq1a7F06VIAUAexY4455u3elrpms0L6+OOPx3e/+1309/ejra0N5XIZixcvxiOPPIKzzz4bP/3pTwP3mzFjht7ESqWCr3zlK7j55pvx7LPPbvJ8//3f/40pU6aYgJ7A7L777vi3f/s3HHfccXjyySdxyimn4OSTT8ahhx66RX2qr68Pq1evxr777qu/7bvvvnjooYd20BWMDic8ufLKKwG4mjHX1uWwFjGHA57DmazPyaBkGNsa9n6+/PLLAQCHHXaY9kl2eOQiGewoWalUtK+efPLJmz3n8ccfj9tuuw2At+Yt69bRaFTfk2q16pxXJgT5fB7/8z//A8DLU/D5z3/+7V34OGOzQrqnpwf/+I//iPvvvx+nnXYaHn30UUyaNAkHHHAADjjgAGdmNhrXXXcdDjroIBxwwAGbFNIrV67EV77ylQm3pmCM5LTTTsPDDz+Mgw46CKFQCD//+c8BeDP9zfUpia2UuEr5LINFvXDppZfq56uuugqAF1cqg04oFHKqD3Hu7hdeeAEA8I1vfGNHNtkwjDpji9akTzrpJCxYsACnnXYa7rrrLpxwwglbfILVq1fjuuuuw5NPPrnJ7davX49PfvKTOPvss8dNCTFhv/32G+smjEtOO+00HHHEEbjpppscb83NIYUt0um0aqXpdFo9SI3xh71DY4tk0DvzzDOx9957A/CSoIg5mlN/1mo1NXH//e9/x9lnn71V5xKNe/78+eqfMWnSJCe7HntsS1KSv/3tb7jhhhu29tLGPaEa5yAchaGhIfT09OD3v/89/uEf/gHPP/88ZsyYgTPPPBN33XVX4D677rornnvuOTz00EM49thj1WVf3PQ7OjqwatUqRCIR9PX14WMf+xgOO+wwfP/739+mF2jUJ5lMBvvuuy8++tGP4pFHHsGzzz6Ljo6OLepTADB16lTcfvvt+MQnPgHAG2Refvll3HvvvTvsGrYGKXd5+OGHqymvUCjo5KRcLuuEY+XKlfjUpz41Ng01DEISoiSTSRWilUoFP/nJTza7r/iXPPbYY5vd9rzzzlNzOuf0lqxpOzNbFIKVTCZx9NFH40tf+hIOPPBATfh/ww03OK7w/E8G00996lN4/fXX8cwzz+CZZ57B5Zdfjve///145plnEIlEkE6ncdhhh+H//J//YwJ6J+K8887DAQccgIULF+Kf/umfcOaZZwLYsj4FACeeeCKuvPJK9PX14cUXX8TNN9+8RWtiY0UikdBiBLVaDbVaTTM6RSIRxGIxDA4OYnBw0EIQDcNQtjgE66STTsLChQuxaNGirTpBIpHQIgLAcAC7fPfggw/iiSeewHPPPacOBQBUWzcmHkuWLMGjjz6q/glz587Ffvvth7vvvhvHHXfcFh3jsssuw1lnnYVdd90VqVQK3/zmN82z2zC2MTtKk90SzXxnZYvM3YAXr/nud78ba9as0XUKwzC2jueff96JjxZP2nK5jGXLlgEAjj766DFpm2FsS7bG3G2MzhZp0tVqFXPnzsWxxx5rAtow3gGzZ8/GhRdeqH/LHLlSqVhUg2EYI9iskM5ms+ju7sauu+6KRx99dEe0yTAMwzAMbIW52zAMwzC2FDN3bxt2aIENwzAMwzC2HBPShmEYhlGnmJA2DMMwjDrFhLRhGIZh1CkmpA3DMMaYu+++G01NTfqvoaEBoVBoszUPjImPCWnDMIwx5rjjjnNS4Erxif3333+sm2aMMSakDcMwtpLFixc7mm8ikdCQo23B7bffjhNPPNGpDGXsnFictGEYxjsgnU7joIMOwvnnn4++vr5NFgrq7+/f7PHeeOMNzJo1C6+++ip22223bdjSHcv5558PALj22mvHtB3jHRPShmEYb5NqtYojjjgC06dP32bFKK644gr853/+pyUBMQCYudswDONt861vfQuDg4O47rrrtnif5cuXO6ZyP3fccQdOOumkbdlMYxyzxaUqDcMwjGHuvfde3HPPPXjiiScQi8UAAFdddRWuuuqqUffJZDKYMWMGMplM4O9//OMfsXr1aquEZihm7jYMw9hKnn76aXzyk5/Er3/9a+y3337b7Linn346hoaGcMcdd2yzYxrjGzN3G4ZhbCVLlixBX18fDjnkEDVbf+pTn3pHxxwaGsJ9991npm7DwTRpwzAMw6hTTJM2DMMwjDrFhLRhGIZh1CkmpA3DMAyjTjEhbRiGYRh1iglpwzAMw6hTTEgbhmEYRp1iQtowDMMw6hQT0oZhGIZRp5iQNgzDMIw6xYS0YRiGYdQpJqQNwzAMo04xIW0YhmEYdYoJacMwDMOoU0xIG4ZhGEad8v8DITcud/IPwW8AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 475.2x187.2 with 4 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# mask VMPFC\n",
"mask_file = '/gpfs/gibbs/pi/levy_ifat/Or/ROI/vmpfc_association-test_z_FDR_0.01.nii.gz'\n",
"mask_file = nilearn.image.math_img(\"a>=3\", a=mask_file)\n",
"%matplotlib inline\n",
"nilearn.plotting.plot_roi(mask_file)\n",
"masker = nilearn.input_data.NiftiMasker(mask_img=mask_file, \n",
" smoothing_fwhm=None, standardize=True, \n",
" detrend=False, verbose=0)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "eee121f0-cee1-4113-8c73-4ef8c05409e1",
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2002/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2005/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2027/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2028/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2032/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2044/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2045/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2052/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2059/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2008/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2021/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2022/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2023/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2025/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2033/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2037/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2039/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2043/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2047/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2050/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2051/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2056/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2001/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2013/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2017/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2038/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2042/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2048/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2053/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2055/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2058/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2004/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2010/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2012/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2015/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2024/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2026/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2034/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2036/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2062/modelestimate/results/zstat1.nii.gz\n",
"Running /gpfs/gibbs/pi/levy_ifat/Or/Chadi_Data/results/modelfit/_subject_id_2064/modelestimate/results/zstat1.nii.gz\n"
]
},
{
"data": {
"text/plain": [
"(41, 4485)"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"responders = []\n",
"for func in resp_func:\n",
" print(f'Running {func}')\n",
" beta = masker.fit_transform(func)\n",
" responders.append(beta)\n",
"\n",
"nonResponders = []\n",
"for func in nonResp_func:\n",
" print(f'Running {func}')\n",
" beta = masker.fit_transform(func)\n",
" nonResponders.append(beta)\n",
"\n",
"respArr = np.array(responders)\n",
"respArr_reshape= np.array(respArr).reshape(respArr.shape[0], respArr.shape[2])\n",
"respArr_reshape.shape\n",
"\n",
"\n",
"nonArr = np.array(nonResponders)\n",
"nonArr_reshape= np.array(nonArr).reshape(nonArr.shape[0], nonArr.shape[2])\n",
"nonArr_reshape.shape\n",
"\n",
"\n",
"## Create condition labels (1 = responder, 0 = non responder)\n",
"label1 = [1] * respArr.shape[0]\n",
"label2 = [0] * nonArr.shape[0]\n",
"condition_label = np.concatenate([label1, label2])\n",
"condition_label\n",
"\n",
"X = np.concatenate([respArr, nonArr])\n",
"X = X.reshape(X.shape[0], nonArr_reshape.shape[1])\n",
"X.shape"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "35979000-60e9-4ad8-917b-ba67dc037de9",
"metadata": {},
"outputs": [],
"source": [
"#from sklearn.model_selection import StratifiedKFold\n",
"model = XGBClassifier(n_jobs=12, \n",
" random_state=None, booster = 'dart')\n",
"\n",
"## Here we use stratified K-fold with shuffling to generate different shuffling of leave one subject out\n",
"cv = StratifiedKFold(n_splits=18, shuffle=True) # running for each subject\n"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "e688a4b4-4f7a-4b5d-b491-f563feecfb79",
"metadata": {},
"outputs": [],
"source": [
"scores = cross_val_score(model,\n",
" X,\n",
" y=condition_label,\n",
" cv=cv,\n",
" groups=condition_label,\n",
" scoring= \"roc_auc\",\n",
" n_jobs=5, # set number of CPUs\n",
" #verbose = 5 # set some details of the activity \n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "ee0210ea-36e5-4335-819f-fffbe0deb8d4",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1. 0. 0.5 1. 0. 1. 0. 1. 0. 1. 0. 1. 1. 0. 0. 1. 1. 1. ]\n",
"0.5833333333333334\n"
]
}
],
"source": [
"print(scores)\n",
"print(np.mean(scores))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "856c3c02-c364-4fbb-90ee-6d5eeb6a3fb4",
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Running 1 iteration\n",
" Running 2 iteration\n",
" Running 3 iteration\n",
" Running 4 iteration\n",
" Running 5 iteration\n",
" Running 6 iteration\n",
" Running 7 iteration\n",
" Running 8 iteration\n",
" Running 9 iteration\n",
" Running 10 iteration\n",
" Running 11 iteration\n",
" Running 12 iteration\n",
" Running 13 iteration\n",
" Running 14 iteration\n",
" Running 15 iteration\n",
" Running 16 iteration\n",
" Running 17 iteration\n",
" Running 18 iteration\n",
" Running 19 iteration\n",
" Running 20 iteration\n",
" Running 21 iteration\n",
" Running 22 iteration\n",
" Running 23 iteration\n",
" Running 24 iteration\n",
" Running 25 iteration\n",
" Running 26 iteration\n",
" Running 27 iteration\n",
" Running 28 iteration\n",
" Running 29 iteration\n",
" Running 30 iteration\n",
" Running 31 iteration\n",
" Running 32 iteration\n",
" Running 33 iteration\n",
" Running 34 iteration\n",
" Running 35 iteration\n",
" Running 36 iteration\n",
" Running 37 iteration\n",
" Running 38 iteration\n",
" Running 39 iteration\n",
" Running 40 iteration\n",
" Running 41 iteration\n",
" Running 42 iteration\n",
" Running 43 iteration\n",
" Running 44 iteration\n",
" Running 45 iteration\n",
" Running 46 iteration\n",
" Running 47 iteration\n",
" Running 48 iteration\n",
" Running 49 iteration\n",
" Running 50 iteration\n",
" Running 51 iteration\n",
" Running 52 iteration\n",
" Running 53 iteration\n",
" Running 54 iteration\n",
" Running 55 iteration\n",
" Running 56 iteration\n"
]
}
],
"source": [
"n_iter = 1000\n",
"rand_score = []\n",
"for i in range(n_iter):\n",
" print(f' Running {i+1} iteration')\n",
" mean_scores = []\n",
" scores = cross_val_score(model,\n",
" X,\n",
" y=condition_label,\n",
" cv=cv,\n",
" groups=condition_label,\n",
" scoring= \"roc_auc\",#\"f1\",#\"accuracy\",\n",
" n_jobs=12, # set number of CPUs\n",
" #verbose = 5 # set some details of the activity \n",
" )\n",
" mean_scores.append(scores.mean())\n",
" rand_score.append(mean_scores)"
]
},
{
"cell_type": "markdown",
"id": "f3a49b7a-f378-4201-b6d8-6ade8477ad67",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Plotting area under ROC curve ditribution and printing average and standard deviation of the distribution"
]
},
{
"cell_type": "code",
"execution_count": 45,
"id": "378ae1af-9cb4-4e62-9ef4-2c5f7483420d",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Area under curve: 0.54 (+/- 0.15)\n",
"95% CI is [0.39958333 0.68333333]\n"
]
},
{
"data": {
"text/plain": [
"<AxesSubplot:ylabel='Density'>"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD4CAYAAADmWv3KAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAn20lEQVR4nO3deXxVd53/8dfnZiX7vieEsAcCAQIU22KttqWtLai0tnUZHS3OqNUZ/Y3Wzm9m6jKLv/k9avuzU7U6U+vWnVql2E3a0k0grIGwBchK9n0hyV2+vz8SKsVAbsI999zl83w87oOQ3NzzPo+bvDmc8z3frxhjUEopFXocdgdQSillDS14pZQKUVrwSikVorTglVIqRGnBK6VUiIq0O8C5MjIyTHFxsd0xlFIqaOzevbvDGJM50dcCquCLi4uprKy0O4ZSSgUNEam70Nf0FI1SSoUoLXillApRWvBKKRWitOCVUipEacErpVSI0oJXSqkQpQWvlFIhSgteKaVClBa8UkqFqIC6k1Upb/xmR73PXuuO1UU+ey2lAo0ewSulVIjSgldKqRBlacGLSIqIPC0iR0TksIissXJ7Siml/szqc/APAC8YYzaKSDQQZ/H2lFJKjbOs4EUkCVgLfAbAGDMKjFq1PaWUUu9l5SmaEqAdeERE9orIz0Qk/vwnicgmEakUkcr29nYL4yilVHixsuAjgeXAj4wxy4BB4O7zn2SMedgYU2GMqcjMnHBREqWUUtNgZcE3Ao3GmB3jf3+ascJXSinlB5adgzfGtIhIg4jMN8YcBT4IVFu1PaUCgd6EpQKJ1aNo7gJ+PT6C5iTwWYu3p5RSapylBW+M2QdUWLkNpZRSE9M7WZVSKkRpwSulVIjSgldKqRClBa+UUiFKC14ppUKUFrxSSoUoLXillApRWvBKKRWitOCVUipEacErpVSI0oJXSqkQpQWvlFIhSgteKaVClBa8UkqFKC14pZQKUVrwSikVorTglVIqRGnBK6VUiLJ6TValglZTzxn2N/RQ3zWE0+UhKymGhblJlOUnIyJ2x1NqUlrwSp3D7TFsOXCaR9+uZU99z4TPKUidwV+tKeZTa2YSGxXh34BKTYEWvFLj9jX0cPczBzjS0s/szHi+uW4BV8zJoCQznqgIBy29w+ys7WLznkb+dethfvGnWh64bRnLi1Ltjq7UhLTgVdjzeAwPvlrD/a8cIycplv93+zI+XJaLw/He0zBF6XEUpcexcUUBbx7v4O7NB7j1x+/wnfWLuWN1kU3plbowLXgV1oadbr7+5H6er2pmQ3ke39mwmKTYqEm/74q5GTx/15V89Ym93PNsFcNON399xSw/JFbKe1rwKmyNujzc+YtK3qzp4J4bFnDnlSVTuniaHBfFw5+q4K7H9vCdLdWkxk/+D4NS/qTDJFVYcnk8/HpHHW/WdPCfG5eyae3saY2MiY508OAdy7msJI1vPlNFfdeQBWmVmh5LC15EakWkSkT2iUilldtSylvGGDbvaeJ42wDf/+gSNq4ouKTXi4pw8KNPrCAnKZbHd9Yz7HT7KKlSl8YfR/AfMMaUG2Mq/LAtpSb1+rF29jX0cE1pNreuLPTJa6bGR/PAbeX0DTv5/f7TPnlNpS6VnqJRYeVYaz8vV7dSXpjCVfMyffray4pSuWp+FnsbejjS3OfT11ZqOqwueAO8JCK7RWTTRE8QkU0iUikile3t7RbHUeGsb9jJk5UNZCfFsqE835K7Ua+an0lmYgxbqppxuj0+f32lpsLqgr/cGLMcuB74koisPf8JxpiHjTEVxpiKzEzfHlEpdZbHGJ7Z3YjT7eG2VYVER1rzox/pcHDz0jy6BkfZflwPWJS9LC14Y8zp8T/bgGeBVVZuT6kLqazt5njbANcvziUrMdbSbc3OTGBxXhJvHOtgYMRl6baUuhjLCl5E4kUk8ezHwLXAQau2p9SF9A87eeFQMyUZ8ayeleaXbV5TmoPL4+G1o21+2Z5SE7HyCD4beFNE9gM7geeNMS9YuD2lJvSHgy043Yb1Fp13n0hmYgzLi1LZcaqLnqFRv2xTqfNZVvDGmJPGmKXjj0XGmH+1altKXUhN2wD7GnpYO3fs4qc/Xb0gC2MMb9R0+HW7Sp2lwyRVyHK6PTy3r4n0+Giumu//C/gpcdGUF6ZQWdvFoJ6LVzbQglch6+2aDjoHR7m5PI+oCHt+1NfOzcTpNrx9otOW7avwpgWvQtLQqIvXj7ezICeRuVmJtuXISoplYW4SO0516rh45Xda8CokbT/WzojTw7WlOXZH4X2z0xkadXOgscfuKCrMaMGrkNN7xsnbJzopL0whJ9naMe/eKMmIJzsphndOdGKMsTuOCiNa8CrkbDvSijHwwYXZdkcBQES4rCSd073DOp2w8isteBVS2vtH2F3XzaqSNNLio+2O867ywhSiIx1U1nXbHUWFES14FVJeO9pGhEP4wPwsu6O8R0xkBGX5yVQ19TLq0outyj+04FXI6BocZX9jD6uK00iICbzVKFcUpTLq8nCwqdfuKCpMaMGrkPHG8XYE4Yq5gTkr6cz0ONLjo9ldr6dplH9owauQ0D/sZHddN8uKUkieEZiLX4sIK2amcqpjkK5BnZ9GWU8LXoWEt2o6cXsMa328SpOvLStKRYDderFV+YEWvAp6Z0bd7DjVyeL8ZDIS/Duh2FQlz4hiTlYCe+u78eiYeGUxLXgV9Hac6mTE5eH9AX70ftaKman0nHFysn3Q7igqxGnBq6Dm9hj+dLKTOZkJ5KXMsDuOVxbmJhET6dCpC5TltOBVUDt4upe+YRfvm51udxSvRUU4KM1N4tDpPlweHROvrKMFr4LaOyc6SYuPZl6OfTNGTkdZfjJnnG5OtOlpGmUdLXgVtBq7h6jvGmJNSToOPy3F5ytzshOIjXJQ1dRjdxQVwrTgVdB6+0QnMZEOVsxMtTvKlEU6HJTmJlPd3IdL54lXFtGCV0Gpb9hJVWMvy2emEhsVYXecaVlSkMyw08PxtgG7o6gQpQWvgtLOU114jOF9JcFzcfV8szMTmBEVQZXOTaMsogWvgo7bY9hV28Xc7ATSA/zGpouJcAiL8pKobu7T5fyUJbTgVdA50tJH/7CLVcXBe/R+VllBMqMuD8da++2OokKQFrwKOjtPdZEUG8n8IBsaOZGSjLHTNIdO99kdRYUgywteRCJEZK+IbLF6Wyr01XcOUdM2QEVxGhGO4BoaOZEIh7AwN5EjLX24PTo3jfItfxzBfxU47IftqDDw2K56AFYWp9mcxHdKc8dG05zq0JuelG9ZWvAiUgDcCPzMyu2o8DDq8vBUZQMLcpMCds736ZiTlUBUhFDdrKNplG9ZfQR/P/AN4IJDBERkk4hUikhle3u7xXFUMHupuoWOgVFWhdDRO0B0pIO5WYlUn+7D6BTCyocsK3gR+TDQZozZfbHnGWMeNsZUGGMqMjODY7pXZY8ndjWQnzKDudkJdkfxudK8JPqGXTT1nLE7igohVh7BXw7cLCK1wOPA1SLyKwu3p0JYU88Z3qzpYOOKgqCbd8YbC7ITcQhU62ga5UOWFbwx5lvGmAJjTDFwG7DNGPNJq7anQtszuxsxBjauKLA7iiXiYiIpTo+nulkLXvmOjoNXAc/jMTy9u5E1JekUpsXZHccypXlJtPWP0NE/YncUFSL8UvDGmNeMMR/2x7ZU6NlZ20V91xC3rgzNo/ezSnOTAPQoXvmMHsGrgPdUZSOJMZGsW5RrdxRLpcRFk5cSqwWvfEYLXgW0gREXW6ua+fDSXGZEB+e0wFNRmptMQ9cQfcNOu6OoEBBpdwClLub5A6c543SzcUWhJa//mx31lrzudJXmJfHK4VaONOvkY+rSeXUELyLPiMiNIqJH/MqvnqpspCQznuVFKXZH8YvsxBjS4qP1rlblE94W9o+AO4DjIvIfIrLAwkxKAXCyfYDKum5uWVGIhODY94mICItykzjRNki/nqZRl8irgjfGvGKM+QSwHKgFXhaRt0XksyISOpOCqIDy9O5GHAIfXZ5vdxS/WpibhNsYXjuqU3eoS+P1KRcRSQc+A3we2As8wFjhv2xJMhXW3B7D5j1NvH9eJtlJsXbH8aui9DjioyN4ubrV7igqyHl7Dn4z8AYQB9xkjLnZGPOEMeYuIPQmBlG2e+N4Oy19w9xaYc3F1UDmEGFhbhKvHmlj1KVL+anp8/YI/mfGmFJjzL8bY5oBRCQGwBhTYVk6FbaeqmwkNS6KDy7MtjuKLRbmJtE/4uJPJzvtjqKCmLcF/70JPveOL4ModVbP0CgvV7eyvjyf6MjwHLg1J2tsKT89TaMuxUV/e0QkR0RWADNEZJmILB9/XMXY6RqlfO65facZdXu4pSK0pya4mKgIB2vnZfBydavOEa+mbbIbna5j7MJqAXDfOZ/vB+6xKJMKc0/tbqA0N4lFecl2R7HVtaU5vHiolaqmXpYUpNgdRwWhixa8MeZR4FER+Zgx5hk/ZVJh7HBzHweb+viXm0rtjmK7qxdkEeEQXjrUqgWvpuWiBS8inzTG/AooFpGvnf91Y8x9E3ybUtP2VGUjURHC+vLwGvs+kdT4aFYWp/JSdQv/67r5dsdRQWiyK1jx438mAIkTPJTymVGXh9/ua+JDC7NJi4+2O05AuLY0h2OtA9R2DNodRQWhyU7R/GT8z2/7J44KZ9uOtNI1OBqWY98v5JrSbL6zpZqXq1u5c22J3XFUkPH2Rqf/IyJJIhIlIn8UkQ4R0eX3lE89VdlIVmIMV87NsDtKwChMi2NhbhIvVbfYHUUFIW8HGV9rjOkDPgw0AvOAf7AslQo7bf3DvHasnY8uLyAyIjzHvl/ItaXZ7K7rpmNAl/JTU+Ptb9LZCcVuAB4zxnRZlEeFqWf3NOH2mLAe+34h15Rm4zGw7XCb3VFUkPG24H8vIkeACuCPIpIJDFsXS4UTYwxP7W5keVEKszN1aqPzLcpLIj9lBi8c0tM0amq8WtHJGHO3iHwf6DPGuEVkEFhvbTQVSi62clJ95yA1bQN8ZFl+wK2wFAhEhOsX5/CLd+roG3aSFKszdCvvTOVk50Lg4yLyaWAjcK01kVS4qazrJjrCwZL88L5z9WKuL8tl1O3hj4d1bhrlPW9H0fwS+L/AFcDK8YfOIqku2YjLzYGmXsryk4mJCv1FtadrWWEKOUmxbK3S0zTKe94uul0BlBqd9Uj5WFVjL6MuDxXFqXZHCWgOh7BucQ6/2VnPwIiLhBhvf3VVOPP2FM1BIGcqLywisSKyU0T2i8ghEdGbpdRf2F3XTUZCDEVpOjnpZG4oy2XU5WHbER1No7zjbcFnANUi8qKI/O7sY5LvGQGuNsYsBcqBdSJy2SVkVSGmrX+Yuq4hKmamhs2i2pdixcxUMhNj+ENVs91RVJDw9v959071hcdP5wyM/zVq/KGneNS7dtd14xBYVpRid5SgEOEQ1i3K4andDQyNuoiL1tM06uK8OoI3xrwO1AJR4x/vAvZM9n0iEiEi+4A24GVjzI7pR1WhxO0x7KnvYX5OEok67M9r15flMOz08NrRdrujqCDg7SiaO4GngZ+Mfyof+O1k32eMcRtjyhlbMGSViCye4LU3iUiliFS2t+sPbbg42tLH4IiLipl6cXUqVhWnkR4fzfMH9DSNmpy35+C/BFwO9AEYY44DWd5uxBjTA7wGrJvgaw8bYyqMMRWZmZnevqQKcpV13STGRjIvW2ednorICAc3lOXyyuFW+oeddsdRAc7bgh8xxoye/YuIRDLJ+XQRyRSRlPGPZwAfAo5MM6cKIX3DTo619rO8KJUIh15cnaoNy/IYcXl48ZDe9KQuztuCf11E7mFs8e1rgKeA30/yPbnAqyJygLFz9i8bY7ZMP6oKFXvquvGYsVEhauqWF6VSmDaD3+5tsjuKCnDeXoa/G/gcUAV8AdgK/Oxi32CMOQAsu6R0KuR4jGFnbRclmfFkJMTYHScoiQjrl+bz0Gs1tPUNk5UUa3ckFaC8HUXjYeyi6heNMRuNMT/Vu1rVdBxr7adnyMnqWel2RwlqG5bl4THwu/2n7Y6iAthFC17G3CsiHYydPz8qIu0i8s/+iadCzc5TXSTGRFKam2R3lKA2JyuRxflJPLdPC15d2GRH8H/H2OiZlcaYdGNMGrAauFxE/t7qcCq0dA+NcrSln4pivbjqCxvK86lq6uVE+8DkT1ZhabKC/zRwuzHm1NlPGGNOAp8c/5pSXttVO7YQ2MriNJuThIabluYhAs/pxVZ1AZMVfJQxpuP8Txpj2vnzMn5KTcrl8VBZ2838nERS4qLtjhMSspNiuWJOBs+ML3eo1PkmK/jRaX5NqfeoPt3HwIhLL6762MdXFtLUc4Y3jutd4OovTVbwS0Wkb4JHP1Dmj4AqNLxzopO0+GjmZuuaq750TWk2afHRPL6zwe4oKgBdtOCNMRHGmKQJHonGGD1Fo7xyoLGHuq4h1pSk49BpgX0qJjKCjSsKeOVwK239w3bHUQFmKmuyKjUtj7xVS0ykQ+9ctcjHVxbi8hie3t1odxQVYLTglaXa+obZcuA0K2amEqtrrlpidmYCq2el8fjOBjx6sVWdQwteWepXf6rD5TGsKdGLq1a6fVUR9V1DvHOy0+4oKoBowSvLDDvd/HpHPR9ckE26zjtjqXWLc0ieEcVvdtTbHUUFEC14ZZnf7TtN5+Aof31Fsd1RQl5sVAS3VhTwwqEWmnrO2B1HBQgteGUJj8fwk+0nKM1N0tMzfvKZy2cB8MibpyZ5pgoXWvDKEi9Vt3KifZC/vWo2okMj/SI/ZQY3luXy+K4G+nS1J4UWvLKAMYYfvX6Cmelx3FCWa3ecsHLnlSUMjLh4Qm98UmjBKwu8c6KT/Q09fGHtbJ010s/KCpJZPSuNR946hdPtsTuOspkWvPK5H71+gszEGD66PN/uKGHpzitLON07zNaqZrujKJtpwSufqmrs5Y3jHXzuill6Y5NNrl6QRUlmPD95/SS68Fp404JXPvXDbcdJjI3kE6uL7I4SthwO4YtXzaG6uY8XD7XaHUfZSAte+UxVYy8vVbdy55UlJMbqXHR22lCeR0lmPPe9fFTnig9jWvDKZ37wyjGSZ0Tx2cuL7Y4S9iIjHPz9h+ZxrHWALQd03dZwpQWvfGJvfTfbjrSxaa0evQeKG8tyWZCTyP2vHMelI2rCUqTdAVRgmuqcJo+8dYq46AjioiJ0PpQA4XAIX7tmHpt+uZtn9zZxS0Wh3ZGUn+kRvLpkdZ2DHG8bYO3cTGJ05ExAuaY0myUFydz/ynGGnW674yg/s6zgRaRQRF4VkcMickhEvmrVtpR9jDG8VN1KfEwkl+mcMwFHRPjW9Qtp6jnDT14/aXcc5WdWHsG7gK8bYxYClwFfEpFSC7enbHC0pZ9THYN8cEEW0ZH6H8JAtGZ2OjcuyeWh12po7B6yO47yI8t+I40xzcaYPeMf9wOHAb21MYS4PYYXDrWQkRDNyuI0u+Ooi7jnhoWIwL9tPWx3FOVHfjnkEpFiYBmwY4KvbRKRShGpbG9v90cc5SN76rtp6x/hukU5OudMgMtPmcEXr5rD1qoW3q7psDuO8hPLC15EEoBngL8zxvSd/3VjzMPGmApjTEVmZqbVcZSPjLo8vHK4laK0OEpzk+yOo7ywaW0JBakzuPf3hxh16bDJcGBpwYtIFGPl/mtjzGYrt6X8682advqHXVy/OEfnew8SsVER3HvTIo61DvDQazV2x1F+YOUoGgH+GzhsjLnPqu0o/+s94+T1Y+0syktiZnq83XHUFHyoNJubl+bxX6/WcLSl3+44ymJW3uh0OfApoEpE9o1/7h5jzFYLt6n84A8HmzEGblisi3kEi3NvPlucn8wrh1v53KO7pjVn/x06kVzQsHIUzZvGGDHGLDHGlI8/tNyD3KmOQQ409rJ2Xiap8dF2x1HTkBATyc1L82jsPsPbJ/SCayjTgcvKax5j2HLgNMkzolg7Vy+IB7Oy/GRKc5N4ubqV1r5hu+Moi2jBK6/tqu2iuXeYG8py9aamICcirC/PIybSwZOVDToZWYjS31LllaFRFy8damVWRjyL83RYZChIjI3iYysKaO4d5sVDLXbHURbQgldeeeVwK8NONzctydNhkSFkQU4Sl5Wk8daJTo636qiaUKMFrybV3HuGHSe7WF2STk5yrN1xlI9dvziXrMQYnt7dyMCIy+44yoe04NVFeYzht3ubiIuO4EMLs+yOoywQFeHg4ysLGXK62bynURfqDiFa8OqidtV20dB9hhvKcomL1vVhQlVu8gzWLcrhSEs/O0512R1H+Yj+xqoL6h928uKhFkoy4ikvTLE7Ttjx98pYa2anc6y1n61VzRRnxJOTpKfjgp0ewasL+sPBFpxuw/ryfL2wGgYcImxcUUBMVASP76zHqUMng54WvJpQTdsA+xp6WDs3k8zEGLvjKD9JjI3ilhUFtPWPsLWq2e446hJpwau/MOx089y+JtLio7lqvt6xGm7mZSdy5ZwMdpzq4tDpXrvjqEugBa/+wo9fP0Hn4Cjrl+YRFaE/IuHomkXZ5KfMYPOeJnqGRu2Oo6ZJf3vVe5zqGOShV0+wpCCZudmJdsdRNol0OLhtZSFuY3iysgG3R4dOBiMtePUuYwz/9NuDxEQ5uLFMpwIOd+kJMaxfmkdt5xCvHW2zO46aBi149a7n9p3mzZoOvnHdfBJjo+yOowLAsqJUlhWmsO1IG6c6Bu2Oo6ZIC14B0DU4yne2VFNemMIdq2faHUcFkJuX5pEWH82TlQ0MjepUBsFEC14B8L0t1fSdcfL9jy2Z8go/KrTFREVw28oiBoZdbN7TpFMZBBEteMXrx9rZvLeJL141m/k5emFV/aX81Blcuyib6uY+fu3nO2zV9GnBh7nBERf3bK5idmY8X7p6jt1xVAC7fE4Gc7MS+O6Wal2wO0howYe5+14+RlPPGf7jY0uIiYywO44KYGenMkiMjeKux/Yw7HTbHUlNQgs+jO1r6OGRt07xycuKWFmcZnccFQQSY6O479alHGsd4HvPV9sdR01CCz5MOd0e7n7mAFmJsXxj3QK746ggsnZeJl9YW8Kv/lTPCwd1qb9ApgUfph7efpIjLf18d8NiknTMu5qir187n7L8ZO7efIDm3jN2x1EXoAUfhk60D/DAH49zY1ku15Rm2x1HBaHoSAcP3FbOqMvD157Yr1MZBCjLCl5E/kdE2kTkoFXbUFPn9hj+4an9xEY6+JebS+2Oo4JYSWYC9968iHdOdvKT7SfsjqMmYOUR/M+BdRa+vpqGn75xkj31PXxn/WKyEnXFHnVpbllRwI1LcrnvpWPsb+ixO446j2UFb4zZDujijgHkWGs/9710jHWLclhfnmd3HBUCRIR/21BGVmIMX318LwMjOpVBILH9HLyIbBKRShGpbG9vtztOyHK6PXztyX0kxkbyvY8s1iX4lM8kx0Vx/23LqO8a4t7fHbI7jjqH7QVvjHnYGFNhjKnIzNTVg6zy0KsnONjUx79+ZDEZCboEn/KtVbPS+PIH5vD07kZ+t/+03XHUONsLXlnvYFMvP9x2nA3leaxbrPO8K2t85YNzWV6Uwj8+W0VD15DdcRRa8CFv2Onm60/uJz0hmm/fvNjuOCqERUY4eOC2ZRgDf//EPlxuj92Rwp6VwyQfA94B5otIo4h8zqptqQv7962HOdraz/c/toTkOL2hSVmrMC2O721YTGVdNz/cVmN3nLAXadULG2Nut+q1lXdePNTCo+/U8fkrZnHV/Cy746gwsWFZPq8fa+eH246zuiSN983OsDtS2NJTNCGqqecM33j6AGX5yTrXjPK7721YzKyMeL7y2F5a+4btjhO2JJBWZ6moqDCVlZV2xwhavxlfiMHtMfzsjZO09A3z5Q/MIV1HzSgbtPYN86PXTpCbEsvnryi54Ephd6wu8nOy0CIiu40xFRN9TY/gQ9DWg83UdQ2xvjxfy13ZJjsplg3L8qnrHOLFQzrrpB0sOwev7LG7rpt3TnRy+ex0ygtT7I6jwlx5YQr1XYO8WdNBbnIsy4pS7Y4UVvQIPoQ0dA3x3L4mSjLjdby7Chg3lOUyKyOezXubqO8ctDtOWNGCDxEtvcP8ekcdibGR3L6y6ILnO5Xyt0iHg0+sKiJ5RhS/3FFP99Co3ZHChhZ8COgdcvJX/7OTYZeHT142k/gYPfOmAktcTCSfXjMTt8fDL9+p48yorufqD1rwQW7Y6ebzv9jFqY5BPrl6JrnJM+yOpNSEshJjuX1VEe39Izz6Ti2jLr3T1Wpa8EHM5fbw5d/spbKumx98vJw5WQl2R1LqouZmJXLrykIauob41Y46nc7AYlrwQWrE5eZLv9nDK4db+fbNi7hxiV5UVcGhLD+Zjy7Pp6ZtgMd3NeDUkreMFnwQGhp18flHK3nxUCv33lTKp9cU2x1JqSlZMTONDy/Jpbq5j88/WsnQqC4UYgUt+CDTO+Tk0/+9k7dqOvjPjUv4zOWz7I6k1LS8b3YGH1mWzxvH27n9pzvoGtTRNb6mBR9EDjf3cdODb7K/sYcH71jOLRWFdkdS6pKsLE7jx59cwZHmPjb++G1OtA/YHSmkaMEHid/vP81HH3qbYaebxzet4YYyPeeuQsO1i3L45edW0z04ys0/fFNXhPIhLfgA1zfs5FubD3DXY3spzUtiy11XsGKm3u6tQsuqWWk8/5UrWZCbxFce28s/Plul5+V9QAs+gL16tI3rfrCdJ3Y1sGltCY/deRlZSbF2x1LKEnkpM3h802V8YW0Jv95RzzX3beclnaTskmjBB6AT7QP8zS9389lHdpEQE8kzf/s+7rlhIdGR+nap0BYV4eBbNyzkyS+sISEmkk2/3M3nfr6LY639dkcLSnpPewBp6jnDg9tqeLKygdhIB1+7Zh5feH8JMZERdkdTyq9WzUpjy1eu4Odv1XL/K8e49gfbWbcohy9fPYfF+cl2xwsaWvAB4EBjDz994xRbq5pxCHzqspl8+eo5ZOhc7iqMRUU4uHNtCRtXFPDIW6d45O1aXjjUwsriVG6tKOTGJbnERWuFXYyu6GST3jNOfr//ND9+/QSN3WeIiXSwsjiNNbPTSY2LtjueUgFn2Olm56kuKuu66BgYJSbSwcLcJEpzk5iXnUh0pCMsV4e62IpO+s+fHw2MuNh2pI0XDjbzx8NtjLg8ZCfFcGNZLitmphIbpadilLqQ2KgI1s7L5Mq5GdR2DrGnvpvq033sa+ghKkKYlRHP0KiLK+ZmMD87ERGdMluP4C3W1j/M9mMdvHCwme3HOxh1echIiOH6xTncUlFAVWOv/iAqNU1uj6G2c5Dq5j6Otw7QMTACQEZCDJfPSWfVrDRWFqcxJzMBR4iukaBH8H404nJTWdvN9uPtbD/WweHmPgDykmP5xOoirl88drR+dkGOg019dsZVKqhFOITZmQnMzhybSfWq+Zm8WdPBWzUdvH2ik+f2jd00lRIXRcXMVFYWp1FRnEZZfnJYjErTgr9EIy43VY297KrtZsepTv50spNhp4eoCGHFzFS+sW4+a+dmsigvSY/UlbJYXsoMbq0o5NaKQowx1HUOsau2i121XVTWdvPK4TYAYiIdlBemsGpWGuWFKSwpSCEzMfQGNWjBT4Fn/L+DB0/3caiplz313exv7H134YLZmfHctrKIK+dmcFlJuq6spJSNRITijHiKM+LfnbepvX+EytoudtV2s6u2i/96tQbP+Fnq/JQZLClIZmlhCksKklmYk0RqfHAPeLC0gURkHfAAEAH8zBjzH1Zuz1dGXG4auoY42T5IbecgpzoGOd46wOHmPgbHlxqLjnBQmpfEZ95XTMXMVFbMTCVdhzUqFdAyE2O4viyX68fnchoadXGwqY8DjT3sb+xlf0MPfzjY8p7nz8tOYG5WInOzEyhMjaMgdQb5qTOC4v4UywpeRCKA/wKuARqBXSLyO2NMta+3ZYzB5TE43R6c7rN/enC5DaPnfDw06qbvjJP+ESf9wy76h130DTvpHhyltW+E1r5h2vpH/mLa0vT4aEoyx44CSvOSWJyXzJyshLA4h6dUKIuLjmTVrDRWzUp793Pdg6McaOrlWEs/x1rHHk/sauCM88/ryIpAdmIs2UkxpCfEkB4fTUZiDGlx0cTHRBIfE0F8dOS7H8dFRxId4SAyQoiMEKIcYx9HRTiIdAgRDrHkFK6VR/CrgBpjzEkAEXkcWA/4vOAX/vMLDDuntypMdISDlLgocpJjKUiNY8XMVLISY5mZHkdxRjyz0uNJjovycWKlVKBKjY/m/fMyef+8zHc/5/EYWvqGaew+Q0PXEA3dQzR0naF9YOzAsPp0H52DIzjd0xuVmJEQQ+X//pCvduFdVhZ8PtBwzt8bgdXnP0lENgGbxv86ICJHLcx0KTKADrtD+JDuT2DT/ZmGT1i9gT/z6f7UAfJP0/72mRf6gpUFP9H/N/7inzdjzMPAwxbm8AkRqbzQWNNgpPsT2HR/Aluw7I+VJ5EbgXOXHCoAdCZ/pZTyEysLfhcwV0RmiUg0cBvwOwu3p5RS6hyWnaIxxrhE5MvAi4wNk/wfY8whq7bnBwF/GmmKdH8Cm+5PYAuK/QmouWiUUkr5jg7kVkqpEKUFr5RSIUoL/jwisk5EjopIjYjcPcHX14vIARHZJyKVInKFHTm9Ndn+nPO8lSLiFpGN/sw3VV68P1eJSO/4+7NPRP7Zjpze8ub9Gd+nfSJySERe93fGqfDi/fmHc96bg+M/c2kTvVYg8GJ/kkXk9yKyf/z9+awdOS/IGKOP8QdjF4NPACVANLAfKD3vOQn8+drFEuCI3bkvZX/Oed42YCuw0e7cl/j+XAVssTurD/cnhbG7v4vG/55ld+5L/Xk75/k3Advszn2J7889wPfHP84EuoBou7OffegR/Hu9O72CMWYUODu9wruMMQNm/N0E4png5q0AMun+jLsLeAZo82e4afB2f4KFN/tzB7DZGFMPYIwJ5Pdoqu/P7cBjfkk2Pd7sjwESZWwimQTGCt7l35gXpgX/XhNNr5B//pNE5CMicgR4HvhrP2Wbjkn3R0TygY8AP/Zjruny6v0B1oz/l/kPIrLIP9GmxZv9mQekishrIrJbRD7tt3RT5+37g4jEAesYO7AIVN7sz4PAQsZu4qwCvmqMmd7EWBbQgn8vb6dXeNYYswDYAHzX6lCXwJv9uR/4pjHGPcFzA403+7MHmGmMWQr8EPit1aEugTf7EwmsAG4ErgP+SUTmWR1smrz6/Rl3E/CWMabLwjyXypv9uQ7YB+QB5cCDIpJkbSzvacG/15SmVzDGbAdmi0iG1cGmyZv9qQAeF5FaYCPwkIhs8Eu6qZt0f4wxfcaYgfGPtwJRQf7+NAIvGGMGjTEdwHZgqZ/yTdVUfn9uI7BPz4B3+/NZxk6hGWNMDXAKWOCnfJOz+yJAID0YO1o6CczizxdVFp33nDn8+SLrcqDp7N8D7eHN/pz3/J8T2BdZvXl/cs55f1YB9cH8/jD23/8/jj83DjgILLY7+6X8vAHJjJ2rjrc7sw/enx8B945/nD3eBxl2Zz/70DXlzmEuML2CiPzN+Nd/DHwM+LSIOIEzwMfN+LsbaLzcn6Dh5f5sBP5WRFyMvT+3BfP7Y4w5LCIvAAcAD2Mrox20L/WFTeHn7SPAS8aYQZuiesXL/fku8HMRqWLslM43zdj/tAKCTlWglFIhSs/BK6VUiNKCV0qpEKUFr5RSIUoLXimlQpQWvFJKhSgteKWUClFa8EopFaL+Pxc01A26bHK5AAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"rand_score = np.array(rand_score)\n",
"print(\"Area under curve: %0.2f (+/- %0.2f)\" % (np.mean(rand_score), np.std(rand_score) * 2))\n",
"print(f'95% CI is {np.quantile(rand_score, [0.025, 0.975])}')\n",
"sns.distplot(rand_score)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "9d81d323-7369-465a-979b-aaa381791025",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Area under curve: 0.32 (+/- 0.15)\n",
"95% CI is [0.18735119 0.46428571]\n"
]
},
{
"data": {
"text/plain": [
"<AxesSubplot:ylabel='Density'>"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD4CAYAAADmWv3KAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAoxElEQVR4nO3dd3xcV5n/8c8zo1HvvViWLEtyL7HlksQlPY5TIQHSYHcDBAi9LD/awi4/Sn7AkoUFNo1dqgmkZ8FpxIlL7DiWuxzLVrEkS7as3utozu8PKYmTuIzsuXOnPO/XSy+rzdzveY30+Ojec58jxhiUUkqFHofdAZRSSllDC7xSSoUoLfBKKRWitMArpVSI0gKvlFIhKsLuACdLT083hYWFdsdQSqmgsXPnzjZjTMapvhZQBb6wsJDy8nK7YyilVNAQkfrTfU1P0SilVIjSAq+UUiFKC7xSSoUoLfBKKRWitMArpVSI0gKvlFIhSgu8UkqFKEvXwYtIHdALjAFuY0yZlcdTSin1Nn/c6HSpMabND8dRSil1koC6k1UFt3XbG+yOMCm3L5tqdwSlLGX1OXgDvCAiO0Xk7lN9g4jcLSLlIlLe2tpqcRyllAofVhf4i40xi4BrgE+LyKp3f4Mx5kFjTJkxpiwj45T9cpRSSp0DSwu8MebYxL8twJPAUiuPp5RS6m2WFXgRiRORhDffB64CKqw6nlJKqXey8iJrFvCkiLx5nHXGmOcsPJ5SSqmTWFbgjTG1wAKrnl8ppdSZ6Z2sSikVorTAK6VUiNICr5RSIUrvZFUhY3TMQ3vfCL1Do0Q4HaTEukiKcTFxoV+psKMFXgU1YwyVzb28fqSDmtY+3B7zjq8nREcwPy+JpdPSyEiIesfXrGqtoC0QVKDQAq+CVlPXIE/ubuRY1xBJMS6WTUslPzWWhGgXYx5Da98wNS19vFbbwdaadhYVpLBmTjZxUfpjr8KD/qSroGOMYdPhVl48eIK4qAhuWTyFBVOScTreeSqmODOeC4vS6Bt2s/FQC6/VdnCouZdbFk+hNCvBpvRK+Y8WeBVUxjyGJ3Y1svtoF/PykrhpYR4xkc4zPiY+KoJr5+eyqCCFR8sb+e3WOtbMzWZFcbqen1chTVfRqKAx5jH8eUcDu492cfmsTG5dkn/W4n6ynKQYPrl6OnNyE3m2opnnKpoxxpz9gUoFKS3wKigYMz5zrzjWw9p5OVw+M+ucZt+REQ5uXTqV5UWpbK5uY/3+41rkVcjSUzQqKGw41MLuo11cMSuTFcXp5/VcDhGun5+LILxa0058tIvVpdqqWoUeLfAq4B083sNLB1u4ID+ZS2dk+uQ5RYRr5+fQP+Lm+QPNpMS6mD8l2SfPrVSg0FM0KqB1Dozw2M5GcpOjuemCPJ9eFHWIcMuiKRSkxfL4rkaOdw/67LmVCgRa4FXA8hjDX3YcxWMMty2Zisvp+x/XCKeD25dOJcbl5I/bGxgeHfP5MZSyixZ4FbC2VrdR3zHA9QtySYuPOvsDzlFCtItbl0yls3+Ev+4/btlxlPI3LfAqILX1DfPCGyeYlZ3ABfnJlh+vMD2O1aUZ7KzvpKKp2/LjKeUPWuBVwDHG8MyeYzgdwo0+Pu9+JpfPyiIvOYYndzfRMzjql2MqZSUt8Crg7Gvsprq1j6vmZJMY7fLbcZ0O4YNl+bg9Hh7f1ajr41XQ0wKvAsqI28OzFcfJS45h2bRUvx8/IyGKq+dkU9XSx75GPVWjgpsWeBVQXq1po2fIzdp5OThs6hOzvCiNvOQY1u8/zpCuqlFBTAu8Chh9w242HW5lVk4i09LjbMvhEOGmhXn0Dbt54Y0TtuVQ6nxpgVcB46WDJxgd87BmTrbdUchLiWFZURrba9tp7BywO45S50QLvAoIrb3D7KjrYElh6nt2XrLLVbOziI+K4Ok9x/DoBVcVhLTAq4DwwhvNuJwOLp+VZXeUt0S7nFwzL5umrkH2HO2yO45Sk6YFXtnuePcgB471cHFxOvEBtp3e/CnJ5CXH8MKBZkbcHrvjKDUpWuCV7TZUthAV4eDi6efXBtgKDhHWzsuhZ8jNluo2u+MoNSla4JWtmruHOHCsh4ump09qdyZ/mpYex+ycRDYdbqV3SO9wVcFDC7yy1YZDE7P34jS7o5zRmrnZuD0e/n5Ql02q4KEFXtnmRM8QB5q6uXB6GrGRgXXu/d3S46NYXpRGeV0nLb1DdsdRyita4JVtXj7UgivCwYoAPPd+KpfMyMQV4eClgy12R1HKK1rglS06+kfY39jNsmmpxAbYypnTiY+K4OLpaexv6uZYl+7+pAKfFnhli1er23CIBOTKmTNZUZxBtMuh5+JVULC8wIuIU0R2i8hfrT6WCg4DI27K6ztYkJ9MYoz/2gH7Qkykk1UlGVQ299LQoS0MVGDzxwz+88BBPxxHBYntRzoYHTOsKAmu2fubLpyeRlxUBC++0Wx3FKXOyNICLyJTgGuBh608jgoeo2MettW0U5oVT3ZitN1xzklUhJNLSjOoae2nprXP7jhKnZbVM/j/AL4KnPYebxG5W0TKRaS8tbXV4jjKbnuOdtE37GZlSYbdUc7L0mmpJEZH8OIbJ3TnJxWwLCvwInId0GKM2Xmm7zPGPGiMKTPGlGVkBPcvvTozjzFsqWojNymaIhv7vfuCy+ng0pmZNHQMcOhEr91xlDolK2fwFwM3iEgd8AhwmYj8wcLjqQB3uLmX1r5hVpZk+G0jbSuVFaSSGhfJ3w/qLF4FJssKvDHm68aYKcaYQuBWYIMx5k6rjqcC36aqNpJjXMzNS7I7ik84HcJlMzI51jXEweM9dsdR6j10Hbzyi6MdA9S193NRcTpOR/DP3t+0ID+ZtLhI/n6wRTcFUQHHLwXeGPOKMeY6fxxLBabN1W1EuxwsKUixO4pPOR3C5bMyae4Z74qpVCDRGbyyXEf/CAeaullamEaUKzBbAp+P+VOSyUiI4qWDJ3QWrwKKFnhluTfbElw0PbBbAp8rhwiXz8ykpXeY/Y3ddsdR6i1a4JWlBoaDty3BZMzNSyIrMYqXKk/gHtOt/VRg0AKvLLW9LrjbEnhrfBafRVvfCM/sPWZ3HKUALfDKQqHQlmAyZucmkpMUzc9eqmJUZ/EqAGiBV5YJlbYE3nKIcMWsLOrbB3hyV5PdcZTSAq+s4TGGzVVt5CYHf1uCyZiZncD8KUn8fEMVI26dxSt7aYFXlqg83ktb3zCrQqQtgbdEhC9eWUpj5yCP7Wy0O44Kc1rglSU2VbWSEutiTm5otCWYjEtKM7hgajK/2FDFsHvM7jgqjGmBVz5X395PQ8cAK0KsLYG3RIQvXVnKse4h/rzjqN1xVBjTAq98btPhVmIjnSwuSLU7im1WFKeztDCVX75czdCozuKVPbTAK59q6RniYHMvy4vSiIwI3x+vN8/Fn+gZZt32BrvjqDAVvr+ByhJbqttwOYXlRaHZlmAyLpyexoVFafzqlRoGR3QWr/xPC7zymZ7BUXYf7WLR1BTioyLsjhMQvnhlKW19w/xuW53dUVQY0gKvfGZrTTsej2FFcWi3JZiMpdNSWVmSzn9trKF7cNTuOCrMaIFXPtEzNMr2I+3MzUsiLT7K7jgB5f+smUn34Ci/eqXa7igqzGiBVz7xu611DLs9rCoNj7YEkzE3L4n3Lczjf16to6lr0O44KoxogVfnrW/YzcNbjjAzO4G85Bi74wSkL189A4B/f/6QzUlUONECr87bH16rp2tglEtnZNodJWDlJcdw18XTeHJPExVNuimI8g8t8Oq8DIy4eWhTLatKM8hPjbU7TkD71CXTSY5x8cNnD2J0az/lB1rg1XlZt72B9v4RPndZsd1RAl5SjIvPXlbCq9XtbDzcanccFQa0wKtzNjQ6xgObarloehplheHblmAy7lxeQEFaLD9cX8mYR2fxylpa4NU5+922Olp7h/n85SV2RwkakREOvnr1TA6d6OUv5dqITFlLC7w6Jz1Do/zqlRpWl2awTNsSTMraedksLUzlR89V0jUwYnccFcK0wKtz8tCmWroGRvnnieV/ynsiwr/dOIeeITc/eUGXTSrraIFXk9baO8zDm49w3fwc5uaF34YevjArJ5EPLy/gj9sbdNmksowWeDVpv9hQxciYhy9fpbP38/HFK0tJi4vkX56uwKMXXJUFtMCrSWloH2Dd6w18sCyfaWG0mbYVkmJcfO2aWexu6OKxXbp/q/I97emqJuV7f3sDl9PBF67QlTO+8P4L8li3vZ7/92wlV8/JJinGddrvtWrjkNuXTbXkeZX9dAavvLalqo0X3jjBpy8tJisx2u44IcHhEL5741w6B0a478XDdsdRIUYLvPKKe8zDd/96gPzUGD66YprdcULK3Lwk7lhWwO+21ekFV+VTWuCVV9a93sDhE318c+1sol1Ou+OEnK9cPYPUuCi+8eR+vcNV+YxlBV5EokXkdRHZKyIHROTfrDqWslZn/wj//sJhLpqextVzsuyOE5KSYlx8+/rZ7Gvs5ve6vZ/yEStn8MPAZcaYBcBCYI2ILLfweMoiP1h/kL5hN9++fjYiYneckHX9/BxWlqTzkxcO09w9ZHccFQK8KvAi8riIXCsiXv+HYMb1TXzomnjTvz2DzNaaNh7d2cjdq4qYmZ1od5yQJiJ876a5jE5c71DqfHlbsP8LuB2oEpF7RWSmNw8SEaeI7AFagBeNMdvPLaayw9DoGN98soKCtFhtKOYnBWlxfPayYtbvb+blyha746gg51WBN8b83RhzB7AIqANeFJGtIvJPInLahbvGmDFjzEJgCrBUROa++3tE5G4RKReR8tZW7ZEdSH75cjVH2vr5/k3z9MKqH929ajrFmfH8y9MVDI6M2R1HBTGvT7mISBrwj8DHgN3Azxgv+C+e7bHGmC7gFWDNKb72oDGmzBhTlpGhGzYHisrmHu7fWMP7F+WxoiTd7jhhJTLCwfdvmktj5yA/e6nK7jgqiHl7Dv4JYDMQC1xvjLnBGPNnY8xngfjTPCZDRJIn3o8BrgAqfZJaWWrYPcYXHtlDUoyLb1072+44YWlZURofWDyFhzfXUtncY3ccFaS8ncE/bIyZbYz5oTHmOICIRAEYY8pO85gc4GUR2QfsYPwc/F/PO7Gy3E9fOExlcy8/umU+qXGRdscJW19fO4uE6Ai+8cR+bUamzom3Bf57p/jctjM9wBizzxhzgTFmvjFmrjHmu5OPp/zttdp2Htxcy+3LpnLZTF3zbqfUuEi+ee1sdjV08cgO3f1JTd4ZC7yIZIvIYiBGRC4QkUUTb5cwfrpGhZCeoVG+/Je9FKbF8a1rZ9kdRwE3L8pjeVEq9z57kN6hUbvjqCBzthn81cBPGF8F81Pg3yfevgR8w9poyt++8/QBmnuGuO9DC4mN1EajgWB8bfw8BkfHeLai2e44Ksic8bfYGPNb4LcicrMx5nE/ZVI2+Ou+Yzy5u4kvXFHCwvxku+OokxRnxvOp1dP5+YZqFk1NoTjzlOsalHqPs52iuXPi3UIR+dK73/yQT/lBc/cQ33yygoX5yXzm0mK746hTuOfSYtLiInl6TxOjYx6746ggcbZTNG9u2RMPJJziTQU5j8fwz4/tZcTt4b4PLSTCqQ1GA1G0y8mNC/No7x/hlUN6Q6DyztlO0Tww8a92ggxRv91Wx+aqNr7/vrm6BV+AK86MZ8GUJDYdbmVBfhKZCbrpijozb290+pGIJIqIS0ReEpG2k07fqCBVdaKXe5+t5PKZmdy+VLdtCwZr5+XgihCe3nMMY3RtvDozb/8ev8oY0wNcBzQCpcA/W5ZKWW7E7eHzj+whPiqCe2+er22Ag0RCtIs1c3I40tbP7oYuu+OoAOdtgX+zodha4E/GmA6L8ig/ue/vh3njeA/33jyfjIQou+OoSSgrTGFqaizrK44zMOy2O44KYN4udv5fEakEBoF7RCQD0B0JLLZue4Mlz3ukrZ+HN9dSVpBCa++wZcdR1nCIcNPCPH7xchXPHmjm5kVT7I6kApS37YK/BlwIlBljRoF+4EYrgylrDLvHeHxXIylxkVw7P8fuOOocZSdFs6I4nZ31nRxp67c7jgpQk1kTNwv4kIh8BLgFuMqaSMpKz1U009k/ws2LphAVoT3eg9llM7NIjnXxzN4mPHrBVZ2Ct6tofs94y4IVwJKJt9N1kVQBqrqlj+1HOrhoepouiQwBkREO1szJ5kTPMLsbOu2OowKQt+fgy4DZRtdlBa2h0TGe2NVIenwUV83JtjuO8pF5eUlsqW7jxTdOMC8vmcgIvVFNvc3bn4YKQKtCEHuuopnuwVE+sHgKLr1bNWSICNfMzaFnyM3Wmja746gA4+0MPh14Q0ReB4bf/KQx5gZLUimfqm/v5/W6DlYUp5Ofql2eQ8209Dhm5SSy8XArZYWpxEdpJ1A1ztufhH+1MoSyzpjH8NSeJpJiXFw+K9PuOMoiV8/J4ucv9bCh8gQ3LMizO44KEN4uk9wI1AGuifd3ALsszKV8ZEt1Gyd6hrlhQa6umglhmQnRLC5IZceRTjoHRuyOowKEt6toPg48Bjww8ak84CmLMikf6egfYUPlCWbnJDIrJ9HuOMpil87IAIGXK1vsjqIChLdX2z4NXAz0ABhjqgD9ez/A/W3/cUSE6xfk2h1F+UFybCRLC1PZ1dBJe9/w2R+gQp63BX7YGPPW330iEgHokskAVtPax8HjPVxamkFSjOvsD1AhYfWMDBwibNBZvML7Ar9RRL7B+ObbVwKPAv9rXSx1PjzGsH7/cZJjXVxUnG53HOVHidEulhelsedoFy292i4q3Hlb4L8GtAL7gU8A64FvWRVKnZ/dDV0c7x7i6jnZuuY9DK0qzcDldOgsXnm3TNIY4xGRp4CnjDG6X1gAG3F7eOGNZvJTYpifl2R3HGWD+KgILpyexqbDrVwxc5h0bQcdts626baIyL+KSBtQCRwSkVYR+bZ/4qnJ2lTVSu+Qm2vn5+omHmHs4uJ0nA5hU5XOx8LZ2f5+/wLjq2eWGGPSjDGpwDLgYhH5otXh1OQMDLt5tbqNubmJTNU7VsNafFQEZYWp7G7ookvXxYetsxX4jwC3GWOOvPkJY0wtcOfE11QA2Vzdxojbw+WzsuyOogLAypJ0DIYt1dqjJlydrcC7jDHv+emYOA+va+8CSP+wm2017cybkkRWYrTdcVQASImNZGF+MjvqOujTrf3C0tkK/Jn+ttO/+wLI5qpWRsc8XDZD7z9Tb1tVkoF7zLBNO02GpbOtolkgIj2n+LwAOk0MEH3DbrbVtrMgP5lMnb2rk2QmRjM7N5Ftte2sLMkg2qX9iMLJGWfwxhinMSbxFG8Jxhg9RRMgNh1uxT1mdPauTumS0kyGRj1sP9JhdxTlZ3oXTJAbHBnj9SMdLMhP1vXO6pTyUmIoyYxnS3Ubo2Meu+MoP7KswItIvoi8LCIHReSAiHzeqmOFs+1H2hkZ87CyRFsSqNNbPSOD/mE3O+t179ZwYuUM3g182RgzC1gOfFpEZlt4vLAzOuZha007pVnx5CTF2B1HBbBpaXFMSYnh1eo2PLq1ctiwrMAbY44bY3ZNvN8LHGS8j7zykT0NXfQNu1lZkmF3FBXgRIQVxem0949QefxU6yZUKPLLOXgRKQQuALb743jhwGMMm6tbyUuOoSg9zu44KgjMyU0iOdbFZr3xKWxYXuBFJB54HPiCMeY9UwcRuVtEykWkvLVV+2Z46+DxHtr6RlhZkq49Z5RXnA7h4unp1LcPcLRjwO44yg8sLfAi4mK8uP/RGPPEqb7HGPOgMabMGFOWkaGnGry1paqNlFgXc3K1Y6TyXllBCtEuh7YvCBNWrqIR4NfAQWPMT606Tjhq6hqkvmOAi6aPdwxUyltRLidLClOpaOqms19vRg91Vs7gLwY+DFwmInsm3tZaeLywsa2mnUing8UFKXZHUUHowqI0RGCrti8IeV5t+HEujDFbGG9poHyof9jNvsYuFhek6G3n6pwkx0YyLy+J8vpO7Twa4vRO1iBTXteB22NYXpRmdxQVxFaUZDDs9rCjTtsXhDLLZvDK98Y8hteOdFCUEactgQPYuu0Ndkc4q7zkGKalx7G1pp3RMY/u3Rui9FUNIpXNPXQPjnKRzt6VD6wsTqd7cJT1+4/bHUVZRAt8ENlW205yjIuZOYl2R1EhoDQ7gfT4KB7aXIvR9gUhSQt8kGjrG6a2tZ+l01Jx6I1NygccE+0LKpp6eK1Wz8WHIi3wQaK8rhOHwCJdGql86IKpyaTFRfLw5lq7oygLaIEPAmMew66GTmZkJ5IYrfusKN9xOR3cubyAlypbqG7pszuO8jEt8EHg4PEe+obdLCnU2bvyvQ9fWEBkhINfbzlidxTlY1rgg0B5fQeJ0RGUZiXYHUWFoPT4KG5elMcTuxpp7xu2O47yIS3wAa5zYISqE30sLtCLq8o6H1tZxLDbw2+31dsdRfmQFvgA9+YWa2V6ekZZaHpGPFfOzuL32+oYGHHbHUf5iBb4AOYxhp31nRRnxpMSG2l3HBXiPrGqiM6BUR4tb7Q7ivIRLfABrOpEL92DoywpTLU7igoDZYWpLJqazMNbanGPeeyOo3xAC3wA21HXSVxUBDNz9OKq8o9PrJ7O0Y5Bnq1otjuK8gEt8AGqpWeIyuYeFk1NJsKhL5PyjytnZVGUHseDm7R9QSjQyhGgHt3ZiMfAkgI9PaP8x+EQPrayiP1N3Wyrbbc7jjpPWuADkMdj+POOo0xLjyM9IcruOCrMvH9RHunxkTywUdsXBDst8AFoW207DR0DeueqskW0y8k/XlTIxsOtHDzeY3ccdR60wAegR3YcJTE6gjm5SXZHUWHqzuUFxEY6eXCTzuKDmRb4ANPRP8LzFc28f9EU3WVH2SY5NpLbl07lmb3HqG/vtzuOOkdaQQLME7saGRnzcOvSfLujqDD38VVFOB3C/Rtr7I6izpEW+ABijOGRHUdZmJ/MzGzdtUnZKysxmg+V5fPYzkaOdQ3aHUedAy3wAWRnfSfVLX3cukRn7yowfPKS6RgDD+gsPihpgQ8gf3r9KHGRTq5fkGt3FKUAyEuO4eZFU/jTjqO09A7ZHUdNkhb4ANE9OMrf9h/jhoV5xEVF2B1Hqbd86pLpuMc8PKQraoKOFvgA8fSeJoZGPdy+dKrdUZR6h8L0OG5amMfvX6unpUdn8cFEC3wAMMawbnsDc3ITmTdF176rwPP5K0oYHTP86hU9Fx9MtMAHgL2N3VQ293Kbzt5VgCpIi+ODZVNYt72BJl1REzS0wAeAR15vIMbl5MaFenFVBa7PXlYCwH++VGVzEuUtLfA26x0a5Zm9x7h+QQ4J0S674yh1WrnJMdy+bCqP7mzkSJve3RoMtMDb7Jm9xxgYGdPTMyoo3HPpdCKdDu578bDdUZQXtMDb7JHXjzIzO4GF+cl2R1HqrDITorlrRSHP7D3G3qNddsdRZ6EF3kYVTd3sb+rmtqVTERG74yjllU+unk56fCTf+9sbuutTgLPsjhoR+W/gOqDFGDPXquMEsz+93kBUhIObLsizO4oKY+u2N0z6MSuKM3hqTxPffLKCuXmnXtp7+zI97Wg3K2fwvwHWWPj8Qa1/2M3Te45x7fwckmL04qoKLosLUshKjOK5A824xzx2x1GnYVmBN8ZsAjqsev5g98TuJvqG3dyxrMDuKEpNmtMhrJ2bQ0f/CK/p3q0By/Zz8CJyt4iUi0h5a2ur3XH8whjD77bWMS8viUVTk+2Oo9Q5KclKoDQrnpcqW+gZGrU7jjoF2wu8MeZBY0yZMaYsIyPD7jh+sa2mnaqWPv7hokK9uKqC2nXzcxnzGNbvP253FHUKthf4cPSbrXWkxkVy3fwcu6ModV7S46NYXZrBvsZuqlp67Y6j3kULvJ81dg7w94MnuHVJPtEup91xlDpvq0ozSIuL5Jk9xxjVC64BxbICLyJ/ArYBM0SkUUQ+atWxgskfXhtfknbncr24qkKDy+ngxoV5tPePsPFweFxHCxaWrYM3xtxm1XMHq6HRMR7Z0cBVs7PJTY6xO45SPlOcGc+CKUlsPNTKnNxEcpL05zsQ6CkaP3p0ZyNdA6P848WFdkdRyueum59LTKSTR8sbdW18gNAC7ydvbnm2MD+ZZdNS7Y6jlM/FRUXw/gvyaO4Z4qXKFrvjKLTA+81zB5pp6Bjgk6un69JIFbJm5iRSVpDCpsOt7KzX+xztpgXeD4wx3L+xhqL0OK6cnWV3HKUstXZeDsmxLr70l716A5TNtMD7wavV7VQ09XD3qiKcDp29q9AW7XLywbJ8GjsH+eqj+7TjpI20wPvB/RtryEyI4n2LtGukCg8FaXF8bc1MnjvQzK+3HLE7TtjSAm+xvUe72FLdxl0rphEVoTc2qfDxsZXTuHpOFvc+W0l5nZ6Pt4MWeIv95IVDpMZF6o1NKuyICD/+wALyUmL49LpdHO8etDtS2NECb6HXatvZXNXGPZdMJz7KsnvKlApYidEuHvjwYvqHx7jrN+X06kVXv9ICbxFjDD95/hBZiVE6e1dhbWZ2Ir+8YxGHT/Ty6XW7tV+NH2mBt8grh1opr+/kc5eXaFMxFfZWl2bw/ZvmsulwK99+ukJX1viJnjewgMdj+PHzh5iaGssHy/LtjqNUQLh16VSOdg7wy5driIuM4JvXztKb/iymBd4CT+1p4o3jPfz0gwtwOfWPJKXe9JWrZtA35ObhLUdwOoWvrZmpRd5CWuB9rHtwlB+sP8jC/GRuWqjr3pU6mYjwrzfMYcwYHthYS4RD+MpVM7TIW0QLvI/99IVDdPSP8Jt/WopD71pV6j1EhO/eMJcxj+GXL9fQO+TmO9fP0bu8LaAF3ocqmrr5/Wv1fHh5AXPzkuyOo1TAcjiE7980b3wZ5aZaTvQM8bNbL9AFCT6mJ4h9xOMxfOupClLjovjSVTPsjqNUwHM4hK+vncW3r5vNC2+c4I6Ht9PWN2x3rJCiBd5Hfv9aPXuOdvGNtTNJinHZHUepoHHXimn88vZFVDR1c93Pt7CzvtPuSCFDC7wP1LT28cNnD3LJjAzed4FeWFVqstbOy+GJey4iMsLBrQ9u47db63StvA9ogT9Pw+4xvvjnPUS7nPzo5vm6GkCpczQnN4n//cwKVpVk8J1nDvDx35XT0jtkd6ygpgX+PP3gbwfZ19jNve+fT2ZitN1xlApqSbEuHvpIGd+6dhabqtq46r5N/HXfMbtjBS0t8Ofhqd1N/HZbPR9dMY01c7PtjqNUSHA4hI+tLGL951ZQkBrLZ9bt5qO/2UFD+4Dd0YKOFvhzVF7XwVcf28eyaal87ZqZdsdRKuQUZybw+Kcu4uvXzGRbbTtX3LeR+148zODImN3RgoYW+HNQdaKXj/+unNzkaO6/c7G2I1DKIhFOB59YPZ0NX76Eq+dk87OXqlj945f54/Z67UrpBb3RaZLq2vq54+HtRDgd/OaflpISF2l3JKUC0rrtDT59vguL0shJjGZXQyfffLKChzbVcs+lxdy4MFd3SzsNnXpOwqHmXj7wwDZGxzz84aPLKEyPszuSUmGlMD2ORz95Ib/+hzJiIiP46mP7WPWjl3lgYw09upnIe+gM3ktbqtq45487iYl08pdPXEhJVoLdkZQKSyLC5bOyuGxmJpuq2nhwUw0/fLaS/9xQzW1L87lt6VSKMuLtjhkQtMCfhcdjeGhzLT96/hDFGfE8/A9l5KfG2h1LqbAnIqwuzWB1aQYVTd3cv7GG/361joc2H2HptFRuXZLPNXNziIkM39M3WuDP4GjHAF9/Yj9bqtu4Zm42P/7AAt1bVakANDcviV/cvoiWniEe39XEn3c08KW/7OU7Tx9g7bwcrpmXzcXF6WG3IEKr1Sl0D4zy6y21PLi5FqcIP3jfPG5bmq93qSoV4DITo/nUJdP55Ooith/p4C87jvK3/cf5c/lRkmJcXDk7i2vmZnPh9DRiI0O//IX+CCeho3+E32yt43+2HKF32M3aedl869rZ5CbH2B1NKTUJIsLyojSWF6UxNDrGlqo21lcc5/kDzTy2s5FIp4PFBSmsKElnVUkGc3ITQ3L/hrAv8AMjbjZUtvDU7iZeOdSK22O4Zm42n7u8hFk5iXbHU0qdp2iXkytmZ3HF7CxG3B62H2lnc1Ubm6va+PHzh/jx84dIiI5gwZRkFuZPvE1NJj0+yu7o583SAi8ia4CfAU7gYWPMvVYe72yMMRzrHuLwiV5213eyrbadPUe7GB0zZCVGcdeKadyyeAqlukJGqZAUGeFgZUkGK0syAGjtHebV6jZ21HWw52gX/7WxhjHPeBfL9PgoijPjKM6MpzgjnqKMeHKTY8hNjg6a0zuWpRQRJ/BL4EqgEdghIs8YY97w9bGOdgzQN+ymf9g98e8YPUOjtPQM09wzREvPECd6h6hrG/8+AIfAvLwk7loxjVUlGSwvStMtw5QKMxkJUdx0QR43TbT5HhwZY39TN/sauzh8opfqlj6e3nOM3iH3Ox6XHOsiJymG9PhIkmMjSY5xkRLrImni/bioCKJcDqIjnO/81+UkKsKBy+HA6RScIjgdQoRDLDlFZOV/Q0uBamNMLYCIPALcCPi8wF/x040Mu09923J6fCSZCdFkJUaxaGoKpVkJlGYlMCsngYRo3ZhDKfW2mEgnS6elsnRa6lufM8bQ2jtMbVs/x7sHOdY19Na/Hf0jNHYO0jUwQvfgKJ5zbGGfHh9J+beu9NEo3mZlgc8Djp70cSOw7N3fJCJ3A3dPfNgnIod8GaL+7XfTgTZfPncAC5ex6jhDz1nHeoefgljsHeOsB+Rfzvm5Ck73BSsL/Kn+3njP/2/GmAeBBy3MMR5GpNwYU2b1cQJBuIxVxxl6wmWs/hqnlav+G4H8kz6eAmjnfqWU8hMrC/wOoEREpolIJHAr8IyFx1NKKXUSy07RGGPcIvIZ4HnGl0n+tzHmgFXH84Llp4ECSLiMVccZesJlrH4Zp+jO5UopFZrCq/OOUkqFES3wSikVokKuwIvIGhE5JCLVIvK1U3x9pohsE5FhEfmKHRl9wYtx3iEi+ybetorIAjty+oIXY71xYpx7RKRcRFbYkfN8nW2cJ33fEhEZE5Fb/JnPV7x4PS8Rke6J13OPiHzbjpy+4M1rOjHePSJyQEQ2+jSAMSZk3hi/mFsDFAGRwF5g9ru+JxNYAnwf+IrdmS0c50VAysT71wDb7c5t4Vjjeft60nyg0u7cVozzpO/bAKwHbrE7t0Wv5yXAX+3O6qexJjN+d//UiY8zfZkh1Gbwb7VHMMaMAG+2R3iLMabFGLMDCOYNHL0Z51ZjTOfEh68xfh9CMPJmrH1m4rcDiOMUN9QFgbOOc8JngceBFn+G8yFvxxkKvBnr7cATxpgGGK9PvgwQagX+VO0R8mzKYqXJjvOjwLOWJrKOV2MVkfeJSCXwN+AuP2XzpbOOU0TygPcB9/sxl695+7N7oYjsFZFnRWSOf6L5nDdjLQVSROQVEdkpIh/xZYDg6HnpPa/aI4QAr8cpIpcyXuCD8rw03re8eBJ4UkRWAf8XuMLqYD7mzTj/A/g/xpixIN5dzJtx7gIKjDF9IrIWeAoosTqYBbwZawSwGLgciAG2ichrxpjDvggQagU+XNojeDVOEZkPPAxcY4xp91M2X5vUa2qM2SQi00Uk3RgTTA26vBlnGfDIRHFPB9aKiNsY85RfEvrGWcdpjOk56f31IvKrIHw9wbvXtBFoM8b0A/0isglYAPikwIfaKZpwaY9w1nGKyFTgCeDDvpoN2MSbsRbLRNUTkUWMX9AKtv/QzjpOY8w0Y0yhMaYQeAy4J8iKO3j3emaf9HouZbxOBdvrCd7Vo6eBlSISISKxjHfcPeirACE1gzenaY8gIp+c+Pr9IpINlAOJgEdEvsD4le2e0z1voPFmnMC3gTTgVxO/K24ThF36vBzrzcBHRGQUGAQ+dNJF16Dg5TiDnpfjvAX4lIi4GX89bw221xO8G6sx5qCIPAfsAzyM73xX4asM2qpAKaVCVKidolFKKTVBC7xSSoUoLfBKKRWitMArpVSI0gKvlFIhSgu8UkqFKC3wSikVov4/NRlSUmmeNJEAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"rand_score = np.array(rand_score)\n",
"print(\"Area under curve: %0.2f (+/- %0.2f)\" % (np.mean(rand_score), np.std(rand_score) * 2))\n",
"print(f'95% CI is {np.quantile(rand_score, [0.025, 0.975])}')\n",
"sns.distplot(rand_score)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "bab61c42-ee7a-4227-b5f7-9f25eae2bca6",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"[Parallel(n_jobs=10)]: Using backend LokyBackend with 10 concurrent workers.\n",
"[Parallel(n_jobs=10)]: Done 21 tasks | elapsed: 12.1min\n",
"[Parallel(n_jobs=10)]: Done 50 out of 50 | elapsed: 25.1min finished\n"
]
}
],
"source": [
"# use scikit learn permutation \n",
"score_rand, perm_scores_rand, pvalue_rand = permutation_test_score(\n",
" model, X, condition_label, scoring=\"roc_auc\", cv=cv, n_jobs=10, verbose=2,\n",
" n_permutations=50)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "f641a75e-0988-4240-a62e-0fcd84f35d3d",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.725 , 0.65 , 0.675 , 0.5 , 0.36666667,\n",
" 0.43333333, 0.41666667, 0.3 , 0.59166667, 0.4 ,\n",
" 0.51666667, 0.675 , 0.55833333, 0.38333333, 0.28333333,\n",
" 0.5 , 0.475 , 0.575 , 0.3 , 0.4 ,\n",
" 0.725 , 0.375 , 0.48333333, 0.325 , 0.675 ,\n",
" 0.39166667, 0.43333333, 0.5 , 0.45833333, 0.30833333,\n",
" 0.35 , 0.49166667, 0.625 , 0.675 , 0.29166667,\n",
" 0.23333333, 0.4 , 0.39166667, 0.66666667, 0.54166667,\n",
" 0.55833333, 0.475 , 0.7 , 0.50833333, 0.45 ,\n",
" 0.5 , 0.78333333, 0.35 , 0.40833333, 0.75833333])"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"perm_scores_rand"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"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.9.9"
}
},
"nbformat": 4,
"nbformat_minor": 5
}