<jupyter_start><jupyter_text>Appendix BThe following contains code for Appendix B of the paper, specifically, evaluating the performance of the models with and without certain covariates.<jupyter_code>from __future__ import print_function

# Import libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import psycopg2
import sys, os
import subprocess

# hack to append the parent path to the system path
# this allows us to import sepsis_utils
sys.path.append(os.path.dirname(os.path.abspath('.')))

from sepsis_utils import sepsis_utils as su
from sepsis_utils import roc_utils as ru

from statsmodels.formula.api import logit

from IPython.display import display, HTML # used to print out pretty pandas dataframes

# default colours for prettier plots
col = [[0.9047, 0.1918, 0.1988],
    [0.2941, 0.5447, 0.7494],
    [0.3718, 0.7176, 0.3612],
    [1.0000, 0.5482, 0.1000],
    [0.4550, 0.4946, 0.4722],
    [0.6859, 0.4035, 0.2412],
    [0.9718, 0.5553, 0.7741],
    [0.5313, 0.3359, 0.6523]];
marker = ['v','o','d','^','s','o','+']
ls = ['-','-','-','-','-','s','--','--']
%matplotlib inline<jupyter_output><empty_output><jupyter_text>First, load the data and specify the scores we are interested in.<jupyter_code># load data
df = pd.read_csv('../sepsis3-df.csv',sep=',')
df_mdl = pd.read_csv('../sepsis3-design-matrix.csv',sep=',')

preds_header = ['sirs','qsofa','sofa','mlods']
target_header = 'hospital_expire_flag'<jupyter_output><empty_output><jupyter_text>For convenience, we reproduce the table of AUROC values from the paper using logistic regression.<jupyter_code># reproduce the AUC table

# print out the baseline AUROC of a model without severity of illness
# in Seymour2016, this was 0.58 (validation)

fn_in = "../sepsis3-design-matrix.csv"
fn_out = "../sepsis3-preds.csv"
formula = target_header + " ~ fp(age) + fp(elixhauser_hospital) + is_male + race_black + race_other"
rcmd = ["Rscript ../r-make-sepsis3-models.R", fn_in, fn_out, '"' + formula + '"']
err = subprocess.call(' '.join(rcmd), shell=True)
if err!=0:
    print('RScript returned error status {}.'.format(err))
else:
    preds_mfp_b = pd.read_csv(fn_out, sep=',', header=0)
    preds_mfp_b = preds_mfp_b.values[:,0]
    auc_mfp_b, ci_mfp_b = ru.calc_auc(preds_mfp_b, df[target_header].values,
                                      with_ci=True,
                                      alpha=0.05)
    print('Performance of an MFP model without severity of illness: {:0.3f} [{:0.3f}, {:0.3f}]'.format(
            auc_mfp_b, ci_mfp_b[0], ci_mfp_b[1]))
reload(su)
print('AUROC table for adjusted scores (using logistic regression):')
preds_logreg = su.calc_predictions(df, preds_header, target_header, model='logreg')
su.print_auc_table(preds_logreg, df[target_header].values, preds_header, with_alpha=False)
print('\n')<jupyter_output>AUROC table for adjusted scores (using logistic regression):
     	sirs                	qsofa               	sofa                	mlods               	
sirs 	0.729 [0.712, 0.746]	                    	                    	                    	
qsofa	0.042               	0.718 [0.701, 0.735]	                    	                    	
sofa 	< 0.001               	< 0.001               	0.783 [0.767, 0.799]	                    	
mlods	< 0.001               	< 0.001               	0.279               	0.788 [0.773, 0.803]<jupyter_text>We examine a logistic regression model using qSOFA in addition to the baseline covariates.<jupyter_code># print a summary of the qsofa model using logistic regression (Table 4)
print('Summary of a logistic regression model using qSOFA:')
_ = su.calc_predictions(df, ['qsofa'], target_header, model='logreg', print_summary=True)
print('\n\n')
print('Summary of a logistic regression model using qSOFA *components*:')
model = logit(formula=target_header + " ~ age + elixhauser_hospital" +
              " + race_black + race_other + is_male + " +
              " + qsofa_sysbp_score + qsofa_gcs_score + qsofa_resprate_score",
              data=df).fit(disp=0)
