<jupyter_start><jupyter_text>Home Credit - ModellingIn this notebook, we will1. select an appropriate __evaluation matrix__2. devise our __modelling and validation strategy__3. select the __best model__ and4. evaluate __business impact__We will use the tidy data with newly derived features from the Feature Engineering notebook.<jupyter_code>import numpy as np 
import pandas as pd 
import matplotlib.pyplot as plt 
from matplotlib.ticker import PercentFormatter
import seaborn as sns
from scipy.stats import loguniform

from sklearn.model_selection import KFold, StratifiedKFold
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.model_selection import RandomizedSearchCV, cross_validate


from sklearn.dummy import DummyClassifier
from sklearn.linear_model import LogisticRegression
import xgboost as xgb
import lightgbm as lgb
from sklearn.metrics import roc_auc_score, roc_curve, confusion_matrix, classification_report

%matplotlib inline
%precision %.2f

plt.style.use('bmh')
pd.set_option('display.max_rows', 30)
pd.set_option('display.min_rows', 10)
pd.set_option('display.max_columns', 100)

# We will work with ~25% data in the process development phase.
process_dev_flag = False
random_seed = 79
lr_iterations = 10
xgb_iterations = 50
lgb_iterations = 50
train_path = '../data/interim/df_train_fe.csv'
test_path = '../data/interim/df_test_fe.csv'
dtype_path = '../data/interim/data_types_fe.csv'
df_train = pd.read_csv(train_path,index_col=0)
df_test = pd.read_csv(test_path,index_col=0)

df_dtype = pd.read_csv(dtype_path,index_col=0)
dict_dtype = df_dtype.dtype.to_dict()

df_train.shape, df_test.shape
df_train = df_train.astype(dict_dtype)
df_test = df_test.astype(dict_dtype)
df_train.info()
df_test.info()
df_train.head()<jupyter_output><empty_output><jupyter_text>Getting the list of categorical features with more than 2 levels.<jupyter_code>cat_object_var_path = '../data/interim/cat_encode_features.csv'
df_cat_object = pd.read_csv(cat_object_var_path,index_col=0)
cat_object_features = list(df_cat_object.cat_object)
len(cat_object_features)<jupyter_output><empty_output><jupyter_text>Getting the list of categorical features with just 2 levels. These features are stored as integers.<jupyter_code>cat_bin_path = '../data/interim/cat_bin_num_features.csv'
df_cat_bin_num = pd.read_csv(cat_bin_path,index_col=0)
cat_bin_num_features = list(df_cat_bin_num.cat_bin_num)
len(cat_bin_num_features)<jupyter_output><empty_output><jupyter_text>Finally, remaining features are numerical<jupyter_code>num_features = set(df_train.columns) - (set(cat_object_features) | set(cat_bin_num_features))
num_features = list(num_features)
num_features.remove('sk_id_curr')
num_features.remove('target')
print(len(num_features))<jupyter_output>124<jupyter_text>Evaluation matrixFor Home Credit, it is equally important to  1. Avoid customers who are likely to default on the loan and2. Increase customers who are likely to pay back the loanGiven this business objective, we will be __using the roc_auc__ matrix as it considers both the classes(default/no default). The pr_auc matrix mainly cares about positive class(default). Hence, we will avoid it.  Other important considerations in model evaluation will be __prediction time and training time__ in that order.  Modelling StrategyHomeCredit is more interested in accurate predictions than interpretations. Hence, we consider following algorithms for our classification task. 1. __Logistic Regression with l1 regularization__ - Logistic regression is a simple linear algorithm with decent performance. As we have 180+ features and many are correlated features, we will use l1 regularization for feature selection. This will also be our base model.2. __XgBoost__ - XgBoost is known to outperform linear models on tabular data in many cases.3. __LightGBM__ - LightGBM is faster than XgBoost and tends to perform well on tabular data.Logistic regression with manual/automatic feature selection resulting in parameter coefficients would be a better choice if transparency and simplicity was the preference.  Handling missing values1. __Categorical__ features - Create a seperate nan category/level for missing value. We will do this at the time of one-hot-encoding using pandas get_dummy function.2. __Numerical__ features - For logistic regression, we will impute missing values with median because numerical features are highly skewed as observed in EDA. XgBoost & LightGBM can handle missing in a supervised way to optimize performance. Hence, we will not impute numerical features for these two algorithms.   Creating dummy variables<jupyter_code>df_train_dummies = pd.get_dummies(df_train,columns = cat_object_features, dummy_na=True)
df_test_dummies = pd.get_dummies(df_test,columns = cat_object_features, dummy_na=True)<jupyter_output><empty_output><jupyter_text>Aligning test data with train data so that it has same feaures.<jupyter_code>#source : https://kanoki.org/2019/08/27/pandas-dataframe-align-function/
df_train_dummies, df_test_dummies = df_train_dummies.align(df_test_dummies, join = 'left', axis = 1)
df_train_dummies.shape, df_test_dummies.shape
new_features = [ col for col in df_train_dummies.columns if col not in df_train.columns]
X_train = df_train_dummies.drop(columns=['target','sk_id_curr'])
y_train = df_train_dummies.target

