{"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":"markdown","source":"# Housing Price for Beginner using Basic Ridge, Lasso, Elastic","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-07-21T22:46:22.473143Z","iopub.execute_input":"2022-07-21T22:46:22.473736Z","iopub.status.idle":"2022-07-21T22:46:22.488558Z","shell.execute_reply.started":"2022-07-21T22:46:22.473620Z","shell.execute_reply":"2022-07-21T22:46:22.487287Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Acknowledgement \n\nNotebooks from which I inspired the most:\n- https://www.kaggle.com/houcembenmansour/house-price-prediction\n- https://www.kaggle.com/rbyron/simple-linear-regression-models\n- https://www.kaggle.com/apapiu/regularized-linear-models","metadata":{}},{"cell_type":"markdown","source":"# Objective:\nThe goal of this work is to build a model that can correctly predict `SalePrice`  \nThis notebook is targeted to give exposure to beginner (myself) to work with continuous target variable, we will only implement basic model such as:\n- Simple Linear Regression (OLS)\n- Ridge Regression (Regression with L2 Regularization)\n- Lasso Regression (Regression with L1 Regularization)\n- Elastic (Regression with combination of L1 and L2)\n\nWe will use RMSE as a metric","metadata":{}},{"cell_type":"code","source":"# Basic\nimport pandas as pd\nimport numpy as np\nimport random\n\n# Plots\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom scipy import stats\n\n# Misc.\nimport warnings\nwarnings.filterwarnings('ignore')","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:46:22.501136Z","iopub.execute_input":"2022-07-21T22:46:22.501678Z","iopub.status.idle":"2022-07-21T22:46:23.334671Z","shell.execute_reply.started":"2022-07-21T22:46:22.501643Z","shell.execute_reply":"2022-07-21T22:46:23.333508Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = pd.read_csv(\"../input/house-prices-advanced-regression-techniques/train.csv\")\ntest  = pd.read_csv(\"../input/house-prices-advanced-regression-techniques/test.csv\")\nsample= pd.read_csv(\"../input/house-prices-advanced-regression-techniques/sample_submission.csv\")","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:46:23.336530Z","iopub.execute_input":"2022-07-21T22:46:23.336978Z","iopub.status.idle":"2022-07-21T22:46:23.430699Z","shell.execute_reply.started":"2022-07-21T22:46:23.336930Z","shell.execute_reply":"2022-07-21T22:46:23.429597Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# EDA\n\nThis section will explore the data, goal:\n- get to know the dataset, how many features etc.\n- quick overview of feature correlation with dependent variable","metadata":{}},{"cell_type":"markdown","source":"## Getting to know data","metadata":{}},{"cell_type":"code","source":"train.head()","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:46:23.433072Z","iopub.execute_input":"2022-07-21T22:46:23.433551Z","iopub.status.idle":"2022-07-21T22:46:23.482497Z","shell.execute_reply.started":"2022-07-21T22:46:23.433506Z","shell.execute_reply":"2022-07-21T22:46:23.481291Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('Data shape: ', train.shape)\nprint('There are %d instances' %train.shape[0])\nprint('There are %d features' %train.shape[1])","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:46:23.484558Z","iopub.execute_input":"2022-07-21T22:46:23.485232Z","iopub.status.idle":"2022-07-21T22:46:23.493540Z","shell.execute_reply.started":"2022-07-21T22:46:23.485155Z","shell.execute_reply":"2022-07-21T22:46:23.492199Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Correlation heatmap\nGoal: see the correlation between each features, correlation range between -1 ant 1, the values near -1 or 1 shows stronge negative correlation or positive correlation respectively, while weak correlation indicated by values closer to zero. The goal is to see which features has high correlation (either positive or negative) with `SalePrice`","metadata":{}},{"cell_type":"code","source":"# Create correlation matrix\ncorrmat = train.corr()\n# corrmat\n\nplt.figure(figsize=(10, 10))\nax = sns.heatmap(corrmat, square=True, vmax=1, vmin=-1)\nax.set_title('Correlation Heatmap of Housing Pricing Train data')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:46:23.495095Z","iopub.execute_input":"2022-07-21T22:46:23.495890Z","iopub.status.idle":"2022-07-21T22:46:24.773667Z","shell.execute_reply.started":"2022-07-21T22:46:23.495841Z","shell.execute_reply":"2022-07-21T22:46:24.772387Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"For now we only need to focus on the last row `SalePrice`:\n- positively correlated : `OverallQual`, `GrLivArea`\n- negatively correlated : `None`\nThere are no variales that has high negative correlation with `SalePrice` \nLet see correlation values of those two variable in positive correlated and plot them to see what patterns they have","metadata":{}},{"cell_type":"code","source":"train['OverallQual']","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:46:24.775050Z","iopub.execute_input":"2022-07-21T22:46:24.775413Z","iopub.status.idle":"2022-07-21T22:46:24.785258Z","shell.execute_reply.started":"2022-07-21T22:46:24.775380Z","shell.execute_reply":"2022-07-21T22:46:24.784064Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Set seaborn theme\nsns.set(style='darkgrid', palette='muted')","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:46:24.786879Z","iopub.execute_input":"2022-07-21T22:46:24.787703Z","iopub.status.idle":"2022-07-21T22:46:24.796668Z","shell.execute_reply.started":"2022-07-21T22:46:24.787527Z","shell.execute_reply":"2022-07-21T22:46:24.795748Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Example of positive correlation\ntarget = 'SalePrice'\nvar1 = 'OverallQual'\nvar2 = 'GrLivArea'\n\n\nfig, (ax1, ax2)  = plt.subplots(1, 2, figsize=(12, 5), sharey=True)\nsns.boxplot(x=var1, y=target, data=train, ax=ax1)\nax1.set_title('Correlation values %.3f' %corrmat.loc[target, var1])\nsns.scatterplot(x=var2, y=target, data=train, ax=ax2)\nax2.set_title('Correlation values %.3f' %corrmat.loc[target, var2])\nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:46:24.798630Z","iopub.execute_input":"2022-07-21T22:46:24.799155Z","iopub.status.idle":"2022-07-21T22:46:25.601343Z","shell.execute_reply.started":"2022-07-21T22:46:24.799112Z","shell.execute_reply":"2022-07-21T22:46:25.600188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This is examples of and definition of positive correlation, value of one variable increase as the other increase, the same logic applies to negative correlation.  \nWe will further observed which feature has high positive correlation and choose them as our features to train model","metadata":{}},{"cell_type":"markdown","source":"## Observe missing values","metadata":{}},{"cell_type":"code","source":"train.shape","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:46:25.602784Z","iopub.execute_input":"2022-07-21T22:46:25.603092Z","iopub.status.idle":"2022-07-21T22:46:25.609997Z","shell.execute_reply.started":"2022-07-21T22:46:25.603062Z","shell.execute_reply":"2022-07-21T22:46:25.608720Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"total = train.isna().sum().sort_values(ascending=False)\npercent = total/len(train)\n\nmissing = pd.concat([total, percent], axis=1)\nmissing.columns = ['total', 'percentage']\n\n# let's also see corresponding correlation values\ncorr_tmp = corrmat.SalePrice\ncorr_tmp.name = 'corrval'\nmissingcorr = missing.merge(corr_tmp, how='outer', left_index=True,right_index=True).sort_values(by='percentage', ascending=False)\nmissingcorr.head(20)","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:46:25.611462Z","iopub.execute_input":"2022-07-21T22:46:25.611865Z","iopub.status.idle":"2022-07-21T22:46:25.647878Z","shell.execute_reply.started":"2022-07-21T22:46:25.611818Z","shell.execute_reply":"2022-07-21T22:46:25.646759Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"What can we can we learn here?:\n- Notice that not all variables has correlation values, this indicates that corresponding features are categorical (non-numeric features)\n- The first four feature contains lots of missing data > 80%, we can delete those feature entirely and pretend they don't exist\n- Same thing goes for `LotFrontage` and `FireplaceQu`, also notice that it has low correlation, so dropping it entirely won'e be a problem\n- Interesting percentage values for features `GarageXXX`, notice they have same number of missing values the belong to same instances, for now let's just drop them, plus the correlation values is not that high\n- The same logic applies to `BsmtXX` and `MasVnrXX`, although `MasVnrXX` has relatively low missing values, but for now let's just drop its column\n- `Electrical` only has one missing values, in this case we can delete the row","metadata":{}},{"cell_type":"code","source":"# Dropping columns\ndel_cols = missingcorr[missingcorr['percentage'] > missingcorr.loc['Electrical', 'percentage']].index\ndel_cols\nprint('Initial data shape:', train.shape)\n\n# The train_nona refers to train data without NaN values\ntrain_nona = train.drop(columns=del_cols)\nprint('After dropping columns:', train_nona.shape)\n\ntrain_nona = train_nona.dropna(axis=0, how='any')\nprint('After dropping instance of `Electrical`: ', train_nona.shape)\n\nprint('Total missing values in data after cleaning: ', train_nona.isna().sum().sum())","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:48:44.552539Z","iopub.execute_input":"2022-07-21T22:48:44.552950Z","iopub.status.idle":"2022-07-21T22:48:44.577581Z","shell.execute_reply.started":"2022-07-21T22:48:44.552915Z","shell.execute_reply":"2022-07-21T22:48:44.576568Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Feature Engineering\nThis section will cover the following:\n- pick which features to build our model, based on highest correlation values\n- check normality of features, normalized them if needed","metadata":{}},{"cell_type":"markdown","source":"## Pick important features\nLet's check the correlation matrix once again and sort them descendingly","metadata":{}},{"cell_type":"code","source":"highcorr = corrmat.SalePrice.sort_values(ascending=False)\n\n# Pick correlation values that are > 0.5\nhighcorr = highcorr[highcorr > 0.5]\nhighcorr","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:49:38.370341Z","iopub.execute_input":"2022-07-21T22:49:38.371067Z","iopub.status.idle":"2022-07-21T22:49:38.381747Z","shell.execute_reply.started":"2022-07-21T22:49:38.371016Z","shell.execute_reply":"2022-07-21T22:49:38.380781Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we got some meaningful features that probably best to train our model, but do we need them all?  \nLet's take a look at correlation heatmap with these features one more time","metadata":{}},{"cell_type":"code","source":"# Correlation matrix with above features\nhighcorrmat = corrmat.loc[highcorr.index, highcorr.index]\n\nhighcorrmat","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:49:57.848988Z","iopub.execute_input":"2022-07-21T22:49:57.849397Z","iopub.status.idle":"2022-07-21T22:49:57.988637Z","shell.execute_reply.started":"2022-07-21T22:49:57.849363Z","shell.execute_reply":"2022-07-21T22:49:57.987425Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Correlation heatmap\nplt.figure(figsize=(8, 5))\nsns.heatmap(highcorrmat, annot=True, fmt='.2f')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:50:14.479822Z","iopub.execute_input":"2022-07-21T22:50:14.480407Z","iopub.status.idle":"2022-07-21T22:50:15.463853Z","shell.execute_reply.started":"2022-07-21T22:50:14.480370Z","shell.execute_reply":"2022-07-21T22:50:15.463037Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now, let me introduce **multicollinearity**, this term means that there are high correlation between two independent variables/predictors/features, this can be cause a problem when we train our model. The example from above would be `GarageCars` and `GarageArea`, those two variables are highly correlated, this means we can drop one variable and keep the other. This makes the train data less redundant and also reduce its dimension, such that making our model less complex.\n\nAnother example of multicollinearity from above heatmap:\n- `TotalBsmtSF` and `1stFlrSF`\n- `GrLivArea` and `TotRmsAbvGrd`\n\nIntuitively, we understand why they have high correlations, because essentially both variables indicates the same thing, for example, the number of cars you can fit in your garage (`GarageCars`) implicitly dictates the area of the garage itself (`GarageArea`). The same goes for other pairs\n\nSo how to decide which feature to drop among those pair? For now, let's pick feature with higher correlation to our target variable, we pick:\n- `GarageCars`\n- `TotalBsmtSF`\n- `GrLivArea`\n\nAnd we drop: `GarageArea`, `1stFlrSF`, `TotRmsAbvGrd`  \nFor now, let's make it simplere and remove those in lowest three as well: `YearBuilt`, `YearRemodAdd`\n\nLet's remove those features","metadata":{}},{"cell_type":"code","source":"# Create new dataframe containing selected features \ndrop_cols = ['GarageArea', '1stFlrSF', 'TotRmsAbvGrd', 'YearBuilt', 'YearRemodAdd']\n\n# Select highest correlated features\nsel_train = train_nona[highcorrmat.index]\n\n# Remove feature with multicollinearity\nsel_train = sel_train.drop(columns=drop_cols)\n\nsel_train.shape","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:52:36.714041Z","iopub.execute_input":"2022-07-21T22:52:36.714729Z","iopub.status.idle":"2022-07-21T22:52:36.725818Z","shell.execute_reply.started":"2022-07-21T22:52:36.714672Z","shell.execute_reply":"2022-07-21T22:52:36.724653Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Pairplot & Outliers\nWe can see the plot between features altogether using pair plot, it helps us to see if patterns or outliers that may exist","metadata":{}},{"cell_type":"code","source":"sns.pairplot(sel_train)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:52:47.380233Z","iopub.execute_input":"2022-07-21T22:52:47.380621Z","iopub.status.idle":"2022-07-21T22:52:55.861691Z","shell.execute_reply.started":"2022-07-21T22:52:47.380587Z","shell.execute_reply":"2022-07-21T22:52:55.860914Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Few observations regarding the plot against `SalePrice`:\n- There are two highest values in `GrLivArea` that doesn't follow the trends\n- There are one instances with highest values in `TotalBsmtSF` that doesn't follow the trends\n\nLet's delete those values\n\nThere some three instance in `GarageCars` when it equals to 4, that doesn't follow the trends, but we will ignore it now","metadata":{}},{"cell_type":"code","source":"# Check index of those instance \nprint(sel_train['GrLivArea'].sort_values()[-2:].index)\nprint(sel_train['TotalBsmtSF'].sort_values()[-1:].index)","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:53:36.375869Z","iopub.execute_input":"2022-07-21T22:53:36.376418Z","iopub.status.idle":"2022-07-21T22:53:36.384744Z","shell.execute_reply.started":"2022-07-21T22:53:36.376383Z","shell.execute_reply":"2022-07-21T22:53:36.383594Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# We can just delete the instance from GrLivArea\nsel_train = sel_train.drop(index=sel_train['GrLivArea'].sort_values()[-2:].index)\n\n# Check to see if thsoe outliers have been removed\nfig, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(12, 5))\nfig.suptitle('After Removing Outliers')\nsns.scatterplot(x='GrLivArea', y='SalePrice', data=sel_train, ax=ax1)\nsns.scatterplot(x='TotalBsmtSF', y='SalePrice', data=sel_train, ax=ax2)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:53:37.113949Z","iopub.execute_input":"2022-07-21T22:53:37.114485Z","iopub.status.idle":"2022-07-21T22:53:37.594268Z","shell.execute_reply.started":"2022-07-21T22:53:37.114451Z","shell.execute_reply":"2022-07-21T22:53:37.593205Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Normality\nI consider this section to be quite advance, since we actually can proceed building model using data from before.\nBut since I've just learned about it, I thought I might as well put it here.\n\nThis section will check whether or not each features with numeric values follows normal distribution.  \nThis is done because data with normal distribution is favorable in machine learning settings.   \nHere we will apply log + 1 transformation to convert non-normal distribution to normal distribution.\n\nWe will also do this for columns with continuos values, i.e. `SalePrice`, `GrLivArea` and `TotalBsmtSF`","metadata":{}},{"cell_type":"markdown","source":"### `SalePrice`","metadata":{}},{"cell_type":"code","source":"from scipy.stats import norm","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:54:20.102141Z","iopub.execute_input":"2022-07-21T22:54:20.102551Z","iopub.status.idle":"2022-07-21T22:54:20.106942Z","shell.execute_reply.started":"2022-07-21T22:54:20.102516Z","shell.execute_reply":"2022-07-21T22:54:20.105568Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 4))\nsns.distplot(sel_train.SalePrice, kde=True, fit=norm, ax=ax1)\n_ = stats.probplot(sel_train.SalePrice, plot = ax2)\nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:54:23.231811Z","iopub.execute_input":"2022-07-21T22:54:23.232498Z","iopub.status.idle":"2022-07-21T22:54:24.027268Z","shell.execute_reply.started":"2022-07-21T22:54:23.232463Z","shell.execute_reply":"2022-07-21T22:54:24.026164Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see `SalePrice` is not normal, the distribution has positive skewness, and qq plot shows it doesn't follows diagonal line. Let's apply log transformation!","metadata":{}},{"cell_type":"code","source":"# Applying log transformation\nsel_train['SalePrice'] = np.log1p(sel_train.SalePrice)","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:55:06.143684Z","iopub.execute_input":"2022-07-21T22:55:06.144113Z","iopub.status.idle":"2022-07-21T22:55:06.150834Z","shell.execute_reply.started":"2022-07-21T22:55:06.144076Z","shell.execute_reply":"2022-07-21T22:55:06.149597Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# See effect of log transformation \nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 4))\nsns.distplot(sel_train.SalePrice, kde=True, fit=norm, ax=ax1)\n_ = stats.probplot(sel_train.SalePrice, plot = ax2)\nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:55:07.770286Z","iopub.execute_input":"2022-07-21T22:55:07.770875Z","iopub.status.idle":"2022-07-21T22:55:08.351656Z","shell.execute_reply.started":"2022-07-21T22:55:07.770810Z","shell.execute_reply":"2022-07-21T22:55:08.350499Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### `GrLivArea`","metadata":{}},{"cell_type":"code","source":"# Plot GrLiveArea\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 4))\nsns.distplot(sel_train.GrLivArea, kde=True, fit=norm, ax=ax1)\n_ = stats.probplot(sel_train.GrLivArea, plot = ax2)\nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:55:15.996638Z","iopub.execute_input":"2022-07-21T22:55:15.997252Z","iopub.status.idle":"2022-07-21T22:55:16.570296Z","shell.execute_reply.started":"2022-07-21T22:55:15.997216Z","shell.execute_reply":"2022-07-21T22:55:16.569263Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Same phenomenon as before, `GrLivArea` also experience positive skewness","metadata":{}},{"cell_type":"code","source":"# Apply log transformation\nsel_train.GrLivArea = np.log(sel_train.GrLivArea)","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:55:20.859553Z","iopub.execute_input":"2022-07-21T22:55:20.859929Z","iopub.status.idle":"2022-07-21T22:55:20.866051Z","shell.execute_reply.started":"2022-07-21T22:55:20.859895Z","shell.execute_reply":"2022-07-21T22:55:20.864814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# See effect of log transformation \nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 4))\nsns.distplot(sel_train.GrLivArea, kde=True, fit=norm, ax=ax1)\n_ = stats.probplot(sel_train.GrLivArea, plot = ax2)\nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:55:24.409279Z","iopub.execute_input":"2022-07-21T22:55:24.409698Z","iopub.status.idle":"2022-07-21T22:55:25.002210Z","shell.execute_reply.started":"2022-07-21T22:55:24.409663Z","shell.execute_reply":"2022-07-21T22:55:25.001387Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### `TotalBsmtSF`","metadata":{}},{"cell_type":"code","source":"# Plot TotalBsmtSF\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 4))\nsns.distplot(sel_train.TotalBsmtSF, kde=True, fit=norm, ax=ax1)\n_ = stats.probplot(sel_train.TotalBsmtSF, plot = ax2)\nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:55:26.986791Z","iopub.execute_input":"2022-07-21T22:55:26.987531Z","iopub.status.idle":"2022-07-21T22:55:27.626701Z","shell.execute_reply.started":"2022-07-21T22:55:26.987483Z","shell.execute_reply":"2022-07-21T22:55:27.625550Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Notice that couple of values are zeros, these values can be transformed using log.   \nTo solve this, we can only apply log transformation for the non-zero values.","metadata":{}},{"cell_type":"code","source":"# transform data\nsel_train['TotalBsmtSF'] = np.log1p(sel_train.TotalBsmtSF)","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:55:29.312249Z","iopub.execute_input":"2022-07-21T22:55:29.312652Z","iopub.status.idle":"2022-07-21T22:55:29.318063Z","shell.execute_reply.started":"2022-07-21T22:55:29.312618Z","shell.execute_reply":"2022-07-21T22:55:29.317256Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# See effect of log transformation\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 4))\nsns.distplot(sel_train.TotalBsmtSF[sel_train.TotalBsmtSF>0], kde=True, fit=norm, ax=ax1)\n_ = stats.probplot(sel_train.TotalBsmtSF[sel_train.TotalBsmtSF>0], plot = ax2)\nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:55:30.701305Z","iopub.execute_input":"2022-07-21T22:55:30.701954Z","iopub.status.idle":"2022-07-21T22:55:31.306830Z","shell.execute_reply.started":"2022-07-21T22:55:30.701898Z","shell.execute_reply":"2022-07-21T22:55:31.305717Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Split the target variable from the predictors\ny = sel_train.SalePrice\nX = sel_train.drop(columns='SalePrice')","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:55:32.985622Z","iopub.execute_input":"2022-07-21T22:55:32.985994Z","iopub.status.idle":"2022-07-21T22:55:32.992543Z","shell.execute_reply.started":"2022-07-21T22:55:32.985961Z","shell.execute_reply":"2022-07-21T22:55:32.991397Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Training","metadata":{}},{"cell_type":"code","source":"# Models\nfrom sklearn.model_selection import train_test_split, cross_val_score\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.linear_model import Lasso, LassoCV, Ridge, RidgeCV\nfrom sklearn.linear_model import ElasticNet, ElasticNetCV\nfrom sklearn.metrics import mean_squared_error","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:55:34.677831Z","iopub.execute_input":"2022-07-21T22:55:34.678245Z","iopub.status.idle":"2022-07-21T22:55:34.954048Z","shell.execute_reply.started":"2022-07-21T22:55:34.678207Z","shell.execute_reply":"2022-07-21T22:55:34.952833Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:55:36.886945Z","iopub.execute_input":"2022-07-21T22:55:36.887692Z","iopub.status.idle":"2022-07-21T22:55:36.894833Z","shell.execute_reply.started":"2022-07-21T22:55:36.887657Z","shell.execute_reply":"2022-07-21T22:55:36.893924Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Ordinary Least Square (OLS), Regression\nFirst is let's use ordinary least square (OLS) model, this is basic linear regression without regularization.   \nThe definition of regularization will not be discussed extensively in this course.","metadata":{}},{"cell_type":"code","source":"# Function to compute RMSE (Root Mean Squared Error), using 5-fold CV\ndef rmse(model, X, y, cv):\n    rmse = np.sqrt(-cross_val_score(model, X, y, scoring='neg_mean_squared_error', cv=cv))\n    return rmse.mean()","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:55:45.602932Z","iopub.execute_input":"2022-07-21T22:55:45.603469Z","iopub.status.idle":"2022-07-21T22:55:45.608668Z","shell.execute_reply.started":"2022-07-21T22:55:45.603434Z","shell.execute_reply":"2022-07-21T22:55:45.607546Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"lin = LinearRegression()","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:56:42.374418Z","iopub.execute_input":"2022-07-21T22:56:42.374977Z","iopub.status.idle":"2022-07-21T22:56:42.379845Z","shell.execute_reply.started":"2022-07-21T22:56:42.374943Z","shell.execute_reply":"2022-07-21T22:56:42.378841Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# For now we will use the root mean square error (rmse) metric to evaluate each model\nrmse_sc = rmse(lin, X, y, 5)\nrmse_sc","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:56:43.143035Z","iopub.execute_input":"2022-07-21T22:56:43.143715Z","iopub.status.idle":"2022-07-21T22:56:43.185331Z","shell.execute_reply.started":"2022-07-21T22:56:43.143652Z","shell.execute_reply":"2022-07-21T22:56:43.183975Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create List to append dictionary of scores and model\nall_scores = []","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:56:43.517720Z","iopub.execute_input":"2022-07-21T22:56:43.518110Z","iopub.status.idle":"2022-07-21T22:56:43.523129Z","shell.execute_reply.started":"2022-07-21T22:56:43.518074Z","shell.execute_reply":"2022-07-21T22:56:43.522158Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_scores.append(dict(model='OLD', score=rmse_sc))","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:56:43.859079Z","iopub.execute_input":"2022-07-21T22:56:43.859508Z","iopub.status.idle":"2022-07-21T22:56:43.863940Z","shell.execute_reply.started":"2022-07-21T22:56:43.859473Z","shell.execute_reply":"2022-07-21T22:56:43.862952Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Ridge Regression\nIf you have studied about L1 and L2 Regularization, it's good thing to know them by different name: \n- L1 is also called Lasso\n- L2 is also called Ridge\n\nSome source I recommend to read on:\n- https://towardsdatascience.com/l1-and-l2-regularization-methods-ce25e7fc831c\n- https://stats.stackexchange.com/questions/866/when-should-i-use-lasso-vs-ridge\n- https://stats.stackexchange.com/questions/200416/is-regression-with-l1-regularization-the-same-as-lasso-and-with-l2-regularizati","metadata":{}},{"cell_type":"code","source":"# let's try a default value alpha =1 \nridge = Ridge(alpha=1)","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:56:44.394399Z","iopub.execute_input":"2022-07-21T22:56:44.395031Z","iopub.status.idle":"2022-07-21T22:56:44.399769Z","shell.execute_reply.started":"2022-07-21T22:56:44.394996Z","shell.execute_reply":"2022-07-21T22:56:44.398605Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rmse(ridge, X, y, cv=5)","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:56:44.921313Z","iopub.execute_input":"2022-07-21T22:56:44.921702Z","iopub.status.idle":"2022-07-21T22:56:44.964494Z","shell.execute_reply.started":"2022-07-21T22:56:44.921667Z","shell.execute_reply":"2022-07-21T22:56:44.963672Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"If you ever heard of term `lambda` as regularization term, in this scikit module it's defined by the variable `alpha`, itgoverns how much we want to regularize the model, it's a hyperparameter, meaning we can set this value according to our will that gives the best metric we concern about, in this case `rmse`, let's create numbers for alpha","metadata":{}},{"cell_type":"code","source":"# Selection for alphas\nalphas = np.logspace(5,-5,50)\nalphas","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:56:45.924295Z","iopub.execute_input":"2022-07-21T22:56:45.924674Z","iopub.status.idle":"2022-07-21T22:56:45.932673Z","shell.execute_reply.started":"2022-07-21T22:56:45.924641Z","shell.execute_reply":"2022-07-21T22:56:45.931318Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Ridge CV is a ridge regression with built-in CV implementation\nridgecv = RidgeCV(alphas=alphas, scoring='neg_mean_squared_error')\nridgecv.fit(X, y)\n\n# RidgeCV gives us the model trained with best alpha values\nridgecv.alpha_","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:56:46.911934Z","iopub.execute_input":"2022-07-21T22:56:46.912296Z","iopub.status.idle":"2022-07-21T22:56:46.949026Z","shell.execute_reply.started":"2022-07-21T22:56:46.912264Z","shell.execute_reply":"2022-07-21T22:56:46.947999Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mod = Ridge(alpha=ridgecv.alpha_)\nmod.fit(X, y)","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:56:47.754534Z","iopub.execute_input":"2022-07-21T22:56:47.755211Z","iopub.status.idle":"2022-07-21T22:56:47.767109Z","shell.execute_reply.started":"2022-07-21T22:56:47.755152Z","shell.execute_reply":"2022-07-21T22:56:47.766132Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rmse(mod, X, y, cv=5)","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:56:48.402660Z","iopub.execute_input":"2022-07-21T22:56:48.403050Z","iopub.status.idle":"2022-07-21T22:56:48.449480Z","shell.execute_reply.started":"2022-07-21T22:56:48.403014Z","shell.execute_reply":"2022-07-21T22:56:48.448153Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Computing rmse\nrmse_sc = rmse(ridgecv, X, y, cv=5)\nrmse_sc","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:56:48.829465Z","iopub.execute_input":"2022-07-21T22:56:48.830126Z","iopub.status.idle":"2022-07-21T22:56:48.973947Z","shell.execute_reply.started":"2022-07-21T22:56:48.830090Z","shell.execute_reply":"2022-07-21T22:56:48.972847Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Oops, it gives the same values as before with OLS, don't worry it happens. Notice that also in this case, using default alpha value or new alpha value doesn't really affect rmse_sc.\n\nAlso one thing to note is that, ridgecv returns the ridge model trained with best alpha, ","metadata":{}},{"cell_type":"code","source":"all_scores.append(dict(model='Ridge', score=rmse_sc))","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:56:49.538852Z","iopub.execute_input":"2022-07-21T22:56:49.539480Z","iopub.status.idle":"2022-07-21T22:56:49.546712Z","shell.execute_reply.started":"2022-07-21T22:56:49.539445Z","shell.execute_reply":"2022-07-21T22:56:49.545542Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Lasso Regression\nThe same logic from before, applies to Lasso model","metadata":{}},{"cell_type":"code","source":"# Try lasso with 1 alpha values\nlasso = Lasso(alpha=1)\nrmse(lasso, X, y, cv=5)","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:56:50.333583Z","iopub.execute_input":"2022-07-21T22:56:50.334012Z","iopub.status.idle":"2022-07-21T22:56:50.380762Z","shell.execute_reply.started":"2022-07-21T22:56:50.333976Z","shell.execute_reply":"2022-07-21T22:56:50.379583Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Find new alphas\nlassocv = LassoCV(alphas=alphas)\nlassocv.fit(X, y)\nlassocv.alpha_","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:56:50.736387Z","iopub.execute_input":"2022-07-21T22:56:50.736753Z","iopub.status.idle":"2022-07-21T22:56:50.793603Z","shell.execute_reply.started":"2022-07-21T22:56:50.736720Z","shell.execute_reply":"2022-07-21T22:56:50.792423Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Compute rmse and add to list\nrmse_sc = rmse(lassocv, X, y, cv=5)\nrmse_sc","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:56:51.523897Z","iopub.execute_input":"2022-07-21T22:56:51.524422Z","iopub.status.idle":"2022-07-21T22:56:51.768612Z","shell.execute_reply.started":"2022-07-21T22:56:51.524388Z","shell.execute_reply":"2022-07-21T22:56:51.767588Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see that rmse with new alpha gives better (lower) rmse than default values.  \nBut still the rmse with new alpha is similar to previous two.","metadata":{}},{"cell_type":"code","source":"all_scores.append(dict(model='Lasso', score=rmse_sc))","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:56:52.603542Z","iopub.execute_input":"2022-07-21T22:56:52.604120Z","iopub.status.idle":"2022-07-21T22:56:52.609760Z","shell.execute_reply.started":"2022-07-21T22:56:52.604067Z","shell.execute_reply":"2022-07-21T22:56:52.608588Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Elastic Net\nElastic net is a combination of both L1 and L2 regularization.  \nWe will skip using the default alpha values, and immediately jump to using ElasticNetCV","metadata":{}},{"cell_type":"code","source":"elasticcv = ElasticNetCV(alphas=alphas)\nelasticcv.fit(X, y)","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:56:54.844819Z","iopub.execute_input":"2022-07-21T22:56:54.845213Z","iopub.status.idle":"2022-07-21T22:56:54.908339Z","shell.execute_reply.started":"2022-07-21T22:56:54.845150Z","shell.execute_reply":"2022-07-21T22:56:54.907056Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Compute rmse and add to list\nrmse_sc = rmse(elasticcv, X, y, cv=5)\nrmse_sc","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:56:55.187008Z","iopub.execute_input":"2022-07-21T22:56:55.187561Z","iopub.status.idle":"2022-07-21T22:56:55.435602Z","shell.execute_reply.started":"2022-07-21T22:56:55.187519Z","shell.execute_reply":"2022-07-21T22:56:55.434269Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_scores.append(dict(model='ElasticNet', score=rmse_sc))","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:56:55.694060Z","iopub.execute_input":"2022-07-21T22:56:55.694699Z","iopub.status.idle":"2022-07-21T22:56:55.701371Z","shell.execute_reply.started":"2022-07-21T22:56:55.694649Z","shell.execute_reply":"2022-07-21T22:56:55.699838Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Results summary ","metadata":{}},{"cell_type":"code","source":"pd.DataFrame(all_scores).sort_values(by='score')","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:56:56.866049Z","iopub.execute_input":"2022-07-21T22:56:56.866595Z","iopub.status.idle":"2022-07-21T22:56:56.880112Z","shell.execute_reply.started":"2022-07-21T22:56:56.866562Z","shell.execute_reply":"2022-07-21T22:56:56.878970Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Current results suggest Ridge Regression gives the lowest RSME score, let's use it to predict test data","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Test","metadata":{}},{"cell_type":"markdown","source":"## Grab desired columns and impute missing values","metadata":{}},{"cell_type":"code","source":"# Desired columns\nwant_cols = X.columns\n\nsel_test = test[want_cols]","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:57:20.553293Z","iopub.execute_input":"2022-07-21T22:57:20.553682Z","iopub.status.idle":"2022-07-21T22:57:20.559870Z","shell.execute_reply.started":"2022-07-21T22:57:20.553649Z","shell.execute_reply":"2022-07-21T22:57:20.558947Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Imput missing values with an average, notice the column with missing values\nsel_test.isna().sum()","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:57:21.772617Z","iopub.execute_input":"2022-07-21T22:57:21.773142Z","iopub.status.idle":"2022-07-21T22:57:21.782000Z","shell.execute_reply.started":"2022-07-21T22:57:21.773109Z","shell.execute_reply":"2022-07-21T22:57:21.780973Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Recall that `GarageCars` only contains integer values, so let's obey that. ","metadata":{}},{"cell_type":"code","source":"sel_test.GarageCars = sel_test.GarageCars.fillna(round(sel_test.GarageCars.mean()))\nsel_test.TotalBsmtSF = sel_test.TotalBsmtSF.fillna(sel_test.TotalBsmtSF.mean())","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:57:40.409881Z","iopub.execute_input":"2022-07-21T22:57:40.410417Z","iopub.status.idle":"2022-07-21T22:57:40.417159Z","shell.execute_reply.started":"2022-07-21T22:57:40.410382Z","shell.execute_reply":"2022-07-21T22:57:40.416272Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# After imputing missing values\nsel_test.isna().sum().sum()","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:57:41.148003Z","iopub.execute_input":"2022-07-21T22:57:41.148532Z","iopub.status.idle":"2022-07-21T22:57:41.156006Z","shell.execute_reply.started":"2022-07-21T22:57:41.148498Z","shell.execute_reply":"2022-07-21T22:57:41.154956Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Predict using Ridge","metadata":{}},{"cell_type":"code","source":"# Predict\nypred = ridgecv.predict(sel_test)\n\n# Creating output csv\nresult = pd.DataFrame({sample.columns[0] : sample['Id'],\n                        sample.columns[1] : ypred\n})\n","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:57:43.026341Z","iopub.execute_input":"2022-07-21T22:57:43.026847Z","iopub.status.idle":"2022-07-21T22:57:43.034947Z","shell.execute_reply.started":"2022-07-21T22:57:43.026813Z","shell.execute_reply":"2022-07-21T22:57:43.033805Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"result.to_csv('./20210629-housing-ridge.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2022-07-21T22:57:47.011728Z","iopub.execute_input":"2022-07-21T22:57:47.012339Z","iopub.status.idle":"2022-07-21T22:57:47.023667Z","shell.execute_reply.started":"2022-07-21T22:57:47.012304Z","shell.execute_reply":"2022-07-21T22:57:47.022266Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}