smry = model.summary2()
smry = smry.tables[1]
smry['Odds Ratio'] = np.exp(smry['Coef.'])

display(HTML(smry[['Coef.','Odds Ratio','P>|z|','[0.025','0.975]']].to_html()))<jupyter_output>Summary of a logistic regression model using qSOFA:
                            Logit Regression Results                            
================================================================================
Dep. Variable:     hospital_expire_flag   No. Observations:                 6781
Model:                            Logit   Df Residuals:                     6774
Method:                             MLE   Df Model:                            6
Date:                  Wed, 23 Nov 2016   Pseudo R-squ.:                 0.09039
Time:                          18:46:54   Log-Likelihood:                -2495.6
converged:                         True   LL-Null:                       -2743.6
                                          LLR p-value:                6.104e-104
=======================================================================================
                          coef    std err          z      P>|z|      [0.025      0.975]
-------------------------------------------[...]<jupyter_text>We now evaluate the performance of the model as we omit covariates.<jupyter_code># score to evaluate
score_name = "qsofa"

# below used for training MFP model... leave the text as is since it is hard coded in the R script
fn_in = "sepsis3-design-matrix" # the function adds "-dev.csv" and "-val.csv"
fn_out = "sepsis3-preds"

model_covariates = ["age","elixhauser_hospital","race_black","race_other","is_male", score_name]

# build full model
formula = target_header + " ~ " + " + ".join([x for x in model_covariates])
model = logit(formula=formula,data=df).fit(disp=0)
pred = model.predict(exog=df)

# remove all but the odds ratio column (with p value)
smry = model.summary2()
smry = smry.tables[1]
smry['Odds Ratio'] = smry.apply(lambda x: '{:1.2f}, p={:1.3f}'.format(np.exp(x['Coef.']),x['P>|z|']),axis=1)
smry.drop(['Coef.','Std.Err.','z','P>|z|','[0.025','0.975]'],inplace=True,axis=1)

# add AUROC to the bottom
auc, ci = ru.calc_auc(pred, df[target_header].values, with_ci=True)
smry.loc['auc','Odds Ratio'] = '{:0.2f}, {:0.2f}-{:0.2f}'.format(auc,ci[0],ci[1])

for m in model_covariates:
    formula = target_header + " ~ " + " + ".join([x for x in model_covariates if x!=m])
    model = logit(formula=formula,data=df).fit(disp=0)
    pred = model.predict(exog=df)
    auc, ci = ru.calc_auc(pred.values,df[target_header].values, with_ci=True)
    
    smry2 = model.summary2()
    smry2 = smry2.tables[1]
    smry2['Odds Ratio'] = smry2.apply(lambda x: '{:1.2f}, p={:1.3f}'.format(np.exp(x['Coef.']),x['P>|z|']),axis=1)
    smry2.loc['auc','Odds Ratio'] = '{:0.2f}, {:0.2f}-{:0.2f}'.format(auc,ci[0],ci[1])
    
    # add this as another column of smry
    smry.loc[smry2.index,m] = smry2['Odds Ratio']
    
display(HTML(smry.to_html().replace('NaN','')))<jupyter_output><empty_output><jupyter_text>We can do the same as the above with development and validation sets.<jupyter_code># score to evaluate
score_name = "qsofa"

# split into development/validation sets
df_dev = df.sample(frac=0.5, replace=False, weights=None, random_state=12875, axis=0)
df_val = df.loc[[x for x in df.index if x not in df_dev.index],:]

# below used for training MFP model... leave the text as is since it is hard coded in the R script
fn_in = "sepsis3-design-matrix" # the function adds "-dev.csv" and "-val.csv"
fn_out = "sepsis3-preds"