if process_dev_flag == True :
    records = 50000
    X_train = X_train[:records]
    y_train = y_train[:records]

X_test = df_test_dummies.drop(columns=['target','sk_id_curr'])
y_test = df_test_dummies.target<jupyter_output><empty_output><jupyter_text>Validation StrategyWe will use a __5-fold__ cross-validation strategy on train data for hyperparameter tuning and model selection. We will train these 3 models on same set train-validation datasets. To select the best hyperparameters, we will use __RandomSearchCV__ technique. And finally, we will test the performance of the best model on test the dataset. This will give us better idea of how well the model performs on unseen data.<jupyter_code>num_folds = 5
cv = KFold(n_splits= num_folds, shuffle=True, random_state=random_seed)<jupyter_output><empty_output><jupyter_text>Logistic regression We will use median imputation to fill missing values for numerical features using column transformer. As we are using l1 regularization and for greater convergence speed, we will scale these features using StandardScaler.We will tune the C (regularization) parameter using RandomSearchCV with 10 iterations.<jupyter_code>numeric_transformer = Pipeline(steps=[
                                ('imputer', SimpleImputer(strategy="median")),
                                ('scaler',StandardScaler())])

preprocessor = ColumnTransformer( 
                    transformers=[('num', numeric_transformer, num_features)], 
                    remainder = 'passthrough',
                    sparse_threshold=0 ) # to get the dense matrix

# saga works with l1 penalty and tol is set to 1e-2 for faster convergence
classifier_lr = LogisticRegression(penalty='l1', solver = 'saga', 
                                   max_iter=500, tol = 1e-2, 
                                   random_state=random_seed)

clf_lr = Pipeline(steps=[
                    ('preprocessor', preprocessor),
                    ('classifier_lr', classifier_lr)])

param_dists = { 'classifier_lr__C': loguniform(1e-2, 1e2) }

%%time
cv_lr = RandomizedSearchCV(clf_lr, 
                           param_distributions=param_dists, 
                           scoring = 'roc_auc',
                           n_iter=lr_iterations, 
                           cv=cv, 
                           n_jobs= -2,
                           random_state = random_seed)

cv_lr.fit(X_train, y_train)
cv_lr_results = pd.DataFrame(cv_lr.cv_results_)
# modified from https://inria.github.io
def shorten_param_names_pipeline(param_name):
    """
    Shorten the hyperparameter names from CrossValidation results of pipeline.
    
    Keyword arguments:
    param_name -- column name of cv_results_ dataframe
    """
    
    if '__' in param_name:
        return param_name.rsplit("__", 1)[1]
    return param_name

