Predicting_Hospitalization_LogReg.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Predicting Utilization with Logistic Regression\n",
    "In this notebook, we will run queries against the data warehouse to create a dataset.  We will then use this dataset to build a model to predict the risk of hospitalization (admission or ER visit) in the next 6 months."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use logistic regression and its variants as our main modeling tool in this notebook.  Through this, we will see some of the strengths and weaknesses of logistic regression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Standard Imports and Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import sqlalchemy as sa\n",
    "import datetime\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "#import xgboost as xgb\n",
    "#import ml_insights as mli\n",
    "\n",
    "pd.set_option(\"display.max_rows\",9999)\n",
    "pd.set_option(\"display.max_columns\",9999)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run queries to get data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cloverdb\n",
    "\n",
    "conn_str = cloverdb.connection_str_for_tunnel('stg-dwh-db')\n",
    "#conn_str = cloverdb.connection_str_for_tunnel('prod-dwh-db')\n",
    "engine = sa.create_engine(conn_str)\n",
    "engine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query_labs = \"\"\"\n",
    "\n",
    "with cmb_meas_persons as\n",
    "(\n",
    "select * from trg_analytics.fact_lab_observations\n",
    "where (('80053'=ANY(test_codes)) or\n",
    "        ('80048'=ANY(test_codes)) or\n",
    "        ('1188'=ANY(test_codes)) or\n",
    "        ('85027'=ANY(test_codes)) or\n",
    "         ('82948'=ANY(test_codes)) or\n",
    "       ('85025'=ANY(test_codes)) ) \n",
    "\n",
    "and test_date between '2017-01-01' and '2017-06-30'\n",
    ")\n",
    "\n",
    "select distinct on (personid, result_description) *\n",
    "from cmb_meas_persons\n",
    "where result_description is not null\n",
    "order by personid, result_description, test_date desc\n",
    " \n",
    "\"\"\"\n",
    "\n",
    "print(query_labs)\n",
    "\n",
    "with engine.begin() as conn:\n",
    "    df_labs = pd.read_sql(sa.text(query_labs), conn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query_hosps = \"\"\"\n",
    "select personid, sum(case when (contact_type = 'hospital'\n",
    "and contact_subtype = 'admission') then 1 else 0 end) as hosp_visits,\n",
    "sum(case when (contact_type = 'hospital'\n",
    "and contact_subtype = 'emergency_department') then 1 else 0 end) as er_visits,\n",
    "sum(case when (contact_type = 'hospital'\n",
    "and contact_subtype = 'observation') then 1 else 0 end) as obs_visits\n",
    "from trg_analytics.fact_contact_events \n",
    "where low_service_date between '2017-07-01' and '2017-12-31'\n",
    "group by personid\n",
    "\"\"\"\n",
    "\n",
    "print(query_hosps)\n",
    "\n",
    "with engine.begin() as conn:\n",
    "    df_hosps = pd.read_sql(sa.text(query_hosps), conn)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query_hosps_pred = \"\"\"\n",
    "\n",
    "\n",
    "\n",
    "select personid, sum(case when (contact_type = 'hospital'\n",
    "and contact_subtype = 'admission') then 1 else 0 end) as hosp_visits_prev,\n",
    "sum(case when (contact_type = 'hospital'\n",
    "and contact_subtype = 'emergency_department') then 1 else 0 end) as er_visits_prev,\n",
    "sum(case when (contact_type = 'hospital'\n",
    "and contact_subtype = 'observation') then 1 else 0 end) as obs_visits_prev\n",
    "from trg_analytics.fact_contact_events \n",
    "where low_service_date between '2017-01-01' and '2017-06-30'\n",
    "group by personid\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "print(query_hosps_pred)\n",
    "\n",
    "with engine.begin() as conn:\n",
    "    df_hosps_pred = pd.read_sql(sa.text(query_hosps_pred), conn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Process and Clean Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_labs=df_labs[df_labs.result_value!='ALT. TEST PERFORMED']\n",
    "df_labs=df_labs[df_labs.result_value!='SEE NOTE']\n",
    "df_labs=df_labs[df_labs.result_value!='NOT MEASURED']\n",
    "df_labs=df_labs[df_labs.result_value!='false']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_labs=df_labs.loc[df_labs.result_description!='glomerular_filtration_rate_1_73_sq_m_predicted_among_blacks_by_creatinine_based_formula_(mdrd)',:]\n",
    "df_labs=df_labs.loc[df_labs.result_description!='glomerular_filtration_rate_1_73_sq_m_predicted_by_creatinine_based_formula_(mdrd)',:]\n",
    "df_labs=df_labs.loc[df_labs.result_description!='glomerular_filtration_rate_1_73_sq_m_predicted_by_creatinine_based_formula_(ckd_epi)',:]\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_labs.result_value=pd.to_numeric(df_labs.result_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Process result_description to replace troublesome characters with '_'\n",
    "\n",
    "temp = df_labs.result_description\n",
    "df_labs.result_description = list(map(lambda x: x.replace('/','_').replace('[','_').replace(']','_').replace('<','_').replace('>','_')\n",
    "                           .replace('.','_').replace(' ','_').replace('-','_').replace(',','_').lower(),temp))\n",
    "\n",
    "# Create a short version of the name \n",
    "# This assumes, perhaps not safely, that if lab names have the same first x charaters, they are the same lab\n",
    "df_labs['result_desc_short'] = np.array(list(map(lambda x: x[:15].lower(), df_labs.result_description)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_labs.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Make a pivot table to get lab values into columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_lab = df_labs.pivot_table(index='personid', values='result_value', columns = 'result_description').reset_index()\n",
    "df_lab.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_lab.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_lab_short = df_labs.pivot_table(index='personid', values='result_value', columns = 'result_desc_short').reset_index()\n",
    "df_lab_short.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_hosps.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Do some sanity checks on the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(np.isnan(df_hosps))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "np.sum(np.isnan(df_hosps_pred))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_full = pd.merge(df_lab_short, df_hosps_pred, how='left', on='personid')\n",
    "#df_full.rename(columns = {'hosp_visits':'hosp_visits_prev'}, inplace = True)\n",
    "\n",
    "df_full = pd.merge(df_full, df_hosps, how='left', on='personid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(np.isnan(df_full))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Fill in zeros for missing values in these columns\n",
    "\n",
    "for col in ['hosp_visits_prev','er_visits_prev','obs_visits_prev','hosp_visits','er_visits','obs_visits']:\n",
    "    df_full[col] = np.nan_to_num(df_full[col])\n",
    "np.sum(np.isnan(df_full))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_full['is_utilizer'] = ((df_full.hosp_visits>0) | (df_full.er_visits>0)).astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Define some sets of variables\n",
    "\n",
    "selected_predictors_labs = [\n",
    "#'hosp_visits_prev',\n",
    "#'er_visits_prev',\n",
    "#'obs_visits_prev',\n",
    "'glucose__mass_v',\n",
    "'glomerular_filt',\n",
    "'erythrocytes__#',  \n",
    "'platelets__#_vo',\n",
    "'leukocytes__#_v', \n",
    "'hematocrit__vol', \n",
    "'urea_nitrogen_c', \n",
    "'hemoglobin__mas']\n",
    "\n",
    "selected_predictors_ut = [\n",
    "'hosp_visits_prev',\n",
    "'er_visits_prev',\n",
    "'obs_visits_prev',\n",
    "]\n",
    "\n",
    "selected_predictors_utlabs = [\n",
    "'hosp_visits_prev',\n",
    "'er_visits_prev',\n",
    "'obs_visits_prev',\n",
    "'glucose__mass_v',\n",
    "'glomerular_filt',\n",
    "'erythrocytes__#',  \n",
    "'platelets__#_vo',\n",
    "'leukocytes__#_v', \n",
    "'hematocrit__vol', \n",
    "'urea_nitrogen_c', \n",
    "'hemoglobin__mas']\n",
    "\n",
    "\n",
    "\n",
    "selected_predictors_big = [\n",
    "'er_visits_prev',\n",
    "'hosp_visits_prev',\n",
    "'obs_visits_prev',\n",
    "'erythrocyte_dis',\n",
    "'glomerular_filt',\n",
    "'erythrocyte_mea', \n",
    "'erythrocytes__#', \n",
    "'lymphocytes_100',\n",
    "'glucose__mass_v',\n",
    "'leukocytes__#_v',\n",
    "'eosinophils_100', \n",
    "'platelets__#_vo',\n",
    "'alkaline_phosph',\n",
    "'neutrophils_100',\n",
    "'urea_nitrogen_c',\n",
    "'urea_nitrogen__',\n",
    "'neutrophils__#_',\n",
    "'platelet_mean_v',\n",
    "'hemoglobin__mas',\n",
    "'hematocrit__vol',\n",
    "'glucose_random',\n",
    "'basophils_100_l',\n",
    "'monocytes_100_l',\n",
    "'calcium__mass_v',\n",
    "'albumin__mass_v',\n",
    "'lymphocytes__#_',\n",
    "'alanine_aminotr',\n",
    "'aspartate_amino',\n",
    "'potassium__mole',\n",
    "'neutrophils_ban',\n",
    "'sodium__moles_v',\n",
    "'creatinine__mas',\n",
    "'monocytes__#_vo',\n",
    "'globulin__mass_',\n",
    "'bicarbonate__mo',\n",
    "'bilirubin_total',\n",
    "'eosinophils__#_',\n",
    "'hematocrit',\n",
    "'protein__mass_v',\n",
    "'mononuclear_cel',\n",
    "'chloride__moles',\n",
    "'albumin_globuli']\n",
    "\n",
    "selected_predictors = selected_predictors_ut"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define some useful functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_reliability_diagram(y,x,bins=np.linspace(0,1,21),size_points=True, show_baseline=True,\n",
    "                                error_bars=True, ax=None, marker='+',c='red', **kwargs):\n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "        fig = ax.get_figure()\n",
    "    digitized_x = np.digitize(x, bins)\n",
    "    mean_count_array = np.array([[np.mean(y[digitized_x == i]),len(y[digitized_x == i]),np.mean(x[digitized_x==i])] for i in np.unique(digitized_x)])\n",
    "    x_pts_to_graph = mean_count_array[:,2]\n",
    "    y_pts_to_graph = mean_count_array[:,0]\n",
    "    pt_sizes = mean_count_array[:,1]\n",
    "    if show_baseline:\n",
    "        ax.plot(np.linspace(0,1,100),(np.linspace(0,1,100)),'k--')\n",
    "    for i in range(len(y_pts_to_graph)):\n",
    "        if size_points:\n",
    "            plt.scatter(x_pts_to_graph,y_pts_to_graph,s=pt_sizes,marker=marker,c=c, **kwargs)\n",
    "        else:\n",
    "            plt.scatter(x_pts_to_graph,y_pts_to_graph, **kwargs)\n",
    "    #plt.axis([-0.1,1.1,-0.1,1.1])\n",
    "    if error_bars:\n",
    "        plt.errorbar(x_pts_to_graph, x_pts_to_graph,\n",
    "                    yerr=1.96*x_pts_to_graph*(1-x_pts_to_graph)/(np.sqrt(pt_sizes)), capsize=5)\n",
    "\n",
    "    return(x_pts_to_graph,y_pts_to_graph,pt_sizes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def histogram_pair(value_vec, binary_vec, **kwargs):\n",
    "    plt.figure(figsize = (10,6))\n",
    "    plt.subplot(3,1,1)\n",
    "    out1 = plt.hist(value_vec[binary_vec==1], **kwargs)\n",
    "    plt.subplot(3,1,2)\n",
    "    out0 = plt.hist(value_vec[binary_vec==0], **kwargs)\n",
    "    bin_leftpts = (out1[1])[:-1]\n",
    "    bin_rightpts = (out1[1])[1:]\n",
    "    bin_centers = (bin_leftpts+bin_rightpts)/2\n",
    "    prob_numer = out1[0]\n",
    "    prob_denom = out1[0]+out0[0]\n",
    "    smoothing_const = .001\n",
    "    probs = (prob_numer+smoothing_const)/(prob_denom+2*smoothing_const)\n",
    "    plt.subplot(3,1,3)\n",
    "    plt.plot(bin_centers, prob_numer/prob_denom)\n",
    "    plt.errorbar(bin_centers,probs, yerr=1.96*probs*(1-probs)/np.sqrt(prob_denom),capsize=3 )\n",
    "    plt.xlim(bin_leftpts[0], bin_rightpts[-1])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_test_evaluate_model(model, X_tr, y_tr, X_te, y_te):\n",
    "    model.fit(X_tr,y_tr)\n",
    "    test_preds = model.predict_proba(X_te)[:,1]\n",
    "    prec, rec, _ = precision_recall_curve(y_test, test_preds)\n",
    "    fpr, tpr, _ = roc_curve(y_test, test_preds)\n",
    "    \n",
    "    plt.figure(figsize=(20,8))\n",
    "    plt.subplot(1,3,1)\n",
    "    plt.plot(prec[:-1], rec[:-1])\n",
    "    plt.xlim(0, 1) \n",
    "    plt.ylim(0,1)\n",
    "\n",
    "    plt.subplot(1,3,2)\n",
    "    plt.plot(fpr, tpr)\n",
    "    plt.xlim(0, 1) \n",
    "    plt.ylim(0,1)\n",
    "    plt.plot([0,1],[0,1], 'k--')\n",
    "\n",
    "    plt.subplot(1,3,3)\n",
    "    plot_reliability_diagram(y_test,test_preds)\n",
    "    \n",
    "    return roc_auc_score(y_test, test_preds), log_loss(y_test, test_preds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Do a little data exploration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "histogram_pair(df_full['hosp_visits_prev'].values, df_full.is_utilizer.values, bins=10, range = (-.5,9.5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "histogram_pair(df_full['er_visits_prev'].values, df_full.is_utilizer, bins=10, range = (-.5,9.5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import precision_recall_curve, roc_auc_score, log_loss, roc_curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = df_full.loc[:,selected_predictors]\n",
    "y = df_full['is_utilizer']\n",
    "#y = ((df_full.hosp_visits>0) | (df_full.er_visits>0) | (df_full.obs_visits>0)).astype(int)\n",
    "\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X,y, test_size = .25, random_state=42)\n",
    "X_train.shape, X_test.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_train.value_counts(), y_test.value_counts(), y_train.mean(), y_test.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegressionCV, LogisticRegression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logistic Regression \"off-the-shelf\"\n",
    "In sklearn, the \"default\" Logistic Regression is actually an L2-Penalized Logistic Regresssion with C=1 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_model0 = LogisticRegression()  \n",
    "train_test_evaluate_model(lr_model0, X_train, y_train, X_test, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(zip(selected_predictors,lr_model0.coef_[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_model0a = LogisticRegression(C=1000000)  \n",
    "train_test_evaluate_model(lr_model0a, X_train, y_train, X_test, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(zip(selected_predictors,lr_model0a.coef_[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_model0b = LogisticRegression(C=.01)  \n",
    "train_test_evaluate_model(lr_model0b, X_train, y_train, X_test, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(zip(selected_predictors,lr_model0b.coef_[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using Cross-Validation to choose best Nuisance Parameter Value\n",
    "Rather than just picking arbitrary values of \"C\", let's use cross-validation to try and find the best value.  Also, let's use an L1 penalty rather than L2, so that we get a more parsimonious model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import make_scorer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_model1 = LogisticRegressionCV(Cs = 10**np.linspace(-5,5,61), penalty='l1', solver = 'liblinear', \n",
    "                                 scoring = make_scorer(log_loss, greater_is_better=False))  \n",
    "train_test_evaluate_model(lr_model1, X_train, y_train, X_test, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(zip(selected_predictors,lr_model1.coef_[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_model1.Cs_, lr_model1.C_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(lr_model1.scores_[1], axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(np.log10(lr_model1.Cs_),np.mean(lr_model1.scores_[1], axis=0), 'x-')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Standardizing Variables Beforehand\n",
    "The issue with regularizing based on coefficient size is that arbitrary differences in the values of the variables (such as units) can have a big effect on how a variable is penalized.  One solution is to \"standardize\" the variables beforehand (i.e. subtract the mean and divide by the standard deviation of the variable)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "std_scaler = StandardScaler()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "std_scaler.fit(X_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "std_scaler.mean_ ,std_scaler.scale_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_std  = std_scaler.transform(X_train)\n",
    "X_test_std = std_scaler.transform(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_model2 = LogisticRegressionCV(Cs = 10**np.linspace(-5,5,61), penalty='l1', solver = 'liblinear',\n",
    "                                scoring = make_scorer(log_loss, greater_is_better=False))    \n",
    "train_test_evaluate_model(lr_model2, X_train_std, y_train, X_test_std, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(zip(selected_predictors,lr_model2.coef_[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Repeat this process using just labs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_predictors = selected_predictors_labs\n",
    "\n",
    "X = df_full.loc[:,selected_predictors]\n",
    "y = ((df_full.hosp_visits>0) | (df_full.er_visits>0)).astype(int)\n",
    "#y = ((df_full.hosp_visits>0) | (df_full.er_visits>0) | (df_full.obs_visits>0)).astype(int)\n",
    "plot_mat = X.copy()\n",
    "plot_mat[\"had_hosp_visit\"] = y\n",
    "\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X,y, test_size = .25, random_state=42)\n",
    "X_train.shape, X_test.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Need to remove missing values - will impute the median"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import Imputer\n",
    "impute_median = Imputer(strategy = 'median')\n",
    "X_train_med = impute_median.fit_transform(X_train)\n",
    "X_test_med = impute_median.transform(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_med = pd.DataFrame(X_train_med, columns = selected_predictors)\n",
    "X_test_med = pd.DataFrame(X_test_med, columns = selected_predictors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "std_scaler = StandardScaler()\n",
    "std_scaler.fit(X_train_med)\n",
    "X_train_med_std  = std_scaler.transform(X_train_med)\n",
    "X_test_med_std = std_scaler.transform(X_test_med)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_model0 = LogisticRegression()\n",
    "train_test_evaluate_model(lr_model0, X_train_med, y_train, X_test_med, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(zip(selected_predictors,lr_model0.coef_[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_model1 = LogisticRegressionCV(Cs = 10**np.linspace(-5,5,61), penalty='l1', solver = 'liblinear', \n",
    "                                 scoring = make_scorer(log_loss, greater_is_better=False))  \n",
    "train_test_evaluate_model(lr_model1, X_train_med, y_train, X_test_med, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(zip(selected_predictors,lr_model1.coef_[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_model2 = LogisticRegressionCV(Cs = 10**np.linspace(-5,5,61), penalty='l1', solver = 'liblinear', \n",
    "                                 scoring = make_scorer(log_loss, greater_is_better=False))  \n",
    "train_test_evaluate_model(lr_model2, X_train_med_std, y_train, X_test_med_std, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(zip(selected_predictors,lr_model2.coef_[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_test_evaluate_model(lr_model2, X_train_med_std, y_train, X_test_med_std, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(zip(selected_predictors,lr_model2.coef_[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "histogram_pair(X_train_med['erythrocytes__#'].values,y_train.values, bins=12, range = (2,8))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "histogram_pair(X_train_med['hematocrit__vol'].values,y_train.values, bins=24, range=(28,52))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "histogram_pair(X_train_med['platelets__#_vo'].values,y_train.values,bins=np.concatenate(([0,100],np.linspace(150,300,16),[350,450,600])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Repeat with both labs and utilization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_predictors = selected_predictors_utlabs\n",
    "\n",
    "X = df_full.loc[:,selected_predictors]\n",
    "y = ((df_full.hosp_visits>0) | (df_full.er_visits>0)).astype(int)\n",
    "#y = ((df_full.hosp_visits>0) | (df_full.er_visits>0) | (df_full.obs_visits>0)).astype(int)\n",
    "plot_mat = X.copy()\n",
    "plot_mat[\"had_hosp_visit\"] = y\n",
    "\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X,y, test_size = .25, random_state=42)\n",
    "X_train.shape, X_test.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import Imputer\n",
    "impute_median = Imputer(strategy = 'median')\n",
    "X_train_med = impute_median.fit_transform(X_train)\n",
    "X_test_med = impute_median.transform(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "std_scaler = StandardScaler()\n",
    "std_scaler.fit(X_train_med)\n",
    "X_train_med_std  = std_scaler.transform(X_train_med)\n",
    "X_test_med_std = std_scaler.transform(X_test_med)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_model0 = LogisticRegression()\n",
    "train_test_evaluate_model(lr_model0, X_train_med, y_train, X_test_med, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(zip(selected_predictors,lr_model0.coef_[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_model1 = LogisticRegressionCV(Cs = 10**np.linspace(-5,5,61), penalty='l1', solver = 'liblinear', \n",
    "                                 scoring = make_scorer(log_loss, greater_is_better=False))  \n",
    "train_test_evaluate_model(lr_model1, X_train_med, y_train, X_test_med, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(zip(selected_predictors,lr_model1.coef_[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_model2 = LogisticRegressionCV(Cs = 10**np.linspace(-5,5,61), penalty='l1', solver = 'liblinear', \n",
    "                                 scoring = make_scorer(log_loss, greater_is_better=False))  \n",
    "train_test_evaluate_model(lr_model2, X_train_med_std, y_train, X_test_med_std, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(zip(selected_predictors,lr_model2.coef_[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}