{"cells":[{"metadata":{"_cell_guid":"768f94f1-8d9c-4e7b-82c8-7f51f98db891","_uuid":"10eef87aa2f378f87177da0d71e6765990c93389"},"cell_type":"markdown","source":"A lot of people have looked at blending kernels and asked where the weights come from.  And the answer is usually something like, \"Intuition, after taking into account their public leaderboard scores and the correlations among their forecasts.\"  Understandably, people don't find that answer very satisfying.  So I made this notebook to show the elegant way of choosing blending weights.  I call it \"stacking\" instead of \"blending\" because it fits an explicit metamodel on top of the base models, but the distinction in meaning isn't clear cut.  Some will say it's still blending becuase the \"stacking model\" is linear and the validation data used to fit the weights are separate from the original training data.  (With k-fold validation, you could use out-of-fold predictions from the original training data, but k-fold validation may be the wrong approach given the time component in the data.)"},{"metadata":{"collapsed":true,"_cell_guid":"f571fc89-4de5-4f43-b9bf-c7c7390292a4","_uuid":"844d76db2685e46d9872df1269a873a2919bd8a1","trusted":true},"cell_type":"code","source":"# File containing validation data\n# (These are selected from the last day of the original training set\n#  to correspond to the times of day used in the test set.)\nVAL_FILE = '../input/training-and-validation-data-pickle/validation.pkl.gz'","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os\nfrom sklearn.linear_model import LogisticRegression\nfrom scipy.special import expit, logit\nfrom sklearn.metrics import roc_auc_score\n\nprint(os.listdir(\"../input\"))","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"18eb2734-ebcc-49ec-9194-4ddaad317c5f","_uuid":"cbc85f52d166fd3e4cb07f75f1d7c24e0fd1f526","trusted":true},"cell_type":"code","source":"almost_zero = 1e-10\nalmost_one = 1 - almost_zero","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"# Just names to identify the models\nbase_models = {\n    'lgb1 ': \"Python LGBM based on Pranav Pandya's R version\",\n    'wbftl': \"anttip's Wordbatch FM-FTRL\",\n    'nngpu': \"Downampled Neural Network run on GPU\"\n    }","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"2d708642-8f1b-4fd3-b33e-500def92647a","_uuid":"87839ab18b15bfc16f71016685a5540ca56bd6dc","trusted":true},"cell_type":"code","source":"# Files with validation set predictions from each of the base models\n# (These were fit on a subset of the training data that ends a day before\n#  the end of the full training set.)\ncvfiles = {\n    'lgb1 ': '../input/validate-pranav-lgb-model/pranav_lgb_val_nostop.csv',\n    'wbftl': '../input/validate-anttip-s-wordbatch-fm-ftrl-9711-version/wordbatch_fm_ftrl_val.csv',\n    'nngpu': '../input/gpu-validation/gpu_val1.csv'\n    }","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"5b5459d8-0ac9-4c8e-9c22-3e19517d5eb6","_uuid":"8e81ec28924b1440ab3604d71a45504f7998e3a4","trusted":true},"cell_type":"code","source":"# Files with test set predictions\n# (These were fit on the full training set\n#  or on a subset at the end, to accommodate memory limits.)\nsubfiles = {\n    'lgb1 ': '../input/try-pranav-s-r-lgbm-in-python/sub_lgbm_r_to_python_nocv.csv',\n    'wbftl': '../input/anttip-s-wordbatch-fm-ftrl-9711-version/wordbatch_fm_ftrl.csv',\n    'nngpu': '../input/talkingdata-gpu-example-with-multiple-runs/gpu_test2.csv'\n    }","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"f89413ac-a2ea-4c4d-8b66-72491bbb0d6c","_uuid":"59de9568978d88cd8823c14a6287a3a256a6a646","trusted":true},"cell_type":"code","source":"# Public leaderbaord scores, for comparison\nlbscores = {\n    'lgb1 ': .9694,\n    'wbftl': .9711,\n    'nngpu': .9678\n}","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"914b72c1-5cfe-4175-9de4-06af6aa38a45","_uuid":"1f7b5e8d323c5582b8a1d87f3da8152dd1fa83e1"},"cell_type":"markdown","source":"You can click on the \"Data\" tab and follow the links to the kernels that generated each of the outputs above. Usually there are \"forked from\" links that you can follow to see where the models originated. (In the case of Pranav's LGBM, though, the link is misleading, because the original was [a separate kernel in R](https://www.kaggle.com/pranav84/single-lightgbm-in-r-with-75-mln-rows-lb-0-9690), and the fork was from a different Python kernel.)"},{"metadata":{"_cell_guid":"141a62bd-772f-4a68-94f9-c8fc14acacc7","_uuid":"d1db97182f30004f0fce30b8b29a1bc757296a94"},"cell_type":"markdown","source":"To train my stacker, I use logit transformations of each base model's validation predictions.  IMO this is the most straightforward way of doing it: my stacking model is a logistic regression, which takes a linear combination of its inputs and does a logistic (inverse logit) transformation on the result. It is essentially interpreting its inputs to be in units of logit (log odds) and then converting the result back into units of probability.\n\nBut you could also experiment with using other kinds of inputs.  Although the base models generate results that are nominally expressed as probabilities, they are mostly optimized for rank (i.e., using AUC) rather than probability value.  So in a sense they aren't *really* probabilities and maybe don't deserve to be treated as probabilities by having log odds taken.  Also, logit takes extreme predictions at face value, but maybe you don't trust the base models enough to believe their extreme predictions: maybe you want some regularization of the base model predicitons, and one way to get that regularization is to take the raw probabilities, or the ranks, instead of the logits."},{"metadata":{"collapsed":true,"_uuid":"f6e911d2bf71e23890211b63860d915ee263bf90","trusted":true},"cell_type":"code","source":"model_order = [m for m in base_models]  # To make sure order is consistent when converting to array","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"81a2a11d-fa38-4aa2-8b28-8c97fd72a1b0","_uuid":"8ba754b121183e83d1542c558150ae1a3bb8313a","trusted":true},"cell_type":"code","source":"cvdata = pd.DataFrame( { \n    m:pd.read_csv(cvfiles[m])['is_attributed'].clip(almost_zero,almost_one).apply(logit) \n    for m in base_models\n    } )\nX_train = np.array(cvdata[model_order])\ny_train = pd.read_pickle(VAL_FILE)['is_attributed']  # relies on validation cases being in same order","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"0c713caf-2ee7-432a-8dbd-d10590ce0fd9","_uuid":"ec8a93d77e9063ab6aad20a6c910b54b9798a425","trusted":true},"cell_type":"code","source":"cvdata.head()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"9765374d-f5df-414b-89fd-8ca186453de9","_uuid":"b071581ff35f4900a110c342b275390d1cd4d225","trusted":true},"cell_type":"code","source":"cvdata.corr()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"b3ce0323-0e5f-409e-b404-7c14b616d769","_uuid":"7b8e096669df5ae976feaef5a3165b60de08105f","trusted":true},"cell_type":"code","source":"stack_model = LogisticRegression()\nstack_model.fit(X_train, y_train)\nstack_model.coef_","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"f1981258-f50c-47cb-9417-7a52f97255be","_uuid":"63faf7a8a6af88edc671ddffa34a20b131dfc946"},"cell_type":"markdown","source":"Note that the evaluation criterion for this competition (AUC) depends only on rank. Therefore:\n1. Any linear transformation applied to the coefficients won't affect the score of the result.\n2. We don't care much about the value of the intercept term (not shown above).\n3. It doesn't matter whether the coefficients sum to 1.\n4. If we normalize the coefficients to sum to 1, the result has the same interpretation as blending weights."},{"metadata":{"_cell_guid":"c4531d65-0675-446c-80bb-2d21c8dd96f2","_uuid":"17bd987ce19a889c637bd27d48e3bafc60582ee2","trusted":true},"cell_type":"code","source":"weights = stack_model.coef_/stack_model.coef_.sum()\ncolumns = cvdata[model_order].columns\nscores = [ roc_auc_score( y_train, expit(cvdata[c]) )  for c in columns ]\nnames = [ base_models[c] for c in columns ]\nlb = [ lbscores[c] for c in columns ]\npd.DataFrame( data={'LB score': lb, 'CV score':scores, 'weight':weights.reshape(-1)}, index=names )","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"2c08fb71-6977-4683-a83b-00be98a9afa5","_uuid":"6551c01befb1d8ce2426f6114576fa34f8d2f9e0","trusted":true},"cell_type":"code","source":"print(  'Stacker score: ', roc_auc_score( y_train, stack_model.predict_proba(X_train)[:,1] )  )","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"e30502b3-9ec5-41cc-8d1c-586b06638cc7","_uuid":"1026a2476bcd384367988ba242d03202cf4311c9"},"cell_type":"markdown","source":"Take the stacker score with a grain of salt. It's based on a directly optimized fit, so it's kind of like training set performance. (The validation set is a training set for the stacker.) You can't expect the leaderboard score to improve as much (relative to single model performance) as the stacker score does.  And when you make changes, you shouldn't expect a close relationship between changes in the stacker score and changes in leaderboard score.  (One would expect them typically to move in the same direction, but there are a lot of exceptions.)  The stacker score will almost always be higher than the individual model CV scores, since the stacker could have chosen to use just one model.  (Not quite always, because the stacker is optimizing log likelihood rather than AUC.)  But the leaderbaord score of the stacked model is not necessarily going to be higher than best individual model leaderbaord score (though it usually is, which is why we do stacking).  "},{"metadata":{"collapsed":true,"_cell_guid":"638c438d-f89a-4f16-862e-1963cab66982","_uuid":"6695317f3b2022328774893364f04263a8236329","trusted":true},"cell_type":"code","source":"final_sub = pd.DataFrame()\nsubs = {m:pd.read_csv(subfiles[m]).rename({'is_attributed':m},axis=1) for m in base_models}\nfirst_model = list(base_models.keys())[0]\nfinal_sub['click_id'] = subs[first_model]['click_id']","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"2aa5ec3e-8a7d-47a9-8e42-15daf7d2f970","_uuid":"d89f2ac0df48446ddd04b62edaa77f8a2a1acef4","trusted":true},"cell_type":"code","source":"df = subs[first_model]\nfor m in base_models:\n    if m != first_model:\n        df = df.merge(subs[m], on='click_id')  # being careful in case clicks are in different order\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"a03efa79-dea6-495d-a071-27543000e1a5","_uuid":"3fa0b547b3435bb5f650a03eda1ff2e49238cef3","trusted":true},"cell_type":"code","source":"X_test = np.array( df.drop(['click_id'],axis=1)[model_order].clip(almost_zero,almost_one).apply(logit) )\nfinal_sub['is_attributed'] = stack_model.predict_proba(X_test)[:,1]\nfinal_sub.head(10)","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"_cell_guid":"c29b10eb-1897-4b68-af20-2f600652d140","_uuid":"bd2cdb0cce3eaf2f6cfee9fe2e3a25d1aacda47c","trusted":true},"cell_type":"code","source":"final_sub.to_csv(\"sub_stacked.csv\", index=False, float_format='%.9f')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"e211cea1aa15e73c57d9980b5bddc26bc66b1d87"},"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}