def shorten_param_names(param_name):
    """
    Shorten the hyperparameter names from CrossValidation results of model.
    
    Keyword arguments:
    param_name -- column name of cv_results_ dataframe
    """
    
    if 'param_' in param_name:
        return param_name.replace("param_", "")
    return param_name
# to get only the parameters columns
col_results = column_results = ['param_{}'.format(name) for name in param_dists.keys()]

# add score columns
column_results += ['mean_test_score', 'std_test_score', 'rank_test_score']

cv_lr_results = cv_lr_results[column_results].sort_values(
                            'mean_test_score', ascending=False)

cv_lr_results.columns = [shorten_param_names_pipeline(col) for col in cv_lr_results.columns]
cv_lr_results<jupyter_output><empty_output><jupyter_text>Top logistic regression model has 0.7643 mean AUC score. Also, top 3 models have same mean & standard deviation indicating the stability of the model. We have reduced tolerance level to 0.01 for faster convergence with less iterations. This will be our base model.  XgBoost   As XgBoost has more hyperparameters to tune, we will use __50 iterations__ in RandomSearch cross-validation. I have enabled GPU support for XgBoost training. 50 iterations are chosen considering the __time and hardware availability__.<jupyter_code>gpu_params = {'tree_method':'gpu_hist', 
              'predictor':'gpu_predictor',
              'n_jobs': -1}

xgb_clf = xgb.XGBClassifier(objective= 'binary:logistic', max_delta_step = 1,
                             **gpu_params, use_label_encoder=False, 
                            seed=random_seed, verbosity=0)

xgb_param_dists = {'learning_rate': loguniform(0.001, 1),
               'max_depth' : [1,4,8,16],
               'min_child_weight' : [1,4,8],
               'colsample_bytree' : [0.6,0.8],
               'subsample' : [0.6,0.8],
               'n_estimators' : [200,400,800,1600,6400],
               'reg_alpha' : loguniform(0.001, 10) }

%%time
xgb_cv = RandomizedSearchCV(xgb_clf, 
                           param_distributions=xgb_param_dists, 
                           scoring = 'roc_auc',
                           n_iter=xgb_iterations, 
                           cv=cv, 
                           random_state = random_seed)

xgb_cv.fit(X_train, y_train)
xgb_cv_results = pd.DataFrame(xgb_cv.cv_results_)

# to get only the parameters columns
col_results = column_results = ['param_{}'.format(name) for name in xgb_param_dists.keys()]

# add score columns
column_results += ['mean_test_score', 'std_test_score', 'rank_test_score']

xgb_cv_results = xgb_cv_results[column_results].sort_values(
                            'mean_test_score', ascending=False)

xgb_cv_results.columns = [shorten_param_names(col) for col in xgb_cv_results.columns]
xgb_cv_results<jupyter_output><empty_output><jupyter_text>Top XgBoost model has __0.7826 mean AUC and 0.0047 standard deviation in AUC__. It is significantly better than Logistic Regression model with 0.7643.   As we can see, XgBoost took 14+ hours to tune.   LightGBM   Similar to XgBoost, LightGBM has more hyperparameters to tune. Hence, we will use __50 iterations__ in RandomSearch cross-validation. LightGBM is faster on multicore CPUs and enabling GPU support for LightGBM was challenging and unsuccessful.<jupyter_code># store original column names as LightGBM does not allow special characters in column names.
import re

X_train = X_train.rename(columns = lambda x:re.sub('[^A-Za-z0-9_]+', '', x))

X_col_names = X_train.columns

X_train = X_train.rename(columns = lambda x:re.sub('[^A-Za-z0-9_]+', '', x))

lgb_clf = lgb.LGBMClassifier(objective= 'binary',
                            random_state =random_seed)


