SubjectSelection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Subject Selection\n",
    "Code to generate final subject lists"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import os\n",
    "import glob\n",
    "import json\n",
    "from scipy.stats import ttest_ind\n",
    "import pingouin as pg\n",
    "\n",
    "# directories\n",
    "PROJECT_DIR = '/projects/f_mc1689_1/ClinicalActFlow/'\n",
    "TASKS = ['rest','bart','scap','stopsignal','taskswitch']\n",
    "FMRIPREP_DIR = PROJECT_DIR+'data/prepro/fmriprep/'\n",
    "data_dir = PROJECT_DIR+'data/prepro/ciftify/'\n",
    "outdir = PROJECT_DIR+'docs/scripts/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data after preprocessing: 261\n"
     ]
    },
    {
     "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>rest</th>\n",
       "      <th>bart</th>\n",
       "      <th>scap</th>\n",
       "      <th>stopsignal</th>\n",
       "      <th>taskswitch</th>\n",
       "      <th>group</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>sub-70055</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>ADHD</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sub-10189</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>CTRL</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sub-60077</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>BPLR</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sub-50032</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>SCHZ</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sub-70015</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>ADHD</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           rest  bart  scap  stopsignal  taskswitch group\n",
       "sub-70055   1.0   1.0   1.0         1.0         1.0  ADHD\n",
       "sub-10189   1.0   1.0   1.0         1.0         1.0  CTRL\n",
       "sub-60077   1.0   1.0   1.0         1.0         1.0  BPLR\n",
       "sub-50032   1.0   0.0   1.0         1.0         1.0  SCHZ\n",
       "sub-70015   1.0   1.0   1.0         1.0         1.0  ADHD"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CTRL    122\n",
      "SCHZ     50\n",
      "BPLR     49\n",
      "ADHD     40\n",
      "Name: group, dtype: int64\n"
     ]
    }
   ],
   "source": [
    "# change directory\n",
    "os.chdir(data_dir)\n",
    "\n",
    "# get list of subjects in ciftdir\n",
    "subjects = glob.glob('sub*')\n",
    "output = np.zeros((len(subjects),len(TASKS)))\n",
    "\n",
    "for i,s in enumerate(subjects):\n",
    "    for j, task in enumerate(TASKS):\n",
    "        #check if the fully preprocessed data exists:\n",
    "        file = s + '/MNINonLinear/Results/task-' + task + '_bold/task-' + task + '_bold_Atlas_s0.dtseries.nii'\n",
    "        output[i,j] = os.path.isfile(file)\n",
    "        \n",
    "df = pd.DataFrame(data=output,index=subjects,columns=TASKS)\n",
    "\n",
    "# generate a group variable\n",
    "df['group'] = np.where(df.index.str.startswith('sub-1'),'CTRL',\n",
    "                       np.where(df.index.str.startswith('sub-5'),'SCHZ',\n",
    "                                np.where(df.index.str.startswith('sub-6'),'BPLR',\n",
    "                                         np.where(df.index.str.startswith('sub-7'),'ADHD',0))))\n",
    "\n",
    "print('Data after preprocessing:',len(df))\n",
    "display(df.head(5))\n",
    "print(df['group'].value_counts())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exclusions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gorgolewski et al., 2017 and missing data\n",
    "Note that other poor subjects from that paper have been excluded by way of preprocessing (no T1 files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data after missing data: 252\n",
      "CTRL    118\n",
      "SCHZ     48\n",
      "BPLR     48\n",
      "ADHD     38\n",
      "Name: group, dtype: int64\n"
     ]
    }
   ],
   "source": [
    "# Gorgolewski et al., 2017 \n",
    "df = df.drop('sub-11067') \n",
    "\n",
    "# remove subjects without all the data\n",
    "df['total'] = df.sum(axis=1)\n",
    "df = df[df['total'] == len(TASKS)]\n",
    "\n",
    "# subject has very odd patterns of response in BART task\n",
    "df = df.drop('sub-50010')\n",
    "df = df.drop('sub-70048') # odd QC (not in study cohort anyway)\n",
    "\n",
    "print('Data after missing data:',len(df))\n",
    "print(df['group'].value_counts())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## relative RMS (Satterthwaite et al. 2013)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dropped data due to RMS head motion\n",
      "SCHZ    12\n",
      "BPLR     7\n",
      "CTRL     6\n",
      "ADHD     3\n",
      "Name: group, dtype: int64\n",
      "Data after RMS: 224\n",
      "CTRL    112\n",
      "BPLR     41\n",
      "SCHZ     36\n",
      "ADHD     35\n",
      "Name: group, dtype: int64\n"
     ]
    }
   ],
   "source": [
    "# relative RMS cutoff in mm (Satterthwaite et al. 2013)\n",
    "RMSCutOff = 0.25\n",
    "\n",
    "# length cutoff (i.e., 20% of data)\n",
    "RMS_length_CutOff = 0.20\n",
    "\n",
    "# preallocate RMS array\n",
    "RMS_data = np.zeros((len(df.index),len(TASKS)+2))\n",
    "labels = []\n",
    "for j, task in enumerate(TASKS):\n",
    "    for s,subj in enumerate(df.index):  \n",
    "        # load the subjects confound data\n",
    "        df_confounds = pd.read_csv(FMRIPREP_DIR + subj + '/func/' + subj + '_task-' + task + '_bold_confounds.tsv', sep='\\t', header=0)\n",
    "\n",
    "        # get 6 motions parameters\n",
    "        motionParams = df_confounds[['X','Y','Z','RotX','RotY','RotZ']].values\n",
    "\n",
    "        # calculate relative RMS\n",
    "        displacement = np.zeros(np.shape(motionParams))\n",
    "        for tr in range(1,len(motionParams)):\n",
    "            for dim in range(0,3):\n",
    "                displacement[tr,dim] = motionParams[tr,dim] - motionParams[tr-1,dim];\n",
    "\n",
    "        relativeRMS = np.sqrt(displacement[:,0]**2 + displacement[:,1]**2+ displacement[:,2]**2)\n",
    "        \n",
    "        # save the data out\n",
    "        RMS_data[s,j] = np.mean(relativeRMS)\n",
    "        if (np.sum(relativeRMS > RMSCutOff) / float(len(relativeRMS))) > RMS_length_CutOff:\n",
    "            RMS_data[s,-1] = 1\n",
    "    labels.append('RMS-mean-'+task)\n",
    "labels.append('RMS-avg')\n",
    "RMS_data[:,-2] = np.mean(RMS_data[:,0:-3],axis=1)\n",
    "labels.append('RMS-exclude')\n",
    "\n",
    "# put the data in the df\n",
    "RMS_df = pd.DataFrame(data=RMS_data,index=df.index,columns=labels)\n",
    "df = pd.concat([df,RMS_df],axis=1)\n",
    "\n",
    "# show excluded subjects\n",
    "print('Dropped data due to RMS head motion')\n",
    "print(df['group'].loc[df['RMS-exclude']==1].value_counts())\n",
    "\n",
    "# drop the data and print\n",
    "df = df[df['RMS-exclude']==0]\n",
    "print('Data after RMS:',len(df))\n",
    "print(df['group'].value_counts())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic demographics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# basic between groups t-test based on the subject dataframe\n",
    "def CTRL_SCHZ_ttest(df,labels=['age','sex','site']):\n",
    "    for l in labels:\n",
    "        x = df[df['group']=='CTRL'][l].values\n",
    "        y = df[df['group']=='SCHZ'][l].values\n",
    "        t,p = ttest_ind(x,y,equal_var=False)\n",
    "        print('t-test,Var=',l)\n",
    "        print('\\tt=',np.round(t,3))\n",
    "        print('\\tp=',np.round(p,4))\n",
    "        \n",
    "def CTRL_SCHZ_stats(df,labels=['age','sex','site']):\n",
    "    \n",
    "    for l in labels:\n",
    "        if l == 'sex' or l == 'site':\n",
    "            new_df = df.loc[(df.group=='CTRL') | (df.group=='SCHZ')]\n",
    "            expected, observed, stats = pg.chi2_independence(new_df, x=l,y='group')\n",
    "            print('chi2 test,Var=',l)\n",
    "            print('\\tchi2=',np.round(stats.chi2[0],3))\n",
    "            print('\\tp=',np.round(stats.p[0],4))\n",
    "            \n",
    "        else:  #do a ttest\n",
    "            x = df[df['group']=='CTRL'][l].values\n",
    "            y = df[df['group']=='SCHZ'][l].values\n",
    "            t,p = ttest_ind(x,y,equal_var=False)\n",
    "            print('t-test,Var=',l)\n",
    "            print('\\tt=',np.round(t,3))\n",
    "            print('\\tp=',np.round(p,4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CTRL    112\n",
      "BPLR     41\n",
      "SCHZ     36\n",
      "ADHD     35\n",
      "Name: group, dtype: int64\n",
      "             age       sex      site\n",
      "group                               \n",
      "ADHD   31.371429  0.514286  0.457143\n",
      "BPLR   34.634146  0.512195  0.536585\n",
      "CTRL   31.392857  0.535714  0.812500\n",
      "SCHZ   35.500000  0.722222  0.472222\n",
      "t-test,Var= age\n",
      "\tt= -2.394\n",
      "\tp= 0.0199\n",
      "chi2 test,Var= sex\n",
      "\tt= 3.165\n",
      "\tp= 0.0753\n",
      "chi2 test,Var= site\n",
      "\tt= 14.316\n",
      "\tp= 0.0002\n"
     ]
    }
   ],
   "source": [
    "# .tsv with age/gender/site information\n",
    "df_participants = pd.read_csv('/projects/f_mc1689_1/ClinicalActFlow/data/ds000030_R105/participants.tsv', sep='\\t', header=0)\n",
    "\n",
    "# .tsv with other demographic info\n",
    "df_demo = pd.read_csv('/projects/f_mc1689_1/ClinicalActFlow/data/ds000030_R105/phenotype/phenotype_demographics.tsv', sep='\\t', header=0)\n",
    "df_wais = pd.read_csv('/projects/f_mc1689_1/ClinicalActFlow/data/ds000030_R105/phenotype/phenotype_wais.tsv', sep='\\t', header=0)\n",
    "\n",
    "#preallocate\n",
    "age = np.zeros((len(df.index)))\n",
    "sex = age.copy()\n",
    "site = age.copy()\n",
    "edu_years = age.copy()\n",
    "hand = age.copy()\n",
    "wais_mr = age.copy() \n",
    "wais_lns = age.copy()\n",
    "wais_voc = age.copy()\n",
    "\n",
    "for subidx,subj in enumerate(df.index):\n",
    "    idx = df_participants.participant_id==subj\n",
    "    \n",
    "    #age\n",
    "    age[subidx] = np.int(df_participants.age[idx].values)\n",
    "    \n",
    "    #sex\n",
    "    _sex = df_participants.gender[idx].values\n",
    "    if _sex == 'F':\n",
    "        sex[subidx] = 0\n",
    "\n",
    "    elif _sex == 'M':\n",
    "        sex[subidx] = 1\n",
    "        \n",
    "    #site\n",
    "    data = []\n",
    "    file = '/projects/f_mc1689_1/ClinicalActFlow/data/ds000030_R105/' + subj + '/func/' + subj + '_task-rest_bold.json'\n",
    "    with open(file) as json_data:\n",
    "        data.append(json.load(json_data))\n",
    "\n",
    "    if data[0]['SoftwareVersions'] == 'syngo MR B15':\n",
    "        site[subidx] = 1\n",
    "    elif data[0]['SoftwareVersions'] == 'syngo MR B17':\n",
    "        site[subidx] = 0\n",
    "    \n",
    "    # education (how many years of schooling complete?)\n",
    "    edu_years[subidx] = np.int(df_demo.school_yrs[idx].values)\n",
    "    \n",
    "    #\n",
    "    #hand = \n",
    "    # wais\n",
    "    wais_mr[subidx] = np.int(df_wais.mr_totalraw[idx].values) \n",
    "    wais_lns[subidx] = np.int(df_wais.lns_totalraw[idx].values)\n",
    "    wais_voc[subidx] = np.int(df_wais.voc_totalraw[idx].values)\n",
    "    \n",
    "# put the data in the df\n",
    "labels = ['age','sex','site','edu_years','wais_mr','wais_lns','wais_voc']\n",
    "demo_data = np.vstack((age,sex,site,edu_years,wais_mr,wais_lns,wais_voc)).T\n",
    "demo_df = pd.DataFrame(data=demo_data,index=df.index,columns=labels)\n",
    "df = pd.concat([df,demo_df],axis=1)\n",
    "\n",
    "# print average demos\n",
    "print(df['group'].value_counts())\n",
    "print(df.groupby(['group'])['age','sex','site'].mean())\n",
    "CTRL_SCHZ_stats(df,labels=['age','sex','site'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Manual HC cohort selection\n",
    "There are substantial differences in sex and age between HC and SZ, so we trim the HC cohort for a better sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CTRL    94\n",
      "BPLR    41\n",
      "SCHZ    36\n",
      "ADHD    35\n",
      "Name: group, dtype: int64\n",
      "             age       sex      site\n",
      "group                               \n",
      "ADHD   31.371429  0.514286  0.457143\n",
      "BPLR   34.634146  0.512195  0.536585\n",
      "CTRL   33.021277  0.627660  0.776596\n",
      "SCHZ   35.500000  0.722222  0.472222\n",
      "t-test,Var= age\n",
      "\tt= -1.419\n",
      "\tp= 0.161\n",
      "chi2 test,Var= sex\n",
      "\tt= 0.653\n",
      "\tp= 0.419\n",
      "chi2 test,Var= site\n",
      "\tt= 9.937\n",
      "\tp= 0.0016\n"
     ]
    }
   ],
   "source": [
    "for subidx,subj in enumerate(df.index[df['group']=='CTRL']):\n",
    "    subj_data = df[df.index==subj]\n",
    "    \n",
    "    age = subj_data['age'].values\n",
    "    sex = subj_data['sex'].values\n",
    "    site = subj_data['site'].values\n",
    "        \n",
    "    if age < 26:\n",
    "        if sex==0:\n",
    "            if site==1:\n",
    "                #print('dropping subj',subj)\n",
    "                df = df.drop(subj)\n",
    "    if age < 22:\n",
    "        if sex==1:\n",
    "            if site==1:\n",
    "                #print('dropping subj',subj)\n",
    "                df = df.drop(subj)\n",
    "\n",
    "# print average demos\n",
    "print(df['group'].value_counts())\n",
    "print(df.groupby(['group'])['age','sex','site'].mean())\n",
    "#CTRL_SCHZ_ttest(df,labels=['age','sex','site'])\n",
    "CTRL_SCHZ_stats(df,labels=['age','sex','site'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Behavioral outliers on the SCAP\n",
    "The analyses are focussed on the SCAP task, just want to make sure there are no odd participants.\n",
    "There are clear outliers based on the histogram of z-scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Poor subjs: Index(['sub-10788'], dtype='object')\n",
      "Odd HC participant scored: [31.25] %\n",
      "Odd HC participant scored: [-5.32460934] (zscore)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEBCAYAAAB2RW6SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAD5tJREFUeJzt3X+onuV9x/H38aibxOFxqZrq0WqL+W4VdYkVZTQycNp2NFirMDNtCltZrZsuVDYhZU03kDqNI7WJJK0UnctktY6oG1ZwIDM4oV0N+AO/2s64+KumsceSFWM4yf547tjHZ+d47ufHOfd5rvN+wc05z3X/uL4kJ59z5bp/jRw8eBBJUrkOa7oASdLsMuglqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9JhTPoJalwBr0kFe7wpjqOiF8DzgVeAyabqkOShswo8EHgB5m5r84OjQU9rZB/rMH+JWmYrQC219mwyaB/DWDr1q0sWbKkwTIkaXi8/vrrXHnllVBlaB1NBv0kwJIlSxgfH2+wDEkaSrWnvD0ZK0mFM+glqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SY16Z39zT0Bpsu+51OQNU5LEkUeMsvL6+xvp+8FbL2mk37nmiF6SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVrtbllRGxDTgNOADsBa7NzB0RsRN4u1oAbsjMh2ehTklSj+peR//5zHwLICIuAb4DLK/WXZ6ZT89GcZKk/tUK+kMhXzmG1si+togYA8Y6mn2tlCTNgdp3xkbEHcDFwAjwybZVWyNihNZLatdm5sQUu68B1vVTqCSpN7VPxmbmFzLzFGAtcEvVvCIzzwbOpfULYOM0u2+gNcffvqzotWhJUn1dP+smM++OiG9FxOLM3FW17YuI24EHptlnAnjPSD8ieqlXktSlGUf0EXF0RJzc9nkl8CbwdkQcU7WNAFcAO2arUElSb+qM6BcB90bEImCSVsivBE4A7ouIUWAUeBa4ZrYKlST1Zsagz8yfAudPs3rZYMuRJA2ad8ZKUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9JhavzcnAiYhtwGnAA2Atcm5k7ImIpcBewGNgDrM7MF2arWElS9+qO6D+fmWdn5jJgPfCdqn0zsCkzlwKbgC2zUKMkqQ+1RvSZ+Vbbx2OAAxFxPLAcuKhqvwfYGBHHZebu9v0jYgwY6zjseG8lS5K6USvoASLiDuBiYAT4JHAy8EpmTgJk5mREvFq17+7YfQ2wbiAVS5K6UvtkbGZ+ITNPAdYCt3TZzwZac/zty4oujyFJ6kHtEf0hmXl3RHwLeBk4KSJGq9H8KHAisGuKfSaAifa2iOixZElSN2Yc0UfE0RFxctvnlcCbwBvADmBVtWoV8GTn/LwkqVl1RvSLgHsjYhEwSSvkV2bmwYi4GrgrIr4K/BxYPXulSpJ6MWPQZ+ZPgfOnWfcccN6gi5IkDY53xkpS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFm/GdsRGxGLgb+AiwD/gx8MXM3B0RB4GngAPV5p/LzKdmq1hJUvdmDHrgIHBzZj4KEBG3ADcBf1Kt/93M3Ds75UmS+jVj0Gfmm8CjbU1PAF/qppOIGAPGOprHuzmGJKk3dUb074qIw2iF/ANtzY9GxOHAQ8DXMnPfFLuuAdb1XKUkqWfdnoz9JrAX2Fh9PiUzPwZcAHwU+Otp9tsAnNaxrOi6WklS12qP6CNiPXA6sDIzDwBk5q7q6y8i4g7gy1Ptm5kTwETH8XqtWZLUhVoj+oi4ETgH+MyhqZmIODYijqq+Pxy4HNgxW4VKknpT5/LKM4C1wPPA49VI/EXgZmBLdYnlEcDjTD91I0lqSJ2rbp4BRqZZfdZgy5EkDZp3xkpS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFq/Ny8MXA3cBHgH3Aj4EvZubuiDgf2AIcBewErsrMN2avXElSt+qM6A8CN2dmZOZZwE+AmyJiBPhH4M8ycynwH8BNs1eqJKkXM47oM/NN4NG2pieALwEfA97OzO1V+2Zao/o/7jxGRIwBYx3N492XK0nqVldz9BFxGK2QfwA4BXjp0LrM/BlwWET85hS7rgFe7Fge67FmSbPgnf2TTZegWTLjiL7DN4G9wEbg0i722wDc2dE2jmEvzRtHHjHKyuvvn/N+H7z1kjnvc6GpHfQRsR44HViZmQci4n+AD7Wt/wBwsJrqeY/MnAAmOo7Xc9GSpPpqTd1ExI3AOcBnMnNf1fxfwFER8fHq89XAdwdfoiSpH3UurzwDWAs8DzxejcRfzMxLI+JzwJaI+HWqyytnsVZJUg/qXHXzDDAyzbrHgTMHXZQkaXC8M1aSCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXtGA19Xyfue6322fdSFIxFsrzfRzRS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9Jhav1rJuIWA9cBpwKnJmZT1ftO4G3qwXghsx8eOBVSpJ6VvehZtuAbwCPTbHu8kPBL0maf2oFfWZuB4iInjqJiDFgrKN5vKeDSZK6MojHFG+NiBFgO7A2Myem2GYNsG4AfUmSutTvydgVmXk2cC4wAmycZrsNwGkdy4o++5Yk1dDXiD4zd1Vf90XE7cAD02w3AbxnpN/rNJAkqTs9j+gjYlFEHFN9PwJcAewYVGGSpMGoe3nlbcBngSXAIxGxB1gJ3BcRo8Ao8CxwzWwVKknqTd2rbq4Drpti1bLBliNJGjTvjJWkwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKN+M7YyNiPXAZcCpwZmY+XbUvBe4CFgN7gNWZ+cLslSpJ6kWdEf024ALgpY72zcCmzFwKbAK2DLg2SdIAzDiiz8ztABHxbltEHA8sBy6qmu4BNkbEcZm5u/MYETEGjHU0j/dYsySpC73O0Z8MvJKZkwDV11er9qmsAV7sWB7rsW9JUhdmHNEPyAbgzo62cQx7SZp1vQb9LuCkiBjNzMmIGAVOrNr/n8ycACba29qngiRJs6enqZvMfAPYAayqmlYBT041Py9JataMQR8Rt0XEy7SmWh6JiGeqVVcD10bE88C11WdJ0jxT56qb64Drpmh/DjhvNoqSJA2Od8ZK88w7+yebLkGFmaurbiTVdOQRo6y8/v457/fBWy+Z8z41NxzRS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpqCL/9QSXzxiDSFpl7+Ab4ARIPXd9BHxE7g7WoBuCEzH+73uJKkwRjUiP7yzHx6QMeSJA3QnEzdRMQYMNbRPD4XfUvSQjeooN8aESPAdmBtZk50rF8DrBtQX5KkLgziqpsVmXk2cC4wAmycYpsNwGkdy4oB9C1JmkHfI/rM3FV93RcRtwMPTLHNBPCeUX5E9Nu1JKmGvkb0EbEoIo6pvh8BrgB2DKIwSdJg9DuiPwG4LyJGgVHgWeCavquSJA1MX0Gfmf8NLBtQLZKkWeAjECSpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHrNa74AROqfLx7RvNbUC0B8+YdK4ohekgpn0EtS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKN7RB39St8d6SL2nYDO0jELw1XpLq6TvoI2IpcBewGNgDrM7MF/o9riRpMAYxdbMZ2JSZS4FNwJYBHFOSNCB9jegj4nhgOXBR1XQPsDEijsvM3W3bjQFjHbt/COD111/vuf/9v3yz53179fLLL895nwtdU3/PTfTbZN8Lrd8m++4nR9oyc7TuPiMHDx7sucOIOAf4h8w8o63tWeCqzPxRW9vXgHU9dyRJ6rQiM7fX2XCuTsZuAO7saDsS+DDwAjCXl7KMA48BK4BhGZ5b89wYxpphOOu25t6NAh8EflB3h36DfhdwUkSMZuZkRIwCJ1bt78rMCWBiiv2f77P/rkXEoW9fzsydc91/L6x5bgxjzTCcdVtz337SzcZ9nYzNzDeAHcCqqmkV8GT7/LwkqVmDmLq5GrgrIr4K/BxYPYBjSpIGpO+gz8zngPMGUIskaRYM7SMQ+jAB/A1TnzOYr6x5bgxjzTCcdVvzHOrr8kpJ0vy3EEf0krSgGPSSVLihfXplvyLiTuD3gZ9VTfdm5o3NVVRfRPwe8O/AX2TmxobLeV8R8RXgD2ndFDcCfD0z/7nZqmYWEZuAC4F9wF5af9Y/bLaq9xcRVwF/BXwUWDNffzaG8UGIEbEeuAw4FTgzM59utqLuLPQR/U2Z+TvVMiwh/xvA3wEPNV1LTRsz86zMXAb8AfDtiDi26aJqeIjWP+izga8D8/6XE617Wq4A/qnpQmYwjA9C3AZcALzUdCG9WOhBP4z+HriFX/1PZF7LzLfaPh4NHGQIfu4y818zc3/18T+B8YiY13Vn5tOZ+SxwoOlaptP2IMR7qqZ7gOURcVxzVc0sM7dn5q6Zt5yf5vUP7hz4ckQ8FRHbIuK3my5mJhHxKWAsM7/XdC3diIirI+I54EngTzNzT9M1denPgX/LzHkboEPkZOCVzJwEqL6+WrVrlhQ7Rx8RPwJOmWb1CcBXgNcy80BErAa+HxEfPvQD2IQZag7gJn71SOh5YaY/58yczMzNwOaIOBPYGhGPNB32dequtrsC+CNa/21vVN2apU5eR1+JiD3A8sycl3NwEfFx4F+AX1ZNH6B1ovAbmfm3jRXWpYj4PvDtzLyv6VpmEhGXAuuBC+fBQ6xqqy40+OF8PBlbTd08DyxuexDiHuD0YXhGVkTsBD7tydghEREntX3/CVpXhbzSXEXvr5ojPD4zT83MU4HvAevme8i3T4lFxGnAMuDZ5iqqJyI+Tet8yCeGKeTnOx+E2IwFO6KPiEdoTeEcAH4B/GVmPtFsVfXN51Fbu4j4LnAGsJ/WL9Obh+Tyyt3AO0B7AF3Y9JTT+4mIVbRO1B9Lq/b/BS6uTtDOGxHxW7QurzyW6kGImZnNVvX+IuI24LPAEloXQuxpf+HSfLdgg16SFooFO3UjSQuFQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuH+D40iUZHuKak0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEBCAYAAACKUEVYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAADi5JREFUeJzt3X+MZfVZx/H3MuzWatWtDe3SDj9qyz6lVKUgWtLQYipVA5stagXKD9s/kKJESbZp1Wh3bYMidttNhQQtUaAiQbRhWaQVSbCyWBsIkAbRpxaXFgpb+eGoILLTYfzj3rGXyezec2fvOXd49v1KJnvnzPfO9zNnZz7znXPPPXfV/Pw8kqSXtoMmHUCStP8sc0kqwDKXpAIsc0kqwDKXpAIsc0kqwDKXpAIsc0kqwDKXpAIsc0kqwDKXpAIObvOTR8TLgBOAx4G5NueSpEKmgEOBuzPz+SZ3aLXM6RX5nS3PIUlVnQTsbDKw7TJ/HOC6665j3bp1LU8lSTXs3r2bs88+G/od2kTbZT4HsG7dOqanp1ueSpLKaXx42gdAJakAy1ySCrDMJakAy1ySCrDMJakAy1ySCrDMJakAy1yagD2zk7m6xaTmVfvaftKQpCWsWT3Fhk3bO593x9aNnc+pbrgyl6QCLHNJKsAyl6QCLHNJKsAyl6QCLHNJKsAyl6QCLHNJKsAyl6QCLHNJKsAyl6QCLHNJKsAyl6QCLHNJKsAyl6QCLHNJKsAyl6QCLHNJKqDRy8ZFxGnAx4FV9H4BbMnMz7UZTJLU3NCVeUSsAj4LnJuZxwLnANdEhKt6SVohmr6g8wvA9/dvrwUez8wXBgdExNr+xwZN7188SVITQ8s8M+cj4heA7RHxLPC9wKlLDL0Y2DzmfJK0X/bMzrFm9VT5eYeWeUQcDPwGsDEz74qItwM3RMSbM/OZgaHbgKsX3X0auHNcYSVpVGtWT7Fh0/bO592xdWOn8zU5zHIs8NrMvAugX+jPAkcDdy8MyswZYGbwjhExxqiSpL1p8iDmo8B09Js5Io4G1gEPtRlMktRck2PmuyPiQuAvI2LhQc8PZObT7UaTJDXV6GyWzLwOuK7lLJKkZfJccUkqwDKXpAIsc0kqwDKXpAIsc0kqwDKXpAIsc0kqwDKXpAIsc0kqwDKXpAIsc0kqwDKXpAIsc0kqwDKXpAIsc0kqwDKXpAIsc0kqwDKXpAIsc0kqwDKXpAIsc0kqwDKXpAIsc0kqwDKXpAIsc0kqwDKXpAIsc0kqwDKXpAIsc0kqwDKXpAIsc0kqwDKXpAIsc0kqwDKXpAIsc0kqwDKXpAIsc0kq4OAmgyLiu4BPAT8J/C/wpcz8pTaDSZKaa1TmwGX0Snx9Zs5HxGtazCRJGtHQMo+IVwDnAdOZOQ+Qmd9aYtxaYO2izdPjCClJ2rcmK/M3AE8BmyPiJ4BngN/KzJ2Lxl0MbB5zPklSA00eAD0Y+EHgvsz8UeAjwOci4vsWjdsGvH7R20ljzCpJ2osmK/OvA98GrgfIzC9HxJPAeuCehUGZOQPMDN4xIsaXVJK0V0NX5pn5JHAHcApARKwHXg18rd1okqSmmp7N8kHgTyJiKzALnNtfiUuSVoBGZZ6Z/wac3G4USdJy+QxQSSrAMpekAixzSSrAMpekAixzSSrAMpekAixzSSrAMpekAixzSSrAMpekAixzSSrAMpekAixzSSrAMpekAixzSSrAMpekAixzSSrAMpfUuj2zc5OOUF7T1wCVpGVbs3qKDZu2T2TuHVs3TmTerrkyl6QCLHNJKsAyl6QCLHNJKsAyl6QCLHNJKsAyl6QCLHNJKsAyl6QCLHNJKsAyl6QCLHNJKsAyl6QCLHNJKsAyl6QCLHNJKsAyl6QCLHNJKmCkMo+IzRExHxFvaSuQJGl0jcs8Io4D3gZ8o704kqTlaPSCzhHxMuAK4H3AHXsZsxZYu2jz9H6lkyQ10nRl/jHgzzJz1z7GXAzsWvR25/7Fk9qzZ3Zu0hGksRm6Mo+IE4ETgF8fMnQbcPWibdNY6Fqh1qyeYsOm7ROZe8fWjROZV3U1OczyTuBNwK6IgF5B/01EfCAzb1sYlJkzwMzgHfvjJUktG1rmmXkpcOnC+xHxMHBaZj7QXixJ0ig8z1ySCmh0NsugzDyyhRySpP3gylySCrDMJakAy1ySCrDMJakAy1ySCrDMJakAy1ySCrDMJakAy1ySCrDMJakAy1ySCrDMJakAy1ySCrDMJakAy1ySCrDMJakAy1ySCljRZb5ndu6Amldqm9/bdY38snFdWrN6ig2btnc+746tGzufU+qCP1N1reiVuSSpGctckgqwzCWpAMtckgqwzCWpAMtckgqwzCWpAMtckgqwzCWpAMtckgqwzCWpAMtckgqwzCWpAMtckgqwzCWpAMtckgqwzCWpAMtckgoY+rJxEfEq4LPAG4Dnga8BF2TmEy1nkyQ11GRlPg9clpmRmT8MPARc2m4sSdIohq7MM/Np4O8GNv0jcOHicRGxFli7aPP0/oSTJDUztMwHRcRB9Ir85iU+fDGweRyhDmR7ZudYs3rqgJlX0niMVObAHwLPAJcv8bFtwNWLtk0Dd44e68C1ZvUUGzZt73zeHVs3dj6npPFpXOYR8QngKGBDZr6w+OOZOQPMLLrPfgeUJA3XqMwj4hLgeODUzHy+3UiSpFE1OTXxGOA3ga8C/9Bfbe/KzNNbziZJaqjJ2Sz/BKzqIIskaZl8BqgkFWCZS1IBlrkkFWCZS1IBlrkkFWCZS1IBlrkkFWCZS1IBlrkkFWCZS1IBlrkkFWCZS1IBlrkkFWCZS1IBlrkkFWCZS1IBlrkkFWCZL2HP7NykI3Rukl/zgbi/pXFr9ILOB5o1q6fYsGn7RObesXXjROY9EL9mqRJX5pJUgGUuSQVY5pJUgGUuSQVY5pJUgGUuSQVY5pJUgGUuSQVY5pJUgGUuSQVY5pJUgGUuSQVY5pJUgGUuSQVY5pJUgGUuSQVY5pJUgGUuSQU0etm4iFgPXAO8CngKOC8z/7XNYJKk5pquzK8ErsjM9cAVwB+1F0mSNKqhK/OIeDVwHHBKf9P1wOURcUhmPjEwbi2wdtHdjwDYvXv3sgPO/s/Ty77vcj366KMTmXeSc0/6a54U/59rzzvJuffn+3qgM6ea3mfV/Pz8PgdExPHAtZl5zMC2B4FzMvPegW1bgM0j5JUk7dtJmbmzycBGx8wb2gZcPfD+NHAn8H5gJzA3xrnasJD3JGByS8VmzNoOs7bDrKObAg4F7m56hyZl/gjwuoiYysy5iJgCXtvf/v8ycwaYWXg/IhZufjEzH24aaFIG8j660vOatR1mbYdZl+2hUQYPfQA0M/8duB84q7/pLOC+wePlkqTJanqY5YPANRHxUeA/gPPaiyRJGlWjMs/MfwF+vOUskqRlavMZoDPA7zBwHH2FeynlNWs7zNoOs3Zg6KmJkqSVz2uzSFIBlrkkFTDOJw0REVcA7wKeB54Bfi0z71li3PvpPcno4f6mXZl5+jizDNM0a3/sb9N78hPA1Zn58U5CvjjDOcCHgTcDF2fm5XsZdzJwK/DV/qbnM7PTB6+bZu2PPR/4CLAK+Dzwq5n5QidBe/N/N/CnwPHAt4EPZeYtS4w7mQns1yYXues/9+PTwE8D88ClmXlV29mWmXUL8MvAY/1Nd2Xmr3SZs5/jE8DPAUcCP5SZDywxZkXs16bGvTL/PL0d8yPA7wE37GPs7Zl5bP+t0yLva5Q1It4BvBd4S//tvf1tXbsfOBP48wZjHxzYt5M4C6lR1oh4Pb1LQJwIHNV/O6f1dC/2IeC/M/ONwAbgqoh4xV7GTmK/NrnI3dnAG+ntvxOBLRFxZEf5BjW9IN+1A/ux8yLvuwl4B/D1fYxZKfu1kbGWeWbekpmz/Xe/BExHxIo8lDNC1jPoffM9l5nPAdf2t3UqMx/IzAeBzlatyzVC1p8HbsrMJ/qr8c/Q/b49g14J0V9F3gP8TMcZljRwkbvr+5uuB46LiEMWDT0D+ExmvtB/Mt9N9BYgnRkh64qQmTsz85Ehwya+X0fRZtFeBPz1Pv5kfmdE3B8Rfx8Rp7aYo4l9ZT2cF//2/gZwWCeplm99RNwbEV+OiF+cdJh9WAn7dpQMXe/Xw4BvZuYcQP/fx5bItxL2Y9OsAGdGxFci4raIOLHLkCNaCfu1sZGOmUfEvfS+wKW8ZuE/MiLOBN5H78+YpdwC3JCZz0XEW4EvRMTJmfnPo+TpKGsnmuZt4F7gsMz8z/5hjNsj4puZeftYgjLWrK0blnWET9X6fj1AXAlckpmzEXEKsD0ijs7MpyYd7KVupDLPzOOGjYmI04FLgHdl5rf28nmeHLh9X0TsBH4MGFuZjysrvd/GRwy8fziLLjI2Dk3yNvw8/zVwe1dE3AS8HRhb6YwrKx3s22FZI2Ihw8K1hg4H7lji87S+X5fQ6CJ3fGc/Llxhb/GKsgtNL8i3e+D230bEI/Qei/pip2mbWQn7tbGxHmaJiNOATwI/ta8rjkXE6wZuHwG8DfjKOLMM0zQrcCNwXkS8PCJeTu+6NH/RQcRliYhDI2JV//YPAO+m94DkSvRXwHsi4pD+4xXn0/2+vRG4ACAijgJOAL6weNAk9usIF7m7ETg/Ig7qH6N+D71925mmWRf97B9L72yS7CjmqCa+X0cx1meARsQTwB6+s8qB3qr3qYi4Crg5M2+OiN8FNtI7FQzgk5l5zdiCjDFrf+wW4Fx6p89dm5lbuszaz3AW8AfAK+nlfhZ4d2Y+GBEfAx7LzCsj4iLgQmCW3l9e12bmZSsxa3/sBfROYwS4Dbioy8M0EfE99K7D/1Z619z/cGZu739s4vs1It5E73S/V9K/yF1mZkTcCnw0M+/pr4Ivp/cLBuD3M/OP2862zKzX0DsNdI7e98bmzLx1Alk/DfwssA54EngqM49Zifu1KZ/OL0kFrMjTBiVJo7HMJakAy1ySCrDMJakAy1ySCrDMJakAy1ySCrDMJamA/wPfk5XjrlNATgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CTRL    93\n",
      "BPLR    41\n",
      "SCHZ    36\n",
      "ADHD    35\n",
      "Name: group, dtype: int64\n",
      "             age       sex      site\n",
      "group                               \n",
      "ADHD   31.371429  0.514286  0.457143\n",
      "BPLR   34.634146  0.512195  0.536585\n",
      "CTRL   33.000000  0.634409  0.774194\n",
      "SCHZ   35.500000  0.722222  0.472222\n",
      "t-test,Var= age\n",
      "\tt= -1.427\n",
      "\tp= 0.1586\n",
      "chi2 test,Var= sex\n",
      "\tt= 0.543\n",
      "\tp= 0.4613\n",
      "chi2 test,Var= site\n",
      "\tt= 9.696\n",
      "\tp= 0.0018\n"
     ]
    }
   ],
   "source": [
    "from functions import load_scan_behaviour\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import zscore\n",
    "accuracy,reaction_time,_,_,task_labels = load_scan_behaviour(df.index,['scap'],df['group'].values)\n",
    "hc_acc = np.mean(np.mean(accuracy['scap']['CTRL'],axis=0),axis=0)\n",
    "sz_acc = np.mean(np.mean(accuracy['scap']['SCHZ'],axis=0),axis=0)\n",
    "\n",
    "idx = zscore(hc_acc) < -4\n",
    "bad_subjs = df.index[df['group']=='CTRL'][idx]\n",
    "print('Poor subjs:',bad_subjs)\n",
    "print('Odd HC participant scored:',hc_acc[idx]*100,'%')\n",
    "print('Odd HC participant scored:',zscore(hc_acc)[idx],'(zscore)')\n",
    "\n",
    "#plot the histograms\n",
    "plt.hist(zscore(hc_acc))\n",
    "plt.show()\n",
    "plt.hist(zscore(sz_acc))\n",
    "plt.show()\n",
    "\n",
    "# # drop the subject\n",
    "df = df.drop(bad_subjs)\n",
    "\n",
    "# print average demos\n",
    "print(df['group'].value_counts())\n",
    "print(df.groupby(['group'])['age','sex','site'].mean())\n",
    "#CTRL_SCHZ_ttest(df,labels=['age','sex','site'])\n",
    "CTRL_SCHZ_stats(df,labels=['age','sex','site'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save the subject list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write text out\n",
    "os.chdir(outdir)\n",
    "df = df.sort_index()\n",
    "df.to_csv(outdir+'subject_list.txt', sep='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "subj_list = list(df.index)\n",
    "with open('batch_subj_list.txt', 'w') as f:\n",
    "    for item in subj_list:\n",
    "        f.write(\"%s\\n\" % item)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demographics table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CTRL    93\n",
      "BPLR    41\n",
      "SCHZ    36\n",
      "ADHD    35\n",
      "Name: group, dtype: int64\n",
      "             age       sex      site   RMS-avg  edu_years    wais_mr  \\\n",
      "group                                                                  \n",
      "ADHD   31.371429  0.514286  0.457143  0.055353  14.571429  20.628571   \n",
      "BPLR   34.634146  0.512195  0.536585  0.070833  14.682927  19.585366   \n",
      "CTRL   33.000000  0.634409  0.774194  0.060001  15.161290  20.430108   \n",
      "SCHZ   35.500000  0.722222  0.472222  0.085636  12.777778  15.777778   \n",
      "\n",
      "        wais_lns   wais_voc  \n",
      "group                        \n",
      "ADHD   19.885714  42.914286  \n",
      "BPLR   19.707317  43.365854  \n",
      "CTRL   21.053763  43.483871  \n",
      "SCHZ   17.750000  32.000000  \n",
      "t-test,Var= age\n",
      "\tchi2= -1.427\n",
      "\tp= 0.1586\n",
      "chi2 test,Var= sex\n",
      "\tt= 0.543\n",
      "\tp= 0.4613\n",
      "chi2 test,Var= site\n",
      "\tt= 9.696\n",
      "\tp= 0.0018\n",
      "t-test,Var= RMS-avg\n",
      "\tchi2= -3.952\n",
      "\tp= 0.0002\n",
      "t-test,Var= edu_years\n",
      "\tchi2= 8.266\n",
      "\tp= 0.0\n",
      "t-test,Var= wais_mr\n",
      "\tchi2= 5.1\n",
      "\tp= 0.0\n",
      "t-test,Var= wais_lns\n",
      "\tchi2= 4.883\n",
      "\tp= 0.0\n",
      "t-test,Var= wais_voc\n",
      "\tchi2= 6.498\n",
      "\tp= 0.0\n"
     ]
    }
   ],
   "source": [
    "print(df['group'].value_counts())\n",
    "labels = ['age','sex','site','RMS-avg','edu_years','wais_mr','wais_lns','wais_voc']\n",
    "print(df.groupby(['group'])[labels].mean())\n",
    "#CTRL_SCHZ_ttest(df,labels=labels)\n",
    "CTRL_SCHZ_stats(df,labels=labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Supplementary: exclude based on site"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dropping subj sub-10171\n",
      "dropping subj sub-10958\n",
      "dropping subj sub-10746\n",
      "dropping subj sub-10329\n",
      "dropping subj sub-10471\n",
      "dropping subj sub-10844\n",
      "dropping subj sub-10506\n",
      "dropping subj sub-10280\n",
      "dropping subj sub-10624\n",
      "dropping subj sub-10686\n",
      "dropping subj sub-10877\n",
      "dropping subj sub-10438\n",
      "dropping subj sub-10708\n",
      "dropping subj sub-10159\n",
      "dropping subj sub-10692\n",
      "dropping subj sub-10530\n",
      "dropping subj sub-10525\n",
      "dropping subj sub-10527\n",
      "dropping subj sub-10523\n",
      "dropping subj sub-10575\n",
      "dropping subj sub-10227\n",
      "dropping subj sub-10707\n",
      "dropping subj sub-10557\n",
      "dropping subj sub-10882\n",
      "dropping subj sub-10325\n",
      "dropping subj sub-10440\n",
      "dropping subj sub-10940\n",
      "dropping subj sub-11052\n",
      "dropping subj sub-10949\n",
      "dropping subj sub-10968\n",
      "dropping subj sub-10987\n",
      "dropping subj sub-10235\n",
      "dropping subj sub-10249\n",
      "dropping subj sub-10273\n",
      "CTRL    59\n",
      "BPLR    41\n",
      "SCHZ    36\n",
      "ADHD    35\n",
      "Name: group, dtype: int64\n",
      "             age       sex      site\n",
      "group                               \n",
      "ADHD   31.371429  0.514286  0.457143\n",
      "BPLR   34.634146  0.512195  0.536585\n",
      "CTRL   37.084746  0.593220  0.644068\n",
      "SCHZ   35.500000  0.722222  0.472222\n",
      "t-test,Var= age\n",
      "\tt= 0.853\n",
      "\tp= 0.3964\n",
      "chi2 test,Var= sex\n",
      "\tchi2= 1.106\n",
      "\tp= 0.2929\n",
      "chi2 test,Var= site\n",
      "\tchi2= 2.049\n",
      "\tp= 0.1523\n"
     ]
    }
   ],
   "source": [
    "for subidx,subj in enumerate(df.index[df['group']=='CTRL']):\n",
    "    subj_data = df[df.index==subj]\n",
    "    \n",
    "    age = subj_data['age'].values\n",
    "    sex = subj_data['sex'].values\n",
    "    site = subj_data['site'].values\n",
    "        \n",
    "    if age < 32:\n",
    "        if sex==0:\n",
    "            if site==1:\n",
    "                print('dropping subj',subj)\n",
    "                df = df.drop(subj)\n",
    "    if age < 29:\n",
    "        if sex==1:\n",
    "            if site==1:\n",
    "                print('dropping subj',subj)\n",
    "                df = df.drop(subj)\n",
    "\n",
    "# print average demos\n",
    "print(df['group'].value_counts())\n",
    "print(df.groupby(['group'])['age','sex','site'].mean())\n",
    "#CTRL_SCHZ_ttest(df,labels=['age','sex','site'])\n",
    "CTRL_SCHZ_stats(df,labels=['age','sex','site'])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save site controlled subject list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write text out\n",
    "os.chdir(outdir)\n",
    "df = df.sort_index()\n",
    "df.to_csv(outdir+'subject_list_siteControl.txt', sep='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Site controlled demographics table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CTRL    59\n",
      "BPLR    41\n",
      "SCHZ    36\n",
      "ADHD    35\n",
      "Name: group, dtype: int64\n",
      "             age       sex      site   RMS-avg  edu_years    wais_mr  \\\n",
      "group                                                                  \n",
      "ADHD   31.371429  0.514286  0.457143  0.055353  14.571429  20.628571   \n",
      "BPLR   34.634146  0.512195  0.536585  0.070833  14.682927  19.585366   \n",
      "CTRL   37.084746  0.593220  0.644068  0.062943  15.118644  19.474576   \n",
      "SCHZ   35.500000  0.722222  0.472222  0.085636  12.777778  15.777778   \n",
      "\n",
      "        wais_lns   wais_voc  \n",
      "group                        \n",
      "ADHD   19.885714  42.914286  \n",
      "BPLR   19.707317  43.365854  \n",
      "CTRL   20.576271  43.440678  \n",
      "SCHZ   17.750000  32.000000  \n",
      "t-test,Var= age\n",
      "\tt= 0.853\n",
      "\tp= 0.3964\n",
      "t-test,Var= sex\n",
      "\tt= -1.297\n",
      "\tp= 0.1984\n",
      "t-test,Var= site\n",
      "\tt= 1.633\n",
      "\tp= 0.1069\n",
      "t-test,Var= RMS-avg\n",
      "\tt= -3.277\n",
      "\tp= 0.0017\n",
      "t-test,Var= edu_years\n",
      "\tt= 7.31\n",
      "\tp= 0.0\n",
      "t-test,Var= wais_mr\n",
      "\tt= 3.644\n",
      "\tp= 0.0005\n",
      "t-test,Var= wais_lns\n",
      "\tt= 4.011\n",
      "\tp= 0.0002\n",
      "t-test,Var= wais_voc\n",
      "\tt= 6.056\n",
      "\tp= 0.0\n"
     ]
    }
   ],
   "source": [
    "print(df['group'].value_counts())\n",
    "labels = ['age','sex','site','RMS-avg','edu_years','wais_mr','wais_lns','wais_voc']\n",
    "print(df.groupby(['group'])[labels].mean())\n",
    "CTRL_SCHZ_ttest(df,labels=labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}