{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport seaborn as sns\nimport xgboost as xgb\nimport matplotlib.pyplot as plt\nfrom sklearn.metrics import roc_auc_score, f1_score, roc_curve, confusion_matrix, accuracy_score\nfrom sklearn.model_selection import train_test_split, GridSearchCV \nfrom sklearn.linear_model import  LogisticRegression\nfrom sklearn.utils import resample\npd.set_option('display.max_columns', None)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We begin by importing our data. As the number of indepedent varibales is large, i.e. 66, we will examine if there is correlation among them. If we find that two predictors have correaltion larger in 0.7 in absolute value, we will drop one of them. In this way we will reduce the dimensionality of the problem. ","metadata":{}},{"cell_type":"code","source":"df = pd.read_csv('train.csv')\ndf = df.apply(pd.to_numeric, errors='coerce', downcast='integer')\nprint('Initial number of predictors is 66')\nX, y = df.iloc[:,1:-1], df.iloc[:,-1]\ncorr_matrix = X.corr().abs()\n\n# Select upper triangle of correlation matrix\nupper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(bool))\nto_drop = [column for column in upper.columns if any(upper[column] > 0.7)]\n\n# Drop features \nX.drop(to_drop, axis=1, inplace=True)\nX, y = X.to_numpy(), y.to_numpy()\nX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.3, random_state=123)\nprint('Number of predictors after dropping the higly correlated ones is', X_train.shape[1])","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We see that almost half of them were dropped. Let us proceed by checking what is the probability of bankruptcy in our sample. ","metadata":{}},{"cell_type":"code","source":"sns.countplot(x = y_train, palette = 'Set3')\nplt.show()\np_estimate = np.mean(y_train)\nprint('The sample estimator for the probability of going bunkrupt is', p_estimate)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We see that it is extremely low, suggesting that we have higly inbalanced data, with the majority class being the banks that did not go bankrupt.This creates a problem, as the minority class instances contribute little to the empirical risk and consequently may be ignored during optimization. To treat this, we upsample the minority class and arrive at balanced data. \n\nThen we verify that the two classes have the same number of labels in the training set.","metadata":{}},{"cell_type":"code","source":"X_train_0 = X_train[y_train==0]\ny_train_0 = y_train[y_train==0]\ny_train_0 = y_train_0.reshape(y_train_0.shape[0],1)\nX_train_1 = X_train[y_train==1]\ny_train_1 = y_train[y_train==1]\ny_train_1 = y_train_1.reshape(y_train_1.shape[0],1)\n\n\nX_0 = np.concatenate((X_train_0,y_train_0),axis=1)\nX_1 = np.concatenate((X_train_1,y_train_1),axis=1)\n\n\n\nupsampled_one = resample(X_1,\n                      replace=True, # sample with replacement\n                      n_samples=X_0.shape[0],\n                      random_state=3)\n                         # match number in majority class)\nupsampled_one = np.array(upsampled_one)\n\n\ny_train = np.concatenate((upsampled_one[:,-1].reshape(X_0.shape[0]),y_train_0.reshape(X_0.shape[0])),axis=0)\ny_train = y_train.astype('int')\nX_train = np.concatenate((upsampled_one[:,:-1],X_train_0),axis=0)\nnr_pos = np.sum(y_train)\nnr_neg = y_train.shape[0]-nr_pos\nprint('The ratio of negative and positive examples is', nr_pos/nr_neg)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We now proceed with creating a grid for the parameters of the gradient boosted tree classifier and then use 3-fold crossvalidation to find the optimal values of this grid. ","metadata":{}},{"cell_type":"code","source":"paramGrid = {\n    'max_depth': range (2, 10, 2),\n    'alpha': range(0, 3, 10),\n    'subsample' : [0.6,  1.0],\n    'learning_rate' : [0.1, 0.01, 0.05],\n    'n_estimators' : range(100, 200, 40),\n    'colsample_bytree' : [0.8, 1.0],\n    'colsample_bylevel': [0.8, 1.0],\n    'scale_pos_weight': [nr_neg/nr_pos],\n    'objective': ['binary:logistic'],\n    'eval_metric': ['logloss']\n}\n\nxgb_clf = xgb.XGBClassifier(tree_method='hist')\nxgb_rs = GridSearchCV(estimator=xgb_clf, param_grid=paramGrid, cv=5, verbose=1, scoring='f1')\nxgb_rs.fit(X_train, y_train)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We pick the best estimator and check how the corresponding Receiver Operator Curve (ROC) looks like. \n\nAs a reminder, the ROC curve is shaped from the True Positive Rate (TPR=TP/(TP+FN)) and the False Positive Rate (FPR=FP/(FP+TN)), when altering the threshold of probability for classifying an observation as positive or negative.  \n(FP=False Positive, TN=True Negative, FN=False Negative)","metadata":{}},{"cell_type":"code","source":"best = xgb_rs.best_estimator_\ny_pred = best.predict(X_test)\ny_pred_scores = best.predict_proba(X_test)[:,1]\nauc = roc_auc_score(y_test, best.predict(X_test))\nfpr, tpr, thresholds = roc_curve(y_test, y_pred_scores)\nplt.figure()\nlw = 2\nplt.plot(\n    fpr,\n    tpr,\n    color=\"darkorange\",\n    lw=lw,\n    label=\"ROC curve (area = %0.2f)\" %  auc,\n)\nplt.plot([0, 1], [0, 1], color=\"navy\", lw=lw, linestyle=\"--\")\nplt.xlim([0.0, 1.0])\nplt.ylim([0.0, 1.05])\nplt.xlabel(\"False Positive Rate\")\nplt.ylabel(\"True Positive Rate\")\nplt.title(\"Receiver operating characteristic example\")\nplt.legend(loc=\"lower right\")\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We see that the curve gives good performance and the area under it is 0.71. This value ranges from 0 to 1, with higher values being better. Finally, to test our classifer, let us do a confusion matrix and compute its test accuracy:","metadata":{}},{"cell_type":"code","source":"def plot_confusion_matrix(cm, classes=None, title='Confusion matrix'):\n    \"\"\"Plots a confusion matrix.\"\"\"\n    if classes is not None:\n        sns.heatmap(cm, cmap=\"YlGnBu\", xticklabels=classes, yticklabels=classes, vmin=0., vmax=1., annot=True, annot_kws={'size':50})\n    else:\n        sns.heatmap(cm, vmin=0., vmax=1.)\n    plt.title(title)\n    plt.ylabel('True label')\n    plt.xlabel('Predicted label')\n# Visualizing cm\n\ncm = confusion_matrix(y_test, best.predict(X_test))\ncm_norm = cm / cm.sum(axis=1).reshape(-1,1)\n\nplot_confusion_matrix(cm_norm, classes = best.classes_, title='Confusion matrix')\nprint('Test accuracy of the model is', accuracy_score(y_test, best.predict(X_test)))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"While the model still could improved, since it has a 57% chance of predicting a postive example as negative, it performs well. The test accuracy is 97%. Finally, we compute the predictions for the test set of the Kaggle competition and create the CSV file to upload it. ","metadata":{}},{"cell_type":"code","source":"df = pd.read_csv('test.csv')\ndf = df.apply(pd.to_numeric, errors='coerce', downcast='integer')\nX_new = df.iloc[:,1:]\nX_new.drop(to_drop, axis=1, inplace=True)\nX_new = X_new.to_numpy()\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub = pd.DataFrame()\ny_new = best.predict(X_new)\nsub['id'] = df['id']\nsub['class'] = y_new\nsub.to_csv('sub.csv', index=False)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}