lgb_param_dists = {
            'learning_rate': loguniform(0.01, 1),
            'n_estimators': [200,400,800,1600,6400],
            'num_leaves': [4,8,16,50],
            'min_child_weight' : [1,4,8,16],
            'colsample_bytree' : [0.6,0.8],
            'subsample' : [0.6,0.8,1],
            'reg_alpha': loguniform(0.001, 1)}
%%time
lgb_cv = RandomizedSearchCV(lgb_clf, 
                           param_distributions=lgb_param_dists, 
                           scoring = 'roc_auc',
                           n_iter=lgb_iterations, 
                           cv=cv, 
                           n_jobs= -2,
                           random_state = random_seed)

lgb_cv.fit(X_train, y_train)
lgb_cv_results = pd.DataFrame(lgb_cv.cv_results_)

# to get only the parameters columns
col_results = column_results = ['param_{}'.format(name) for name in lgb_param_dists.keys()]

# add score columns
column_results += ['mean_test_score', 'std_test_score', 'rank_test_score']

lgb_cv_results = lgb_cv_results[column_results].sort_values(
                            'mean_test_score', ascending=False)

lgb_cv_results.columns = [shorten_param_names(col) for col in lgb_cv_results.columns]
lgb_cv_results<jupyter_output><empty_output><jupyter_text>Top LightGBM model has __0.7814 mean AUC and 0.0049 standard deviation in AUC__. It is slightly less than top XgBoost model. __Both the scores are within 1 standard deviation of each other.__LightGBM took 1 hour 50 minutes to tune. LightGBM with 5 CPU cores is 7 times faster than XgBoost and 6GB GPU.   Feature Importance :Let's look at feature importance plot of XgBoost & LightGBM models. We will use Gains to calculate & compare feature importance.<jupyter_code>best_xgboost_model = xgb_cv.best_estimator_
df_xgboost_feature_imp = pd.DataFrame(X_train.columns,columns=['feature'])
df_xgboost_feature_imp['importance_gains'] = best_xgboost_model.feature_importances_

df_xgboost_feature_imp.sort_values(by='importance_gains', ascending=False, inplace=True)
df_xgboost_feature_imp_top = df_xgboost_feature_imp.head(20)

fig, ax = plt.subplots(figsize=(8,12))
sns.barplot(x ='importance_gains', y='feature', data = df_xgboost_feature_imp_top)
plt.title('XgBoost- Feature Importance(Gains) - Top 20')
plt.show()
best_lightgbm_model = lgb_cv.best_estimator_
df_feature_imp = pd.DataFrame(X_train.columns,columns=['feature'])
df_feature_imp['importance_gains'] = best_lightgbm_model.booster_.feature_importance(importance_type='gain')
df_feature_imp['importance_split'] = best_lightgbm_model.booster_.feature_importance(importance_type='split')

df_feature_imp.sort_values(by='importance_gains', ascending=False, inplace=True)
df_featture_imp_top = df_feature_imp.head(20)

