{ "cells": [ { "cell_type": "markdown", "metadata": { "_cell_guid": "e6335977-1dcb-a5bc-4856-184dd7bce3f9" }, "source": [ "## Trying out a linear model: \n", "\n", "Author: Alexandru Papiu \n", "If you use parts of this notebook in your own scripts, please give some sort of credit (for example link back to this). Thanks!\n", "\n", "\n", "There have been a few [great](https://www.kaggle.com/comartel/house-prices-advanced-regression-techniques/house-price-xgboost-starter/run/348739) [scripts](https://www.kaggle.com/zoupet/house-prices-advanced-regression-techniques/xgboost-10-kfolds-with-scikit-learn/run/357561) on [xgboost](https://www.kaggle.com/tadepalli/house-prices-advanced-regression-techniques/xgboost-with-n-trees-autostop-0-12638/run/353049) already so I'd figured I'd try something simpler: a regularized linear regression model. Surprisingly it does really well with very little feature engineering. The key point is to to log_transform the numeric variables since most of them are skewed." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "0d706811-b70c-aeab-a78b-3c7abd9978d3" }, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "import seaborn as sns\n", "import matplotlib\n", "\n", "import matplotlib.pyplot as plt\n", "from scipy.stats import skew\n", "from scipy.stats.stats import pearsonr\n", "\n", "\n", "%config InlineBackend.figure_format = 'retina' #set 'png' here when working on notebook\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "603292c1-44b7-d72a-5468-e6782f311603" }, "outputs": [], "source": [ "train = pd.read_csv(\"../input/train.csv\")\n", "test = pd.read_csv(\"../input/test.csv\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "d646bb1b-56c4-9b45-d5d4-27095f61b1c0" }, "outputs": [], "source": [ "train.head()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "cb2d88d7-7f76-4b04-d28b-d2c315ae4346" }, "outputs": [], "source": [ "all_data = pd.concat((train.loc[:,'MSSubClass':'SaleCondition'],\n", " test.loc[:,'MSSubClass':'SaleCondition']))" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "29fa13df-61e8-b0c2-b3a7-ea92bffd4396" }, "source": [ "###Data preprocessing: \n", "We're not going to do anything fancy here: \n", " \n", "- First I'll transform the skewed numeric features by taking log(feature + 1) - this will make the features more normal \n", "- Create Dummy variables for the categorical features \n", "- Replace the numeric missing values (NaN's) with the mean of their respective columns" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "9b5a3e5b-f683-3fd2-7269-4068975bbe42" }, "outputs": [], "source": [ "matplotlib.rcParams['figure.figsize'] = (12.0, 6.0)\n", "prices = pd.DataFrame({\"price\":train[\"SalePrice\"], \"log(price + 1)\":np.log1p(train[\"SalePrice\"])})\n", "prices.hist()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "4ed54771-95c4-00e7-b2cd-569d17862878" }, "outputs": [], "source": [ "#log transform the target:\n", "train[\"SalePrice\"] = np.log1p(train[\"SalePrice\"])\n", "\n", "#log transform skewed numeric features:\n", "numeric_feats = all_data.dtypes[all_data.dtypes != \"object\"].index\n", "\n", "skewed_feats = train[numeric_feats].apply(lambda x: skew(x.dropna())) #compute skewness\n", "skewed_feats = skewed_feats[skewed_feats > 0.75]\n", "skewed_feats = skewed_feats.index\n", "\n", "all_data[skewed_feats] = np.log1p(all_data[skewed_feats])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "3854ab12-a4f3-4c88-fe6e-1fee08e18af2" }, "outputs": [], "source": [ "all_data = pd.get_dummies(all_data)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "5d417300-0deb-3353-cabf-95f75af62678" }, "outputs": [], "source": [ "#filling NA's with the mean of the column:\n", "all_data = all_data.fillna(all_data.mean())" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "fe687685-cdac-0a89-4d71-af2d11d87a81" }, "outputs": [], "source": [ "#creating matrices for sklearn:\n", "X_train = all_data[:train.shape[0]]\n", "X_test = all_data[train.shape[0]:]\n", "y = train.SalePrice" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "cc4e3014-23b7-2971-ddb0-f67b03f83558" }, "source": [ "###Models\n", "\n", "Now we are going to use regularized linear regression models from the scikit learn module. I'm going to try both l_1(Lasso) and l_2(Ridge) regularization. I'll also define a function that returns the cross-validation rmse error so we can evaluate our models and pick the best tuning par" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "82886739-eee6-5d7a-4be9-e1fe6ac059f1" }, "outputs": [], "source": [ "from sklearn.linear_model import Ridge, RidgeCV, ElasticNet, LassoCV, LassoLarsCV\n", "from sklearn.model_selection import cross_val_score\n", "\n", "def rmse_cv(model):\n", " rmse= np.sqrt(-cross_val_score(model, X_train, y, scoring=\"neg_mean_squared_error\", cv = 5))\n", " return(rmse)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "436ce6e8-917f-8c88-3d7b-245e82a1619f" }, "outputs": [], "source": [ "model_ridge = Ridge()" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "69ff958c-dbbb-4750-3fb0-d0ac17ff6363" }, "source": [ "The main tuning parameter for the Ridge model is alpha - a regularization parameter that measures how flexible our model is. The higher the regularization the less prone our model will be to overfit. However it will also lose flexibility and might not capture all of the signal in the data." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "f6b86166-f581-6e05-5274-d3d3516ebaf3" }, "outputs": [], "source": [ "alphas = [0.05, 0.1, 0.3, 1, 3, 5, 10, 15, 30, 50, 75]\n", "cv_ridge = [rmse_cv(Ridge(alpha = alpha)).mean() \n", " for alpha in alphas]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "f8cf53ba-8441-9233-b7f5-a851d270b770" }, "outputs": [], "source": [ "cv_ridge = pd.Series(cv_ridge, index = alphas)\n", "cv_ridge.plot(title = \"Validation - Just Do It\")\n", "plt.xlabel(\"alpha\")\n", "plt.ylabel(\"rmse\")" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "37486402-4a48-912f-84ee-a3611334b133" }, "source": [ "Note the U-ish shaped curve above. When alpha is too large the regularization is too strong and the model cannot capture all the complexities in the data. If however we let the model be too flexible (alpha small) the model begins to overfit. A value of alpha = 10 is about right based on the plot above." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "d42c18c9-ee70-929f-ce63-aac7f77796cc" }, "outputs": [], "source": [ "cv_ridge.min()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "52a42dcc-9636-3dd5-5cb8-88947b5c77d6" }, "outputs": [], "source": [ "model_ridge = Ridge(alpha = 5).fit(X_train, y)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "fa64c37d-9d5f-5cd0-dd72-99d87a466a3b" }, "outputs": [], "source": [ "#let's look at the residuals as well:\n", "matplotlib.rcParams['figure.figsize'] = (6.0, 6.0)\n", "\n", "preds_ridge = pd.DataFrame({\"preds Ridge\":model_ridge.predict(X_train), \"true\":y})\n", "preds_ridge[\"residuals\"] = preds_ridge[\"true\"] - preds_ridge[\"preds Ridge\"]\n", "preds_ridge.plot(x = \"preds Ridge\", y = \"residuals\",kind = \"scatter\")" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "863fb699-7bcd-3748-3dbb-1c9b18afee9b" }, "source": [ "So for the Ridge regression we get a rmsle of about 0.127\n", "\n", "Let' try out the Lasso model. We will do a slightly different approach here and use the built in Lasso CV to figure out the best alpha for us. For some reason the alphas in Lasso CV are really the inverse or the alphas in Ridge." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "8204520c-a595-2ad2-4685-0b84cc662b84" }, "outputs": [], "source": [ "model_lasso = LassoCV(alphas = [1, 0.1, 0.001, 0.0005]).fit(X_train, y)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "e78e6126-4de0-08ad-250b-46a3f0f48de0" }, "outputs": [], "source": [ "rmse_cv(model_lasso).mean()" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "abc5f43e-1c38-4c1e-cb70-a95c8d9be8de" }, "source": [ "Nice! The lasso performs even better so we'll just use this one to predict on the test set. Another neat thing about the Lasso is that it does feature selection for you - setting coefficients of features it deems unimportant to zero. Let's take a look at the coefficients:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "c7be87ca-412a-cb19-1524-cd94cf698d44" }, "outputs": [], "source": [ "coef = pd.Series(model_lasso.coef_, index = X_train.columns)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "14be641e-bbe0-824d-d90f-f47698c8b5c5" }, "outputs": [], "source": [ "print(\"Lasso picked \" + str(sum(coef != 0)) + \" variables and eliminated the other \" + str(sum(coef == 0)) + \" variables\")" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "ca153134-b109-1afc-e066-44273f65d44c" }, "source": [ "Good job Lasso. One thing to note here however is that the features selected are not necessarily the \"correct\" ones - especially since there are a lot of collinear features in this dataset. One idea to try here is run Lasso a few times on boostrapped samples and see how stable the feature selection is." ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "632e23a8-948f-c692-f5e9-1aa8a75f21d5" }, "source": [ "We can also take a look directly at what the most important coefficients are:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "3efc02df-c877-b1fe-1807-5dd93c896c63" }, "outputs": [], "source": [ "imp_coef = pd.concat([coef.sort_values().head(10),\n", " coef.sort_values().tail(10)])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "87317789-6e7e-d57f-0b54-d8ba0ee26abf" }, "outputs": [], "source": [ "matplotlib.rcParams['figure.figsize'] = (8.0, 10.0)\n", "imp_coef.plot(kind = \"barh\")\n", "plt.title(\"Coefficients in the Lasso Model\")" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "f6e6f820-1ec6-4a69-c309-9992b80652da" }, "source": [ "The most important positive feature is `GrLivArea` - the above ground area by area square feet. This definitely sense. Then a few other location and quality features contributed positively. Some of the negative features make less sense and would be worth looking into more - it seems like they might come from unbalanced categorical variables.\n", "\n", " Also note that unlike the feature importance you'd get from a random forest these are _actual_ coefficients in your model - so you can say precisely why the predicted price is what it is. The only issue here is that we log_transformed both the target and the numeric features so the actual magnitudes are a bit hard to interpret. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "cdeaa3d3-f9ad-2e06-1339-61b4425a43f8" }, "outputs": [], "source": [ "#let's look at the residuals as well:\n", "matplotlib.rcParams['figure.figsize'] = (6.0, 6.0)\n", "\n", "preds = pd.DataFrame({\"preds\":model_lasso.predict(X_train), \"true\":y})\n", "preds[\"residuals\"] = preds[\"true\"] - preds[\"preds\"]\n", "preds.plot(x = \"preds\", y = \"residuals\",kind = \"scatter\")" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "4780532e-2815-e355-9a96-fb1c598f6984" }, "source": [ "The residual plot looks pretty good.To wrap it up let's predict on the test set and submit on the leaderboard:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "7c39c06f-c929-422e-52bc-a1ab762c55ff" }, "outputs": [], "source": [ "model_elas = ElasticNet(alpha=0.001, l1_ratio=0.5, fit_intercept=True, normalize=False, precompute=False, max_iter=1000, copy_X=True, tol=0.0001, warm_start=False, positive=False, random_state=None, selection='cyclic').fit(X_train, y)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "30d40f46-f885-741e-4d48-f91a7eb5d1b2" }, "outputs": [], "source": [ "rmse_cv(model_elas).mean()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "baef7ed8-1270-70d0-d8c3-a7d3d0ca9060" }, "outputs": [], "source": [ "#let's look at the residuals as well:\n", "matplotlib.rcParams['figure.figsize'] = (6.0, 6.0)\n", "\n", "preds_elas = pd.DataFrame({\"preds Elas\":model_elas.predict(X_train), \"true\":y})\n", "preds_elas[\"residuals\"] = preds_elas[\"true\"] - preds_elas[\"preds Elas\"]\n", "preds_elas.plot(x = \"preds Elas\", y = \"residuals\",kind = \"scatter\")" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "f8da43e0-fd51-a4c9-d9b2-364d9911699a" }, "source": [ "### Adding an xgboost model:" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "ae9bcc1a-5106-0909-ecf7-0abc2d2ca386" }, "source": [ "Let's add an xgboost model to our linear model to see if we can improve our score:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "654e4fcf-a049-921a-4783-3c6d6dcca673" }, "outputs": [], "source": [ "import xgboost as xgb" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "be53a9f8-d88b-05fb-734d-3a1388d39864" }, "outputs": [], "source": [ "\n", "dtrain = xgb.DMatrix(X_train, label = y)\n", "dtest = xgb.DMatrix(X_test)\n", "\n", "params = {\"max_depth\":6, \"eta\":0.1}\n", "model = xgb.cv(params, dtrain, num_boost_round=500, early_stopping_rounds=100)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "c9d5bfe5-a0a8-0b10-d54a-1dfbb6123bd0" }, "outputs": [], "source": [ "model.loc[30:,[\"test-rmse-mean\", \"train-rmse-mean\"]].plot()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "00b8a271-0f93-c757-7e33-516c3a297628" }, "outputs": [], "source": [ "model_xgb = xgb.XGBRegressor(n_estimators=360, max_depth=6, learning_rate=0.1) #the params were tuned using xgb.cv\n", "model_xgb.fit(X_train, y)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "1c9c640b-9e6c-a350-0691-7f6a7dc19c41" }, "outputs": [], "source": [ "predictions = pd.DataFrame({\"xgb\":xgb_preds, \"lasso\":lasso_preds})\n", "predictions.plot(x = \"xgb\", y = \"lasso\", kind = \"scatter\")" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "74c9bdd2-afbb-5fdf-a776-a8eebfa30d12" }, "source": [ "Many times it makes sense to take a weighted average of uncorrelated results - this usually imporoves the score although in this case it doesn't help that much." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "2bba690f-ad55-1dd5-19fa-2dccd695780b" }, "outputs": [], "source": [ "#gb_preds = np.expm1(model_gb.predict(X_test))\n", "lasso_preds = np.expm1(model_lasso.predict(X_test))\n", "ridge_preds = np.expm1(model_ridge.predict(X_test))\n", "xgb_preds = np.expm1(model_xgb.predict(X_test))\n", "elas_preds = np.expm1(model_elas.predict(X_test))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "623ed0fe-0150-5226-db27-90a321061d52" }, "outputs": [], "source": [ "preds = 0.15*xgb_preds + 0.75*elas_preds + 0.10*elas_preds" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "569d7154-e3b5-84ab-1d28-57bdc02955d9" }, "outputs": [], "source": [ "solution = pd.DataFrame({\"id\":test.Id, \"SalePrice\":preds})\n", "solution.to_csv(\"ridge_sol.csv\", index = False)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "80ccfbd7-8c93-59fd-5adb-60d84100a746" }, "outputs": [], "source": [ "print(solution)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "9bf8f306-7d6a-30bd-9738-d23a4fd8b3d6" }, "outputs": [], "source": [ "" ] } ], "metadata": { "_change_revision": 0, "_is_fork": false, "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.0" } }, "nbformat": 4, "nbformat_minor": 0 }