code/brady_2020_classification_analyses.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "### load packages\n",
    "import itertools\n",
    "from itertools import compress\n",
    "import pdb\n",
    "\n",
    "import time\n",
    "from datetime import timedelta, datetime\n",
    "from os.path import join\n",
    "import os\n",
    "from os import listdir\n",
    "from os.path import isfile, join\n",
    "import csv\n",
    "from timeit import default_timer as timer\n",
    "\n",
    "import numpy as np\n",
    "import numpy.matlib\n",
    "import pandas as pd\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import figure\n",
    "from matplotlib.ticker import PercentFormatter\n",
    "%matplotlib inline\n",
    "\n",
    "from mne.filter import filter_data\n",
    "\n",
    "import scipy\n",
    "from scipy.signal import welch, decimate, periodogram, find_peaks\n",
    "from scipy import signal\n",
    "from scipy import stats\n",
    "from scipy.stats import pearsonr, mannwhitneyu, spearmanr, ranksums, ttest_ind\n",
    "from statsmodels.stats.multitest import multipletests\n",
    "from scipy.ndimage.filters import uniform_filter1d\n",
    "from statsmodels.stats.anova import AnovaRM\n",
    "from statistics import mode\n",
    "import math  \n",
    "import random \n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.preprocessing import StandardScaler \n",
    "from sklearn.kernel_ridge import KernelRidge\n",
    "from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier, GradientBoostingClassifier\n",
    "from sklearn import tree\n",
    "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis\n",
    "from sklearn.linear_model import LogisticRegression, LinearRegression\n",
    "from sklearn.svm import LinearSVC, SVC\n",
    "from sklearn.model_selection import LeaveOneGroupOut, cross_val_predict, cross_val_score, KFold, train_test_split, StratifiedKFold, GridSearchCV \n",
    "from sklearn.model_selection import permutation_test_score\n",
    "from sklearn.metrics import roc_auc_score, roc_curve, classification_report, confusion_matrix, accuracy_score\n",
    "from sklearn.metrics import precision_score, recall_score, auc, f1_score, precision_recall_curve\n",
    "\n",
    "import sys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pandas: 1.1.3\n",
      "Numpy: 1.19.1\n",
      "Scipy: 1.5.2\n"
     ]
    }
   ],
   "source": [
    "print('Pandas:',pd.__version__)\n",
    "print('Numpy:',np.__version__)\n",
    "print('Scipy:',scipy.__version__)\n",
    "# Pandas: 1.1.3\n",
    "# Numpy: 1.19.1\n",
    "# Scipy: 1.5.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notebook to Perform Supervised Classification Analysis\n",
    "\n",
    "<p> <strong>Part 1:</strong> Loads in extracted features in .csv files per patient. Uses clinical scores from updrs .csv table to select patients to include in analysis based on clinical fluctuations pre vs post medication. </p>\n",
    "<p> <strong>Part 2:</strong> Executes classification analysis.  </p>\n",
    "<p> <strong>Part 3:</strong> Executes classification analysis. Different functions for individual vs group models, optional to include activity filter to remove inactive feature-windows. Seperate function for true predictions and permutation creation. Significancy testing in seperate Results notebook.  </p>\n",
    "<p> <strong>Part 4:</strong> Sub analysis analyzing the influence of the number of patients used as training data in the group model.  </p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Path and working directory defined as /Users/jeroenhabets/Research/pd_sensors/brady_reallife\n"
     ]
    }
   ],
   "source": [
    "# define path, working directory which contains filtered accelerometer files\n",
    "notebook_path = os.getcwd()\n",
    "path = os.path.dirname(notebook_path)\n",
    "images_path = \"/Users/jeroenhabets/Starr Lab Dropbox/jeroen habets/PHD werkmap/UCSF time/Sensor-project/analysis images/\"\n",
    "os.chdir(path)\n",
    "print('Path and working directory defined as %s' %path)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: Feature Loading\n",
    "Previsouly filtered data is loaded in to dataframes per patient-side."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loadUPDRSscores():\n",
    "    '''\n",
    "    Input:\n",
    "    Load in existing updrsScores table\n",
    "    \n",
    "    Returns:\n",
    "    - updrsScores: dataframe with clinical scores of all patients\n",
    "    \n",
    "    '''\n",
    "    # Read in updrs file\n",
    "    updrsScores = pd.read_csv(os.path.join(path,'data','updrsScores.csv'))\n",
    "    # convert int to strings for PtId's\n",
    "    updrsScores['PtId'] = updrsScores['PtId'].astype(str)\n",
    "    # add missing zeros in front of PtId\n",
    "    for row in np.arange(updrsScores.shape[0]):\n",
    "        updrsScores.at[row,'PtId'] = updrsScores.loc[row]['Record Id'][3:]\n",
    "    del(updrsScores['Record Id'])\n",
    "    updrsScores = updrsScores.set_index('PtId')\n",
    "    \n",
    "    # create list for ON and OFF for subscore-lists\n",
    "    OFF_sublists = {\n",
    "        'leftHandBradyOFF': ['OFF_UPDRS_3_3c','OFF_UPDRS_3_4b','OFF_UPDRS_3_5b','OFF_UPDRS_3_6b',],\n",
    "        'rightHandBradyOFF': ['OFF_UPDRS_3_3b','OFF_UPDRS_3_4a','OFF_UPDRS_3_5a','OFF_UPDRS_3_6a',],\n",
    "        'bradyBodyOFF': ['OFF_UPDRS_3_14'],\n",
    "        'leftHandTremorOFF': ['OFF_UPDRS_3_15b','OFF_UPDRS_3_16b','OFF_UPDRS_3_17b',],\n",
    "        'rightHandTremorOFF': ['OFF_UPDRS_3_15a','OFF_UPDRS_3_16a','OFF_UPDRS_3_17a',],\n",
    "        'gaitLegsOFF': ['OFF_UPDRS_3_7a','OFF_UPDRS_3_7b','OFF_UPDRS_3_8a','OFF_UPDRS_3_8b','OFF_UPDRS_3_9','OFF_UPDRS_3_10','OFF_UPDRS_3_11'],\n",
    "        'postureOFF': ['OFF_UPDRS_3_12','OFF_UPDRS_3_13'],\n",
    "        'facialOFF': ['OFF_UPDRS_3_1','OFF_UPDRS_3_2',]}\n",
    "    \n",
    "    ON_sublists = {\n",
    "        'leftHandBradyON': ['ON_UPDRS_3_3c', 'ON_UPDRS_3_4b','ON_UPDRS_3_5b','ON_UPDRS_3_6b',],\n",
    "        'rightHandBradyON': ['ON_UPDRS_3_3b','ON_UPDRS_3_4a','ON_UPDRS_3_5a','ON_UPDRS_3_6a',],\n",
    "        'bradyBodyON': ['ON_UPDRS_3_14'],\n",
    "        'leftHandTremorON': ['ON_UPDRS_3_15b','ON_UPDRS_3_16b','ON_UPDRS_3_17b',],\n",
    "        'rightHandTremorON': ['ON_UPDRS_3_15a','ON_UPDRS_3_16a','ON_UPDRS_3_17a',],\n",
    "        'gaitLegsON': ['ON_UPDRS_3_7a','ON_UPDRS_3_7b','ON_UPDRS_3_8a','ON_UPDRS_3_8b','ON_UPDRS_3_9','ON_UPDRS_3_10','ON_UPDRS_3_11'],\n",
    "        'postureON': ['ON_UPDRS_3_12','ON_UPDRS_3_13'],\n",
    "        'facialON': ['ON_UPDRS_3_1','ON_UPDRS_3_2',]}\n",
    "\n",
    "    for l,(off, on) in enumerate(zip(OFF_sublists.keys(),ON_sublists.keys())):\n",
    "        diffList = []\n",
    "        relative_diffList = []\n",
    "        for pt in updrsScores.index:\n",
    "            diffScore = np.sum(updrsScores.loc[pt][OFF_sublists[off]]) - np.sum(updrsScores.loc[pt][ON_sublists[on]]) # calculate difference between off and on scores\n",
    "            diffList.append(diffScore)\n",
    "            relative_diffList.append(diffScore / (len(OFF_sublists[off])*4)) # take sum-difference, normalize to list-length by dividing by potential total score\n",
    "        # create name for column\n",
    "        rel_colName = off[:-3] + '_%diff'\n",
    "        colName = off[:-3] + '_diff'\n",
    "        updrsScores.insert(loc=l*2, value=diffList, column=colName)\n",
    "        updrsScores.insert(loc=l*2 +1, value=relative_diffList, column=rel_colName)\n",
    "    \n",
    "    # create list with patients who have more bradykinesia fluctuation on left side\n",
    "    affectedSides = []\n",
    "    for pt in updrsScores.index:\n",
    "        if updrsScores.loc[pt]['leftHandBrady_diff'] > updrsScores.loc[pt]['rightHandBrady_diff']:\n",
    "            affectedSides.append('Left')\n",
    "        else:\n",
    "            affectedSides.append('Right')\n",
    "    updrsScores.insert(loc=0, value=affectedSides, column='Side')\n",
    "    \n",
    "    return updrsScores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loadFeatures(minBradyDiff, updrsScores, windowLen=60):\n",
    "    '''\n",
    "    Input:\n",
    "    - minBradyDiff = minimal difference in brady-updrs-subscores to include in analysis\n",
    "    - updrsScores: dataframe with clinical scores of all patients\n",
    "    \n",
    "    Returns:\n",
    "    - accData: one dictionary with accData, each patient has a dictionary per side, \n",
    "    containing pre and post session.\n",
    "    \n",
    "    '''\n",
    "    \n",
    "    # select patients to involve in analysis\n",
    "    selectedIDs = []\n",
    "    for pt in updrsScores.index:\n",
    "        if updrsScores.loc[pt]['leftHandBrady_diff'] > minBradyDiff:\n",
    "            selectedIDs.append(pt)\n",
    "        elif updrsScores.loc[pt]['rightHandBrady_diff'] > minBradyDiff:\n",
    "            selectedIDs.append(pt)\n",
    "        else:\n",
    "            continue\n",
    "    # removing 022 and 080 if included based on clinical difference, because data is not large enough for holdout\n",
    "    for pt in ['022','080']:\n",
    "        if pt in selectedIDs:\n",
    "            selectedIDs.remove(pt)\n",
    "    # dictionary for all acc data\n",
    "    features = {}\n",
    "    for pt in selectedIDs:\n",
    "        features[pt] = pd.read_csv(os.path.join(path,'data','features_oct20/%s_%isec_features.csv' % (pt,windowLen)))\n",
    "    \n",
    "    return features\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "updrsScores = loadUPDRSscores()\n",
    "features = loadFeatures(minBradyDiff=0.5, updrsScores=updrsScores, windowLen=60)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Individual Model Classification (incl. activity-filter sub analysis)\n",
    "\n",
    "Splits data in to x and y dataframes per patient.\n",
    "Uses always same hyperparameters, no hyperparameter optimization.\n",
    "\n",
    "Includes:\n",
    "- Balanced data regarding pre and post medication in training and test data\n",
    "- Stratified, continuous selection of hold-out data: Takes timeblocks of 10% in pre-medication timespan, and 10% in post-medication timespan as holdout-validation data.\n",
    "- Alternative data splitting method: <em>'Distributed data splitting'</em>, uses 2% of every 10-% data block as hold-out, the other 8% as trainin data. (0-8% is training, 8-10% is test, 10-18% is training, 18-20 is test, etc.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Function to perform true predictions for individual model (02DEC20)\n",
    "\n",
    "def indiv_classification_true(clsf, features, feat_select, act_filter, save):\n",
    "    '''\n",
    "    Function calculates classification-outcomes for individual model approach with SVM classifier.\n",
    "    Input: \n",
    "    - clsf: type of classification function: aupport vector: 'SV' or random forest: 'RF'\n",
    "    - features-dataframe containing all patients, all 60-sec features, and medication-states.\n",
    "    - feat_select: 'all' to include all features as predictors, use 'svm' to only include signal\n",
    "    vector magnitude features, '4' only includes the selected 4 main features.\n",
    "    - act_filter: defines activity filter or not\n",
    "    - save: has to be binary, determining whether outcomes are written to results-folder.\n",
    "    \n",
    "    Output: table containing all 41 outcomes from different splits, per patient.\n",
    "    \n",
    "    Data is balanced 50/50 between pre- and post-medication moments.\n",
    "    At least 40 minutes of pre and post are present in features (selected in feature extraction).\n",
    "    Number of folds for cross-validation is individually adjusted to create folds of +/- 10 minutes.\n",
    "    \n",
    "    Data-splitting is adjusted: 2 minutes before and after test-data are not included in training data.\n",
    "    To increase statistical power, 41 different datasplits are performed per patient. Data-split #0 takes % 0-10,\n",
    "    and % 50-60 of features as test-data, #1 takes % 1-11, and % 51-61, ..., #40 takes % 40-50, and % 90-100 as\n",
    "    test-data. The remaining data is always trainings-data (except for the neighbouring 2 minutes).\n",
    "    \n",
    "    SVM classifier using hyperparameters: C 10, gamma 0.001, kernel rbf.\n",
    "\n",
    "    Fits once the model on the real X_train[pt] and y_train[pt]. This is fitting on individualized data.\n",
    "    Calculates for this one real trial different performance metrics.\n",
    "\n",
    "    '''\n",
    "    starttime = timer()\n",
    "\n",
    "## DATA PREP for individual-holdout validation\n",
    "    pts = list(features.keys())\n",
    "    keys = list(features[pts[0]].keys())\n",
    "    if feat_select == 'all':\n",
    "        sel = [k != 'medState' for k in keys]\n",
    "        preds = list(compress(keys,sel))  \n",
    "    elif feat_select == 'SVM':\n",
    "        sel = [k[:3]=='SVM' for k in keys]\n",
    "        preds = list(compress(keys,sel))\n",
    "    elif feat_select == '4':\n",
    "        sel = [k in ['SVM_maxAcc','SVM_coefVar','SVM_RMS','SVM_specPow_totalu4Hz'] for k in keys]\n",
    "        preds = list(compress(keys,sel))\n",
    "    \n",
    "    metricNames = ['Accuracy', 'AUROC','F1_score', 'Precision', 'Recall', 'FPR'] # define metrics to check significancy for\n",
    "    allpreds = pd.DataFrame(index=np.arange(41), columns=[pt+'_'+m for pt in pts for m in metricNames]) # 2d-df to save every pred per pt, per split\n",
    "    ptDataSizes = pd.DataFrame(index=pts, columns=['n_ft_rows'])\n",
    "    for pt in pts:\n",
    "        trueLists = {} # list to store values during calculation\n",
    "        for metric in metricNames:\n",
    "            trueLists[metric]=[]\n",
    "        \n",
    "        if act_filter:\n",
    "            presize = np.sum(features[pt][features[pt]['SVM_variance'] > 0.3]['medState']==0)\n",
    "            postsize = np.sum(features[pt][features[pt]['SVM_variance'] > 0.3]['medState']==1)\n",
    "            minSize=np.min([presize,postsize]) # minimal number of rows per med-state, after actvity filter\n",
    "            # selecting only rows with SVM_variance above threshold\n",
    "            actData = features[pt][features[pt]['SVM_variance'] > 0.3].reset_index(drop=True)\n",
    "            # create predData per pt, per side, with balanced rows pre and post med (minSize)\n",
    "            predData = pd.concat([actData[actData['medState']==0][:minSize],actData[actData['medState']==1][:minSize]]).reset_index(drop=True)\n",
    "        \n",
    "        else:\n",
    "            # minSize is shortest length of samples(rows) per pt-side of dataframe, on or off samples\n",
    "            minSize = min(sum(features[pt]['medState']==0),sum(features[pt]['medState']==1))\n",
    "            # create predData per pt, per side, with balanced rows pre and post med (minSize)\n",
    "            predData = pd.concat([features[pt][features[pt]['medState']==0][:minSize],features[pt][features[pt]['medState']==1][:minSize]]).reset_index(drop=True)\n",
    "        ptDataSizes.loc[pt]['n_ft_rows'] = predData.shape[0]    \n",
    "        # standarize data \n",
    "        scaler = StandardScaler()\n",
    "        scaler.fit(predData[predData['medState']==0][preds]) # standardized only on pre-medication data\n",
    "\n",
    "        # create x and y dataset for this patient\n",
    "        x = pd.DataFrame(data = scaler.transform(predData[preds]), columns = preds)\n",
    "        y = predData['medState']\n",
    "        \n",
    "        ### real prediction with accurate medication-labels (y) ###\n",
    "        '''\n",
    "        leading to 41 different outcome values for true-labels\n",
    "        data splitting; creating 41 separate distributions of training/test splits with more temporal independence \n",
    "        first create individual dictionaries with row-numbers corresponding with percentage of data to split data\n",
    "        \n",
    "        there will be 41 different data-splits using the row-numbers corresponding to the percentile of data\n",
    "        the data-splits will start at row 0, 1, 2, ..., until 40\n",
    "        for every data-split the real-prediction will be performed, and a n-permutation-test will be performed to determine the significancy\n",
    "        ''' \n",
    "        # dictionary of rownumbers corresponding with % of data\n",
    "        rows = {}\n",
    "        for r in np.arange(101):\n",
    "            rows[r] = int(r/100*x.shape[0])\n",
    "        \n",
    "        ## First, Every iteration starts with data splitting, thereafter prediction\n",
    "        for n_split in np.arange(41): # test with 20-fold cv, every second percentage new data split\n",
    "            n_split = n_split # *2for test with 20-fold cv\n",
    "            X_test = pd.concat([x[rows[n_split]:rows[n_split+10]],x[rows[n_split+50]:rows[n_split+60]]], ignore_index=True)\n",
    "            y_test = pd.concat([y[rows[n_split]:rows[n_split+10]],y[rows[n_split+50]:rows[n_split+60]]], ignore_index=True)\n",
    "\n",
    "            if n_split < 3: # for splits 0-1-2, no rows before the test-split are selected for training\n",
    "                X_train = pd.concat([x[rows[n_split+12]:rows[50]],x[rows[n_split+62]:rows[100]]], ignore_index=True)\n",
    "                y_train = pd.concat([y[rows[n_split+12]:rows[50]],y[rows[n_split+62]:rows[100]]], ignore_index=True)\n",
    "                \n",
    "            elif n_split > 37: # for splits 38-39-40, no minutes after the test-split are selected for training\n",
    "                X_train = pd.concat([x[rows[0]:rows[n_split-2]],x[rows[50]:rows[n_split+48]]], ignore_index=True)\n",
    "                y_train = pd.concat([y[rows[0]:rows[n_split-2]],y[rows[50]:rows[n_split+48]]], ignore_index=True)\n",
    "                \n",
    "            else: # minutes before and after the test-split are included in training-split\n",
    "                X_train_a = pd.concat([x[rows[n_split+12]:rows[50]],x[rows[n_split+62]:rows[100]]], ignore_index=True)\n",
    "                X_train_b = pd.concat([x[rows[0]:rows[n_split-2]],x[rows[50]:rows[n_split+48]]], ignore_index=True)\n",
    "                X_train = pd.concat([X_train_a,X_train_b], ignore_index=True)\n",
    "                y_train_a = pd.concat([y[rows[n_split+12]:rows[50]],y[rows[n_split+62]:rows[100]]], ignore_index=True)\n",
    "                y_train_b = pd.concat([y[rows[0]:rows[n_split-2]],y[rows[50]:rows[n_split+48]]], ignore_index=True)\n",
    "                y_train = pd.concat([y_train_a,y_train_b], ignore_index=True)\n",
    "\n",
    "            ## Second, prediction part, define classifier\n",
    "            if clsf == 'SV': # check which part can be used for randomforest\n",
    "                sv = SVC(C=10, kernel='rbf', gamma=0.001, probability=True) \n",
    "                sv.fit(X_train, y_train) # fitting on 80% individual data, with true labels\n",
    "                y_preds = sv.predict(X_test) # predicting\n",
    "                y_probas = pd.DataFrame(data = sv.predict_proba(X_test), columns = sv.classes_) \n",
    "            elif clsf == 'RF':\n",
    "                rf = RandomForestClassifier(max_depth=5, min_samples_leaf=1, min_samples_split=2,n_estimators=100)\n",
    "                rf.fit(X_train, y_train) # fitting on 80% individual data, with true labels\n",
    "                y_preds = rf.predict(X_test) # predicting\n",
    "                y_probas = pd.DataFrame(data = rf.predict_proba(X_test), columns = rf.classes_) \n",
    "            y_true = y_test\n",
    "            tn, fp, fn, tp = confusion_matrix(y_true,y_preds).ravel()\n",
    "            # add outcomes to dedicated lists\n",
    "            trueLists['Accuracy'].append(accuracy_score(y_true,y_preds))\n",
    "            trueLists['AUROC'].append(roc_auc_score(y_true, y_probas[1]))\n",
    "            trueLists['F1_score'].append(f1_score(y_true,y_preds))\n",
    "            trueLists['Precision'].append(precision_score(y_true, y_preds)) # precision/PPV\n",
    "            trueLists['Recall'].append(recall_score(y_true, y_preds)) # sensitivity/recall/TPR\n",
    "            trueLists['FPR'].append(fp / (fp+tn)) # false-positive, false-alarm rate\n",
    "\n",
    "        for metric in metricNames: # all data-splits for real data performed, save values in dataframe\n",
    "            allpreds[pt+'_'+metric] = trueLists[metric]\n",
    "            \n",
    "        print(pt,'all predictions performed')\n",
    "    \n",
    "    endtime = timer()\n",
    "    timePassed = (endtime-starttime)/60 # gives time in minutes\n",
    "    \n",
    "    if save == True:\n",
    "        if act_filter == False:\n",
    "            allpreds.to_csv(os.path.join(path,'results','indiv_preds_41splits_%s_%sfts_c10g001.csv' % (clsf,feat_select)), header=True, index=True)\n",
    "#             ptDataSizes.to_csv(os.path.join(path,'results','indiv_preds_dataSizes_all.csv'), header=True, index=True)\n",
    "        elif act_filter == True:\n",
    "            allpreds.to_csv(os.path.join(path,'results','indiv_preds_41splits_%s_%sfts_c10g001_ACTFILTER.csv' % (clsf,feat_select)), header=True, index=True)\n",
    "#             ptDataSizes.to_csv(os.path.join(path,'results','indiv_preds_dataSizes_ACTFILTER.csv'), header=True, index=True)    \n",
    "    print('Minutes passed: %.1f' % timePassed)\n",
    "    \n",
    "    \n",
    "    return allpreds\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "002 all predictions performed\n",
      "012 all predictions performed\n",
      "013 all predictions performed\n",
      "014 all predictions performed\n",
      "015 all predictions performed\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "017 all predictions performed\n",
      "018 all predictions performed\n",
      "023 all predictions performed\n",
      "024 all predictions performed\n",
      "038 all predictions performed\n",
      "039 all predictions performed\n",
      "043 all predictions performed\n",
      "047 all predictions performed\n",
      "051 all predictions performed\n",
      "054 all predictions performed\n",
      "058 all predictions performed\n",
      "063 all predictions performed\n",
      "065 all predictions performed\n",
      "079 all predictions performed\n",
      "090 all predictions performed\n",
      "Minutes passed: 2.1\n",
      "002 all predictions performed\n",
      "012 all predictions performed\n",
      "013 all predictions performed\n",
      "014 all predictions performed\n",
      "015 all predictions performed\n",
      "017 all predictions performed\n",
      "018 all predictions performed\n",
      "023 all predictions performed\n",
      "024 all predictions performed\n",
      "038 all predictions performed\n",
      "039 all predictions performed\n",
      "043 all predictions performed\n",
      "047 all predictions performed\n",
      "051 all predictions performed\n",
      "054 all predictions performed\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "058 all predictions performed\n",
      "063 all predictions performed\n",
      "065 all predictions performed\n",
      "079 all predictions performed\n",
      "090 all predictions performed\n",
      "Minutes passed: 2.0\n",
      "002 all predictions performed\n",
      "012 all predictions performed\n",
      "013 all predictions performed\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "014 all predictions performed\n",
      "015 all predictions performed\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "017 all predictions performed\n",
      "018 all predictions performed\n",
      "023 all predictions performed\n",
      "024 all predictions performed\n",
      "038 all predictions performed\n",
      "039 all predictions performed\n",
      "043 all predictions performed\n",
      "047 all predictions performed\n",
      "051 all predictions performed\n",
      "054 all predictions performed\n",
      "058 all predictions performed\n",
      "063 all predictions performed\n",
      "065 all predictions performed\n",
      "079 all predictions performed\n",
      "090 all predictions performed\n",
      "Minutes passed: 2.0\n",
      "002 all predictions performed\n",
      "012 all predictions performed\n",
      "013 all predictions performed\n",
      "014 all predictions performed\n",
      "015 all predictions performed\n",
      "017 all predictions performed\n",
      "018 all predictions performed\n",
      "023 all predictions performed\n",
      "024 all predictions performed\n",
      "038 all predictions performed\n",
      "039 all predictions performed\n",
      "043 all predictions performed\n",
      "047 all predictions performed\n",
      "051 all predictions performed\n",
      "054 all predictions performed\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "058 all predictions performed\n",
      "063 all predictions performed\n",
      "065 all predictions performed\n",
      "079 all predictions performed\n",
      "090 all predictions performed\n",
      "Minutes passed: 1.9\n"
     ]
    }
   ],
   "source": [
    "### EXECUTE REAL PREDICTIONS FOR INDIVIDUAL MODEL\n",
    "# repeat for all 4 combinations (features and activity filter)\n",
    "for f_sel in ['4','all']:\n",
    "    for act in [False, True]:\n",
    "        indiv_classification_true(clsf='RF', features=features, feat_select=f_sel, \n",
    "                                  act_filter=act, save=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Function to perform permutations for individual model\n",
    "\n",
    "def indiv_classification_perm(clsf, features, act_filter, n_perm, save):\n",
    "    '''\n",
    "    Function calculates classification-outcomes for individual model approach with SVM classifier.\n",
    "    Input: \n",
    "    Input: \n",
    "    - clsf: type of classification function: aupport vector: 'SV' or random forest: 'RF'\n",
    "    - features: dataframe containing all patients, all 60-sec features, and medication-states.\n",
    "    - feat_select: 'all' to include all features as predictors, use 'svm' to only include signal\n",
    "    vector magnitude features, '4' only includes the selected 4 main features.\n",
    "    - act_filter: defines activity filter or not\n",
    "    - n_perm defines number of permutations to perform.\n",
    "    - save: has to be binary, determining whether outcomes are written to results-folder.\n",
    "    \n",
    "    Output: table with mean outcome of individual model per patient (mean is over the 41 possible data-splits),\n",
    "    including p-values of permutation-test. \n",
    "    \n",
    "    Data is balanced 50/50 between pre- and post-medication moments.\n",
    "    At least 40 minutes of pre and post are present in features (selected in feature extraction).\n",
    "    Number of folds for cross-validation is individually adjusted to create folds of +/- 10 minutes.\n",
    "    \n",
    "    Data-splitting is adjusted: 2 minutes before and after test-data are not included in training data.\n",
    "    To increase statistical power, 41 different datasplits are performed per patient. Data-split #0 takes % 0-10,\n",
    "    and % 50-60 of features as test-data, #1 takes % 1-11, and % 51-61, ..., #40 takes % 40-50, and % 90-100 as\n",
    "    test-data. The remaining data is always trainings-data (except for the neighbouring 2 minutes).\n",
    "    \n",
    "    SVM classifier using hyperparameters: C 10, gamma 0.001, kernel rbf.\n",
    "\n",
    "    Fits once the model on the real X_train[pt] and y_train[pt]. This is fitting on individualized data.\n",
    "    Calculates this n_permutations times, with randomized y_train labels. This is model\n",
    "     training with random input. For every random trial the outcome scores are calculated \n",
    "     and stored in the 'bins'. These bins represent all random performance values.\n",
    "     P-value can be calculated by: \n",
    "     (# of permutations performing better then actual score + 1) / (# of permutations + 1)\n",
    "     Based on sci-kit permutation_test_score, source: http://www.jmlr.org/papers/volume11/ojala10a/ojala10a.pdf\n",
    "\n",
    "     Overview of predictive metrics: https://en.wikipedia.org/wiki/F1_score\n",
    "\n",
    "    '''\n",
    "    starttime = timer()\n",
    "\n",
    "    pts = list(features.keys())\n",
    "    keys = list(features[pts[0]].keys())\n",
    "    ft_sel = [k != 'medState' for k in keys]\n",
    "    preds = list(compress(keys,ft_sel))      \n",
    "    metricNames = ['Accuracy', 'AUROC','F1_score', 'Precision', 'Recall', 'FPR'] # define metrics to check significancy for\n",
    "    allperms = pd.DataFrame(index=np.arange(n_perm), columns=[pt+'_'+m for pt in pts for m in metricNames]) # 2d-df to save every pred per pt, per split\n",
    "    for pt in pts:        \n",
    "        if act_filter:\n",
    "            presize = np.sum(features[pt][features[pt]['SVM_variance'] > 0.3]['medState']==0)\n",
    "            postsize = np.sum(features[pt][features[pt]['SVM_variance'] > 0.3]['medState']==1)\n",
    "            minSize=np.min([presize,postsize]) # minimal number of rows per med-state, after actvity filter\n",
    "            # selecting only rows with SVM_variance above threshold\n",
    "            actData = features[pt][features[pt]['SVM_variance'] > 0.3].reset_index(drop=True)\n",
    "            # create predData per pt, per side, with balanced rows pre and post med (minSize)\n",
    "            predData = pd.concat([actData[actData['medState']==0][:minSize],actData[actData['medState']==1][:minSize]]).reset_index(drop=True) \n",
    "        else:\n",
    "            # minSize is shortest length of samples(rows) per pt-side of dataframe, on or off samples\n",
    "            minSize = min(sum(features[pt]['medState']==0),sum(features[pt]['medState']==1))\n",
    "            # create predData per pt, per side, with balanced rows pre and post med (minSize)\n",
    "            predData = pd.concat([features[pt][features[pt]['medState']==0][:minSize],features[pt][features[pt]['medState']==1][:minSize]]).reset_index(drop=True)\n",
    "\n",
    "        # standarize data \n",
    "        scaler = StandardScaler()\n",
    "        scaler.fit(predData[predData['medState']==0][preds]) # standardized only on pre-medication data\n",
    "        # create x and y dataset for this patient\n",
    "        x = pd.DataFrame(data = scaler.transform(predData[preds]), columns = preds)\n",
    "        y = predData['medState']\n",
    "        '''\n",
    "        For every round of permutation, 41 different data-splits will be predicted,\n",
    "        leading to 41 different outcome values per permutation round, in total: n_datasplits*n_perm.\n",
    "        Original x and y, and rows-dictionary are unchanged and can be used.\n",
    "        ''' \n",
    "        rows = {} # dictionary of rownumbers corresponding with % of data\n",
    "        for r in np.arange(101):\n",
    "            rows[r] = int(r/100*x.shape[0])\n",
    "        for p in np.arange(n_perm): ## Looping over number of permutation to perform, first creating permuted data-set\n",
    "            perm_scores = {}\n",
    "            for metric in metricNames:\n",
    "                perm_scores[metric]=[]\n",
    "            random.seed(p) # set structured different random seed for every shuffle loop\n",
    "            y_random = list(y.values) # set true medication-labels to shuffle\n",
    "            random.shuffle(y_random) # shuffles labels, x and y_random are now data to use\n",
    "            y_rand = pd.DataFrame(data=y_random, columns=['med'])\n",
    "            y_random = y_rand['med']\n",
    "            ## Now, first create data-split, then perform prediction and add values to lists\n",
    "            for n_split in np.arange(41): # test with 20-fold cv, every second percentage new data split\n",
    "                n_split = n_split # *2 for test with 20-fold cv\n",
    "                X_test = pd.concat([x[rows[n_split]:rows[n_split+10]],x[rows[n_split+50]:rows[n_split+60]]], ignore_index=True)\n",
    "                y_test = pd.concat([y_random[rows[n_split]:rows[n_split+10]],y_random[rows[n_split+50]:rows[n_split+60]]], ignore_index=True)\n",
    "                if n_split < 3: # for splits 0-1-2, no rows before the test-split are selected for training\n",
    "                    X_train = pd.concat([x[rows[n_split+12]:rows[50]],x[rows[n_split+62]:rows[100]]], ignore_index=True)\n",
    "                    y_train = pd.concat([y_random[rows[n_split+12]:rows[50]],y_random[rows[n_split+62]:rows[100]]], ignore_index=True)\n",
    "                elif n_split > 37: # for splits 38-39-40, no minutes after the test-split are selected for training\n",
    "                    X_train = pd.concat([x[rows[0]:rows[n_split-2]],x[rows[50]:rows[n_split+48]]], ignore_index=True)\n",
    "                    y_train = pd.concat([y_random[rows[0]:rows[n_split-2]],y_random[rows[50]:rows[n_split+48]]], ignore_index=True)\n",
    "                else: # minutes before and after the test-split are included in training-split\n",
    "                    X_train_a = pd.concat([x[rows[n_split+12]:rows[50]],x[rows[n_split+62]:rows[100]]], ignore_index=True)\n",
    "                    X_train_b = pd.concat([x[rows[0]:rows[n_split-2]],x[rows[50]:rows[n_split+48]]], ignore_index=True)\n",
    "                    X_train = pd.concat([X_train_a,X_train_b], ignore_index=True)\n",
    "                    y_train_a = pd.concat([y_random[rows[n_split+12]:rows[50]],y_random[rows[n_split+62]:rows[100]]], ignore_index=True)\n",
    "                    y_train_b = pd.concat([y_random[rows[0]:rows[n_split-2]],y_random[rows[50]:rows[n_split+48]]], ignore_index=True)\n",
    "                    y_train = pd.concat([y_train_a,y_train_b], ignore_index=True)\n",
    "                ## Second, prediction part\n",
    "                if clsf == 'SV':\n",
    "                    sv = SVC(C=10, kernel='rbf', gamma=0.001, probability=True) # define classifier as indiv-optimized class\n",
    "                    sv.fit(X_train, y_train) # fitting on 80% individual data, with true labels\n",
    "                    y_preds = sv.predict(X_test) # predicting\n",
    "                    y_probas = pd.DataFrame(data = sv.predict_proba(X_test), columns = sv.classes_) \n",
    "                elif clsf == 'RF':\n",
    "                    rf = RandomForestClassifier(max_depth=5, min_samples_leaf=1, min_samples_split=2,n_estimators=100)\n",
    "                    rf.fit(X_train, y_train) # fitting on 80% individual data, with true labels\n",
    "                    y_preds = rf.predict(X_test) # predicting\n",
    "                    y_probas = pd.DataFrame(data = rf.predict_proba(X_test), columns = rf.classes_) \n",
    "                y_true = y_test\n",
    "                tn, fp, fn, tp = confusion_matrix(y_true,y_preds).ravel()\n",
    "                # Third, add outcomes to dedicated lists\n",
    "                perm_scores['Accuracy'].append(accuracy_score(y_true,y_preds))\n",
    "                perm_scores['AUROC'].append(roc_auc_score(y_true, y_probas[1]))\n",
    "                perm_scores['F1_score'].append(f1_score(y_true,y_preds))\n",
    "                perm_scores['Precision'].append(precision_score(y_true, y_preds)) # precision/PPV\n",
    "                perm_scores['Recall'].append(recall_score(y_true, y_preds)) # sensitivity/recall/TPR\n",
    "                perm_scores['FPR'].append(fp / (fp+tn)) # false-positive, false-alarm rate\n",
    "                # end of one data-split, in one permuted data-split\n",
    "            for metric in metricNames: # at end of 41-splits in one permutations, save mean scores in outcome df\n",
    "                allperms.loc[p][pt+'_'+metric] = np.mean(perm_scores[metric])\n",
    "            \n",
    "        print(pt,'all predictions ready')\n",
    "        ptcols = list(compress(allperms.keys(),[pt in colname for colname in list(allperms.keys())]))\n",
    "        ptperms = allperms[ptcols]\n",
    "        if save == True:\n",
    "            if act_filter == False:\n",
    "                ptperms.to_csv(os.path.join(path,'results','perms','temp',\n",
    "                                             '%s_indiv_%s_allfts_%iperms.csv' % (pt,clsf,n_perm)), header=True, index=True)\n",
    "            elif act_filter == True:\n",
    "                ptperms.to_csv(os.path.join(path,'results','perms','temp',\n",
    "                                             '%s_indiv_%s_allfts_actfilter_%iperms.csv' % (pt,cls,n_perm)), header=True, index=True)\n",
    "    \n",
    "    endtime = timer()\n",
    "    timePassed = (endtime-starttime)/60 # gives time in minutes\n",
    "    \n",
    "    if save == True:\n",
    "        if act_filter == False:\n",
    "            allperms.to_csv(os.path.join(path,'results','indiv_%iperms_41splits_%s.csv' % (n_perm,clsf)), header=True, index=True)\n",
    "        elif act_filter == True:\n",
    "            allperms.to_csv(os.path.join(path,'results','indiv_%iperms_41splits_%s_ACTFILTER.csv' % (n_perm,clsf)), header=True, index=True)\n",
    "        \n",
    "    print('Minutes passed: %.1f' % timePassed)    \n",
    "    \n",
    "    return allperms\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "002 all predictions ready\n",
      "012 all predictions ready\n",
      "013 all predictions ready\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "014 all predictions ready\n",
      "015 all predictions ready\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "017 all predictions ready\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "018 all predictions ready\n",
      "023 all predictions ready\n",
      "024 all predictions ready\n",
      "038 all predictions ready\n",
      "039 all predictions ready\n",
      "043 all predictions ready\n",
      "047 all predictions ready\n",
      "051 all predictions ready\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "054 all predictions ready\n",
      "058 all predictions ready\n",
      "063 all predictions ready\n",
      "065 all predictions ready\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "079 all predictions ready\n",
      "090 all predictions ready\n",
      "Minutes passed: 10.0\n",
      "002 all predictions ready\n",
      "012 all predictions ready\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "013 all predictions ready\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "014 all predictions ready\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "015 all predictions ready\n",
      "017 all predictions ready\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "018 all predictions ready\n",
      "023 all predictions ready\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "024 all predictions ready\n",
      "038 all predictions ready\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "039 all predictions ready\n",
      "043 all predictions ready\n",
      "047 all predictions ready\n",
      "051 all predictions ready\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "054 all predictions ready\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "058 all predictions ready\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "063 all predictions ready\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "065 all predictions ready\n",
      "079 all predictions ready\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n",
      "/Users/jeroenhabets/Research/pd_sensors/env/lib/python3.7/site-packages/sklearn/metrics/_classification.py:1221: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.\n",
      "  _warn_prf(average, modifier, msg_start, len(result))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "090 all predictions ready\n",
      "Minutes passed: 9.6\n"
     ]
    }
   ],
   "source": [
    "permall = indiv_classification_perm('RF', features, False, 5, True)\n",
    "permact = indiv_classification_perm('RF', features, True, 5, True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 3: Group Model Classification (incl. activity-filter sub analysis)\n",
    "Trains model for one patient on data from nineteen remaining patients, tests on data from the one specific patient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "def group_classification_true(clsf, features, feat_select, act_filter, save):\n",
    "    \n",
    "    starttime = timer()\n",
    "    ## Data preparation for group classification\n",
    "    pts = list(features.keys())\n",
    "    keys = list(features[pts[0]].keys())\n",
    "    if feat_select == 'all':\n",
    "        sel = [k != 'medState' for k in keys]\n",
    "        preds = list(compress(keys,sel))  \n",
    "    elif feat_select == 'SVM':\n",
    "        sel = [k[:3]=='SVM' for k in keys]\n",
    "        preds = list(compress(keys,sel))\n",
    "    elif feat_select == '4':\n",
    "        sel = [k in ['SVM_maxAcc','SVM_coefVar','SVM_RMS','SVM_specPow_totalu4Hz'] for k in keys]\n",
    "        preds = list(compress(keys,sel))\n",
    "    \n",
    "    metricNames = ['Accuracy', 'AUROC','F1_score', 'Precision', 'Recall', 'FPR'] # define metrics to check significancy for\n",
    "    group_preds = pd.DataFrame(index=['pred'],columns=[pt+'_'+m for pt in pts for m in metricNames]) # 2d-df to save every pred per pt, per split\n",
    "    predDict = {} # store selected features in dictionary for test/training selection based on patients\n",
    "    for pt in pts:\n",
    "        if act_filter:\n",
    "            presize = np.sum(features[pt][features[pt]['SVM_variance'] > 0.3]['medState']==0)\n",
    "            postsize = np.sum(features[pt][features[pt]['SVM_variance'] > 0.3]['medState']==1)\n",
    "            minSize=np.min([presize,postsize]) # minimal number of rows per med-state, after actvity filter\n",
    "            # selecting only rows with SVM_variance above threshold\n",
    "            actData = features[pt][features[pt]['SVM_variance'] > 0.3].reset_index(drop=True)\n",
    "            # create predData per pt, per side, with balanced rows pre and post med (minSize)\n",
    "            predDict[pt] = pd.concat([actData[actData['medState']==0][:minSize],actData[actData['medState']==1][:minSize]]).reset_index(drop=True)\n",
    "        \n",
    "        else:\n",
    "            # minSize is shortest length of samples(rows) per pt-side of dataframe, on or off samples\n",
    "            minSize = min(sum(features[pt]['medState']==0),sum(features[pt]['medState']==1))\n",
    "            # create predData per pt, per side, with balanced rows pre and post med (minSize)\n",
    "            predDict[pt] = pd.concat([features[pt][features[pt]['medState']==0][:minSize],features[pt][features[pt]['medState']==1][:minSize]]).reset_index(drop=True)\n",
    "    # data prepared: one dict with balanced features per patient (activity filtered or not)\n",
    "\n",
    "    for pt in pts: ## Perform classification per patient that is tested\n",
    "        testData = predDict[pt] # pt in loop is test-patient\n",
    "        trainPts = pts.copy() # copy list to leave original pts list unchanged\n",
    "        trainPts.remove(pt) # list with 19 training-patients\n",
    "        trainData = pd.DataFrame(columns=predDict[pt].keys()) #create trainings df with other 19 patients\n",
    "        for trainPt in trainPts:\n",
    "            trainData = pd.concat([trainData,predDict[trainPt]]).reset_index(drop=True)\n",
    "        scaler = StandardScaler() # standarize data, fit on trainingsdata (pre med)\n",
    "        scaler.fit(trainData[trainData['medState']==0][preds]) # standardized only on pre-medication data\n",
    "        # create x and y dataset \n",
    "        x_train = pd.DataFrame(data = scaler.transform(trainData[preds]), columns = preds)\n",
    "        y_train = trainData['medState']\n",
    "        x_test = pd.DataFrame(data = scaler.transform(testData[preds]), columns = preds)\n",
    "        y_test = testData['medState']\n",
    "        ## Performing prediction\n",
    "        if clsf == 'SV':\n",
    "            sv = SVC(C=10, kernel='rbf', gamma=0.001, probability=True) # define classifier as indiv-optimized class\n",
    "            sv.fit(x_train, y_train) \n",
    "            y_preds = sv.predict(x_test) # predicting\n",
    "            y_probas = pd.DataFrame(data = sv.predict_proba(x_test), columns = sv.classes_) \n",
    "        elif clsf == 'RF':\n",
    "            rf = RandomForestClassifier(max_depth=5, min_samples_leaf=1, min_samples_split=2,n_estimators=100)\n",
    "            rf.fit(x_train, y_train) # fitting on 19 training patients\n",
    "            y_preds = rf.predict(x_test) # predicting\n",
    "            y_probas = pd.DataFrame(data = rf.predict_proba(x_test), columns = rf.classes_)\n",
    "        y_true = y_test\n",
    "        tn, fp, fn, tp = confusion_matrix(y_true,y_preds).ravel()\n",
    "        group_preds.loc['pred'][pt+'_'+'Accuracy'] = (accuracy_score(y_true,y_preds))\n",
    "        group_preds.loc['pred'][pt+'_'+'AUROC'] = (roc_auc_score(y_true, y_probas[1]))\n",
    "        group_preds.loc['pred'][pt+'_'+'F1_score'] = (f1_score(y_true,y_preds))\n",
    "        group_preds.loc['pred'][pt+'_'+'Precision'] = (precision_score(y_true, y_preds)) # precision/PPV\n",
    "        group_preds.loc['pred'][pt+'_'+'Recall'] = (recall_score(y_true, y_preds)) # sensitivity/recall/TPR\n",
    "        group_preds.loc['pred'][pt+'_'+'FPR'] = (fp / (fp+tn)) # false-positive, false-alarm rate\n",
    "#         print(pt,'is performed')\n",
    "    endtime = timer()\n",
    "    proc_time = endtime-starttime\n",
    "    print('%s, %s fts, act %s, ready, %.1f minutes passed.' % (clsf,feat_select, \n",
    "                                                               str(act_filter),proc_time/60))\n",
    "    \n",
    "    if save:\n",
    "        if act_filter == False:\n",
    "            group_preds.to_csv(os.path.join(path,'results','group_pred_%s_%sfts.csv' % (clsf, feat_select)), header=True, index=True)\n",
    "        elif act_filter == True:\n",
    "            group_preds.to_csv(os.path.join(path,'results','group_pred_%s_%sfts_ACTFILTER.csv' % (clsf,feat_select)), header=True, index=True)\n",
    "            \n",
    "    return group_preds\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RF, all fts, act True, ready, 0.2 minutes passed.\n",
      "RF, all fts, act False, ready, 0.3 minutes passed.\n",
      "RF, 4 fts, act True, ready, 0.1 minutes passed.\n",
      "RF, 4 fts, act False, ready, 0.1 minutes passed.\n"
     ]
    }
   ],
   "source": [
    "for cls in ['RF']:\n",
    "    for ft_sel in ['all','4']:\n",
    "        for act in [True,False]:\n",
    "            group_classification_true(cls, features, ft_sel, act, True)\n",
    "# SV, all fts, act True, ready, 0.6 minutes passed.\n",
    "# SV, all fts, act False, ready, 1.1 minutes passed.\n",
    "# SV, 4 fts, act True, ready, 0.1 minutes passed.\n",
    "# SV, 4 fts, act False, ready, 0.2 minutes passed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "def group_classification_perm(clsf, features, act_filter, n_perm, save):\n",
    "    \n",
    "    starttime = timer()\n",
    "    ## Data preparation for group classification\n",
    "    pts = list(features.keys())\n",
    "    keys = list(features[pts[0]].keys())\n",
    "    sel = [k != 'medState' for k in keys]\n",
    "    preds = list(compress(keys,sel))  \n",
    "    \n",
    "    metricNames = ['Accuracy', 'AUROC','F1_score', 'Precision', 'Recall', 'FPR'] # define metrics to check significancy for\n",
    "    group_perms = pd.DataFrame(index=np.arange(n_perm),columns=[pt+'_'+m for pt in pts for m in metricNames]) # 2d-df to save every pred per pt, per split\n",
    "    predDict = {} # store selected features in dictionary for test/training selection based on patients\n",
    "    for pt in pts:\n",
    "        if act_filter:\n",
    "            presize = np.sum(features[pt][features[pt]['SVM_variance'] > 0.3]['medState']==0)\n",
    "            postsize = np.sum(features[pt][features[pt]['SVM_variance'] > 0.3]['medState']==1)\n",
    "            minSize=np.min([presize,postsize]) # minimal number of rows per med-state, after actvity filter\n",
    "            # selecting only rows with SVM_variance above threshold\n",
    "            actData = features[pt][features[pt]['SVM_variance'] > 0.3].reset_index(drop=True)\n",
    "            # create predData per pt, per side, with balanced rows pre and post med (minSize)\n",
    "            predDict[pt] = pd.concat([actData[actData['medState']==0][:minSize],actData[actData['medState']==1][:minSize]]).reset_index(drop=True)\n",
    "        \n",
    "        else:\n",
    "            # minSize is shortest length of samples(rows) per pt-side of dataframe, on or off samples\n",
    "            minSize = min(sum(features[pt]['medState']==0),sum(features[pt]['medState']==1))\n",
    "            # create predData per pt, per side, with balanced rows pre and post med (minSize)\n",
    "            predDict[pt] = pd.concat([features[pt][features[pt]['medState']==0][:minSize],features[pt][features[pt]['medState']==1][:minSize]]).reset_index(drop=True)\n",
    "    # data prepared: one dict with balanced features per patient (activity filtered or not)\n",
    "\n",
    "    for pt in pts: ## Perform classification per patient that is tested\n",
    "        testData = predDict[pt] # pt in loop is test-patient\n",
    "        trainPts = pts.copy() # copy list to leave original pts list unchanged\n",
    "        trainPts.remove(pt) # list with 19 training-patients\n",
    "        trainData = pd.DataFrame(columns=predDict[pt].keys()) #create trainings df with other 19 patients\n",
    "        for trainPt in trainPts:\n",
    "            trainData = pd.concat([trainData,predDict[trainPt]]).reset_index(drop=True)\n",
    "        scaler = StandardScaler() # standarize data, fit on trainingsdata (pre med)\n",
    "        scaler.fit(trainData[trainData['medState']==0][preds]) # standardized only on pre-medication data\n",
    "        x_train = pd.DataFrame(data = scaler.transform(trainData[preds]), columns = preds)\n",
    "        y_train = trainData['medState']\n",
    "        x_test = pd.DataFrame(data = scaler.transform(testData[preds]), columns = preds)\n",
    "        y_test = testData['medState']\n",
    "        ## data prepared for randomisation and permuted classification\n",
    "        permlist={} # dict to store lists for perm-scores per metric (for every test-pt again)\n",
    "        for m in metricNames:\n",
    "            permlist[m] = []\n",
    "        ## Start actual permutations, shuffling of y-labels\n",
    "        for p in np.arange(n_perm): # loops over number of permutations\n",
    "            y_random_train = list(y_train.values) # set true labels again as starting list, y_train itself stays unchanged\n",
    "            y_random_test = list(y_test.values)\n",
    "            random.seed(p) # set different random seed for every shuffle loop\n",
    "            random.shuffle(y_random_train) # shuffles labels in y_random\n",
    "            random.shuffle(y_random_test)\n",
    "            ## Performing prediction\n",
    "            if clsf == 'SV':\n",
    "                sv = SVC(C=10, kernel='rbf', gamma=0.001, probability=True) # define classifier as indiv-optimized class\n",
    "                sv.fit(x_train, y_random_train) \n",
    "                y_preds = sv.predict(x_test) # predicting\n",
    "                y_probas = pd.DataFrame(data = sv.predict_proba(x_test), columns = sv.classes_) \n",
    "            elif clsf == 'RF':\n",
    "                rf = RandomForestClassifier(max_depth=5, min_samples_leaf=1, min_samples_split=2,n_estimators=100)\n",
    "                rf.fit(x_train, y_random_train) # fitting on 19 training patients\n",
    "                y_preds = rf.predict(x_test) # predicting\n",
    "                y_probas = pd.DataFrame(data = rf.predict_proba(x_test), columns = rf.classes_)\n",
    "            y_true = y_random_test\n",
    "            tn, fp, fn, tp = confusion_matrix(y_true,y_preds).ravel()\n",
    "            permlist['Accuracy'].append(accuracy_score(y_true,y_preds))\n",
    "            permlist['AUROC'].append(roc_auc_score(y_true, y_probas[1]))\n",
    "            permlist['F1_score'].append(f1_score(y_true,y_preds))\n",
    "            permlist['Precision'].append(precision_score(y_true, y_preds)) # precision/PPV\n",
    "            permlist['Recall'].append(recall_score(y_true, y_preds)) # sensitivity/recall/TPR\n",
    "            permlist['FPR'].append(fp / (fp+tn)) # false-positive, false-alarm rate\n",
    "        \n",
    "        for metric in metricNames: # at end of 41-splits in one permutations, save mean scores in outcome df\n",
    "            group_perms[pt+'_'+metric] = permlist[metric]    \n",
    "        print(pt,'all predictions ready')\n",
    "    \n",
    "    endtime = timer()\n",
    "    timePassed = (endtime-starttime)/60 # gives time in minutes\n",
    "    if save == True:\n",
    "        if act_filter == False:\n",
    "            group_perms.to_csv(os.path.join(path,'results','group_%s_%iperms.csv' % (clsf,n_perm)), header=True, index=True)\n",
    "        elif act_filter == True:\n",
    "            group_perms.to_csv(os.path.join(path,'results','group_%s_%iperms_ACTFILTER.csv' % (clsf,n_perm)), header=True, index=True)\n",
    "        \n",
    "    print('Minutes passed: %.1f' % timePassed)   \n",
    "\n",
    "    return group_perms\n",
    "\n",
    "## PM activity filter litearture\n",
    "#actThresholdCheck(accData, features, 'SVM_variance', .3, save=True)\n",
    "# eye-balling of episode without movement: variance +/- 0.005\n",
    "# Van Hilten '91: > 0.1 g\n",
    "# Keijsers '06: > 0.05 m/s\n",
    "# Mahadevan '20: Coef of Var > 0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "002 all predictions ready\n",
      "012 all predictions ready\n",
      "013 all predictions ready\n",
      "014 all predictions ready\n",
      "015 all predictions ready\n",
      "017 all predictions ready\n",
      "018 all predictions ready\n",
      "023 all predictions ready\n",
      "024 all predictions ready\n",
      "038 all predictions ready\n",
      "039 all predictions ready\n",
      "043 all predictions ready\n",
      "047 all predictions ready\n",
      "051 all predictions ready\n",
      "054 all predictions ready\n",
      "058 all predictions ready\n",
      "063 all predictions ready\n",
      "065 all predictions ready\n",
      "079 all predictions ready\n",
      "090 all predictions ready\n",
      "Minutes passed: 5.5\n",
      "002 all predictions ready\n",
      "012 all predictions ready\n",
      "013 all predictions ready\n",
      "014 all predictions ready\n",
      "015 all predictions ready\n",
      "017 all predictions ready\n",
      "018 all predictions ready\n",
      "023 all predictions ready\n",
      "024 all predictions ready\n",
      "038 all predictions ready\n",
      "039 all predictions ready\n",
      "043 all predictions ready\n",
      "047 all predictions ready\n",
      "051 all predictions ready\n",
      "054 all predictions ready\n",
      "058 all predictions ready\n",
      "063 all predictions ready\n",
      "065 all predictions ready\n",
      "079 all predictions ready\n",
      "090 all predictions ready\n",
      "Minutes passed: 3.1\n",
      "002 all predictions ready\n",
      "012 all predictions ready\n",
      "013 all predictions ready\n",
      "014 all predictions ready\n",
      "015 all predictions ready\n",
      "017 all predictions ready\n",
      "018 all predictions ready\n",
      "023 all predictions ready\n",
      "024 all predictions ready\n",
      "038 all predictions ready\n",
      "039 all predictions ready\n",
      "043 all predictions ready\n",
      "047 all predictions ready\n",
      "051 all predictions ready\n",
      "054 all predictions ready\n",
      "058 all predictions ready\n",
      "063 all predictions ready\n",
      "065 all predictions ready\n",
      "079 all predictions ready\n",
      "090 all predictions ready\n",
      "Minutes passed: 1.2\n",
      "002 all predictions ready\n",
      "012 all predictions ready\n",
      "013 all predictions ready\n",
      "014 all predictions ready\n",
      "015 all predictions ready\n",
      "017 all predictions ready\n",
      "018 all predictions ready\n",
      "023 all predictions ready\n",
      "024 all predictions ready\n",
      "038 all predictions ready\n",
      "039 all predictions ready\n",
      "043 all predictions ready\n",
      "047 all predictions ready\n",
      "051 all predictions ready\n",
      "054 all predictions ready\n",
      "058 all predictions ready\n",
      "063 all predictions ready\n",
      "065 all predictions ready\n",
      "079 all predictions ready\n",
      "090 all predictions ready\n",
      "Minutes passed: 0.9\n"
     ]
    }
   ],
   "source": [
    "for cls in ['SV','RF']:\n",
    "    for act in [False,True]:\n",
    "        group_classification_perm(cls, features, act, 5, True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 4: Sub analyses on influence on amount of training data for group model \n",
    "\n",
    "\n",
    "Only Group Model analyzed, individual data is not large enough to decrease amount of training data used per paraticipant.\n",
    "\n",
    "Test data: \n",
    "\n",
    "always data from one same patient\n",
    "\n",
    "Training data:\n",
    "\n",
    "1: start with data from 1 other patient (randomly selected); 5 random repititions\n",
    "\n",
    "2: data from 2 other patients; 5 random rep's\n",
    "\n",
    "3: data from 3 other patients; 5 random rep's\n",
    "\n",
    "..\n",
    "\n",
    "19: data from all 19 other patients; only one repitition possible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def amount_training_data(features, clsf, act_filter, test_pts, save):\n",
    "    '''\n",
    "    Function calculates classification-outcomes for individual model approach with SVM classifier.\n",
    "    Input: \n",
    "    - features-dataframe containing all patients, all 60-sec features, and medication-states.\n",
    "    - ptList: contains pt names to analyse.\n",
    "    - save: has to be binary, determining whether outcomes are written to results-folder.\n",
    "    \n",
    "    Output: table with mean outcome of individual model per patient (mean is over the 41 possible data-splits),\n",
    "    including p-values of permutation-test. \n",
    "        \n",
    "    SVM classifier using hyperparameters: C 10, gamma 0.001, kernel rbf.\n",
    "    \n",
    "    Discussionpoints: standarisation, per patient, or per combi? Now: per combi.\n",
    "\n",
    "    '''\n",
    "    starttime = timer()\n",
    "    pts = list(features.keys())\n",
    "    preds = list(features[pts[0]].keys()) # define pred's\n",
    "    preds.remove('medState')\n",
    "    \n",
    "    metrics = ['Accuracy', 'AUROC', ] \n",
    "    outcomeCols = [metr + '_mean' for metr in metrics]\n",
    "    outcomeCols.extend([metr + '_sd' for metr in metrics])\n",
    "    outcomeCols = [col + '_n%i' %npt for npt in np.arange(1,20) for col in outcomeCols]\n",
    "    outcomes = pd.DataFrame(index=pts, columns=outcomeCols) # store final results\n",
    "    ## Preparing data\n",
    "    predDict={} # dict to store prediction-ready data per pt\n",
    "    for pt in pts: \n",
    "        if act_filter:\n",
    "            presize = np.sum(features[pt][features[pt]['SVM_variance'] > 0.3]['medState']==0)\n",
    "            postsize = np.sum(features[pt][features[pt]['SVM_variance'] > 0.3]['medState']==1)\n",
    "            minSize=np.min([presize,postsize]) # minimal number of rows per med-state, after actvity filter\n",
    "            # selecting only rows with SVM_variance above threshold\n",
    "            actData = features[pt][features[pt]['SVM_variance'] > 0.3].reset_index(drop=True)\n",
    "            # create predData per pt, per side, with balanced rows pre and post med (minSize)\n",
    "            predDict[pt] = pd.concat([actData[actData['medState']==0][:minSize],actData[actData['medState']==1][:minSize]]).reset_index(drop=True)\n",
    "        else:\n",
    "            minSize = min(sum(features[pt]['medState']==0),sum(features[pt]['medState']==1)) # minSize is shortest length of samples(rows) per pt-side of dataframe, on or off samples\n",
    "            predDict[pt] = pd.concat([features[pt][features[pt]['medState']==0][:minSize],features[pt][features[pt]['medState']==1][:minSize]]).reset_index(drop=True)\n",
    "    ## data splitting with increasing number of training patients\n",
    "    for test_pt in test_pts: # loop over inserted patients to test\n",
    "        for nPts in np.arange(1,20): # loop over number of training pt's [1-19]\n",
    "            n_reps = 5 # number of repetitions, for every number of training pt's if possible\n",
    "            if nPts == 19: # exceptions where 5 different repetitions are not possible\n",
    "                n_reps = 1\n",
    "            elif nPts == 18:\n",
    "                n_reps = 2\n",
    "            elif nPts == 17:\n",
    "                n_reps = 3\n",
    "            elif nPts == 16:\n",
    "                n_reps = 4\n",
    "            other_pts = list(features.keys()).copy() # list with all potential training pt's\n",
    "            other_pts.remove(test_pt)\n",
    "            acc_list = [] # list to store results (max5) per number of training pts\n",
    "            auc_list = []\n",
    "            for rep in np.arange(n_reps):\n",
    "                random.seed(rep) # make sure every repetitions has always the same random.seed\n",
    "                trainPts = random.sample(other_pts,nPts) # takes random sample of n=nPts, from other other_pts\n",
    "                # make one df of data of trainings patients\n",
    "                for n,train_pt in enumerate(trainPts):\n",
    "                    if n == 0:\n",
    "                        trainData = predDict[train_pt]\n",
    "                    else:\n",
    "                        trainData = pd.concat([trainData,predDict[train_pt]]).reset_index(drop=True)\n",
    "                # standarize data \n",
    "                scaler = StandardScaler()\n",
    "                scaler.fit(trainData[trainData['medState']==0][preds]) # standardized only on pre-medication data\n",
    "                # create x and y dataset for this repetition, of this nr of trainingpts, for this test-patient\n",
    "                x_train = pd.DataFrame(data = scaler.transform(trainData[preds]), columns = preds)\n",
    "                y_train = trainData['medState']\n",
    "                x_test = pd.DataFrame(data = scaler.transform(predDict[test_pt][preds]), columns = preds)\n",
    "                y_test = predDict[test_pt]['medState']\n",
    "                ## Performing prediction            \n",
    "                if clsf == 'SV':\n",
    "                    sv = SVC(C=10, kernel='rbf', gamma=0.001, probability=True) # define classifier as indiv-optimized class\n",
    "                    sv.fit(x_train, y_train) \n",
    "                    y_preds = sv.predict(x_test) # predicting\n",
    "                    y_probas = pd.DataFrame(data = sv.predict_proba(x_test), columns = sv.classes_) \n",
    "                elif clsf == 'RF':\n",
    "                    rf = RandomForestClassifier(max_depth=5, min_samples_leaf=1, min_samples_split=2,n_estimators=100)\n",
    "                    rf.fit(x_train, y_train) # fitting on 19 training patients\n",
    "                    y_preds = rf.predict(x_test) # predicting\n",
    "                    y_probas = pd.DataFrame(data = rf.predict_proba(x_test), columns = rf.classes_)\n",
    "                y_true = y_test\n",
    "    #             tn, fp, fn, tp = confusion_matrix(y_true,y_preds).ravel()\n",
    "                acc_list.append(accuracy_score(y_true,y_preds))\n",
    "                auc_list.append(roc_auc_score(y_true, y_probas[1]))\n",
    "    #             indiv_outcomes.iloc[n_split]['F1_score'] = f1_score(y_true,y_preds)\n",
    "    #             indiv_outcomes.iloc[n_split]['Precision'] = precision_score(y_true, y_preds) # precision/PPV\n",
    "    #             indiv_outcomes.iloc[n_split]['Recall'] = recall_score(y_true, y_preds) # sensitivity/recall/TPR\n",
    "    #             indiv_outcomes.iloc[n_split]['FPR'] = fp / (fp+tn) # false-alarm rate\n",
    "            outcomes.loc[test_pt]['Accuracy_mean_n%i' % nPts] = np.mean(acc_list) # at end of 5 rep's            \n",
    "            outcomes.loc[test_pt]['Accuracy_sd_n%i' % nPts] = np.std(acc_list)\n",
    "            outcomes.loc[test_pt]['AUROC_mean_n%i' % nPts] = np.mean(auc_list)\n",
    "            outcomes.loc[test_pt]['AUROC_sd_n%i' % nPts] = np.std(auc_list)\n",
    "        print('Test patient',test_pt,'ready')    \n",
    "    endtime = timer()\n",
    "    print('Time passed: %.1f minutes' % ((endtime-starttime)/60)) # gives time in minutes\n",
    "    \n",
    "    if save == True:\n",
    "        if act_filter:\n",
    "            outcomes.to_csv(os.path.join(path,'results/preds_%s_nr_train_pts_actfilter.csv' % clsf), header=True, index=True)\n",
    "        else:\n",
    "            outcomes.to_csv(os.path.join(path,'results/preds_%s_nr_train_pts.csv' % clsf), header=True, index=True)\n",
    "    \n",
    "    return outcomes\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test patient 002 ready\n",
      "Test patient 012 ready\n",
      "Test patient 013 ready\n",
      "Test patient 014 ready\n",
      "Test patient 015 ready\n",
      "Test patient 017 ready\n",
      "Test patient 018 ready\n",
      "Test patient 023 ready\n",
      "Test patient 024 ready\n",
      "Test patient 038 ready\n",
      "Test patient 039 ready\n",
      "Test patient 043 ready\n",
      "Test patient 047 ready\n",
      "Test patient 051 ready\n",
      "Test patient 054 ready\n",
      "Test patient 058 ready\n",
      "Test patient 063 ready\n",
      "Test patient 065 ready\n",
      "Test patient 079 ready\n",
      "Test patient 090 ready\n",
      "Time passed: 16.4 minutes\n"
     ]
    }
   ],
   "source": [
    "allPts = list(features.keys())\n",
    "test_pts = allPts\n",
    "for cls in ['SV']:#'RF',\n",
    "    amount_training_data(features, cls, True, test_pts, True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['043_indiv_SV_allfts_actfilter_5000perms.csv',\n",
       " '039_indiv_SV_allfts_actfilter_5000perms.csv',\n",
       " '002_indiv_SV_allfts_actfilter_5000perms.csv',\n",
       " '012_group_SV_allfts_actfilter_5000perms.csv',\n",
       " '013_group_SV_allfts_actfilter_5000perms.csv',\n",
       " '038_indiv_SV_allfts_actfilter_5000perms.csv',\n",
       " '079_indiv_SV_allfts_actfilter_5000perms.csv',\n",
       " '014_group_SV_allfts_actfilter_5000perms.csv',\n",
       " '047_indiv_SV_allfts_actfilter_5000perms.csv',\n",
       " '051_indiv_SV_allfts_actfilter_5000perms.csv',\n",
       " '065_indiv_SV_allfts_actfilter_5000perms.csv',\n",
       " '024_indiv_SV_allfts_actfilter_5000perms.csv',\n",
       " '013_indiv_SV_allfts_actfilter_5000perms.csv',\n",
       " '090_indiv_SV_allfts_actfilter_5000perms.csv',\n",
       " '012_indiv_SV_allfts_actfilter_5000perms.csv',\n",
       " '002_group_SV_allfts_actfilter_5000perms.csv',\n",
       " '018_indiv_SV_allfts_actfilter_5000perms.csv',\n",
       " '023_indiv_SV_allfts_actfilter_5000perms.csv',\n",
       " '017_indiv_SV_allfts_actfilter_5000perms.csv',\n",
       " '058_indiv_SV_allfts_actfilter_5000perms.csv',\n",
       " '063_indiv_SV_allfts_actfilter_5000perms.csv',\n",
       " '015_indiv_SV_allfts_actfilter_5000perms.csv',\n",
       " '054_indiv_SV_allfts_actfilter_5000perms.csv',\n",
       " '014_indiv_SV_allfts_actfilter_5000perms.csv']"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## ADDING INDIVIDUAL PERMUTATIONS TO GROUP FILE\n",
    "\n",
    "listdir(os.path.join(path,'results/perms/temp'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'015',\n",
       " '017',\n",
       " '018',\n",
       " '023',\n",
       " '024',\n",
       " '038',\n",
       " '039',\n",
       " '043',\n",
       " '047',\n",
       " '051',\n",
       " '054',\n",
       " '058',\n",
       " '063',\n",
       " '065',\n",
       " '079',\n",
       " '090'}"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "f = pd.read_csv(os.path.join(path,'results/perms/indiv_SV_4fts_actfilter_5000perms.csv'), index_col=0)\n",
    "set([k[:3] for k in f.keys()])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# main = pd.read_csv(os.path.join(path,'results/perms/002_043_indiv_SV_allfts_5000perms.csv'), index_col=0)\n",
    "# main\n",
    "\n",
    "# for pt in ['047','051','054','058','063','065','079','090']:\n",
    "#     temp = pd.read_csv(os.path.join(path,'results/perms/%s_indiv_SV_allfts_5000perms.csv'%pt), index_col=0)\n",
    "#     for col in list(main.keys()):\n",
    "#         if col[:3] == pt:\n",
    "#             main[col] = temp[col]\n",
    "# main.to_csv(os.path.join(path,'results/perms/indiv_SV_allfts_5000perms.csv'), index=True)\n",
    "\n",
    "# model = 'indiv_SV_allfts_actfilter_5000perms'\n",
    "# files = listdir(os.path.join(path,'results/perms/temp'))\n",
    "# pts = np.sort(list(set([file[:3] for file in files])))\n",
    "# first = pd.read_csv(os.path.join(path,'results/perms/temp/002_%s.csv' % model), index_col=0)\n",
    "# main = pd.DataFrame(columns=[pt+var[3:] for pt in pts for var in first.keys() ])\n",
    "# for file in files:\n",
    "#     if file[4:-4] == model:\n",
    "#         dat = pd.read_csv(os.path.join(path,'results/perms/temp/%s' % file), index_col=0)\n",
    "#         for col in dat.keys():\n",
    "#             main[col] = dat[col]\n",
    "# main.to_csv(os.path.join(path,'results/perms/%s.csv' % model), index=True)\n",
    "# main\n",
    "\n",
    "# ## STILL TO DO: [002-014] for indiv SV 4fts actfilter (now filled with 0.5's)\n",
    "# main = pd.read_csv(os.path.join(path,'results/perms/indiv_SV_4fts_actfilter_5000perms.csv'), index_col=0)\n",
    "# model = 'indiv_SV_4fts_actfilter_5000perms'\n",
    "# pts = set([k[:3] for k in main.keys()])\n",
    "# pts_todo = ['002','012','013','014']\n",
    "# for pt in pts_todo:\n",
    "#     dat = pd.read_csv(os.path.join(path,'results/perms/temp/%s_%s.csv' % (pt,model)), index_col=0)\n",
    "#     for c in dat.keys():\n",
    "#         main[c] = dat[c]\n",
    "# #     for col in set([c[3:] for c in main.keys()]):\n",
    "# #         main[pt+col] = [0.5]*5000\n",
    "           \n",
    "# main.to_csv(os.path.join(path,'results/perms/%s.csv' % model), index=True)\n"
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}