fig, ax = plt.subplots(figsize=(8,12))
sns.barplot(x ='importance_gains', y='feature', data = df_featture_imp_top)
plt.title('LightGBM- Feature Importance(Gains) - Top 20')
plt.show()
# Common features in top 20
set(df_xgboost_feature_imp_top.feature) & set(df_featture_imp_top.feature)<jupyter_output><empty_output><jupyter_text>External feature 2 & 3 are top two features in both the models. It is assuring that we observed trends in many of the top 20 features during EDA. Derived __Ratio features__ are coming out to be important predictors.   ComparisonWe have 8 common features in top 20 between XgBoost & LightGBM. XgBoost has more categorical features in top 20 whereas LightGBM has more continuous features in top 20. This may need further investigation.<jupyter_code>%%time
y_train_pred_probs = best_xgboost_model.predict_proba(X_train)[:,1]
%%time
y_train_pred_probs = best_lightgbm_model.predict_proba(X_train)[:,1]<jupyter_output>Wall time: 36.8 s<jupyter_text>XgBoost prediction is 1.8 times faster than LightGBM prediction.Hence, we will select __XgBoost  as the best model__ because1. slightly higher accuracy (0.7826 vs 0.7814)2. faster prediction time (20.4 s vs 36.8 s) Let us analyize features with null importance.<jupyter_code>null_imp_features = df_xgboost_feature_imp.loc[df_xgboost_feature_imp.importance_gains ==0, 'feature'].values
print('Number of Features with NULL Importance : {} '.format(len(null_imp_features)))
print("Features : ")
null_imp_features<jupyter_output>Number of Features with NULL Importance : 111 
Features :<jupyter_text>Organization type, last loan status from balance bureau, payment type of last loans and some flags features seem to be not at all importance. Let's remove these features and cross-validate the XgBoost model.<jupyter_code>X_train_dimR = X_train.drop(columns = list(null_imp_features))
print(X_train_dimR.shape)
best_model_DimR = xgb.XGBClassifier()
best_model_DimR.set_params(**best_xgboost_model.get_params())

cv_dimR = cross_validate(best_model_DimR, X_train_dimR, y_train, cv = cv, scoring = 'roc_auc')
test_scores = cv_dimR["test_score"]
print('AUC score mean :  {:.4f} std : {:.4f} '.format(test_scores.mean(), test_scores.std()))<jupyter_output>AUC score mean :  0.7827 std : 0.0052<jupyter_text>We have almost same score after removing these null importance features. But as we have noticed from feature importance camparison, ranking can somewhat differ by algorithms even when the algorithm performances are close. Assumig 423 features are manageable, let us not remove these features untill further investigation.  Prediction on Test Data<jupyter_code># changing column names for X_test as we removed special characters from column names for X_train while training LightGBM
X_test.columns = X_train.columns
y_pred_probs = best_xgboost_model.predict_proba(X_test)[:,1]
roc_auc_score(y_test.values, y_pred_probs)<jupyter_output><empty_output><jupyter_text>AUC of __0.7841__ is slightly higher than __0.7814__ mean AUC of cross-validation but well within the range __(0.7861, 0.7767)__ of 1 standard deviation. XgBoost model has well generalized on test set.   Let us look at model performance in detail using ROC curve and confusion matrix.<jupyter_code>fpr, tpr, thresholds = roc_curve(y_test.values, y_pred_probs)
# To get the best threshold, we will use J statistics
# source : https://machinelearningmastery.com/threshold-moving-for-imbalanced-classification/
# https://en.wikipedia.org/wiki/Youden%27s_J_statistic
j_stat = tpr - fpr
ix = np.argmax(j_stat)
best_thresh = thresholds[ix]
print('Best Threshold=%f' % (best_thresh))<jupyter_output>Best Threshold=0.073093<jupyter_text>J statistics achieves optimum balance between tpr & fpr by giving equal weight to false positive and false negative pay-offs.<jupyter_code># plot the roc curve for the model
fig, ax = plt.subplots(figsize=(6,6))
plt.plot([0,1], [0,1], linestyle='--', label='Chance', color='k', alpha= 0.5)
plt.plot(fpr, tpr, linestyle='-',label='Light gbm', color='brown', alpha= 0.7)
plt.scatter(fpr[ix], tpr[ix], marker='o', s=200, color='orange', label='Best Threshold')
# axis labels
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.legend()
# show the plot
plt.show()
# Corresponding fpr & tpr
print('Best Threshold: {}'.format(round(best_thresh,4)))
print('Corresponding fpr : {} and tpr: {}'.format(round(fpr[ix],4), round(tpr[ix],4)))
def get_confusion_matrix(y_probabilities, y_true, threshold):
    """
    Get the confusion matrix using specified threshold
    
    Keyword arguments:
    y_probabilities -- predicted probabilities of positive class
    y_true -- true value of y
    threshold -- threshold to use to get 0/1 prediction
    
    Returns :
    2*2 confusion_matrix
    """
    y_pred = np.zeros(y_true.shape)
    y_pred = np.where( y_probabilities > threshold ,1, 0 )
    return confusion_matrix(y_test, y_pred)