model_covariates = ["age","elixhauser_hospital","race_black","race_other","is_male",score_name]

# build full model
formula = target_header + " ~ " + " + ".join([x for x in model_covariates])
model = logit(formula=formula,data=df_dev).fit(disp=0)
pred = model.predict(exog=df_val)

# remove all but the odds ratio column (with p value)
smry = model.summary2()
smry = smry.tables[1]
smry['Odds Ratio'] = smry.apply(lambda x: '{:1.2f}, p={:1.3f}'.format(np.exp(x['Coef.']),x['P>|z|']),axis=1)
smry.drop(['Coef.','Std.Err.','z','P>|z|','[0.025','0.975]'],inplace=True,axis=1)

# add AUROC to the bottom
auc, ci = ru.calc_auc(pred, df_val[target_header].values, with_ci=True)
smry.loc['auc','Odds Ratio'] = '{:0.2f}, {:0.2f}-{:0.2f}'.format(auc,ci[0],ci[1])

for m in model_covariates:
    formula = target_header + " ~ " + " + ".join([x for x in model_covariates if x!=m])
    model = logit(formula=formula,data=df_dev).fit(disp=0)
    pred = model.predict(exog=df_val)
    auc, ci = ru.calc_auc(pred.values,df_val[target_header].values, with_ci=True)
    
    smry2 = model.summary2()
    smry2 = smry2.tables[1]
    smry2['Odds Ratio'] = smry2.apply(lambda x: '{:1.2f}, p={:1.3f}'.format(np.exp(x['Coef.']),x['P>|z|']),axis=1)
    smry2.loc['auc','Odds Ratio'] = '{:0.2f}, {:0.2f}-{:0.2f}'.format(auc,ci[0],ci[1])
    
    # add this as another column of smry
    smry.loc[smry2.index,m] = smry2['Odds Ratio']
    
    
display(HTML(smry.to_html().replace('NaN','')))<jupyter_output><empty_output><jupyter_text>We can further verify that this likely occurs with the MFP model as well by comparing the AUROCs between the MFP model and the logistic regression model (this training can take some time).<jupyter_code># score to evaluate
score_name = "qsofa"

# split into development/validation sets
df_dev = df.sample(frac=0.5, replace=False, weights=None, random_state=12875, axis=0)
df_val = df.loc[[x for x in df.index if x not in df_dev.index],:]

# below used for training MFP model... leave the text as is since it is hard coded in the R script
fn_in = "sepsis3-design-matrix" # the function adds "-dev.csv" and "-val.csv"
fn_out = "sepsis3-preds"

model_covariates = ["age","elixhauser_hospital","race_black","race_other","is_male",score_name]


# build full model
formula = target_header + " ~ " + " + ".join([x for x in model_covariates])
model = logit(formula=formula,data=df_dev).fit(disp=0)
pred = model.predict(exog=df_val)

# remove all but the odds ratio column (with p value)
smry = model.summary2()
smry = smry.tables[1]
smry['Odds Ratio'] = smry.apply(lambda x: '{:1.2f}, p={:1.3f}'.format(np.exp(x['Coef.']),x['P>|z|']),axis=1)
smry.drop(['Coef.','Std.Err.','z','P>|z|','[0.025','0.975]'],inplace=True,axis=1)

# add AUROC to the bottom
auc, ci = ru.calc_auc(pred, df_val[target_header].values, with_ci=True)
smry.loc['auc','Odds Ratio'] = '{:0.2f}, {:0.2f}-{:0.2f}'.format(auc,ci[0],ci[1])

# also get performance if we were using an MFP model
formula = formula.replace(" age ", " fp(age) ").replace(" elixhauser_hospital "," fp(elixhauser_hospital) ")
formula = formula.replace(" " + score_name, " fp(" + score_name + ")")
rcmd = ["Rscript ../r-make-sepsis3-models-dev-val.R", fn_in, fn_out, '"' + formula + '"']
err = subprocess.call(' '.join(rcmd), shell=True)
if err!=0:
    print('RScript returned error status {}.'.format(err))

