{"cells":[{"metadata":{},"cell_type":"markdown","source":"[Elo](https://elo.com.br/) is a Brazillian debit and credit card brand.  They offer credit and prepaid transactions, and have paired up with merchants in order offer promotions to cardholders.  In order to offer more relevant and personalized promotions, in a [recent Kaggle competition](https://www.kaggle.com/c/elo-merchant-category-recommendation), Elo challenged Kagglers to predict customer loyalty based on transaction history.  Presumably they plan to use a loyalty-predicting model in order to determine what promotions to offer to customers based on how certain offers are predicted to effect card owners' card loyalty.\n\nIn previous kernels, we performed some [exploratory data analysis](http://www.kaggle.com/brendanhasz/elo-data-cleaning-and-eda), and then did some [feature engineering and selection](http://www.kaggle.com/brendanhasz/elo-feature-engineering-and-feature-selection). In this kernel, we'll build models to predict the customer loyalty given the features we engineered.  First we'll run some baseline models, then create some models which make use of outlier prediction, and finally construct ensemble and stacking models.\n\n**Outline**\n\n- [Loading the Data](#loading-the-data)\n- [Evaluation Metric](#evaluation-metric)\n- [Baseline](#baseline)\n- [Gradient-boosted decision tree regression](#gradient-boosted-decision-tree-regression)\n- [Predictions without outliers](#predictions-without-outliers)\n- [Outlier prediction model](#outlier-prediction-model)\n- [Hyperparameter Optimization](#hyperparameter-optimization)\n- [Performance of Individual Models](#performance-of-individual-models)\n- [Ensemble Model](#ensemble-model)\n- [Stacking Model](#stacking-model)"},{"metadata":{"trusted":true,"_uuid":"a24c2b06292c72efd602b460c800d7bcfec7e087"},"cell_type":"code","source":"import gc\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.impute import SimpleImputer\nfrom sklearn.preprocessing import RobustScaler\nfrom sklearn.metrics import make_scorer\nfrom sklearn.model_selection import cross_val_score\nfrom sklearn.model_selection import KFold\nfrom sklearn.linear_model import BayesianRidge\n\nfrom catboost import CatBoostRegressor\nfrom xgboost import XGBClassifier, XGBRegressor\nfrom lightgbm import LGBMRegressor\n\n!pip install git+http://github.com/brendanhasz/dsutils.git\nfrom dsutils.encoding import TargetEncoderCV\nfrom dsutils.optimization import optimize_params\nfrom dsutils.ensembling import EnsembleRegressor, StackedRegressor\n\n# Plot settings\n%matplotlib inline\n%config InlineBackend.figure_format = 'svg'\nsns.set()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0f49e76fcb13fe760a037b14babd877d2707a1be"},"cell_type":"markdown","source":"<a id='loading-the-data'></a>\n## Loading the Data\n\nFirst, let's load the data which includes the features we engineered in the [last kernel](http://www.kaggle.com/brendanhasz/elo-feature-engineering-and-feature-selection)."},{"metadata":{"trusted":true,"_uuid":"2818cb72b6c4fc94becf5eaf64419b9fc6d9c7b3"},"cell_type":"code","source":"# Load data containing all the features\nfname = '../input/elo-feature-selection/card_features_top100.feather'\ncards = pd.read_feather(fname)\ncards.set_index('card_id', inplace=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Then, we can split the dataset into test vs train, and features (X) vs the target (y)."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Test indexes\ntest_ix = cards['target'].isnull()\n\n# Test data\nX_test = cards[test_ix]#TODO: .copy()\ndel X_test['target']\n\n# Training data\ncards_train = cards[~test_ix]#TODO: .copy()\ny_train = cards_train['target']#TODO: .copy()\nX_train = cards_train.copy()\ndel X_train['target']\n\n# Clean up\ndel cards_train\ndel cards\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"812d71c29dd9940d068510a1ebd2a3283b300a94"},"cell_type":"markdown","source":"<a id='evaluation-metric'></a>\n## Evaluation Metric\n\nWe'll be using the root mean squared error as our evaluation metric:\n\n$$\nRMSE(y, \\hat{y}) = \\sqrt{ \\frac{1}{N} \\sum_{i=1}^N (y_i - \\hat{y}_i)^2 }\n$$\n\nUnfortunately, sklearn doesn't have a built-in RMSE metric function, so we'll build our own."},{"metadata":{"trusted":true,"_uuid":"db8ad5d26726d3290d4f638b6b30f767686aeccc"},"cell_type":"code","source":"def root_mean_squared_error(y_true, y_pred):\n    \"\"\"Root mean squared error regression loss\"\"\"\n    return np.sqrt(np.mean(np.square(y_true-y_pred)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And we'll also create a sklearn scorer corresponding to that metric:"},{"metadata":{"trusted":true},"cell_type":"code","source":"rmse_scorer = make_scorer(root_mean_squared_error)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"848e49db4c208d155339a0be381d7ef2efb97626"},"cell_type":"markdown","source":"<a id='baseline'></a>\n## Baseline\n\nBefore we start creating predictive models, it's good to have a baseline with which to compare our models' performance.  That way, if our models perform worse than the baseline, we'll know something's broken somewhere!  \n\nOur predictive models should *definitely* be able to do better than if we just predict the mean loyalty score for each sample.  What RMSE can we acheive if we just predict the mean?"},{"metadata":{"trusted":true,"_uuid":"24e7b013b3b42bee0d3a30298b888fd7c4b46422"},"cell_type":"code","source":"root_mean_squared_error(np.mean(y_train), y_train)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"OK, so our models should *for sure* be getting RMSE values lower than 3.85."},{"metadata":{"_uuid":"14729d532ea68b83851b99059cb4b47278e2ee18"},"cell_type":"markdown","source":"<a id='gradient-boosted-decision-tree-regression'></a>\n## Gradient-boosted decision tree regression\n\nAs a slightly more realistic baseline, let's first just use [CatBoost](http://catboost.ai) by itself, without any parameter tuning or anything fancy.  \n\nWell, one fancy thing: we'll also target-encode the categorical columns.  The only categorical columns in the dataset were created by taking the most common category during an aggregation, and so they all have 'mode' in their column name."},{"metadata":{"trusted":true,"_uuid":"1202ba7547b5284e11550cf6a98145420ec340d8"},"cell_type":"code","source":"%%time\n\n# Categorical columns\ncat_cols = [c for c in X_train if 'mode' in c] \n\n# Regression pipeline\nmodel = Pipeline([\n    ('targ_enc',  TargetEncoderCV(cols=cat_cols)),\n    ('scaler',    RobustScaler()),\n    ('imputer',   SimpleImputer(strategy='median')),\n    ('regressor', CatBoostRegressor(verbose=False))\n])\n\n# Cross-validated performance\nscores = cross_val_score(model, X_train, y_train, \n                         cv=3, scoring=rmse_scorer)\nprint('Cross-validated MSE: %0.3f +/- %0.3f'\n      % (scores.mean(), scores.std()))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So just with CatBoost, we can do better than just predicting the mean.  But not a whole lot better!  Let's see if we can improve our predictive performance by creating more complex models."},{"metadata":{"_uuid":"400d86a54a038dca5781f6d6211a091a607bdb33"},"cell_type":"markdown","source":"<a id='predictions-without-outliers'></a>\n## Predictions without outliers\n\nOne weird thing about the target is that there's a lot of outliers:"},{"metadata":{"trusted":true,"_uuid":"415155721134ffe255a96fdd1201471505296afe"},"cell_type":"code","source":"# Show histogram of target\ny_train.hist(bins=30)\nplt.xlabel('target')\nplt.ylabel('count')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"All the outliers have exactly the same value:"},{"metadata":{"trusted":true,"_uuid":"02e2e0b524e02bd6f822477229e3607a17f274d0"},"cell_type":"code","source":"y_train[y_train<-20].unique()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And around 1% of the samples have these outlier values:"},{"metadata":{"trusted":true,"_uuid":"23f800eda180cfb743d60f6ea04e2c9c432cc00a"},"cell_type":"code","source":"print('Percent of targets which are outliers:', 100*np.mean(y_train<-20))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"47a8ecef21e78003cd56b5f54426b4f87b7e45f6"},"cell_type":"markdown","source":"Those outliers are probably having an outsized effect on the training of the model.  It might be best to just train the model w/o including the outliers.  To try that, we'll build a function which lets us specify what samples to use for testing and training.  Then, we can use that function to compute the cross-validated performance of our regressor when we train using only the non-outlier data.\n"},{"metadata":{"trusted":true,"_uuid":"8f2a0557e9a9e8d683585724030395d13a2d0604"},"cell_type":"code","source":"def cross_val_metric(model, X, y, cv=3, \n                     metric=root_mean_squared_error, \n                     train_subset=None, test_subset=None, \n                     shuffle=False, display=None):\n    \"\"\"Compute a cross-validated metric for a model.\n    \n    Parameters\n    ----------\n    model : sklearn estimator or callable\n        Model to use for prediction.  Either an sklearn estimator (e.g. a \n        Pipeline), or a function which takes 3 arguments: \n        (X_train, y_train, X_test), and returns y_pred.  X_train and X_test\n        should be pandas DataFrames, and y_train and y_pred should be \n        pandas Series.\n    X : pandas DataFrame\n        Features.\n    y : pandas Series\n        Target variable.\n    cv : int\n        Number of cross-validation folds\n    metric : sklearn.metrics.Metric\n        Metric to evaluate.\n    train_subset : pandas Series (boolean)\n        Subset of the data to train on. \n        Must be same size as y, with same index as X and y.\n    test_subset : pandas Series (boolean)\n        Subset of the data to test on.  \n        Must be same size as y, with same index as X and y.\n    shuffle : bool\n        Whether to shuffle the data. Default = False\n    display : None or str\n        Whether to print the cross-validated metric.\n        If None, doesn't print.\n    \n    Returns\n    -------\n    metrics : list\n        List of metrics for each test fold (length cv)\n    preds : pandas Series\n        Cross-validated predictions\n    \"\"\"\n    \n    # Use all samples if not specified\n    if train_subset is None:\n        train_subset = y.copy()\n        train_subset[:] = True\n    if test_subset is None:\n        test_subset = y.copy()\n        test_subset[:] = True\n    \n    # Perform the cross-fold evaluation\n    metrics = []\n    TRix = y.copy()\n    TEix = y.copy()\n    all_preds = y.copy()\n    kf = KFold(n_splits=cv, shuffle=shuffle)\n    for train_ix, test_ix in kf.split(X):\n        \n        # Indexes for samples in training fold and in train_subset\n        TRix[:] = False\n        TRix.iloc[train_ix] = True\n        TRix = TRix & train_subset\n        \n        # Indexes for samples in test fold and in test_subset\n        TEix[:] = False\n        TEix.iloc[test_ix] = True\n        TEix = TEix & test_subset\n        \n        # Predict using a function\n        if callable(model):\n            preds = model(X.loc[TRix,:], y[TRix], X.loc[TEix,:])\n        else:\n            model.fit(X.loc[TRix,:], y[TRix])\n            preds = model.predict(X.loc[TEix,:])\n        \n        # Store metrics for this fold\n        metrics.append(metric(y[TEix], preds))\n        \n        # Store predictions for this fold\n        all_preds[TEix] = preds\n\n    # Print the metric\n    metrics = np.array(metrics)\n    if display is not None:\n        print('Cross-validated %s: %0.3f +/- %0.3f'\n              % (display, metrics.mean(), metrics.std()))\n        \n    # Return a list of metrics for each fold\n    return metrics, all_preds","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"49cd31bdfb98181123f4ef8834235b970f31268e"},"cell_type":"markdown","source":"Now, we can train the same model as before, but only on datapoints which are *not* outliers, and then test its cross-validated performance on all the data."},{"metadata":{"trusted":true,"_uuid":"fe5995f23a818297cee639a2280cc23a252a579b"},"cell_type":"code","source":"%%time\n\n# Compute which samples are outliers\nnonoutliers = y_train>-20\n\n# Cross-validated performance training only on non-outliers\ncross_val_metric(model, X_train, y_train, cv=3, \n                 metric=root_mean_squared_error, \n                 train_subset=nonoutliers,\n                 display='RMSE')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Hmm, that didn't seem to help at all.  The model still performs better than just guessing the mean (a RMSE of 3.85), but doesn't beat simply using CatBoost on all the data by a long shot (which results in a RMSE of 3.67)."},{"metadata":{"_uuid":"e01791212c89fbef0b575d4ba9567437256aad0f"},"cell_type":"markdown","source":"<a id='outlier-prediction-model'></a>\n## Outlier prediction model\n\nIt may work better to train one model to try and predict whether a sample will be an outlier or not, and then train another model on the non-outlier data.\n\nLet's create a classification model to classify outlier vs. non-outlier, and a regression model to estimate the values of non-outlier data.  And we'll use XGBoost for the classifier instead of CatBoost, because what the hey."},{"metadata":{"trusted":true,"_uuid":"4ffc3fc361b2a4d3ee06a1febbafb0d3c50e4c3e"},"cell_type":"code","source":"# Classification pipeline\nclassifier = Pipeline([\n    ('targ_enc',   TargetEncoderCV(cols=cat_cols)),\n    ('scaler',     RobustScaler()),\n    ('imputer',    SimpleImputer(strategy='median')),\n    ('classifier', XGBClassifier())\n])\n\n# Regression pipeline\nregressor = Pipeline([\n    ('targ_enc',  TargetEncoderCV(cols=cat_cols)),\n    ('scaler',    RobustScaler()),\n    ('imputer',   SimpleImputer(strategy='median')),\n    ('regressor', CatBoostRegressor())\n])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0ba9bbb515533b32f328a862776e95e7cda4f5a3"},"cell_type":"markdown","source":"Then, we can use the classification model to predict the probability that each test sample is an outlier, and the regression model to estimate the value for non-outlier values.  For samples which the classifier predicts are outliers, we'll just set the predicted value to be the outlier value (`-33.2`...).  For samples which the classifier predicts are *not* outliers, we'll use the regression model to predict the value.\n\nIn order to evaluate our model with cross-fold validation, we'll create a function which performs the whole (somewhat non-orthodox) pipeline."},{"metadata":{"trusted":true,"_uuid":"3b6362bdb2b436c062ea5361c773abba7bcb347c"},"cell_type":"code","source":"def classifier_regressor(X_tr, y_tr, X_te):\n    \"\"\"Classify outliers, and set to outlier value if a predicted outlier\n    else use a regressor trained on non-outliers\n    \n    Parameters\n    ----------\n    X_tr : pandas DataFrame\n        Training features\n    y_tr : pandas Series\n        Training target\n    X_te : pandas DataFrame\n        Test features\n        \n    Returns\n    -------\n    y_pred : pandas Series\n        Predicted y values for samples in X_te\n    \"\"\"\n    \n    # Fit the classifier to predict outliers\n    outliers = y_tr<-20\n    fit_classifier = classifier.fit(X_tr, outliers)\n    \n    # Samples which do not have outlier target values\n    X_nonoutlier = X_tr.loc[~outliers, :]\n    y_nonoutlier = y_tr[~outliers]\n\n    # Fit the regressor to estimate non-outlier targets\n    fit_regressor = regressor.fit(X_nonoutlier, y_nonoutlier)\n    \n    # Predict outlier probability\n    pred_outlier = fit_classifier.predict_proba(X_te)[:,1]\n\n    # Estimate target\n    y_pred = fit_regressor.predict(X_te)\n    \n    # Estimate number of outliers in test data\n    outlier_num = int(X_te.shape[0]*np.mean(outliers))\n\n    # Set that proportion of top estimated outliers to outlier value\n    thresh = np.sort(pred_outlier)[-outlier_num]\n    y_pred[pred_outlier>thresh] = -33.21928024\n    \n    # Return predictions\n    return y_pred","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b426c77642bb55b311e1d0b794e41692548c5594"},"cell_type":"markdown","source":"Now we can evaluate the model's performance:"},{"metadata":{"trusted":true,"_uuid":"eab6ee71b62e1df190255a9a31f6121763610397"},"cell_type":"code","source":"%%time\n\n# Performance of mixed classifier + regressor\ncross_val_metric(classifier_regressor, \n                 X_train, y_train,\n                 metric=root_mean_squared_error, \n                 cv=3, display='RMSE')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Yikes, that made things even worse than the baseline!  It might work better to *weight* the predictions by the outlier probability according to the classifier model (instead of using a binary threshold, as we did above).  However, at that point one may as well just use a regression model for all datapoints, so that's what we'll do."},{"metadata":{"_uuid":"7a806961f1301e73761b67438fd28ecbd1277d9f"},"cell_type":"markdown","source":"<a id='hyperparameter-optimization'></a>\n## Hyperparameter Optimization\n\nLater on, we'll be using an ensemble of multiple models to make predictions.  Before we do that, however, we should find the hyperparameters for each model which work best with this dataset.  To do that we'll use Bayesian hyperparameter optimization, which uses Gaussian processes to find the best set of parameters efficiently (see my post on [Bayesian hyperparameter optimization](http://brendanhasz.github.io/2019/03/28/hyperparameter-optimization.html))."},{"metadata":{},"cell_type":"markdown","source":"### Bayesian Ridge Regression\n\nThe first model we'll be using is a [Bayesian ridge regression](http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.BayesianRidge.html).  This model has several hyperparameters, including:\n\n* `alpha_1` and `alpha_2` - control the prior on the regression weight priors,\n* `lambda_1` and `lambda_2` - control the prior on the noise standard deviation prior.\n\nTo find the settings of these hyperparameters which allow the model to make the best predictions with this dataset, we first need to define a processing pipeline."},{"metadata":{"trusted":true},"cell_type":"code","source":"# XGBoost pipeline\nxgb_pipeline = Pipeline([\n    ('targ_enc',  TargetEncoderCV(cols=cat_cols)),\n    ('scaler',    RobustScaler()),\n    ('imputer',   SimpleImputer(strategy='median')),\n    ('regressor', BayesianRidge(alpha_1=1e-6,\n                                alpha_2=1e-6,\n                                lambda_1=1e-6,\n                                lambda_2=1e-6))\n])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Then, we'll define a dictionary containing the parameters whose values we want to optimize.  The keys contain a string with the element of the pipeline (in this case all our parameters are for the `regressor`), and the parameter of that pipeline element (the parameter name), separated by a double underscore.  The values of the bounds dictionary contain a list with three elements: the first element is the lower bound, the second element contains the upper bound, and the last element is the datatype of the parameter.  The optimizer will test parameter values between the lower and upper bounds which match the corresponding datatype."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Parameter bounds\nbounds = {\n    'regressor__alpha_1':  [1e-7, 1e-5, float],\n    'regressor__alpha_2':  [1e-7, 1e-7, float],\n    'regressor__lambda_1': [1e-7, 1e-7, float],\n    'regressor__lambda_2': [1e-7, 1e-7, float],\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Finally, we can use Bayesian optimization to find the values of these parameters which cause the model to generate the most accurate predictions.  Note that I've run the hyperparameter optimization for each model in separate runs of the kernel, because running them all would exceed Kaggle kernel's runtime limits!"},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"\n%%time\n# Find the optimal parameters\nopt_params, gpo = optimize_params(\n    X_train, y_train, xgb_pipeline, bounds,\n    n_splits=2, max_evals=300, n_random=100,\n    metric=root_mean_squared_error)\n    \n# Show the expected best parameters\nprint(opt_params)\n\"\"\"\n\n# Output:\n#    {'regressor__alpha_1': 2.76e-6,\n#     'regressor__alpha_2': 7.96e-7,\n#     'regressor__lambda_1': 7.32e-7,\n#     'regressor__lambda_2': 9.69e-6}\n#    Wall time: 5h 35min 34s","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ea5bf353a5525fb9cae7dba0211accefc3710eff"},"cell_type":"markdown","source":"### XGBoost\n\nThe second model we'll be using is [XGBoost](http://github.com/dmlc/xgboost), a gradient-boosted decision tree model.  This model has several [parameters](http://github.com/dmlc/xgboost/blob/master/doc/parameter.rst#learning-task-parameters) which control how the decision trees are constructed, including:\n\n* `max_depth` - controls the maximum number of splits in the decision trees\n* `learning_rate` - the learning rate used for boosting\n* `n_estimators` - how many decision trees to use\n\nSide note: the [default loss function for XGBoost is the squared error](http://github.com/dmlc/xgboost/blob/master/doc/parameter.rst#learning-task-parameters), unlike [LightGBM](http://lightgbm.readthedocs.io/en/latest/Parameters.html#core-parameters) and [CatBoost](http://catboost.ai/docs/concepts/python-reference_catboostregressor.html), which use a default loss (or objective) of the root mean squared error, which is our metric for this competition.  In theory we could create a custom objective function and pass that as the `objective` argument to XGBRegressor, but we'd need to compute the gradient and the hessian of the predictions with respect to the RMSE.  So: the reason XGBoost isn't performing quite as well as the other two doesn't necessarily mean it's a worse algorithm - it could just be because we're not using quite the correct loss function.  Using the MSE instead of the RMSE will cause outliers to have a larger effect on the fit.  And as we saw above, there's a reasonable amount of outliers.\n\nAfter defining a pipeline which uses XGBoost, and bounds for the parameters, we can again use Bayesian optimization to find good parameter choices."},{"metadata":{"trusted":true,"_uuid":"833f6b6399d5517bbd0ca9f28f988f4edd666ffa"},"cell_type":"code","source":"\"\"\"\n%%time\n\n# XGBoost pipeline\nxgb_pipeline = Pipeline([\n    ('targ_enc',  TargetEncoderCV(cols=cat_cols)),\n    ('scaler',    RobustScaler()),\n    ('imputer',   SimpleImputer(strategy='median')),\n    ('regressor', XGBRegressor(max_depth=3,\n                               learning_rate=0.1,\n                               n_estimators=100))\n])\n\n# Parameter bounds\nbounds = {\n    'regressor__max_depth': [2, 3, int],\n    'regressor__learning_rate': [0.01, 0.5, float],\n    'regressor__n_estimators': [50, 200, int],\n}\n\n# Find the optimal parameters\nopt_params, gpo = optimize_params(\n    X_train, y_train, xgb_pipeline, bounds,\n    n_splits=2, max_evals=100, n_random=20,\n    metric=root_mean_squared_error)\n    \n# Show the expected best parameters\nprint(opt_params)\n\"\"\"\n\n# Output:\n#    {'regressor__max_depth': 5, \n#     'regressor__learning_rate': 0.0739,\n#     'regressor__n_estimators': 158}\n#    Wall time: 5h 20min 3s","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"135fcc7c6003931135d63aa43ae57c29a4484b5b"},"cell_type":"markdown","source":"### LightGBM\n\nAnother tree-based model we'll use is Microsoft's [LightGBM](http://github.com/microsoft/LightGBM).  Its [parameters](http://lightgbm.readthedocs.io/en/latest/Parameters.html#core-parameters) include:\n\n* `num_leaves` - the maximum number of leaves in a tree\n* `max_depth` - controls the maximum number of splits in the decision trees\n* `learning_rate` - the learning rate used for boosting\n* `n_estimators` - how many decision trees to use\n\nAgain we'll use Bayesian optimization to find a good combination of parameters:"},{"metadata":{"trusted":true,"_uuid":"0820160411fdf01a856d53ff1cce09ae9c707387"},"cell_type":"code","source":"\"\"\"\n%%time\n\n# LightGBM pipeline\nlgbm_pipeline = Pipeline([\n    ('targ_enc',  TargetEncoderCV(cols=cat_cols)),\n    ('scaler',    RobustScaler()),\n    ('imputer',   SimpleImputer(strategy='median')),\n    ('regressor', LGBMRegressor(num_leaves=31,\n                                max_depth=-1,\n                                learning_rate=0.1,\n                                n_estimators=100))\n])\n\n# Parameter bounds\nbounds = {\n    'regressor__num_leaves': [10, 50, int],\n    'regressor__max_depth': [2, 8, int],\n    'regressor__learning_rate': [0.01, 0.2, float],\n    'regressor__n_estimators': [50, 200, int],\n}\n\n# Find the optimal parameters\nopt_params, gpo = optimize_params(\n    X_train, y_train, lgbm_pipeline, bounds,\n    n_splits=2, max_evals=100, n_random=20,\n    metric=root_mean_squared_error)\n    \n# Show the expected best parameters\nprint(opt_params)\n\"\"\"\n\n# Output\n#    {'regressor__num_leaves': 30,\n#     'regressor__max_depth': 6,\n#     'regressor__learning_rate': 0.0439,\n#     'regressor__n_estimators': 160}\n#    Wall time: 5h 19min 45s","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### CatBoost\n\nFinally, we'll also use Yandex's [CatBoost](https://github.com/catboost/catboost), whose [parameters](http://catboost.ai/docs/concepts/python-reference_parameters-list.html#python-reference_parameters-list) include:\n\n* `depth` - the depth of the decision trees\n* `learning_rate` - the learning rate used for boosting\n* `iterations` - the maximum number of decision trees to use\n* `l2_leaf_reg` - regularization parameter\n\nAgain let's use Bayesian optimization to find a good combination of parameters for the CatBoost model:"},{"metadata":{"trusted":true},"cell_type":"code","source":"\"\"\"\n%%time\n\n# CatBoost pipeline\ncb_pipeline = Pipeline([\n    ('targ_enc', TargetEncoderCV(cols=cat_cols)),\n    ('scaler', RobustScaler()),\n    ('imputer', SimpleImputer(strategy='median')),\n    ('regressor', CatBoostRegressor(depth=6,\n                                    learning_rate=0.03,\n                                    iterations=1000,\n                                    l2_leaf_reg=3.0))\n])\n\n# Parameter bounds\nbounds = {\n    'regressor__depth': [4, 10, int],\n    'regressor__learning_rate': [0.01, 0.1, float],\n    'regressor__iterations': [500, 1500, int],\n    'regressor__l2_leaf_reg': [1.0, 5.0, float],\n}\n\n# Find the optimal parameters\nopt_params, gpo = optimize_params(\n    X_train, y_train, cb_pipeline, bounds,\n    n_splits=2, max_evals=100, n_random=20,\n    metric=root_mean_squared_error)\n        \n# Show the expected best parameters\nprint(opt_params)\n\"\"\"\n\n# Output\n#    {'regressor__depth': 9, \n#     'regressor__learning_rate': 0.0240, \n#     'regressor__iterations': 650, \n#     'regressor__l2_leaf_reg': 1.47}\n#    Wall time: 5h 25min 22s","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<a id='performance-of-individual-models'></a>\n## Performance of Individual Models\n\nNow that we've found hyperparameter values which work well for each model, let's test the performance of each model individually before creating ensembles.  We'll define processing pipelines for each models, using the Bayesian-optimized hyperparameter values."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Bayesian ridge regression\nridge = Pipeline([\n    ('targ_enc',  TargetEncoderCV(cols=cat_cols)),\n    ('scaler',    RobustScaler()),\n    ('imputer',   SimpleImputer(strategy='mean')),\n    ('regressor', BayesianRidge(alpha_1=3e-6,\n                                alpha_2=1e-6,\n                                lambda_1=1e-6,\n                                lambda_2=1e-5))\n])\n\n# XGBoost\nxgboost = Pipeline([\n    ('targ_enc',  TargetEncoderCV(cols=cat_cols)),\n    ('scaler',    RobustScaler()),\n    ('imputer',   SimpleImputer(strategy='median')),\n    ('regressor', XGBRegressor(max_depth=5,\n                               learning_rate=0.07,\n                               n_estimators=150, \n                               n_jobs=2))\n])\n\n# LightGBM\nlgbm = Pipeline([\n    ('targ_enc',  TargetEncoderCV(cols=cat_cols)),\n    ('scaler',    RobustScaler()),\n    ('imputer',   SimpleImputer(strategy='mean')),\n    ('regressor', LGBMRegressor(num_leaves=30,\n                                max_depth=6,\n                                learning_rate=0.05,\n                                n_estimators=160))\n])\n\n# CatBoost\ncatboost = Pipeline([\n    ('targ_enc',  TargetEncoderCV(cols=cat_cols)),\n    ('scaler',    RobustScaler()),\n    ('imputer',   SimpleImputer(strategy='median')),\n    ('regressor', CatBoostRegressor(verbose=False, \n                                    depth=7,\n                                    learning_rate=0.02,\n                                    iterations=1000,\n                                    l2_leaf_reg=2.0))\n])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"How well does each model perform individually?"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\n# Bayesian Ridge Regression\nrmse_br, preds_br = cross_val_metric(\n    ridge, X_train, y_train,\n    metric=root_mean_squared_error,\n    cv=3, display='RMSE')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\n# XGBoost\nrmse_xgb, preds_xgb = cross_val_metric(\n    xgboost, X_train, y_train,\n    metric=root_mean_squared_error,\n    cv=3, display='RMSE')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\n# LightGBM\nrmse_lgb, preds_lgb = cross_val_metric(\n    lgbm, X_train, y_train,\n    metric=root_mean_squared_error,\n    cv=3, display='RMSE')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\n# CatBoost\nrmse_cb, preds_cb = cross_val_metric(\n    catboost, X_train, y_train,\n    metric=root_mean_squared_error,\n    cv=3, display='RMSE')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6ee4466ed153c3f559f3a22c7f04c1b0234b9706"},"cell_type":"markdown","source":"<a id='ensemble-model'></a>\n## Ensemble Model\n\nOften a way to get better predictions is to combine the predictions from multiple models, a technique known as [ensembling](https://en.wikipedia.org/wiki/Ensemble_learning).  This is because combining predictions from multiple models can reduce overfitting.  The simplest way to do this is to just average the predictions of several base learners (the models which are a part of the ensemble).\n\nHowever, ensembles tend to perform better when the predictions of their base learners aren't highly correlated.  How correlated are the predictions of our four individual models?\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Construct a DataFrame with each model's predictions\npdf = pd.DataFrame()\npdf['Ridge'] = preds_br\npdf['XGB'] = preds_xgb\npdf['LGBM'] = preds_lgb\npdf['CatBoost'] = preds_cb\n\n# Plot the correlation matrix\ncorr = pdf.corr()\nmask = np.zeros_like(corr, dtype=np.bool)\nmask[np.triu_indices_from(mask)] = True\nsns.set(style=\"white\")\nsns.heatmap(corr, mask=mask, annot=True,\n            square=True, linewidths=.5,\n            cbar_kws={'label': 'Correlation coefficient'})","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"On one hand, as we saw in the previous section the performance of the tree-based models is far better than the ridge regression.  On the other hand, the ridge regression's predictions aren't very highly correlated with the tree models', and including models whose predictions are not highly correlated in an ensemble usually leads to better performance.\n\nHow well can we predict customer loyalty if we average the predictions from all 4 models?"},{"metadata":{"trusted":true,"_uuid":"9f4761af0a4772a3bbe5e050f485abe70fbb8821"},"cell_type":"code","source":"# Compute the averaged predictions\nmean_preds = pdf.mean(axis=1)\n\n# Compute RMSE of averaged predictions\nroot_mean_squared_error(y_train, mean_preds)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And if we only use the tree-based models?"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Compute the averaged predictions\nmean_preds = pdf[['XGB', 'LGBM', 'CatBoost']].mean(axis=1)\n\n# Compute RMSE of averaged predictions\nroot_mean_squared_error(y_train, mean_preds)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Looks like adding the ridge regression doesn't really help the performance of the model when ensembling."},{"metadata":{"_uuid":"02c9fbaabd30fc3832fdc56bd08983c39468f115"},"cell_type":"markdown","source":"<a id='stacking-model'></a>\n## Stacking Model\n\nAnother ensembling method is called [stacking](https://en.wikipedia.org/wiki/Ensemble_learning#Stacking), which often performs even better than just averaging the predictions of the base learners.  With stacking, we first make predictions with the base learners.  Then, we use a separate \"meta\" model to make predictions for the target based of the predictions of the base learners.  That is, we train the meta learner with the features being the predictions of the base learners (and the target still being the target).\n\nWhen stacking it's important to ensure you train the meta-learner on out-of-fold predictions of the base learners (i.e. make cross-fold predictions with the base learners and train the meta learner on those).  The `StackedRegressor` model we'll use below handles this for us (see the code for it [on my github](http://github.com/brendanhasz/dsutils/blob/master/src/dsutils/ensembling.py#L144)).\n\nFirst let's try using all four models as base learners, and then a Bayesian ridge regression as the meta-learner:"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\n# Create the ensemble regressor\nmodel = StackedRegressor([ridge, xgboost, catboost, lgbm],\n                         meta_learner=BayesianRidge())\n\n# Performance of ensemble\ncross_val_metric(model, X_train, y_train,\n                 metric=root_mean_squared_error,\n                 cv=3, display='RMSE')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And how well do we do if we just use the tree models?"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\n# Create the ensemble regressor\nmodel = StackedRegressor([xgboost, catboost, lgbm],\n                         meta_learner=BayesianRidge())\n\n# Performance of ensemble\ncross_val_metric(model, X_train, y_train,\n                 metric=root_mean_squared_error,\n                 cv=3, display='RMSE')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Again, adding the ridge regression to the stacking ensemble didn't really seem to help either.  Since it's a bit faster, we'll use the stacking model which uses only the tree models.  Note that the ensemble model in this case is just barely better than the tree-based models by themselves!  \n\nFinally, we just need to save our predictions to file:"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Fit model on training data\nfit_model = model.fit(X_train, y_train)\n\n# Predict on test data\npredictions = fit_model.predict(X_test)\n\n# Write predictions to file\ndf_out = pd.DataFrame()\ndf_out['card_id'] = X_test.index\ndf_out['target'] = predictions\ndf_out.to_csv('predictions.csv', index=False)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}