{"cells":[{"metadata":{},"cell_type":"markdown","source":"This kernel shows how to use [Shapley values](https://www.google.com/search?q=Shapley+values) for feature importance scoring. This is not a new thing and has been used in several competitions in the past e.g. [a old kernel.](https://www.kaggle.com/hmendonca/lightgbm-predictions-explained-with-shap-0-796). Shapley values are not the default metric from LGBM but can easily be implemented with the [SHAP library](https://github.com/slundberg/shap), which is also very useful for model explainability.   \nOnce you start experimenting and adding many features to your model, these techniques became crucial to select only the best features and avoid over-fitting."},{"metadata":{},"cell_type":"markdown","source":"### ASHRAE - Great Energy Predictor III\n\nI'll summarise below the insights from my brief data analysis and whatever is shared in public discussions.\nThis is still highly under progress and will be updated when possible.\n\nQuoting Chris Balbach (our Competition Host):\n> Consider the scenario laid out in this diagram. This competition simulates the modelling challenge presented at the end of the timeline when measured energy and weather conditions are known and the adjusted baseline energy must be calculated. \n![Energy Consumption/Demand vs Time](https://www.mdpi.com/make/make-01-00056/article_deploy/html/images/make-01-00056-g001-550.jpg)\n> For further details, I recommend reading [this paper](https://www.mdpi.com/2504-4990/1/3/56) by Clayton Miller.\n\nWe are actually using data from 2016 to predict the demand for both 2017 and 2018, which is a hard task. However, other public kernels have showed that you can give a relatively good estimate even with a very simple linear model."},{"metadata":{"_kg_hide-input":true,"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os, gc\nimport random\nimport datetime\n\nfrom tqdm import tqdm_notebook as tqdm\n\n# matplotlib and seaborn for plotting\nimport matplotlib\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\n# import plotly.offline as py\n# py.init_notebook_mode(connected=True)\n# from plotly.offline import init_notebook_mode, iplot\n# init_notebook_mode(connected=True)\n# import plotly.graph_objs as go\n# import plotly.offline as offline\n# offline.init_notebook_mode()\n\nfrom pandas.plotting import register_matplotlib_converters\nregister_matplotlib_converters()\n\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.model_selection import KFold, StratifiedKFold\nfrom sklearn import preprocessing\nfrom sklearn.metrics import mean_squared_error, mean_squared_log_error\n\nimport lightgbm as lgb\nimport shap","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"_kg_hide-output":true,"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"path = '../input/ashrae-energy-prediction'\n# Input data files are available in the \"../input/\" directory.\nfor dirname, _, filenames in os.walk(path):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load data"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n# unimportant features (see importance below)\nunimportant_cols = ['wind_direction', 'wind_speed', 'sea_level_pressure']\ntarget = 'meter_reading'\n\ndef load_data(source='train', path=path):\n    ''' load and merge all tables '''\n    assert source in ['train', 'test']\n    \n    building = pd.read_csv(f'{path}/building_metadata.csv', dtype={'building_id':np.uint16, 'site_id':np.uint8})\n    weather  = pd.read_csv(f'{path}/weather_{source}.csv', parse_dates=['timestamp'],\n                                                           dtype={'site_id':np.uint8, 'air_temperature':np.float16,\n                                                                  'cloud_coverage':np.float16, 'dew_temperature':np.float16,\n                                                                  'precip_depth_1_hr':np.float16},\n                                                           usecols=lambda c: c not in unimportant_cols)\n    df = pd.read_csv(f'{path}/{source}.csv', dtype={'building_id':np.uint16, 'meter':np.uint8}, parse_dates=['timestamp'])\n    df = df.merge(building, on='building_id', how='left')\n    df = df.merge(weather, on=['site_id', 'timestamp'], how='left')\n    return df\n\n# load and display some samples\ntrain = load_data('train')\ntrain.sample(5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Test data"},{"metadata":{"trusted":true},"cell_type":"code","source":"test = load_data('test')\ntest.sample(5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# EDA and sample statistics"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"print(f'Training from {train.timestamp.min()} to {train.timestamp.max()}, and predicting from {test.timestamp.min()} to {test.timestamp.max()}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Building types on each site location"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"building_counts = train.groupby(['primary_use', 'site_id']).building_id.nunique().to_frame('counts')\nbuilding_counts = building_counts.reset_index().pivot(index='primary_use', columns='site_id', values='counts')\n\nfig, ax = plt.subplots(figsize=(16,8))\n_ = sns.heatmap(building_counts, annot=True, cmap='Reds',\n                xticklabels=building_counts.columns.values,\n                yticklabels=building_counts.index.values)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Distribution of Meter readings"},{"metadata":{"trusted":true},"cell_type":"code","source":"# target's log-log histogram:\nax = np.log1p(train.meter_reading).hist()\nax.set_yscale('log')\n\n# describe raw values first\ntrain.meter_reading.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### distribution in the types of meters"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# check the distribution in the types of meters\nmeters = train.groupby('building_id').meter.nunique().to_frame()\nax = meters.hist()\n_ = ax[0][0].set_title('Distribution of types of meters\\n{0:electricity, 1:water, 2:steam, 3:hotwater}') # from the official starter kernel\n# from the graphs it looks like steam and hotwater are reversed (e.g.: 3:steam, 2:hotwater) but that shouldn't make any difference to the model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### display a single time series (notice measurement errors and discontinuities)"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"building_id = 1258  # a building with all 4 meters:  meters[meters.meter == 4]\nmeters = train[train['building_id'] == building_id].meter.nunique()\n\nfor meter in range(meters):\n    fig, ax = plt.subplots()\n    plt.title(f'Building {building_id} Meter {meter}')\n    ax2 = ax.twinx()\n    # plot meter_reading\n    idx = (train['building_id'] == building_id) & (train['meter'] == meter)\n    dates = matplotlib.dates.date2num(train.loc[idx, 'timestamp'])\n    ax2.plot_date(dates, train.loc[idx, 'meter_reading'], '-', label='meter_reading', alpha=0.8)\n    # plot air_temperature\n    dates = matplotlib.dates.date2num(train.loc[train['building_id'] == building_id, 'timestamp'])\n    ax.plot_date(dates, train.loc[train['building_id'] == building_id, 'air_temperature'], '.', color='tab:cyan', label='air_temperature')\n    ax.set_ylabel('air_temperature'); ax2.set_ylabel('meter_reading')\n    ax.legend(loc='upper left'); ax2.legend(loc='upper right')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### now let's see what's the expected prediction in the test set for the same building"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"meter = 1 # pick a meter\n\ntrain_sample = train[(train['building_id'] == building_id) & (train['meter'] == meter)]  # same train sample as above\n\ntest['meter_reading'] = 0.0\ntest_sample = test[(test['building_id'] == building_id) & (test['meter'] == meter)]  # and the same meter in the test set\n\nfig, ax = plt.subplots(figsize=(16,4))\nplt.title(f'Building {building_id} Meter {meter}')\nax.xaxis.set_tick_params(rotation=30, labelsize=10)\nax2 = ax.twinx()\n\n# plot training sample\ndates = matplotlib.dates.date2num(train_sample['timestamp'])\nax2.plot_date(dates, train_sample['meter_reading'], '-', label='train', alpha=0.8)\nax.plot_date(dates, train_sample['air_temperature'], '.', color='tab:cyan', label='air_temperature')\n\n# plot test sample\ndates = matplotlib.dates.date2num(test_sample['timestamp'])\nax2.plot_date(dates, test_sample['meter_reading'], '*', label='test', alpha=0.8)\nax.plot_date(dates, test_sample['air_temperature'], '.', color='tab:cyan', label='air_temperature')\n\nax.set_ylabel('air_temperature'); ax2.set_ylabel('meter_reading')\nax.legend(loc='upper left'); ax2.legend(loc='upper right')\n\ndel train_sample; del test_sample; del dates","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Check for missing data"},{"metadata":{"trusted":true},"cell_type":"code","source":"# some feature stats\ntrain.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# the counts above expose the missing data (Should we drop or refill the missing data?)\nprint(\"Ratio of available data (not NAN's):\")\ndata_ratios = train.count()/len(train)\ndata_ratios","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Is the same happening in the test set? Yes\nprint(\"Ratio of available data (not NAN's):\")\ntest.count()/len(test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# we can refill with averages\ntrain.loc[:, data_ratios < 1.0].mean()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Preprocess data"},{"metadata":{"trusted":true},"cell_type":"code","source":"class ASHRAE3Preprocessor(object):\n    @classmethod\n    def fit(cls, df, data_ratios=data_ratios):\n        cls.avgs = df.loc[:,data_ratios < 1.0].mean()\n        cls.pu_le = LabelEncoder()\n        cls.pu_le.fit(df[\"primary_use\"])\n\n    @classmethod\n    def transform(cls, df):\n        df = df.fillna(cls.avgs) # refill NAN with averages\n        df['primary_use'] = np.uint8(cls.pu_le.transform(df['primary_use']))  # encode labels\n\n        # expand datetime into its components\n        df['hour'] = np.uint8(df['timestamp'].dt.hour)\n        df['day'] = np.uint8(df['timestamp'].dt.day)\n        df['weekday'] = np.uint8(df['timestamp'].dt.weekday)\n        df['month'] = np.uint8(df['timestamp'].dt.month)\n        df['year'] = np.uint8(df['timestamp'].dt.year-2000)\n        \n        # parse and cast columns to a smaller type\n        df.rename(columns={\"square_feet\": \"log_square_feet\"}, inplace=True)\n        df['log_square_feet'] = np.float16(np.log(df['log_square_feet']))\n        df['year_built'] = np.uint8(df['year_built']-1900)\n        df['floor_count'] = np.uint8(df['floor_count'])\n        \n        # remove redundant columns\n        for col in df.columns:\n            if col in ['timestamp', 'row_id']:\n                del df[col]\n    \n        # extract target column\n        if 'meter_reading' in df.columns:\n            df['meter_reading'] = np.log1p(df['meter_reading']).astype(np.float32) # comp metric uses log errors\n\n        return df\n        \nASHRAE3Preprocessor.fit(train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = ASHRAE3Preprocessor.transform(train)\ntrain.sample(7)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.dtypes","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Feature ranked correlation"},{"metadata":{"trusted":true,"_kg_hide-output":false,"_kg_hide-input":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(16,8))\n# use a ranked correlation to catch nonlinearities (linear correlation is not important for boosted  trees e.g. LightGBM)\ncorr = train[[col for col in train.columns if col != 'year']].sample(100100).corr(method='spearman')\n_ = sns.heatmap(corr, annot=True,\n                xticklabels=corr.columns.values,\n                yticklabels=corr.columns.values)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Train K folds"},{"metadata":{"trusted":true},"cell_type":"code","source":"# force the model to use the weather data instead of dates, to avoid overfitting to the past history\nfeatures = [col for col in train.columns if col not in [target, 'year', 'month', 'day']]\n# sample features\ntrain[features].sample(5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"folds = 4\nseed = 42\nkf = StratifiedKFold(n_splits=folds, shuffle=True, random_state=seed)\nmodels = []\nshap_values = np.zeros(train[features].shape)\nshap_sampling = 125000  # reduce compute cost\noof_pred = np.zeros(train.shape[0])  # out of fold predictions\n\n## stratify data by building_id\nfor i, (tr_idx, val_idx) in tqdm(enumerate(kf.split(train, train['building_id'])), total=folds):\n    def fit_regressor(tr_idx, val_idx): # memory closure\n        tr_x, tr_y = train[features].iloc[tr_idx],  train[target].iloc[tr_idx]\n        vl_x, vl_y = train[features].iloc[val_idx], train[target].iloc[val_idx]\n        print({'fold':i, 'train size':len(tr_x), 'eval size':len(vl_x)})\n\n        tr_data = lgb.Dataset(tr_x, label=tr_y)\n        vl_data = lgb.Dataset(vl_x, label=vl_y)  \n        clf = lgb.LGBMRegressor(n_estimators=900,\n                                learning_rate=0.33,\n                                feature_fraction=0.9,\n                                subsample=0.25,  # batches of 25% of the data\n                                subsample_freq=1,\n                                num_leaves=20,\n                                lambda_l1=0.5,  # regularisation\n                                lambda_l2=0.5,\n                                seed=i,   # seed diversification\n                                metric='rmse')\n        clf.fit(tr_x, tr_y,\n                eval_set=[(vl_x, vl_y)],\n#                 early_stopping_rounds=50,\n                verbose=150)\n        # sample shapley values\n        fold_importance = shap.TreeExplainer(clf).shap_values(vl_x[:shap_sampling])\n        # out of fold predictions\n        valid_prediticion = clf.predict(vl_x, num_iteration=clf.best_iteration_)\n        oof_loss = np.sqrt(mean_squared_error(vl_y, valid_prediticion)) # target is already in log scale\n        print(f'Fold:{i} RMSLE: {oof_loss:.4f}')\n        return clf, fold_importance, valid_prediticion\n\n    clf, shap_values[val_idx[:shap_sampling]], oof_pred[val_idx] = fit_regressor(tr_idx, val_idx)\n    models.append(clf)\n    \ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"oof_loss = np.sqrt(mean_squared_error(train[target], oof_pred)) # target is already in log scale\nprint(f'OOF RMSLE: {oof_loss:.4f}')\n\n# save out of fold predictions\npd.DataFrame(oof_pred).to_csv(f'oof_preds{oof_loss:.4f}.csv.bz')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Feature importance"},{"metadata":{},"cell_type":"markdown","source":"### Classical LGBM feature importance (gain on first fold)"},{"metadata":{"trusted":true},"cell_type":"code","source":"_ = lgb.plot_importance(models[0], importance_type='gain')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Shapley Values feature importance across all folds"},{"metadata":{"trusted":true},"cell_type":"code","source":"shap.summary_plot(shap_values, train[features], plot_type=\"bar\")","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"ma_shap = pd.DataFrame(sorted(zip(abs(shap_values).mean(axis=0), features), reverse=True),\n                       columns=['Mean Abs Shapley', 'Feature']).set_index('Feature')\n# fig, ax = plt.subplots(figsize=(2,6))\n# _ = sns.heatmap(ma_shap, annot=True, cmap='Blues', fmt='.06f')\nma_shap","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Model explainability with SHAP"},{"metadata":{"trusted":true},"cell_type":"code","source":"# visualize the first prediction's explanation (use matplotlib=True to avoid Javascript)\nshap.force_plot(shap.TreeExplainer(models[0]).expected_value, shap_values[0,:], train[features].iloc[0,:], matplotlib=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Check prediction"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# load and pre-process test data\ntest = ASHRAE3Preprocessor.transform(test)\ntest.sample(5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### now let's revisit the same meter we had initially looked at, and check the predictions for each fold"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def recover_timestamp(x):\n    ''' reassemble timestamp using date components '''\n    return datetime.datetime.strptime(f'{x.year}-{x.month}-{x.day} {x.hour}', '%y-%m-%d %H')\n\nfig, ax = plt.subplots(figsize=(16,4))\nplt.title(f'Building {building_id} Meter {meter} on all {folds} prediction folds')\nax.xaxis.set_tick_params(rotation=30, labelsize=10)\nax2 = ax.twinx()\n\ntrain_sample = train[(train['building_id'] == building_id) & (train['meter'] == meter)]  # same training sample as before\ntest_sample = test[(test['building_id'] == building_id) & (test['meter'] == meter)]   # and the same meter in the test set\n\n# plot training sample\ndates = matplotlib.dates.date2num(train_sample[['year', 'month', 'day', 'hour']].apply(recover_timestamp, axis=1))\nax.plot_date(dates, train_sample['air_temperature'], '.', color='tab:cyan', label='air_temperature')\nax2.plot_date(dates, np.expm1(train_sample['meter_reading']), '-', color='tab:blue', label='train')\n\n# plot prediction sample\ndates = matplotlib.dates.date2num(test_sample[['year', 'month', 'day', 'hour']].apply(recover_timestamp, axis=1))\nax.plot_date(dates, test_sample['air_temperature'], '.', color='tab:cyan', label='air_temperature')\nfor i,model in enumerate(models):\n    ax2.plot_date(dates, np.expm1(model.predict(test_sample[features])), '-', label=f'prediction{i}', alpha=0.4)\n\nax.set_ylabel('air_temperature'); ax2.set_ylabel('meter_reading (+prediction)')\nax.legend(loc='upper left'); ax2.legend(loc='upper right')\n_ = plt.show()\n\ndel test_sample; del train_sample\n_ = gc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# # Check if all test buildings and meters are the same as in the training data\n# train_buildings = np.unique(train[['building_id', 'meter']].values, axis=0)\n# # del train; gc.collect()\n\n# test_buildings  = np.unique(test[['building_id', 'meter']].values, axis=0)\n\n# print(len(train_buildings), len(test_buildings))\n# [b for b in test_buildings if b not in train_buildings]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Test Inference and Submission"},{"metadata":{"trusted":true},"cell_type":"code","source":"# split test data into batches\nset_size = len(test)\niterations = 50\nbatch_size = set_size // iterations\n\nprint(set_size, iterations, batch_size)\nassert set_size == iterations * batch_size","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"meter_reading = []\nfor i in tqdm(range(iterations)):\n    pos = i*batch_size\n    fold_preds = [np.expm1(model.predict(test[features].iloc[pos : pos+batch_size])) for model in models]\n    meter_reading.extend(np.mean(fold_preds, axis=0))\n\nprint(len(meter_reading))\nassert len(meter_reading) == set_size","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Save submission"},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.read_csv(f'{path}/sample_submission.csv')\nsubmission['meter_reading'] = np.clip(meter_reading, a_min=0, a_max=None) # clip min at zero","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.to_csv('submission.csv', index=False)\nsubmission.head(9)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# prediction's log-log histogram:\nplt.yscale('log')\n_ = np.log1p(submission['meter_reading']).hist()\n\n# describe raw submission\nsubmission['meter_reading'].describe()","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.5"}},"nbformat":4,"nbformat_minor":1}