# load in the predictions
pred_val = pd.read_csv(fn_out + '-val.csv', sep=',', header=0)
pred_val = pred_val.values[:,0]
tar_val = pd.read_csv(fn_out + '-val-tar.csv', sep=',', header=0)
tar_val = tar_val.values[:,0]

auc, ci = ru.calc_auc(pred_val, tar_val, with_ci=True)

# add AUROC to the bottom
smry.loc['auc_mfp','Odds Ratio'] = '{:0.2f}, {:0.2f}-{:0.2f}'.format(auc,ci[0],ci[1])

for m in model_covariates:
    formula = target_header + " ~ " + " + ".join([x for x in model_covariates if x!=m])
    model = logit(formula=formula,data=df_dev).fit(disp=0)
    pred = model.predict(exog=df_val)
    auc, ci = ru.calc_auc(pred.values,df_val[target_header].values, with_ci=True)
    
    smry2 = model.summary2()
    smry2 = smry2.tables[1]
    smry2['Odds Ratio'] = smry2.apply(lambda x: '{:1.2f}, p={:1.3f}'.format(np.exp(x['Coef.']),x['P>|z|']),axis=1)
    smry2.loc['auc','Odds Ratio'] = '{:0.2f}, {:0.2f}-{:0.2f}'.format(auc,ci[0],ci[1])
    

    # also get performance if we were using an MFP model
    formula = formula.replace(" age ", " fp(age) ").replace(" elixhauser_hospital "," fp(elixhauser_hospital) ")
    formula = formula.replace(" " + score_name, " fp(" + score_name + ")")
    rcmd = ["Rscript ../r-make-sepsis3-models-dev-val.R", fn_in, fn_out, '"' + formula + '"']
    err = subprocess.call(' '.join(rcmd), shell=True)
    if err!=0:
        print('RScript returned error status {}.'.format(err))

    # load in the predictions
    pred_val = pd.read_csv(fn_out + '-val.csv', sep=',', header=0)
    pred_val = pred_val.values[:,0]
    tar_val = pd.read_csv(fn_out + '-val-tar.csv', sep=',', header=0)
    tar_val = tar_val.values[:,0]

    smry2.loc['auc_mfp','Odds Ratio'] = '{:0.2f}, {:0.2f}-{:0.2f}'.format(auc,ci[0],ci[1])
    
    # add this as another column of smry
    smry.loc[smry2.index,m] = smry2['Odds Ratio']
    
    
display(HTML(smry.to_html().replace('NaN','')))<jupyter_output><empty_output><jupyter_text>The above tells us:* omitting race and gender have very little effect on model performance* omitting qsofa reduces model performance by a small fraction* omitting elixhauser reduces model AUROC by 0.05For comparison, we can repeat this with SOFA.<jupyter_code># score to evaluate
score_name = "sofa"

# split into development/validation sets
df_dev = df.sample(frac=0.5, replace=False, weights=None, random_state=12875, axis=0)
df_val = df.loc[[x for x in df.index if x not in df_dev.index],:]

# below used for training MFP model... leave the text as is since it is hard coded in the R script
fn_in = "sepsis3-design-matrix" # the function adds "-dev.csv" and "-val.csv"
fn_out = "sepsis3-preds"

model_covariates = ["age","elixhauser_hospital","race_black","race_other","is_male",score_name]


# build full model
formula = target_header + " ~ " + " + ".join([x for x in model_covariates])
model = logit(formula=formula,data=df_dev).fit(disp=0)
pred = model.predict(exog=df_val)

