{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"markdown","source":"# Results analysis of simple Lgbm model: visualization and thoughts\n\n# [Update] What is the impact of negative downsampling ? \n\nThis notebook investigates the discrepancies between actual and predicted revenues based on a simple Light GBM model. \n\n- The first section considers a simple LGBM model trained on the entire data set at the session level. \n- The second section investigates the effect of negative downsampling on the predictions and distribution of results\n\n Of interest to some users may also be the helper functions used to visualize the model's outputs.  \n \nI have hidden most of the code for the model as the primary interest is the discussion of the errors and visualizations. Feel free to post your comments / thoughts !\n\n--------------------------------\n\n##  Light GBM model\nThe model is a pretty standard one, largely based on popular public kernels. It does not perform great on LB but has reasonable user level CV = 1.5821\n\nI have added / removed the following features:\n\n**Time features**\n- Local hour of the day (time zone aware) \n- Time since last visit\n- Day of week, month of year, day of year\n*Removed*: visitStartTime (as this decreases CV but increases LB)\n\n**Objective and CV**\n- Predicting at session level on log1p of revenues\n- 5 fold cross validation, with model prediction on test and averaging over each fold \n- RMSE and gbdt with standard hyperparameters\n\n**Encoding**: label encoding of categorical features (note, this is done elsewhere)\n\n## Negative downsampling\nA good recamp is given [here][1]. Quoting the author:\n\n    **Negative downsampling (imbalanced)**: different samples sizes are used in this procedure. In all these samples all the observations from the minority class are kept and we take different number of observations from the majority class by performing sampling without replacement.\n\nHere the difficulty is deermining the ratio of positive to negative downsamples as well as the cross validation scheme to use. \n\n# Preparation code\n\n[1]: https://medium.com/bluekiri/dealing-with-highly-imbalanced-classes-7e36330250bc"},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true,"_uuid":"14321ca71665c2dc685760a0277fc74c3056ec95"},"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nimport time\nimport warnings\nimport datetime\n\n#Plotting\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nplt.style.use('ggplot')\n\n#Sklearn\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.model_selection import KFold\n\n#lgm and graph viz\nimport graphviz \nimport lightgbm as lgb\n\nwarnings.filterwarnings('ignore')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"b46292b12693eb224ceb7261a73bc1c17159f5e6","_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"os.listdir('../input/kernel-for-file-processing-2')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5caf2bc414a15cfa635b035f52cc9ba1313b8d0d","_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# Extract target values and Ids\ncat_cols = ['channelGrouping','device.browser',\n       'device.deviceCategory', 'device.isMobile', 'device.operatingSystem',\n       'geoNetwork.city', 'geoNetwork.continent', 'geoNetwork.country',\n       'geoNetwork.metro', 'geoNetwork.networkDomain', 'geoNetwork.region',\n       'geoNetwork.subContinent','trafficSource.adContent',\n       #'trafficSource.adwordsClickInfo.adNetworkType',\n       'trafficSource.adwordsClickInfo.gclId',\n       #'trafficSource.adwordsClickInfo.isVideoAd',\n       #'trafficSource.adwordsClickInfo.page',\n       #'trafficSource.adwordsClickInfo.slot', #Drop as only 3 values and always poor\n        'trafficSource.campaign',\n       'trafficSource.isTrueDirect', 'trafficSource.keyword',\n       'trafficSource.medium', 'trafficSource.referralPath',\n       'trafficSource.source'  ]\n\nto_drop = ['trafficSource.adwordsClickInfo.adNetworkType','trafficSource.adwordsClickInfo.isVideoAd',\n          'trafficSource.adwordsClickInfo.page','trafficSource.adwordsClickInfo.slot']\n\nnum_cols = ['visitNumber', 'totals.bounces', 'totals.hits',\n            'totals.newVisits', 'totals.pageviews', ]\n\n#interaction_cols = ['totals.hits / totals.pageviews']\n\nvisitStartTime = ['visitStartTime']\n\ntime_cols = ['_dayofweek', '_monthofyear', '_dayofyear', '_local_hourofday', '_time_since_last_visit']\n\nID_cols = ['date', 'fullVisitorId', 'sessionId', 'visitId']\n\ntarget_col = ['totals.transactionRevenue']\n\n","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true,"_uuid":"25537f11c7332e6e6cc55d22036d74862fcc1eb6"},"cell_type":"code","source":"%%time\n\ntrain_df = pd.read_pickle('../input/kernel-for-file-processing-2/train_flat_FE.pkl')\ntest_df = pd.read_pickle('../input/kernel-for-file-processing-2/test_flat_FE.pkl')\n\ntrain_df.drop(to_drop, axis = 1, inplace = True)\ntest_df.drop(to_drop, axis = 1, inplace = True)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6661664721784e12a4d62a3f6680b20a88f34055"},"cell_type":"markdown","source":"## Features "},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true,"_uuid":"dc906c8b468438f9457612410750b858b2d2faa7"},"cell_type":"code","source":"#Time features\ntrain_df['_dayofweek'] = train_df['visitStartTime'].dt.dayofweek\ntrain_df['_monthofyear'] = train_df['visitStartTime'].dt.month\ntrain_df['_dayofyear'] = train_df['visitStartTime'].dt.dayofyear\n#train_df['_dayofmonth'] = train_df['visitStartTime'].dt.day\n\ntest_df['_dayofweek'] = test_df['visitStartTime'].dt.dayofweek\ntest_df['_monthofyear'] = test_df['visitStartTime'].dt.month\ntest_df['_dayofyear'] = test_df['visitStartTime'].dt.dayofyear\n#test_df['_dayofmonth'] = test_df['visitStartTime'].dt.day","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true,"_uuid":"6a2c522f41f290251656f0623b99224f3cdfc5b6"},"cell_type":"code","source":"%%time\n\n#Numeric as float ##\nfor n in [num_cols + time_cols]:\n    train_df[n] = train_df[n].fillna(0).astype('int')\n    test_df[n] = test_df[n].fillna(0).astype('int')\n    \n#Time as float\ntrain_df['_time_since_last_visit'] = pd.to_numeric(train_df['_time_since_last_visit'])\ntest_df['_time_since_last_visit'] = pd.to_numeric(test_df['_time_since_last_visit'])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"74f9083151de8cde26b8dfa62477028d54fddcfc"},"cell_type":"markdown","source":"## Label encoding "},{"metadata":{"trusted":true,"_uuid":"f52a16bdb952688ab471f2cc42a95ab7301fa042","_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"%%time\n\ntrain_df[cat_cols] = train_df[cat_cols].fillna('unknown')\ntest_df[cat_cols] = test_df[cat_cols].fillna('unknown')\ntest_df['device.isMobile'] = test_df['device.isMobile'].astype('int')\n\n#Factorize cats\nfor f in cat_cols:\n    train_df[f], indexer = pd.factorize(train_df[f])\n    test_df[f] = indexer.get_indexer(test_df[f])","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"train_df['totals.transactionRevenue'] = train_df['totals.transactionRevenue'].fillna(0).astype('float')\n\n#Index\ntrain_idx = train_df['fullVisitorId']\ntest_idx = test_df['fullVisitorId']\n\n#Targets\ntrain_target = np.log1p(train_df.groupby(\"fullVisitorId\")[\"totals.transactionRevenue\"].sum())\ntrain_y = np.log1p(train_df[\"totals.transactionRevenue\"])\n\n#Datasets\ntrain_X = train_df[cat_cols + num_cols + time_cols].copy()\ntest_X = test_df[cat_cols + num_cols + time_cols ].copy()\n\nprint(train_X.shape)\nprint(test_X.shape)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e833490e9f8d6d3b8b763bc355dcf080c2252dda"},"cell_type":"markdown","source":"# Light GBM\n## Parameters"},{"metadata":{"trusted":true,"_uuid":"b82effd29b1ef85335e9f2b9a58eb37989bf84f4"},"cell_type":"code","source":"from lightgbm import LGBMRegressor\n\n#Initialize LGBM\ngbm = LGBMRegressor(objective = 'regression', \n                     boosting_type = 'gbdt', \n                     metric = 'rmse',\n                     n_estimators = 10000, #10000\n                     num_leaves = 10, #10\n                     learning_rate = 0.08, #0.01\n                     bagging_fraction = 0.9,\n                     feature_fraction = 0.3,\n                     bagging_seed = 0,\n                     max_depth = 10,\n                                         )\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fd12cc707cb2bccae633813ee5eb0d1357c8c938"},"cell_type":"markdown","source":"## Fit the model\n**In a nutshell**: K-fold training where each fold is used once for early stopping validation. At each fold, a test prediction is made using the trained model. Final prediction is an average of the K predictions\n\n**Steps:**\n- Fit the LGBM model K times on the dataset - the Kth fold\n - For each fitted model, predict on validation set (oof_pred) and on test set (sub_preds)\n - Average presub_preds for final predictions\n\nIdea for averaged models comes from Olivier and: https://www.kaggle.com/sz8416/lb-1-4439-gacr-prediction-eda-lgb-baseline"},{"metadata":{"trusted":true,"_uuid":"df4147ca5ac8b84997abb78631d16667977fa88f","_kg_hide-output":true,"_kg_hide-input":true},"cell_type":"code","source":"%%time\n#Initilization\nall_K_fold_results = []\nkf = KFold(n_splits=5, shuffle = True)\noof_preds = np.zeros(train_X.shape[0])\nsub_preds = np.zeros(test_X.shape[0])\n\n\nfor dev_index, val_index in kf.split(train_X):\n    X_dev, X_val = train_X.iloc[dev_index], train_X.iloc[val_index]\n    y_dev, y_val = train_y[dev_index], train_y[val_index]\n\n    #Fit the model\n    model = gbm.fit(X_dev,y_dev, eval_set=[(X_val, y_val)],verbose = 100, \n                    eval_metric = 'rmse', early_stopping_rounds = 100) ##100\n    \n    #Predict out of fold \n    oof_preds[val_index] = gbm.predict(X_val, num_iteration= model.best_iteration_)\n    \n    oof_preds[oof_preds < 0] = 0\n    \n    #Predict on test set based on current fold model. Average results\n    sub_prediction = gbm.predict(test_X, num_iteration= model.best_iteration_) / kf.n_splits\n    sub_prediction[sub_prediction<0] = 0\n    sub_preds = sub_preds + sub_prediction\n    \n    #Save current fold values\n    fold_results = {'best_iteration_' : model.best_iteration_, \n                   'best_score_' : model.best_score_['valid_0']['rmse'], \n                   'evals_result_': model.evals_result_['valid_0']['rmse'],\n                   'feature_importances_' : model.feature_importances_}\n\n    all_K_fold_results.append(fold_results.copy())\n    \n\nresults = pd.DataFrame(all_K_fold_results)\n\n\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f193c2e3f0c793609059e8188b25f89556e45684"},"cell_type":"markdown","source":"# Utility functions\n## Visualization, RMSE and saving"},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true,"_uuid":"bad61b5037669b13b4b06197f95339341ccce86f"},"cell_type":"code","source":"def RMSE_log_sum(pred_val, val_df):\n    #set negative values to zero\n    pred_val[pred_val < 0] = 0\n    \n    #Build new dataframe\n    val_pred_df = pd.DataFrame(data = {'fullVisitorId': val_df['fullVisitorId'].values, \n                                       'transactionRevenue': val_df['totals.transactionRevenue'].values,\n                                      'predictedRevenue':np.expm1(pred_val) })\n    #Compute sum\n    val_pred_df = val_pred_df.groupby('fullVisitorId').sum().reset_index()\n\n    mse_log_sum = mean_squared_error( np.log1p(val_pred_df['transactionRevenue'].values), \n                             np.log1p(val_pred_df['predictedRevenue'].values)  )\n\n    #print('log (sum + 1): ',np.sqrt(mse_log_sum))\n    return np.sqrt(mse_log_sum)\n\n\ndef save_submission(pred_test, test_df, file_name):\n    #Zero negative predictions\n    pred_test[pred_test < 0] = 0\n    \n    #Create temporary dataframe\n    sub_df = pd.DataFrame(data = {'fullVisitorId':test_df['fullVisitorId'], \n                             'predictedRevenue':np.expm1(pred_test)})\n    sub_df = sub_df.groupby('fullVisitorId').sum().reset_index()\n    sub_df.columns = ['fullVisitorId', 'predictedLogRevenue']\n    sub_df['predictedLogRevenue'] = np.log1p(sub_df['predictedLogRevenue'])\n    sub_df.to_csv(file_name, index = False)\n\n    \ndef visualize_results(results):\n#Utility function to plot fold loss and best model feature importance\n    plt.figure(figsize=(16, 12))\n\n    #----------------------------------------\n    # Plot validation loss\n    plt.subplot(2,2,1)\n\n    for K in range(results.shape[0]):\n        plt.plot(np.arange(len(results.evals_result_[K])), results.evals_result_[K], label = 'fold {}'.format(K))\n\n    plt.xlabel('Boosting iterations')\n    plt.ylabel('RMSE')\n    plt.title('Validation loss vs boosting iterations')\n    plt.legend()\n\n    #----------------------------------------\n    # Plot box plot of RMSE\n    plt.subplot(2, 2, 2)    \n    scores = results.best_score_\n    plt.boxplot(scores)\n    rmse_mean = np.mean(scores)\n    rmse_std = np.std(scores)\n    plt.title('RMSE Mean:{:.3f} Std: {:.4f}'.format(rmse_mean,rmse_std ))\n    \n    #----------------------------------------\n    # Plot feature importance\n    #feature_importance = results.sort_values('best_score_').feature_importances_[0]\n    df_feature_importance = pd.DataFrame.from_records(results.feature_importances_)\n    feature_importance = df_feature_importance.mean()\n    std_feature_importance = df_feature_importance.std()\n    \n    # make importances relative to max importance\n    #feature_importance = 100.0 * (mean_feature_importance / mean_feature_importance.sum())\n    sorted_idx = np.argsort(feature_importance)\n    pos = np.arange(sorted_idx.shape[0]) + .5\n    plt.subplot(2, 1, 2)\n    plt.bar(pos, feature_importance[sorted_idx], align='center', yerr = std_feature_importance)\n    xlabels = [ train_X.columns.values[i] for i in sorted_idx]\n    plt.xticks(pos, xlabels, rotation = 90)\n    plt.xlabel('Feature')\n    plt.ylabel('Avg Importance score')\n    plt.title('Mean Feature Importance over K folds') \n    \n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"bd4106c976bed4d1f6bfb1b4b528f2408ac3c85b"},"cell_type":"code","source":"print('Session level CV score: ', np.mean(results.best_score_))\nprint('User level CV score: ', RMSE_log_sum(oof_preds, train_df))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0699f9308757485537d05fc042f8aaec52d5bded"},"cell_type":"markdown","source":"## Model visualization \n- Feature importance, validation RMSE and box plot\n- View of first decision tree"},{"metadata":{"trusted":true,"_uuid":"f6391e0169800dbe2df8318892b8ed7b8147aec9","_kg_hide-output":true,"_kg_hide-input":true},"cell_type":"code","source":"import graphviz \ndot_data = lgb.create_tree_digraph(model, tree_index = 1,show_info=['split_gain'])\n\ngraph = graphviz.Source(dot_data)  \ngraph ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"188d9553f6a1e937a0d6507830554e6352f40162"},"cell_type":"markdown","source":"# Error analysis\n\nThis section is the main purpose of this kernel, your thoughts and comments are welcome. \n\nNote that unless stated otherwise, plots are of **log of revenues, at the session level**. \n\n### Utility dataframe\nIt appears that the log revenue of the first 10 sessions are well predicted. Zero true revenues correspond to zero predicted revenues. Lets look at the overall distributions"},{"metadata":{"trusted":true,"_uuid":"2a158042c7f4baaccd0d83d5537124ce99ba9ca8","_kg_hide-input":true},"cell_type":"code","source":"error_df = pd.DataFrame(data = {'visitStartTime':train_df['visitStartTime'],'fullVisitorId':train_df['sessionId'], \n                                'True_log_revenue' : np.log1p(train_df['totals.transactionRevenue']), \n                                'Predicted_log_revenue':oof_preds  })\n\nerror_df['Difference'] = error_df['True_log_revenue'] - error_df['Predicted_log_revenue']\nerror_df['True_is_non_zero'] = error_df['True_log_revenue'] > 0\n#temp_df.columns = ['fullVisitorId', 'predictedLogRevenue']\n#sub_df['predictedLogRevenue'] = np.log1p(sub_df['predictedLogRevenue'])\n#sub_df.to_csv(file_name, index = False)\nerror_df.sort_values('visitStartTime').head(10)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"dbd24c9a2ba4cd08bca6deafc39cdf676c3c9746"},"cell_type":"markdown","source":"## Distributions of true and predicted log revenues\n\n- **Left hand side**: distribution of predicted revenues for sessions which have *zero true revenues*\n- **Right hand side**: distribution of true and predicted revenues for sessions which have *non zero true revenues* i.e. totals.transactionRevenue > 0\n\n### Comments\n\nThe distribution of predictions of log revenues for true zero sessions is strongly centered around the zero value, with a long right tail. I would argue that this isn't too bad for a start.\n\nFor true non zero sessions, however, the picture is very different:\n- The mean of the predicted values is clearly smaller than true values (bias) \n- The variance of the predicted value is larger \n- The shape is somewhat bell shaped (manually truncated at zero) \n- The overlap between the predicted and true values is very small\n"},{"metadata":{"trusted":true,"_uuid":"29f2d0e6d4152e99aa8c6e4dbc42cc4f67e83720","_kg_hide-input":true},"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize = (20,7))\n\n\nsns.distplot(error_df[error_df['True_is_non_zero'] == False]['True_log_revenue'], ax = ax1, label = 'true')\nsns.distplot(error_df[error_df['True_is_non_zero'] == False ]['Predicted_log_revenue'], ax = ax1, label = 'pred')\nax1.legend()\nax1.set_ylim(0,.1)\nax1.set_xlabel('Log revenue (session)')\nax1.set_title('Distribution of log revenues for sessions with zero true revenue ')\n\nsns.distplot(error_df[error_df['True_is_non_zero'] == True]['True_log_revenue'], ax = ax2, label = 'true')\nsns.distplot(error_df[error_df['True_is_non_zero'] == True ]['Predicted_log_revenue'], ax = ax2, label = 'pred')\nax2.legend()\nax2.set_ylim(0,.5)\nax2.set_xlabel('Log revenue (session)')\nax2.set_title('Distribution of log revenues for sessions with non zero true revenue ')\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9eefcc214289b419d5d351606518e306606251e2"},"cell_type":"markdown","source":"### Predicted distribution for all data points\n\nIn effect, what we have here is a strongly *bi-modal* true distribution which the lgbm is trying to predict. It does faily well on the zero sessions, but poorly on the non zero sessions. Note that these are log values at session level, and thus need to be exponentiated, summed and then log again...\n"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"842adba8750054b2d2a304349b88ff87779bd87e"},"cell_type":"code","source":"fig, ax = plt.subplots(nrows=1, ncols=1, figsize = (15,7))\n\nsns.distplot(error_df['True_log_revenue'], ax = ax, label = 'true')\nsns.distplot(error_df['Predicted_log_revenue'], ax = ax, label = 'pred')\nax.legend()\nax.set_ylim(0,.04)\nax.set_xlabel('Log revenue (session)')\nax.set_title('Distribution of log revenues for all sessions')\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7c7bf37ef58dfdbbae2e33bfbdb6d767955df864"},"cell_type":"markdown","source":"## Predictions and errors over time  \n\nLooking at the plot of true and predicted log revenues for true non zero sessions only (top plots) we can clearly see the high bias and variance in our predictions. Time doesn't seem to play a particular role and there are no clear trends or seasonality. Checks on auto-correlation and stationarity could be performed here (to do) . Clearly, the model is doing poorly for these data points\n\nWhen looking at true and predicted values for *all sessions* (i.e. both zero and non zero true values), the model performs much better. Monthly and daily average patterns are well predicted, the main issues is that extreme values for daily revenues are poorly predicted. See plot on bottom right corner"},{"metadata":{"trusted":true,"_uuid":"0e4619247dea80a37901a06389e6441aa309ab68","_kg_hide-input":true},"cell_type":"code","source":"sorted_non_zero = error_df[error_df['True_is_non_zero'] == True].sort_values('visitStartTime')\n\nplt.figure(figsize = (20,15))\nplt.subplot(2,2,1)\nplt.plot(sorted_non_zero.visitStartTime, sorted_non_zero.True_log_revenue , label = 'True')\nplt.plot(sorted_non_zero.visitStartTime, sorted_non_zero.Predicted_log_revenue , alpha = .5, label = 'Pred')\nplt.title('Log revenue over time (non zero true sessions only)')\nplt.legend()\nplt.xlabel('Time: sessions')\n\nplt.subplot(2,2,2)\ndaily_error_non_zero_df = sorted_non_zero.set_index('visitStartTime', drop = True).resample('D').mean()\nplt.plot(daily_error_non_zero_df.index, daily_error_non_zero_df.True_log_revenue , label = 'True')\nplt.plot(daily_error_non_zero_df.index, daily_error_non_zero_df.Predicted_log_revenue , label = 'Pred')\nplt.title('Daily average log revenue (non zero true sessions only)')\n\nplt.subplot(2,2,3)\nweekly_error_df = error_df.set_index('visitStartTime', drop = True).resample('W').mean()\nplt.plot(weekly_error_df.index, weekly_error_df.True_log_revenue , label = 'True')\nplt.plot(weekly_error_df.index, weekly_error_df.Predicted_log_revenue , label = 'Pred')\nplt.title('Weekly average log revenue (all session)')\n\n\nplt.subplot(2,2,4)\ndaily_error_df = error_df.set_index('visitStartTime', drop = True).resample('D').mean()\nplt.plot(daily_error_df.index, daily_error_df.True_log_revenue , label = 'True')\nplt.plot(daily_error_df.index, daily_error_df.Predicted_log_revenue , label = 'Pred')\nplt.title('Daily average log revenue (all session)')\n\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c7a96fff1b64a8c55e77d770ec8620fb89e36900"},"cell_type":"markdown","source":"## Residuals over time\n\n- Residuals = True - Predicted\n- No obvious patterns over time\n- Clear bias for true non zero sessions\n- Manual truncation of negative values helps the true zero session predictions\n"},{"metadata":{"trusted":true,"_uuid":"ada25f70841a12b497c964fa70c0475bdf59be25","_kg_hide-input":true},"cell_type":"code","source":"sorted_non_zero = error_df[error_df['True_is_non_zero'] == True].sort_values('visitStartTime')\nsorted_zero = error_df[error_df['True_is_non_zero'] == False].sort_values('visitStartTime')\n\n\nplt.figure(figsize = (20,5))\nplt.subplot(1,3,1)\nts_error_df = error_df.set_index('visitStartTime', drop = True)\ndifference_rev_df = error_df.sort_values('visitStartTime')\nplt.plot(error_df.visitStartTime, error_df.Difference , label = 'True - predicted', color = 'grey')\nplt.title('Train - Pred (log rev) for all sessions')\n\nplt.subplot(1,3,2)\nplt.plot(sorted_non_zero.visitStartTime, sorted_non_zero.Difference , label = 'True - predicted',\n         color = 'grey')\nplt.title('Train - Pred for non zero sessions only')\n\nplt.subplot(1,3,3)\nplt.plot(sorted_zero.visitStartTime, sorted_zero.Difference,\n         color = 'grey')\nplt.title('Train - Pred for zero sessions only')\n\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e714667d931c99fb4d4315a49ee4c7c44a10d1af"},"cell_type":"markdown","source":"## Joint distribution of true and predicted revenues\nFor true non zero sessions only"},{"metadata":{"trusted":true,"_uuid":"a48de8af0a395f8c24584d1edd8216b7ec4c741f","_kg_hide-input":true},"cell_type":"code","source":"\nsns.jointplot(x=\"True_log_revenue\", y=\"Predicted_log_revenue\", data=sorted_non_zero)\ndisplay('Joint distribution of log rev for non zero sessions only')\n\nplt.show()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"00979940c00037743cdcf6c64a57c3ebd33df356"},"cell_type":"markdown","source":"# Part 2) Negative downsampling\n\nLets investivate an alternative approach using negative downsampling. In other words, artificially adjusting the proportion of positive and negative samples. \n\nTo start with, we will set this proportion at 1/2, i.e. 50% positive (with true revenue >0) and %)% negative (with true revenue = 0) \n\n### Open question:\n\nWhat cross validation approach would work in this case ? Current approach is to use a single dev/val split, where the validation set is used to perform early stopping and the RMSE and error analysis is also performed on the validation set.\n\nA cross validated approach would be useful here but I still need to think about how to make this work with the ramdom sampling of the negative samples... *help wanted*"},{"metadata":{"_uuid":"a33641019255646ed637ae41754dc1d4ba37d9be"},"cell_type":"markdown","source":"## Lgbm with negative downsampling\n\n"},{"metadata":{"trusted":true,"_uuid":"9621c938ef2b22ef06d990532cb2bcb1bc240135","_kg_hide-output":true},"cell_type":"code","source":"%%time\nall_K_fold_results = []\n\n#Downsampling preparation\nidx_pos = train_df[train_df['totals.transactionRevenue'] > 0].index\nidx_neg = train_df[train_df['totals.transactionRevenue'] == 0].index\ntrain_pos_X, train_pos_y = train_X.loc[idx_pos], train_y[idx_pos]\ntrain_neg_X, train_neg_y = train_X.loc[idx_neg], train_y[idx_neg]\n\noof_preds = [] \nsub_preds = np.zeros(test_X.shape[0])\nsize_dev = round( (4/5) * len(idx_pos) )  #4/5 for dev and 1/5 for val\n\n#Indices\nidx_dev_pos = np.random.choice(idx_pos, size = size_dev, replace = False)\nidx_val_pos = np.setdiff1d(idx_pos,idx_dev_pos)\n\nidx_dev_neg = np.random.choice(idx_neg, size = size_dev, replace = False)\nidx_val_neg = np.random.choice(idx_neg, size = len(idx_val_pos), replace = False)\n\n#Datasets\nX_dev_pos, X_val_pos = train_pos_X.loc[idx_dev_pos], train_pos_X.loc[idx_val_pos]\ny_dev_pos, y_val_pos = train_pos_y[idx_dev_pos], train_pos_y[idx_val_pos]\n\nX_dev_neg, X_val_neg = train_neg_X.loc[idx_dev_neg], train_neg_X.loc[idx_val_neg]\ny_dev_neg, y_val_neg = train_neg_y[idx_dev_neg], train_neg_y[idx_val_neg]\n\n#Concatenate\nX_dev, X_val = pd.concat([X_dev_pos, X_dev_neg], axis = 0), pd.concat([X_val_pos, X_val_neg], axis = 0)\ny_dev, y_val = pd.concat([y_dev_pos, y_dev_neg]), pd.concat([y_val_pos, y_val_neg])\n\n#Fit the model\nmodel = gbm.fit(X_dev,y_dev, eval_set=[(X_val, y_val)],verbose = 100, \n                eval_metric = 'rmse', early_stopping_rounds = 100) ##100\n\n#Predict on val\noof_preds.append( gbm.predict(X_val, num_iteration= model.best_iteration_).copy() )\n\n#Predict on test set based on current fold model. \nsub_prediction = gbm.predict(test_X, num_iteration= model.best_iteration_)\nsub_prediction[sub_prediction<0] = 0\nsub_preds = sub_preds + sub_prediction #not needed\n\n#Save current fold values\nfold_results = {'best_iteration_' : model.best_iteration_, \n               'best_score_' : model.best_score_['valid_0']['rmse'], \n               'evals_result_': model.evals_result_['valid_0']['rmse'],\n               'feature_importances_' : model.feature_importances_}\n\nall_K_fold_results.append(fold_results.copy())\n\n\nresults2 = pd.DataFrame(all_K_fold_results)\n\n#Save as array and flatten\noof_preds = np.asarray(oof_preds).flatten()\noof_preds[oof_preds < 0] = 0\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"706ba8b746d81c1da6f44825d11284340e494dfe"},"cell_type":"code","source":"print('Session level validation score: ', np.mean(results2.best_score_))\nprint('User level validation score: ', RMSE_log_sum(oof_preds, train_df.iloc[X_val.index]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8f02aedd862edead19ab2dfb50dc9b9e00f559df"},"cell_type":"code","source":"X_val.index","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e28a65efdf619e5d9226bd139877b516c0139b3e"},"cell_type":"code","source":"visualize_results(results2)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b4fef10d3ca95dfacdaa51f4871ad1215f478df0"},"cell_type":"markdown","source":"# Error analysis for downsampled model \n\n ### Inspecting the first 20 rows"},{"metadata":{"trusted":true,"_uuid":"0b5fdc802426fdf655d6665aad180b1c46843cb3","_kg_hide-input":true},"cell_type":"code","source":"idx_oof_pred= X_val.index\n\nerror_df = pd.DataFrame(data = {'visitStartTime':train_df.loc[idx_oof_pred]['visitStartTime'],\n                                'fullVisitorId':train_df.loc[idx_oof_pred]['sessionId'], \n                                'True_log_revenue' : np.log1p(train_df.loc[idx_oof_pred]['totals.transactionRevenue']), \n                                'Predicted_log_revenue':oof_preds  })\n\nerror_df['Difference'] = error_df['True_log_revenue'] - error_df['Predicted_log_revenue']\nerror_df['True_is_non_zero'] = error_df['True_log_revenue'] > 0\nerror_df.head(10)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"539250cc3f973f7344045265ba852d4aac4f0b3c"},"cell_type":"markdown","source":"### Breakdown of error: Sum of Squares of Error (SSE)\nPredictions on true zero revenue sessions account for the majority of the error. Lets confirm this visually\n"},{"metadata":{"trusted":true,"_uuid":"43033e7f57d8df4ceaf9309f00a9fa1ad3734971"},"cell_type":"code","source":"Total_SSE = np.square(error_df.Difference).sum()\nZero_SSE = np.square(error_df[error_df['True_is_non_zero'] == False].Difference).sum()\nNon_zero_SSE = np.square(error_df[error_df['True_is_non_zero'] == True].Difference).sum() \n\nprint('Total SSE: ',round(Total_SSE), ' 100%'  )\nprint('Zero true revenue SSE: ', round(Zero_SSE), round(Zero_SSE / Total_SSE , 3) * 100, '%'  )\nprint('Non zero true revenue SSE: ',round(Non_zero_SSE), round(Non_zero_SSE / Total_SSE, 3) * 100, '%' )","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9b52c206065cab652db69668f1437fab36321314"},"cell_type":"markdown","source":"## Distributions of true and predicted log revenues (downsampled model)\n\n**Left hand side**: distribution of predicted revenues for sessions which have *zero true revenues*\n - As previously, a long right hand side tail for predicted revenues\n\n**Right hand side**: distribution of true and predicted revenues for sessions which have *non zero true revenues* i.e. totals.transactionRevenue > 0\n- Much improved distribution ! The shape, size, mean and variance of the predicted distribution fits much closer to reality. Predictions are slightly skewed with a thick left tail. \n"},{"metadata":{"trusted":true,"_uuid":"12b6bc860c8f01ae9e99f9753f20cef7c1948df0"},"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize = (20,7))\n\n\nsns.distplot(error_df[error_df['True_is_non_zero'] == False]['True_log_revenue'], ax = ax1, label = 'true')\nsns.distplot(error_df[error_df['True_is_non_zero'] == False ]['Predicted_log_revenue'], ax = ax1, label = 'pred')\nax1.legend()\nax1.set_ylim(0,.1)\nax1.set_xlabel('Log revenue (session)')\nax1.set_title('Distribution of log revenues for sessions with zero true revenue ')\n\nsns.distplot(error_df[error_df['True_is_non_zero'] == True]['True_log_revenue'], ax = ax2, label = 'true')\nsns.distplot(error_df[error_df['True_is_non_zero'] == True ]['Predicted_log_revenue'], ax = ax2, label = 'pred')\nax2.legend()\nax2.set_ylim(0,.5)\nax2.set_xlabel('Log revenue (session)')\nax2.set_title('Distribution of log revenues for sessions with non zero true revenue ')\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"90940ae95d6358a58169d984f43315d251079c58"},"cell_type":"markdown","source":"## Distribution on all datapoints (downsampled set only) \n\nHere the *bi-modal* nature of the prediction is much more visible. The error however is larger, mostly due to incorrectly predeicted revenues for zero, or small true revenues. "},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"d5507231f2b7f464ee0ffaaa7d58096183f9d9de"},"cell_type":"code","source":"fig, ax = plt.subplots(nrows=1, ncols=1, figsize = (15,7))\n\nsns.distplot(error_df['True_log_revenue'], ax = ax, label = 'true', bins = 100)\nsns.distplot(error_df['Predicted_log_revenue'], ax = ax, label = 'pred', bins = 100)\nax.legend()\nax.set_ylim(0,.1)\nax.set_xlabel('Log revenue (session)')\nax.set_title('Distribution of log revenues for all sessions')\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8b30c1dcedb8189b4c2b0428cefb1f19a56a5c54"},"cell_type":"markdown","source":"## Over time "},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"649a69fecb6d7408bdd5ab0c3a43f2e6cc270a6e"},"cell_type":"code","source":"sorted_non_zero = error_df[error_df['True_is_non_zero'] == True].sort_values('visitStartTime')\n\nplt.figure(figsize = (20,15))\nplt.subplot(2,2,1)\nplt.plot(sorted_non_zero.visitStartTime, sorted_non_zero.True_log_revenue , label = 'True')\nplt.plot(sorted_non_zero.visitStartTime, sorted_non_zero.Predicted_log_revenue , alpha = .5, label = 'Pred')\nplt.title('Log revenue over time (non zero true sessions only)')\nplt.legend()\nplt.xlabel('Time: sessions')\n\nplt.subplot(2,2,2)\ndaily_error_non_zero_df = sorted_non_zero.set_index('visitStartTime', drop = True).resample('D').mean()\nplt.plot(daily_error_non_zero_df.index, daily_error_non_zero_df.True_log_revenue , label = 'True')\nplt.plot(daily_error_non_zero_df.index, daily_error_non_zero_df.Predicted_log_revenue , label = 'Pred')\nplt.title('Daily average log revenue (non zero true sessions only)')\n\nplt.subplot(2,2,3)\nweekly_error_df = error_df.set_index('visitStartTime', drop = True).resample('W').mean()\nplt.plot(weekly_error_df.index, weekly_error_df.True_log_revenue , label = 'True')\nplt.plot(weekly_error_df.index, weekly_error_df.Predicted_log_revenue , label = 'Pred')\nplt.title('Weekly average log revenue (all session)')\n\n\nplt.subplot(2,2,4)\ndaily_error_df = error_df.set_index('visitStartTime', drop = True).resample('D').mean()\nplt.plot(daily_error_df.index, daily_error_df.True_log_revenue , label = 'True')\nplt.plot(daily_error_df.index, daily_error_df.Predicted_log_revenue , label = 'Pred')\nplt.title('Daily average log revenue (all session)')\n\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e81ffe6a345e17e02768e05b47314c7c254b7513"},"cell_type":"markdown","source":"## Residuals over time\nWhereas previously the residuals over time did not seem to follow any particular pattern, we see here (by inspection only) a shape more representative of time series fluctuations and stochastic processes. \n"},{"metadata":{"trusted":true,"_uuid":"28901e49a77a909604596a2e9f1a3d06a5c87844"},"cell_type":"code","source":"sorted_non_zero = error_df[error_df['True_is_non_zero'] == True].sort_values('visitStartTime')\nsorted_zero = error_df[error_df['True_is_non_zero'] == False].sort_values('visitStartTime')\nsorted_all = error_df.sort_values('visitStartTime')\n\n\nplt.figure(figsize = (20,5))\nplt.subplot(1,3,1)\nplt.plot(sorted_all.visitStartTime, sorted_all.Difference , label = 'True - predicted', color = 'grey')\nplt.title('Train - Pred (log rev) for all sessions')\n\nplt.subplot(1,3,2)\nplt.plot(sorted_non_zero.visitStartTime, sorted_non_zero.Difference , label = 'True - predicted',\n         color = 'grey')\nplt.title('Train - Pred for non zero sessions only')\n\nplt.subplot(1,3,3)\nplt.plot(sorted_zero.visitStartTime, sorted_zero.Difference,\n         color = 'grey')\nplt.title('Train - Pred for zero sessions only')\n\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f75ac98370d07c53154da893e4b32eb97f14a9ef"},"cell_type":"markdown","source":"### Joint distribution "},{"metadata":{"trusted":true,"_uuid":"9b42caf65729a3b6b7179a5b21468d8da8458164"},"cell_type":"code","source":"\nsns.jointplot(x=\"True_log_revenue\", y=\"Predicted_log_revenue\", data=sorted_non_zero)\ndisplay('Joint distribution of log rev for non zero sessions only')\n\nplt.show()\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1848cf8acf54bb17b929ade442b59ba63bbe96a7"},"cell_type":"markdown","source":"# Next steps\n\n- Repeat the negative downsampled modelling, and average results over many repetitions ? \n- Hyper parameter tuning (as these are the parameters used for full data set) \n- Implement cross validation \n..."},{"metadata":{"trusted":true,"_uuid":"84defbf32040123848fdd96e99ac22c28775cd10"},"cell_type":"code","source":"#Save and submit\nsave_submission(sub_preds, test_df, 'submission.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"af9d2df9d73ef5ebd0d94da549d42ef9286a09e4"},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"750de54a4f49440bb5e14f12798557a6f2bed052"},"cell_type":"markdown","source":"## Appendix - code snips and pieces"},{"metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true,"_uuid":"4f3ac1ed41745d3257519fc7df3be75fe851a4f3"},"cell_type":"code","source":"'''%%time\n\n#Downsampling preparation\nidx_pos = train_df[train_df['totals.transactionRevenue'] > 0].index\nidx_neg = train_df[train_df['totals.transactionRevenue'] == 0].index\ntrain_pos_X, train_pos_y = train_X.loc[idx_pos], train_y[idx_pos]\ntrain_neg_X, train_neg_y = train_X.loc[idx_neg], train_y[idx_neg]\n\n#Initilization\nall_K_fold_results = []\nkf = KFold(n_splits=5, shuffle = False)\n\noof_preds = [] #np.zeros(train_pos_X.shape[0] * 2)\nsub_preds = np.zeros(test_X.shape[0])\n\n\nfor i, v in kf.split(train_pos_X):\n    #Positive samples\n    idx_dev_pos, idx_val_pos  = train_pos_X.index[i], train_pos_X.index[v]\n    X_dev_pos, X_val_pos = train_pos_X.loc[idx_dev_pos], train_pos_X.loc[idx_val_pos]\n    y_dev_pos, y_val_pos = train_pos_y[idx_dev_pos], train_pos_y[idx_val_pos]\n\n    #Negative samples\n    idx_dev_neg = np.random.choice(idx_neg, size = len(idx_dev_pos), replace = False)\n    idx_val_neg = np.random.choice(idx_neg, size = len(idx_val_pos), replace = False)\n    \n    X_dev_neg, X_val_neg = train_neg_X.loc[idx_dev_neg], train_neg_X.loc[idx_val_neg]\n    y_dev_neg, y_val_neg = train_neg_y[idx_dev_neg], train_neg_y[idx_val_neg]\n\n    #Concatenate\n    X_dev, X_val = pd.concat([X_dev_pos, X_dev_neg], axis = 0), pd.concat([X_val_pos, X_val_neg], axis = 0)\n    y_dev, y_val = pd.concat([y_dev_pos, y_dev_neg]), pd.concat([y_val_pos, y_val_neg])\n    \n     #Fit the model\n    model = gbm.fit(X_dev,y_dev, eval_set=[(X_val, y_val)],verbose = 100, \n                    eval_metric = 'rmse', early_stopping_rounds = 100) ##100\n    \n    #Predict out of fold \n    oof_preds.append( gbm.predict(X_val, num_iteration= model.best_iteration_).copy() )\n    \n    #Predict on test set based on current fold model. Average results\n    sub_prediction = gbm.predict(test_X, num_iteration= model.best_iteration_) / kf.n_splits\n    sub_prediction[sub_prediction<0] = 0\n    sub_preds = sub_preds + sub_prediction\n\n    #Save current fold values\n    fold_results = {'best_iteration_' : model.best_iteration_, \n                   'best_score_' : model.best_score_['valid_0']['rmse'], \n                   'evals_result_': model.evals_result_['valid_0']['rmse'],\n                   'feature_importances_' : model.feature_importances_}\n\n    all_K_fold_results.append(fold_results.copy())\n\n\nresults = pd.DataFrame(all_K_fold_results)\n\n#Save as array and flatten\noof_preds = np.asarray(oof_preds).flatten()\noof_preds[oof_preds < 0] = 0 '''","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d866a606efae1564c7ce3dba8fd1e720d6b70eb0"},"cell_type":"code","source":"","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}