<jupyter_start><jupyter_text>Predicting Length-of-Stay at Hospitals (Project Notebook)Author: <NAME>   Email: <EMAIL>   The goal of this project is offer a predicted Length-of-Stay (LOS) for each patient at time of hospital admission. The inputs to the prediction model will include features such as patient’s diagnosis category (e.g. heart disease, delivery, injury/poisoning, etc.), age, gender, and race.  Table of Contents1. [**Project Overview**](Overview)1.1 [Problem Statement](Problem_statement)1.2 [Metrics](Metrics)2. [**Data Exploration and Feature Engineering**](Data-Exploration)   2.1 [ADMISSIONS.csv](ADMISSIONS.csv)   2.1.1 [*Length-of-Stay (Feature Engineering)*](LOS_FE)    2.1.2 [*Deathtime*](Deathtime)  2.1.3 [*Ethnicity*](Ethnicity)  2.1.4 [*Religion*](Religion)   2.1.5 [*Admission Type*](Adm_type)   2.1.6 [*Insurance*](Insurance)  2.1.7 [*Marital Status*](Marital)  2.2 [DIAGNOSES_ICD.csv](DIAGNOSES_ICD.csv)      2.2.1 [ICD-9 Code (Feature Engineering)](ICD_FE)   2.3 [PATIENTS.csv](PATIENTS.csv)     2.3.1 [*Age (Feature Engineering)*](Age_FE)   2.3.2 [*Gender*](Gender)   2.4 [ICUSTAYS.csv](ICUSTAYS.csv)  2.4.1 [*ICU (Feature Engineering)*](ICU_FE)  3. [**Data Preprocessing**](Data-Preprocessing)   3.1 [Data Cleaning Function](Data-cleanup)4. [**Length-of-Stay Prediction Model**](Prediction)   4.1 [Model Implementation](Implementation)  4.2 [Model Refinement](Refinement)      5. [**Results**](Results)6. [**Conclusion**](Conclusion) 1. Project OverviewPredictive analytics is an increasingly important tool in the healthcare field since modern machine learning (ML) methods can use large amounts of available data to predict individual outcomes for patients. For example, ML predictions can help healthcare providers determine likelihoods of disease, aid in diagnosis, recommend treatment, and predict future wellness. For this project, I chose to focus on a more logistical metric of healthcare, hospital length-of-stay (LOS). LOS is defined as the time between hospital admission and discharge measured in days.U.S. hospital stays cost the health system at least $377.5 billion per year and recent Medicare legislation standardizes payments for procedures performed, regardless of the number of days a patient spends in the hospital. This incentivizes hospitals to identify patients of high LOS risk at the time of admission. Once identified, patients with high LOS risk can have their treatment plan optimized to minimize LOS and lower the chance of getting a hospital-acquired condition such as staph infection. Another benefit is that prior knowledge of LOS can aid in logistics such as room and bed allocation planning. 1.1 Problem StatementThe goal of this project is to create a model that predicts the length-of-stay for each patient at time of admission. My first intuition was that the inputs to the prediction model should ideally include features such as patient’s diagnosis category (e.g. heart disease, delivery, injury/poisoning, etc.), age, gender, marital status, and race. After looking for a useful medical database, I ended up choosing the MIT MIMIC-III database due to the robust amount of information it held. MIMIC is described as the following [(source)](https://mimic.physionet.org/). “MIMIC is an openly available dataset developed by the MIT Lab for Computational Physiology, comprising de-identified health data associated with ~40,000 critical care patients. It includes demographics, vital signs, laboratory tests, medications, and more.” MIMIC is a relational database containing tables of data relating to patients who stayed within the intensive care units at Beth Israel Deaconess Medical Center in Boston, MA. As a side note, MIMIC-III requires taking a research ethics and compliance training [course]( https://about.citiprogram.org/en/homepage/) and filling out a research application form for MIMIC on https://physionet.org. Additionally, the database does not include pediatric information (ages 2-13). Due to the restrictive nature of the dataset, I will not be able to provide a copy in the GitHub repo associated with this notebook.  1.2 MetricsThe expected outcome of this project is to develop a model that will be better at predicting hospital LOS than the industry standards of median and average LOS. The median LOS is simply the median LOS of past admissions to a hospital. Similarly, a second commonly used metric in healthcare is the average, or mean LOS. To measure performance, I’ll compare the prediction model against the median and average LOS using the root-mean-square error (RMSE). The RMSE is a commonly used measure of the differences between values predicted by a model and the values observed, where a lower score implies better accuracy. For example, a perfect prediction model would have an RMSE of 0. The RMSE equation for this work is given as follows, where (n) is the number of hospital admission records, ($\hat{y}$) the prediction LOS, and (y) is the actual LOS.The ultimate goal is to develop a prediction model that results in a lower RMSE than the average or median models. There is a multitude of regression models available for predicting LOS. To determine the best regression model for this work (of the subset of models that will be evaluated), the R2 (R-squared) score will be used. The R2 is a measure of the goodness of the fit of a model. In other words, it is the proportion of the variance in the dependent variable that is predictable from the independent variables. R2 is defined as the following equation where ($y_i$) is an observed data point, ($\hat{y}$) is the mean of the observed data, and ($f_i$) the predicted model value.Best possible R2 score is 1.0 and a negative value means it is worse than a constant model, average or median in this case.<jupyter_code># Imports
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score, mean_squared_error
from sklearn.preprocessing import MinMaxScaler
from sklearn.neighbors import KNeighborsRegressor
from sklearn.linear_model import LinearRegression
from sklearn.svm import SVR
from sklearn.ensemble import RandomForestRegressor
from sklearn.tree import DecisionTreeRegressor
from scipy.stats import pearsonr
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.linear_model import SGDRegressor
import statsmodels.api as sm
from sklearn.model_selection import GridSearchCV
import seaborn as sns<jupyter_output><empty_output><jupyter_text>The following csv files were downloaded from the MIMIC-III database [source](https://mimic.physionet.org/).<jupyter_code># Primary Admissions information
df = pd.read_csv('data/ADMISSIONS.csv')

# Patient specific info such as gender
df_pat = pd.read_csv('data/PATIENTS.csv')

# Diagnosis for each admission to hospital
df_diagcode = pd.read_csv('data/DIAGNOSES_ICD.csv')

# Intensive Care Unit (ICU) for each admission to hospital
df_icu = pd.read_csv('data/ICUSTAYS.csv')<jupyter_output><empty_output><jupyter_text>2. Data Exploration and Feature EngineeringIn this section, I'll examine the various imported MIMIC DataFrames to understand how the data is distributed. Additionaly, I need to figure out a strategy to extract the target Length-of-Stay (LOS) values and understand what features (independent variables) may be useful in predicting LOS. **NOTE**: I have to perform a good amount of feature engineering in this section to enable the visualization of certain data categories. Any additional data 'tidyness' problems will be addressed in the 'Data Preprocessing' section including a comprehensive preprocessing function.  2.1 ADMISSIONS.csv Exploration From [MIMIC](https://mimic.physionet.org/mimictables/admissions/): The ADMISSIONS table gives information regarding a patient’s admission to the hospital. Since each unique hospital visit for a patient is assigned a unique HADM_ID, the ADMISSIONS table can be considered as a definition table for HADM_ID. Information available includes timing information for admission and discharge, demographic information, the source of the admission, and so on.<jupyter_code>df.info()
print('Dataset has {} number of unique admission events.'.format(df['HADM_ID'].nunique()))
print('Dataset has {} number of unique patients.'.format(df['SUBJECT_ID'].nunique()))<jupyter_output>Dataset has 58976 number of unique admission events.
Dataset has 46520 number of unique patients.<jupyter_text>2.1.1 Length-of-Stay (Feature Engineering)The first task is to figure out a way to calculate the LOS. LOS is defined as the time between admission and discharge from the hospital.<jupyter_code>df.head()
# Convert admission and discharge times to datatime type
df['ADMITTIME'] = pd.to_datetime(df['ADMITTIME'])
df['DISCHTIME'] = pd.to_datetime(df['DISCHTIME'])

# Convert timedelta type into float 'days', 86400 seconds in a day
df['LOS'] = (df['DISCHTIME'] - df['ADMITTIME']).dt.total_seconds()/86400
# Verify
df[['ADMITTIME', 'DISCHTIME', 'LOS']].head()
df['LOS'].describe()
# Look at what is happening with negative LOS values
df[df['LOS'] < 0]
# Drop rows with negative LOS, usually related to a time of death before admission
df['LOS'][df['LOS'] > 0].describe()
# Drop LOS < 0 
df = df[df['LOS'] > 0]
# Plot LOS Distribution
plt.hist(df['LOS'], bins=200, color = '#55a868')
plt.xlim(0, 50)
plt.title('Distribution of LOS for all hospital admissions \n incl. deceased')
plt.ylabel('Count')
plt.xlabel('Length-of-Stay (days)')
plt.tick_params(top=False, right=False) 
plt.show();
# Pre-emptively drop some columns that I don't need anymore
df.drop(columns=['DISCHTIME', 'ROW_ID', 
                'EDREGTIME', 'EDOUTTIME', 'HOSPITAL_EXPIRE_FLAG',
                'HAS_CHARTEVENTS_DATA'], inplace=True)<jupyter_output><empty_output><jupyter_text>2.1.2 DeathtimePatients who died at the hospital will later be **dropped** from the DataFrame as they would bias the LOS since they never were discharged 'alive' from the hospital. In other words, LOS would be shorter for this group.<jupyter_code># Mark admissions where patients died in boolean column
df['DECEASED'] = df['DEATHTIME'].notnull().map({True:1, False:0})
print("{} of {} patients died in the hospital".format(df['DECEASED'].sum(), 
                                                     df['SUBJECT_ID'].nunique()))
# Look at statistics less admissions resulting in death
df['LOS'].loc[df['DECEASED'] == 0].describe()
# Hospital LOS metrics for later comparison
actual_mean_los = df['LOS'].loc[df['DECEASED'] == 0].mean() 
actual_median_los = df['LOS'].loc[df['DECEASED'] == 0].median() 

print(actual_mean_los)
print(actual_median_los)
plt.hist(df['LOS'].loc[df['DECEASED'] == 0], bins=200, color = '#55a868')
plt.xlim(0, 50)
plt.title('Distribution of LOS for hospital admissions')
plt.ylabel('Count')
plt.xlabel('Length-of-Stay (days)')
plt.tick_params(top=False, right=False) 
plt.show();<jupyter_output><empty_output><jupyter_text>2.1.3 Ethnicity Given the large amount of categories in the 'ETHNICITY' series, my aim was to reduce the column to a few categories. Interestly, asians have the lowest median LOS.<jupyter_code>df['ETHNICITY'].value_counts()
# Compress the number of ethnicity categories
df['ETHNICITY'].replace(regex=r'^ASIAN\D*', value='ASIAN', inplace=True)
df['ETHNICITY'].replace(regex=r'^WHITE\D*', value='WHITE', inplace=True)
df['ETHNICITY'].replace(regex=r'^HISPANIC\D*', value='HISPANIC/LATINO', inplace=True)
df['ETHNICITY'].replace(regex=r'^BLACK\D*', value='BLACK/AFRICAN AMERICAN', inplace=True)
df['ETHNICITY'].replace(['UNABLE TO OBTAIN', 'OTHER', 'PATIENT DECLINED TO ANSWER', 
                         'UNKNOWN/NOT SPECIFIED'], value='OTHER/UNKNOWN', inplace=True)
df['ETHNICITY'].loc[~df['ETHNICITY'].isin(df['ETHNICITY'].value_counts().nlargest(5).index.tolist())] = 'OTHER/UNKNOWN'
df['ETHNICITY'].value_counts()
# Re-usable plotting function
def plot_los_groupby(variable, size=(7,4)):
    ''' 
    Plot Median LOS by df categorical series name 
    '''
    results = df[[variable, 'LOS']].groupby(variable).median().reset_index()
    values = list(results['LOS'].values)
    labels = list(results[variable].values)

    fig, ax = plt.subplots(figsize=size)
    ind = range(len(results))
    ax.barh(ind, values, align='center', height=0.6, color = '#55a868', alpha=0.8)
    ax.set_yticks(ind)
    ax.set_yticklabels(labels)
    ax.set_xlabel('Median Length of Stay (days)')
    ax.tick_params(left=False, top=False, right=False) 
    ax.set_title('Comparison of {} labels'.format(variable))
    
    plt.tight_layout()
    plt.show();
    
# Re-usable boxplot function
def boxplot_los_groupby(variable, los_range=(-1, 30), size=(8,4)):
    ''' 
    Boxplot of LOS by df categorical series name 
    '''
    results = df[[variable, 'LOS']].groupby(variable).median().reset_index()

    categories = results[variable].values.tolist()

    hist_data = []
    for cat in categories:
        hist_data.append(df['LOS'].loc[df[variable]==cat].values)

    fig, ax = plt.subplots(figsize=size)
    ax.boxplot(hist_data, 0, '', vert=False)
    ax.set_xlim(los_range)
    ax.set_yticklabels(categories)
    ax.set_xlabel('Length of Stay (days)')
    ax.tick_params(left=False, right=False) 
    ax.set_title('Comparison of {} categories'.format(variable))
    plt.tight_layout()
    plt.show();

boxplot_los_groupby('ETHNICITY', los_range=(-1, 30))<jupyter_output><empty_output><jupyter_text>2.1.4 ReligionI thought it was interesting that this was tracked in MIMIC. BUDDHIST, HINDU, and MUSLIM trend with lower LOS medians, but are a very small proportion of the dataset.<jupyter_code>df['RELIGION'].value_counts()
boxplot_los_groupby('RELIGION', los_range=(-5, 30), size=(7, 7))
# Reduce categories to terms of religious or not
# I tested with and without category reduction, with little change in R2 score
df['RELIGION'].loc[~df['RELIGION'].isin(['NOT SPECIFIED', 'UNOBTAINABLE'])] = 'RELIGIOUS'

print(df['RELIGION'].value_counts())
print(df['RELIGION'].value_counts()[0]/len(df['RELIGION']))
print(df['RELIGION'].value_counts()[1]/len(df['RELIGION']))
print(df['RELIGION'].value_counts()[2]/len(df['RELIGION']))
# Look at median LOS for groups
plot_los_groupby('RELIGION', size=(5,4))<jupyter_output><empty_output><jupyter_text>2.1.5 Admission TypeNewborns have the lowest median LOS followed by elective admissions. This is expected since these are often somewhat planned for and with the risks being understood.<jupyter_code>df['ADMISSION_TYPE'].value_counts()
boxplot_los_groupby('ADMISSION_TYPE', los_range=(-5, 35), size=(7, 4))<jupyter_output><empty_output><jupyter_text>2.1.6 InsuranceIf a patient is 'Self-Pay', typically means they can't or didn't pay. Self-pay patients have the lowest LOS.<jupyter_code>df['INSURANCE'].value_counts()
boxplot_los_groupby('INSURANCE', los_range=(-5, 30), size=(7, 4))<jupyter_output><empty_output><jupyter_text>2.1.7 Marital status<jupyter_code>df['MARITAL_STATUS'].value_counts(dropna=False)
# Fix NaNs and file under 'UNKNOWN'
df['MARITAL_STATUS'] = df['MARITAL_STATUS'].fillna('UNKNOWN (DEFAULT)')
df['MARITAL_STATUS'].value_counts(dropna=False)
plot_los_groupby('MARITAL_STATUS')
boxplot_los_groupby('MARITAL_STATUS')<jupyter_output><empty_output><jupyter_text>2.2 DIAGNOSES_ICD.csv Exploration This section explore the ICUSTAYS.csv table of the MIMIC-III dataset.  2.2.1 ICD-9 Code Feature EngineeringFrom [Wiki](https://en.wikipedia.org/wiki/International_Statistical_Classification_of_Diseases_and_Related_Health_ProblemsICD-9): International Classification of Diseases, Clinical Modification (ICD-9-CM) is an adaption created by the U.S. National Center for Health Statistics (NCHS) and used in assigning diagnostic and procedure codes associated with inpatient, outpatient, and physician office utilization in the United States.<jupyter_code># Imports
#import pandas as pd
#import numpy as np
#import matplotlib.pyplot as plt
#df_diagcode = pd.read_csv('DIAGNOSES_ICD.csv')
df_diagcode.info()
print('There are {} unique ICD9 codes in this dataset.'.format(df_diagcode['ICD9_CODE'].value_counts().count()))<jupyter_output>There are 6984 unique ICD9 codes in this dataset.<jupyter_text>Because it's not feasible to have 6984 unique values to use as features for predicting LOS, I need to reduce the diagnosis into more general categories. After researching the ICD9 approach, I discovery that they arranged into super categories as the following [(source)](https://en.wikipedia.org/wiki/List_of_ICD-9_codes):* 001–139: infectious and parasitic diseases* 140–239: neoplasms* 240–279: endocrine, nutritional and metabolic diseases, and immunity disorders* 280–289: diseases of the blood and blood-forming organs* 290–319: mental disorders* 320–389: diseases of the nervous system and sense organs* 390–459: diseases of the circulatory system* 460–519: diseases of the respiratory system* 520–579: diseases of the digestive system* 580–629: diseases of the genitourinary system* 630–679: complications of pregnancy, childbirth, and the puerperium* 680–709: diseases of the skin and subcutaneous tissue* 710–739: diseases of the musculoskeletal system and connective tissue* 740–759: congenital anomalies* 760–779: certain conditions originating in the perinatal period* 780–799: symptoms, signs, and ill-defined conditions* 800–999: injury and poisoning* E and V codes: external causes of injury and supplemental classification, *using 999 as placeholder even though it overlaps with complications of medical care*Now the task is to re-code the in terms of these category groups.<jupyter_code># Filter out E and V codes since processing will be done on the numeric first 3 values
df_diagcode['recode'] = df_diagcode['ICD9_CODE']
df_diagcode['recode'] = df_diagcode['recode'][~df_diagcode['recode'].str.contains("[a-zA-Z]").fillna(False)]
df_diagcode['recode'].fillna(value='999', inplace=True)
# https://stackoverflow.com/questions/46168450/replace-specific-range-of-values-in-data-frame-pandas
df_diagcode['recode'] = df_diagcode['recode'].str.slice(start=0, stop=3, step=1)
df_diagcode['recode'] = df_diagcode['recode'].astype(int)
# ICD-9 Main Category ranges
icd9_ranges = [(1, 140), (140, 240), (240, 280), (280, 290), (290, 320), (320, 390), 
               (390, 460), (460, 520), (520, 580), (580, 630), (630, 680), (680, 710),
               (710, 740), (740, 760), (760, 780), (780, 800), (800, 1000), (1000, 2000)]

# Associated category names
diag_dict = {0: 'infectious', 1: 'neoplasms', 2: 'endocrine', 3: 'blood',
             4: 'mental', 5: 'nervous', 6: 'circulatory', 7: 'respiratory',
             8: 'digestive', 9: 'genitourinary', 10: 'pregnancy', 11: 'skin', 
             12: 'muscular', 13: 'congenital', 14: 'prenatal', 15: 'misc',
             16: 'injury', 17: 'misc'}

# Re-code in terms of integer
for num, cat_range in enumerate(icd9_ranges):
    df_diagcode['recode'] = np.where(df_diagcode['recode'].between(cat_range[0],cat_range[1]), 
            num, df_diagcode['recode'])
    
# Convert integer to category name using diag_dict
df_diagcode['recode'] = df_diagcode['recode']
df_diagcode['cat'] = df_diagcode['recode'].replace(diag_dict)
# Verify
df_diagcode.head()<jupyter_output><empty_output><jupyter_text>For each admission, there could be (and usually is) more than one diagnosis. Often, there are more than 1 diagnoses for 1 category. Therefore, I need to create a dummy matrix that highlights all the diagnoses for each admission. This should not be done on the SUBJECT_ID since each patient could have different diagnoses for each admission.<jupyter_code># Create list of diagnoses for each admission
hadm_list = df_diagcode.groupby('HADM_ID')['cat'].apply(list).reset_index()
hadm_list.head()
# Convert diagnoses list into hospital admission-item matrix
hadm_item = pd.get_dummies(hadm_list['cat'].apply(pd.Series).stack()).sum(level=0)
hadm_item.head()
# Join back with HADM_ID, will merge with main admissions DF later
hadm_item = hadm_item.join(hadm_list['HADM_ID'], how="outer")
hadm_item.head()
# Merge with main admissions df
df = df.merge(hadm_item, how='inner', on='HADM_ID')
# Verify Merge
df.head()
# Look at the median LOS by diagnosis category
diag_cat_list = ['skin', 'infectious',  'misc', 'genitourinary', 'neoplasms', 'blood', 'respiratory', 
                  'congenital','nervous', 'muscular', 'digestive', 'mental', 'endocrine', 'injury',
                 'circulatory', 'prenatal',  'pregnancy']

results = []
for variable in diag_cat_list:
    results.append(df[[variable, 'LOS']].groupby(variable).median().reset_index().values[1][1])

sns.set(style="whitegrid")
#sns.set_style("ticks")
fig, ax = plt.subplots(figsize=(7,5))
ind = range(len(results))
ax.barh(ind, results, align='edge', color = '#55a868', alpha=0.8)
ax.set_yticks(ind)
ax.set_yticklabels(diag_cat_list)
ax.set_xlabel('Median Length of Stay (days)')
ax.tick_params(left=False, right=False, top=False) 
ax.set_title('Comparison of Diagnoses'.format(variable))
plt.show();<jupyter_output><empty_output><jupyter_text>2.3 Patients.csv Exploration The PATIENTS table provides age and gender information. To protect identity, a patient’s age is given by the difference between their 'DOB' date of birth and the date of their first admission. Therefore, subsequent admissions for the same patient need to be ignored in the calculation. The only things that need to be done with this table are to extract the DOB and gender information and merge them with the admissions dataframe.<jupyter_code>df_pat.head()
df_pat['GENDER'].value_counts()
# Convert to datetime type
df_pat['DOB'] = pd.to_datetime(df_pat['DOB'])
df_pat = df_pat[['SUBJECT_ID', 'DOB', 'GENDER']]
df_pat.head()
df = df.merge(df_pat, how='inner', on='SUBJECT_ID')<jupyter_output><empty_output><jupyter_text>2.3.1 Age (Feature Engineering)Now that the admissions dataFrame has DOB information, age can be calculated.<jupyter_code># Find the first admission time for each patient
df_age_min = df[['SUBJECT_ID', 'ADMITTIME']].groupby('SUBJECT_ID').min().reset_index()
df_age_min.columns = ['SUBJECT_ID', 'ADMIT_MIN']
df_age_min.head()
df = df.merge(df_age_min, how='outer', on='SUBJECT_ID')
# Verify merge
df.head()
# Age is decode by finding the difference in admission date and date of birth
df['age'] = (df['ADMIT_MIN'] - df['DOB']).dt.days // 365
df['age'] = np.where(df['age'] < 0, 90, df['age'])
#df['age'] = np.where(df['age'] == -0, 0, df['age'])
df['age'].isnull().sum()
# Note that no ‘middle’ patients show up - this reflects the fact that MIMIC-III does not contain data from pediatric patients.
plt.hist(df['age'], bins=20, color='#c44e52')
plt.ylabel('Count')
plt.xlabel('Age (years)')
plt.title('Distribution of Age in MIMIC-III')
plt.tick_params(left=False, bottom=False, top=False, right=False) 
plt.show();
plt.scatter(df['age'], df['LOS'], alpha=0.005)
#plt.yscale('sqrt')
plt.ylabel('LOS (days)')
plt.xlabel('Age (years)')
plt.title('Age versus Length-of-stay')
plt.ylim(1, 50)
# https://en.wikipedia.org/wiki/List_of_ICD-9_codes
age_ranges = [(0, 13), (13, 36), (36, 56), (56, 100)]
for num, cat_range in enumerate(age_ranges):
    df['age'] = np.where(df['age'].between(cat_range[0],cat_range[1]), 
            num, df['age'])
    
age_dict = {0: 'newborn', 1: 'young_adult', 2: 'middle_adult', 3: 'senior'}
df['age'] = df['age'].replace(age_dict)
df.age.value_counts()<jupyter_output><empty_output><jupyter_text>2.3.2 Gender<jupyter_code>boxplot_los_groupby('GENDER', los_range=(0, 30))
df['GENDER'].replace({'M': 0, 'F':1}, inplace=True)<jupyter_output><empty_output><jupyter_text>2.4 ICUSTAYS.csv Exploration<jupyter_code># Intensive Care Unit (ICU) for each admission to hospital
df_icu.info()
df_icu['HADM_ID'].nunique()
df_icu.groupby('FIRST_CAREUNIT').median()<jupyter_output><empty_output><jupyter_text>2.4.1 ICU (Feature Engineering)<jupyter_code># Based on above statistics, reduce to just ICU and NICU groups
df_icu['FIRST_CAREUNIT'].replace({'CCU': 'ICU', 'CSRU': 'ICU', 'MICU': 'ICU',
                                  'SICU': 'ICU', 'TSICU': 'ICU'}, inplace=True)
df_icu['cat'] = df_icu['FIRST_CAREUNIT']
icu_list = df_icu.groupby('HADM_ID')['cat'].apply(list).reset_index()
icu_list.head()
df_icu['FIRST_CAREUNIT'].value_counts()
# Create admission-ICU matrix
icu_item = pd.get_dummies(icu_list['cat'].apply(pd.Series).stack()).sum(level=0)
icu_item[icu_item >= 1] = 1
icu_item = icu_item.join(icu_list['HADM_ID'], how="outer")
icu_item.head()
print("Number of admissions to ICU {}.".format(icu_item.ICU.sum()))
print("Number of admissions to NICU {}.".format(icu_item.NICU.sum()))
# Merge ICU data with main dataFrame
df = df.merge(icu_item, how='outer', on='HADM_ID')
# Replace NaNs with 0
df['ICU'].fillna(value=0, inplace=True)
df['NICU'].fillna(value=0, inplace=True)
# Verify NaN fix
print(df.ICU.value_counts(dropna=False))
print(df.NICU.value_counts(dropna=False))<jupyter_output>1.0    49794
0.0     9168
Name: ICU, dtype: int64
0.0    50970
1.0     7992
Name: NICU, dtype: int64<jupyter_text>3. Data PreprocessingEven after completing the feature engineering for age and ICD-9, there were some loose ends that needed tidying up before the data could be used for the prediction model. First, I ensured that no admissions resulting in death were part the cleaned dataset. I dropped all unused columns and verified that no NaNs existed in the data. For the admission type, insurance type, religion, ethnicity, age, and marital status columns, I performed the Pandas get_dummies command to convert these categorical variables into dummy/indicator variables. The final DataFrame size resulted in 48 feature columns and 1 target column with an entry count of 53,104.<jupyter_code># Look at what is no longer needed in the DataFrame
df.info()
# Remove deceased persons as they will skew LOS result
df = df[df['DECEASED'] == 0]

# Remove LOS with negative number, likely entry form error
df = df[df['LOS'] > 0]
# Drop unused or no longer needed columns
df.drop(columns=['SUBJECT_ID', 'HADM_ID', 'ADMITTIME', 'ADMISSION_LOCATION',
                'DISCHARGE_LOCATION', 'LANGUAGE', 'ADMIT_MIN', 'DOB',
                'DIAGNOSIS', 'DECEASED',  'DEATHTIME'], inplace=True)

df.info()
# Create dummy columns for categorical variables
prefix_cols = ['ADM', 'INS', 'REL', 'ETH', 'AGE', 'MAR', 'RELIGION']
dummy_cols = ['ADMISSION_TYPE', 'INSURANCE', 'RELIGION',
             'ETHNICITY', 'age', 'MARITAL_STATUS', 'RELIGION']
df = pd.get_dummies(df, prefix=prefix_cols, columns=dummy_cols)
df.info()
# Verify
df.head()<jupyter_output><empty_output><jupyter_text>3.1 Data Cleanup FunctionThe `mimic_los_cleanup` performs all of the above data cleanup and feature engineering steps.<jupyter_code>def mimic_los_cleanup(adm_csv='ADMISSIONS.csv', patients_csv='PATIENTS.csv',
                      diagcode_csv='DIAGNOSES_ICD.csv', icu_csv='ICUSTAYS.csv',
                      verbose=True):
    '''
    This function take 4 csv files from the MIMIC-III database, converts them to
    DataFrames for cleanup and feature engineering for use in a Length-of-Stay
    regression model such as the sklearn GradientBoostingRegressor. 
    
    INPUT:
    adm_csv - Primary Admissions information
    patients_csv - Patient specific info such as gender and DOB
    diagcode_csv - ICD9 Diagnosis for each admission to hospital
    icu_csv - Intensive Care Unit (ICU) data for each admission 
    
    OUTPUT:
    df - clean DataFrame for use in an regression model
    actual_median_los - Median LOS for all admissions 
    actual_mean_los - Average LOS for all admissions
    '''
    
    # Import CSV tables
    df = pd.read_csv(adm_csv)
    df_pat = pd.read_csv(patients_csv)
    df_diagcode = pd.read_csv(diagcode_csv)
    df_icu = pd.read_csv(icu_csv)
    if verbose: 
        print('(1/5) Completed .csv imports')
        
    # Feature Engineering for Length of Stay (LOS) target variable
    # Convert admission and discharge times to datatime type
    df['ADMITTIME'] = pd.to_datetime(df['ADMITTIME'])
    df['DISCHTIME'] = pd.to_datetime(df['DISCHTIME'])
    # Convert timedelta type into float 'days', 86400 seconds in a day
    df['LOS'] = (df['DISCHTIME'] - df['ADMITTIME']).dt.total_seconds()/86400
    
    # Drop columns that are not needed for next steps
    df.drop(columns=['DISCHTIME', 'ROW_ID', 
                    'EDREGTIME', 'EDOUTTIME', 'HOSPITAL_EXPIRE_FLAG',
                    'HAS_CHARTEVENTS_DATA'], inplace=True)
    
    # Track patients who died at the hospital by admission event
    df['DECEASED'] = df['DEATHTIME'].notnull().map({True:1, False:0})
    
    # Hospital LOS metrics
    actual_mean_los = df['LOS'].loc[df['DECEASED'] == 0].mean() 
    actual_median_los = df['LOS'].loc[df['DECEASED'] == 0].median() 
    
    # Compress the number of ethnicity categories
    df['ETHNICITY'].replace(regex=r'^ASIAN\D*', value='ASIAN', inplace=True)
    df['ETHNICITY'].replace(regex=r'^WHITE\D*', value='WHITE', inplace=True)
    df['ETHNICITY'].replace(regex=r'^HISPANIC\D*', value='HISPANIC/LATINO', inplace=True)
    df['ETHNICITY'].replace(regex=r'^BLACK\D*', value='BLACK/AFRICAN AMERICAN', inplace=True)
    df['ETHNICITY'].replace(['UNABLE TO OBTAIN', 'OTHER', 'PATIENT DECLINED TO ANSWER', 
                             'UNKNOWN/NOT SPECIFIED'], value='OTHER/UNKNOWN', inplace=True)
    df['ETHNICITY'].loc[~df['ETHNICITY'].isin(df['ETHNICITY'].value_counts().nlargest(5).index.tolist())] = 'OTHER/UNKNOWN'

    # Reduce categories to terms of religious or not
    df['RELIGION'].loc[~df['RELIGION'].isin(['NOT SPECIFIED', 'UNOBTAINABLE'])] = 'RELIGIOUS'

    # Re-categorize NaNs into 'Unknown'
    df['MARITAL_STATUS'] = df['MARITAL_STATUS'].fillna('UNKNOWN (DEFAULT)')
    
    if verbose: 
        print('(2/5) Completed ADMISSIONS.csv cleanup and feature engineering.')
        
    # Feature Engineering for ICD9 code categories
    # Filter out E and V codes since processing will be done on the numeric first 3 values
    df_diagcode['recode'] = df_diagcode['ICD9_CODE']
    df_diagcode['recode'] = df_diagcode['recode'][~df_diagcode['recode'].str.contains("[a-zA-Z]").fillna(False)]
    df_diagcode['recode'].fillna(value='999', inplace=True)
    df_diagcode['recode'] = df_diagcode['recode'].str.slice(start=0, stop=3, step=1)
    df_diagcode['recode'] = df_diagcode['recode'].astype(int)
    
    # ICD-9 Main Category ranges
    icd9_ranges = [(1, 140), (140, 240), (240, 280), (280, 290), (290, 320), (320, 390), 
                   (390, 460), (460, 520), (520, 580), (580, 630), (630, 680), (680, 710),
                   (710, 740), (740, 760), (760, 780), (780, 800), (800, 1000), (1000, 2000)]

    # Associated category names
    diag_dict = {0: 'infectious', 1: 'neoplasms', 2: 'endocrine', 3: 'blood',
                 4: 'mental', 5: 'nervous', 6: 'circulatory', 7: 'respiratory',
                 8: 'digestive', 9: 'genitourinary', 10: 'pregnancy', 11: 'skin', 
                 12: 'muscular', 13: 'congenital', 14: 'prenatal', 15: 'misc',
                 16: 'injury', 17: 'misc'}

    # Re-code in terms of integer
    for num, cat_range in enumerate(icd9_ranges):
        df_diagcode['recode'] = np.where(df_diagcode['recode'].between(cat_range[0],cat_range[1]), 
                num, df_diagcode['recode'])

    # Convert integer to category name using diag_dict
    df_diagcode['recode'] = df_diagcode['recode']
    df_diagcode['cat'] = df_diagcode['recode'].replace(diag_dict)
    
    # Create list of diagnoses for each admission
    hadm_list = df_diagcode.groupby('HADM_ID')['cat'].apply(list).reset_index()
    
    # Convert diagnoses list into hospital admission-item matrix
    hadm_item = pd.get_dummies(hadm_list['cat'].apply(pd.Series).stack()).sum(level=0)
    
    # Join back with HADM_ID, will merge with main admissions DF later
    hadm_item = hadm_item.join(hadm_list['HADM_ID'], how="outer")

    # Merge with main admissions df
    df = df.merge(hadm_item, how='inner', on='HADM_ID')
    
    if verbose: 
        print('(3/5) Completed DIAGNOSES_ICD.csv cleanup and feature engineering.')
    
    # Feature Engineering for Age and Gender
    # Convert to datetime type
    df_pat['DOB'] = pd.to_datetime(df_pat['DOB'])
    df_pat = df_pat[['SUBJECT_ID', 'DOB', 'GENDER']]
    df = df.merge(df_pat, how='inner', on='SUBJECT_ID')
    
    # Find the first admission time for each patient
    df_age_min = df[['SUBJECT_ID', 'ADMITTIME']].groupby('SUBJECT_ID').min().reset_index()
    df_age_min.columns = ['SUBJECT_ID', 'ADMIT_MIN']
    df = df.merge(df_age_min, how='outer', on='SUBJECT_ID')
    
    # Age is decode by finding the difference in admission date and date of birth
    df['age'] = (df['ADMIT_MIN'] - df['DOB']).dt.days // 365
    df['age'] = np.where(df['age'] < 0, 90, df['age'])
    
    # Create age categories
    age_ranges = [(0, 13), (13, 36), (36, 56), (56, 100)]
    for num, cat_range in enumerate(age_ranges):
        df['age'] = np.where(df['age'].between(cat_range[0],cat_range[1]), 
                num, df['age'])
    age_dict = {0: 'newborn', 1: 'young_adult', 2: 'middle_adult', 3: 'senior'}
    df['age'] = df['age'].replace(age_dict)
    
    # Re-map Gender to boolean type
    df['GENDER'].replace({'M': 0, 'F':1}, inplace=True)
    
    if verbose: 
        print('(4/5) Completed PATIENT.csv cleanup and feature engineering.')
    
    # Feature engineering for Intensive Care Unit (ICU) category
    # Reduce ICU categories to just ICU or NICU
    df_icu['FIRST_CAREUNIT'].replace({'CCU': 'ICU', 'CSRU': 'ICU', 'MICU': 'ICU',
                                  'SICU': 'ICU', 'TSICU': 'ICU'}, inplace=True)
    df_icu['cat'] = df_icu['FIRST_CAREUNIT']
    icu_list = df_icu.groupby('HADM_ID')['cat'].apply(list).reset_index()
    icu_item = pd.get_dummies(icu_list['cat'].apply(pd.Series).stack()).sum(level=0)
    icu_item[icu_item >= 1] = 1
    icu_item = icu_item.join(icu_list['HADM_ID'], how="outer")
    df = df.merge(icu_item, how='outer', on='HADM_ID')
    
    # Cleanup NaNs
    df['ICU'].fillna(value=0, inplace=True)
    df['NICU'].fillna(value=0, inplace=True)
    
    if verbose: 
        print('(5/5) Completed ICUSTAYS.csv cleanup and feature engineering.')
        
    # Remove deceased persons as they will skew LOS result
    df = df[df['DECEASED'] == 0]

    # Remove LOS with negative number, likely entry form error
    df = df[df['LOS'] > 0]    
    
    # Drop unused columns, e.g. not used to predict LOS
    df.drop(columns=['SUBJECT_ID', 'HADM_ID', 'ADMITTIME', 'ADMISSION_LOCATION',
                'DISCHARGE_LOCATION', 'LANGUAGE', 'ADMIT_MIN', 'DOB',
                'DIAGNOSIS', 'DECEASED', 'DEATHTIME'], inplace=True)
    
    prefix_cols = ['ADM', 'INS', 'REL', 'ETH', 'AGE', 'MAR']
    dummy_cols = ['ADMISSION_TYPE', 'INSURANCE', 'RELIGION',
                 'ETHNICITY', 'age', 'MARITAL_STATUS']
    df = pd.get_dummies(df, prefix=prefix_cols, columns=dummy_cols)
    
    if verbose: 
        print('Data Preprocessing complete.')
    
    return df, actual_median_los, actual_mean_los

df_clean, actual_median_los, actual_mean_los = mimic_los_cleanup(adm_csv='data/ADMISSIONS.csv', 
                                                                 patients_csv='data/PATIENTS.csv',
                                                                 diagcode_csv='data/DIAGNOSES_ICD.csv', 
                                                                 icu_csv='data/ICUSTAYS.csv')
# Verify output of cleaning function
df_clean.info(verbose=1)
# Check for any remaining NaNs
df_clean.isnull().values.sum()<jupyter_output><empty_output><jupyter_text>4. Length-of-Stay Prediction ModelIn this section I'll walkthrough the implementation and selection of the supervised learning prediction model using the Scikit-Learn machine learning library. 4.1 Model ImplementationTo implement the prediction model, I split the LOS target variable and features into training and test sets at an 80:20 ratio using the scikit-learn train_test_split function. Using the training set, I'll fit five different regression models (from the scikit-learn library) using default settings to see what the R2 score comparison looked like.<jupyter_code># Target Variable (Length-of-Stay)
LOS = df_clean['LOS'].values
# Prediction Features
features = df_clean.drop(columns=['LOS'])
# Split into train 80% and test 20%
X_train, X_test, y_train, y_test = train_test_split(features, 
                                                    LOS, 
                                                    test_size = .20, 
                                                    random_state = 0)

# Show the results of the split
print("Training set has {} samples.".format(X_train.shape[0]))
print("Testing set has {} samples.".format(X_test.shape[0]))
# Regression models for comparison
models = [SGDRegressor(random_state = 0), 
          GradientBoostingRegressor(random_state = 0), 
          LinearRegression(),
          KNeighborsRegressor(),
          RandomForestRegressor(random_state = 0)]

results = {}

for model in models:
    
    # Instantiate and fit Regressor Model
    reg_model = model
    reg_model.fit(X_train, y_train)
    
    # Make predictions with model
    y_test_preds = reg_model.predict(X_test)
    
    # Grab model name and store results associated with model
    name = str(model).split("(")[0]

    results[name] = r2_score(y_test, y_test_preds)
    print('{} done.'.format(name))
# R2 score results
fig, ax = plt.subplots()
ind = range(len(results))
ax.barh(ind, list(results.values()), align='center', 
        color = '#55a868', alpha=0.8)
ax.set_yticks(ind)
ax.set_yticklabels(results.keys())
ax.set_xlabel('R-squared score')
ax.tick_params(left=False, top=False, right=False) 
ax.set_title('Comparison of Regression Models')
fig.savefig('images/compare_models.png', bbox_inches = 'tight')<jupyter_output><empty_output><jupyter_text>The GradientBoostingRegressor took the win with an R2 score of ~37% so I decided to focus on refining this particular model. Because of past success with the RandomForestRegressor, I played with that model’s parameters but was never able to exceed the GradientBoostingRegressor score.<jupyter_code># GradientBoostingRegressor will be used as the LOS prediction model
reg_model = GradientBoostingRegressor(random_state=0)   
reg_model.fit(X_train, y_train)
y_test_preds = reg_model.predict(X_test)
r2_not_refined = r2_score(y_test, y_test_preds)
print("R2 score is: {:2f}".format(r2_not_refined))<jupyter_output>R2 score is: 0.367296<jupyter_text>4.2 Model Refinement The GradientBoostingRegressor performed well versus the other regression model. To refine the GradientBoostingRegressor model, I used the GridSearchCV function from scikit-learn to test out various permutations of parameters such as n_estimators, max_depth, and loss. The best estimator result from GridSearchCV was n_estimators=200, max_depth=4, and loss=ls.<jupyter_code># Split into train 80% and test 20%
X_train, X_test, y_train, y_test = train_test_split(features, 
                                                    LOS, 
                                                    test_size = .20, 
                                                    random_state = 42)

# Set the parameters by cross-validation
#tuned_parameters = [{'n_estimators': [100, 200, 300],
#                     'max_depth' : [2, 3, 4],
#                     'loss': ['ls', 'lad', 'huber']}]
tuned_parameters = [{'n_estimators': [200, 300],
                     'max_depth' : [3, 4],
                     'loss': ['ls', 'lad']}]

# create and fit a ridge regression model, testing each alpha
reg_model = GradientBoostingRegressor()  
grid = GridSearchCV(reg_model, tuned_parameters)
grid.fit(X_train, y_train)
reg_model_optimized = grid.best_estimator_

# summarize the results of the grid search
print(grid.best_score_)
print(grid.best_estimator_)
#reg_model = GradientBoostingRegressor(n_estimators = 200, max_depth=4, random_state=0)   
#reg_model.fit(X_train, y_train)
y_test_preds = reg_model_optimized.predict(X_test)
r2_optimized = r2_score(y_test, y_test_preds)
print("Optimized R2 score is: {:2f}".format(r2_optimized)) 
print('Model refinement improved R2 score by {:.4f}'.format(r2_optimized-r2_not_refined))<jupyter_output>Model refinement improved R2 score by 0.0190<jupyter_text>5. ResultsBefore looking at the RMSE benchmark, I wanted to investigate what features were most important in predicting hospital length-of-stay when using the gradient boosting regression model. Diagnoses related to prenatal issues have the highest feature importance coefficient followed by respiratory and injury. As I alluded to earlier, the ICD-9 diagnoses categories are by far the most important features. In fact, in the top 20 top features, only emergency admission type, gender, and Medicaid insurance showed any importance outside of diagnoses groups.<jupyter_code># https://towardsdatascience.com/running-random-forests-inspect-the-feature-importances-with-this-code-2b00dd72b92e
feature_imp = pd.DataFrame(reg_model_optimized.feature_importances_,
                                   index = X_train.columns,
                                   columns=['importance']).sort_values('importance', ascending=False)

feature_imp.head(20)
feature_imp.index[0:10].tolist()
# Plot feature importance
fig, ax = plt.subplots(figsize=(7, 5))
ind = range(0,10)
ax.barh(ind, feature_imp['importance'].values[0:10], 
        align='center', color='#c44e52', alpha=0.9)
ax.set_yticks(ind)
ax.set_yticklabels(feature_imp.index[0:10].tolist())
ax.tick_params(left=False, top=False, right=False) 
ax.set_title("Top 10 features for predicting LOS")
ax.set_xlabel('Feature Importance Coefficient \n(GradientBoostingRegressor)')
plt.gca().invert_yaxis()
fig.savefig('images/feature_importance.png', bbox_inches = 'tight')
fig, ax = plt.subplots(figsize=(10, 4.5))

# https://matplotlib.org/examples/api/barchart_demo.html
ind = np.arange(0,20)
pad = 0.15
width = 0.35
set_actual = ax.bar(pad+ind, LOS_actual, width, color='#4c72b0', alpha=0.8)
set_predict = ax.bar(pad+ind+width, LOS_predict, width, color='#dd8452')

ax.set_ylabel('Length of Stay (days)')
ax.set_xlabel('Admission Event')
ax.set_title('Prediction Model for Hospital Length of Stay (LOS)')
ax.text(0.5, 11, 'Average LOS', fontdict=None, withdash=False, color='red', alpha=0.8)
ax.text(0.5, 7.5, 'Median LOS', fontdict=None, withdash=False, color='gray')
ax.set_xticks(pad + ind + width)
ax.set_ylim(0, 30)
ax.set_xticklabels(list(range(20)))
ax.axhline(y=actual_median_los, xmin=0, xmax=20, ls='--', color='gray')
ax.axhline(y=actual_mean_los, xmin=0, xmax=20, ls='--', color='red', alpha=0.8)
ax.legend( (set_actual, set_predict), ('Actual LOS', 'Predicted LOS'), 
          loc='upper left')
ax.tick_params(bottom=False, top=False, right=False)<jupyter_output><empty_output><jupyter_text>While the RMSE trend is promising, I also wanted to evaluate the model from a few other perspectives. The following figure takes the first 20 admissions from the test set and directly compares the actual, predicted (gradient boosting model), average, and median LOS values. This gives a more convoluted picture of the prediction model; in some admissions, it predicts well but not as well in others.<jupyter_code>#y_test_preds = reg_model.predict(X_test)

ml_count, md_count, avg_count  = 0, 0, 0
ml_days, md_days, avg_days  = 0, 0, 0
ml_days_rms, md_days_rms, avg_days_rms  = 0, 0, 0

for i in range(y_test_preds.shape[0]):
    ml_model = abs(y_test_preds[i] - y_test[i])
    median_model = abs(actual_median_los - y_test[i])
    average_model = abs(actual_avg_los - y_test[i])
    
    ml_days += ml_model
    md_days += median_model
    avg_days += average_model
    
    ml_model_rms = (y_test_preds[i] - y_test[i])**2
    median_model_rms = (actual_median_los - y_test[i])**2
    average_model_rms = (actual_avg_los - y_test[i])**2
    
    ml_days_rms += ml_model_rms
    md_days_rms += median_model_rms
    avg_days_rms += average_model_rms

print("Prediction Model days {}".format(ml_days/y_test_preds.shape[0]))
print("Median Model days {}".format(md_days/y_test_preds.shape[0]))
print("Average Model days {}".format(avg_days/y_test_preds.shape[0]))

print("Prediction Model RMS {}".format((ml_days_rms**0.5)/y_test_preds.shape[0]))
print("Median Model RMS {}".format((md_days_rms**0.5)/y_test_preds.shape[0]))
print("Average Model RMS {}".format((avg_days_rms**0.5)/y_test_preds.shape[0]))<jupyter_output>Prediction Model days 5.443111108662211
Median Model days 6.378620461036283
Average Model days 7.238142660243353
Prediction Model RMS 0.09538079881644249
Median Model RMS 0.12679346060609392
Average Model RMS 0.12180964220883837
Prediction Model r2 score: 0.38681380810675525
Median LOS Model r2 score: -0.08358794923765545
Average LOS Model r2 score: -7.78096559899133e-05<jupyter_text>In the metrics section, I stated that the RMSE would be used to compare the prediction model versus the industry standard average and median LOS metrics. The gradient boosting model RMSE is better by more than 24% (percent difference) versus the constant average or median models.<jupyter_code># RMSE plot for writeup
data = pd.DataFrame({'RMSE': [(ml_days_rms**0.5)/y_test_preds.shape[0], 
                             (avg_days_rms**0.5)/y_test_preds.shape[0],
                             (md_days_rms**0.5)/y_test_preds.shape[0]],
                     'LOS Model Type': ['Gradient Boosting', 'Average', 'Median'] })

fig, ax = plt.subplots()
ax = sns.barplot(x='RMSE', y='LOS Model Type', data=data)
ax.set_title('RMSE comparison of Length-of-Stay models')
ax.tick_params(top=False, left=False, right=False)

fig.savefig('images/rms_comparison.png', bbox_inches = 'tight')<jupyter_output><empty_output><jupyter_text>The final way I wanted to look at the model was to plot the proportion of accurate predictions in the test set versus an allowed margin of error. Other studies qualify a LOS prediction as correct if it falls within a certain margin of error. It follows that as the margin of error allowance increases, so should the proportion of accurate predictions for all models. The gradient boosting prediction model performs better than the other constant models across the margin of error range up to 50%.<jupyter_code># Calculate Proportion of 'accurate' prediction as a function of allowed margin of error
reg_array = []
median_array = []
average_array = []

for i in list(range(6)):
    reg_count, median_count, average_count = 0, 0, 0
    
    for j in range(y_test_preds.shape[0]):
        # Percent Difference
        reg_model = (y_test_preds[j] - y_test[j])/y_test[j]
        median_model = (actual_median_los - y_test[j])/y_test[j]
        average_model = (actual_avg_los - y_test[j])/y_test[j]
        if abs(reg_model) < i/10:
            reg_count += 1
        if abs(median_model) < i/10:
            median_count += 1
        if abs(average_model) < i/10:
            average_count += 1
            
    reg_array.append((reg_count/y_test_preds.shape[0]))
    median_array.append((median_count/y_test_preds.shape[0]))
    average_array.append((average_count/y_test_preds.shape[0]))

# Plot proportion of 'accurate' prediction as a function of allowed margin of error  
fig, ax = plt.subplots()
ax.plot(reg_array, label='Gradient Boosting')
ax.plot(median_array, label='Median LOS model')
ax.plot(average_array, label='Average LOS model')
ax.set_title('Proportion of Accurate Predictions vs. Percent Error')
ax.set_xlabel('Allowed Margin of Error (Percent Error)')
ax.set_ylabel('Proportion of Accurate Predictions')
ax.set_xticklabels(['0%', '10%', '20%', '30%', '40%', '50%'])
ax.legend(loc='lower right');
ax.tick_params(top=False, right=False)
fig.savefig('images/rms_comparison.png', bbox_inches = 'tight')<jupyter_output>6.46
10.2
{'fraction_non_alphanumeric': 0.12201386568660207, 'fraction_numerical': 0.025851987971934515, 'mean_word_length': 4.573024554870243, 'pattern_counts': {'":': 0, '<': 111, '<?xml version=': 0, '>': 110, 'https://': 11, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 25, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '27', 'dir': 'jupyter-structured-clean-dedup', 'id': '240660', 'license': 'MIT', 'n_tokens_mistral': 16624, 'n_tokens_neox': 15136, 'n_words': 4522, 'path': 'hospital_los_prediction.ipynb', 'repo_name': 'satyaki-bhargava/Neonate-los-predictor', 'size': '47914'}