conf_matx = get_confusion_matrix(y_pred_probs, y_test, best_thresh)
def custom_plot_confusion_matrix(confusion_matrix_array):
    """
    Plot the confusion matrix
    
    Keyword arguments:
    confusion_matrix -- 2*2 confusion matrix
    """
    
    fig, ax = plt.subplots(1,2, figsize=(16,4))
    display_labels = ['Paid', 'Default']
    sns.heatmap(confusion_matrix_array, xticklabels = display_labels, yticklabels = display_labels, 
                annot=True, cmap='Greens', fmt='d', annot_kws={"size": 15}, ax=ax[0])
    ax[0].set_ylabel('True label')
    ax[0].set_xlabel('Predicted label')
    ax[0].set_title('Confusion Matrix')
    
    confusion_matrix_norm = confusion_matrix_array/confusion_matrix_array.sum(axis=1).reshape(-1,1)
    sns.heatmap(confusion_matrix_norm, xticklabels = display_labels, yticklabels = display_labels, 
                annot=True, cmap='Greens', fmt=".2%", annot_kws={"size": 15}, ax=ax[1])
    ax[1].set_ylabel('True label')
    ax[1].set_xlabel('Predicted label')
    ax[1].set_title('Confusion Matrix - Normalized')
    
    
custom_plot_confusion_matrix(conf_matx)<jupyter_output><empty_output><jupyter_text>Normalization is done as a % of true labels. We are able to correctly identify 73.23% of truely default loans and 26.77% of the default loans are missed with the selected thresold. Also, we are missing out in lending to 30.44% of truely payable loans by identifying them as potentially default loans. We are giving loans to 60981(58987+1994) applications which is 66.10% of total applications. Let us extend this analysis by incorperating pay-offs associated with each class.   Business PerformanceAssume that HomeCredit wants to __maximize profit__ as a business objective. We have following parameters from business.1. __5% of credit as a final profit__ on loans which are fully paid - this is the average % profit after removing commission, expenses, cost of capital etc.2. __40% of credit as a loss__ on loans which are defaulted - this is the average % loss after considering outstanding amount, recoveries etc. 90% of loans are cash loans and losses tend to be higher on cash loans. With this in mind, let us optimize the threshold for maximum profit and calculate confusion matrix. Now we will1. calculate average loan size for paid and default loans2. create a class for business parameters3. create a function to calculate expected profit at given threshold4. simulate different scenarios and plot profitability curve.<jupyter_code>mean_credit_paid = df_test.loc[df_test.target==0, ].amt_credit.mean()
mean_credit_default = df_test.loc[df_test.target==1, ].amt_credit.mean()
mean_credit_paid,mean_credit_default
profit_on_paid_loans = 0.05
loss_on_default_loans = 0.40

class cls_business_params:
    """
    Class to store business information regarding loan size, profit/loss
    set_business_params - to set the values of business parameters
    get_business_params - to get the values of business parameters
    """

    def __init__(self,mean_credit_paid,mean_credit_default,profit_on_paid_loans,loss_on_default_loans): 
        # instance as passed as a first agrument. Referred as self
        self.mean_credit_paid = mean_credit_paid
        self.mean_credit_default = mean_credit_default
        self.profit_on_paid_loans = profit_on_paid_loans
        self.loss_on_default_loans = loss_on_default_loans
        
    def set_business_params(mean_credit_paid,mean_credit_default,profit_on_paid_loans,loss_on_default_loans):
        self.mean_credit_paid = mean_credit_paid
        self.mean_credit_default = mean_credit_default
        self.profit_on_paid_loans = profit_on_paid_loans
        self.loss_on_default_loans = loss_on_default_loans        
    
    def get_business_params(self):
        return mean_credit_paid, mean_credit_default, profit_on_paid_loans, loss_on_default_loans
    
