supervised_distance_feature_elimination.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25ca4ba2-058e-41dc-a6be-e6f1049e22c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import xgboost\n",
    "import shap\n",
    "import matplotlib.pylab as pl\n",
    "import pandas as pd\n",
    "from sklearn.model_selection import train_test_split\n",
    "import numpy as np\n",
    "import os\n",
    "from sklearn.metrics import roc_auc_score, average_precision_score\n",
    "import argparse\n",
    "from feature_selector import FeatureSelector\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "import pickle\n",
    "import sklearn\n",
    "import scipy as sp\n",
    "from scipy.cluster.hierarchy import complete, fcluster\n",
    "from collections import defaultdict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fec600fc-ada1-4eb6-85ae-8200d49b79b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def label(permth, mortstat, month):\n",
    "    if permth > month:\n",
    "        return 0\n",
    "    else:\n",
    "        if mortstat == 1:\n",
    "            return 1\n",
    "        else:\n",
    "            return 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee5b9732-aa8a-474a-aa11-528364f09820",
   "metadata": {},
   "outputs": [],
   "source": [
    "def bootstrap_ci(y_pre, y_label, sample_size, repetitions = 1000, alpha = 0.05): \n",
    "    y_pre = np.array(y_pre)\n",
    "    y_label = np.array(y_label)\n",
    "    \n",
    "    auc = []\n",
    "    ap = []\n",
    "    for i in range(repetitions):\n",
    "        np.random.seed(i)\n",
    "        idx = list(np.random.choice(len(y_pre), replace = True, size = sample_size))\n",
    "        y_pre_bootstrap = y_pre[idx]\n",
    "        y_label_bootstrap = y_label[idx]\n",
    "        auc.append(roc_auc_score(y_label_bootstrap, y_pre_bootstrap))\n",
    "        ap.append(average_precision_score(y_label_bootstrap, y_pre_bootstrap))\n",
    "    # confidence interval\n",
    "    left_auc = np.percentile(auc, alpha/2*100)\n",
    "    right_auc = np.percentile(auc, 100-alpha/2*100)\n",
    "    left_ap = np.percentile(ap, alpha/2*100)\n",
    "    right_ap = np.percentile(ap, 100-alpha/2*100)\n",
    "    # point estimate\n",
    "    print('average AUROC', np.mean(auc))\n",
    "    print((1-alpha)*100,'%','confidence interval for the AUROC:', (round(left_auc,4), round(right_auc,4)))\n",
    "    print('average AP', np.mean(ap))\n",
    "    print((1-alpha)*100,'%','confidence interval for the AP:', (round(left_ap,4), round(right_ap,4)))\n",
    "    return auc, left_auc, right_auc, ap, left_ap, right_ap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63ecdc12-4baa-4ee5-be9c-d971885aea2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "year_num = 5\n",
    "path = './model/supervised_distance_feature_elimination/'\n",
    "if not os.path.isdir(path):\n",
    "    os.mkdir(path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d858c3d-8406-4937-9b30-6d890d4654a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = pd.read_csv('./data/NHANES/NHANES.csv')\n",
    "if str(year_num)+'_year_label' not in X.columns:\n",
    "    X[str(year_num)+'_year_label'] = X.apply(lambda x: label(x['permth_int'], x['mortstat'], 12*int(year_num)), axis=1)\n",
    "    \n",
    "X = X[X[str(year_num)+'_year_label']!=2]\n",
    "y = X[str(year_num)+'_year_label']\n",
    "\n",
    "if int(year_num) not in [1,2,3,4,5]:\n",
    "    X = X.drop([str(year_num)+'_year_label'], axis=1)\n",
    "\n",
    "mortstat = X['mortstat']\n",
    "permth_int = X['permth_int']\n",
    "drop_list = [\"mortstat\", \"permth_int\", '1_year_label', '2_year_label', '3_year_label', '4_year_label', '5_year_label']\n",
    "X = X.drop(drop_list, axis=1)\n",
    "\n",
    "X = X.drop(['Demographics_ReleaseCycle'], axis=1)\n",
    "\n",
    "fea_list = pd.read_csv('./data/NHANES/NHANES_feature_list.csv')\n",
    "nominal_fea = fea_list[fea_list['Nominal']==1]['Type_Short_Name'].tolist()\n",
    "nominal_fea = list(set(nominal_fea) & set(X.columns))\n",
    "X = pd.get_dummies(X, columns=nominal_fea, drop_first=True)\n",
    "\n",
    "display_name = pd.read_csv('./data/NHANES/NHANES_feature_list_Display_name.csv')\n",
    "display_col=[]\n",
    "for col in X.columns:\n",
    "    display_col.append(list(display_name.loc[display_name['Type_Short_Name']==col, 'Display_Name'])[0])\n",
    "col_dict = dict(zip(X.columns, display_col))\n",
    "\n",
    "print(X.columns)\n",
    "print('After encoding', X.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bed81de-120b-4637-b806-ac8b749940e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(X.columns)\n",
    "print(X.shape)\n",
    "print('# samples: ', X.shape[0])\n",
    "print('# positive samples: ', sum(y==1))\n",
    "print('# negative samples: ', sum(y==0))\n",
    "print('# features: ', X.shape[1])  \n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=7)\n",
    "X_train, X_val, y_train, y_val = train_test_split(X_train, y_train, test_size=0.2, random_state=7)\n",
    "\n",
    "y_train = np.array(y_train); y_test = np.array(y_test); y_val = np.array(y_val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88622dc8-469f-4723-9bf8-d770cef48053",
   "metadata": {},
   "outputs": [],
   "source": [
    "clustering = shap.utils.hclust(X_train, y_train, random_state=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42a8fdb6-cb36-4b0d-aec0-dd65877c8921",
   "metadata": {},
   "outputs": [],
   "source": [
    "maxclust_list = [151] + [i for i in range(145, 50, -5)] + [53, 1]\n",
    "X_train_all = X_train.copy()\n",
    "X_test_all = X_test.copy()\n",
    "X_val_all = X_val.copy()\n",
    "\n",
    "features_ranking_dict = {}\n",
    "auc_dict = {}\n",
    "ap_dict = {}\n",
    "left_auc_dict = {}\n",
    "right_auc_dict = {}\n",
    "left_ap_dict = {}\n",
    "right_ap_dict = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "635afd9f-c8a5-431d-bb4d-2bcc5188c569",
   "metadata": {},
   "outputs": [],
   "source": [
    "features_input = {}\n",
    "ranked_features = list(X_train_all.columns)\n",
    "feature_temp = ranked_features\n",
    "for maxclust in maxclust_list:\n",
    "    print('maxclust: ', maxclust)\n",
    "    if maxclust == X_train_all.shape[1]:\n",
    "        print('# features: ', maxclust)\n",
    "        X_train = X_train_all\n",
    "        X_test = X_test_all\n",
    "        X_val = X_val_all\n",
    "        cluster_num = maxclust\n",
    "        features_input[cluster_num] = feature_temp\n",
    "    else:\n",
    "        cluster_temp = fcluster(clustering, maxclust, criterion='maxclust')\n",
    "        cluster_num = max(cluster_temp)\n",
    "        if (max(cluster_temp) == len(ranked_features)):\n",
    "            features_input[cluster_num] = feature_temp\n",
    "            continue\n",
    "        cluster_dict = defaultdict(list)\n",
    "        cluster_rank_dict = defaultdict(list)\n",
    "        for fea, cluster_id in zip(X_train_all.columns, cluster_temp):\n",
    "            if fea in ranked_features:\n",
    "                cluster_dict[cluster_id].append(fea)\n",
    "                cluster_rank_dict[cluster_id].append(ranked_features.index(fea))\n",
    "        feature_temp = []\n",
    "        for cluster_id in cluster_dict:\n",
    "            feature_temp.append(cluster_dict[cluster_id][np.argmin(cluster_rank_dict[cluster_id])])\n",
    "        print('# features: ', len(feature_temp))\n",
    "        X_train = X_train_all.loc[:, feature_temp]\n",
    "        X_test = X_test_all.loc[:, feature_temp]\n",
    "        X_val = X_val_all.loc[:, feature_temp]\n",
    "        features_input[cluster_num] = feature_temp\n",
    "    xlf = xgboost.XGBClassifier(n_estimators=1000, max_depth=4, subsample=0.5, min_child_weight=3, objective='binary:logistic', random_state=7)\n",
    "    xlf.fit(X_train, y_train, eval_set = [(X_val, y_val)], early_stopping_rounds=100, verbose=False)\n",
    "    model_train = xlf\n",
    "    pickle.dump(model_train, open(path+\"model_\"+str(cluster_num)+\".pickle.dat\", \"wb\"))\n",
    "    y_pre = model_train.predict_proba(X_test)[:, 1]\n",
    "    auc, left_auc, right_auc, ap, left_ap, right_ap = bootstrap_ci(y_pre, y_test, len(y_test), repetitions = 1000, alpha = 0.05)\n",
    "    auc_dict[cluster_num] = auc\n",
    "    ap_dict[cluster_num] = ap\n",
    "    left_auc_dict[cluster_num] = left_auc\n",
    "    right_auc_dict[cluster_num] = right_auc\n",
    "    left_ap_dict[cluster_num] = left_ap\n",
    "    right_ap_dict[cluster_num] = right_ap    \n",
    "\n",
    "\n",
    "    if len(X_train)>=5000:\n",
    "        back_data = X_train.sample(n=5000, random_state=428)\n",
    "    else:\n",
    "        back_data = X_train\n",
    "    if len(X_test)>=2000:\n",
    "        fore_data = X_test.sample(n=2000, random_state=528)\n",
    "        fore_data_label = pd.DataFrame(y_test).sample(n=20, random_state=528)\n",
    "    else:\n",
    "        fore_data = X_test\n",
    "        fore_data_label = pd.DataFrame(y_test)\n",
    "\n",
    "    explainer = shap.TreeExplainer(model_train, data=back_data)\n",
    "    shap_values = explainer.shap_values(fore_data)\n",
    "    ranked_features = list(X_train.columns[np.argsort(-np.sum(np.abs(shap_values), axis=0))])\n",
    "\n",
    "    features_ranking_dict[cluster_num] = ranked_features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ff74cc5-c036-4f91-8ed8-d1465fbccaf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "pickle.dump(auc_dict, open(path+\"auc_dict.pickle.dat\", \"wb\"))\n",
    "pickle.dump(ap_dict, open(path+\"ap_dict.pickle.dat\", \"wb\"))\n",
    "pickle.dump(left_auc_dict, open(path+\"left_auc_dict.pickle.dat\", \"wb\"))\n",
    "pickle.dump(right_auc_dict, open(path+\"right_auc_dict.pickle.dat\", \"wb\"))\n",
    "pickle.dump(left_ap_dict, open(path+\"left_ap_dict.pickle.dat\", \"wb\"))\n",
    "pickle.dump(right_ap_dict, open(path+\"right_ap_dict.pickle.dat\", \"wb\"))\n",
    "pickle.dump(features_ranking_dict, open(path+\"features_ranking_dict.pickle.dat\", \"wb\"))\n",
    "pickle.dump(features_input, open(path+\"features_input.pickle.dat\", \"wb\"))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}