{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport shap\nfrom sklearn.pipeline import Pipeline \nfrom sklearn.ensemble import GradientBoostingRegressor\nfrom sklearn.model_selection import cross_val_score, GridSearchCV\nfrom sklearn.impute import SimpleImputer\nfrom sklearn.linear_model import LinearRegression\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-07-01T03:54:27.7161Z","iopub.execute_input":"2022-07-01T03:54:27.71655Z","iopub.status.idle":"2022-07-01T03:54:31.990433Z","shell.execute_reply.started":"2022-07-01T03:54:27.71646Z","shell.execute_reply":"2022-07-01T03:54:31.989221Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# <center>SHAP: all you need to know.</center>","metadata":{}},{"cell_type":"markdown","source":"The idea is to build a simple dummy model on this housing data set and use SHAP of this model to understand how much the features impact the model output. This will enable us to understand how SHAP values can be useful in explaining why we get such results when using black box ML algorithms. \n\nThis notebook is divided in 4 parts described below\n\n<a id=\"toc\"></a>\n* [1. Load data ](#1)<br>\n* [2. A simple model](#2)<br>\n    * pd.get_dummies → sklearn SimpleImputer → GradientBoostingRegressor\n* [3. Explain Model with SHAP](#3)<br>\n    * [3.1. partial dependence plots: how a specific feature impacts the model output](#3.1.)<br>\n    * [3.2. Understand a single prediction](#3.2.)<br>\n* [4. Feature Selection with SHAP](#4)<br>\n    * [4.1. feature importance plots](#4.1.)\n    * [4.2. final model with most important features](#4.2.)\n\n\n\n\n\n\n\n    ","metadata":{}},{"cell_type":"markdown","source":"<a id=\"1\"></a>\n# **<center><span style=\"color:#FF7B5F;\">1. Load data</span></center>**\nLet's first start by loading the data using pandas.","metadata":{}},{"cell_type":"code","source":"# first, let's increase the number of columns that pandas can display since the data set has a lot of columns \npd.options.display.max_columns=200\n\ntrain = pd.read_csv('/kaggle/input/house-prices-advanced-regression-techniques/train.csv')\ntest = pd.read_csv('/kaggle/input/house-prices-advanced-regression-techniques/test.csv')\n\n#and set the random_state for later\nrandom_state = 0\n\n# ... and print the shape and the first lines of the training set\nprint(train.shape)\nprint(test.shape)\ntrain.head()","metadata":{"execution":{"iopub.status.busy":"2022-07-01T03:54:31.99209Z","iopub.execute_input":"2022-07-01T03:54:31.992728Z","iopub.status.idle":"2022-07-01T03:54:32.13212Z","shell.execute_reply.started":"2022-07-01T03:54:31.992694Z","shell.execute_reply":"2022-07-01T03:54:32.130899Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train = train.drop(['Id', 'SalePrice'], axis=1)\ny_train = train['SalePrice']\nX_test = test.drop('Id', axis=1)","metadata":{"execution":{"iopub.status.busy":"2022-07-01T03:54:32.13366Z","iopub.execute_input":"2022-07-01T03:54:32.134074Z","iopub.status.idle":"2022-07-01T03:54:32.150963Z","shell.execute_reply.started":"2022-07-01T03:54:32.134032Z","shell.execute_reply":"2022-07-01T03:54:32.149831Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"2\"></a>\n# **<center><span style=\"color:#FF7B5F;\">2. A simple model</span></center>**","metadata":{}},{"cell_type":"markdown","source":"Let's create a simple model which cn predict the price of a house based on the data available. We'll then use shap to understand why the model predicts such prices. \n\nLets begin by creating dummies with our data: ","metadata":{}},{"cell_type":"code","source":"X_train_dummified = pd.get_dummies(X_train)","metadata":{"execution":{"iopub.status.busy":"2022-07-01T03:54:32.154718Z","iopub.execute_input":"2022-07-01T03:54:32.155343Z","iopub.status.idle":"2022-07-01T03:54:32.209449Z","shell.execute_reply.started":"2022-07-01T03:54:32.155307Z","shell.execute_reply":"2022-07-01T03:54:32.208447Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Then, let's use sklearn SimpleImputer so we can get rid of null values.","metadata":{}},{"cell_type":"code","source":"imputer = SimpleImputer()\nX_train_dummified_imputed = pd.DataFrame(imputer.fit_transform(X_train_dummified), columns = X_train_dummified.columns)","metadata":{"execution":{"iopub.status.busy":"2022-07-01T03:54:32.211081Z","iopub.execute_input":"2022-07-01T03:54:32.211782Z","iopub.status.idle":"2022-07-01T03:54:32.247325Z","shell.execute_reply.started":"2022-07-01T03:54:32.211739Z","shell.execute_reply":"2022-07-01T03:54:32.246247Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Finally, let's use sklearn's GradientBoostingRegressor as our estimator. This will enable us to then use the shap library to understand this model and see which features are important for the prediction. ","metadata":{}},{"cell_type":"code","source":"# estimator = GradientBoostingRegressor(random_state = random_state) # LinearRegression()\nestimator = GradientBoostingRegressor(random_state = random_state)\nestimator.fit(X_train_dummified_imputed, y_train)","metadata":{"execution":{"iopub.status.busy":"2022-07-01T03:54:32.249125Z","iopub.execute_input":"2022-07-01T03:54:32.249512Z","iopub.status.idle":"2022-07-01T03:54:33.342863Z","shell.execute_reply.started":"2022-07-01T03:54:32.24948Z","shell.execute_reply":"2022-07-01T03:54:33.341893Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scores = cross_val_score(estimator, X_train_dummified_imputed, y_train, cv=5, scoring= 'r2')\nprint(\"Average CV score:\", scores.mean())","metadata":{"execution":{"iopub.status.busy":"2022-07-01T03:54:33.344122Z","iopub.execute_input":"2022-07-01T03:54:33.344457Z","iopub.status.idle":"2022-07-01T03:54:37.650253Z","shell.execute_reply.started":"2022-07-01T03:54:33.344429Z","shell.execute_reply":"2022-07-01T03:54:37.649036Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"3\"></a>\n# **<center><span style=\"color:#FF7B5F;\">3. Model explanation with SHAP</span></center>**","metadata":{}},{"cell_type":"markdown","source":"<a id=\"3.1.\"></a>\n## **<center><span style=\"color:#FF7B5F;\">3.1. See how a specific feature impacts the output</span></center>**","metadata":{}},{"cell_type":"markdown","source":"One question that is often considered when building ML models is how a specific feature impacts our model. An answer to this are Partial Dependence Plots (aka PDPs). PDPs enable us to understand how much the output changes when a feature changes. \n\nHere, for example: why did the model predict such a price estimation?   ","metadata":{}},{"cell_type":"code","source":"shap.partial_dependence_plot(\n    \"OverallQual\"\n    , estimator.predict\n    , X_train_dummified_imputed\n    , ice=False\n    , model_expected_value=True\n    , feature_expected_value=True\n)","metadata":{"execution":{"iopub.status.busy":"2022-07-01T03:54:37.651508Z","iopub.execute_input":"2022-07-01T03:54:37.651812Z","iopub.status.idle":"2022-07-01T03:54:39.000521Z","shell.execute_reply.started":"2022-07-01T03:54:37.651784Z","shell.execute_reply":"2022-07-01T03:54:38.999203Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here, the partial dependence plot indicates that as OverallQual increases the output of the ML model (i.e. Estimated Price) also increases. This can help us see if our model overfits too, just like the example described below: ","metadata":{}},{"cell_type":"code","source":"shap.partial_dependence_plot(\n    \"BsmtUnfSF\"\n    , estimator.predict\n    , X_train_dummified_imputed\n    , ice=False\n    , model_expected_value=True\n    , feature_expected_value=True\n)\n","metadata":{"execution":{"iopub.status.busy":"2022-07-01T03:54:39.002581Z","iopub.execute_input":"2022-07-01T03:54:39.004455Z","iopub.status.idle":"2022-07-01T03:54:40.377246Z","shell.execute_reply.started":"2022-07-01T03:54:39.004374Z","shell.execute_reply":"2022-07-01T03:54:40.37597Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We see here that there might be something weird with this value: why is our model decreasing the estimated price when BsmtUnfSF is between 1700 and 1800 and then re-increasing it. This is an indication of overfitting. ","metadata":{}},{"cell_type":"markdown","source":"<a id=\"3.2.\"></a>\n## **<center><span style=\"color:#FF7B5F;\">3.2. Visualizing SHAP values</span></center>**","metadata":{}},{"cell_type":"markdown","source":"With the help of the shap library, we can also visualize why we had such results for a single row of the data. ","metadata":{}},{"cell_type":"code","source":"X100 = shap.utils.sample(X_train_dummified_imputed, 100)","metadata":{"execution":{"iopub.status.busy":"2022-07-01T03:54:40.382927Z","iopub.execute_input":"2022-07-01T03:54:40.383993Z","iopub.status.idle":"2022-07-01T03:54:40.392926Z","shell.execute_reply.started":"2022-07-01T03:54:40.383933Z","shell.execute_reply":"2022-07-01T03:54:40.391453Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"explainer = shap.Explainer(estimator.predict, X100)\nshap_values = explainer(X100)","metadata":{"execution":{"iopub.status.busy":"2022-07-01T03:54:40.396006Z","iopub.execute_input":"2022-07-01T03:54:40.39725Z","iopub.status.idle":"2022-07-01T03:55:00.598478Z","shell.execute_reply.started":"2022-07-01T03:54:40.39715Z","shell.execute_reply":"2022-07-01T03:55:00.597136Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_index = 12\nshap.partial_dependence_plot(\n    \"OverallQual\"\n    , estimator.predict\n    , X100\n    , ice=False\n    , model_expected_value=True\n    , feature_expected_value=True\n    , shap_values=shap_values[sample_index:sample_index+1,:]\n)","metadata":{"execution":{"iopub.status.busy":"2022-07-01T03:55:00.600136Z","iopub.execute_input":"2022-07-01T03:55:00.600609Z","iopub.status.idle":"2022-07-01T03:55:01.600304Z","shell.execute_reply.started":"2022-07-01T03:55:00.600563Z","shell.execute_reply":"2022-07-01T03:55:01.598757Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shap.plots.waterfall(shap_values[sample_index], max_display=14)","metadata":{"execution":{"iopub.status.busy":"2022-07-01T03:55:01.602559Z","iopub.execute_input":"2022-07-01T03:55:01.603177Z","iopub.status.idle":"2022-07-01T03:55:02.75416Z","shell.execute_reply.started":"2022-07-01T03:55:01.603079Z","shell.execute_reply":"2022-07-01T03:55:02.752763Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"With these graphs you can now understand why the model preicted such an output: here the building was built in 1925 which decreases the model prediction but the overall quality is high (7) which increases the model prediction etc. ","metadata":{}},{"cell_type":"markdown","source":"<a id=\"4\"></a>\n# **<center><span style=\"color:#FF7B5F;\">4. Feature Selection with SHAP</span></center>**","metadata":{}},{"cell_type":"markdown","source":"<a id=\"4.1.\"></a>\n## **<center><span style=\"color:#FF7B5F;\">4.1. Feature importance plots</span></center>**\n\n\nThe two most useful feature importance plots that come with the shap library are the barplot and the beeswarm plot. \n- barplot plot: gives the importance of features in order \n- beeswarm plot: gives the importance of feature + color on how the value of each feature impacts the output","metadata":{}},{"cell_type":"code","source":"shap.summary_plot(shap_values, max_display=15, show=False)","metadata":{"execution":{"iopub.status.busy":"2022-07-01T03:55:02.755625Z","iopub.execute_input":"2022-07-01T03:55:02.756912Z","iopub.status.idle":"2022-07-01T03:55:03.431572Z","shell.execute_reply.started":"2022-07-01T03:55:02.756855Z","shell.execute_reply":"2022-07-01T03:55:03.430449Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shap.summary_plot(shap_values, max_display=15, show=False, plot_type='bar')","metadata":{"execution":{"iopub.status.busy":"2022-07-01T03:55:03.43327Z","iopub.execute_input":"2022-07-01T03:55:03.433931Z","iopub.status.idle":"2022-07-01T03:55:03.757633Z","shell.execute_reply.started":"2022-07-01T03:55:03.433889Z","shell.execute_reply":"2022-07-01T03:55:03.756438Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"feature_importance = pd.DataFrame({'name': X100.columns, 'importance': shap_values.abs.sum(0).values})\nfeature_importance= feature_importance.sort_values(by='importance', ascending=False).reset_index(drop=True)\nfeature_importance[feature_importance['importance']>0]","metadata":{"execution":{"iopub.status.busy":"2022-07-01T03:55:03.758835Z","iopub.execute_input":"2022-07-01T03:55:03.759128Z","iopub.status.idle":"2022-07-01T03:55:03.781501Z","shell.execute_reply.started":"2022-07-01T03:55:03.759101Z","shell.execute_reply":"2022-07-01T03:55:03.780424Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"4.2.\"></a>\n## **<center><span style=\"color:#FF7B5F;\"> 4.2. Only keeping top features</span></center>**\n","metadata":{}},{"cell_type":"code","source":"X_train_final = X_train_dummified_imputed[feature_importance.head(25)[\"name\"]]\n\nfinal_estimator = GradientBoostingRegressor(random_state = random_state)\nfinal_estimator.fit(X_train_final, y_train)","metadata":{"execution":{"iopub.status.busy":"2022-07-01T03:55:03.78309Z","iopub.execute_input":"2022-07-01T03:55:03.783401Z","iopub.status.idle":"2022-07-01T03:55:04.220018Z","shell.execute_reply.started":"2022-07-01T03:55:03.783357Z","shell.execute_reply":"2022-07-01T03:55:04.218741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scores = cross_val_score(final_estimator, X_train_final, y_train, cv=5, scoring= 'r2')\nprint(\"Average CV score:\", scores.mean())","metadata":{"execution":{"iopub.status.busy":"2022-07-01T03:55:04.221614Z","iopub.execute_input":"2022-07-01T03:55:04.221967Z","iopub.status.idle":"2022-07-01T03:55:06.017727Z","shell.execute_reply.started":"2022-07-01T03:55:04.221937Z","shell.execute_reply":"2022-07-01T03:55:06.016458Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We see here that only keeping the top features allowed us to increase our CV score while dividing by 10 the number of features of the model (from 280 to 25!!!). Of course, hyperparameter tuning can now be used to further improve this performance, now that we know which features we want to keep! ","metadata":{}},{"cell_type":"markdown","source":"Hope you liked the article. If you need, feel free to upvote this notebook! Happy learning!!","metadata":{}},{"cell_type":"markdown","source":"## Appendix: useful functions to handle null values if you don't want to use sklearn imputer","metadata":{}},{"cell_type":"code","source":"def calculate_missing_percentage(df):\n    percent_missing = np.round(df.isnull().sum() * 100 / len(df),2)\n    missing_value_df = pd.DataFrame({'percent_missing': percent_missing})\n    return missing_value_df.sort_values(by='percent_missing', ascending=False)\n\nstats_missing = calculate_missing_percentage(train)\nstats_missing[stats_missing['percent_missing']>30]","metadata":{"execution":{"iopub.status.busy":"2022-07-01T03:55:06.019284Z","iopub.execute_input":"2022-07-01T03:55:06.020288Z","iopub.status.idle":"2022-07-01T03:55:06.044732Z","shell.execute_reply.started":"2022-07-01T03:55:06.020239Z","shell.execute_reply":"2022-07-01T03:55:06.043339Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def handle_null_values(df, drop_null_column_threshold=30):\n    stats_missing = calculate_missing_percentage(df)\n    df = df.drop(stats_missing[stats_missing['percent_missing']>drop_null_column_threshold].index, axis=1)\n    df = df.fillna(df.median())\n    return df\ndf_train_preprocessed = handle_null_values(train)","metadata":{"execution":{"iopub.status.busy":"2022-07-01T03:55:06.046684Z","iopub.execute_input":"2022-07-01T03:55:06.047097Z","iopub.status.idle":"2022-07-01T03:55:06.086344Z","shell.execute_reply.started":"2022-07-01T03:55:06.047056Z","shell.execute_reply":"2022-07-01T03:55:06.084037Z"},"trusted":true},"execution_count":null,"outputs":[]}]}