business_params = cls_business_params(mean_credit_paid,mean_credit_default,profit_on_paid_loans,loss_on_default_loans)
def calculate_expectancy(confusion_matrix, business_params):
    """
    Calculates  business expectancy/profit 
    
    Keyword arguments:
    confusion_matrix -- to set the values of business parameters
    business_params -- instance of  business_params containing profit,/loss, loan size values
    Returns :
    expectancy -- profit/loss based on confusion matrix and business parameters
    """
    tn, fp, fn, tp = confusion_matrix.flatten()
    mean_credit_paid, mean_credit_default, profit_on_paid_loans, loss_on_default_loans = business_params.get_business_params()

    expectancy = (tn * mean_credit_paid * profit_on_paid_loans         # profit from fully paid loans
                + tp * mean_credit_default * loss_on_default_loans     # notional profit from avoiding defaults 
                - fp * mean_credit_paid * profit_on_paid_loans         # notional loss from avoiding potentially paid back loans
                - fn * mean_credit_default * loss_on_default_loans)    # loss from deafults on loans

    return expectancy
def plot_profitability_curve(y_probabilities, y_true, business_params):
    """
    Plots profit against thresholds for positive class
    
    Keyword arguments:
    y_probabilities -- predicted probabilities of positive class
    y_true -- true value of y
    business_params -- instance of  business_params containing profit,/loss, loan size values
    Returns :
    max_profit -- maximum profit achievable with business constraints and model capabilities
    max_profit_threshold -- positive class threshold associated with maximum profit
    """

    thresholds = np.linspace(0,1,30)
    expectancy = []
    for threshold in thresholds:
        conf_matx = get_confusion_matrix(y_probabilities, y_true, threshold)
        expt = calculate_expectancy(conf_matx, business_params)
        expectancy.append(expt)
    
    idx = np.argmax(expectancy)
    max_profit = np.max(expectancy)
    max_profit_threshold = thresholds[idx]

    fig, ax = plt.subplots(figsize=(7,5))
    plt.plot(thresholds, expectancy, linestyle='-',label='profit', color='brown', alpha= 0.7)
    #plt.scatter(max_profit_threshold, max_profit, marker='o', s=200, color='orange', label='Best Threshold')
    max_profit_fmt = 'Max profit: {0:.1f} billion'.format(max_profit/1e9)
    max_profit_fmt  
    ax.annotate(max_profit_fmt, xy=(max_profit_threshold, max_profit),  xycoords='data',
            xytext=(0.6, 0.92), textcoords='axes fraction',
            arrowprops=dict(facecolor='black', shrink=0.05),
            horizontalalignment='right', verticalalignment='top',
            fontsize=13)
    
    plt.xlabel('Thresholds')
    plt.ylabel('Profits')
    plt.legend()
    plt.title('Profitability Curve')
    plt.show()
    
    return max_profit, max_profit_threshold
