analysis/hospitalisation_prediction_urti.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Predcition of Hospital Admission Related to Infections\n",
    "\n",
    "In this notebook, we develop a Cox proportional hazard regression model to predict the risk of hospital admission for common infections including urinary tract infection (urti), upper respiratory tract infection (URTI), lower respiratory tract infection (LRTI), sinusitis, otitis media or middle ear infection (ot media), and ear infection or otitis externa (ot externa)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 244,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from matplotlib import pyplot\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import math\n",
    "import numpy as np\n",
    "import os\n",
    "import glob\n",
    "import gzip\n",
    "from matplotlib.ticker import PercentFormatter\n",
    "from patsy import dmatrices\n",
    "from lifelines import CoxPHFitter\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "# import sklearn\n",
    "# from sklearn.model_selection import train_test_split\n",
    "from sklearn.impute import SimpleImputer\n",
    "from sklearn.experimental import enable_iterative_imputer\n",
    "from sklearn.impute import IterativeImputer\n",
    "from sklearn.impute import KNNImputer\n",
    "import pickle\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "from lifelines.utils import k_fold_cross_validation\n",
    "from lifelines.utils import concordance_index\n",
    "from lifelines.calibration import survival_probability_calibration\n",
    "import io\n",
    "import sys\n",
    "from contextlib import redirect_stdout\n",
    "# import miceforest as mf\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 245,
   "metadata": {},
   "outputs": [],
   "source": [
    "#os.chdir(\"..\")\n",
    "current_working_dir = os.getcwd()\n",
    "# print(current_working_dir)\n",
    "\n",
    "# !python --version\n",
    "# import sklearn\n",
    "# print(sklearn.__version__)\n",
    "# import lifelines\n",
    "# print(lifelines.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 246,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No Warning Shown\n"
     ]
    }
   ],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "warnings.warn('DelftStack')\n",
    "warnings.warn('Do not show this message')\n",
    "print(\"No Warning Shown\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0- Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 247,
   "metadata": {},
   "outputs": [],
   "source": [
    "# def makeFilePath(pathFromWorkspace, fileName):\n",
    "#     workingdir = \"/workspace/\" if (os.name != \"nt\") else \"\"\n",
    "#     filepath = workingdir+pathFromWorkspace+fileName\n",
    "\n",
    "#     return filepath"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 248,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_histplot(title, df, ax):\n",
    "    ax.hist(df, color = 'purple', edgecolor = 'black', bins=20)\n",
    "    ax.set_title(title)\n",
    "    ax.set_ylabel('Frequency')\n",
    "    ax.xaxis.set_tick_params(which='both', labelbottom=True, labeltop=False, rotation=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 249,
   "metadata": {},
   "outputs": [],
   "source": [
    "# def create_lineplot(data, var, title, legend_title, ax1):\n",
    "#     a = data.groupby(['date', var]).size().reset_index()\n",
    "#     a.columns = ['date', var, 'Freq']\n",
    "#     sns.lineplot(x='date', y='Freq', \n",
    "#             hue=var, \n",
    "#             data=a,\n",
    "#             ax=ax1)\n",
    "#     ax1.set_title(title)\n",
    "#     ax1.set_ylabel('Frequency')\n",
    "#     ax1.xaxis.set_tick_params(which='both', labelbottom=True, labeltop=False, rotation=90)\n",
    "#     ax1.legend(title=legend_title)\n",
    "# #     ax1.axvline(x=\"2020-01\", color='black', ls='--', lw=1.5)\n",
    "# #     ax1.axvline(x=\"2020-04\", color='black', ls='--', lw=1.5)\n",
    "# #     ax1.axvline(x=\"2021-04\", color='black', ls='--', lw=1.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 250,
   "metadata": {},
   "outputs": [],
   "source": [
    "# def create_lineplot_hue(var, hue, title, legend_title, ax1):\n",
    "#     a = data_inclusions_exlusions.groupby(['date', var, hue]).size().reset_index()\n",
    "#     a.columns = ['date', var, hue, 'Freq']\n",
    "#     a = a[a[var] != 0]\n",
    "#     sns.lineplot(x='date', y='Freq', \n",
    "#             hue=hue, \n",
    "#             data=a,\n",
    "#             ax=ax1)\n",
    "#     ax1.set_title(title)\n",
    "#     ax1.set_ylabel('Frequency')\n",
    "#     ax1.xaxis.set_tick_params(which='both', labelbottom=True, labeltop=False, rotation=90)\n",
    "#     ax1.legend(title=legend_title)\n",
    "#     ax1.axvline(x=\"2020-01\", color='black', ls='--', lw=1.5)\n",
    "#     ax1.axvline(x=\"2021-03\", color='black', ls='--', lw=1.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 251,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to transform cph model summary and save it as html\n",
    "def GetPrintSummary(model):\n",
    "    output = \"\"\n",
    "    with io.StringIO() as buf, redirect_stdout(buf):\n",
    "        model.print_summary(style=\"html\")\n",
    "        output = buf.getvalue()\n",
    "    return output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1- Data Preparation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 252,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data2019_month = {}\n",
    "# data2020_month = {}\n",
    "# data2021_month = {}\n",
    "# data2019 = []\n",
    "# data2020 = []\n",
    "# data2021 = []\n",
    "# for  i in range(1, 13):\n",
    "#     # data of 2019\n",
    "#     input_file_path = makeFilePath(\"output/hospitalisation_data/\", \"input_hospitalisation_2019-\"+str(i).zfill(2)+\"-01.csv.gz\")\n",
    "#     data2019_month[\"%s\" %i] = pd.read_csv(input_file_path, compression='gzip')\n",
    "#     data2019_month[\"%s\" %i]['date'] = \"2019-\"+str(i).zfill(2)\n",
    "#     if (i<3) or (i>11):\n",
    "#         data2019_month[\"%s\" %i]['season'] = \"Winter\"\n",
    "#     elif (i>2) and (i<6):\n",
    "#         data2019_month[\"%s\" %i]['season'] = \"Spring\"\n",
    "#     elif (i>5) and (i<9):\n",
    "#         data2019_month[\"%s\" %i]['season'] = \"Summer\"\n",
    "#     elif (i>8) and (i<12):\n",
    "#         data2019_month[\"%s\" %i]['season'] = \"Autumn\"\n",
    "#     data2019.append(data2019_month[\"%s\" %i])\n",
    "#     # data of 2020\n",
    "#     input_file_path = makeFilePath(\"output/hospitalisation_data/\", \"input_hospitalisation_2020-\"+str(i).zfill(2)+\"-01.csv.gz\")\n",
    "#     data2020_month[\"%s\" %i] = pd.read_csv(input_file_path, compression='gzip')\n",
    "#     data2020_month[\"%s\" %i]['date'] = \"2020-\"+str(i).zfill(2)\n",
    "#     if (i<3) or (i>11):\n",
    "#         data2020_month[\"%s\" %i]['season'] = \"Winter\"\n",
    "#     elif (i>2) and (i<6):\n",
    "#         data2020_month[\"%s\" %i]['season'] = \"Spring\"\n",
    "#     elif (i>5) and (i<9):\n",
    "#         data2020_month[\"%s\" %i]['season'] = \"Summer\"\n",
    "#     elif (i>8) and (i<12):\n",
    "#         data2020_month[\"%s\" %i]['season'] = \"Autumn\"\n",
    "#     data2020.append(data2020_month[\"%s\" %i])\n",
    "#     # data of 2021\n",
    "#     input_file_path = makeFilePath(\"output/hospitalisation_data/\", \"input_hospitalisation_2021-\"+str(i).zfill(2)+\"-01.csv.gz\")\n",
    "#     data2021_month[\"%s\" %i] = pd.read_csv(input_file_path, compression='gzip')\n",
    "#     data2021_month[\"%s\" %i]['date'] = \"2021-\"+str(i).zfill(2)\n",
    "#     if (i<3) or (i>11):\n",
    "#         data2021_month[\"%s\" %i]['season'] = \"Winter\"\n",
    "#     elif (i>2) and (i<6):\n",
    "#         data2021_month[\"%s\" %i]['season'] = \"Spring\"\n",
    "#     elif (i>5) and (i<9):\n",
    "#         data2021_month[\"%s\" %i]['season'] = \"Summer\"\n",
    "#     elif (i>8) and (i<12):\n",
    "#         data2021_month[\"%s\" %i]['season'] = \"Autumn\"\n",
    "#     data2021.append(data2021_month[\"%s\" %i])\n",
    "    \n",
    "# data2019 = pd.concat(data2019_month, axis=0, ignore_index=True)\n",
    "# data2020 = pd.concat(data2020_month, axis=0, ignore_index=True)\n",
    "# data2021 = pd.concat(data2021_month, axis=0, ignore_index=True)\n",
    "\n",
    "# data = data2019.append([data2020, data2021], ignore_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 253,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data2019_month = {}\n",
    "# # data2020_month = {}\n",
    "# # data2021_month = {}\n",
    "# # data2022_month = {}\n",
    "# data2019 = []\n",
    "# # data2020 = []\n",
    "# # data2021 = []\n",
    "# # data2022 = []\n",
    "# i=1\n",
    "# while  i < 13:\n",
    "#     # data of 2019\n",
    "#     input_file_path = makeFilePath(\"output/hospitalisation_data/\", \"input_hospitalisation_2019-\"+str(i).zfill(2)+\"-01.csv.gz\")\n",
    "#     data2019_month[\"%s\" %i] = pd.read_csv(input_file_path, compression='gzip')\n",
    "#     data2019_month[\"%s\" %i]['date'] = \"2019-\"+str(i).zfill(2)\n",
    "#     if (i<3) or (i>11):\n",
    "#         data2019_month[\"%s\" %i]['season'] = \"Winter\"\n",
    "#     elif (i>2) and (i<6):\n",
    "#         data2019_month[\"%s\" %i]['season'] = \"Spring\"\n",
    "#     elif (i>5) and (i<9):\n",
    "#         data2019_month[\"%s\" %i]['season'] = \"Summer\"\n",
    "#     elif (i>8) and (i<12):\n",
    "#         data2019_month[\"%s\" %i]['season'] = \"Autumn\"\n",
    "#     data2019.append(data2019_month[\"%s\" %i])\n",
    "#     i=i+1\n",
    "    \n",
    "# # i=1\n",
    "# # while  i < 13:    \n",
    "# #     # data of 2020\n",
    "# #     input_file_path = makeFilePath(\"output/hospitalisation_data/\", \"input_hospitalisation_2020-\"+str(i).zfill(2)+\"-01.csv.gz\")\n",
    "# #     data2020_month[\"%s\" %i] = pd.read_csv(input_file_path, compression='gzip')\n",
    "# #     data2020_month[\"%s\" %i]['date'] = \"2020-\"+str(i).zfill(2)\n",
    "# #     if (i<3) or (i>11):\n",
    "# #         data2020_month[\"%s\" %i]['season'] = \"Winter\"\n",
    "# #     elif (i>2) and (i<6):\n",
    "# #         data2020_month[\"%s\" %i]['season'] = \"Spring\"\n",
    "# #     elif (i>5) and (i<9):\n",
    "# #         data2020_month[\"%s\" %i]['season'] = \"Summer\"\n",
    "# #     elif (i>8) and (i<12):\n",
    "# #         data2020_month[\"%s\" %i]['season'] = \"Autumn\"\n",
    "# #     data2020.append(data2020_month[\"%s\" %i])\n",
    "# #     i=i+1\n",
    "\n",
    "# # i=1\n",
    "# # while  i < 13:\n",
    "# #     # data of 2021\n",
    "# #     input_file_path = makeFilePath(\"output/hospitalisation_data/\", \"input_hospitalisation_2021-\"+str(i).zfill(2)+\"-01.csv.gz\")\n",
    "# #     data2021_month[\"%s\" %i] = pd.read_csv(input_file_path, compression='gzip')\n",
    "# #     data2021_month[\"%s\" %i]['date'] = \"2021-\"+str(i).zfill(2)\n",
    "# #     if (i<3) or (i>11):\n",
    "# #         data2021_month[\"%s\" %i]['season'] = \"Winter\"\n",
    "# #     elif (i>2) and (i<6):\n",
    "# #         data2021_month[\"%s\" %i]['season'] = \"Spring\"\n",
    "# #     elif (i>5) and (i<9):\n",
    "# #         data2021_month[\"%s\" %i]['season'] = \"Summer\"\n",
    "# #     elif (i>8) and (i<12):\n",
    "# #         data2021_month[\"%s\" %i]['season'] = \"Autumn\"\n",
    "# #     data2021.append(data2021_month[\"%s\" %i])\n",
    "# #     i=i+1\n",
    "\n",
    "# # i=1\n",
    "# # while  i < 4:\n",
    "# #     # data of 2022\n",
    "# #     input_file_path = makeFilePath(\"output/hospitalisation_data/\", \"input_hospitalisation_2022-\"+str(i).zfill(2)+\"-01.csv.gz\")\n",
    "# #     data2022_month[\"%s\" %i] = pd.read_csv(input_file_path, compression='gzip')\n",
    "# #     data2022_month[\"%s\" %i]['date'] = \"2022-\"+str(i).zfill(2)\n",
    "# #     if (i<3) or (i>11):\n",
    "# #         data2022_month[\"%s\" %i]['season'] = \"Winter\"\n",
    "# #     elif (i>2) and (i<6):\n",
    "# #         data2022_month[\"%s\" %i]['season'] = \"Spring\"\n",
    "# #     elif (i>5) and (i<9):\n",
    "# #         data2022_month[\"%s\" %i]['season'] = \"Summer\"\n",
    "# #     elif (i>8) and (i<12):\n",
    "# #         data2022_month[\"%s\" %i]['season'] = \"Autumn\"\n",
    "# #     data2022.append(data2022_month[\"%s\" %i])\n",
    "# #     i=i+1\n",
    "    \n",
    "# data2019 = pd.concat(data2019_month, axis=0, ignore_index=True)\n",
    "# # data2020 = pd.concat(data2020_month, axis=0, ignore_index=True)\n",
    "# # data2021 = pd.concat(data2021_month, axis=0, ignore_index=True)\n",
    "# # data2022 = pd.concat(data2022_month, axis=0, ignore_index=True)\n",
    "\n",
    "# # data = data2019.append([data2020, data2021, data2022], ignore_index=True)\n",
    "\n",
    "# data = data2019"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 254,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data2019_month = {}\n",
    "# data2019 = []\n",
    "\n",
    "# i=1\n",
    "# while  i < 13:\n",
    "#     # data of 2019\n",
    "#     # if (i<4):\n",
    "#     data2019_month[\"%s\" %i] = pd.read_csv(f'../output/hospitalisation_data/input_hospitalisation_2019-'+str(i).zfill(2)+'-01.csv.gz')\n",
    "#     data2019_month[\"%s\" %i]['date'] = \"2019-\"+str(i).zfill(2)\n",
    "#     # if (i>3) and (i<7):\n",
    "#     #     data2019_month[\"%s\" %i] = pd.read_csv(f'output/hospitalisation_data/input_hospitalisatio2019-'+str(i).zfill(2)+'-01.csv.gz')\n",
    "#     #     data2019_month[\"%s\" %i]['date'] = \"2019-\"+str(i).zfill(2)\n",
    "#     if (i<3) or (i>11):\n",
    "#         data2019_month[\"%s\" %i]['season'] = \"Winter\"\n",
    "#     elif (i>2) and (i<6):\n",
    "#         data2019_month[\"%s\" %i]['season'] = \"Spring\"\n",
    "#     elif (i>5) and (i<9):\n",
    "#         data2019_month[\"%s\" %i]['season'] = \"Summer\"\n",
    "#     elif (i>8) and (i<12):\n",
    "#         data2019_month[\"%s\" %i]['season'] = \"Autumn\"\n",
    "\n",
    "#     data2019.append(data2019_month[\"%s\" %i])\n",
    "#     i=i+1\n",
    "    \n",
    "# data2019 = pd.concat(data2019_month, axis=0, ignore_index=True)\n",
    "\n",
    "# data = data2019"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 255,
   "metadata": {},
   "outputs": [],
   "source": [
    "data201901 = pd.read_csv(f'../output/hospitalisation_data/input_hospitalisation_2019-01-01.csv.gz')\n",
    "# data201901['date'] = \"2019-01\"\n",
    "# data201901['season'] = \"Winter\"\n",
    "data201902 = pd.read_csv(f'../output/hospitalisation_data/input_hospitalisation_2019-02-01.csv.gz')\n",
    "# data201902['date'] = \"2019-02\"\n",
    "# data201902['season'] = \"Winter\"\n",
    "data201903 = pd.read_csv(f'../output/hospitalisation_data/input_hospitalisation_2019-03-01.csv.gz')\n",
    "# data201903['date'] = \"2019-03\"\n",
    "# data201903['season'] = \"Spring\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 256,
   "metadata": {},
   "outputs": [],
   "source": [
    "data2019 = [data201901, data201902, data201903]\n",
    "data = pd.concat(data2019)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 257,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data = data201901"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 258,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data2019 = [data201901, data201902]\n",
    "# data = pd.concat(data2019)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 259,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data = data.loc[data['sex'] == 'F']\n",
    "# data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 260,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>bmi_date_measured</th>\n",
       "      <th>smoking_status_date</th>\n",
       "      <th>most_recent_unclear_smoking_cat_date</th>\n",
       "      <th>flu_vaccine_med</th>\n",
       "      <th>flu_vaccine_clinical</th>\n",
       "      <th>sgss_positive</th>\n",
       "      <th>primary_care_covid</th>\n",
       "      <th>admitted_date</th>\n",
       "      <th>deregistered_date</th>\n",
       "      <th>died_date</th>\n",
       "      <th>...</th>\n",
       "      <th>hemiplegia_comor</th>\n",
       "      <th>hiv_comor</th>\n",
       "      <th>metastatic_cancer_comor</th>\n",
       "      <th>mild_liver_comor</th>\n",
       "      <th>mod_severe_liver_comor</th>\n",
       "      <th>mod_severe_renal_comor</th>\n",
       "      <th>mi_comor</th>\n",
       "      <th>peptic_ulcer_comor</th>\n",
       "      <th>peripheral_vascular_comor</th>\n",
       "      <th>patient_id</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>22500</td>\n",
       "      <td>3000</td>\n",
       "      <td>3000</td>\n",
       "      <td>3000</td>\n",
       "      <td>3000</td>\n",
       "      <td>3000</td>\n",
       "      <td>3000</td>\n",
       "      <td>9000</td>\n",
       "      <td>1500</td>\n",
       "      <td>3000</td>\n",
       "      <td>...</td>\n",
       "      <td>30000.000000</td>\n",
       "      <td>30000.000000</td>\n",
       "      <td>30000.000000</td>\n",
       "      <td>30000.000000</td>\n",
       "      <td>30000.000000</td>\n",
       "      <td>30000.000000</td>\n",
       "      <td>30000.000000</td>\n",
       "      <td>30000.000000</td>\n",
       "      <td>30000.000000</td>\n",
       "      <td>30000.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>unique</th>\n",
       "      <td>40</td>\n",
       "      <td>40</td>\n",
       "      <td>40</td>\n",
       "      <td>13</td>\n",
       "      <td>13</td>\n",
       "      <td>1104</td>\n",
       "      <td>1102</td>\n",
       "      <td>1194</td>\n",
       "      <td>830</td>\n",
       "      <td>461</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>top</th>\n",
       "      <td>2021-12</td>\n",
       "      <td>2021-12</td>\n",
       "      <td>2020-10</td>\n",
       "      <td>2018-03</td>\n",
       "      <td>2018-05</td>\n",
       "      <td>2020-10-27</td>\n",
       "      <td>2020-07-22</td>\n",
       "      <td>2021-06-01</td>\n",
       "      <td>2020-03-15</td>\n",
       "      <td>2022-04-10</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>freq</th>\n",
       "      <td>627</td>\n",
       "      <td>105</td>\n",
       "      <td>90</td>\n",
       "      <td>278</td>\n",
       "      <td>276</td>\n",
       "      <td>8</td>\n",
       "      <td>9</td>\n",
       "      <td>20</td>\n",
       "      <td>7</td>\n",
       "      <td>35</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>0.100000</td>\n",
       "      <td>0.100000</td>\n",
       "      <td>0.100000</td>\n",
       "      <td>0.100000</td>\n",
       "      <td>0.100000</td>\n",
       "      <td>0.100000</td>\n",
       "      <td>0.100000</td>\n",
       "      <td>0.100000</td>\n",
       "      <td>0.100000</td>\n",
       "      <td>50120.305433</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>0.300005</td>\n",
       "      <td>0.300005</td>\n",
       "      <td>0.300005</td>\n",
       "      <td>0.300005</td>\n",
       "      <td>0.300005</td>\n",
       "      <td>0.300005</td>\n",
       "      <td>0.300005</td>\n",
       "      <td>0.300005</td>\n",
       "      <td>0.300005</td>\n",
       "      <td>28941.362264</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>25133.750000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>50142.500000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>75270.250000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>99994.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>11 rows × 85 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       bmi_date_measured smoking_status_date  \\\n",
       "count              22500                3000   \n",
       "unique                40                  40   \n",
       "top              2021-12             2021-12   \n",
       "freq                 627                 105   \n",
       "mean                 NaN                 NaN   \n",
       "std                  NaN                 NaN   \n",
       "min                  NaN                 NaN   \n",
       "25%                  NaN                 NaN   \n",
       "50%                  NaN                 NaN   \n",
       "75%                  NaN                 NaN   \n",
       "max                  NaN                 NaN   \n",
       "\n",
       "       most_recent_unclear_smoking_cat_date flu_vaccine_med  \\\n",
       "count                                  3000            3000   \n",
       "unique                                   40              13   \n",
       "top                                 2020-10         2018-03   \n",
       "freq                                     90             278   \n",
       "mean                                    NaN             NaN   \n",
       "std                                     NaN             NaN   \n",
       "min                                     NaN             NaN   \n",
       "25%                                     NaN             NaN   \n",
       "50%                                     NaN             NaN   \n",
       "75%                                     NaN             NaN   \n",
       "max                                     NaN             NaN   \n",
       "\n",
       "       flu_vaccine_clinical sgss_positive primary_care_covid admitted_date  \\\n",
       "count                  3000          3000               3000          9000   \n",
       "unique                   13          1104               1102          1194   \n",
       "top                 2018-05    2020-10-27         2020-07-22    2021-06-01   \n",
       "freq                    276             8                  9            20   \n",
       "mean                    NaN           NaN                NaN           NaN   \n",
       "std                     NaN           NaN                NaN           NaN   \n",
       "min                     NaN           NaN                NaN           NaN   \n",
       "25%                     NaN           NaN                NaN           NaN   \n",
       "50%                     NaN           NaN                NaN           NaN   \n",
       "75%                     NaN           NaN                NaN           NaN   \n",
       "max                     NaN           NaN                NaN           NaN   \n",
       "\n",
       "       deregistered_date   died_date  ... hemiplegia_comor     hiv_comor  \\\n",
       "count               1500        3000  ...     30000.000000  30000.000000   \n",
       "unique               830         461  ...              NaN           NaN   \n",
       "top           2020-03-15  2022-04-10  ...              NaN           NaN   \n",
       "freq                   7          35  ...              NaN           NaN   \n",
       "mean                 NaN         NaN  ...         0.100000      0.100000   \n",
       "std                  NaN         NaN  ...         0.300005      0.300005   \n",
       "min                  NaN         NaN  ...         0.000000      0.000000   \n",
       "25%                  NaN         NaN  ...         0.000000      0.000000   \n",
       "50%                  NaN         NaN  ...         0.000000      0.000000   \n",
       "75%                  NaN         NaN  ...         0.000000      0.000000   \n",
       "max                  NaN         NaN  ...         1.000000      1.000000   \n",
       "\n",
       "       metastatic_cancer_comor mild_liver_comor mod_severe_liver_comor  \\\n",
       "count             30000.000000     30000.000000           30000.000000   \n",
       "unique                     NaN              NaN                    NaN   \n",
       "top                        NaN              NaN                    NaN   \n",
       "freq                       NaN              NaN                    NaN   \n",
       "mean                  0.100000         0.100000               0.100000   \n",
       "std                   0.300005         0.300005               0.300005   \n",
       "min                   0.000000         0.000000               0.000000   \n",
       "25%                   0.000000         0.000000               0.000000   \n",
       "50%                   0.000000         0.000000               0.000000   \n",
       "75%                   0.000000         0.000000               0.000000   \n",
       "max                   1.000000         1.000000               1.000000   \n",
       "\n",
       "       mod_severe_renal_comor      mi_comor peptic_ulcer_comor  \\\n",
       "count            30000.000000  30000.000000       30000.000000   \n",
       "unique                    NaN           NaN                NaN   \n",
       "top                       NaN           NaN                NaN   \n",
       "freq                      NaN           NaN                NaN   \n",
       "mean                 0.100000      0.100000           0.100000   \n",
       "std                  0.300005      0.300005           0.300005   \n",
       "min                  0.000000      0.000000           0.000000   \n",
       "25%                  0.000000      0.000000           0.000000   \n",
       "50%                  0.000000      0.000000           0.000000   \n",
       "75%                  0.000000      0.000000           0.000000   \n",
       "max                  1.000000      1.000000           1.000000   \n",
       "\n",
       "       peripheral_vascular_comor    patient_id  \n",
       "count               30000.000000  30000.000000  \n",
       "unique                       NaN           NaN  \n",
       "top                          NaN           NaN  \n",
       "freq                         NaN           NaN  \n",
       "mean                    0.100000  50120.305433  \n",
       "std                     0.300005  28941.362264  \n",
       "min                     0.000000      0.000000  \n",
       "25%                     0.000000  25133.750000  \n",
       "50%                     0.000000  50142.500000  \n",
       "75%                     0.000000  75270.250000  \n",
       "max                     1.000000  99994.000000  \n",
       "\n",
       "[11 rows x 85 columns]"
      ]
     },
     "execution_count": 260,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.describe(include='all')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1-3-2- Imputation of missing values\n",
    "We impute variables with missing values that at most 50 pecent of their values are missing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 261,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prep before plot\n",
    "data['age'][data['age'] < 0] = np.nan\n",
    "\n",
    "# ## remove died cases\n",
    "# #data.died_date.value_counts()\n",
    "\n",
    "# # remove cases with bmi equal 0\n",
    "# data['bmi'][data['bmi'] <= 0] = 0\n",
    "\n",
    "# # remove cases with negative gp_count\n",
    "# data['gp_count'][data['gp_count'] < 0] = np.nan\n",
    "\n",
    "# # remove cases with negative antibacterial_brit\n",
    "# data['antibacterial_brit'][data['antibacterial_brit'] < 0] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 262,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['Age_cat'] = pd.cut(data['age'], right=False, bins=[0,18,60,200], labels=['0-17 years', '18-59 years', '>59 years'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3- Data preparation\n",
    "This includes data type conversion and exclusion of cases that are not considered in the infection-related hospital admission, i.e. those who were diagnosed with Covid19 three months (=90 days) before and one month (=30 days) after hospital admission.\n",
    "\n",
    "### 3-1- URTI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 263,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ConvertColumns(InfectionKey, NumSets):\n",
    "    for i in range(1, NumSets+1):        \n",
    "        arrayKey1 = 'admitted_'+InfectionKey+'_date_'+str(i)\n",
    "        #arrayKey2 = InfectionKey+'_date_'+str(i)\n",
    "        data[arrayKey1] = data[arrayKey1].astype('datetime64[ns]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### URTI-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 264,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(81, 9)"
      ]
     },
     "execution_count": 264,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_urti_1 = data[['sex', #'msoa', #'date', \n",
    "                #     'age', 'age_cat', 'sex', 'bmi_imp', 'BMI_cat', 'gp_count', 'gp_count_cat', 'gp_count_imp', 'flu_vaccine', \n",
    "                #    'region', 'smoking_status_imp', 'antibacterial_brit', 'imd_imp', 'ethnicity_imp', 'antibacterial_brit_imp',\n",
    "                #    'hx_hosp', 'CCI_cat', 'CCI', 'meds_nonabs', #'season', \n",
    "                   'died_date', 'deregistered_date', \n",
    "                    'gp_cons_urti_1', 'gp_cons_urti_ab_1',\n",
    "                    'urti_date_1', 'incdt_urti_date_1', 'admitted_urti_date_1',\n",
    "                    'sgss_gp_cov_urti_date_1']]\n",
    "\n",
    "data_urti_1 = data_urti_1[data_urti_1['urti_date_1'].notnull()]\n",
    "\n",
    "data_urti_1.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 265,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(67, 9)"
      ]
     },
     "execution_count": 265,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# exclusion of covid positive 90 days before and 30 days after dx with urti_1\n",
    "data_urti_1 = data_urti_1[data_urti_1['sgss_gp_cov_urti_date_1'] == 0]\n",
    "data_urti_1.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 266,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['age', 'died_date', 'deregistered_date', 'gp_cons_urti_1',\n",
       "       'gp_cons_urti_ab_1', 'urti_date_1', 'incdt_urti_date_1',\n",
       "       'admitted_urti_date_1', 'sgss_gp_cov_urti_date_1'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 266,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_urti_1.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### URTI-2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 267,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(81, 9)"
      ]
     },
     "execution_count": 267,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_urti_2 = data[['sex', #'msoa', #'date',\n",
    "                #     'age', 'age_cat', 'sex', 'bmi_imp', 'BMI_cat', 'gp_count', 'gp_count_cat', 'gp_count_imp', 'flu_vaccine', \n",
    "                #    'region', 'smoking_status_imp', 'antibacterial_brit', 'imd_imp', 'ethnicity_imp', 'antibacterial_brit_imp',\n",
    "                #    'hx_hosp', 'CCI_cat', 'CCI', 'meds_nonabs',  #'season',\n",
    "                   'died_date', 'deregistered_date', \n",
    "                    'gp_cons_urti_2', 'gp_cons_urti_ab_2',\n",
    "                    'urti_date_2', 'incdt_urti_date_2', 'admitted_urti_date_2',\n",
    "                    'sgss_gp_cov_urti_date_2']]\n",
    "\n",
    "data_urti_2 = data_urti_2[data_urti_2['urti_date_2'].notnull()]\n",
    "\n",
    "data_urti_2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 268,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(75, 9)"
      ]
     },
     "execution_count": 268,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# exclusion of covid positive 90 days before and 30 days after dx with urti_2\n",
    "data_urti_2 = data_urti_2[data_urti_2['sgss_gp_cov_urti_date_2'] == 0]\n",
    "data_urti_2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 269,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['age', 'died_date', 'deregistered_date', 'gp_cons_urti_2',\n",
       "       'gp_cons_urti_ab_2', 'urti_date_2', 'incdt_urti_date_2',\n",
       "       'admitted_urti_date_2', 'sgss_gp_cov_urti_date_2'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 269,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_urti_2.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### URTI_3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 270,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(79, 9)"
      ]
     },
     "execution_count": 270,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_urti_3 = data[['sex', #'msoa', #'date',\n",
    "                #     'age', 'age_cat', 'sex', 'bmi_imp', 'BMI_cat', 'gp_count', 'gp_count_cat', 'gp_count_imp', 'flu_vaccine', \n",
    "                #    'region', 'smoking_status_imp', 'antibacterial_brit', 'imd_imp', 'ethnicity_imp', 'antibacterial_brit_imp',\n",
    "                #    'hx_hosp', 'CCI_cat', 'CCI', 'meds_nonabs',  #'season',\n",
    "                   'died_date', 'deregistered_date', \n",
    "                    'gp_cons_urti_3', 'gp_cons_urti_ab_3',\n",
    "                    'urti_date_3', 'incdt_urti_date_3', 'admitted_urti_date_3',\n",
    "                    'sgss_gp_cov_urti_date_3']]\n",
    "\n",
    "data_urti_3 = data_urti_3[data_urti_3['urti_date_3'].notnull()]\n",
    "\n",
    "data_urti_3.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 271,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(74, 9)"
      ]
     },
     "execution_count": 271,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# exclusion of covid positive 90 days before and 30 days after dx with urti_3\n",
    "data_urti_3 = data_urti_3[data_urti_3['sgss_gp_cov_urti_date_3'] == 0]\n",
    "data_urti_3.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 272,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['age', 'died_date', 'deregistered_date', 'gp_cons_urti_3',\n",
       "       'gp_cons_urti_ab_3', 'urti_date_3', 'incdt_urti_date_3',\n",
       "       'admitted_urti_date_3', 'sgss_gp_cov_urti_date_3'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 272,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_urti_3.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### URTI_4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 273,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(82, 9)"
      ]
     },
     "execution_count": 273,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_urti_4 = data[['sex', #'msoa', #'date',\n",
    "                #     'age', 'age_cat', 'sex', 'bmi_imp', 'BMI_cat', 'gp_count', 'gp_count_cat', 'gp_count_imp', 'flu_vaccine', \n",
    "                #    'region', 'smoking_status_imp', 'antibacterial_brit', 'imd_imp', 'ethnicity_imp', 'antibacterial_brit_imp',\n",
    "                #    'hx_hosp', 'CCI_cat', 'CCI', 'meds_nonabs',  #'season',\n",
    "                   'died_date', 'deregistered_date', \n",
    "                   'gp_cons_urti_4', 'gp_cons_urti_ab_4',\n",
    "                   'urti_date_4', 'incdt_urti_date_4', 'admitted_urti_date_4',\n",
    "                   'sgss_gp_cov_urti_date_4']]\n",
    "\n",
    "data_urti_4 = data_urti_4[data_urti_4['urti_date_4'].notnull()]\n",
    "\n",
    "data_urti_4.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 274,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(73, 9)"
      ]
     },
     "execution_count": 274,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# exclusion of covid positive 90 days before and 30 days after dx with urti_4\n",
    "data_urti_4 = data_urti_4[data_urti_4['sgss_gp_cov_urti_date_4'] == 0]\n",
    "data_urti_4.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 275,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['age', 'died_date', 'deregistered_date', 'gp_cons_urti_4',\n",
       "       'gp_cons_urti_ab_4', 'urti_date_4', 'incdt_urti_date_4',\n",
       "       'admitted_urti_date_4', 'sgss_gp_cov_urti_date_4'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 275,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_urti_4.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Combining URTI_1, URTI_2, URTI_3, and URTI_4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 276,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(289, 9)"
      ]
     },
     "execution_count": 276,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_urti_1.rename(columns={'urti_date_1': 'urti_date', 'incdt_urti_date_1': 'incdt_urti_date', \n",
    "                           'gp_cons_urti_1': 'gp_cons_urti', 'gp_cons_urti_ab_1': 'gp_cons_urti_ab',\n",
    "                           'incdt_urti_date_1': 'incdt_urti_date', #'prevl_urti_date_1': 'prevl_urti_date',\n",
    "                           'admitted_urti_date_1': 'admitted_urti_date',\n",
    "                           'sgss_gp_cov_urti_date_1': 'sgss_gp_cov_urti_date'},\n",
    "                inplace=True)\n",
    "data_urti_2.rename(columns={'urti_date_2': 'urti_date', 'incdt_urti_date_2': 'incdt_urti_date',\n",
    "                           'gp_cons_urti_2': 'gp_cons_urti', 'gp_cons_urti_ab_2': 'gp_cons_urti_ab',\n",
    "                           'incdt_urti_date_2': 'incdt_urti_date', #'prevl_urti_date_2': 'prevl_urti_date', \n",
    "                           'admitted_urti_date_2': 'admitted_urti_date',\n",
    "                           'sgss_gp_cov_urti_date_2': 'sgss_gp_cov_urti_date'},\n",
    "                inplace=True)\n",
    "data_urti_3.rename(columns={'urti_date_3': 'urti_date', 'incdt_urti_date_3': 'incdt_urti_date',\n",
    "                           'gp_cons_urti_3': 'gp_cons_urti', 'gp_cons_urti_ab_3': 'gp_cons_urti_ab',\n",
    "                           'incdt_urti_date_3': 'incdt_urti_date', #'prevl_urti_date_3': 'prevl_urti_date', \n",
    "                           'admitted_urti_date_3': 'admitted_urti_date',\n",
    "                           'sgss_gp_cov_urti_date_3': 'sgss_gp_cov_urti_date'},\n",
    "                 inplace=True)\n",
    "data_urti_4.rename(columns={'urti_date_4': 'urti_date', 'incdt_urti_date_4': 'incdt_urti_date', \n",
    "                           'gp_cons_urti_4': 'gp_cons_urti', 'gp_cons_urti_ab_4': 'gp_cons_urti_ab',\n",
    "                           'incdt_urti_date_4': 'incdt_urti_date', #'prevl_urti_date_4': 'prevl_urti_date',\n",
    "                           'admitted_urti_date_4': 'admitted_urti_date',\n",
    "                           'sgss_gp_cov_urti_date_4': 'sgss_gp_cov_urti_date'}, \n",
    "                 inplace=True)\n",
    "\n",
    "data_urti = pd.concat([data_urti_1, data_urti_2, data_urti_3, data_urti_4])\n",
    "data_urti.reset_index(inplace=True, drop=True)\n",
    "data_urti.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 277,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig, axes = plt.subplots(7,2,\n",
    "#                         figsize=(15, 33)\n",
    "#                         )\n",
    "# create_lineplot(data_urti, 'sex', 'Sex', 'Sex', axes[0,0])\n",
    "# create_lineplot(data_urti, 'Age_cat', 'Age category', 'Age category', axes[0,1])\n",
    "# create_lineplot(data_urti, 'imd_imp', 'IMD', 'IMD', axes[1,0])\n",
    "# create_lineplot(data_urti, 'ethnicity_imp', 'Ethnicity', 'Ethnicity', axes[1,1])\n",
    "# create_lineplot(data_urti, 'smoking_status_imp', 'Smoking status', 'Smoking status', axes[2,0])\n",
    "# create_lineplot(data_urti, 'BMI_cat', 'BMI', 'BMI', axes[2,1])\n",
    "# create_lineplot(data_urti, 'CCI_cat', 'CCI', 'CCI', axes[3,0])\n",
    "# create_lineplot(data_urti, 'gp_count_cat', 'GP consultation', 'GP consultation', axes[3,1])\n",
    "# create_lineplot(data_urti, 'region', 'Region', 'Region', axes[4,0])\n",
    "# create_lineplot(data_urti, 'flu_vaccine', 'Flu vaccine', 'Flu vaccine', axes[4,1])\n",
    "# create_lineplot(data_urti, 'antibacterial_brit', 'Antibacterial prescriptions', 'Antibacterial prescriptions', axes[5,0])\n",
    "# create_lineplot(data_urti, 'msoa', 'msoa', 'MSOA', axes[5,1])\n",
    "# create_lineplot(data_urti, 'meds_nonabs', 'Medicines excpet antibacterials', 'Medicines excpet antibacterials', axes[6,0])\n",
    "# fig.delaxes(ax=axes[6,1])\n",
    "\n",
    "# fig.suptitle('Risk Factors of URTI patients', fontsize=16)\n",
    "# fig.tight_layout(rect=[0, 0.03, 1, 0.98])\n",
    "# #plt.tight_layout()\n",
    "# plt.savefig('../output/hospitalisation_prediction_urti/risk_factors_cat_urti.jpg')\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Adding event and duration columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 278,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(289, 9)"
      ]
     },
     "execution_count": 278,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_urti['urti_date'] = data_urti['urti_date'].astype('datetime64[ns]')\n",
    "data_urti['admitted_urti_date'] = data_urti['admitted_urti_date'].astype('datetime64[ns]')\n",
    "data_urti.died_date = data_urti.died_date.astype('datetime64[ns]')\n",
    "data_urti.deregistered_date = data_urti.deregistered_date.astype('datetime64[ns]')\n",
    "data_urti.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 279,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0    197\n",
       "Name: event_urti_admitted, dtype: int64"
      ]
     },
     "execution_count": 279,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#scenario 1\n",
    "#not hosped (nothing happened)\n",
    "#data_urti = data_urti[data_urti['admitted_urti_date'].notnull()]\n",
    "data_urti.loc[data_urti['admitted_urti_date'].isnull(), 'event_urti_admitted'] = 0\n",
    "data_urti['event_urti_admitted'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 281,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0    197\n",
       "1.0     13\n",
       "Name: event_urti_admitted, dtype: int64"
      ]
     },
     "execution_count": 281,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#scenario 2\n",
    "# become a case (uncensoring)\n",
    "#calculating days between infection dx and hosp\n",
    "data_urti['delta_urti_admitted'] = (data_urti['admitted_urti_date'] - data_urti['urti_date']).dt.days\n",
    "data_urti.loc[((data_urti['delta_urti_admitted'] > 0) & (data_urti['delta_urti_admitted'] < 30)), 'event_urti_admitted'] = 1\n",
    "data_urti['event_urti_admitted'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 282,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "506.0    2\n",
       "495.0    2\n",
       "735.0    2\n",
       "505.0    2\n",
       "524.0    2\n",
       "        ..\n",
       "470.0    1\n",
       "618.0    1\n",
       "31.0     1\n",
       "492.0    1\n",
       "208.0    1\n",
       "Name: delta_urti_admitted, Length: 85, dtype: int64"
      ]
     },
     "execution_count": 282,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#scenario 2\n",
    "#drop whoever was admitted before urti consultation, i.e. negative value for delta_urti_admitted\n",
    "data_urti = data_urti[~(data_urti['delta_urti_admitted'] < 0)]\n",
    "data_urti['delta_urti_admitted'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 283,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0    207\n",
       "1.0     11\n",
       "Name: event_urti_admitted, dtype: int64"
      ]
     },
     "execution_count": 283,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#scenario 3\n",
    "#censor died patients\n",
    "data_urti['delta_admitted_died'] = (data_urti['died_date'] - data_urti['admitted_urti_date']).dt.days\n",
    "data_urti.loc[data_urti['delta_admitted_died'] < 0, 'delta_admitted_died'] = np.NaN\n",
    "data_urti.loc[data_urti['delta_admitted_died'] > 0, 'event_urti_admitted'] = 0\n",
    "data_urti['event_urti_admitted'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 284,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0    207\n",
       "1.0     11\n",
       "Name: event_urti_admitted, dtype: int64"
      ]
     },
     "execution_count": 284,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#scenario 3\n",
    "#censor deregistered patients\n",
    "data_urti['delta_admitted_deregistered'] = (data_urti['deregistered_date'] - data_urti['admitted_urti_date']).dt.days\n",
    "data_urti.loc[data_urti['delta_admitted_deregistered'] < 0, 'delta_admitted_deregistered'] = np.NaN\n",
    "data_urti.loc[((data_urti['delta_admitted_deregistered'] > 0) & (data_urti['delta_admitted_deregistered'] < 30)), 'event_urti_admitted'] = 0\n",
    "data_urti['event_urti_admitted'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 285,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "279"
      ]
     },
     "execution_count": 285,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#agg scenario 3s\n",
    "data_urti['delta_admitted_died_deregistered'] = data_urti['delta_admitted_deregistered'].combine_first(data_urti['delta_admitted_died'])\n",
    "data_urti.loc[data_urti['delta_admitted_died_deregistered'] < 0, 'delta_admitted_died_deregistered'] = np.NaN\n",
    "data_urti['delta_admitted_died_deregistered'].isnull().sum()#.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 286,
   "metadata": {},
   "outputs": [],
   "source": [
    "#censor patients that had a history of hospital admission\n",
    "# data_urti.loc[data_urti['hx_hosp'] == 1, 'event_urti_admitted'] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 287,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0    278\n",
       "1.0     11\n",
       "Name: event_urti_admitted, dtype: int64"
      ]
     },
     "execution_count": 287,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#scenario 1\n",
    "#any other patients (nothing happened)\n",
    "data_urti['event_urti_admitted'] = data_urti['event_urti_admitted'].replace(np.NaN, 0)\n",
    "data_urti['event_urti_admitted'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 288,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "230.0    267\n",
       "226.0      1\n",
       "32.0       1\n",
       "105.0      1\n",
       "194.0      1\n",
       "209.0      1\n",
       "68.0       1\n",
       "65.0       1\n",
       "137.0      1\n",
       "112.0      1\n",
       "38.0       1\n",
       "223.0      1\n",
       "87.0       1\n",
       "132.0      1\n",
       "125.0      1\n",
       "220.0      1\n",
       "106.0      1\n",
       "31.0       1\n",
       "16.0       1\n",
       "118.0      1\n",
       "180.0      1\n",
       "123.0      1\n",
       "208.0      1\n",
       "Name: duration_urti_admitted, dtype: int64"
      ]
     },
     "execution_count": 288,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#assign values for duration column\n",
    "data_urti['duration_urti_admitted'] = data_urti['delta_urti_admitted'].combine_first(data_urti['delta_admitted_died_deregistered'])\n",
    "data_urti['duration_urti_admitted'] = data_urti['duration_urti_admitted'].replace(np.NaN, 30)\n",
    "data_urti.loc[(data_urti['duration_urti_admitted'] > 30), 'duration_urti_admitted'] = 30\n",
    "data_urti['duration_urti_admitted'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>delta_urti_admitted</th>\n",
       "      <th>urti_date</th>\n",
       "      <th>admitted_urti_date</th>\n",
       "      <th>delta_urti_admitted</th>\n",
       "      <th>delta_admitted_died_deregistered</th>\n",
       "      <th>duration_urti_admitted</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [delta_urti_admitted, urti_date, admitted_urti_date, delta_urti_admitted, delta_admitted_died_deregistered, duration_urti_admitted]\n",
       "Index: []"
      ]
     },
     "execution_count": 164,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_urti[data_urti['duration_urti_admitted']<0][['delta_urti_admitted', 'urti_date', 'admitted_urti_date', 'delta_urti_admitted', 'delta_admitted_died_deregistered', 'duration_urti_admitted']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {},
   "outputs": [],
   "source": [
    "#scenario2 (uncensoring) again to prevent conflict with other scenarios\n",
    "data_urti.loc[((data_urti['delta_urti_admitted'] > 0) & (data_urti['delta_urti_admitted'] < 30)), 'event_urti_admitted'] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Separate data of URTI\n",
    "We separate data of urti to form four sub-datasets for 1) no antibiotics and incident admission, 2) with antibiotics and incident admission, 3) no antibiotics and prevalent admission, and 4) with antiobiotics and prevalent admission."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shape of data_urti_no_abs_incdt: (239, 15)\n"
     ]
    }
   ],
   "source": [
    "# no antibiotics and incident hospital admission\n",
    "data_urti_no_abs_incdt = data_urti[data_urti['gp_cons_urti_ab'].isnull()]\n",
    "data_urti_no_abs_incdt = data_urti_no_abs_incdt[data_urti_no_abs_incdt['incdt_urti_date'] == 0]\n",
    "\n",
    "# # with antibiotics and incident hospital admission\n",
    "# data_urti_abs_incdt = data_urti[data_urti['gp_cons_urti_ab'].notnull()]\n",
    "# data_urti_abs_incdt = data_urti_abs_incdt[data_urti_abs_incdt['incdt_urti_date'] == 0]\n",
    "\n",
    "# # no antibiotics and prevalent hospital admission\n",
    "# data_urti_no_abs_prevl = data_urti[data_urti['gp_cons_urti_ab'].isnull()]\n",
    "# data_urti_no_abs_prevl = data_urti_no_abs_prevl[data_urti_no_abs_prevl['incdt_urti_date'] == 1]\n",
    "\n",
    "# # with antibiotics and prevalent hospital admission\n",
    "# data_urti_abs_prevl = data_urti[data_urti['gp_cons_urti_ab'].notnull()]\n",
    "# data_urti_abs_prevl = data_urti_abs_prevl[data_urti_abs_prevl['incdt_urti_date'] == 1]\n",
    "\n",
    "print(f'shape of data_urti_no_abs_incdt: {data_urti_no_abs_incdt.shape}')\n",
    "# print(f'shape of data_urti_abs_incdt: {data_urti_abs_incdt.shape}')\n",
    "# print(f'shape of data_urti_no_abs_prevl: {data_urti_no_abs_prevl.shape}')\n",
    "# print(f'shape of data_urti_abs_prevl: {data_urti_abs_prevl.shape}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>age</th>\n",
       "      <th>died_date</th>\n",
       "      <th>deregistered_date</th>\n",
       "      <th>gp_cons_urti</th>\n",
       "      <th>gp_cons_urti_ab</th>\n",
       "      <th>urti_date</th>\n",
       "      <th>incdt_urti_date</th>\n",
       "      <th>admitted_urti_date</th>\n",
       "      <th>sgss_gp_cov_urti_date</th>\n",
       "      <th>event_urti_admitted</th>\n",
       "      <th>delta_urti_admitted</th>\n",
       "      <th>delta_admitted_died</th>\n",
       "      <th>delta_admitted_deregistered</th>\n",
       "      <th>delta_admitted_died_deregistered</th>\n",
       "      <th>duration_urti_admitted</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>59</td>\n",
       "      <td>NaT</td>\n",
       "      <td>2021-04-13</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2019-01-03</td>\n",
       "      <td>0</td>\n",
       "      <td>NaT</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>30.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>49</td>\n",
       "      <td>NaT</td>\n",
       "      <td>NaT</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2019-01-23</td>\n",
       "      <td>0</td>\n",
       "      <td>NaT</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>30.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>70</td>\n",
       "      <td>NaT</td>\n",
       "      <td>NaT</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2019-01-28</td>\n",
       "      <td>0</td>\n",
       "      <td>2021-04-21</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>814.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>30.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>8</td>\n",
       "      <td>NaT</td>\n",
       "      <td>NaT</td>\n",
       "      <td>2019-08-03</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2019-01-24</td>\n",
       "      <td>0</td>\n",
       "      <td>2020-08-22</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>576.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>30.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>5</td>\n",
       "      <td>NaT</td>\n",
       "      <td>NaT</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2019-01-05</td>\n",
       "      <td>0</td>\n",
       "      <td>NaT</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>30.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>284</th>\n",
       "      <td>41</td>\n",
       "      <td>NaT</td>\n",
       "      <td>NaT</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2019-01-28</td>\n",
       "      <td>0</td>\n",
       "      <td>2022-03-19</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1146.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>30.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>285</th>\n",
       "      <td>67</td>\n",
       "      <td>2022-04-02</td>\n",
       "      <td>NaT</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2019-01-13</td>\n",
       "      <td>0</td>\n",
       "      <td>NaT</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>30.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>286</th>\n",
       "      <td>34</td>\n",
       "      <td>NaT</td>\n",
       "      <td>NaT</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2019-01-24</td>\n",
       "      <td>0</td>\n",
       "      <td>NaT</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>30.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>287</th>\n",
       "      <td>13</td>\n",
       "      <td>NaT</td>\n",
       "      <td>NaT</td>\n",
       "      <td>2021-12-26</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2019-01-14</td>\n",
       "      <td>0</td>\n",
       "      <td>NaT</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>30.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>288</th>\n",
       "      <td>70</td>\n",
       "      <td>NaT</td>\n",
       "      <td>NaT</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>2019-01-08</td>\n",
       "      <td>0</td>\n",
       "      <td>NaT</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>30.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>239 rows × 15 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     age  died_date deregistered_date gp_cons_urti gp_cons_urti_ab  urti_date  \\\n",
       "2     59        NaT        2021-04-13          NaN             NaN 2019-01-03   \n",
       "3     49        NaT               NaT          NaN             NaN 2019-01-23   \n",
       "5     70        NaT               NaT          NaN             NaN 2019-01-28   \n",
       "6      8        NaT               NaT   2019-08-03             NaN 2019-01-24   \n",
       "7      5        NaT               NaT          NaN             NaN 2019-01-05   \n",
       "..   ...        ...               ...          ...             ...        ...   \n",
       "284   41        NaT               NaT          NaN             NaN 2019-01-28   \n",
       "285   67 2022-04-02               NaT          NaN             NaN 2019-01-13   \n",
       "286   34        NaT               NaT          NaN             NaN 2019-01-24   \n",
       "287   13        NaT               NaT   2021-12-26             NaN 2019-01-14   \n",
       "288   70        NaT               NaT          NaN             NaN 2019-01-08   \n",
       "\n",
       "     incdt_urti_date admitted_urti_date  sgss_gp_cov_urti_date  \\\n",
       "2                  0                NaT                      0   \n",
       "3                  0                NaT                      0   \n",
       "5                  0         2021-04-21                      0   \n",
       "6                  0         2020-08-22                      0   \n",
       "7                  0                NaT                      0   \n",
       "..               ...                ...                    ...   \n",
       "284                0         2022-03-19                      0   \n",
       "285                0                NaT                      0   \n",
       "286                0                NaT                      0   \n",
       "287                0                NaT                      0   \n",
       "288                0                NaT                      0   \n",
       "\n",
       "     event_urti_admitted  delta_urti_admitted  delta_admitted_died  \\\n",
       "2                    0.0                  NaN                  NaN   \n",
       "3                    0.0                  NaN                  NaN   \n",
       "5                    0.0                814.0                  NaN   \n",
       "6                    0.0                576.0                  NaN   \n",
       "7                    0.0                  NaN                  NaN   \n",
       "..                   ...                  ...                  ...   \n",
       "284                  0.0               1146.0                  NaN   \n",
       "285                  0.0                  NaN                  NaN   \n",
       "286                  0.0                  NaN                  NaN   \n",
       "287                  0.0                  NaN                  NaN   \n",
       "288                  0.0                  NaN                  NaN   \n",
       "\n",
       "     delta_admitted_deregistered  delta_admitted_died_deregistered  \\\n",
       "2                            NaN                               NaN   \n",
       "3                            NaN                               NaN   \n",
       "5                            NaN                               NaN   \n",
       "6                            NaN                               NaN   \n",
       "7                            NaN                               NaN   \n",
       "..                           ...                               ...   \n",
       "284                          NaN                               NaN   \n",
       "285                          NaN                               NaN   \n",
       "286                          NaN                               NaN   \n",
       "287                          NaN                               NaN   \n",
       "288                          NaN                               NaN   \n",
       "\n",
       "     duration_urti_admitted  \n",
       "2                      30.0  \n",
       "3                      30.0  \n",
       "5                      30.0  \n",
       "6                      30.0  \n",
       "7                      30.0  \n",
       "..                      ...  \n",
       "284                    30.0  \n",
       "285                    30.0  \n",
       "286                    30.0  \n",
       "287                    30.0  \n",
       "288                    30.0  \n",
       "\n",
       "[239 rows x 15 columns]"
      ]
     },
     "execution_count": 167,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_urti_no_abs_incdt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4- Cox modelling for hospital admissions related to infections \n",
    "### 4-1- URTI\n",
    "#### 4-1-1- Incident hospital admission related to URTI with no antibiotics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['age', 'died_date', 'deregistered_date', 'gp_cons_urti',\n",
       "       'gp_cons_urti_ab', 'urti_date', 'incdt_urti_date', 'admitted_urti_date',\n",
       "       'sgss_gp_cov_urti_date', 'event_urti_admitted', 'delta_urti_admitted',\n",
       "       'delta_admitted_died', 'delta_admitted_deregistered',\n",
       "       'delta_admitted_died_deregistered', 'duration_urti_admitted'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 168,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_urti_no_abs_incdt.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(239, 15)"
      ]
     },
     "execution_count": 169,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_urti_no_abs_incdt.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [],
   "source": [
    "#keep relevant variables\n",
    "# data_urti_no_abs_incdt = data_urti_no_abs_incdt[['age', 'age_cat', 'sex', 'bmi_imp', 'BMI_cat', 'gp_count_imp', 'flu_vaccine',\n",
    "#                                                'region', 'antibacterial_brit_imp', 'imd_imp', 'ethnicity_imp', 'CCI_cat', 'CCI',\n",
    "#                                                'event_urti_admitted', 'duration_urti_admitted']]\n",
    "\n",
    "data_urti_no_abs_incdt = data_urti_no_abs_incdt[['sex',\n",
    "                                                #  'age', 'bmi_imp', 'gp_count_imp', 'flu_vaccine', 'sex', 'region', #'season', \n",
    "                                                #  'imd_imp', 'antibacterial_brit_imp', 'ethnicity_imp', 'smoking_status_imp',\n",
    "                                                #  'meds_nonabs',\n",
    "                                                 'event_urti_admitted', 'duration_urti_admitted']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "age                       0\n",
       "event_urti_admitted       0\n",
       "duration_urti_admitted    0\n",
       "dtype: int64"
      ]
     },
     "execution_count": 171,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# printing number of missing values in each variable \n",
    "data_urti_no_abs_incdt.isnull().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #drop ros with missing values\n",
    "# data_urti_no_abs_incdt = data_urti_no_abs_incdt.dropna(axis='rows')\n",
    "# # data_urti = data_urti.drop(columns=['smoking_status'])\n",
    "# data_urti_no_abs_incdt.isnull().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['age', 'event_urti_admitted', 'duration_urti_admitted'], dtype='object')"
      ]
     },
     "execution_count": 173,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_urti_no_abs_incdt.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shape of data_urti_no_abs_incdt for training: (179, 3)\n",
      "shape of data_urti_no_abs_incdt for testing: (60, 3)\n"
     ]
    }
   ],
   "source": [
    "fold_thre = len(data_urti_no_abs_incdt)-round(len(data_urti_no_abs_incdt)/4)\n",
    "data_urti_no_abs_incdt_train = data_urti_no_abs_incdt.iloc[:fold_thre]\n",
    "data_urti_no_abs_incdt_test = data_urti_no_abs_incdt.iloc[fold_thre:]\n",
    "\n",
    "print(f'shape of data_urti_no_abs_incdt for training: {data_urti_no_abs_incdt_train.shape}')\n",
    "print(f'shape of data_urti_no_abs_incdt for testing: {data_urti_no_abs_incdt_test.shape}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "30.0    235\n",
       "28.0      1\n",
       "15.0      1\n",
       "13.0      1\n",
       "6.0       1\n",
       "Name: duration_urti_admitted, dtype: int64"
      ]
     },
     "execution_count": 177,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_urti_no_abs_incdt.duration_urti_admitted.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>model</th>\n",
       "      <td>lifelines.CoxPHFitter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>duration col</th>\n",
       "      <td>'duration_urti_admitted'</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>event col</th>\n",
       "      <td>'event_urti_admitted'</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>penalizer</th>\n",
       "      <td>0.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>l1 ratio</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>baseline estimation</th>\n",
       "      <td>breslow</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>number of observations</th>\n",
       "      <td>179</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>number of events observed</th>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>partial log-likelihood</th>\n",
       "      <td>-20.69</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>time fit was run</th>\n",
       "      <td>2022-04-11 13:49:32 UTC</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div><table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th style=\"min-width: 12px;\"></th>\n",
       "      <th style=\"min-width: 12px;\">coef</th>\n",
       "      <th style=\"min-width: 12px;\">exp(coef)</th>\n",
       "      <th style=\"min-width: 12px;\">se(coef)</th>\n",
       "      <th style=\"min-width: 12px;\">coef lower 95%</th>\n",
       "      <th style=\"min-width: 12px;\">coef upper 95%</th>\n",
       "      <th style=\"min-width: 12px;\">exp(coef) lower 95%</th>\n",
       "      <th style=\"min-width: 12px;\">exp(coef) upper 95%</th>\n",
       "      <th style=\"min-width: 12px;\">cmp to</th>\n",
       "      <th style=\"min-width: 12px;\">z</th>\n",
       "      <th style=\"min-width: 12px;\">p</th>\n",
       "      <th style=\"min-width: 12px;\">-log2(p)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>age</th>\n",
       "      <td>0.00</td>\n",
       "      <td>1.00</td>\n",
       "      <td>0.01</td>\n",
       "      <td>-0.01</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.99</td>\n",
       "      <td>1.01</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.21</td>\n",
       "      <td>0.83</td>\n",
       "      <td>0.26</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table><br><div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Concordance</th>\n",
       "      <td>0.62</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Partial AIC</th>\n",
       "      <td>43.39</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>log-likelihood ratio test</th>\n",
       "      <td>0.04 on 1 df</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>-log2(p) of ll-ratio test</th>\n",
       "      <td>0.26</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/latex": [
       "\\begin{tabular}{lrrrrrrrrrrr}\n",
       "\\toprule\n",
       "{} &  coef &  exp(coef) &  se(coef) &  coef lower 95\\% &  coef upper 95\\% &  exp(coef) lower 95\\% &  exp(coef) upper 95\\% &  cmp to &    z &    p &  -log2(p) \\\\\n",
       "covariate &       &            &           &                 &                 &                      &                      &         &      &      &           \\\\\n",
       "\\midrule\n",
       "age       &  0.00 &       1.00 &      0.01 &           -0.01 &            0.01 &                 0.99 &                 1.01 &    0.00 & 0.21 & 0.83 &      0.26 \\\\\n",
       "\\bottomrule\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "<lifelines.CoxPHFitter: fitted with 179 total observations, 175 right-censored observations>\n",
       "             duration col = 'duration_urti_admitted'\n",
       "                event col = 'event_urti_admitted'\n",
       "                penalizer = 0.2\n",
       "                 l1 ratio = 0.0\n",
       "      baseline estimation = breslow\n",
       "   number of observations = 179\n",
       "number of events observed = 4\n",
       "   partial log-likelihood = -20.69\n",
       "         time fit was run = 2022-04-11 13:49:32 UTC\n",
       "\n",
       "---\n",
       "            coef  exp(coef)   se(coef)   coef lower 95%   coef upper 95%  exp(coef) lower 95%  exp(coef) upper 95%\n",
       "covariate                                                                                                         \n",
       "age         0.00       1.00       0.01            -0.01             0.01                 0.99                 1.01\n",
       "\n",
       "            cmp to    z    p   -log2(p)\n",
       "covariate                              \n",
       "age           0.00 0.21 0.83       0.26\n",
       "---\n",
       "Concordance = 0.62\n",
       "Partial AIC = 43.39\n",
       "log-likelihood ratio test = 0.04 on 1 df\n",
       "-log2(p) of ll-ratio test = 0.26"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cph_urti_no_abs_incdt = CoxPHFitter(penalizer=.2)\n",
    "cph_urti_no_abs_incdt.fit(data_urti_no_abs_incdt_train,\n",
    "        duration_col='duration_urti_admitted',\n",
    "        event_col='event_urti_admitted',\n",
    "        formula=\"sex\"\n",
    "        # \"age + sex + bmi_imp + gp_count_imp + flu_vaccine + region + antibacterial_brit_imp + imd_imp + ethnicity_imp + smoking_status_imp + meds_nonabs\" # + season\n",
    "        )\n",
    "cph_urti_no_abs_incdt.print_summary()\n",
    "# cph_urti_no_abs_incdt_summary = cph_urti_no_abs_incdt.print_summary(style=\"html\")\n",
    "# cph_urti_no_abs_incdt.print_summary(columns=[\"coef\", \"se(coef)\", \"p\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2870\n"
     ]
    }
   ],
   "source": [
    "output = GetPrintSummary(cph_urti_no_abs_incdt)\n",
    "\n",
    "with open('../output/hospitalisation_prediction_urti/cph_urti_no_abs_incdt_summary.html', \"w\") as fo:\n",
    "    print(fo.write(output))\n",
    "    fo.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = cph_urti_no_abs_incdt.baseline_hazard_\n",
    "# cph_urti_no_abs_incdt.baseline_survival_\n",
    "# cph_urti_no_abs_incdt.predict_cumulative_hazard(data_urti_no_abs_incdt, data_urti_no_abs_incdt['duration_urti_admitted']<=30)\n",
    "# cph_urti_no_abs_incdt.predict_survival_function(data_urti_no_abs_incdt)\n",
    "\n",
    "# save \n",
    "np.savetxt('../output/hospitalisation_prediction_urti/sum_betas_urti_no_abs_incdt.txt', a, delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAd8AAAEdCAYAAABNIG3nAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA21klEQVR4nO3debyVZbn/8c9XcCrFDc6JBkWTaWFqPz2mIZmVR1ND07ABjcqgU3qSk0YlZjRplh218pBiKQ5pWmlWJFpYaeI8K7YpcSIZRM0Q9fr9cd8LHh7W2nvtzV4D2+/79Vqvvde9nuF6xuu572dSRGBmZmbNs06rAzAzM3u5cfI1MzNrMidfMzOzJnPyNTMzazInXzMzsyZz8jUzM2uyupOvpN0lXSLpUUnPS1ooaaakj0ka0IjgJK0j6XuSHpP0kqQrcvkbJc2StFRSSDpI0hRJPbpvStKo3P+oBoRfGcc4SUc1avhdjHd3STdKejZP48g6uu9y+RbmV+XznKR7JH1F0oaFYU2XNL/GePZp9DxfU3mZhaRhveg3JE2po7sDJN0p6d+5n45ehFpPPMPy8Mc1aPjTJc2ro7tez9MexHKMpA/U2W0lnhGNiqe3JM2TNL3wfVTev/WqstSM/VweT13rfg+Gt8p8aNdh9tbAejqSdAxwGjAL+ALwd2AwsC/wA2AJ8IsGxHcI8Dng88BfgIW5/DTgNcAH87jvB+YAv+nh8G8Bdgfu6YNYaxlHms/nNHAc1fwYeA44APgX8ECtDnuxfD8L3AS8AngPcCIwAvho305C/yRpIHAB8GdgIvA88HSDRvcYaR1/qEHDPxk4vUHD7qljgOuBn7c4jjV1MLC08H0UaRv7GvBSKwJqkfJ8aNdh9kq3yVfSXqQd8xkR8dnSz7+QdBrwykYEB7wp//1eRLxUKv9jRBST7WKgam2rlohYCtywZiG2n3yE/AZgakTM6qbb3izfeyOiMt9mSdoCGCfpmIhY1AeT0N9tA2wMXBIRf1zTgeWWCUXEC+XfImIZDVzHI6JRSf1lKyJubXUM7aAR86Gt5m1EdPkBrgKeBDbortvc/duB3wPPAM8C1wBvr9LdO/NvT+fufgvsUPh9HhClz7gqZZG7n1L5vzCMgaSa3D3Av4F/kmrHb8y/j8rDGFXq7wOkHda/SLW+nwHblbqZB5wPHA7cm6dhDvCOQjfXVYn3uvzbVsB5wKPAMlIN5Upgi27m7yDgjEJ/9wPHkna+1JhH8/pi+Rbm1z6l8gm5/O35+3Rgfo1h7FNtnlfpbkru7o153XgW+AdwZP79I8B9eT27Fnhtqf91STWFeaSa5bz8fd1Sd6/J8+Bfef04HfhUHvewUrefBG7P69KTpNaFIaVuAphSx3RVWyeUl+X9OebH8rIeVGUcU4HjgU7gRWCnGuMblrsfVyibTjpQ3QmYnaf9QeDoKv0PB34KPJ7Xt78Bp5eGNa8J8/RrpBaXTtI+4w/Am7vZX0zvYjmMy93sRmqFWErapr5PaVsAtgZ+kuNbBtwBfLjUTZfbc2E5TCAd7C7I8+fKKvNkXiX2GutLFLo9idSCtzTHNwvYrcZ22902tytwaV43niOth18HNix1NyAvj8fyNFwHvJnSus+ab8Mr5kOd83ggqSXmocL6dD2r7pNXGWYu6zZn0YNtpt5PdzvbAXkkM+oaGLwlL7SbSU3GY0jNk88Bby1095/AC6SmzAPz58+k2uu2uZudgHNZuYHsBrw6/11A2rh3q6xoVE++l+bxnAq8FziItOLvXWulBI7OZecA+wGHkZJrJ7BxaSH+PU/fIcD+wK2kZN2Ru9metGHcXpiG7fNvM0lNwUcAewGHAj+ktCGWpmedvOCfJTXF70vasQXw9dzN5sAeuWxaHmetHXNPl29lfpWT7ym5/PWlFXVglc++5XleY1xTcnd3kna67wYur0xrXl8OyvPtUeDGUv8z8rL/ah7nFGB5cVqB9Ugb6qPAkaT18pfAw5QSBfDN3P938vCOBB4BbgQGFLrrLvkOzetLkHYUxXXi67n8DFJz/rGkHcJsYJ3SOB7J5WNI6/aWNcY3jOrJdylpvf5Unrczcnd7F7obTkqefyclyb2BjwEXlIY1rwnzdB5pB/7+PP86gbnAwML+4jHSwXVlW3ttF8thXB7ug3kd2Qf4MulA5qRCd68kbaf/zPPgfaRkHcAnC911uT0XlsPDwK/yfDkyx/wAhYNCVk2+Q0nbcZC26xX7vPz7NFIS25u0D7qIdOC2Y5XttrttbgzwpTycd5IOFB4HLip1dzKp+fvUvNy+mJd5reTb2214xXyocx5PJm0vn8vxH0A6OHl/F8OsN2dNp45tpief7na2W+aBf6POnfOlFJJPLhsELAJ+XiibC1xT6ncQ6Ujle4Wyr1FKqLl8PqsfvUwpdguMzrF/tot4V1kpgY2Ap4BzSt0Nzyv0MaWFuBgYXCjbJQ9vbKHsOuD6KuN+pqvYasS7P6UdaWEDXAZslr8PpJsk0MvlW5lf++ZxDMor7DPAraUVNbr5dLcjmJK7+2ihbDApoS6kUBskbdgBvDp/36Ha9JN2LAG8JX//RP5e3JmtA9xNIVGQdpwvAl8pDa9ykHNQoaye+T6ivByBIXkZltfrD+du318ax6OUaiQ1xjWsyrgqy6eYaNfP8/XsQtlP8rJ9VRfDn86qybdR8/RBVk1QlQOY/yhtk+fXuS6Py/2fVCq/Enig8P0z1dZXUk1pAfkggW6258JyuIdVD6Qq0/vx0nRML3yfkrsZ2M00DSBtl/ezauvEqGrT0M2wlIf1YVKi3bSwDT4D/LDU/ReonXx7vA3XmA/dzeMrKeSZGt2Uh1lvzppOHdtMTz59favRXsCVEbGkUhDpvOovSUciSHod8FrgAkkDKx9SDewveRh9oVLD+r8e9LM7acaXY3uY1DxSju0vEbG48P3O/He7OsZ1EzBJ0uck7ShJdfSzF2lDmFEqP59U49i9jmH0hd+SaixPkZrkryUdwRYtIDVjlT8Teziuqyv/5Hm9ALghr1cV9+W/2+a/leV0fmlYle/vzH93Bx6OleeviXRtwSWl/t5NSiDl9eJGUhNoX6yzu5GWYTnmi0g7q3eWyn8TEc+twfj+FRHXVr5EOjf8AKuuu/uStudHezDcRs3TmRGxvPC9J9taV64qfb+zNMy9gEci4rpSd+eTWpm2z9/r3Z4vjcL1KxHxJ1Jlolfbbr574FpJC0nryXLg9aRrPno6rEGSviXpIdKB4HLSKQcBr8ud7UhqDSgvz4u6GHRvtuFqupvHNwH7SZoq6R2S1utiWBXd5qyCeraZunWXfBeSqt+vrnN4Q0jNKGWPk454ALbIf39MWrjFz/7ApnWOqzubAot6uIOqxPb7KrHtWCW2VS4uygsDYIM6xnUYaQH/D+kc0iP5lp2ulskQ0jQ9Xyp/vPB7T/R0+VZMJCXSHYCNIuKAiPh7qZvlETGn/KGLq65rWFz6/nyNMlg53yvzobwulufT1sATVcZZLqusF3NZfb3YmL5ZZ6vGHOkiqoWsvmyrbWc9UZ6HkHa4xXV3U3p4ESONm6flC/l6sq11pdpw1y9872qfVvkd6t+ea82bbXoSNICktwG/JtUIP046gNuVdJqrN/PlXNJpt++TDo6KB8uV4W1diLmo2nRV9GYbrqa7efx10lXh7yedklko6VxJm3UxzHpyVq3pgNW3mbp1ebVzRLwg6Trg3ZLWLySXWhaRToqXbcXKwCu3C51ASnJl5cTSW08CQyRt2IMEXIltHKmZrKzPbgeJiAWkFXuipDeQzqWdRDq39IMavS0iTdN6pQS8VeH3nsTQ0+Vb8UBOpO2qMh+2YtVbbMrz6THShSJlW5a+V9aLfam+AS6sUtZTxZhXrHu5Nrgpqy/b6INxdudJep4U2mme9oVFVK9FrrIu9WB7Ls+HStltvYhtDKm2+4Fiq4CkwaSm1LpJ2oB07c2UiDi9UL5jqdNKotqSVfeR1aarT3U3j/M8+BbwLUlbkSpzp5FuiTysxmDryVkNUU+z8zdJG/+3q/0oabikt+SvfyBV+zcu/L4x6cT3dbnoflK7+5ur1Ywi4o7eTcpqfkdqLhnfg37+TEqwI2rEdn8v4lgGbNhVBxFxf0R8kbSwd+ii0z+QltmhpfIjSActf+lFfD1ZvmuLyu07h5fKj8h/r8t//wJsK2m3Sgf5KPqDpf5mkpr7t6uxXnT2Qcw3kJZhOebDSAfJ15V7aILfAftL2rrbLldq5TztdlvrhT8AQyXtUSofS2o+Xe0ZAd1sz4cUa8N5uEPpetutHBSXp+0VpPPmKw7EJI2md82g65POGS8vlY8rfb+DdMFneXmW19uG6m6fGRGPR8Q0UgWvu31qdzmrIbq9zzci/ijpv4HTJG1POvH8D1KV/F2k5DaWtFBOJh1tXCPpW6SV4gukleSreXghaSLpHtL1SOcOniQdOf0H8I+IOG1NJywirpV0WY57W9Il+OuS2vivqnIOh4hYKmkScKakzUnnKp4iHf2/k3RLSPl8a3fuASZIOoxUC3ua1KTxe9JVk/eRVvgDSfP0d10M62rSpfM/zPHdTboiezzpoqknexhbT5fvWiEi7pJ0ITAl1xz/TDqn9mXgwoionC88j3S7zs8lfZG0Mz2adN6/OLyH8vp8Rj7i/gPpVoZtSc1z04rngnoZ8yJJ3wFOkPQsqTnxTaSLDq9n9XOTzXAiaf36s6Svk5qItwHeGxEfrtFPK+fpPcCekvYnbWNPRsS8Hg6jbDrp6tmfS5pMaoY/Isf4qYh4UdIm1L89bwxcIelHpHPG3yBdTPaTbqYL4POSrgZezC1PvyE9WGS6pHNJ53q/TLpivEci4ilJN+RxPEbaJx9FqeUjIpZI+i4wWdLTefp2JTV7N0w981jSL0hN7reQkvJOpDsBftTFoLvNWQ1T75VZpMT4M1Kzw3JSdf13pKvhilfv/T/qu893d9LVaYtJG9080kn73Qvd9Ppq51w2kHT5+QOkWsU/STu1N0QXVwGSdjjXki4tr9zPdQ75lpBYedXcaldWsvoVf1vlcT6df7uOdJT5I1LyfCaP5yYKV0l3sRwq9/k+lqfpAQr3+Rame5U4+mL5UuNWoyrDmk7f3ec7sFS+2nyvFhfp4qWvkW6TWZ7/1rrP99fUd0/qR0g11Gfzcrs3L4uhtZZ/jWlb7WrnXF7tPt8zqX6f79fqXK7DyuOqtXzyunldqey1wIWknfG/SQeQp5WGNa8J8/RrdUzXG1l5D2ZQ332+I+rYj2xNuvCo6n2+1LE9s/p9vv/McV4FDK+yjk8vfB+Q14MFpNaCKPz2X6Tbrp7L49ynvByp/1ajYaQD/KfzuM4g3RK1Sr+svM/38Tze60gXnpX3fVNYs214xXyocx5/Pq9LlWtZ7s8xVL2Nq1DWbc6iB9tMvZ/KgxnMzKxBlJ5p3Ql8IlJzqL3M+a1GZmZmTebka2Zm1mRudjYzM2sy13zNzMyazMnXzMysyZx8X0YkDZMUksb1ot9jJH2gATHNkzS9m25GSpoiqaePz+xqmKPyvNini7jOL3wfl7uvfJ6X9JCkr+enA1H6vdZnXu52Sv7e7b32vZi26ZXx5O/D8vhe09109nddrUt5eUwpfK9rGRXWpVE9jKUjj+NtVX67Lj99zvqpPt/wrd86hvSwh5+3YNwjSQ98OJ8ePkKzAQ4l3We+MXAw6TGpG5Putyw/HP9y0k3/Uwpl9T7Cc02cTLq3tmIYaf5dT3of78vZSGqvS7vT82dZQ3qow+5UedpVNzpyLPPzMIom9CIOW4s4+a7levhMZltzt0XE3Pz/TKW3dB0l6XNReJMPgKRlpKcs3bDaUBooIh7qvisr6+1yivQWnD5dxhHR00Ruaxk3O69FCk2fe0n6maQlpFewofRKthMk3SdpmaRHJX2n0iTaxTB3lXSppPmSnpN0f25K3bDQzTzSm4+OKDSfTi/8/lZJv5S0OA/jT5L2rDKuz+Vmzn9LmlOtm2rTTHrbCsCDhfEPy78PknRGnt5lOf5jpbpe0dgXbiE9iq6rN6f0iqTPS/qXCq9Gk3RZublc0ickvSBpUP6+otk5N4VWHtM4szD/RpXGdbikeyU9m5fNO+qIr9Is+zpJV0l6RtLfVeVtPpLeIOlySUvyOnKDpPfWMY7NJf1I0gN5XjwsaYakbUrddRtLHevSKs3OBW9Sem3fvyQ9JumrWvX5zKs1Oys5Nq+Pz+f+zigso2Gkh24A/F8hlnH599WanfO8OCvPg2X5708lrZ9/f32exwvyNvYPpf2EK1ltyMl37XQBacM9hPQcXUjNaF8ivev3P0nPjP147rYr25HeqHI06Tmop5Oe6XpuoZuDSY+S+y2peW13UtNm5bVmfya9musTpDetLAR+L2nnygAkfRz4Hivf/Tud9NjC8mu7yq4iPcoOUpNvZfyP5R3gVcCRwHdID0P/DenxfVO7GW5fGUZ6/ncj3sJzLelh+rtB2qGTHsP3HDC60N1o4OZY9f2oFbew8rVwn2Xl/Cs2c+5JejTfl0kvchgAXCmpo844Lyc9O/0g4ArSm2Y+VvlR0qtITd5vJb2c/oOkt+5cJel93Qx7COmxlieQ1s9JpHfL/knVDyy7iqXmutRNDFeQHj94EGn7+jLwlW76mUpaD2eS1stvkx5peVVebx8DKtdQfKMQS9VneCu9qejPpOVzGukRuP9Del595eDsKtKzmD8NvIe0b1iG9/PtqTfPpPSnNR9WPo/2u6XyPXP5R0vlR+Tykfn7MKo8U7jQvUinIj5MeobspoXf5lH9WdbXkJ7Hu16hbEAuuyJ/Xwd4mPQC+GK/h9HNM3hL011+Du/+1aYHmEba6WzWxTBH0cVzqsvTW4jhDXkeDSYdpLwAfKaeYZR+m0KV596WulmHdF7yxPx9ZF4upwF/KXT3GPDNwvfpFJ633NW05hgXA4MLZbvk7rt81nhhGo4sld8J/K7w/dQ8n0YUygaQnr17Sw+3gQGkly8EcHAvYqm6LuXfaj2b+PhSd/9Hev5xR2n+jsrfh+T1b3qpvw/n7t5f2h7HV4nlOlZ9PvNXSW8w2qnGfNmsOGx/2v/jI6K10+Wl7+8lPYj/UqXm54G5qanyRpW9ag0oN9t+S9JDpB3GctJD5EWqYdSk1DT9TtILGV4qjFekmkJlvEPz55LSIC4j7ZR7ay9SMiq/aep8Um2gfAFUX6i8UWUR8GPgRxFxRgPGQ0S8RHrbT6WWO5r0UP+fAbtI2ljpTVRbsbJpuTf+EhGLC98rb32q99V05draXaV+9wJuiJXnyomIF0ktHyMrTbG1SPq0pNslPUNaX/6Rf6r2nt3uYumN8np7EbARtV9Vtxtp/StfRX4RKf539iKGfYGbIuLWGr8vJF1M9818GqLLbddaz8l37VRuJtuCtLE/S0oMlc+C/PumXQzrXFKT8/dJr0nblZXNlF2eLyYd4Q8gNcMtL30+AwzOTWyV98E+Uew5Il5gzZprhwCLIuL5Uvnjhd9rqST9ATV+H0D1A4ODSfNoP9IBxgRJH60v3F65FtgtH+jsnb/fRGqK3TOXLSc16/bWKlf9xsoL+Lpb/lX7Jx3EFfsdQvWm3cdJB2o1Tz1I+i/gLNK8/gDwdnIzfI34uoulN56o8X2bcodZZb1bZZoL63tvbpnblC6uxI5U/X03MIfUjP2ApL9J+nQvxmVN4BPxa6fyM0EXsnJnXM2j1QrzObMDSU1tpxfKd6wzjiWkmueZ1HgfaUS8pPR+UEjvbC6OfyBdHxh0ZxEwRNJ6pQS8VeH3WioHJq8q/5Dj2oLVd7oAd1VqcJJmkWqip0i6LCKe7ekE1OFa0oHVXvlzdkS8IGk2qSY8HPhrg8bdVxaxcpkUbUValxdX+a3icOCaiPh8pUDS8L4Nr1tbsuotWpX1uNZ7cyvr3VakV+ABq6zvvbld7klqJ3sAIuJvwEfztQGV8+tnSZoXEVf3YpzWQK759g+/IR3dbxIRc6p8qiZf0jsyB5BqTkXjqnS7jHTxzwp5hz+btKHfUm3cudP5pHO+HywNcwz1HQBWamIblsr/QFqHDy2VH0Fqhv9LF8N8MMdV7cEh+5MSXpdNubmGOImUqBt1X+ZdpHe/TgJeSZpmSBcVvYvUhNldk3Ot+dcsfyDV3odVCiQNIJ3zvzWqXyhW8QpWXz+PXINYejMvyuvt4aR3v95ZpVtItx09n7srOoy0vl/Xi1h+B7xd0lu76zCS24D/zkW1msethVzz7Qci4jpJF5LO+Z4G/JVUIx1Gah79QkQ8UKW/pyTdAHw+106fJF1EVO0I+x5gT0n7k5oLn4yIeaQN/I/AbyX9mNTUthnwNmBARByfa78nAdMknUs69zWCdDVmVzve4rgBJko6j7QzvoP04u/rgR9K2pxUy9gPGA98IyKe7GKehaQTgJ9Kuox03ngpqUn5i6Tk9tvuAouIX0q6iTQPz4iI5+qYnrrlOK8jHWDcVEhU1wKn5P9ndTOYB0hN6EdJWkTa6d8fEU/3Zaxd+C7pgG6mpBNJ83kC8HrSlfld+Q3wBUlfJK3Xo0lX+fdW1XWpyqmLok/k0yc3ka4iHk9qLXqqWscRsUjSd4ATJD0L/Bp4E+lK6+tZeV76CVKr1eGS7iCdNuqMiGqnYr4LjCXdRfA1UuLfjNRydTSpBeR04GJgLumgehxpuXe3flgrtPqKL3/q/9D1lZrrAJ8jPVHp36TbX24n3eKwSe5mGKWrg3PZ1aSrNxcAZ5B2iCuu3szdvZFUy/0XpSuUSTuWi3L/y0g1yl8C+5Vi/Bzw9xzfHOAdpKttp9cx7SeSmvlezOMflssH5ZgfI9U2HgCOJb+xq47hHkCqmT2d+38Q+DqwYQ/m/b75t2NL5fNYg6udC91+OndbvKK5ciX0v4ENSt1Pp3C1cy77FKnp9IXisq0VI6Urf3syDTXG/wbSLTtP5ZhvAN5bx7RvCPyAVPt/GriSlGhqXZlcTyy11qVaw9yBdLDzHOnA82RgnUJ3o1h9e1FeD+/P69VjpNMzg0qxHEQ6IFhOYdukdLVzLtsCOJuV6/rDwHmkFqwt8v8PkLbRRaT1+j1rut/xpzEfv1LQzMysyXzO18zMrMmcfM3MzJrMydfMzKzJnHzNzMyazMnXzMysyZx8zczMmszJ18zMrMmcfM3MzJrMydfMzKzJnHzNzMyazMnXzMysyZx8zczMmszJ18zMrMmcfM3MzJrMydfMzKzJnHzNzMyazMnXzMysyQa2OgCzVttss81i2LBhrQ7DCjo7Oxk+fHirw7CCm2+++cmI2LzVcfQXTr72sjds2DDmzJnT6jCsYOzYscyYMaPVYViBpL+3Oob+xM3OZmZmTebka2ZtZ/To0a0Owayh3OxsZm1h0aJFLF26FIB99tmHefPmATBo0CCGDBnSwsjM+p5rvmYNtmTJEqZMmcKSJUtaHUrbWrRoESNGjGD48OGrfUaMGMGiRYtaHWJb8zq29nHyrYOkgySFpDe2OI4OSRMK30dJekrSryXtKOm2/FkkqTP///vSMM6RtEDSXaXykZJuyP3MkfT2GjFMLwz7Nkkjc/lgSZdLukPSXyXt0IBZ0COSDpM0V9KVrYxjyZIlnHTSSd4xdmHp0qUsXryY2bNn09nZueIze/ZsFi9evKJGbNV5HVv7uNm5Ph8Crs9/T2xhHB3ABOCsQtnsiNg//z8SUoIEroyIS6sMYzpwBvCTUvm3gZMi4mpJ++Xvo2rEManKsL8I3BYRB+eDlDOBd3U/SWtO0oCIeLFcHhEXS3oCOK4ZcXRn/vz5rQ6hbVXmzdChQ6l225fnXdc8f9Y+Tr7dkLQR8A5gb+BX5OQraR1SEhsNPAwsB86JiEsl7QycBmwEPAmMi4jHagx/BPBDYHPgReBQ4AngF8BgYF3gSxHxC+CbwGsl3QbMBK7q6fRExB8lDav2EzAo/78J8GgPB719jo+IuE/SMElbRsQTlQ4kHQW8JSKOyd8/AWwfEcdK+jDwWWA94EZgQkS8KOkHwK7AhsClEVGZ//OAi4F3A9+WtAVwNPACcE9EHN5VsJI+CXwSYLvttuvhpPbOnnvu2ZTx9Eeed9bfOPl270DgNxHxgKSFknaOiJuBDwDDSElnC+Be4BxJ6wL/CxwYEf+UdBgwFTiqxvAvAL4ZEZdL2oB0KuB54OCIWCppM+AGSb8Ejgd2iIiRkJqd+3A6jwF+K+nUHMN/dNHtVElfAa4Bjo+IZcDtpHkyOzdZvxoYSjqQqLgEmCxpUkQsB44EPiXpTcBhwB4RsVzSWcARpNr55IhYJGkAcI2kt0TEHXl4CyPibQCSHgWGR8QySR3dTWxEnA2cDbDLLrtEt3OnD8yePZuhQ4c2Y1Rrnfnz53eZYD3vutbd/LP24+TbvQ8Bp+f/L8rfbybVhn8WES8Bj0u6NnfzBmAHYKYkgAFArVrvxsA2EXE5QET8O5evC3xd0l7AS8A2wJZ9P2mr+DRwbERcJumDwI+Bfap0dwLwOKmGejbwBeCrpFrv6blWfidwK6kmv0JEPCNpFrC/pHuBdSPiTkmfAXYGbsrzbENgQe7tg7mWOhDYmnSwU0m+FxcGfwdwgaQrgCt6OQ8aqlaTqq1Ubj7trjnabG3l5NsFSUNIzco7SgpSIg1Jk7rqDbg7InZfg1EfQWqG3jnXBOcBG6zB8OrxMeBz+f+fAdMAJP2WlPjnRMT4QvP5Mknnks+nRsRSUk0WpQzaCfytynimkc4P3wecm8sEnBcRJxQ7lDQ8D3/XiFicz2UX58Ozhf//E9gLOIBUu94xIl7o0RxokI6ODk488UQ6OjpaHUrbGjRoEIMHD65aexs8eDCDBg2q0pdVeB1b+/hq564dAvw0Il4dEcMiYltSUtkT+BMwRtI6krZk5cVJ9wObS9odUi1W0purDTwingbmSzood7u+pFeQzrkuyIl3b1ITLsDTwMaNmFDSOd535v9HAw/mGN8TESMjYnyOcev8V8BBwF35e4ek9XL/44E/5oS8ioi4EdgWGAtcmIuvAQ7J522RNETSq0nnoJ8Fnsrz+H3VAs/n37eNiGtJNfFNSOfb20JHRwdTpkzxjrELQ4YMYe7cuSuucj7ggANW/D937lzf59sNr2NrH9d8u/Yh4Fulssty+UTS1bz3kC64ugV4KiKel3QI8H1Jm5Dm8feAu2uM4yPAjyR9lXTR1qGk88C/knQnMIdUSyQiFkr6U75N6GrquOBK0quAaRGxX/5+IelAYTNJ84ETI+LHwCdIzcYDgX+TL0aq4gJJm5Nqq7eRLnICeBNwXm4huBv4eBdhXQKMjIjFebrukfQl4Hc5kS4HJkbEDZJuzdP/MOmAp5oBwPl5fgv4fkQs6WL81oaGDBmyIslutNFGbma2fk0RTbnWpF+StFE+j7kp8FfSBUOPN3H8o4DjCrcarRXyfbffjYhrmjCuUXQzj3bZZZfwixXai99q1H4k3RwRu7Q6jv7Czc5r5sp8gdFs4ORmJt7seWAHSb9u8nh7JTdNPwA816TEexjpnujFjR6XmVlPuObbJJLOBPYoFZ8eEedW696axzXf9uNXCrYf13z7ls/5NklETGx1DGZm1h7c7GxmZtZkTr5m1nbGjBnT6hDMGsrJ18zajpOv9XdOvmbWdiZMmNB9R2ZrMSdfM2s7fi+t9XdOvmZmZk3m5GtmbcdPt7L+zsnXzNrO1KlTWx2CWUM5+ZpZ25k2bVqrQzBrKCdfM2s7s2bNanUIZg3l5GtmZtZkTr5mZmZN5uRrZm3nzDPPbHUIZg3l5Gtmbaezs7PVIZg1lJOvmbWdU089tdUhmDWUk6+ZmVmTOfmamZk1mZOvmbWd8ePHtzoEs4Zy8jWztjN69OhWh2DWUE6+ZtZ2xo4d2+oQzBrKydfMzKzJnHzNzMyazMnXzNrOTjvt1OoQzBrKydfM2s6kSZNaHYJZQzn5mlnbOeWUU1odgllDOfmaWdu59dZbWx2CWUM5+ZqZmTWZk6+ZmVmTOfmaWduZMWNGq0MwaygnXzNrO7NmzWp1CGYN5eRrZm1n2rRprQ7BrKGcfM3MzJrMydfMzKzJnHzNrO0cd9xxrQ7BrKGcfM2s7QwfPrzVIZg1lJOvmbWdiRMntjoEs4Zy8jUzM2syJ18zM7Mmc/I1s7YzevToVodg1lBOvmbWdsaPH9/qEMwaysnXzNrO5MmTWx2CWUM5+ZpZ2+ns7Gx1CGYN5eRrZmbWZE6+ZtZ2Ojo6Wh2CWUM5+ZpZ2znrrLNaHYJZQzn5mlnbueyyy1odgllDOfmaWdtx8rX+zsnXzMysyZx8zczMmszJ18zaztSpU1sdgllDOfmamZk1mZOvmbUdP17S+jsnXzMzsyZz8jUzM2syJ18zaztjxoxpdQhmDeXka2Ztx8nX+jsnXzNrOxMmTGh1CGYN5eRrZm1nyZIlrQ7BrKGcfM3MzJrMydfM2s7w4cNbHYJZQzn5mlnb8eMlrb9z8jWztjNt2rRWh2DWUE6+ZtZ2Zs2a1eoQzBrKydfMzKzJnHzNzMyazMnXzNrOmWee2eoQzBrKydfM2k5nZ2erQzBrKCdfM2s7p556aqtDMGsoJ18zM7Mmc/I1MzNrMidfM2s748ePb3UIZg3l5GtmbWf06NGtDsGsoZx8zaztjB07ttUhmDWUk6+ZmVmTOfmamZk1mZOvmbWdnXbaqdUhmDWUk6+ZtZ1Jkya1OgSzhnLyNbO2c8opp7Q6BLOGcvI1s7Zz6623tjoEs4Zy8jUzM2syJ18zM7Mmc/I1s7YzY8aMVodg1lBOvmbWdmbNmtXqEMwaysnXzNrOtGnTWh2CWUM5+ZqZmTWZk6+ZmVmTOfmaWds57rjjWh2CWUM5+ZpZ2xk+fHirQzBrKCdfM2s7EydObHUIZg3l5GtmZtZkTr5mZmZN5uRrZm1n9OjRrQ7BrKGcfM2s7YwfP77VIZg1lJOvmbWdyZMntzoEs4Zy8jWzttPZ2dnqEMwaysnXzMysyZx8zaztdHR0tDoEs4Zy8jWztnPWWWe1OgSzhnLyNbO2c9lll7U6BLOGcvI1s7bj5Gv9nZOvmZlZkzn5mpmZNZmTr5m1nalTp7Y6BLOGcvI1MzNrMidfM2s7fryk9XdOvmZmZk3m5GtmZtZkTr5m1nbGjBnT6hDMGsrJ18zajpOv9XdOvmbWdiZMmNDqEMwaysnXzNrOkiVLWh2CWUM5+ZqZmTWZk6+ZtZ3hw4e3OgSzhnLyNbO248dLWn/n5GtmbWfatGmtDsGsoZx8zaztzJo1q9UhmDWUk6+ZmVmTOfmamZk1mZOvmbWdM888s9UhmDWUk6+ZtZ3Ozs5Wh2DWUE6+ZtZ2Tj311FaHYNZQTr5mZmZN5uRrZmbWZE6+ZtZ2xo8f3+oQzBrKydfM2s7o0aNbHYJZQzn5mlnbGTt2bKtDMGsoJ18zM7Mmc/I1MzNrMidfM2s7O+20U6tDMGsoJ18zazuTJk1qdQhmDeXka2Zt55RTTml1CGYN5eRrZm3n1ltvbXUIZg3l5GtmZtZkTr5mZmZN5uRrZm1nxowZrQ7BrKGcfM2s7cyaNavVIZg1lJOvmbWdadOmtToEs4Zy8jUzM2syJ18zM7Mmc/I1s7Zz3HHHtToEs4Zy8jWztjN8+PBWh2DWUE6+ZtZ2Jk6c2OoQzBrKydfMzKzJnHzNzMyazMnXzNrO6NGjWx2CWUM5+ZpZ2xk/fnyrQzBrKCdfM2s7kydPbnUIZg3l5Gtmbaezs7PVIZg1lJOvmZlZkzn5mlnb6ejoaHUIZg3l5Gtmbeess85qdQhmDeXka2Zt57LLLmt1CGYN5eRrZm3Hydf6OydfMzOzJnPyNTMzazInXzNrO1OnTm11CGYN5eRrZmbWZE6+ZtZ2/HhJ6++cfM3MzJrMydfMzKzJnHzNrO2MGTOm1SGYNdTAVgdgZgawaNEili5dCsDOO+/MvHnzABg0aBBDhgxpYWRmfc81X7NuLFmyhClTprBkyZJWh9JvLVq0iBEjRjB8+PDVPiNGjGDRokWtDrHf8vrdGk6+fUDSQZJC0htbHEeHpAmF76MkPSXp15J2lHRb/iyS1Jn//30PxzG90O9tkkbm8sGSLpd0h6S/StqhjyevxyQdJmmupCvXZDhLlizhpJNO8s6pgZYuXcrixYuZPXs2nZ2dKz6zZ89m8eLFK2rE1ve8freGm537xoeA6/PfE1sYRwcwASi+EmZ2ROyf/x8JKYECV0bEpeUBSBocEYu7Gc+kKv1+EbgtIg7OByFnAu/q8RT0gqQBEfFiuTwiLpb0BHBcX4xn/vz5fTEYq6Iyb4cOHcqwYcNq/m59z/O2NZx815CkjYB3AHsDvyInX0nrAGcAo4GHgeXAORFxqaSdgdOAjYAngXER8ViN4Y8AfghsDrwIHAo8AfwCGAysC3wpIn4BfBN4raTbgJnAVb2YpDmSbgB+DFwbEVFnf9vn8RMR90kaJmnLiHiiMC1HAW+JiGPy908A20fEsZI+DHwWWA+4EZgQES9K+gGwK7AhcGlEVObvPOBi4N3AtyVtARwNvADcExGHdxWspE8CnwTYbrvt6prAPffcs67urO953lt/4+S75g4EfhMRD0haKGnniLgZ+AAwjJSUtgDuBc6RtC7wv8CBEfFPSYcBU4Gjagz/AuCbEXG5pA1IpwqeBw6OiKWSNgNukPRL4Hhgh4gYCanZuRfT83rgfcBngDMl/RSYHhGPFrqZKukrwDXA8RGxDLg9T/NsSW8HXg0MJR0oVFwCTJY0KSKWA0cCn5L0JuAwYI+IWC7pLOAI4CfA5IhYJGkAcI2kt0TEHXl4CyPibXlaHwWGR8QySR3dTWREnA2cDbDLLrvUdYAxe/Zshg4dWk+n1kPz58/vMsF63jdOd/PeGsPJd819CDg9/39R/n4zqTb8s4h4CXhc0rW5mzcAOwAzJQEMAGrVejcGtomIywEi4t+5fF3g65L2Al4CtgG27IuJyc23VwJXStoc+AbwD0n/ERF/BU4AHifVUM8GvgB8lVTrPT3Xuu8EbiXV1IvDfkbSLGB/SfcC60bEnZI+A+wM3JTnyYbAgtzbB3MtdSCwNelgppJ8Ly4M/g7gAklXAFf0xbwoq9Ukan2n3ATaXXO02drKyXcNSBpCalbeUVKQEmlImtRVb8DdEbH7Goz6CFIz9M65pjgP2GANhrcKSZsAhwPjSLXso8gJr9A8vkzSueTzqRGxlFSTRSmDdgJ/qzL4aaTzw/cB51ZGCZwXESeU4hieh79rRCzO56qL0/ls4f//BPYCDiDVrneMiBd6Ou3VdHR0cOKJJ9LR0dEXg7MqBg0axODBg6vWwAYPHsygQYNaENXLg9fv1vDVzmvmEOCnEfHqiBgWEduSks6ewJ+AMZLWkbQlMCr3cz+wuaTdIdViJb252sAj4mlgvqSDcrfrS3oFsAmwICfevUlNvABPAxuvyQRJOh+4BRgOfDQi3hkRPynUurfOfwUcBNyVv3dIWi8PZjzwx5yQy9N0I7AtMBa4MBdfAxySz9siaYikVwODSAn2qTwP31cj5nWAbSPiWlJNfBPS+fQ+0dHRwZQpU7xzaqAhQ4Ywd+7cFVc5H3DAASv+nzt3ru/zbSCv363hmu+a+RDwrVLZZbl8Iulq33tIF1zdAjwVEc9LOgT4fq5hDgS+B9xdYxwfAX4k6auki7YOJZ0H/pWkO4E5pFokEbFQ0p8k3QVcTR0XXEl6FTAtIvbLRZeQLgCrVWu8IDdHC7iNdJETwJuA83ILwN3Ax7sY7SXAyMpV1RFxj6QvAb/LiXQ5MDEibpB0a56+h0kHNNUMAM7P81PA9yNiSTeTbm1myJAhK5LsRhtt5GZm69dU/8Ws1lOSNsrnOTcF/kq6oOjxJo5/FHBc4VajtpDvu/1uRFzThHGNopt5sMsuu8ScOXMaHYr1wOLFixk8eHCrw7ACSTdHxC6tjqO/cLNzY12ZL0CaDZzczMSbPQ/sIOnXTR5vVblp+gHguSYl3sNI9zx3d9+ytZnOzs5Wh2DWUK75tglJZwJ7lIpPj4hzq3Vvfcc13/YzduxYZsyY0eowrMA1377l5Gsve5L+Cfwd2Iz00JO1WX+YBvB0tJvNgFdGxOatDqS/cPI1yyTNWduP7PvDNICno930l+loJz7na2Zm1mROvmZmZk3m5Gu20tmtDqAP9IdpAE9Hu+kv09E2fM7XzMysyVzzNTMzazInXzMzsyZz8rV+S9Khku6W9JKkXUq/nSBprqT7Jb2nUP7eXDZX0vGF8uGSbszlF1deIpFfdnFxLr9R0rAGT9MUSY9Iui1/9iv81ifT1Gq14m0XkuZJujPP/zm5bIikmZIezH8H53JJ+n6eljskva0wnI/l7h+U9LEmxH2OpAX52e+Vsj6LW9LOeb7Mzf2q0dO0VosIf/zplx/Syx7eAFwH7FIo3x64HVif9Pamh0gvZxiQ/38N6X3FtwPb534uAQ7P//8Q+HT+fwLww/z/4cDFDZ6mKaRnVZfL+2yaWrzMasbbLh9gHrBZqezbwPH5/+OBb+X/9yO95ETAbsCNuXwI6ZWbQ4DB+f/BDY57L+BtwF2NiJv0/Prdcj9XA+9r9bJq549rvtZvRcS9EXF/lZ8OBC6KiGUR0QnMBd6eP3Mj4m8R8TxwEXBgPoIfDVya+z+P9DrFyrDOy/9fCryrRUf8fTlNrVQ13hbHVI/ielBeP34SyQ1Ah9JrOd8DzIyIRZHe7jUTeG8jA4yIPwKLGhF3/m1QRNwQKRP/hPZYn9qWk6+9HG1DekVhxfxcVqt8U2BJrHzNYqV8lWHl35/K3TfSZ3JT4DmVZsIuYu/NNLVSrXjbSZBef3mzpE/msi0j4rH8/+PAlvn/ni6XZuuruLfJ/5fLrQa/z9fWapJ+D2xV5afJEfGLZsfTF7qaJuAHwMmkBHAy8B3gqOZFZ8A7IuIRSVsAMyXdV/wxIkLpvdZrlbU17rWVk6+t1SJin1709giwbeH70FxGjfKFpGa3gbmmWOy+Mqz5kgYCm+Tue63eaZL0f8CVpTjKsVOjvKtpaqWupqMtRMQj+e8CSZeTmsqfkLR1RDyWm2AX5M5rTc8jwKhS+XUNDr2avor7kfx/uXurwc3O9nL0S+DwfKXycOB1pItFbgJel68CXo90AdUv8zmsa4FDcv8fA35RGFblis9DgFm5+4bIO8iKg4HKlat9OU2tVDXeFse0gqRXStq48j+wL2kZFNeD8vrx0Xz18G7AU7mZ97fAvpIG51MH++ayZuuTuPNvSyXtlq8n+CjtsT61r1Zf8eWPP436kJLTfGAZ8ARpJ1H5bTLpqtr7KVyVSbrK84H82+RC+WtIyWwu8DNg/Vy+Qf4+N//+mgZP00+BO4E7SDvIrft6mlr9qRVvO3zyPLs9f+6uxEc6h34N8CDwe2BILhdwZp6WO1n1qvuj8ryfCxzZhNgvBB4Dluft4uN9GTewC+lA5CHgDPITFP2p/vHjJc3MzJrMzc5mZmZN5uRrZmbWZE6+ZmZmTebka2Zm1mROvmZmZk3m5Gv2MiPpmTXs/1JJr8n/z5O0WeG3UZKuzP+Pk/TP/Paf+yQdW+juM5JqPplL0jGSPpr/f6ukv+Q35vxK0qBcPkzSc1r5hqcf5vL1Jf1G0l2SJhSGeXbx7TxVxvk+SXMk3SPpVknfyeVTJB2X/z9V0ujezTmzlZx8zaxukt4MDIiIv9XZy8URMRLYA5gsqfLUpHOA/6oxjoGke0ln5KJppDfv7AhcDkwqdP5QRIzMn6Nz2XuA64G3AB/Jw3xrjvuWGuPcgXRv6ocjYnvSPatzq3T6v6S3/5itESdfs5ep/PSiU3IN8U5Jh+XydSSdlWurMyX9WlLlSVhH0IsnF0XEQlIy2zp//xcwT9Lbq3Q+GrglVr704fXAH/P/M4Ex3YxuOfAKYF3SwyIgPQf7y1308z/A1Ii4L8f3YkT8oMp0/B3YVFK1Z2+b1c3J1+zl6wPASOCtwD7AKfnxlR8AhpHeEfwRYPdCP3sAN5eGc22l6ZdUS12NpO1ITwO7o1A8B9izSuflcdzNytcKHsqqzxwenpuI/yCpMqyZOf4bgO9Lej8pmT9aLbZshyrTVcstOUazXvOLFcxevt4BXBgRL5IesP8HYNdc/rOIeAl4XNK1hX62Bv5ZGs7eEfEkpHO+wHGF3w6TtBfwRuAzEfHvwm8LcnnZ1sC9he9HkZLol0mP1Hw+lz8GbBcRCyXtDFwh6c0RsRQYm+NZl/Q84gMlnQZsR3pP7Zo8L3oB8Ko16N/MNV8z65HnSDXYel0cEW8B/gP4Zqm5doM8vC7HERH3RcS+EbEz6fnED+XyZbk5m4i4OZe/vjSsCaQXu+9GetfyYcDnq4zzbmDnOqepVtxmdXPyNXv5mk2qmQ6QtDmwF+lFC38CxuRzv1uy6ivk7gVG9HREETGH9FKIzxWKX8/KtzIVrTIOpffmImkd4EtA5armzSUNyP+/hvQmp78V+hsM7E9Kvq8AXiK9B3nDKuM8BfiipNdXxiXp6CrddRW3Wd2cfM1evi4nnYO9HZgF/E9EPA5cRnrrzT3A+aRznE/lfq5i1WTcE98Cjqy8ko903nRmle6uJh0IVHxI0gPAfcCjwLm5fC/gjnyu+VLg6IhYVOjvK6SLqF4iNT3vSXpDz0/LI4yIO4BjgAsl3UtKrq8pd5ebsUeQzleb9ZrfamRmq5G0UUQ8I2lTUm14j4h4XNKGpPcA75HPFfd2+DsB/x0RH6nx++Wkg4EHezuORpB0MPC2iOjqymmzbvmCKzOr5kpJHcB6wMm5RkxEPCfpRGAb4B9rMPzN6PrWn+NJF161VfIl7TO/0+ogbO3nmq+ZmVmT+ZyvmZlZkzn5mpmZNZmTr5mZWZM5+ZqZmTWZk6+ZmVmT/X8zyd/tXfAgFQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 360x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(5, 4)\n",
    "                        )\n",
    "cph_urti_no_abs_incdt.plot()\n",
    "plt.suptitle('Coefficients of CPH model for incident hospital admission \\n related to URTI with no antibiotics', fontsize=16)\n",
    "#_ = ax.set_yticklabels(['Age', 'BMI', 'GP consultations', 'CCI', 'IMD', 'Flu vaccine', 'Antibacterial prescriptions'], \n",
    "#                        rotation=45, fontsize=12, ha='right')\n",
    "plt.savefig('../output/hospitalisation_prediction_urti/coefficients_cph_urti_no_abs_incdt.jpg', bbox_inches='tight')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # concordance index of cph model with testing data\n",
    "# concord_index_urti_no_abs_incdt_test = concordance_index(data_urti_no_abs_incdt_test['event_urti_admitted'], -cph_urti_no_abs_incdt.predict_partial_hazard(data_urti_no_abs_incdt_test), data_urti_no_abs_incdt_test['duration_urti_admitted'])\n",
    "# print(f'Concordance index for urti_no_abs_incdt with testing data: {concord_index_urti_no_abs_incdt_test}')\n",
    "# # print(concordance_index(data_urti_no_abs_incdt_test['event_urti_admitted'], -cph_urti_no_abs_incdt.predict_partial_hazard(data_urti_no_abs_incdt_test), data_urti_no_abs_incdt_test['duration_urti_admitted']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a = [concord_index_urti_no_abs_incdt_test]\n",
    "# np.savetxt('../output/hospitalisation_prediction_urti/concord_index_urti_no_abs_incdt_test.txt', a, delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig, ax = plt.subplots(figsize=(8, 5)\n",
    "#                         )\n",
    "# survival_probability_calibration(cph_urti_no_abs_incdt, data_urti_no_abs_incdt_train, t0=30)\n",
    "# # plt.suptitle('Calibration of CPH model for incident hospital admission related to urti with no antibiotics', fontsize=16)\n",
    "# plt.savefig('../output/hospitalisation_prediction_urti/calibration_cph_urti_no_abs_incdt_train.jpg', bbox_inches='tight')\n",
    "# plt.show()\n",
    "\n",
    "# # t0 (float) – the time to evaluate the probability of event occurring prior at.\n",
    "# # ICI – mean absolute difference between predicted and observed\n",
    "# # E50 – median absolute difference between predicted and observed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig, ax = plt.subplots(figsize=(8, 5)\n",
    "#                         )\n",
    "# survival_probability_calibration(cph_urti_no_abs_incdt, data_urti_no_abs_incdt_test, t0=30)\n",
    "# # plt.suptitle('Calibration of CPH model for incident hospital admission related to URTI with no antibiotics', fontsize=16)\n",
    "# plt.savefig('../output/hospitalisation_prediction_urti/calibration_cph_urti_no_abs_incdt_test.jpg', bbox_inches='tight')\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cph_urti_no_abs_incdt_predicted = cph_urti_no_abs_incdt.predict_partial_hazard(data_urti_no_abs_incdt)\n",
    "\n",
    "# fig, ax = plt.subplots(figsize=(8, 5))\n",
    "# fig = sm.qqplot(cph_urti_no_abs_incdt_predicted, line='45', ax=ax)\n",
    "# plt.suptitle('Q-Q plot for calibration of CPH model for incident hospital \\n admission related to URTI with no antibiotics', fontsize=16)\n",
    "# plt.savefig('../output/hospitalisation_prediction_urti/qq_calibration_cph_urti_no_abs_incdt.jpg', bbox_inches='tight')\n",
    "# # pylab.show()\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print(cph_urti_no_abs_incdt.score(data_urti_no_abs_incdt_1_test))\n",
    "# print(cph_urti_no_abs_incdt.concordance_index_)\n",
    "# print(cph_urti_no_abs_incdt.AIC_partial_)\n",
    "# print(cph_urti_no_abs_incdt.concordance_index_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4-1-2- Incident hospital admission related to URTI with antibiotics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data_urti_abs_incdt.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data_urti_abs_incdt.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #keep relevant variables\n",
    "# data_urti_abs_incdt = data_urti_abs_incdt[['age', 'bmi_imp', 'gp_count_imp', 'flu_vaccine', 'sex', 'region', #'season', \n",
    "#                                          'imd_imp', 'antibacterial_brit_imp', 'ethnicity_imp', 'smoking_status_imp',\n",
    "#                                          'meds_nonabs',\n",
    "#                                          'event_urti_admitted', 'duration_urti_admitted']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # printing number of missing values in each variable \n",
    "# data_urti_abs_incdt.isnull().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #drop ros with missing values\n",
    "# data_urti_abs_incdt = data_urti_abs_incdt.dropna(axis='rows')\n",
    "# #data_urti_abs_incdt = data_urti_abs_incdt.drop(columns=['smoking_status'])\n",
    "# data_urti_abs_incdt.isnull().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fold_thre = len(data_urti_abs_incdt)-round(len(data_urti_abs_incdt)/4)\n",
    "# data_urti_abs_incdt_train = data_urti_abs_incdt.iloc[:fold_thre]\n",
    "# data_urti_abs_incdt_test = data_urti_abs_incdt.iloc[fold_thre:]\n",
    "\n",
    "# print(f'shape of data_urti_abs_incdt for training: {data_urti_abs_incdt_train.shape}')\n",
    "# print(f'shape of data_urti_abs_incdt for testing: {data_urti_abs_incdt_test.shape}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cph_urti_abs_incdt = CoxPHFitter().fit(data_urti_abs_incdt_train,\n",
    "#         duration_col='duration_urti_admitted',\n",
    "#         event_col='event_urti_admitted',\n",
    "#         formula=\"age + sex + bmi_imp + gp_count_imp + flu_vaccine + region + antibacterial_brit_imp + imd_imp + ethnicity_imp + smoking_status_imp + meds_nonabs\" # season +\n",
    "#         )\n",
    "# cph_urti_abs_incdt.print_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8608\n"
     ]
    }
   ],
   "source": [
    "# output = GetPrintSummary(cph_urti_abs_incdt)\n",
    "\n",
    "# with open('../output/hospitalisation_prediction_urti/cph_urti_abs_incdt_summary.html', \"w\") as fo:\n",
    "#     print(fo.write(output))\n",
    "#     fo.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a = cph_urti_abs_incdt.baseline_hazard_\n",
    "# # cph_urti_abs_incdt.baseline_survival_\n",
    "# # cph_urti_abs_incdt.predict_cumulative_hazard(data_urti_abs_incdt, data_urti_abs_incdt['duration_urti_admitted']<=30)\n",
    "# # cph_urti_abs_incdt.predict_survival_function(data_urti_abs_incdt)\n",
    "\n",
    "# # save \n",
    "# np.savetxt('../output/hospitalisation_prediction_urti/sum_betas_urti_abs_incdt.txt', a, delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig, ax = plt.subplots(figsize=(5, 9)\n",
    "#                         )\n",
    "# cph_urti_abs_incdt.plot()\n",
    "# plt.suptitle('Coefficients of CPH model for incident hospital admission \\n related to URTI with antibiotics', fontsize=16)\n",
    "# #_ = ax.set_yticklabels(['Age', 'BMI', 'GP consultations', 'CCI', 'IMD', 'Flu vaccine', 'Antibacterial prescriptions'], \n",
    "# #                        rotation=45, fontsize=12, ha='right')\n",
    "# plt.savefig('../output/hospitalisation_prediction_urti/coefficients_cph_urti_abs_incdt.jpg', bbox_inches='tight')\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # concordance index of cph model with testing data\n",
    "# concord_index_urti_abs_incdt_test = concordance_index(data_urti_abs_incdt_test['event_urti_admitted'], -cph_urti_abs_incdt.predict_partial_hazard(data_urti_abs_incdt_test), data_urti_abs_incdt_test['duration_urti_admitted'])\n",
    "# print(f'Concordance index for urti_abs_incdt with testing data: {concord_index_urti_abs_incdt_test}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a = [concord_index_urti_abs_incdt_test]\n",
    "# np.savetxt('../output/hospitalisation_prediction_urti/concord_index_urti_abs_incdt_test.txt', a, delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig, ax = plt.subplots(figsize=(8, 5)\n",
    "#                         )\n",
    "# survival_probability_calibration(cph_urti_abs_incdt, data_urti_abs_incdt_train, t0=30)\n",
    "# # plt.suptitle('Calibration of CPH model for incident hospital admission related to URTI with antibiotics', fontsize=16)\n",
    "# plt.savefig('../output/hospitalisation_prediction_urti/calibration_cph_urti_abs_incdt_train.jpg', bbox_inches='tight')\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig, ax = plt.subplots(figsize=(8, 5)\n",
    "#                         )\n",
    "# survival_probability_calibration(cph_urti_abs_incdt, data_urti_abs_incdt_test, t0=30)\n",
    "# # plt.suptitle('Calibration of CPH model for incident hospital admission related to urti with antibiotics', fontsize=16)\n",
    "# plt.savefig('../output/hospitalisation_prediction_urti/calibration_cph_urti_abs_incdt_test.jpg', bbox_inches='tight')\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cph_urti_abs_incdt_predicted = cph_urti_abs_incdt.predict_partial_hazard(data_urti_abs_incdt)\n",
    "\n",
    "# fig, ax = plt.subplots(figsize=(8, 5))\n",
    "# fig = sm.qqplot(cph_urti_abs_incdt_predicted, line='45', ax=ax)\n",
    "# plt.suptitle('Q-Q plot for calibration of CPH model for incident hospital \\n admission related to URTI with antibiotics', fontsize=16)\n",
    "# plt.savefig('../output/hospitalisation_prediction_urti/qq_calibration_cph_urti_abs_incdt.jpg', bbox_inches='tight')\n",
    "# # pylab.show()\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4-1-3- Prevalent hospital admission related to URTI with no antibiotics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data_urti_no_abs_prevl.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data_urti_no_abs_prevl.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #keep relevant variables\n",
    "# data_urti_no_abs_prevl = data_urti_no_abs_prevl[['age', 'bmi_imp', 'gp_count_imp', 'flu_vaccine', 'sex', 'region', #'season', \n",
    "#                                                 'imd_imp', 'antibacterial_brit_imp', 'ethnicity_imp', 'smoking_status_imp',\n",
    "#                                                 'meds_nonabs',\n",
    "#                                                 'event_urti_admitted', 'duration_urti_admitted']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # printing number of missing values in each variable \n",
    "# data_urti_no_abs_prevl.isnull().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #drop rows with missing values\n",
    "# data_urti_no_abs_prevl = data_urti_no_abs_prevl.dropna(axis='rows')\n",
    "# #data_urti_no_abs_prevl = data_urti_no_abs_prevl.drop(columns=['smoking_status'])\n",
    "# data_urti_no_abs_prevl.isnull().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data_urti_no_abs_prevl.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fold_thre = len(data_urti_no_abs_prevl)-round(len(data_urti_no_abs_prevl)/4)\n",
    "# data_urti_no_abs_prevl_train = data_urti_no_abs_prevl.iloc[:fold_thre]\n",
    "# data_urti_no_abs_prevl_test = data_urti_no_abs_prevl.iloc[fold_thre:]\n",
    "\n",
    "# print(f'shape of data_urti_no_abs_prevl for training: {data_urti_no_abs_prevl_train.shape}')\n",
    "# print(f'shape of data_urti_no_abs_prevl for testing: {data_urti_no_abs_prevl_test.shape}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cph_urti_no_abs_prevl = CoxPHFitter().fit(data_urti_no_abs_prevl_train,\n",
    "#         duration_col='duration_urti_admitted',\n",
    "#         event_col='event_urti_admitted',\n",
    "#         formula=\"age + sex + bmi_imp + gp_count_imp + flu_vaccine + region + antibacterial_brit_imp + imd_imp + ethnicity_imp + smoking_status_imp + meds_nonabs\" # season +\n",
    "#         )\n",
    "# cph_urti_no_abs_prevl.print_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# output = GetPrintSummary(cph_urti_no_abs_prevl)\n",
    "\n",
    "# with open('../output/hospitalisation_prediction_urti/cph_urti_no_abs_prevl_summary.html', \"w\") as fo:\n",
    "#     print(fo.write(output))\n",
    "#     fo.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a = cph_urti_no_abs_prevl.baseline_hazard_\n",
    "# # cph_urti_no_abs_prevl.baseline_survival_\n",
    "# # cph_urti_no_abs_prevl.predict_cumulative_hazard(data_urti_no_abs_prevl, data_urti_no_abs_prevl['duration_urti_admitted']<=30)\n",
    "# # cph_urti_no_abs_prevl.predict_survival_function(data_urti_no_abs_prevl)\n",
    "\n",
    "# # save \n",
    "# np.savetxt('../output/hospitalisation_prediction_urti/sum_betas_urti_no_abs_prevl.txt', a, delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig, ax = plt.subplots(figsize=(5, 9)\n",
    "#                         )\n",
    "# cph_urti_no_abs_prevl.plot()\n",
    "# plt.suptitle('Coefficients of CPH model for prevalent hospital admission \\n related to URTI with no antibiotics', fontsize=16)\n",
    "# #_ = ax.set_yticklabels(['Age', 'BMI', 'GP consultations', 'CCI', 'IMD', 'Flu vaccine', 'Antibacterial prescriptions'], \n",
    "# #                        rotation=45, fontsize=12, ha='right')\n",
    "# plt.savefig('../output/hospitalisation_prediction_urti/coefficients_cph_urti_no_abs_prevl.jpg', bbox_inches='tight')\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data_urti_no_abs_prevl_test.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # concordance index of cph model with testing data\n",
    "# concord_index_urti_no_abs_prevl_test = concordance_index(data_urti_no_abs_prevl_test['event_urti_admitted'], -cph_urti_no_abs_prevl.predict_partial_hazard(data_urti_no_abs_prevl_test), data_urti_no_abs_prevl_test['duration_urti_admitted'])\n",
    "# print(f'Concordance index for cph_urti_no_abs_prevl with testing data: {concord_index_urti_no_abs_prevl_test}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a = [concord_index_urti_no_abs_prevl_test]\n",
    "# np.savetxt('../output/hospitalisation_prediction_urti/concord_index_urti_no_abs_prevl_test.txt', a, delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig, ax = plt.subplots(figsize=(8, 5)\n",
    "#                         )\n",
    "# survival_probability_calibration(cph_urti_no_abs_prevl, data_urti_no_abs_prevl_train, t0=30)\n",
    "# plt.suptitle('Calibration of CPH model for prevalent hospital admission related to URTI with no antibiotics', fontsize=16)\n",
    "# plt.savefig('../output/hospitalisation_prediction_urti/calibration_cph_urti_no_abs_prevl_train.jpg', bbox_inches='tight')\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig, ax = plt.subplots(figsize=(8, 5)\n",
    "#                         )\n",
    "# survival_probability_calibration(cph_urti_no_abs_prevl, data_urti_no_abs_prevl_test, t0=30)\n",
    "# plt.suptitle('Calibration of CPH model for prevalent hospital admission related to URTI with no antibiotics', fontsize=16)\n",
    "# plt.savefig('../output/hospitalisation_prediction_urti/calibration_cph_urti_no_abs_prevl_test.jpg', bbox_inches='tight')\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cph_urti_no_abs_prevl_predicted = cph_urti_no_abs_prevl.predict_partial_hazard(data_urti_no_abs_prevl)\n",
    "\n",
    "# fig, ax = plt.subplots(figsize=(8, 5))\n",
    "# fig = sm.qqplot(cph_urti_no_abs_prevl_predicted, line='45', ax=ax)\n",
    "# plt.suptitle('Q-Q plot for calibration of CPH model for prevalent hospital \\n admission related to URTI with no antibiotics', fontsize=16)\n",
    "# plt.savefig('../output/hospitalisation_prediction_urti/qq_calibration_cph_urti_no_abs_prevl.jpg', bbox_inches='tight')\n",
    "# # pylab.show()\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4-1-4- Prevalent hospital admission related to URTI with antibiotics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 385,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data_urti_abs_prevl.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 386,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #keep relevant variables\n",
    "# data_urti_abs_prevl = data_urti_abs_prevl[['age', 'age_cat', 'sex', 'bmi', 'BMI_cat', 'gp_count', 'flu_vaccine',\n",
    "#                                         'region', 'antibacterial_brit', 'imd', 'ethnicity', 'CCI_cat', 'CCI',\n",
    "#                                         'event_urti_admitted', 'duration_urti_admitted']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 387,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # print number of missing values in each variable \n",
    "# data_urti_abs_prevl.isnull().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 388,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #drop ros with missing values\n",
    "# data_urti_abs_prevl = data_urti_abs_prevl.dropna(axis='rows')\n",
    "# #data_urti_abs_prevl = data_urti_abs_prevl.drop(columns=['smoking_status'])\n",
    "# data_urti_abs_prevl.isnull().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 389,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cph_urti_abs_prevl = CoxPHFitter(penalizer=0.1)\n",
    "# cph_urti_abs_prevl.fit(data_urti_abs_prevl, \n",
    "#         duration_col='duration_urti_admitted',\n",
    "#         event_col='event_urti_admitted',\n",
    "#         formula=\"age_cat + sex + gp_count + flu_vaccine + region + antibacterial_brit + imd + ethnicity\"\n",
    "#         )\n",
    "# cph_urti_abs_prevl.print_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 390,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig, ax = plt.subplots(figsize=(5, 9)\n",
    "#                         )\n",
    "# cph_urti_abs_prevl.plot()\n",
    "# plt.suptitle('Coefficients of CPH model for prevalent hospital admission related to URTI with antibiotics', fontsize=16)\n",
    "# #_ = ax.set_yticklabels(['Age', 'BMI', 'GP consultations', 'CCI', 'IMD', 'Flu vaccine', 'Antibacterial prescriptions'], \n",
    "# #                        rotation=45, fontsize=12, ha='right')\n",
    "# #plt.savefig('output/repeat_ABs_outputs/coefficients_cph_urti_abs_prevl.jpg', bbox_inches='tight')\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [],
   "source": [
    "##save cph model with pickle package\n",
    "\n",
    "# current_working_dir = os.getcwd()\n",
    "# print(current_working_dir)\n",
    "\n",
    "# with open('output/hospitalisation_prediction_urti/cph_urti_abs_prevl.pickle', 'wb') as f:\n",
    "#     pickle.dump(cph_urti_abs_prevl, f)\n",
    "\n",
    "# with open('output/hospitalisation_prediction_urti/cph_urti_abs_prevl.pickle', 'rb') as f:\n",
    "#     cph_urti_abs_prevl = pickle.load(f)\n",
    "\n",
    "# cph_urti_abs_prevl.summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {},
   "outputs": [],
   "source": [
    "# hosp history, ab type history, "
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "7a2d6c6560a78692f44e6c684bf0cb5541339d681bcf02c80a9e36dbc5e3249d"
  },
  "kernelspec": {
   "display_name": "Python 3.9.6 64-bit",
   "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.10.4"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}