notebooks/Linear_Regression.ipynb

{"cells":[{"cell_type":"markdown","source":["# Linear Regression"],"metadata":{"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"517adce9-82d7-43b2-851a-77c27b8ddc90"}}},{"cell_type":"markdown","source":["## Import useful Python packages"],"metadata":{"id":"YOrmY8FiOMwa","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"f378aa8c-ba33-4b12-ba87-02960091c207"}}},{"cell_type":"code","source":["import requests\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n%matplotlib inline\n\nimport pyspark\nfrom pyspark.sql import *\nfrom pyspark.sql.types import *\nfrom pyspark.sql.functions import *\nfrom pyspark import SparkContext, SparkConf"],"metadata":{"id":"Fh8zPg5APmYv","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"7f41c6cf-8a94-41aa-ad5b-7e3d781bfee7"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["## Check everything is ok"],"metadata":{"id":"_eRvyl2xwCV6","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"56878782-5b05-46fc-85c5-b981d46cd1ce"}}},{"cell_type":"code","source":["spark"],"metadata":{"id":"4fqJ5f0JE3BL","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"01a94540-34c4-45e8-8393-0be7feea668d"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["sc._conf.getAll()"],"metadata":{"id":"qhTN342EEOYZ","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"4dcd187c-d9c6-4248-ab03-c242dee2bb3f"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["# **The Prediction Task**\n\nIn this notebook, we will be using a dataset from [Kaggle](https://www.kaggle.com/burhanykiyakoglu/predicting-house-prices/data) containing house sale records from [King County, WA, USA](https://en.wikipedia.org/wiki/King_County,_Washington) to predict house prices.\n\nIntuitively, each sold house is described by a set of properties (i.e., _features_) and the price it has been eventually sold. Our goal is to use those features (or, at least, a subset of them) to accurately predict the price which an _unseen_ house should be sold.\n\nIn the following, we implement each step of the pipeline which we have discussed in our classes."],"metadata":{"id":"Mu77QbD2vC_o","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"0ff2a8f2-3090-48b8-882e-2940b2c9f85f"}}},{"cell_type":"markdown","source":["## **1. Data Acquisition**\n\nThis is the first step we need to accomplish before going any further. The dataset will be downloaded and loaded to DBFS, as usual."],"metadata":{"id":"g-vLczwIuUhG","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"d1492548-37e6-4d4c-9dd9-616c34686968"}}},{"cell_type":"markdown","source":["### Download the dataset to the local driver node's ```/tmp``` folder using ```wget```"],"metadata":{"id":"j5fI3wiGvKBl","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"9bb8cf40-3ff7-4fe9-b0a7-464737238f32"}}},{"cell_type":"code","source":["%sh wget -P /tmp https://github.com/gtolomei/big-data-computing/raw/master/datasets/king-county-house-sales.csv.bz2"],"metadata":{"id":"kxrLDE_4e7KH","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"4cb180f1-2487-4305-b343-35d47b5e07e7"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["%fs ls file:/tmp/"],"metadata":{"id":"quiuGbfyv8vT","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"42723f28-7880-4535-942b-b82262e512ca"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### Move the file from local driver node's file system to DBFS"],"metadata":{"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"ee5b2ecd-374f-4614-8e53-1d788e34db37"}}},{"cell_type":"code","source":["dbutils.fs.mv(\"file:/tmp/king-county-house-sales.csv.bz2\", \"dbfs:/bdc-2020-21/datasets/king-county-house-sales.csv.bz2\")"],"metadata":{"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"e6408890-ad7e-4f55-b4af-ae1747d19390"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["%fs ls /bdc-2020-21/datasets/"],"metadata":{"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"1e0c9331-4e88-42e1-9e96-f6647f11328c"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Read dataset file into a Spark Dataframe**"],"metadata":{"id":"DevlrMcPw1ZI","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"b20dc7da-9117-411e-8398-55cb42548406"}}},{"cell_type":"code","source":["house_df = spark.read.load(\"dbfs:/bdc-2020-21/datasets/king-county-house-sales.csv.bz2\", \n                         format=\"csv\", \n                         sep=\",\", \n                         inferSchema=\"true\", \n                         header=\"true\"\n                         )"],"metadata":{"id":"qKi5Hd60FFcX","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"552f4965-db9d-4382-b8b0-661488c86f90"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Check the shape of the loaded dataset, i.e., number of rows and columns**"],"metadata":{"id":"IPTd8ep9x74H","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"6334ec72-2922-4b52-b92e-3be9948ce596"}}},{"cell_type":"code","source":["print(\"The shape of the dataset is {:d} rows by {:d} columns\".format(house_df.count(), len(house_df.columns)))"],"metadata":{"id":"JyRyYqeXGA4l","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"c0a0a5c2-310d-424d-8443-f41fd49a82bc"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Print out the schema of the loaded dataset**"],"metadata":{"id":"-WC4RPQgyEsB","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"50463729-da22-4159-bd9a-b0404b45b69d"}}},{"cell_type":"code","source":["house_df.printSchema()"],"metadata":{"id":"Q3KhtSnvGIwG","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"5efddbc0-b446-4860-965f-606c9909e09e"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Cast `price` type to `double` to avoid compatibility issues**"],"metadata":{"id":"YzkeDrSUAuRo","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"b4055440-d82b-4df7-868a-0105a2e81565"}}},{"cell_type":"code","source":["house_df = house_df.withColumn(\"price\", house_df[\"price\"].cast(\"double\"))"],"metadata":{"id":"xpe-xgZuAFCD","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"a867d221-13c7-4787-a7e4-c2f450d08a3d"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Dataset Shape and Schema**\n\nThe dataset contains **21,613** historical records of house sales; each record, contains the following set of **21** columns:\n- `id`: unique ID for each house sold (_numerical_, _discrete_);\n- `date`: Date when the house was sold (_datetime_);\n- **`price`**: The price a home has been sold (_numerical_, _continuous_); **[This is the _target_ variable we want to predict]**\n- `bedrooms`: Number of bedrooms (_numerical_, _discrete_);\n- `bathrooms`: Number of bathrooms (_numerical_, _continuous_) [.5 accounts for a room with a toilet but no shower];\n- `sqft_living`: Square footage of the house's interior living space (_numerical_, _continuous_);\n- `sqft_lot`: Square footage of the house's land space (_numerical_, _continuous_);\n- `floors`: Number of floors (_numerical_, _discrete_);\n- `waterfront`:  A dummy (i.e., binary) variable for whether the house is overlooking the waterfront or not (_categorical_, _nominal_);\n- `view`: An index from 0 to 4 of how good the view of the property was (_categorical_, _ordinal_);\n- `condition`: An index from 1 to 5 on the condition of the house (_categorical_, _ordinal_);\n- `grade`: An index from 1 to 13, where 1-3 falls short of building construction and design, 7 has an average level of construction and design, and 11-13 have a high quality level of construction and design (_categorical_, _ordinal_);\n- `sqft_above`: The square footage of the interior housing space that is above ground level (_numerical_, _continuous_);\n- `sqft_basement`: The square footage of the interior housing space that is below ground level (_numerical_, _continuous_);\n- `yr_built`: The year the house was initially built (_numerical_, _discrete_);\n- `yr_renovated`: The year of the house's last renovation (_numerical_, _discrete_);\n- `zipcode`: The zipcode area where the house is located (_categorical_, _nominal_);\n- `lat`: Latitude (_numerical_, _continuous_);\n- `long`: Longitude (_numerical_, _continuous_);\n- `sqft_living15`: The square footage of interior housing living space for the nearest **15** neighbors (_numerical_, _continuous_);\n- `sqft_lot15`: The square footage of the land lots of the nearest **15** neighbors (_numerical_, _continuous_)."],"metadata":{"id":"dXlcRfzmNNCv","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"4432942b-2b75-4553-9d4c-73032ea55d9a"}}},{"cell_type":"code","source":["# Let's define some constants which we will use throughout this notebook\nINDEX_FEATURE = \"id\"\nNUMERICAL_FEATURES = [\"bedrooms\", \n                      \"bathrooms\",\n                      \"floors\",\n                      \"lat\",\n                      \"long\",\n                      \"sqft_above\", \n                      \"sqft_basement\",\n                      \"sqft_living\", \n                      \"sqft_living15\", \n                      \"sqft_lot\", \n                      \"sqft_lot15\",\n                      \"yr_built\",\n                      \"yr_renovated\"\n                      ]\nCATEGORICAL_FEATURES = [\"condition\", \"grade\", \"view\", \"waterfront\", \"zipcode\"]\nTARGET_VARIABLE = \"price\""],"metadata":{"id":"xsHSUHEzDS1-","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"659e5da7-8967-410a-a0a4-446fa37aade1"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["print(\"1 Index feature = `{:s}`\".format(INDEX_FEATURE))\nprint(\"{:d} Numerical features = [{:s}]\".format(len(NUMERICAL_FEATURES), \", \".join([\"`{:s}`\".format(nf) for nf in NUMERICAL_FEATURES])))\nprint(\"{:d} Categorical features = [{:s}]\".format(len(CATEGORICAL_FEATURES), \", \".join([\"`{:s}`\".format(nf) for nf in CATEGORICAL_FEATURES])))\nprint(\"1 Target variable = `{:s}`\".format(TARGET_VARIABLE))"],"metadata":{"id":"pME8ArdkEt_M","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"a6fea43b-5e99-4714-aaca-38becfeeaa7a"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Display the first 5 rows of the dataset**"],"metadata":{"id":"0Pjib1fiylb6","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"1a9c86e0-b9b3-4752-a53d-bfa2f6e63e9a"}}},{"cell_type":"code","source":["house_df.show(5)"],"metadata":{"id":"a4v-Z92rGXoe","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"3a9f7289-e4c7-48f8-acb2-a5453a174ba9"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Check for any missing values**"],"metadata":{"id":"1Df7cSAKmqQi","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"21b53ca0-926f-421a-9e1a-b7f865156da3"}}},{"cell_type":"code","source":["for c in house_df.columns:\n  print(\"N. of missing values of column `{:s}` = {:d}\".format(c, house_df.where(col(c).isNull()).count()))"],"metadata":{"id":"ubkWy6w6mtvG","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"aa5b2f21-ae75-4dd1-938d-c115c7adfd7c"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["## **2. Data Exploration**\n\nBefore starting with the application of any (linear regression) modeling technique, a good practice is to first take a look at a few statistics computed from the data. In addition, drawing specific plots may help us spot interesting facts (e.g., the presence of outliers), which in turn could indicate us what steps to do next (e.g., outliers removal/winsorization)."],"metadata":{"id":"DHTM4wS2ztOJ","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"e6229f14-6127-490f-8f37-fe8dc40f547d"}}},{"cell_type":"markdown","source":["### **Summary of Descriptive Statistics**"],"metadata":{"id":"9TKWRYgl3jPi","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"28d7fdbc-381d-428e-b6d1-638140bc58f6"}}},{"cell_type":"code","source":["house_df.describe().show()"],"metadata":{"id":"FvTv38i73pVh","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"9162ed63-57b9-4226-8a8c-d7535cdbce11"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["## Built-in <code>display</code> function\n\nNote that, somehow, it is also possible to plot and visualize data using PySpark's built-in <code>display</code> function. In this way, you won't need to convert your PySpark DataFrame into a Pandas one; however, accessing to <code>seaborn</code>/<code>matplotlib</code> libraries after converting your DataFrame to Pandas seems still the most flexible choice."],"metadata":{"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"a22a02f0-1690-4ea0-a7df-a2a3fb1bf8d6"}}},{"cell_type":"code","source":["display(house_df)"],"metadata":{"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"63bee0e1-22c7-4979-bab1-440a70b26f68"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# To access plotting libraries, we need to first transform our PySpark DataFrame into a Pandas DataFrame\nhouse_pdf = house_df.toPandas() "],"metadata":{"id":"CK1VjwrvN4BV","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"3161c4d7-3728-4f41-8e00-013d0f462967"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# Set some default plotting configuration using seaborn properties\nsns.set_style(\"darkgrid\")\nsns.set_context(\"notebook\", rc={\"lines.linewidth\": 2, \n                                \"xtick.labelsize\":14, \n                                \"ytick.labelsize\":14,\n                                \"axes.labelsize\": 18\n                                })"],"metadata":{"id":"okVR9377jPTf","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"b0db2b92-541e-4a09-889f-871546560f37"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Distribution Plots**\n\nThis set of plots will show the distribution of values of a subset of **19** columns of interest. The following feature will **not** be included:\n- `id`\n- `date`"],"metadata":{"id":"DlRGHd9WOr5i","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"9f5d8c1f-4b2c-45a4-b560-b9d9c35f1c24"}}},{"cell_type":"code","source":["# Select only the columns which will be included in the distribution plots\nPLOT_COLUMNS = sorted(NUMERICAL_FEATURES + CATEGORICAL_FEATURES + [TARGET_VARIABLE])\nprint(\"Plotting {:d} columns from the dataset:\\n[{:s}]\".format(len(PLOT_COLUMNS), \n                                                              \", \".join([\"`{:s}`\".format(pc) for pc in PLOT_COLUMNS])\n                                                              ))"],"metadata":{"id":"j76rXSw6PtOc","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"9641d206-a1d9-4c5c-bf3f-d2f1a35388b5"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# Plot the distribution of values of each column of interest\nn_rows = 5\nn_cols = 4\n\nfig, axes = plt.subplots(n_rows, n_cols, figsize=(20,20))\nfig.delaxes(axes[4][3]) # the last ax will be empty as there are 19 columns for 20 slots, so just delete it!\n\nfor i,f in enumerate(PLOT_COLUMNS):\n  _ = sns.distplot(house_pdf[f],\n                   kde_kws={\"color\": \"#ca0020\", \"lw\": 1}, \n                   hist_kws={\"histtype\": \"bar\", \"edgecolor\": \"k\", \"linewidth\": 1,\"alpha\": 0.8, \"color\": \"#92c5de\"},\n                   ax=axes[i//n_cols, i%n_cols]\n                   )\n\nfig.tight_layout(pad=1.5)"],"metadata":{"id":"8esOySR-Oh1Q","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"1e74bbb5-3f09-425e-9463-e08f74245b23"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Relationship between _numerical continuous_ features and the _target variable_ (`price`)**"],"metadata":{"id":"__pFmdPp2cBJ","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"024fb539-3ca8-446e-93f2-2509086b22cc"}}},{"cell_type":"code","source":["NUMERICAL_CONTINUOUS_FEATURES = sorted(list(set(NUMERICAL_FEATURES) - set([\"bedrooms\", \"floors\", \"yr_renovated\"])))\nprint(\"Regression plots of {:d} numerical continuous features:\\n[{:s}]\".format(len(NUMERICAL_CONTINUOUS_FEATURES),\n                                                                              \", \".join([\"`{:s}`\".format(ncf) for ncf in NUMERICAL_CONTINUOUS_FEATURES])\n                                                                              ))"],"metadata":{"id":"fp5l-4kUfZQg","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"a7f63668-800e-4e9c-aff5-2a9e7e463711"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# Plot the relationship between each continuous feature (i.e., independent variable) with the target (i.e., dependent) variable\nn_rows = 2\nn_cols = 5\n\nfig, axes = plt.subplots(n_rows, n_cols, figsize=(20,14))\n\nfor i,f in enumerate(NUMERICAL_CONTINUOUS_FEATURES):\n  _ = sns.regplot(data=house_pdf, \n                  x=f, \n                  y=\"price\",\n                  color=\"#00b159\",\n                  ax=axes[i//n_cols, i%n_cols])\n\nfig.tight_layout(pad=1.5)"],"metadata":{"id":"YtNI7UoE6dw7","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"97f42a66-4b7e-4ff9-8022-5fac562ec25e"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["#### **Observations**\n\nIt seems there is a moderate linear relationship between `sqft_living` and `price`. Still, the former feature _alone_ might not be enough to correctly explain observed values of `price`."],"metadata":{"id":"nGBJyDm2KPVA","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"070c9d11-2b47-46d8-bdb2-6a28c9a1a048"}}},{"cell_type":"markdown","source":["### **Relationship between _numerical discrete_/_categorical_ features and the _target variable_ (`price`)**"],"metadata":{"id":"YBlrbZpqFsYI","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"b09a0a9a-5eb1-4f3b-8525-595e2d0b0a03"}}},{"cell_type":"code","source":["DISCRETE_FEATURES = sorted(list(set(CATEGORICAL_FEATURES + [\"bedrooms\", \"bathrooms\", \"floors\"]) - set([\"zipcode\"])))\nprint(\"Box plots of {:d} discrete features:\\n[{:s}]\".format(len(DISCRETE_FEATURES),\n                                                                              \", \".join([\"`{:s}`\".format(df) for df in DISCRETE_FEATURES])\n                                                                              ))"],"metadata":{"id":"DlFiw8ywF5xV","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"238ec6df-7d7d-41e3-8f85-77900c4c5deb"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# Plot the relationship between each discrete feature (i.e., independent variable) with the target (i.e., dependent) variable\nn_rows = 3\nn_cols = 2\n\nfig, axes = plt.subplots(n_rows, n_cols, figsize=(16,16))\n\nfor i,f in enumerate(DISCRETE_FEATURES[1:]):\n  _ = sns.boxplot(data=house_pdf,\n                  x=f, \n                  y=\"price\", \n                  ax=axes[i//n_cols, i%n_cols])\n  \nfig.tight_layout(pad=1.5)\n\n# Plot `bathrooms` on a dedicated chart\nfig, ax = plt.subplots(1, 1, figsize=(16,8))\n_ = sns.boxplot(data=house_pdf, x=\"bathrooms\", y=\"price\")"],"metadata":{"id":"Ow45r2cG2r9e","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"634ac0cb-aab1-456f-a3dc-732f32e8bf9b"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Correlation between features**\n\nWe know that having too many features in a model is not always a good thing because it might cause _overfitting_ and therefore bad generalization performance to unseen examples. Thus, if a feature does not improve our model a lot, removing it may indeed be a better choice.\n\nThings get even worst if there exists very high **correlation** between a subset of features, as keeping all of them most of the time will again cause overfitting. For instance, if we figure out that `sqt_above` is highly correlated with `sqt_living` there is no point for us to keep them both, as one can be easily derived from the other and together they do not add any valuable predictive power to our model.\n\nThe presence of high correlation (whether it be positive or negative) can be checked by computing and visualizing the **correlation matrix**. However, this does not mean that we must _always_ remove one of the highly correlated features. In some cases, even if two features are highly correlated they still may capture different aspects of the domain objects."],"metadata":{"id":"c7LJAWBLQcf6","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"b3dec80b-f98b-4809-a777-f40afe8b77ed"}}},{"cell_type":"markdown","source":["### **Pearson Correlation Coefficient and Correlation Matrix**\n\nRemember that given 2 random variables _X_ and _Y_, their _Pearson Correlation Coefficient_ can be computed as follows:\n$$ \\rho(X, Y) = \\frac{E \\Big[(X-E[X]) \\times (Y - E[Y])\\Big]}{\\sigma_X \\sigma_Y} = \\frac{E \\Big[(X-\\mu_X) \\times (Y - \\mu_Y)\\Big]}{\\sigma_X \\sigma_Y} $$\n\n$$ =Cov(X_{\\text{std}}, Y_{\\text{std}}) = E\\Big[ (X_{\\text{std}} - E[X_{\\text{std}}])  \\times (Y_{\\text{std}} - E[Y_{\\text{std}}])\\Big] = $$\n\n$$ =E\\Big[ (X_{\\text{std}} - 0)  \\times (Y_{\\text{std}} - 0)\\Big] = $$\n\n$$ =E \\Bigg[ \\frac{(X - \\mu_X)}{\\sigma_X} \\times \\frac{(Y - \\mu_Y)}{\\sigma_Y}\\Bigg ] = $$\n\n$$ =\\frac{E\\Big[(X - \\mu_X) \\times (Y - \\mu_Y)\\Big]}{\\sigma_X \\sigma_Y} = \\frac{Cov(X,Y)}{\\sigma_X \\sigma_Y} $$\n\nEventually, the **correlation matrix** _M_ is the triangular square matrix such that:\n\n$$ M[i][j] = \\rho(X_i, X_j) $$"],"metadata":{"id":"uUOs38v3kCuY","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"b0f86636-98d9-4031-bb88-e9e9efdc633a"}}},{"cell_type":"code","source":["# Select the features we want to use to compute the correlation matrix (i.e., everything except `id` and `date`)\nfeatures = [TARGET_VARIABLE] + sorted(NUMERICAL_FEATURES + CATEGORICAL_FEATURES)\n\nmask = np.zeros_like(house_pdf[features].corr(), dtype=np.bool) \nmask[np.triu_indices_from(mask)] = True \n\nwith sns.axes_style(\"white\"): # Temporarily set the background to white\n  fig, ax = plt.subplots(figsize=(16, 12))\n  plt.title('Pearson Correlation Matrix', fontsize=24)\n\n  cmap = sns.diverging_palette(220, 10, as_cmap=True)\n\n  _ = sns.heatmap(house_pdf[features].corr(), \n              linewidths=0.25, \n              vmax=0.7, \n              square=True,\n              ax=ax, \n              cmap=cmap, \n              linecolor='w', \n              annot=True, \n              annot_kws={\"size\":8}, \n              mask=mask, \n              cbar_kws={\"shrink\": .9});"],"metadata":{"id":"KYfav2GJS0yv","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"b92ea235-2dac-46d2-8cb9-6c3e3745d6c4"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["#### **Observations**\n\nAs expected, `sqft_living` and `sqft_above` are high-positively correlated with each other ($\\rho = 0.88$). Similarly, also `bathrooms`/`bedrooms` and `sqft_living` are high-positively correlated with each other ($\\rho$=0.75 and 0.70, respectively). Again, this seems totally legitimate, as the number of bathrooms and bedrooms depends on the square footage of a house.\n\nLater on, we will see how to deal with this multicollinearity between features, and what will be the effect on the quality of the learned model when those features are retained or removed."],"metadata":{"id":"qTeXo6eFL3Em","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"fe16e99b-7d67-42b4-ac5d-e837259e4960"}}},{"cell_type":"markdown","source":["## ** The Learning Pipeline**"],"metadata":{"id":"66pFGdch-BU_","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"7ed3bf71-b094-4693-8f03-bb866566696e"}}},{"cell_type":"markdown","source":["### **Important Remark on Dataset Transformations**\n\nAny dataset transformation we want to perform must be **derived** from the _training set_ portion of the original dataset, and then applied to (_validation_ and) _test set_.\n\nJust to make an example, if we need to standardize our feature values we will also need to compute the **mean** and the **standard deviation** of each feature. Now, those statistics **cannot** be computed by looking at the whole dataset. In other words, we first need to split our dataset into at least two portions: _training_ and _test set_. Then, we will compute the mean and standard deviation of all the feature values from the observations contained in the training set, and standardize those values accordingly. Of course, we will need also to standardize features appearing in the test set, and in order to do that you will use the means and standard deviations **as previoulsly computed from the training set**.\n\nThis is just following the general principle: anything we learn must be learned from the model's training data."],"metadata":{"id":"Wv3q6FKlWgbu","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"9e8d9b27-1d2f-4964-b59a-393d1aa064c0"}}},{"cell_type":"markdown","source":["### **Dataset Splitting: Training vs. Test Set**\n\nBefore moving along with any preprocessing involving data transformations, we will split our dataset into **2** portions:\n- _training set_ (e.g., accounting for 90% of the total number of instances);\n- _test set_ (e.g., accounting for the remaining 10% of instances)"],"metadata":{"colab_type":"text","id":"xlUz0-TNb4Jo","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"9f86b3e6-c951-41ce-ae4a-75611b4f9370"}}},{"cell_type":"code","source":["RANDOM_SEED = 42"],"metadata":{"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"2d93f772-a336-4fb1-8897-b13fba10161c"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# Randomly split our original dataset `house_df` into 90÷10 for training and test, respectively\ntrain_df, test_df = house_df.randomSplit([0.9, 0.1], seed=RANDOM_SEED)"],"metadata":{"id":"ZSzZLA9QcA_P","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"9cc9f1ee-b754-4295-aafd-ce53a4c9096d"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Working on the Training Set only**\n\nFrom now on, we will be working on the training set portion only. The test set will come back into play when we evaluate our learned model."],"metadata":{"id":"Kh179KdCcS81","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"97cdd8e7-4bb6-40a6-a56d-220195d14f3d"}}},{"cell_type":"markdown","source":["## **The Quickest Linear Regression Model**\n\nInitially, we try to learn a simple (i.e., _univariate_) linear regression model. To this end, we select only one feature from the original set, which we believe it is mostly related to our target variable `price`.\nAs an example, suppose we choose `sqft_living`.\n\nOur linear regression model would thus look like the following:\n\n$$ \\underbrace{h_{\\theta}(x)}_{\\texttt{price}} = \\theta_0 + \\theta_1 \\cdot \\underbrace{x}_{\\texttt{sqft\\\\_living}} $$"],"metadata":{"id":"tdIyqhU89eju","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"f8cb9612-8723-4d5b-87f0-54808eebb201"}}},{"cell_type":"markdown","source":["### **Extract the single predictor variable (`sqft_living`) and the target variable (`price`)**"],"metadata":{"id":"-VmTwlYp-cUL","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"d151cfb3-a4b4-43c0-acb3-ce5059792200"}}},{"cell_type":"code","source":["# Select `sqft_living` feature only together with the target variable `price` from the original `train_df`\nfrom pyspark.ml.feature import VectorAssembler\n\nassembler = VectorAssembler(inputCols=[\"sqft_living\"], \n                            outputCol=\"features\")\n\nsimple_train_df = assembler.transform(train_df)"],"metadata":{"id":"IZc7OWNK-E9h","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"228c871f-6aaa-4de2-916c-f214bbd190b8"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# Select only the two columns which we are interested in, i.e., `features` and `price`\nsimple_train_df = simple_train_df.select([\"features\", \"price\"])\nsimple_train_df.show(5, truncate=False)"],"metadata":{"id":"5S0LFBk5-XMN","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"31b042bf-14ef-478e-a1b4-1a5842962551"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **The General Optimization Problem behind Linear Regression**\nGiven an _m_-by-_n+1_ **feature matrix** **_X_**, an _m_-by-1 **target vector** **_y_**, and an _n+1_-by-1 **parameter vector** __theta__ associated with a hypothesis function defined as:\n\n$$ h_{\\boldsymbol{\\theta}}(\\mathbf{x}): \\theta_0 + \\theta_1 x_1 + \\ldots + \\theta_n x_n = \\mathbf{x}^T\\cdot \\boldsymbol{\\theta} $$\n\nIn its simplest form (a.k.a. **Ordinary Least Squares** or **OLS**), the linear regression problem aims to find the **parameter vector** __theta^*__  which minimizes the **Mean Squared Error**.\n\n$$\n\\boldsymbol{\\theta}^* = \\text{argmin}_{\\boldsymbol{\\theta}\\in \\mathbb{R}^n} \\underbrace{\\frac{1}{m} ||\\mathbf{X}\\cdot \\boldsymbol{\\theta} -\\mathbf{y}||^2}_{\\text{Mean Squared Error}}\n$$\n\nOLS is a special case of a more general optimization framework, which is known as **Elastic Net** and is defined as follows:\n\n$$\n\\boldsymbol{\\theta}^* = \\text{argmin}_{\\boldsymbol{\\theta}\\in \\mathbb{R}^n} \\frac{1}{m} ||\\mathbf{X}\\cdot \\boldsymbol{\\theta} -\\mathbf{y}||^2 + \\lambda\\Big(\\alpha |\\boldsymbol{\\theta}| + (1-\\alpha)||\\boldsymbol{\\theta}||^2\\Big)\n$$\nwhere:\n\n$$ \\lambda \\geq 0 \\text{ is the {\\bf regularization parameter}} $$\n$$ \\alpha \\in [0,1] \\text{ is the {\\bf tradeoff parameter} for the regularization penalties} $$\n$$ |\\boldsymbol{\\theta}| \\text{ is the L1-norm of the {\\bf parameter vector}} $$\n$$ ||\\boldsymbol{\\theta}||^2 \\text{ is the squared L2-norm (i.e., squared Euclidean norm) of the {\\bf parameter vector}} $$\n\nDepending on the values that are plugged in the general optimization framework above, we have the following cases:\n\n$$ \\lambda = 0 \\Rightarrow \\text{ {\\bf OLS} (no regularization at all)} $$\n$$ \\lambda > 0,~\\alpha = 1 \\Rightarrow \\text{ {\\bf Least Absolute Shrinkage and Selection Operator} or {\\bf LASSO} (L1-regularization only)} $$\n$$ \\lambda > 0,~\\alpha = 0 \\Rightarrow \\text{ {\\bf Ridge} (L2-regularization only)} $$\n$$ \\text{Any other case where } \\lambda > 0,~0 < \\alpha < 1 \\text{ balances between L1- and L2-regularization} $$"],"metadata":{"id":"zPF-F3d0H0yh","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"ffcd78ec-57c1-47d3-9511-443c6817e983"}}},{"cell_type":"markdown","source":["### **Learn the model using the `LinearRegression` object**\n\nWe use the `LinearRegression` object provided within the `pyspark.ml.regression` package. Any further information can be found in the [API documentation](https://spark.apache.org/docs/latest/api/python/pyspark.ml.html#pyspark.ml.regression.LinearRegression).\n\nAmong the parameters we can specify, let's have a look at the most important ones:\n- `regParam` is the regularization parameter (or $\\lambda$);\n- `elasticNetParam` is the tradeoff parameter for regularization penalties (or $\\alpha$);\n  - `regParam = 0` and `elasticNetParam = 0` means there is no regularization (i.e., OLS);\n  - `regParam > 0` and `elasticNetParam = 0` means there is only L2-regularization (Ridge); \n  - `regParam > 0` and `elasticNetParam = 1` means there is only L1-regularization (LASSO);\n  - `regParam > 0` and `0 < elasticNetParam < 1` means there is both L1- and L2-regularization (Elastic Net);\n- `solver` can take one of the following values: `l-bfgs`, `normal` and `auto`. \n  - `l-bfgs` stands for **Limited-memory BFGS** which is a limited-memory quasi-Newton optimization _iterative_ method; \n  - `normal` refers to **Normal Equation** _analytical_ method, which computes the solution to the optimization problem _directly_ using the **pseudo-inverse** of the feature matrix and multiplying it by the target vector;\n  - `auto` leaves PySpark to select the solver algorithm automatically."],"metadata":{"id":"MJCWzadb-qL-","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"250c46dd-5edb-43a6-b3bc-372443d37c9d"}}},{"cell_type":"markdown","source":["### **Training a simple univariate OLS linear regressor (with no regularization)**"],"metadata":{"id":"_7YsIr07i5RK","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"c4af2594-b0da-43a2-81b0-dfdfdd606872"}}},{"cell_type":"code","source":["from pyspark.ml.regression import LinearRegression\n\n# First of all, let's just setup a very basic OLS linear regressor (i.e., univariate with no regularizatio term)\nlr = LinearRegression(featuresCol=\"features\", labelCol=\"price\")\n# Train the actual model on our training set `simple_train_df`\nlr_model = lr.fit(simple_train_df)"],"metadata":{"id":"k2YVAZcO9u1c","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"50c70415-b728-4314-ad3e-4d5a22f6f6f4"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["#### **Intercept ($\\theta_0$) and Coefficient ($\\theta_1$)**"],"metadata":{"id":"J_3Hblc4jCIP","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"1bd735eb-84df-4e2e-a613-a41247825fe1"}}},{"cell_type":"code","source":["print(\"Intercept: \" + str(lr_model.intercept))\nprint(\"Coefficient: \" + str(lr_model.coefficients))"],"metadata":{"id":"__i2fdx69w5d","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"806b1bc6-a52d-48a1-89ac-0f07c0dbde89"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["#### **Measuring performance on the Training Set: $\\text{RMSE}$ and $\\text{R}^2$ statistic**"],"metadata":{"id":"eWBZCzw6jjJm","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"cc80a0d3-f419-4e3c-81bd-f0bf0296814c"}}},{"cell_type":"code","source":["training_result = lr_model.summary\nprint(\"***** Training Set *****\")\nprint(\"RMSE: {:.3f}\".format(training_result.rootMeanSquaredError))\nprint(\"R2: {:.3f}\".format(training_result.r2))\nprint(\"Adjusted R2: {:.3f}\".format(training_result.r2adj))\nprint(\"***** Training Set *****\")"],"metadata":{"id":"JPRYCw3y90OI","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"056e80bd-bfc2-459b-811d-f251641309a8"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Testing the simple univariate OLS linear regressor learned above on the Test Set**"],"metadata":{"id":"a8acPJn5j9rc","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"e1e3ce81-7b35-4812-9b46-5e635dfc78d9"}}},{"cell_type":"code","source":["# Select `sqft_living` feature only together with the target variable `price` from the original `train_df`\nfrom pyspark.ml.feature import VectorAssembler\n\nassembler = VectorAssembler(inputCols=[\"sqft_living\"], \n                            outputCol=\"features\")\n\nsimple_test_df = assembler.transform(test_df)"],"metadata":{"id":"pcMHUfGb-NVe","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"d46eff1c-86a9-4fe2-b122-fd638f7076dc"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# `lr_model` is a Transformer which can be used to \"transform\" our test set\nlr_predictions = lr_model.transform(simple_test_df) "],"metadata":{"id":"WYLcBk8R-AI3","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"c9cdbb73-4429-461e-af79-6e5ad56fa8ab"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# `lr_predictions` is a dataframe containing (among other things) the predictions made by `lr_model` on the test set\nlr_predictions.select(\"features\", \"prediction\", \"price\").show(5)"],"metadata":{"id":"igZ7TcCPqqxM","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"bce01a5e-df41-40cd-8307-5b08269db0ba"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["#### **Measuring performance on the Test Set: $\\text{RMSE}$ and $\\text{R}^2$ statistic**"],"metadata":{"id":"Dy4MzTWckKxV","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"78a55519-bedb-421f-94b1-da6b81f199b1"}}},{"cell_type":"code","source":["test_result = lr_model.evaluate(simple_test_df)\nprint(\"***** Test Set *****\")\nprint(\"RMSE: {:.3f}\".format(test_result.rootMeanSquaredError))\nprint(\"R2: {:.3f}\".format(test_result.r2))\nprint(\"Adjusted R2: {:.3f}\".format(test_result.r2adj))\nprint(\"***** Test Set *****\")"],"metadata":{"id":"JlOofISU-dmw","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"919626a2-ea01-4cd9-92a4-35f09339e918"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["#### **Observations**\n\nIt turns out that the very simple univariate linear regressor does not work well. This may be due to the fact that a single explanatory variable (feature) alone like `sqft_living` could be not enough to nicely predict our target variable `price`.\n\nAlso, the fact that both $R^2$ and $R^2_{\\text{adj}}$ are slightly better on the test set than on the training set may be simply due to a \"lucky\" random split of the initial dataset."],"metadata":{"id":"LSh-GRysnVad","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"ddb95ab6-5cec-4f0d-bc82-b6c8e15e31ce"}}},{"cell_type":"markdown","source":["## **Summary of the main findings from Data Exploration step**\n\nLet's now go back to the result of our exploratory data analysis above, where we were able to observe the following:\n- Both continuous and categorical features;\n- Different scales of feature values;\n- Several outliers on many features (see the boxplots above)."],"metadata":{"id":"Uzic8uSpl8y3","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"c48f1f76-30eb-40a0-81e7-8030a5207680"}}},{"cell_type":"markdown","source":["### **Transform Categorical features into Numerical using One-Hot Encoding**\n\nLet's address the first issue indicated above.\n\nTo transform _categorical_ features into _numerical_ ones we proceed as follows.\nWe setup a pipeline which is composed of the following steps:\n- [`StringIndexer`](https://spark.apache.org/docs/latest/ml-features#stringindexer): encodes a string column of labels to a column of label indices. The indices are in `[0, numLabels)`, and 4 ordering options are supported (default `frequencyDesc`, which assigns the most frequent label the index `0`, and so on and so forth).\n- [`OneHotEncoderEstimator`](https://spark.apache.org/docs/latest/ml-features#onehotencoderestimator): maps a categorical feature, represented as a label index, to a binary vector with at most a single one-value indicating the presence of a specific feature value from among the set of all feature values. An important parameter is `handleInvalid`, which indicates how to deal with previously unseen labels. By default this raises an error but it can be set to as `keep` to assign previously unseen labels a fallback value.\n- [`VectorAssembler`](https://spark.apache.org/docs/latest/ml-features#vectorassembler): is a transformer that combines a given list of columns into a single vector column."],"metadata":{"id":"11HcGwvdg90A","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"c5ab233c-b4d3-408d-bbf9-314711f8d2af"}}},{"cell_type":"code","source":["# This function is responsible to implement the pipeline above for transforming categorical features into numerical ones\ndef to_numerical(df, numerical_features, categorical_features, target_variable):\n\n    \"\"\"\n    Args:\n        - df: the input dataframe\n        - numerical_features: the list of column names in `df` corresponding to numerical features\n        - categorical_features: the list of column names in `df` corresponding to categorical features\n        - target_variable: the column name in `df` corresponding to the target variable\n\n    Return:\n        - transformer: the pipeline of transformation fit to `df` (for future usage)\n        - df_transformed: the dataframe transformed according to the pipeline\n    \"\"\"\n    \n    from pyspark.ml import Pipeline\n    from pyspark.ml.feature import StringIndexer, OneHotEncoder, VectorAssembler\n\n    # 1. Create a list of indexers, i.e., one for each categorical feature\n    indexers = [StringIndexer(inputCol=c, outputCol=\"{0}_indexed\".format(c), handleInvalid=\"keep\") for c in categorical_features]\n\n    # 2. Create the one-hot encoder for the list of features just indexed (this encoder will keep any unseen label in the future)\n    encoder = OneHotEncoder(inputCols=[indexer.getOutputCol() for indexer in indexers], \n                                    outputCols=[\"{0}_encoded\".format(indexer.getOutputCol()) for indexer in indexers], \n                                    handleInvalid=\"keep\")\n    \n    # 3. Assemble all the features into a single vector\n    assembler = VectorAssembler(inputCols=encoder.getOutputCols() + numerical_features, outputCol=\"features\")\n\n    # 4. Setup the pipeline with the stages above\n    pipeline = Pipeline(stages=indexers + [encoder] + [assembler])\n\n    # 5. Transform the input dataframe accordingly\n    transformer = pipeline.fit(df)\n    df_transformed = transformer.transform(df)\n\n    # 6. Optionally, change the name of the target variable to `label` (if that is different), as PySpark implicitly assumes this is the name of the column to predict\n    if target_variable != \"label\":\n        df_transformed = df_transformed.withColumnRenamed(target_variable, \"label\")\n\n    # 7. Eventually, return both the transformed dataframe and the transformer object for future transformations\n    return transformer, df_transformed "],"metadata":{"id":"BEHqXC-yvTft","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"e4c06e64-6e84-46e0-8163-72f79d20f515"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# Transform the training set and get back both the transformer and the new dataset\noh_transformer, oh_train_df = to_numerical(train_df, NUMERICAL_FEATURES, CATEGORICAL_FEATURES, TARGET_VARIABLE)"],"metadata":{"id":"ueRiKbiTxJsr","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"d0b16822-b5df-4d2e-a4a4-d97545fc8c8b"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# Show the result of numerical transformation\noh_train_df.show(5)"],"metadata":{"id":"vTix1SBg71Xr","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"d54ca44a-9985-4d8f-9315-403e783219b6"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# Select `features` and `label` (i.e., formerly `price`) target variable only\ntrain = oh_train_df.select([\"features\", \"label\"])"],"metadata":{"id":"yYpzQRWqyxUx","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"ec65843b-825e-44a8-bb3d-1ce9fdfd51b4"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["train.show(5, truncate=False)"],"metadata":{"id":"K2pZg_Pey-yP","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"6882384e-1a9d-4700-bde1-3990c8d1e960"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Training an OLS Multivariate Linear Regressor (with no regularization) using all the features**"],"metadata":{"id":"teSVf-WFzvs-","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"6cb7bc78-819a-4b67-b31e-fbdbf2808de0"}}},{"cell_type":"code","source":["# `lr` is still the LinearRegressor object we set up above (the only difference here is that now `features` refers to the whole feature vector)\nlr = LinearRegression(featuresCol=\"features\", labelCol=\"label\") # Analogous to lr = LinearRegression()\n# Train the actual model on our training set `train`\nlr_model = lr.fit(train)"],"metadata":{"id":"Mmgq2uMK0AMF","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"8f4da453-2039-4af4-92ea-7103bd4edcae"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Intercept ($\\theta_0$) and Coefficients ($\\theta_1, \\ldots, \\theta_n$)**"],"metadata":{"id":"vf8I0bal0hsi","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"95cf2f47-807e-43b8-b592-671eead0965a"}}},{"cell_type":"code","source":["print(\"Intercept: \" + str(lr_model.intercept))\nprint(\"Coefficients: \" + str(lr_model.coefficients))"],"metadata":{"id":"NyWx968q0zdk","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"50b89191-4995-4daf-ad9f-cc3e7044ca08"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["#### **Measuring performance on the Training Set: $\\text{RMSE}$ and $\\text{R}^2$ statistic**"],"metadata":{"id":"_QEPpDoV1JnK","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"463b752d-c75a-4c57-b422-d235b06cc0e2"}}},{"cell_type":"code","source":["training_result = lr_model.summary\nprint(\"***** Training Set *****\")\nprint(\"RMSE: {:.3f}\".format(training_result.rootMeanSquaredError))\nprint(\"R2: {:.3f}\".format(training_result.r2))\nprint(\"Adjusted R2: {:.3f}\".format(training_result.r2adj))\nprint(\"***** Training Set *****\")"],"metadata":{"id":"OsID3E0r07i5","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"70b1cc42-9647-4fd6-a175-90944ae3ada7"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Use the One-Hot encoding pipeline to transform the Test Set**"],"metadata":{"id":"ena6EJ0Y5qhx","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"532a99c0-b650-4f9f-b017-c1dae1ada35b"}}},{"cell_type":"code","source":["# Here, we use the same transformer as the one returned by the `to_numerical` function above yet applied to the test set\noh_test_df = oh_transformer.transform(test_df)"],"metadata":{"id":"KPgeFLQC1S1t","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"565f836c-5f08-4982-a898-8d1ef6fcf7d5"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["oh_test_df.show(5)"],"metadata":{"id":"YoU_fGuA2ivF","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"8792e321-5ca3-44df-8197-dba1c769693c"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# Rename the target variable `price` to `label`\noh_test_df = oh_test_df.withColumnRenamed(TARGET_VARIABLE, \"label\")\noh_test_df.show(5)"],"metadata":{"id":"slDBka0G1fzo","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"05f9d3ed-9eab-469a-b461-e2e5b65b4669"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# Select `features` and `label` only\ntest = oh_test_df.select([\"features\", \"label\"])\ntest.show(5)"],"metadata":{"id":"oMNJ1NdtwXB8","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"019c9bfe-f1ef-4699-9a19-9b5b22670ac1"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Compute predictions on the Test Set according to the model learned on the Training Set**"],"metadata":{"id":"LokzdSNP6T48","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"3926f05b-82b5-499b-90a4-8fea0c78f097"}}},{"cell_type":"code","source":["# `lr_model` is a Transformer which can be used to \"transform\" our test set\nlr_predictions = lr_model.transform(test)"],"metadata":{"id":"Y-tiLfEA1tml","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"f475022f-97c7-42a6-ab9d-add1fc5ecab6"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# `lr_predictions` is a dataframe containing (among other things) the predictions made by `lr_model` on the test set\nlr_predictions.select(\"features\", \"prediction\", \"label\").show(5)"],"metadata":{"id":"KjjN-q_dw8CA","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"b65becca-61b1-49ed-818d-fa7a2f6c4334"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["#### **Measuring performance on the Test Set: $\\text{RMSE}$ and $\\text{R}^2$ statistic**"],"metadata":{"id":"9fue7UUz5mom","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"9997244d-ac50-45d0-95f5-421226d5fbde"}}},{"cell_type":"code","source":["test_result = lr_model.evaluate(test)\nprint(\"***** Test Set *****\")\nprint(\"RMSE: {:.3f}\".format(test_result.rootMeanSquaredError))\nprint(\"R2: {:.3f}\".format(test_result.r2))\nprint(\"Adjusted R2: {:.3f}\".format(test_result.r2adj))\nprint(\"***** Test Set *****\")"],"metadata":{"id":"0GySC6TS5c_i","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"13b67ee4-f087-476c-8cbe-ec6c9b34bf10"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["## **Summary of the main findings from Data Exploration step**\n\nFrom our exploratory data analysis above, we observed:\n- Both continuous and categorical features (**DONE**);\n- Different scales of feature values;\n- Several outliers on many features (see the boxplots above)."],"metadata":{"id":"kv5bwWvO9HsE","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"b03983ac-a36b-4497-8501-a0f8d651586b"}}},{"cell_type":"markdown","source":["### **Standardize Features**\n\nLet's now consider the second bullet of the list above: feature scales. We can use a feature scaling transformer to standardize all our features to 0-mean and 1-unit of standard deviation. \n\nIn this way, the coefficient $\\theta_0$ (i.e., the _intercept_) will have a nice interpretation, as it can be thought of as the prediction the model would make when _all_ the input features take on their respective mean values (i.e., 0)."],"metadata":{"id":"jAsM9CLq4wU3","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"c5cda769-e213-4f37-993c-fa741ba372e7"}}},{"cell_type":"code","source":["# Let's define a function to standardize all the features\ndef standardize_features(df, input_col=\"features\", with_std=True, with_mean=True):\n\n    from pyspark.ml.feature import StandardScaler\n\n    # 1. Create the StandardScaler\n    scaler = StandardScaler(inputCol=input_col, outputCol=\"std_\"+input_col, withStd=with_std, withMean=with_mean)\n\n    # 2. Compute summary statistics by fitting the StandardScaler\n    scaler_model = scaler.fit(df)\n\n    # 3. Normalize each feature to have 0-mean and unit standard deviation\n    scaled_data = scaler_model.transform(df)\n\n    # 4. Eventually, return both the scaler_model (for future transformations) and the scaled data\n    return scaler_model, scaled_data"],"metadata":{"id":"D_-X9amgfcVJ","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"f15444aa-5ecb-48a1-8dde-dbdf4fad5f8b"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# Call the function above\nscaler_model, scaled_train = standardize_features(train)"],"metadata":{"id":"qzl8-aMo-fzc","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"c701ef61-d5cf-4d42-a97d-25a9c919b12e"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["scaled_train.select([\"features\", \"std_features\", \"label\"]).show(5, truncate=False)"],"metadata":{"id":"QbDvfMuv-pc1","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"940d62f2-dfe4-4cc1-a892-dc72df00f592"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# First of all, let's just setup an OLS linear regressor with no regularization term\nlr = LinearRegression(featuresCol=\"std_features\", labelCol=\"label\")\n# Train the actual model on our training set `scaled_train`\nlr_model = lr.fit(scaled_train)"],"metadata":{"id":"49BjURng_BSf","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"3a7592e5-9b58-45f2-b7ac-4b763c8c98c4"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Intercept ($\\theta_0$) and Coefficients ($\\theta_1, \\ldots, \\theta_n$)**"],"metadata":{"id":"lhOw4gmWyRag","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"61ae0e95-0ec1-45a0-9bd3-8b5a3a8a900c"}}},{"cell_type":"code","source":["print(\"Intercept: \" + str(lr_model.intercept))\nprint(\"Coefficients: \" + str(lr_model.coefficients))"],"metadata":{"id":"KYxq3K2s_dAZ","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"a37ff747-8636-428f-9c6b-c3f0deb53fdf"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["#### **Measuring performance on the Training Set: $\\text{RMSE}$ and $\\text{R}^2$ statistic**"],"metadata":{"id":"GKgUXrfwyVJ3","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"3d5c6b2d-8c54-47d9-b887-d876c1927df6"}}},{"cell_type":"code","source":["training_result = lr_model.summary\nprint(\"***** Training Set *****\")\nprint(\"RMSE: {:.3f}\".format(training_result.rootMeanSquaredError))\nprint(\"R2: {:.3f}\".format(training_result.r2))\nprint(\"Adjusted R2: {:.3f}\".format(training_result.r2adj))\nprint(\"***** Training Set *****\")"],"metadata":{"id":"3m2R-Y9g_gmA","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"9c803801-7311-415f-959c-c7eaeccbb5d5"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# Let's apply the same feature scaling to the test set\nscaled_test = scaler_model.transform(test)"],"metadata":{"id":"_kMA5fuRA9OF","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"9c0652fd-a776-4310-8ed5-987e53c7cddd"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# Compute predictions on the test set and select only the columns of interest\nlr_predictions = lr_model.transform(scaled_test)\nlr_predictions.select(\"features\", \"std_features\", \"prediction\", \"label\").show(5)"],"metadata":{"id":"L3Wvc8u8BJZA","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"1638c864-2795-4837-b1f8-7318419fffb8"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["#### **Measuring performance on the Test Set: $\\text{RMSE}$ and $\\text{R}^2$ statistic**"],"metadata":{"id":"H6-9_TKWylcz","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"ef73519c-2e9a-4f2e-bdd3-01b8c3595fb7"}}},{"cell_type":"code","source":["test_result = lr_model.evaluate(scaled_test)\nprint(\"***** Test Set *****\")\nprint(\"RMSE: {:.3f}\".format(test_result.rootMeanSquaredError))\nprint(\"R2: {:.3f}\".format(test_result.r2))\nprint(\"Adjusted R2: {:.3f}\".format(test_result.r2adj))\nprint(\"***** Test Set *****\")"],"metadata":{"id":"ng34uXbdBSc1","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"d6a6e9b9-99c2-443f-90aa-0e659da01cbe"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Observations**\n\nSince we are using simple OLS linear regression (i.e., no regularization term is included), we obtain *exactly* the same result whether we do feature scaling or not. Of course, the interpretation of the coefficients are much more convenient when features are standardized.\n\nGenerally speaking, the regression coefficient $\\theta_i$ indicates the effect of a change in $x_i$ on the target variable $y$ with all of the other features $x_j$ unchanged. The measurement units of regression coefficient $\\theta_i$ are units of $y$ per unit of $x_i$. For example, if $y$ is the dollar amount of sales and $x_1$ is the number of people in the sales force, then $\\theta_1$ is expressed as units of dollars of sales per person. Now, suppose that the next regression coefficient, $\\theta_2$, is in units of dollars of sales per number of total miles traveled by the sales force. The question of which is more important to sales, staffing level or travel budget, cannot be answered by comparing $\\theta_1$ to $\\theta_2$ because dollars per person and dollars per mile are not directly comparable. Feature standardization solves this problem by expressing each \"unit\" of each feature $x_i$ as a statistical unit equal to one standard deviation.\n\nAnother reason why feature standardization may be helpful is when the optimizer we use to solve OLS is **not** analyitical (i.e., Normal Equation) but iterative (e.g., L-BFGS or Gradient Descent); in this case standardization helps the numerical method to converge faster.\n\nFinally, the effect of feature standardization is instead appreciable whenever we aim to solve a regularized version of the regularization problem (e.g., L1- or L2-norm regularization)."],"metadata":{"id":"aOmr8jNiz9KJ","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"c72f9870-0f7b-4e14-8ce1-0cf06e96295b"}}},{"cell_type":"markdown","source":["## **Summary of the main findings from Data Exploration step**\n\nFrom our exploratory data analysis above, we observed:\n- Both continuous and categorical features (**DONE**);\n- Different scales of feature values (**NOT NECESSARY UNLESS WE USE REGULARIZATION OR NEED A BETTER INTERPRETATION OF COEFFICIENTS**);\n- Several outliers on many features (see the boxplots above)."],"metadata":{"id":"QndOB860B2dm","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"56de2a5c-2af9-4485-bbef-42388f1ebe0a"}}},{"cell_type":"markdown","source":["### **Handling Outliers**\n\nIn this notebook, we will not deal with this specific issue. However, I strongly suggest you to try apply one of the techniques we mentioned in class:\n- outliers removal (_hard_ approach);\n- winsorization (_soft_ approach).\n\nRecall that feature winsorization works as follow: you choose a certain percentile value (e.g., 90) and for each feature $f_i$ you replace any value $x_i$ which falls either below the 5-th percentile or above the 95-th percentile with the 5-th and the 95-th percentile, respectively."],"metadata":{"id":"MZYG0cx4UuEO","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"51e5205f-16d9-4504-acd9-809317637fa9"}}},{"cell_type":"markdown","source":["## **Putting All Together**\n\nIn the following, we try to summarize the whole pipeline making use also of $K$-fold cross validation to get a better estimate of the generalization performance of our model."],"metadata":{"id":"mKWC5LfR6ibL","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"4564e33c-3a24-4954-b3b2-059832abe52b"}}},{"cell_type":"code","source":["# This function defines the general pipeline for linear regression\ndef linear_regression_pipeline(train, \n                               numerical_features, \n                               categorical_features, \n                               target_variable, \n                               with_std=True,\n                               with_mean=True,\n                               k_fold=5):\n\n    from pyspark.ml.feature import StringIndexer, OneHotEncoder, VectorAssembler, StandardScaler\n    from pyspark.ml.tuning import CrossValidator, ParamGridBuilder\n    from pyspark.ml.evaluation import RegressionEvaluator\n    from pyspark.ml.regression import LinearRegression\n    from pyspark.ml import Pipeline\n\n    # Configure a linear regression pipeline, which consists of the following stages: \n    # 1) convert categorical features to numerical ones\n    # 2) standardize feature values (optional)\n    # ... add any other custom transformation here ...\n    # n) fit a linear regressor model\n\n    # 1.a Create a list of indexers, i.e., one for each categorical feature\n    indexers = [StringIndexer(inputCol=c, outputCol=\"{0}_indexed\".format(c), handleInvalid=\"keep\") for c in categorical_features]\n\n    # 1.b Create the one-hot encoder for the list of features just indexed (this encoder will keep any unseen label in the future)\n    encoder = OneHotEncoder(inputCols=[indexer.getOutputCol() for indexer in indexers], \n                                    outputCols=[\"{0}_encoded\".format(indexer.getOutputCol()) for indexer in indexers], \n                                    handleInvalid=\"keep\")\n\n    # 1.c Assemble all the features into a single vector\n    assembler = VectorAssembler(inputCols=encoder.getOutputCols() + numerical_features, outputCol=\"features\")\n\n    # 2.a Create the StandardScaler\n    # scaler = StandardScaler(inputCol=assembler.getOutputCol(), outputCol=\"std_\"+assembler.getOutputCol(), withStd=with_std, withMean=with_mean)\n\n    # ...\n\n    if TARGET_VARIABLE != \"label\":\n        train = train.withColumnRenamed(TARGET_VARIABLE, \"label\")\n\n    lr = LinearRegression(featuresCol=\"features\", labelCol=\"label\") # change `featuresCol=std_features` if scaler is used\n\n    pipeline = Pipeline(stages=indexers + [encoder] + [assembler] + [lr]) # add/remove `[scaler]` to the pipeline if needed\n\n\n    # We use a ParamGridBuilder to construct a grid of parameters to search over.\n    # A CrossValidator requires an Estimator, a set of Estimator ParamMaps, and an Evaluator.\n    # We use a ParamGridBuilder to construct a grid of parameters to search over.\n    # With 3 values for lr.regParam ($\\lambda$) and 3 values for lr.elasticNetParam ($\\alpha$),\n    # this grid will have 3 x 3 = 9 parameter settings for CrossValidator to choose from.\n    param_grid = ParamGridBuilder()\\\n    .addGrid(lr.regParam, [0.0, 0.05, 0.1]) \\\n    .addGrid(lr.elasticNetParam, [0.0, 0.5, 1.0])\\\n    .build()\n    \n    cross_val = CrossValidator(estimator=pipeline, \n                               estimatorParamMaps=param_grid,\n                               evaluator=RegressionEvaluator(metricName=\"rmse\"),\n                               numFolds=k_fold,\n                               collectSubModels=True # this flag allows us to store ALL the models trained during k-fold cross validation\n                               )\n\n    # Run cross-validation, and choose the best set of parameters.\n    cv_model = cross_val.fit(train)\n\n    return cv_model"],"metadata":{"id":"Ux02QmTg7_Ps","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"baf43aa9-49fb-4310-a6d1-98cdbfe4873c"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["cv_model = linear_regression_pipeline(train_df, NUMERICAL_FEATURES, CATEGORICAL_FEATURES, TARGET_VARIABLE)"],"metadata":{"id":"xjnj2jXeRSEJ","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"b0ac96c3-93f4-42d2-87fe-c2189b37927b"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# This function summarizes all the models trained during k-fold cross validation\ndef summarize_all_models(cv_models):\n    for k, models in enumerate(cv_models):\n        print(\"*************** Fold #{:d} ***************\\n\".format(k+1))\n        for i, m in enumerate(models):\n            print(\"--- Model #{:d} out of {:d} ---\".format(i+1, len(models)))\n            print(\"\\tParameters: lambda=[{:.3f}]; alpha=[{:.3f}] \".format(m.stages[-1]._java_obj.getRegParam(), m.stages[-1]._java_obj.getElasticNetParam()))\n            print(\"\\tModel summary: {}\\n\".format(m.stages[-1]))\n        print(\"***************************************\\n\")"],"metadata":{"id":"IHvPkEyEBUOf","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"43d02107-54b7-4e99-8f77-cf3a52770249"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["# Call the function above|\nsummarize_all_models(cv_model.subModels)"],"metadata":{"id":"m50kfzDTRhsd","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"545d6bed-0c65-4a94-8351-741733486fae"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["for i, avg_rmse in enumerate(cv_model.avgMetrics):\n    print(\"Avg. RMSE computed across k-fold cross validation for model setting #{:d}: {:3f}\".format(i+1, avg_rmse))"],"metadata":{"id":"jJyI7BSC8cmA","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"cd72b3e7-86d1-4857-a95d-4d52594d5ee5"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Measuring performance on the Training Set: $\\text{RMSE}$ and $\\text{R}^2$ statistic**"],"metadata":{"id":"BzVrvU9kJOZj","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"eb85a69c-8624-4040-b935-fdfed5e8cfd8"}}},{"cell_type":"code","source":["# `bestModel` is the best resulting model according to K-fold cross validation, which is also entirely retrained on the whole `train_df`\ntraining_result = cv_model.bestModel.stages[-1].summary\nprint(\"***** Training Set *****\")\nprint(\"RMSE: {:.3f}\".format(training_result.rootMeanSquaredError))\nprint(\"R2: {:.3f}\".format(training_result.r2))\nprint(\"Adjusted R2: {:.3f}\".format(training_result.r2adj))\nprint(\"***** Training Set *****\")"],"metadata":{"id":"fLgNpJP0JVtw","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"53f1d247-f22a-4bc5-961e-cbc41f60e6d8"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Using the best model from $K$-fold cross validation to make predictions**"],"metadata":{"id":"yadIWG4pIcrr","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"dfbd6f0e-7f6f-490d-9519-bbcecc8bffb7"}}},{"cell_type":"code","source":["# Make predictions on the test set (`cv_model` contains the best model according to the result of k-fold cross validation)\n# `test_df` will follow exactly the same pipeline defined above, and already fit to `train_df`\ntest_df = test_df.withColumnRenamed(\"price\", \"label\")\npredictions = cv_model.transform(test_df)"],"metadata":{"id":"b01W9KqbSz6c","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"2f069066-4c10-49e1-9b3d-bd2f9ed992c1"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["predictions.select(\"features\", \"prediction\", \"label\").show(5)"],"metadata":{"id":"2Tlh8G3STxGh","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"beae054f-a0df-4ca8-94d1-f9fb1e39fe03"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["def evaluate_model(predictions, metric=\"rmse\"):\n    \n    from pyspark.ml.evaluation import RegressionEvaluator\n\n    evaluator = RegressionEvaluator(labelCol=\"label\",\n                                    predictionCol=\"prediction\",\n                                    metricName=metric)\n\n    return evaluator.evaluate(predictions)"],"metadata":{"id":"FOA-EKGxU16J","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"c4280880-db40-4cb4-813f-6e4e55e3bbca"}},"outputs":[],"execution_count":0},{"cell_type":"code","source":["def r2_adj(predictions):\n    \n    r2 = evaluate_model(predictions, metric=\"r2\")\n    r2_adj_score = (1 - (1 - r2) * ((predictions.count() - 1) / (predictions.count() - predictions.select('features').first()[0].size - 1)))\n\n    return r2_adj_score"],"metadata":{"id":"bsZ2IQ7lOtbQ","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"76e699c0-2b1d-4095-8cd2-db55f68f00f2"}},"outputs":[],"execution_count":0},{"cell_type":"markdown","source":["### **Measuring performance on the Test Set: $\\text{RMSE}$ and $\\text{R}^2$ statistic**"],"metadata":{"id":"7ulP5vi_VpnP","colab_type":"text","application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"05996c75-e1cc-4e08-aca7-d03756c26ad8"}}},{"cell_type":"code","source":["print(\"***** Test Set *****\")\nprint(\"RMSE: {:.3f}\".format(evaluate_model(predictions)))\nprint(\"R2: {:.3f}\".format(evaluate_model(predictions, metric=\"r2\")))\nprint(\"Adjusted R2: {:.3f}\".format(r2_adj(predictions)))\nprint(\"***** Test Set *****\")"],"metadata":{"id":"WDEu2HuaaueC","colab_type":"code","colab":{},"application/vnd.databricks.v1+cell":{"title":"","showTitle":false,"inputWidgets":{},"nuid":"6715704e-2b16-4228-bf60-7a0591aaf83c"}},"outputs":[],"execution_count":0}],"metadata":{"colab":{"name":"Lecture_09_Linear_Regression.ipynb","provenance":[],"collapsed_sections":[],"machine_shape":"hm"},"kernelspec":{"name":"python3","display_name":"Python 3"},"application/vnd.databricks.v1+notebook":{"notebookName":"Linear_Regression","dashboards":[],"notebookMetadata":{"pythonIndentUnit":2},"language":"python","widgets":{},"notebookOrigID":1470792745644131}},"nbformat":4,"nbformat_minor":0}