max_profit, max_profit_threshold = plot_profitability_curve(y_pred_probs, y_test, business_params)
print('Profit Maximizing Threshold: {}'.format(round(max_profit_threshold,4)))
print('Corresponding Profit : {}'.format(round(max_profit)))<jupyter_output><empty_output><jupyter_text>We see that maximum profit of 1.9 billion(test data) will be achived at 0.1034 threshold. This threshold is higher than J statistics based threshold. Let us plot confusion matrix using this threshold.<jupyter_code>max_profit_conf_matx = get_confusion_matrix(y_pred_probs, y_test, max_profit_threshold)
custom_plot_confusion_matrix(max_profit_conf_matx)<jupyter_output><empty_output><jupyter_text>For profit maximization, we relax the threshold and extend the loans to 71147 (77.12%) applications. This is 11 % points or 10166 applications higher than previous analysis. Given the pay-offs, we can accommodate 935 (2929-1994) additional defaults  which are more than compensated for by approving additional 9231 (68218-58987) payable loans.<jupyter_code>tpr_profit = max_profit_conf_matx[1][1] / (max_profit_conf_matx[1][0] + max_profit_conf_matx[1][1])
fpr_profit = max_profit_conf_matx[0][1] / (max_profit_conf_matx[0][0] + max_profit_conf_matx[0][1])
print('Max profit fpr : {} and tpr: {}'.format(round(fpr_profit,4), round(tpr_profit,4)))<jupyter_output>Max profit fpr : 0.1956 and tpr: 0.6067<jupyter_text>Profit %Assuming that we will accept all the applications that are predicted as paid and reject all the application that are predicted as default, let us calculate the expected booked profit %.<jupyter_code>mean_credit_paid, mean_credit_default, profit_on_paid_loans, loss_on_default_loans = business_params.get_business_params()

potential_credit_given = (max_profit_conf_matx[0][0] * mean_credit_paid        # loan given and paid back
                         + max_profit_conf_matx[1][0] * mean_credit_default    # loan given and defaulted
                         )

book_profit = (max_profit_conf_matx[0][0] * mean_credit_paid * profit_on_paid_loans        # profit on paid loans
              - max_profit_conf_matx[1][0] * mean_credit_default * loss_on_default_loans)  # loss on default loans

print('*****   Expected booked profit : {:.2%}   *****'.format(book_profit/potential_credit_given))<jupyter_output>*****   Expected booked profit : 3.29%   *****<jupyter_text>With this recommendation, we conclude HomeCredit default prediction project.  Next Steps  Following ideas can be explored in the future  1. Modelling applicants with and without credit card separately - We have just about 28% applicants who had credit cards with HomeCredit but the records are credible(65k). This leads to missing values in most of the records. Also, applicants without credit cards have a 1.2% points higher default rate. Separating these two groups may lead to better accuracy.2. Baysian Hyperparameter Optimization - We can explore Bayesian Hyperparameter Optimization which is a targeted search strategy based on model improvements. This may lead to improved tuning at a lesser time.3. Feature Importance - We had seen that feature importance can differ by algorithms with almost same performance. This needs further investigation.4. Model Interpretation - HomeCredit is more focused on prediction. But this is a good to have feature where we can demonstrate impact of the individual feature for easy buy-in of the model by business.<jupyter_code>import os
import pickle
import datetime
from sklearn import __version__ as sklearn_version


final_model = best_xgboost_model
final_model.version = '1.0'
final_model.pandas_version = pd.__version__
final_model.numpy_version = np.__version__
final_model.sklearn_version = sklearn_version
final_model.xgboost_version = xgb.__version__
final_model.X_columns = [col for col in X_train.columns]
final_model.build_datetime = datetime.datetime.now()

modelpath = '../models'
model_name = 'home_credit_loan_default_model.pkl'
fpath = os.path.join(modelpath, model_name)

with open(fpath, 'wb') as f:
            pickle.dump(final_model, f)

xgb_tune_results = '../reports/model_tunning/xgb_tune_results.csv'
xgb_cv_results.to_csv(xgb_tune_results)

lgb_tune_results = '../reports/model_tunning/lgb_tune_results.csv'
lgb_cv_results.to_csv(lgb_tune_results)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08278285770069312, 'fraction_numerical': 0.017767075591422342, 'mean_word_length': 4.50259623992838, 'pattern_counts': {'":': 2, '<': 82, '<?xml version=': 0, '>': 83, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '554957', 'license': 'MIT', 'n_tokens_mistral': 9907, 'n_tokens_neox': 9366, 'n_words': 2653, 'path': 'notebooks/04_Modelling.ipynb', 'repo_name': 'rohanML/home_credit_default_risk', 'size': '30731'}