# remove all but the odds ratio column (with p value)
smry = model.summary2()
smry = smry.tables[1]
smry['Odds Ratio'] = smry.apply(lambda x: '{:1.2f}, p={:1.3f}'.format(np.exp(x['Coef.']),x['P>|z|']),axis=1)
smry.drop(['Coef.','Std.Err.','z','P>|z|','[0.025','0.975]'],inplace=True,axis=1)

# add AUROC to the bottom
auc, ci = ru.calc_auc(pred, df_val[target_header].values, with_ci=True)
smry.loc['auc','Odds Ratio'] = '{:0.2f}, {:0.2f}-{:0.2f}'.format(auc,ci[0],ci[1])

# also get performance if we were using an MFP model
formula = formula.replace(" age ", " fp(age) ").replace(" elixhauser_hospital "," fp(elixhauser_hospital) ")
formula = formula.replace(" " + score_name, " fp(" + score_name + ")")
rcmd = ["Rscript ../r-make-sepsis3-models-dev-val.R", fn_in, fn_out, '"' + formula + '"']
err = subprocess.call(' '.join(rcmd), shell=True)
if err!=0:
    print('RScript returned error status {}.'.format(err))

# load in the predictions
pred_val = pd.read_csv(fn_out + '-val.csv', sep=',', header=0)
pred_val = pred_val.values[:,0]
tar_val = pd.read_csv(fn_out + '-val-tar.csv', sep=',', header=0)
tar_val = tar_val.values[:,0]

auc, ci = ru.calc_auc(pred_val, tar_val, with_ci=True)

# add AUROC to the bottom
smry.loc['auc_mfp','Odds Ratio'] = '{:0.2f}, {:0.2f}-{:0.2f}'.format(auc,ci[0],ci[1])

for m in model_covariates:
    formula = target_header + " ~ " + " + ".join([x for x in model_covariates if x!=m])
    model = logit(formula=formula,data=df_dev).fit(disp=0)
    pred = model.predict(exog=df_val)
    auc, ci = ru.calc_auc(pred.values,df_val[target_header].values, with_ci=True)
    
    smry2 = model.summary2()
    smry2 = smry2.tables[1]
    smry2['Odds Ratio'] = smry2.apply(lambda x: '{:1.2f}, p={:1.3f}'.format(np.exp(x['Coef.']),x['P>|z|']),axis=1)
    smry2.loc['auc','Odds Ratio'] = '{:0.2f}, {:0.2f}-{:0.2f}'.format(auc,ci[0],ci[1])
    

    # also get performance if we were using an MFP model
    formula = formula.replace(" age ", " fp(age) ").replace(" elixhauser_hospital "," fp(elixhauser_hospital) ")
    formula = formula.replace(" " + score_name, " fp(" + score_name + ")")
    rcmd = ["Rscript ../r-make-sepsis3-models-dev-val.R", fn_in, fn_out, '"' + formula + '"']
    err = subprocess.call(' '.join(rcmd), shell=True)
    if err!=0:
        print('RScript returned error status {}.'.format(err))

    # load in the predictions
    pred_val = pd.read_csv(fn_out + '-val.csv', sep=',', header=0)
    pred_val = pred_val.values[:,0]
    tar_val = pd.read_csv(fn_out + '-val-tar.csv', sep=',', header=0)
    tar_val = tar_val.values[:,0]

    smry2.loc['auc_mfp','Odds Ratio'] = '{:0.2f}, {:0.2f}-{:0.2f}'.format(auc,ci[0],ci[1])
    
    # add this as another column of smry
    smry.loc[smry2.index,m] = smry2['Odds Ratio']
    
    
display(HTML(smry.to_html().replace('NaN','')))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18501227030585982, 'fraction_numerical': 0.03423714610642246, 'mean_word_length': 3.6566332218506132, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 45, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '140139', 'license': 'MIT', 'n_tokens_mistral': 6603, 'n_tokens_neox': 6192, 'n_words': 1362, 'path': 'appendix/sepsis-3-appendix-b.ipynb', 'repo_name': 'sreered/sepsis3-mimic', 'size': '16707'}