3_INPUTS_age.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5ee962c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import warnings\n",
    "#from MyUtils import *\n",
    "import tensorflow as tf\n",
    "from tensorflow import keras\n",
    "from tensorflow.keras import layers\n",
    "import tensorflow_addons as tfa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72200df7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "ts = str(int(time.time()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2e5d975",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select the time-to-event time\n",
    "def selected_time(selected_time, time, event):\n",
    "    if (time > selected_time): return 0\n",
    "    else:\n",
    "        if (event == 1.0): return 1\n",
    "        else: return None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1446268",
   "metadata": {},
   "outputs": [],
   "source": [
    "def numerical_stages(stage):\n",
    "    if stage=='Stage I' or stage=='Stage IA':\n",
    "        numerical_stage=0.0\n",
    "    if stage=='Stage IB':\n",
    "        numerical_stage=5.0\n",
    "    if stage=='Stage II' or stage=='Stage IIA':\n",
    "        numerical_stage=10.0\n",
    "    if stage=='Stage IIB':\n",
    "        numerical_stage=15.0\n",
    "    if stage=='Stage III' or stage=='Stage IIIA':\n",
    "        numerical_stage=20.0\n",
    "    if stage=='Stage IIIB':\n",
    "        numerical_stage=25.0\n",
    "    if stage=='Stage IV':\n",
    "        numerical_stage=30.0\n",
    "    return numerical_stage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "515edd40",
   "metadata": {},
   "outputs": [],
   "source": [
    "survival_metric='OS'\n",
    "survival_time='OS.time'\n",
    "time=1800\n",
    "model_save_path='./saved_models/lung_model_cnn_{}_{}_{}.h5'.format(survival_metric,time,ts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75662c94",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "Y_type_L=pd.read_hdf('data/Preprocessed_Data/Lung.h5',key='sample_type')\n",
    "Y_clinical_L=pd.read_hdf('data/Preprocessed_Data/Lung.h5',key='clinical')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5e5e9fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_clinical_L.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dbdf8cb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_clinical_L=Y_clinical_L.loc[Y_type_L.pathology_status=='Tumor']\n",
    "Y_clinical_L.dropna(subset=[survival_metric, survival_time], inplace=True)\n",
    "Y_clinical_L.dropna(subset=['age_at_initial_pathologic_diagnosis'],inplace=True)\n",
    "Y_clinical_L.dropna(subset=['ajcc_pathologic_tumor_stage'],inplace=True)\n",
    "Y_clinical_L=Y_clinical_L[~Y_clinical_L['ajcc_pathologic_tumor_stage'].isin(['[Discrepancy]'])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ac1732d",
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_clinical_L['ajcc_pathologic_tumor_stage'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf867c30",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_stage = Y_clinical_L[['ajcc_pathologic_tumor_stage']].apply(lambda row: numerical_stages(row['ajcc_pathologic_tumor_stage']), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43110cdb",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_stage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e043437b",
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_clinical_categorical_L = Y_clinical_L[[survival_metric, survival_time]].apply(lambda row: selected_time(time, row[survival_time], row[survival_metric]), axis=1)\n",
    "Y_clinical_categorical_L.dropna(inplace=True)\n",
    "Y_clinical_categorical_L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac501762",
   "metadata": {},
   "outputs": [],
   "source": [
    "age_os=Y_clinical_L[['age_at_initial_pathologic_diagnosis','OS']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e0d6123",
   "metadata": {},
   "outputs": [],
   "source": [
    "age_os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "baa44965",
   "metadata": {},
   "outputs": [],
   "source": [
    "age=Y_clinical_L['age_at_initial_pathologic_diagnosis']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dce9bcb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "age.value_counts(ascending=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40fec32d",
   "metadata": {},
   "outputs": [],
   "source": [
    "age=np.array(age.loc[[s for s in Y_clinical_categorical_L.index]])-65"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bbff04f",
   "metadata": {},
   "outputs": [],
   "source": [
    "age=age.reshape(age.shape[0],1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fa96443",
   "metadata": {},
   "outputs": [],
   "source": [
    "age.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1aadc10",
   "metadata": {},
   "outputs": [],
   "source": [
    "stage=np.array(num_stage.loc[[s for s in Y_clinical_categorical_L.index]])\n",
    "stage=stage.reshape(stage.shape[0],1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58eba337",
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_clinical_categorical_L.value_counts(normalize=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94e4f286",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "#读取图像数据\n",
    "dir_name_p = 'images/genes_selected/log2tpm_MinMax/lung_with_keggpathway_1k/'\n",
    "dir_name_f = 'images/genes_selected/log2tpm_MinMax/lung_with_function_1k/'\n",
    "#dir_name = 'images/genes_selected/log2tpm_MinMax/lung_no_pathway/'\n",
    "X_p = np.array([np.load(dir_name_p+s+\".npy\") for s in Y_clinical_categorical_L.index]) \n",
    "X_f = np.array([np.load(dir_name_f+s+\".npy\") for s in Y_clinical_categorical_L.index]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35384a73",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(X_p.shape)\n",
    "print(X_f.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c41f356b",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "from sklearn.preprocessing import LabelEncoder\n",
    "y_L = LabelEncoder().fit_transform(Y_clinical_categorical_L)\n",
    "y_L=y_L.reshape(y_L.shape[0],1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90b85c1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "#划分训练集测试集\n",
    "from sklearn.model_selection import train_test_split\n",
    "X_p_train,X_p_test,X_f_train,X_f_test,age_train,age_test,stage_train,stage_test,Y_clinical_categorical_L_train,Y_clinical_categorical_L_test,y_train,y_test=train_test_split(X_p,X_f,age,stage,Y_clinical_categorical_L,y_L,test_size=0.2,random_state=126,stratify=y_L)\n",
    "print(X_p_train.shape)\n",
    "print(X_p_test.shape)\n",
    "print(X_f_train.shape)\n",
    "print(X_f_test.shape)\n",
    "print(age_train.shape)\n",
    "print(age_test.shape)\n",
    "print(stage_train.shape)\n",
    "print(stage_test.shape)\n",
    "print(y_train.shape)\n",
    "print(y_test.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8e19495",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "X_random_p_train=X_p_train.reshape(X_p_train.shape[0],26*26)\n",
    "X_random_f_train=X_f_train.reshape(X_f_train.shape[0],27*27)\n",
    "age_random_train=age_train\n",
    "stage_random_train=stage_train\n",
    "X_i_train=[[i] for i in range(len(y_train))]\n",
    "#X_smote.shape\n",
    "from imblearn.over_sampling import SMOTE as os\n",
    "\n",
    "\n",
    "X_ri_train, y_r_train = os(random_state=321).fit_resample(X_i_train,y_train)\n",
    "X_p_r_train=X_random_p_train[[i[0] for i in X_ri_train]]\n",
    "X_f_r_train=X_random_f_train[[i[0] for i in X_ri_train]]\n",
    "age_r_train=age_random_train[[i[0] for i in X_ri_train]]\n",
    "stage_r_train=stage_random_train[[i[0] for i in X_ri_train]]\n",
    "\n",
    "X_p_r_train=X_p_r_train.reshape(X_p_r_train.shape[0],26,26,1)\n",
    "X_f_r_train=X_f_r_train.reshape(X_f_r_train.shape[0],27,27,1)\n",
    "y_r_train=y_r_train.reshape(y_r_train.shape[0],1)\n",
    "print(X_p_r_train.shape)\n",
    "print(X_f_r_train.shape)\n",
    "print(age_r_train.shape)\n",
    "print(stage_r_train.shape)\n",
    "print(y_r_train.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd807540",
   "metadata": {},
   "outputs": [],
   "source": [
    "import keras_tuner as kt\n",
    "def build_dual_inputs_cnn(hp):\n",
    "    \n",
    "    function_input=keras.layers.Input(shape=(27,27,1))\n",
    "    pathway_input=keras.layers.Input(shape=(26,26,1))\n",
    "    age_input=keras.layers.Input(shape=(1,))\n",
    "    #stage_input=keras.layers.Input(shape=(1,))\n",
    "    \n",
    "    f_conv_1=layers.Conv2D(hp.Choice('filters1',[32,40,48,56,64]), \n",
    "                           kernel_size=(3,3), \n",
    "                           strides=1,\n",
    "                           padding='same',\n",
    "                           use_bias=False, \n",
    "                           activation='relu',\n",
    "                           kernel_initializer='he_uniform', \n",
    "                           name='f_conv_1',\n",
    "                           )(function_input)\n",
    "    \n",
    "    f_pooling1=layers.MaxPooling2D(pool_size=(3,3),\n",
    "                         strides=2, \n",
    "                         padding='same')(f_conv_1)\n",
    "    \n",
    "    f_bn1=layers.BatchNormalization()(f_pooling1)\n",
    "    \n",
    "    f_conv_2=layers.Conv2D(hp.Choice('filters2',[80,96,112,128]), \n",
    "                           kernel_size=(3,3), \n",
    "                           strides=1,\n",
    "                           padding='same',\n",
    "                           use_bias=False, \n",
    "                           activation='relu',\n",
    "                           kernel_initializer='he_uniform', \n",
    "                           name='f_conv_2',\n",
    "                           )(f_bn1)\n",
    "    \n",
    "    #f_padding_1=layers.ZeroPadding2D(3)(f_conv_1)\n",
    "    \n",
    "    f_pooling2=layers.MaxPooling2D(pool_size=(3,3),\n",
    "                         strides=2, \n",
    "                         padding='same')(f_conv_2) \n",
    "    \n",
    "    f_bn2=layers.BatchNormalization()(f_pooling2)\n",
    "    \n",
    "    f_fn1=layers.Dense(hp.Choice('units_f',[128,144,192,256]),activation='relu')(f_bn2)\n",
    "    \n",
    "    f_fn1=layers.Dropout(hp.Choice('drf',[0.1,0.2,0.3]))(f_fn1)\n",
    "    \n",
    "    p_conv_1=layers.Conv2D(hp.Choice('filters3',[32,40,48,56,64]), \n",
    "                           kernel_size=(3,3), \n",
    "                           strides=1,\n",
    "                           padding='same',\n",
    "                           use_bias=False, \n",
    "                           activation='relu',\n",
    "                           kernel_initializer='he_uniform', \n",
    "                           name='p_conv_1',\n",
    "                           )(pathway_input)\n",
    "    \n",
    "    p_pooling1=layers.MaxPooling2D(pool_size=(3,3),\n",
    "                         strides=2, \n",
    "                         padding='same')(p_conv_1)\n",
    "    \n",
    "    p_bn1=layers.BatchNormalization()(p_pooling1)\n",
    "    \n",
    "    p_conv_2=layers.Conv2D(hp.Choice('filters4',[80,96,112,128]), \n",
    "                           kernel_size=(3,3), \n",
    "                           strides=1,\n",
    "                           padding='same',\n",
    "                           use_bias=False, \n",
    "                           activation='relu',\n",
    "                           kernel_initializer='he_uniform', \n",
    "                           name='p_conv_2',\n",
    "                           )(p_bn1)\n",
    "    \n",
    "    #f_padding_1=layers.ZeroPadding2D(3)(f_conv_1)\n",
    "    \n",
    "    p_pooling2=layers.MaxPooling2D(pool_size=(3,3),\n",
    "                         strides=2, \n",
    "                         padding='same')(p_conv_2)\n",
    "    \n",
    "    p_bn2=layers.BatchNormalization()(p_pooling2)\n",
    "    \n",
    "    p_fn1=layers.Dense(hp.Choice('units_p',[128,144,192,256]),activation='relu')(p_bn2)\n",
    "    \n",
    "    p_fn1=layers.Dropout(hp.Choice('drp',[0.1,0.2,0.3]))(p_fn1)\n",
    "    \n",
    "    x=layers.concatenate([f_fn1,p_fn1])\n",
    "    \n",
    "    x=layers.Flatten()(x)\n",
    "    \n",
    "    #x=layers.concatenate([x,age_input])\n",
    "    \n",
    "    #x=layers.GlobalAveragePooling2D()(x)\n",
    "    \n",
    "    #x=layers.Dense(hp.Choice('units1',[32,64,96]),activation='relu')(x)\n",
    "    \n",
    "    #x=layers.Dropout(hp.Choice('dr1',[0.1,0.2,0.3]))(x)\n",
    "    \n",
    "    x=layers.Dense(hp.Choice('units2',[64,128,144,192,256]),activation='relu')(x)\n",
    "    \n",
    "    x=layers.Dropout(hp.Choice('dr2',[0.3,0.4,0.5]))(x)\n",
    "    \n",
    "    #x=layers.Flatten()(x)\n",
    "    \n",
    "    x=layers.concatenate([x,age_input])\n",
    "    \n",
    "    x=layers.Dense(hp.Choice('units3',[32,64,128]),activation='relu')(x)\n",
    "    \n",
    "    x=layers.Dropout(hp.Choice('dr3',[0.3,0.4,0.5]))(x)\n",
    "    \n",
    "    output=layers.Dense(1,activation='sigmoid')(x)\n",
    "    \n",
    "    model=keras.Model(inputs=[function_input,pathway_input,age_input],outputs=output)\n",
    "    \n",
    "    optimizer = tfa.optimizers.AdamW(learning_rate=hp.Choice('lr',[0.0005,0.001,0.002,0.003]), weight_decay=0.0001)\n",
    "\n",
    "    model.compile(\n",
    "        optimizer=optimizer,\n",
    "        loss=keras.losses.BinaryCrossentropy(\n",
    "            from_logits=False, label_smoothing=0.1\n",
    "        ),\n",
    "        metrics=[\n",
    "#            keras.metrics.BinaryAccuracy(name=\"accuracy\"),\n",
    "#            keras.metrics.BinaryAccuracy(\n",
    "#                name=\"binary_accuracy\", dtype=None, threshold=0.5\n",
    "#            ),\n",
    "            keras.metrics.AUC(name='auc',from_logits=False),\n",
    "#            keras.metrics.Precision(name='precision'),\n",
    "#            keras.metrics.Recall(name='recall'),  \n",
    "#            tfa.metrics.F1Score(num_classes,threshold=0.5,name='f1_score')\n",
    "        ],\n",
    "    )\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d849f716",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "tuner = kt.BayesianOptimization(\n",
    "        build_dual_inputs_cnn,\n",
    "        objective=kt.Objective(\"val_auc\", direction=\"max\"),\n",
    "        max_trials=100,\n",
    "        overwrite=True)\n",
    "stop_early = tf.keras.callbacks.EarlyStopping(monitor='val_auc',mode='max', patience=6)\n",
    "tuner.search([X_f_r_train,X_p_r_train,age_r_train], y_r_train,epochs=50, validation_data=([X_f_test,X_p_test,age_test], y_test),callbacks=[stop_early])\n",
    "best_model = tuner.get_best_models()[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42c8d4a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_hyperparameters = tuner.get_best_hyperparameters(1)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9d6a097",
   "metadata": {},
   "outputs": [],
   "source": [
    "model=best_model\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e96e8b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "from tensorflow.keras.utils import plot_model\n",
    "plot_model(best_model, to_file='model1.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11b5157b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import roc_curve\n",
    "import sklearn.metrics\n",
    "def train_model(model):\n",
    "    #set dynamic learning rate\n",
    "    initial_learning_rate = 0.0001\n",
    "    '''\n",
    "    lr_schedule = tf.keras.optimizers.schedules.ExponentialDecay(\n",
    "        initial_learning_rate,\n",
    "        decay_steps=100000,\n",
    "        decay_rate=0.96,\n",
    "        staircase=True)\n",
    "    '''\n",
    "    optimizer = tfa.optimizers.AdamW(learning_rate=0.00001, weight_decay=0.0001)\n",
    "    #optimizer=keras.optimizers.Adam(learning_rate=0.001)\n",
    "    model.compile(\n",
    "        optimizer=optimizer,\n",
    "        loss=keras.losses.BinaryCrossentropy(\n",
    "            from_logits=False, label_smoothing=0.1\n",
    "        ),\n",
    "        metrics=[\n",
    "#            keras.metrics.BinaryAccuracy(name=\"accuracy\"),\n",
    "#            keras.metrics.BinaryAccuracy(\n",
    "#                name=\"binary_accuracy\", dtype=None, threshold=0.5\n",
    "#            ),\n",
    "            keras.metrics.AUC(name='auc',from_logits=False,num_thresholds=500),\n",
    "#            keras.metrics.Precision(name='precision'),\n",
    "#            keras.metrics.Recall(name='recall'),  \n",
    "#            tfa.metrics.F1Score(num_classes,threshold=0.5,name='f1_score')\n",
    "        ],\n",
    "    )\n",
    "\n",
    "    checkpoint_filepath = \"./tmp/checkpoint\"\n",
    "    checkpoint_callback = keras.callbacks.ModelCheckpoint(\n",
    "        checkpoint_filepath,\n",
    "        monitor=\"val_auc\",\n",
    "        mode='max',\n",
    "        save_best_only=True,\n",
    "        save_weights_only=True,\n",
    "    )\n",
    "    \n",
    "    early_stop = keras.callbacks.EarlyStopping(monitor='val_auc',mode='max', patience=10)\n",
    "\n",
    "    history = model.fit(\n",
    "        x=([X_f_r_train,X_p_r_train,age_r_train]),\n",
    "        y=y_r_train,\n",
    "        batch_size=32,\n",
    "        epochs=150,\n",
    "        validation_data=([X_f_test,X_p_test,age_test], y_test),\n",
    "        callbacks=[checkpoint_callback,early_stop],\n",
    "    )\n",
    "\n",
    "    model.load_weights(checkpoint_filepath)\n",
    "    y_pred=model.predict([X_f_test,X_p_test,age_test])\n",
    "    auc=sklearn.metrics.roc_auc_score(y_test,y_pred)\n",
    "    Y_clinical_categorical_pred=pd.DataFrame(Y_clinical_categorical_L_test)\n",
    "    Y_clinical_categorical_pred['prediction']=y_pred\n",
    "    fpr, tpr, thresholds = roc_curve(y_test, y_pred)\n",
    "    J = tpr - fpr\n",
    "    index_max = np.argmax(J)\n",
    "    best_thresh = thresholds[index_max]\n",
    "    \n",
    "    tp=Y_clinical_categorical_pred[(Y_clinical_categorical_pred[0]==1.0) & (Y_clinical_categorical_pred['prediction']>=best_thresh)].shape[0]\n",
    "    fp=Y_clinical_categorical_pred[(Y_clinical_categorical_pred[0]==1.0) & (Y_clinical_categorical_pred['prediction']<best_thresh)].shape[0]\n",
    "    tn=Y_clinical_categorical_pred[(Y_clinical_categorical_pred[0]==0.0) & (Y_clinical_categorical_pred['prediction']<=best_thresh)].shape[0]\n",
    "    fn=Y_clinical_categorical_pred[(Y_clinical_categorical_pred[0]==0.0) & (Y_clinical_categorical_pred['prediction']>best_thresh)].shape[0]\n",
    "    accuracy=(tp+tn)/(tp+fp+tn+fn)\n",
    "    precision=tp/(tp+fp)\n",
    "    recall=tp/(tp+fn)\n",
    "    f1_score=2*(precision*recall)/(precision+recall)\n",
    "    \n",
    "    #    _, accuracy, binary_accuracy,auc_roc = model.evaluate(X_test, y_test)\n",
    "    #_, accuracy,auc,precision,recall = model.evaluate([X_f_test,X_p_test,age_test,stage_test], y_test)\n",
    "   \n",
    "    print('Best Threshold=%f' % (best_thresh))\n",
    "    print(f\"Test accuracy: {round(accuracy,4)}\")\n",
    "    print(f\"Test AUC: {round(auc, 4)}\")\n",
    "    print(f\"Test precision: {round(precision, 4)}\")\n",
    "    print(f\"Test recall: {round(recall, 4)}\")\n",
    "    print(f\"Test f1_score: {round(f1_score,4)}\")\n",
    "    \n",
    "    return history "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03dfd409",
   "metadata": {},
   "outputs": [],
   "source": [
    "model=best_model\n",
    "history=train_model(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c42350ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "from sklearn.model_selection import train_test_split\n",
    "from imblearn.over_sampling import SMOTE as os\n",
    "from sklearn.model_selection import StratifiedKFold\n",
    "from sklearn.metrics import roc_curve\n",
    "import sklearn.metrics\n",
    "accuracy_list=[]\n",
    "auc_list=[]\n",
    "precision_list=[]\n",
    "recall_list=[]\n",
    "f1_score_list=[]\n",
    "init_random_state=255\n",
    "for i in range(0,50):\n",
    "    X_p_train,X_p_test,X_f_train,X_f_test,age_train,age_test,stage_train,stage_test,Y_clinical_categorical_L_train,Y_clinical_categorical_L_test,y_train,y_test=train_test_split(X_p,X_f,age,stage,Y_clinical_categorical_L,y_L,test_size=0.2,random_state=init_random_state+i,stratify=y_L)\n",
    "    print('=============================================================================================================')\n",
    "    print('* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *')\n",
    "    print(' * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ')\n",
    "    print('* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *')\n",
    "    print(' * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ')\n",
    "    print('* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *')\n",
    "    print('=============================================================================================================')\n",
    "    print('current random state is: {}'.format(init_random_state+i))\n",
    "    print('This is the {}th fold'.format(i+1))\n",
    "    \n",
    "    X_random_p_train=X_p_train.reshape(X_p_train.shape[0],26*26)\n",
    "    X_random_f_train=X_f_train.reshape(X_f_train.shape[0],27*27)\n",
    "    age_random_train=age_train\n",
    "    stage_random_train=stage_train\n",
    "    X_i_train=[[i] for i in range(len(y_train))]\n",
    "    X_ri_train, y_r_train = os(random_state=321+i).fit_resample(X_i_train,y_train)\n",
    "    X_p_r_train=X_random_p_train[[i[0] for i in X_ri_train]]\n",
    "    X_f_r_train=X_random_f_train[[i[0] for i in X_ri_train]]\n",
    "    age_r_train=age_random_train[[i[0] for i in X_ri_train]]\n",
    "    stage_r_train=stage_random_train[[i[0] for i in X_ri_train]]\n",
    "\n",
    "    X_p_r_train=X_p_r_train.reshape(X_p_r_train.shape[0],26,26,1)\n",
    "    X_f_r_train=X_f_r_train.reshape(X_f_r_train.shape[0],27,27,1)\n",
    "    y_r_train=y_r_train.reshape(y_r_train.shape[0],1)\n",
    "    \n",
    "    model=build_dual_inputs_cnn(best_hyperparameters)\n",
    "    \n",
    "    lr_schedule = keras.optimizers.schedules.ExponentialDecay(\n",
    "    initial_learning_rate=3e-3,\n",
    "    decay_steps=10,\n",
    "    decay_rate=0.99)\n",
    "\n",
    "    \n",
    "    optimizer = tfa.optimizers.AdamW(learning_rate=lr_schedule, weight_decay=0.0001)\n",
    "    #optimizer=keras.optimizers.Adam(learning_rate=0.001)\n",
    "    model.compile(\n",
    "        optimizer=optimizer,\n",
    "        loss=keras.losses.BinaryCrossentropy(\n",
    "            from_logits=False, label_smoothing=0.1\n",
    "        ),\n",
    "        metrics=[\n",
    "#            keras.metrics.BinaryAccuracy(name=\"accuracy\"),\n",
    "#            keras.metrics.BinaryAccuracy(\n",
    "#                name=\"binary_accuracy\", dtype=None, threshold=0.5\n",
    "#            ),\n",
    "            keras.metrics.AUC(name='auc',from_logits=False,num_thresholds=500),\n",
    "#            keras.metrics.Precision(name='precision'),\n",
    "#            keras.metrics.Recall(name='recall'),  \n",
    "#            tfa.metrics.F1Score(num_classes,threshold=0.5,name='f1_score')\n",
    "        ],\n",
    "    )\n",
    "\n",
    "    checkpoint_filepath = \"./tmp/checkpoint\"\n",
    "    checkpoint_callback = keras.callbacks.ModelCheckpoint(\n",
    "        checkpoint_filepath,\n",
    "        monitor=\"val_auc\",\n",
    "        mode='max',\n",
    "        save_best_only=True,\n",
    "        save_weights_only=True,\n",
    "    )\n",
    "    \n",
    "    early_stop = keras.callbacks.EarlyStopping(monitor='val_auc',mode='max', patience=30)\n",
    "\n",
    "    history = model.fit(\n",
    "        x=([X_f_r_train,X_p_r_train,age_r_train]),\n",
    "        y=y_r_train,\n",
    "        batch_size=32,\n",
    "        epochs=150,\n",
    "        validation_data=([X_f_test,X_p_test,age_test], y_test),\n",
    "        callbacks=[checkpoint_callback,early_stop],\n",
    "    )\n",
    "\n",
    "    model.load_weights(checkpoint_filepath)\n",
    "#    _, accuracy, binary_accuracy,auc_roc = model.evaluate(X_test, y_test)\n",
    "    #_, accuracy,auc,precision,recall = model.evaluate([X_f_test,X_p_test,age_test,stage_test], y_test)\n",
    "    y_pred=model.predict([X_f_test,X_p_test,age_test])\n",
    "    auc=sklearn.metrics.roc_auc_score(y_test,y_pred)\n",
    "    Y_clinical_categorical_pred=pd.DataFrame(Y_clinical_categorical_L_test)\n",
    "    Y_clinical_categorical_pred['prediction']=y_pred\n",
    "    \n",
    "    fpr, tpr, thresholds = roc_curve(y_test, y_pred)\n",
    "    J = tpr - fpr\n",
    "    index_max = np.argmax(J)\n",
    "    best_thresh = thresholds[index_max]\n",
    "    \n",
    "    tp=Y_clinical_categorical_pred[(Y_clinical_categorical_pred[0]==1.0) & (Y_clinical_categorical_pred['prediction']>=best_thresh)].shape[0]\n",
    "    fp=Y_clinical_categorical_pred[(Y_clinical_categorical_pred[0]==1.0) & (Y_clinical_categorical_pred['prediction']<best_thresh)].shape[0]\n",
    "    tn=Y_clinical_categorical_pred[(Y_clinical_categorical_pred[0]==0.0) & (Y_clinical_categorical_pred['prediction']<=best_thresh)].shape[0]\n",
    "    fn=Y_clinical_categorical_pred[(Y_clinical_categorical_pred[0]==0.0) & (Y_clinical_categorical_pred['prediction']>best_thresh)].shape[0]\n",
    "    \n",
    "    accuracy=(tp+tn)/(tp+fp+tn+fn)\n",
    "    precision=tp/(tp+fp)\n",
    "    recall=tp/(tp+fn)\n",
    "    f1_score=2*(precision*recall)/(precision+recall) \n",
    "    \n",
    "   \n",
    "    #print(f\"Test accuracy: {round(accuracy * 100, 2)}%\")\n",
    "#    print(f\"Test binary accuracy: {round(binary_accuracy * 100, 2)}%\")\n",
    "    #print(f\"Test AUC: {round(auc, 4)}\")\n",
    "    #print(f\"Test precision: {round(precision * 100, 2)}%\")\n",
    "    #print(f\"Test recall: {round(recall * 100, 2)}%\")\n",
    "    #print(f\"Test f1_score: {round(f1_score*100,2)}%\")\n",
    "    accuracy_list.append(accuracy)\n",
    "    auc_list.append(auc)\n",
    "    precision_list.append(precision)\n",
    "    recall_list.append(recall)\n",
    "    f1_score_list.append(f1_score)\n",
    "    \n",
    "print('mean accuracy={}±{}'.format(np.mean(accuracy_list),np.std(accuracy_list)))\n",
    "print('mean auc={}±{}'.format(np.mean(auc_list),np.std(auc_list)))\n",
    "print('mean precision={}±{}'.format(np.mean(precision_list),np.std(precision_list)))\n",
    "print('mean recall={}±{}'.format(np.mean(recall_list),np.std(recall_list)))\n",
    "print('mean f1_score={}±{}'.format(np.mean(f1_score_list),np.std(f1_score_list)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c311886",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "f1_score_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3c216e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sklearn.metrics\n",
    "y_pred=model.predict([X_f_test,X_p_test,age_test,stage_test])\n",
    "auc=sklearn.metrics.roc_auc_score(y_test,y_pred)\n",
    "print(auc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2af08682",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3c4f35e",
   "metadata": {},
   "outputs": [],
   "source": [
    "f1_score_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c7db19e",
   "metadata": {},
   "outputs": [],
   "source": [
    "predictions=model.predict([X_f_test,X_p_test,age_test,stage_test])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "793476e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "#perform the youden's J statistics best threshold selection\n",
    "from sklearn.metrics import roc_curve\n",
    "fpr, tpr, thresholds = roc_curve(y_test, y_pred)\n",
    "J = tpr - fpr\n",
    "index_max = np.argmax(J)\n",
    "best_thresh = thresholds[index_max]\n",
    "print('Best Threshold=%f' % (best_thresh))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "940906a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_clinical_categorical_pred=pd.DataFrame(Y_clinical_categorical_L_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05ca3e17",
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_clinical_categorical_pred['prediction']=predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aaaa5fb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_clinical_categorical_pred[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53e47ec4",
   "metadata": {},
   "outputs": [],
   "source": [
    "tp=Y_clinical_categorical_pred[(Y_clinical_categorical_pred[0]==1.0) & (Y_clinical_categorical_pred['prediction']>=best_thresh)].shape[0]\n",
    "fp=Y_clinical_categorical_pred[(Y_clinical_categorical_pred[0]==1.0) & (Y_clinical_categorical_pred['prediction']<best_thresh)].shape[0]\n",
    "tn=Y_clinical_categorical_pred[(Y_clinical_categorical_pred[0]==0.0) & (Y_clinical_categorical_pred['prediction']<=best_thresh)].shape[0]\n",
    "fn=Y_clinical_categorical_pred[(Y_clinical_categorical_pred[0]==0.0) & (Y_clinical_categorical_pred['prediction']>best_thresh)].shape[0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7ca2811",
   "metadata": {},
   "outputs": [],
   "source": [
    "real_acc=(tp+tn)/95\n",
    "real_precision=tp/(tp+fp)\n",
    "real_recall=tp/(tp+fn)\n",
    "real_f1_score=2*(real_precision*real_recall)/(real_precision+real_recall)\n",
    "real_acc,real_precision,real_recall,real_f1_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "edd0966f",
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_positive=Y_clinical_categorical_pred.loc[Y_clinical_categorical_pred.prediction>=0.626445]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63fab812",
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_negative=Y_clinical_categorical_pred.loc[Y_clinical_categorical_pred.prediction<0.626445]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4cd902d",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinical_positive=Y_clinical_L.loc[[s for s in pred_positive.index]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e563ca51",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinical_negative=Y_clinical_L.loc[[s for s in pred_negative.index]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad3b4e43",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinical_positive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d143c5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinical_negative"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb3b3e5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from lifelines import KaplanMeierFitter\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "timeline=[]\n",
    "for i in range(0,1830):\n",
    "    timeline.append(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22aab540",
   "metadata": {},
   "outputs": [],
   "source": [
    "from lifelines.statistics import logrank_test\n",
    "from lifelines.utils import restricted_mean_survival_time\n",
    "from lifelines.plotting import rmst_plot\n",
    "\n",
    "T_p=clinical_positive['OS.time']\n",
    "E_p=clinical_positive['OS']\n",
    "kmf_positive = KaplanMeierFitter()\n",
    "kmf_positive.fit(durations=T_p, event_observed=E_p,timeline=timeline,label='Predicted Positive')\n",
    "ax=kmf_positive.plot()\n",
    "\n",
    "T_n=clinical_negative['OS.time']\n",
    "E_n=clinical_negative['OS']\n",
    "kmf_negative=KaplanMeierFitter()\n",
    "kmf_negative.fit(durations=T_n, event_observed=E_n,timeline=timeline,label='Predicted Negative')\n",
    "ax=kmf_negative.plot(ax=ax)\n",
    "\n",
    "results = logrank_test(T_p, T_n, event_observed_A=E_p, event_observed_B=E_n)\n",
    "'''\n",
    "rmst_positive= restricted_mean_survival_time(kmf_positive, t=900)\n",
    "rmst_negative= restricted_mean_survival_time(kmf_negative, t=900)\n",
    "\n",
    "ax = plt.subplot(311)\n",
    "rmst_plot(kmf_negative, t=900, ax=ax)\n",
    "\n",
    "\n",
    "ax = plt.subplot(312)\n",
    "rmst_plot(kmf_positive, t=900, ax=ax)\n",
    "\n",
    "\n",
    "ax = plt.subplot(313)\n",
    "rmst_plot(kmf_negative, model2=kmf_positive, t=900, ax=ax)\n",
    "'''\n",
    "plt.legend(loc='lower left', bbox_to_anchor=(0.01,0.01),ncol=1,fancybox=False,shadow=False)\n",
    "plt.title('SiameseSurvNet p-value: {}'.format(results.p_value))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7fb48b61",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax.get_figure().savefig(\"figures/SiameseSurvNet_kpcurve.svg\",dpi=600,format='svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a54f89d1",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python3.8(deepSurvival)",
   "language": "python",
   "name": "deepsurvival"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}