.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/plot_california_housing.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_plot_california_housing.py: Improve your data science workflow with skops ============================================= .. GENERATED FROM PYTHON SOURCE LINES 7-9 Introduction ------------ .. GENERATED FROM PYTHON SOURCE LINES 11-14 The goal of this exercise is to go through a semi-realistic data science and machine learning task and develop a practical solution for it. We will learn about the following topics: .. GENERATED FROM PYTHON SOURCE LINES 16-24 - Perform *exploratory data analysis* - Do some non-trivial *feature engineering* - Explain how the feature engineering informs the *choice of machine learning model* and vice versa - Show how to make use of a couple of *advanced scikit-learn* features and explain why we use them - Create a *model card* that provides useful information about the model - Share the model by uploading it to the *Hugging Face Hub* .. GENERATED FROM PYTHON SOURCE LINES 26-28 Imports ------- .. GENERATED FROM PYTHON SOURCE LINES 30-33 Before we start, we need to import a couple of packages. In particular, to run this code, we need to have the following 3rd party packages: jupyter, matplotlib, pandas, scikit-learn, skops. .. GENERATED FROM PYTHON SOURCE LINES 35-37 So if you want to run this exercise yourself and these packages are not installed yet in your Python environment, you should run: .. GENERATED FROM PYTHON SOURCE LINES 39-40 ``python -m pip install jupyter matplotlib pandas scikit-learn skops`` .. GENERATED FROM PYTHON SOURCE LINES 42-73 .. code-block:: Python import os from operator import itemgetter from pathlib import Path from tempfile import mkdtemp import matplotlib.pyplot as plt import numpy as np import pandas as pd import sklearn from matplotlib.patches import Rectangle from sklearn.compose import ColumnTransformer from sklearn.datasets import fetch_california_housing from sklearn.dummy import DummyRegressor from sklearn.ensemble import ( GradientBoostingRegressor, RandomForestRegressor, StackingRegressor, ) from sklearn.inspection import DecisionBoundaryDisplay, permutation_importance from sklearn.linear_model import LinearRegression from sklearn.metrics import get_scorer from sklearn.model_selection import GridSearchCV, cross_val_predict, train_test_split from sklearn.neighbors import KNeighborsRegressor from sklearn.pipeline import Pipeline from sklearn.preprocessing import FunctionTransformer from sklearn.tree import DecisionTreeRegressor import skops from skops import card, hub_utils from skops import io as sio .. GENERATED FROM PYTHON SOURCE LINES 74-76 .. code-block:: Python plt.style.use("seaborn-v0_8") .. GENERATED FROM PYTHON SOURCE LINES 77-79 Analyzing the dataset --------------------- .. GENERATED FROM PYTHON SOURCE LINES 81-83 Fetch the data ~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 85-89 First of all, let’s load our dataset. For this exercise, we use the California Housing dataset. It can be downloaded using the function from scikit-learn. If called the first time, this function will download the dataset, on subsequent calls, it will load the cached version. .. GENERATED FROM PYTHON SOURCE LINES 91-94 We will make use of the option to set , which will return the data as a pandas. This will make our much easier than working with a numpy array, which it would return otherwise. .. GENERATED FROM PYTHON SOURCE LINES 96-98 .. code-block:: Python data = fetch_california_housing(as_frame=True) .. GENERATED FROM PYTHON SOURCE LINES 99-101 The dataset comes with a description. If not already familiar with the dataset, it’s always a good idea to read the included description. .. GENERATED FROM PYTHON SOURCE LINES 103-105 .. code-block:: Python print(data.DESCR) .. rst-class:: sphx-glr-script-out .. code-block:: none .. _california_housing_dataset: California Housing dataset -------------------------- **Data Set Characteristics:** :Number of Instances: 20640 :Number of Attributes: 8 numeric, predictive attributes and the target :Attribute Information: - MedInc median income in block group - HouseAge median house age in block group - AveRooms average number of rooms per household - AveBedrms average number of bedrooms per household - Population block group population - AveOccup average number of household members - Latitude block group latitude - Longitude block group longitude :Missing Attribute Values: None This dataset was obtained from the StatLib repository. https://www.dcc.fc.up.pt/~ltorgo/Regression/cal_housing.html The target variable is the median house value for California districts, expressed in hundreds of thousands of dollars ($100,000). This dataset was derived from the 1990 U.S. census, using one row per census block group. A block group is the smallest geographical unit for which the U.S. Census Bureau publishes sample data (a block group typically has a population of 600 to 3,000 people). A household is a group of people residing within a home. Since the average number of rooms and bedrooms in this dataset are provided per household, these columns may take surprisingly large values for block groups with few households and many empty houses, such as vacation resorts. It can be downloaded/loaded using the :func:`sklearn.datasets.fetch_california_housing` function. .. rubric:: References - Pace, R. Kelley and Ronald Barry, Sparse Spatial Autoregressions, Statistics and Probability Letters, 33 (1997) 291-297 .. GENERATED FROM PYTHON SOURCE LINES 106-108 Exploratory data analysis ~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 110-112 Now it’s time to start exploring the dataset. First of all, let’s determine what the target for this task is: .. GENERATED FROM PYTHON SOURCE LINES 114-117 .. code-block:: Python target_col = data.target_names[0] print(target_col) .. rst-class:: sphx-glr-script-out .. code-block:: none MedHouseVal .. GENERATED FROM PYTHON SOURCE LINES 118-122 The target column is called "MedHouseVal" and from the description, we know it designates "the median house value for California districts, expressed in hundreds of thousands of dollars". .. GENERATED FROM PYTHON SOURCE LINES 124-126 Next let’s extract the actual data, which, as mentioned, is contained in a pandas ``DataFrame``: .. GENERATED FROM PYTHON SOURCE LINES 128-130 .. code-block:: Python df = data["frame"] .. GENERATED FROM PYTHON SOURCE LINES 131-135 For now, we leave the target variable inside the ``DataFrame``, as this will facilitate the upcoming analysis. Once we get to modeling, we should of course separate the target data to avoid accidentally training on the target. .. GENERATED FROM PYTHON SOURCE LINES 137-138 Let’s peak at some properties of the data. .. GENERATED FROM PYTHON SOURCE LINES 140-142 .. code-block:: Python df.shape .. rst-class:: sphx-glr-script-out .. code-block:: none (20640, 9) .. GENERATED FROM PYTHON SOURCE LINES 143-145 .. code-block:: Python df.head() .. raw:: html
MedInc HouseAge AveRooms AveBedrms Population AveOccup Latitude Longitude MedHouseVal
0 8.3252 41.0 6.984127 1.023810 322.0 2.555556 37.88 -122.23 4.526
1 8.3014 21.0 6.238137 0.971880 2401.0 2.109842 37.86 -122.22 3.585
2 7.2574 52.0 8.288136 1.073446 496.0 2.802260 37.85 -122.24 3.521
3 5.6431 52.0 5.817352 1.073059 558.0 2.547945 37.85 -122.25 3.413
4 3.8462 52.0 6.281853 1.081081 565.0 2.181467 37.85 -122.25 3.422


.. GENERATED FROM PYTHON SOURCE LINES 146-148 .. code-block:: Python df.describe() .. raw:: html
MedInc HouseAge AveRooms AveBedrms Population AveOccup Latitude Longitude MedHouseVal
count 20640.000000 20640.000000 20640.000000 20640.000000 20640.000000 20640.000000 20640.000000 20640.000000 20640.000000
mean 3.870671 28.639486 5.429000 1.096675 1425.476744 3.070655 35.631861 -119.569704 2.068558
std 1.899822 12.585558 2.474173 0.473911 1132.462122 10.386050 2.135952 2.003532 1.153956
min 0.499900 1.000000 0.846154 0.333333 3.000000 0.692308 32.540000 -124.350000 0.149990
25% 2.563400 18.000000 4.440716 1.006079 787.000000 2.429741 33.930000 -121.800000 1.196000
50% 3.534800 29.000000 5.229129 1.048780 1166.000000 2.818116 34.260000 -118.490000 1.797000
75% 4.743250 37.000000 6.052381 1.099526 1725.000000 3.282261 37.710000 -118.010000 2.647250
max 15.000100 52.000000 141.909091 34.066667 35682.000000 1243.333333 41.950000 -114.310000 5.000010


.. GENERATED FROM PYTHON SOURCE LINES 149-151 Scaling the target ^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 153-158 Before we continue working with the data, let’s do a small adjustment. From the description, we know that the target variable, the median house price, is expressed in units of $100,000. For our first row, that means the actual price is $52,600, not $52.6 (which would be very cheap, even for 1990). .. GENERATED FROM PYTHON SOURCE LINES 160-164 In theory, the unit should not matter for our work, but let’s still convert it to $. This is because when we search for general solutions to this task, most people work with $ values. If we use a different unit here, it makes comparison to these results unnecessarily difficult. .. GENERATED FROM PYTHON SOURCE LINES 166-168 .. code-block:: Python df["MedHouseVal"] = 100_000 * df["MedHouseVal"] .. GENERATED FROM PYTHON SOURCE LINES 169-171 Differences to other versions of the dataset ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 173-181 Another notable difference in this particular dataset is that some columns are already averages. E.g. we find the column AveRooms, which is the average number of rooms per household of this group of houses. In other versions of the dataset (like the one on `kaggle `__), we will, however, find the total number of rooms and the number of households. So here, some feature engineering was already performed by calculating the averages. This is fine and we can keep it like this. .. GENERATED FROM PYTHON SOURCE LINES 183-185 Missing values ^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 187-190 Furthermore, in the description, we find: Missing Attribute Values: None. This probably means that there are no missing values, but let’s check ourselves just to be certain: .. GENERATED FROM PYTHON SOURCE LINES 193-195 .. code-block:: Python df.isna().any() .. rst-class:: sphx-glr-script-out .. code-block:: none MedInc False HouseAge False AveRooms False AveBedrms False Population False AveOccup False Latitude False Longitude False MedHouseVal False dtype: bool .. GENERATED FROM PYTHON SOURCE LINES 196-199 Indeed, the dataset contains no missing values. If it did, we could have made use of the `imputation features of sklearn `__. .. GENERATED FROM PYTHON SOURCE LINES 201-203 Distributions of variables ^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 205-209 It’s always useful to take a look at the distributions of our variables. This is best achieved by visual inspection, which is why we will create a couple of plots. For this exercise, we use matplotlib for plotting, since it’s very powerful and popular. It also works well with pandas. .. GENERATED FROM PYTHON SOURCE LINES 211-212 First we want to focus on the following columns: .. GENERATED FROM PYTHON SOURCE LINES 215-217 .. code-block:: Python cols = ["MedInc", "HouseAge", "AveRooms", "AveBedrms", "Population", "AveOccup"] .. GENERATED FROM PYTHON SOURCE LINES 218-221 These are all our feature variables except for the geospatial features, "Longitude" and "Latitude". We will analyze those more closely later. .. GENERATED FROM PYTHON SOURCE LINES 223-225 The most basic plot we can create for analyzing the distribution is the histogram. So let’s plot the histograms for each of those variables. .. GENERATED FROM PYTHON SOURCE LINES 227-232 Before we take a closer look at the data, just a few words on how we create the plots. Since we have 6 variables, it would be convenient to plot the data in a 3x2 grid. That’s why createa matplotlib figure with 6 subplots, using 3 rows and 2 columns. The resulting ``axes`` variable is a 3x2 numpy array that contains the individual subplots. .. GENERATED FROM PYTHON SOURCE LINES 234-238 We also want to make use of the pandas plotting method, which we can call using ``df.plot``. This uses matplotlib under the hood, so it’s not strictly needed. But the nice thing is that pandas provides some extra convenience, e.g. by automatically labeling the axes. .. GENERATED FROM PYTHON SOURCE LINES 240-245 In order for pandas to plot onto our created figure with its subplots, we pass the ``X`` argument to ``df.plot``. This tells pandas to plot onto this subplot, instead of creating a new plot. Another little trick is to ``flatten`` the subplot array while looping over it. That way, we don’t need to take care of looping over its two dimensions separately. .. GENERATED FROM PYTHON SOURCE LINES 247-249 Finally, we should also call ``plt.tight_layout()`` to prevent the subplots from overlapping. .. GENERATED FROM PYTHON SOURCE LINES 251-257 .. code-block:: Python fig, axes = plt.subplots(3, 2, figsize=(8, 12)) for ax, col in zip(axes.flatten(), cols): df.plot(kind="hist", y=col, bins=100, title=col, ax=ax, legend=None) plt.tight_layout() .. image-sg:: /auto_examples/images/sphx_glr_plot_california_housing_001.png :alt: MedInc, HouseAge, AveRooms, AveBedrms, Population, AveOccup :srcset: /auto_examples/images/sphx_glr_plot_california_housing_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 258-263 Already, we can make some interesting observations about the data. For "MedInc", but especially for "HouseAge", we see a large bin to the right side of the distribution. This could be an indicator that values might have been clipped. Let’s look at "HouseAge" more closely: .. GENERATED FROM PYTHON SOURCE LINES 265-267 .. code-block:: Python df["HouseAge"].describe() .. rst-class:: sphx-glr-script-out .. code-block:: none count 20640.000000 mean 28.639486 std 12.585558 min 1.000000 25% 18.000000 50% 29.000000 75% 37.000000 max 52.000000 Name: HouseAge, dtype: float64 .. GENERATED FROM PYTHON SOURCE LINES 268-270 .. code-block:: Python df["HouseAge"].value_counts().head().to_frame("count") .. raw:: html
count
HouseAge
52.0 1273
36.0 862
35.0 824
16.0 771
17.0 698


.. GENERATED FROM PYTHON SOURCE LINES 271-277 So we see that there are 1273 samples with a house age of 52, which also happens to be the maximum value. Could this be coincidence? Perhaps, but it’s unlikely. The more likely explanation is that houses that were older than 52 years were just clipped at 52. In the context of the problem we’re trying to solve, this doesn’t make a big difference, so we will just accept this peculiarity. .. GENERATED FROM PYTHON SOURCE LINES 279-283 Next we can see in the histograms above that for "AveRooms", "AveBedrms", "Population", and "AveOccup", the bins are squished to the left. This means that there is a fat right tail in the distribution, which might be problematic. When looking at the description, we find a potential explanation: .. GENERATED FROM PYTHON SOURCE LINES 285-290 An household is a group of people residing within a home. Since the average number of rooms and bedrooms in this dataset are provided per household, these columns may take surpinsingly large values for block groups with few households and many empty houses, such as vacation resorts. .. GENERATED FROM PYTHON SOURCE LINES 293-304 So what should we do about this? For the purpose of plotting the values, we should certainly think about removing these extreme values. For the machine learning model we will train later, the answer is: it depends. When we use a model like linear regressions or a neural network, these extreme values can be problematic and it would make sense to scale the values to have a more uniform or normal distribution. When using decision tree-based models, these exteme values are not problematic though. A decision tree will split the data into those samples that are less than or greater than a certain value – it doesn’t matter how much smaller or greater the values are. Since we will actually rely on tree-based models, let’s leave the data as is (except for plotting). .. GENERATED FROM PYTHON SOURCE LINES 306-314 When it comes to plotting, how can we deal with these extreme values? We could scale the data, e.g. by taking a ``log``. This can be achieved by passing ``logx=True`` to the plotting method. However, the log scale makes it harder to read the actual values of the data. Instead, let’s use a more brute force approach of simply excluding that 1% largest values. For this, we calculate the 99th percentile of the values and exclude all values that exceed that percentile. Apart from that change, the plots are the same as above: .. GENERATED FROM PYTHON SOURCE LINES 317-324 .. code-block:: Python fig, axes = plt.subplots(3, 2, figsize=(8, 12)) for ax, col in zip(axes.flatten(), cols): quantile = df[col].quantile(0.99) mask = df[col] < quantile df[mask].plot(kind="hist", y=col, bins=100, title=col, ax=ax, legend=None) plt.tight_layout() .. image-sg:: /auto_examples/images/sphx_glr_plot_california_housing_002.png :alt: MedInc, HouseAge, AveRooms, AveBedrms, Population, AveOccup :srcset: /auto_examples/images/sphx_glr_plot_california_housing_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 325-328 With extremely large values excluded, we can see that the variables are distributed in a very reasonable fashion. With a bit of squinting, the distributions almost look Gaussian, which is what we like to see. .. GENERATED FROM PYTHON SOURCE LINES 330-332 Correlations with the target ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 334-339 Now it’s time to look at how our features are correlated with our target. After all, we plan on predicting the target based on the features, and even though correlation with a target is not necessary for a feature to be helpful, it’s a strong indicator. As before, we will filter out extreme values. .. GENERATED FROM PYTHON SOURCE LINES 341-353 .. code-block:: Python fig, axes = plt.subplots(3, 2, figsize=(8, 12)) for ax, col in zip(axes.flatten(), cols): quantile = df[col].quantile(0.99) df_subset = df[df[col] < quantile] correlation = df_subset[[col, target_col]].corr().iloc[0, 1] title = f"Pearson correlation coefficient: {correlation:.3f}" df_subset.plot( kind="scatter", x=col, y=target_col, s=1.5, alpha=0.05, ax=ax, title=title ) plt.tight_layout() .. image-sg:: /auto_examples/images/sphx_glr_plot_california_housing_003.png :alt: Pearson correlation coefficient: 0.672, Pearson correlation coefficient: 0.040, Pearson correlation coefficient: 0.329, Pearson correlation coefficient: -0.090, Pearson correlation coefficient: -0.033, Pearson correlation coefficient: -0.280 :srcset: /auto_examples/images/sphx_glr_plot_california_housing_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 354-358 Again, let’s try to improve our understanding of the dataset by visually inspecting the outcomes. The first thing to notice is perhaps that our target, the "MedHouseVal", seems to be clipped at $500,000. Let’s remember that and return to it later. .. GENERATED FROM PYTHON SOURCE LINES 360-364 Next we should notice that apart from "MedInc", none of the variables seem to strongly correlate with the target. If this were a real business problem to solve at a company, now would be a good time to get ahold of a domain expert and verify that this is expected. .. GENERATED FROM PYTHON SOURCE LINES 366-369 We also calculate the Pearson correlation coefficient and show it in the plot titles, but honestly, it cannot tell us much we can’t already determine by visual inspection. .. GENERATED FROM PYTHON SOURCE LINES 371-373 Geospatial features ^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 375-379 Now it’s time to take a look at the geospatial data, namely "Longitude" and "Latitude". We already know that our dataset is limited to the US state of California, so we should expect the data to be exclusively in that area. .. GENERATED FROM PYTHON SOURCE LINES 381-385 Before actually plotting the data, let’s take a quick brake and form a hypothesis. We can reasonably expect that housing prices should be high in metropolitan areas. For California, this could be around Los Angeles and the bay area. Does our data reflect that? .. GENERATED FROM PYTHON SOURCE LINES 387-393 To answer this question, we can plot the target variable as a function of its coordinates. Since we deal with longitude and latitude, and since it’s reasonably likely that the Earth is not flat, it is not quite correct to just plot the data as is. It would be more accurate to use a projection to map the coordinates to 2 dimensions, e.g. by using `geopandas `__. .. GENERATED FROM PYTHON SOURCE LINES 395-398 For our purposes, however, despite the size of California, we should be able to get away without using any projections. So let’s simplify our life by just using raw longitude and latitude. .. GENERATED FROM PYTHON SOURCE LINES 400-417 .. code-block:: Python fig, ax = plt.subplots(figsize=(10, 8)) df.plot( kind="scatter", x="Longitude", y="Latitude", c=target_col, title="House value by location", cmap="coolwarm", s=2.5, ax=ax, ) inset = (-122.5, 37.5) rect = Rectangle( inset, 0.5, 0.5, linewidth=1, edgecolor="k", facecolor="none", alpha=0.5 ) ax.add_patch(rect) .. image-sg:: /auto_examples/images/sphx_glr_plot_california_housing_004.png :alt: House value by location :srcset: /auto_examples/images/sphx_glr_plot_california_housing_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 418-419 This is an interesting plot. Let’s try to draw some conclusions. .. GENERATED FROM PYTHON SOURCE LINES 421-426 First of all, we see that the locations are not evenly distributed across the state. There are big patches without any data, especially on the eastern side of the state, whereas the western coast is more highly populated. Data is also more sparse in the mountainous area, where we can expect population density to be lower. .. GENERATED FROM PYTHON SOURCE LINES 428-433 Regarding our initial hypothesis about house prices, we can indeed find clusters of high prices exceeding $400,000, whereas the majority of the data points seem to fall below $50,000. After checking on a map, the high priced areas seem to be indeed around Los Angeles and the bay area, but there are other high priced areas as well, e.g. in San Diego. .. GENERATED FROM PYTHON SOURCE LINES 435-443 From this, we can already start thinking about some possibly interesting features. Some variations of this dataset contain, for instance, a variable indicating the closeness to the Ocean, since it seems that areas on the coast are more expensive on average. Another interesting feature could be the distance to key cities like Los Angeles and San Francisco. Here distance should probably be measured in terms of how long it takes to drive there by car, not just purely in terms of spatial distance – after all, most people only have a car and not a helicopter. .. GENERATED FROM PYTHON SOURCE LINES 445-450 As you can imagine, getting this type of data requires additional data sources and probably a lot of extra work. Therefore, we won’t take this route. However, we will develop specific geospatial features below, which will probably capture most of the information we need for this task. .. GENERATED FROM PYTHON SOURCE LINES 452-456 Before advancing further, we should zoom into the geospatial data, given how difficult it is to see details on the plot above. For that, we take a small inset of that figure (marked with a rectangle), which seems to be an interesting spot, and plot it below: .. GENERATED FROM PYTHON SOURCE LINES 459-472 .. code-block:: Python fig, ax = plt.subplots(figsize=(10, 8)) df.plot( kind="scatter", x="Longitude", y="Latitude", c=target_col, title="House value by location", cmap="coolwarm", ax=ax, ) ax.set_xlim([inset[0], inset[0] + 0.5]) ax.set_ylim([inset[1], inset[1] + 0.5]) .. image-sg:: /auto_examples/images/sphx_glr_plot_california_housing_005.png :alt: House value by location :srcset: /auto_examples/images/sphx_glr_plot_california_housing_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none (37.5, 38.0) .. GENERATED FROM PYTHON SOURCE LINES 473-479 What can we take away from this plot? First of all, it shows that even though there are clear patterns of high and low price areas, they can also be situated very close to each other. This could already give us the idea that prices in the *neighborhood* could be a strong predictor for our target variable, but the size of the neighborhood should not be too large. .. GENERATED FROM PYTHON SOURCE LINES 481-489 Another interesting pattern we see is that the data points are distributed on a regular grid. This is not what we would expect if, say, each point represented a community, town, or city, since those are not evenly spaced. Instead, this indicates that the data was aggregated with a certain spatial resolution. Also, given the gaps, it could be reasonable to assume that data points with too few houses were removed from the dataset, maybe for privacy concerns. Again, talking to a domain expert would help us better understand the reason. .. GENERATED FROM PYTHON SOURCE LINES 491-492 Anyway, we should keep this regular spatial distribution in mind for later. .. GENERATED FROM PYTHON SOURCE LINES 494-498 A final observation about the coordinates. We might believe that for a given longitude and latitude, there is exactly one sample (or 0, if not present). However, there are duplicates when it comes to coordinates, as shown below: .. GENERATED FROM PYTHON SOURCE LINES 500-502 .. code-block:: Python df.duplicated(["Longitude", "Latitude"]).sum() .. rst-class:: sphx-glr-script-out .. code-block:: none np.int64(8050) .. GENERATED FROM PYTHON SOURCE LINES 503-505 These rows are not completely duplicated, however, as there are no duplicates when considering all columns: .. GENERATED FROM PYTHON SOURCE LINES 507-509 .. code-block:: Python df.duplicated().sum() .. rst-class:: sphx-glr-script-out .. code-block:: none np.int64(0) .. GENERATED FROM PYTHON SOURCE LINES 510-511 At the most extreme end, we find coordinates with 15 samples: .. GENERATED FROM PYTHON SOURCE LINES 513-515 .. code-block:: Python df[["Longitude", "Latitude"]].value_counts().to_frame("count").reset_index().head() .. raw:: html
Longitude Latitude count
0 -122.41 37.80 15
1 -122.42 37.80 11
2 -122.44 37.78 11
3 -122.41 37.75 10
4 -122.43 37.77 10


.. GENERATED FROM PYTHON SOURCE LINES 516-518 Overall, for more than a third of coordinates, we find more than one data point: .. GENERATED FROM PYTHON SOURCE LINES 520-522 .. code-block:: Python (df[["Longitude", "Latitude"]].value_counts() > 1).mean() .. rst-class:: sphx-glr-script-out .. code-block:: none np.float64(0.3457505957108816) .. GENERATED FROM PYTHON SOURCE LINES 523-531 Again, if this were a real business case, we should investigate this further and find someone who can explain this peculiarity in the data. One possible explanation could be that those are samples for the same location but from different points in time. If this were true, we would have to make adjustments, e.g. by considering the time when we make a train/test split. But we have no possibility to verify that time is the explanation. As is, we just accept this fact and keep it in mind for later. .. GENERATED FROM PYTHON SOURCE LINES 533-535 Target variable ^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 537-539 Finally, we should not forget to take a look at the target variable itself. Again, let’s start with plotting its disstribution: .. GENERATED FROM PYTHON SOURCE LINES 541-544 .. code-block:: Python df.plot(kind="hist", y=target_col, bins=100) .. image-sg:: /auto_examples/images/sphx_glr_plot_california_housing_006.png :alt: plot california housing :srcset: /auto_examples/images/sphx_glr_plot_california_housing_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 545-547 As we have already established earlier, we find that the target data seems to be clipped. Let’s take a closer look: .. GENERATED FROM PYTHON SOURCE LINES 549-551 .. code-block:: Python df[target_col].value_counts().head() .. rst-class:: sphx-glr-script-out .. code-block:: none MedHouseVal 500001.0 965 137500.0 122 162500.0 117 112500.0 103 187500.0 93 Name: count, dtype: int64 .. GENERATED FROM PYTHON SOURCE LINES 552-558 Is it possible that this would occur naturally in the data? Sometimes, we may find strange patterns. To give an example, if there was a law that sales above a certain value are taxed differently, we could expect prices to cluster at this value. But this appears to be very unlikely here, especially since prices seem to be rounded to the closest $100, as we can see from the following probe: .. GENERATED FROM PYTHON SOURCE LINES 560-569 .. code-block:: Python ( (df[target_col] % 100) .round() .value_counts() .to_frame() .reset_index() .rename(columns={"index": f"{target_col} % $100", target_col: "count"}) ) .. raw:: html
count count
0 0.0 18396
1 100.0 1275
2 1.0 965
3 99.0 4


.. GENERATED FROM PYTHON SOURCE LINES 570-575 So if we take the modulo of the house price to $100, we find that it it’s almost always 0 (well, we see a few 100’s, but that’s just a rounding issue). Then we see 965 prices ending in , which are exactly those 965 samples we found with a price of $100,001. Finally, we see 4 prices ending in 9. Let’s take a look: .. GENERATED FROM PYTHON SOURCE LINES 577-579 .. code-block:: Python df[np.isclose(df[target_col] % 100, 99)][target_col].to_frame() .. raw:: html
MedHouseVal
2521 14999.0
2799 14999.0
9188 14999.0
19802 14999.0


.. GENERATED FROM PYTHON SOURCE LINES 580-582 .. code-block:: Python df[target_col].min() .. rst-class:: sphx-glr-script-out .. code-block:: none np.float64(14999.000000000002) .. GENERATED FROM PYTHON SOURCE LINES 583-588 So these four samples actually correspond to the lowest prices in the dataset. Therefore, it is very reasonable to assume that the dataset creators decided to set a maximum price of $500,000 and a minimum price of $15,000, with all prices falling outside that range being set to the max/min price +/- . For the prices within the range, they decided to round to $100. .. GENERATED FROM PYTHON SOURCE LINES 590-599 When it comes to our task of predicting the target, this clipping could be dangerous. We cannot know by how much the actual price was clipped, especially when it comes to high prices. For a machine learning model, it could become extraordinarily hard to predict these high prices, because even though the features may, for instance, indicate a price of $1,000,000 for one sample and $500,100 for another, the model is supposed to predict the same price for both. For this reason, let us remove the clipped data from the dataset. Whether that’s a good idea will in reality depend on the use case that we’re trying to solve. .. GENERATED FROM PYTHON SOURCE LINES 601-603 Training a machine learning model --------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 605-608 Now that we have gained a good understanding of our data, we move to the next phase, which is training a machine learning model to predict the target. .. GENERATED FROM PYTHON SOURCE LINES 610-612 Remove samples with clipped target data ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 614-616 As a first step, we will remove the samples where the target data has been clipped, as discussed above. .. GENERATED FROM PYTHON SOURCE LINES 619-625 .. code-block:: Python mask = (15_000 <= df[target_col]) & (df[target_col] <= 500_000) print( f"Discarding {(1 - mask).sum()} ({100 * (1 - mask.mean()):.1f}%) of rows because" " the target is clipped." ) .. rst-class:: sphx-glr-script-out .. code-block:: none Discarding 969 (4.7%) of rows because the target is clipped. .. GENERATED FROM PYTHON SOURCE LINES 626-628 Train/test split ~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 630-638 Then we should split our data into a training and a test set. As every data scientist should know, we need to evaluate the data on a different set than we use for training. It would be even better if we created three sets of data, train/valid/test, and only used the test set at the very end. For the training part, we could use cross-validation to get more reliable results, given that our dataset is not that big. For the purpose of this exercise, we make our lifes simple by only using a single train/test split though. .. GENERATED FROM PYTHON SOURCE LINES 640-645 As to the split itself, we just split the data randomly using sklearn’s ``train_test_split``. There is nothing in the data that would suggest we need to perform a non-random split but again, this will vary from use case to use case. Note that we set the ``random_state`` to make the results reproducible. .. GENERATED FROM PYTHON SOURCE LINES 648-650 .. code-block:: Python df_train, df_test = train_test_split(df[mask], random_state=0) .. GENERATED FROM PYTHON SOURCE LINES 651-653 .. code-block:: Python df_train.shape, df_test.shape .. rst-class:: sphx-glr-script-out .. code-block:: none ((14753, 9), (4918, 9)) .. GENERATED FROM PYTHON SOURCE LINES 654-658 After performing the split, it is now a good time to remove the target data from the rest of the ``DataFrame``, so that we don’t run the risk of accidentally training on the target. We use the ``pop`` method for that. .. GENERATED FROM PYTHON SOURCE LINES 661-664 .. code-block:: Python y_train = df_train.pop(target_col).values y_test = df_test.pop(target_col).values .. GENERATED FROM PYTHON SOURCE LINES 665-667 Feature engineering ~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 669-672 Now let’s get to the feature engineering part. Here it’s worth it to think a little bit about what type of ML model we plan to use and how that decision should inform the feature engineering. .. GENERATED FROM PYTHON SOURCE LINES 674-680 In general, we know that for tabular tasks as this one, ensembles of decision trees perform exceptionally well, without the need for a lot of tuning. Examples would be random forest or gradient boosting. We should not try to challenge this conventional wisdom, this class of models is an excellent choice for this task as well. There is, however, a caveat. Let’s explore it more closely. .. GENERATED FROM PYTHON SOURCE LINES 682-688 We know that the geospatial features are crucial for our task, since we saw that there is a very strong relationship between the house price of a given sample and the house price of its neighbors. On the other hand, we saw that all the other variables, safe for "MedInc", are probably not strong predictors. We thus need to ensure we can make the best use of "Longitude" and "Latitude". .. GENERATED FROM PYTHON SOURCE LINES 690-699 If we use a tree-based model, can we just input the longitude and latitude as features and be done with it? In theory, this would indeed work. However, let’s remember how decision trees work. According to a certain criterion, they split the data at a certain value. As an example, we could find that the first split of a tree is to split the data into samples with latitude less than 34 and latitude greater than 34, which is roughly the median. Each time we split the data along longitude and latitude, we can imagine the map being devided into four quadrants. So far, so good. .. GENERATED FROM PYTHON SOURCE LINES 701-705 The problem now comes with the patchiness of the data. Looking again at the house value by location plot further above, we can immediately see that we would need *a lot of splits* to capture the neighborhood relationship we find in the data. How many splits would we need? Let’s take a look. .. GENERATED FROM PYTHON SOURCE LINES 707-713 To study this question, we will take the ``DecisionTreeRegressor`` from sklearn and fit it on the longitude and latitude features. The number of splits is bounded by the ``max_depth`` parameter. So for each level of depth, the tree splits the data once. That is, for a depth of 10, we get ``2**10=1024`` splits (in reality, the number could be lower, depending on the other parameters of the tree). .. GENERATED FROM PYTHON SOURCE LINES 715-721 To get a feeling of what that means in practice, let us first fit 4 decision trees, using ``max_depth`` values of 1, 2, 5, and 10. Then we plot the *decision boundary* of the tree after fitting it. We use sklearn’s `DecisionBoundaryDisplay `__ to plot the data. Here is what we get: .. GENERATED FROM PYTHON SOURCE LINES 723-741 .. code-block:: Python _, axes = plt.subplots(2, 2, figsize=(8, 8)) max_depths = [1, 2, 5, 10] for max_depth, ax in zip(max_depths, axes.flatten()): dt = DecisionTreeRegressor(random_state=0, max_depth=max_depth) dt.fit(df_train[["Longitude", "Latitude"]], y_train) DecisionBoundaryDisplay.from_estimator( dt, df_train[["Longitude", "Latitude"]], cmap="coolwarm", response_method="predict", ax=ax, xlabel="Longitude", ylabel="Latitude", grid_resolution=1000, ) ax.set_title(f"Decision boundary for a depth of {max_depth}") plt.tight_layout() .. image-sg:: /auto_examples/images/sphx_glr_plot_california_housing_007.png :alt: Decision boundary for a depth of 1, Decision boundary for a depth of 2, Decision boundary for a depth of 5, Decision boundary for a depth of 10 :srcset: /auto_examples/images/sphx_glr_plot_california_housing_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 742-744 As we have described, with a depth of 1, we split our map in half, with two splits, we get 4 quadrants. This is shown in the first row. .. GENERATED FROM PYTHON SOURCE LINES 746-751 In the second row, we see the outcomes for depths of 5 and 10. Especially for 10, we see this resulting in quite a few "boxes" that represent neighborhoods with similar prices. But even with this relatively high number of splits, the "resolution" of the resulting map is quite bad, lumping together many areas with quite different prices. .. GENERATED FROM PYTHON SOURCE LINES 753-762 So what can we do about that? The easiest solution would be to increase the ``max_depth`` parameter to such a high value that we can actually model the spatial relationship well enough. But this has some disadvantages. First of all, a high ``max_depth`` parameter can result in overfitting on the training data. Remember that we also have other variables that we want to include, the tree will also split on those. Second of all, a high ``max_depth`` makes the model bigger and slower. Depending on the application, this could be a problem for productionizing the model. .. GENERATED FROM PYTHON SOURCE LINES 764-770 Another solution is that we could use an ensemble of decision trees. The result of that will be multiple maps as those above layered on top of each other. Although this will certainly help, it’s still not a perfect solution, since it would still require quite a lot of trees to achieve a good resolution. Imagine a dataset with much more samples and much more fine-grained data, we would need a giant ensemble to fit it. .. GENERATED FROM PYTHON SOURCE LINES 772-780 At the end of the day, we have to admit that decision tree-based models are just not the best fit when it comes to modeling geospatial relationships. Can we think of a model that is better able to model this type of data? Why, of course we can. The k-nearest neighbor (KNN) family of models should be a perfect fit for this, since we want to model *neighborhood* relationships. To see this in action, let’s again plot the decision boundaries, this time using the ``KNeighborsRegressor`` from sklearn: .. GENERATED FROM PYTHON SOURCE LINES 783-784 this controls the level of parallelism, feel free to set to a higher number .. GENERATED FROM PYTHON SOURCE LINES 784-786 .. code-block:: Python N_JOBS = 1 .. GENERATED FROM PYTHON SOURCE LINES 787-802 .. code-block:: Python _, ax = plt.subplots(figsize=(5, 5)) knn = KNeighborsRegressor(n_jobs=N_JOBS) knn.fit(df_train[["Longitude", "Latitude"]], y_train) DecisionBoundaryDisplay.from_estimator( knn, df_train[["Longitude", "Latitude"]], cmap="coolwarm", response_method="predict", ax=ax, xlabel="Longitude", ylabel="Latitude", grid_resolution=1000, ) ax.set_title("Decision boundary of KNN") .. image-sg:: /auto_examples/images/sphx_glr_plot_california_housing_008.png :alt: Decision boundary of KNN :srcset: /auto_examples/images/sphx_glr_plot_california_housing_008.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 1.0, 'Decision boundary of KNN') .. GENERATED FROM PYTHON SOURCE LINES 803-812 Now if we compare this to the decision boundaries of the decision trees, even with a depth of 10, it’s a completely different picture. The KNN model, by its very nature, can easily model very fine grained spatial differences. The granularity will depend on the ``k`` part of the name, which indicates the number of neighbors that are used to make the prediction. At the one extreme, for ``k=1``, we would only consider a single data point, resulting in a very spotty map. At the other extreme, when ``k`` is the size of the total dataset, we would simply average across all data points. Choosing a good value for ``k`` is thus important. .. GENERATED FROM PYTHON SOURCE LINES 814-818 We will now go into more details of the KNN model. If you’re wondering how this is related to feature engineering, it will become clear later, as will actually use the KNN model for the purpose of creating a new feature. So please be patient. .. GENERATED FROM PYTHON SOURCE LINES 820-822 Aside: distance metrics ^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 824-828 Before we explore this problem further, let’s talk about distance metrics. For a KNN model to work, we need to define the distance metric it uses to determine the closest neighbors. By default, ``KNeighborsRegressor`` uses the Euclidean distance. .. GENERATED FROM PYTHON SOURCE LINES 830-836 Earlier, we saw that our data points are distributed on a regular grid. This means, if we take a specific sample, and if we assume that it’s neighboring spots are not empty, there should be 4 neighbors at exactly the same distance, namely the neighbors directly north, east, south, and west. Similarly, neighbors 5 to 8, in the directions NE, SE, SW, and NW, should also have the exact same distances. .. GENERATED FROM PYTHON SOURCE LINES 838-840 Let’s check if this is true. For this, we fit a KNN and then use the ``kneighbors`` method to return the closest neighbors. .. GENERATED FROM PYTHON SOURCE LINES 843-846 .. code-block:: Python knn = KNeighborsRegressor(20) knn.fit(df[["Longitude", "Latitude"]], df[target_col]) .. raw:: html 
KNeighborsRegressor(n_neighbors=20)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.


.. GENERATED FROM PYTHON SOURCE LINES 847-852 .. code-block:: Python distances = knn.kneighbors( df[["Longitude", "Latitude"]].iloc[[123]], return_distance=True )[0].round(5) print(distances) .. rst-class:: sphx-glr-script-out .. code-block:: none [[0. 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01414 0.01414 0.01414 0.01414 0.01414 0.01414 0.01414]] .. GENERATED FROM PYTHON SOURCE LINES 853-860 So we find indeed that the distances are discrete, with many samples having exactly the same distance like 0.01. The closest neighbor has a distance of 0 – this is not surprising, as this the sample itself. However, unlike what we expected, we don’t find exactly four equally distant closest neighbors. Instead, for this sample we find 12 samples with a distance of exaclty 0.01 (i.e. exactly one step to the north, east, south, or west). How come? .. GENERATED FROM PYTHON SOURCE LINES 862-865 Remember from earlier that we found duplicate coordinates in our data? This is most likely such a case. So let’s remove duplicates and determine the distances again: .. GENERATED FROM PYTHON SOURCE LINES 868-872 .. code-block:: Python knn = KNeighborsRegressor(20) df_no_dup = df.drop_duplicates(["Longitude", "Latitude"]).copy() knn.fit(df_no_dup[["Longitude", "Latitude"]], df_no_dup[target_col]) .. raw:: html 
KNeighborsRegressor(n_neighbors=20)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.


.. GENERATED FROM PYTHON SOURCE LINES 873-878 .. code-block:: Python distances = knn.kneighbors( df[["Longitude", "Latitude"]].iloc[[123]], return_distance=True )[0].round(5) print(distances) .. rst-class:: sphx-glr-script-out .. code-block:: none [[0. 0.01 0.01 0.01 0.01 0.01414 0.01414 0.01414 0.01414 0.02 0.02 0.02 0.02 0.02236 0.02236 0.02236 0.02236 0.02236 0.02236 0.02828]] .. GENERATED FROM PYTHON SOURCE LINES 879-881 This looks much better. Now we find exactly 4 neighbors tied for the closest distance, 4 tied for the second closest distance, etc. .. GENERATED FROM PYTHON SOURCE LINES 883-889 The same thing should become even clearer when we change the metric from Euclidean to Manhattan distance. In this context, the Manhattan distance basically means how far two points are from each other, if we can only take steps along the north-south or east-west axis. To calculate this metric, we can set the ``p`` parameter of ``KNeighborsRegressor`` to 1: .. GENERATED FROM PYTHON SOURCE LINES 891-894 .. code-block:: Python knn = KNeighborsRegressor(20, p=1) knn.fit(df_no_dup[["Longitude", "Latitude"]], df_no_dup[target_col]) .. raw:: html 
KNeighborsRegressor(n_neighbors=20, p=1)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.


.. GENERATED FROM PYTHON SOURCE LINES 895-899 .. code-block:: Python distances = knn.kneighbors( df[["Longitude", "Latitude"]].iloc[[123]], return_distance=True )[0].round(5) print(distances) .. rst-class:: sphx-glr-script-out .. code-block:: none [[0. 0.01 0.01 0.01 0.01 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.03 0.03 0.03 0.03 0.03 0.03 0.03]] .. GENERATED FROM PYTHON SOURCE LINES 900-905 Again, we find 4 neighbors with a distance of exactly 0.01, as expected. For distance 0.02, we actually find 8 neighbors. This makes sense, because in Manhattan space, the neighbor north-north of the sample has the same distance as the neighbor north-east, etc., as both require two steps to be reached. .. GENERATED FROM PYTHON SOURCE LINES 907-908 What does this mean for us? There are two important considerations: .. GENERATED FROM PYTHON SOURCE LINES 910-925 1. When we train a KNN model with a ``k`` that’s not exactly equal to 4, 8, etc., we run into some problems. If we take, for instance, ``k=3``, and the 4 closest neighbors are equally close, ``KNeighborsRegressor`` will actually pick an arbitrary set of 3 from those 4 possible candidates. This generally something we want to avoid in our ML models. In practice, however, the problem isn’t so bad, because we have duplicates and because some neighbor spots are empty, as we saw earlier. This makes it less likely that many data points exhibit different behavior for a very specific ``k``. 2. The second consideration is that we should figure out which metric is best for our use case, Euclidean or Manhattan. If people were traveling by air, Euclidean makes most sense. If they traveled on a rectangular grid (like the streets in Manhattan, hence the name), Manhattan makes more sense. In practice, we have neither of those. So let’s just use what works better! .. GENERATED FROM PYTHON SOURCE LINES 927-929 Determining the best hyper-parameters for the KNN model ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 931-935 With all that in mind, let’s try to find good hyper-parameters for our KNN regressor. As mentioned, we should check the value for ``k`` and we should check ``p=1`` and ``p=2`` (i.e. Manhattan and Euclidean distance). # .. GENERATED FROM PYTHON SOURCE LINES 937-944 On top of those, let’s check one more hyper-parameter, namely ``weights``. What this means is that when the model finds the ``k`` nearest neighbors, should it simply predict the average target of those values or should closer neighbors have a higher weight than more remote neighbors? To use the former approach, we can set ``weights=’uniform’``, for the latter, we set it to ``’distance’``. .. GENERATED FROM PYTHON SOURCE LINES 946-949 As for the ``k`` parameter, in sklearn, it’s called ``n_neighbors``. Let’s check the values from 1 to 25, as well as some higher values, 50, 75, 100, 150, and 200. .. GENERATED FROM PYTHON SOURCE LINES 951-956 Regarding the metrics, we use root mean squared error (RMSE). This could also be replaced with other metrics, it really depends on the use case. Here we choose it because this is the most common one used for this dataset, so if we wanted to compare our results with the results from others, it makes sense to use the same metrics. .. GENERATED FROM PYTHON SOURCE LINES 958-962 (Note: For scoring, we use the ``’neg_root_mean_squared_error’``, i.e. the negative RMSE. This is because by convention, sklearn considers higher values to be better. For RMSE, however, lower values are better. To circumvent that, we just the negative RMSE.) .. GENERATED FROM PYTHON SOURCE LINES 964-969 To check all the different parameter combinations and to calculate the RMSE out of fold, we rely on sklearn’s ``GridSearchCV``. We won’t explain what it does here, since there already are so many tutorials out there that discuss grid search. Suffice it to say, this is exactly what we need for our problem. .. GENERATED FROM PYTHON SOURCE LINES 971-980 .. code-block:: Python knn = KNeighborsRegressor() params = { "weights": ["uniform", "distance"], "p": [1, 2], "n_neighbors": list(range(1, 26)) + [50, 75, 100, 150, 200], } search = GridSearchCV(knn, params, scoring="neg_root_mean_squared_error") search.fit(df_train[["Longitude", "Latitude"]], y_train) .. raw:: html 
GridSearchCV(estimator=KNeighborsRegressor(),
                 param_grid={'n_neighbors': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
                                             13, 14, 15, 16, 17, 18, 19, 20, 21, 22,
                                             23, 24, 25, 50, 75, 100, 150, 200],
                             'p': [1, 2], 'weights': ['uniform', 'distance']},
                 scoring='neg_root_mean_squared_error')
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.


.. GENERATED FROM PYTHON SOURCE LINES 981-986 Fortunately, the grid search only takes a few seconds because the model is quick, the dataset size is small, and we only test 120 different parameter combinations. If the grid search would take too long, we could switch to ``RandomizedSearchCV`` or ``HalvingGridSearchCV`` from sklearn, or use a Bayesian optimization method from other packages. .. GENERATED FROM PYTHON SOURCE LINES 988-990 Now, let’s put the results of the grid search into a pandas ``DataFrame`` and inspect the top results. .. GENERATED FROM PYTHON SOURCE LINES 993-998 .. code-block:: Python df_cv = pd.DataFrame(search.cv_results_) df_cv.sort_values("rank_test_score")[ ["param_n_neighbors", "param_weights", "param_p", "mean_test_score"] ].head(10) .. raw:: html
param_n_neighbors param_weights param_p mean_test_score
22 6 uniform 2 -48792.277422
26 7 uniform 2 -48908.663856
18 5 uniform 2 -48939.112420
20 6 uniform 1 -48960.566567
24 7 uniform 1 -49049.354866
30 8 uniform 2 -49125.740942
16 5 uniform 1 -49158.915396
28 8 uniform 1 -49264.347864
34 9 uniform 2 -49299.733688
32 9 uniform 1 -49463.991846


.. GENERATED FROM PYTHON SOURCE LINES 999-1003 Okay, so there is a clear trend for small ``k`` values in the range of 5 to 9 to fare better. Also ``’uniform’`` weights seem to be clearly superior. When it comes to the distance metric, no clear trend is discernable. .. GENERATED FROM PYTHON SOURCE LINES 1005-1012 To get a better picture, it is often a good idea to plot the grid search results, instead of just picking the best result blindly, especially if the outcome is noisy. For this purpose we will plot four lines, one for each combination of ``weights`` and ``p``, all showing the RMSE as a function of ``n_neighbors``. Note that we plot ``n_neighbors`` on a log scale, because we want to have higher resolution for small values. .. GENERATED FROM PYTHON SOURCE LINES 1014-1031 .. code-block:: Python fig, ax = plt.subplots() for weight in params["weights"]: # type: ignore for p in params["p"]: # type: ignore query = f"param_weights=='{weight}' & param_p=={p}" df_subset = df_cv.query(query) df_subset.plot( x="param_n_neighbors", y="mean_test_score", xlabel="Log of the n_neighbors parameter", ylabel="negative RMSE", label=query, ax=ax, marker="o", ms=5, logx=True, ) .. image-sg:: /auto_examples/images/sphx_glr_plot_california_housing_009.png :alt: plot california housing :srcset: /auto_examples/images/sphx_glr_plot_california_housing_009.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 1032-1036 So both for ``weights==’uniform’`` and for ``weights==’distance’``, we see that ``p==2`` is slightly better than, or equally good as, ``p==1``. It will thus be safe to use ``p==2``. .. GENERATED FROM PYTHON SOURCE LINES 1038-1045 Furthermore, we see that for low values of ``n_neighbors``, it is better to have ``weights==’uniform’``, while for large values of ``n_neighbors``, ``weights==’distance’`` is better. This is perhaps not too surprising: If we consider a lot of neighbors, many of which are already quite far away, we should give lower weight to those remote neighbors. When we only look at 5 neighbors, this isn’t necessary – their distances will be quite similar anway. .. GENERATED FROM PYTHON SOURCE LINES 1047-1052 Another observation is that the curve for the RMSE as a function of ``n_neighbors`` seems to be quite smooth. This tells us two things: First of all, the metric isn’t very noisy, at least using the 5-fold cross validation that the grid seach uses by default. This is good, we don’t want the outcome to be too dependent on randomness. .. GENERATED FROM PYTHON SOURCE LINES 1054-1059 Moreover, we don’t see any strange jumps at ``n_neighbors`` values like 4 or 8. If you remember the problem we discussed earlier of arbitrary points being chosen when neighbors have the exact same distance, this could have been a concern. Since we don’t see it, it’s probably not as bad as we might have feared. .. GENERATED FROM PYTHON SOURCE LINES 1061-1065 So now that we have a good idea about what good hyper-parameters for KNN are, let’s see how the good the predictions from the model are. For this, we plot the true target as a function of the predictions, which we calculate out of fold using sklearn’s ``cross_val_predict``. .. GENERATED FROM PYTHON SOURCE LINES 1068-1077 .. code-block:: Python knn = KNeighborsRegressor(**search.best_params_) avg_neighbor_val = cross_val_predict(knn, df_train[["Longitude", "Latitude"]], y_train) fig, ax = plt.subplots(figsize=(5, 5)) ax.scatter(avg_neighbor_val, y_train, alpha=0.05, s=1.5) ax.set_xlabel( f"median house price of {search.best_params_['n_neighbors']} closest neighbors" ) ax.set_ylabel(target_col) .. image-sg:: /auto_examples/images/sphx_glr_plot_california_housing_010.png :alt: plot california housing :srcset: /auto_examples/images/sphx_glr_plot_california_housing_010.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(-5.402777777777777, 0.5, 'MedHouseVal') .. GENERATED FROM PYTHON SOURCE LINES 1078-1085 As we can see, the correlation is quite good between the two. This means that our KNN model can explain a lot of the variability in the target using only the coordinates, none of the other features. In particular, we can easily see that this correlation is much better than the correlation of any of the features that we plotted above. **Would it not be nice if we could use the KNN predictions as a feature for a more powerful model?** .. GENERATED FROM PYTHON SOURCE LINES 1087-1093 Why would we want to use the KNN predictions as a feature for another model, instead of just being content with using the KNN as the final predictor? The problem is that with KNN, it is very hard to incorporate the remaining features. Although we can be sure that those features are not as important, they should still contain useful signal to improve the model. .. GENERATED FROM PYTHON SOURCE LINES 1095-1101 The reason why it is difficult to include arbitrary features like "MedInc" or "HouseAge" into a KNN model is that they are on a completely different scale than the coordinates, so we would need to normalize all the data. Even then, we know that our other features are not very uniformely distributed, which in general is something that the KNN benefits from. .. GENERATED FROM PYTHON SOURCE LINES 1103-1108 But even if everything was on the same scale and evenly distributed, should we weight a distance of 0.1 in the "MedInc" space the same as a distance of 0.1 in "Longitude" space? Probably not. It just doesn’t make sense to put two completely measurements into the same KNN model. .. GENERATED FROM PYTHON SOURCE LINES 1110-1115 This is why we will use the KNN predictions as a *feature* for more appropriate models. This also explains why this is part of the feature engineering section. Doing this correctly is not quite trivial, but below we will see how we can use the tools that sklearn provides to achieve this goal. .. GENERATED FROM PYTHON SOURCE LINES 1117-1119 Training ML models ~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 1121-1124 Now that we have gained important insights into the data and also formed a plan when it comes to feature engineering, we can start with the training of the predictive models. .. GENERATED FROM PYTHON SOURCE LINES 1126-1128 Dummy model ^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 1130-1135 As a first step, we start with a dummy model, i.e. a model that isn’t allowed to learn anything from the input data. Often, it’s a good idea to try a dummy model first, as it will give us an idea what the worst score is we should expect. If our proper models cannot beat the dummy model, it most likely means we have done something seriously wrong. .. GENERATED FROM PYTHON SOURCE LINES 1137-1144 For this dataset, what is the appropriate dummy model? We are trying to minimize the root mean squared error. This is the same as trying to minimize the mean squared error (the root is just a monotonic transformation). And to minimize the mean squared error, if we don’t know anything about the input data, requires us to predict the *mean* of the target. For our convenience, sklearn provides a ``DummyRegressor`` that will do exactly this for us: .. GENERATED FROM PYTHON SOURCE LINES 1146-1149 .. code-block:: Python dummy = DummyRegressor() dummy.fit(df_train, y_train) .. raw:: html 
DummyRegressor()
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.


.. GENERATED FROM PYTHON SOURCE LINES 1150-1153 Note: Even though we pass ``df_train``, the ``DummyRegressor`` does not make use of it. We could pass an empty ``df`` and the result would be the same. .. GENERATED FROM PYTHON SOURCE LINES 1156-1158 .. code-block:: Python -get_scorer("neg_root_mean_squared_error")(dummy, df_test, y_test) .. rst-class:: sphx-glr-script-out .. code-block:: none 100528.43327419003 .. GENERATED FROM PYTHON SOURCE LINES 1159-1161 As a reminder, sklearn only comes with the negative RMSE, so we make it positive again at the end. .. GENERATED FROM PYTHON SOURCE LINES 1163-1165 Adding KNN predictions as features ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 1167-1174 Let’s get to the part where we want to add the predictions of the KNN model as features for our final regression model. As we have discussed above, we have some good evidence that the KNN predictions are well suited to capture the spatial relationship between data points, but KNN itself is not good with dealing with the other features, so we want to use a model that can deal with those features on top of the KNN predictions. .. GENERATED FROM PYTHON SOURCE LINES 1176-1179 Does sklearn offer us a way of adding the predictions of one model as features for another model? Indeed it does, this is achieved by using the :class:`sklearn.ensemble.StackingRegressor`. To quote from the docs: .. GENERATED FROM PYTHON SOURCE LINES 1181-1185 *Stacked generalization consists in stacking the output of individual estimator [sic] and use a regressor to compute the final prediction. Stacking allows to use the strength of each individual estimator by using their output as input of a final estimator.* .. GENERATED FROM PYTHON SOURCE LINES 1187-1190 In our case, that means that we use KNN for the "individual estimator" part and then a different model, like an ensemble of decision trees, for the "final esimator". .. GENERATED FROM PYTHON SOURCE LINES 1192-1193 The documentation also states that: .. GENERATED FROM PYTHON SOURCE LINES 1195-1198 *Note that estimators\_ are fitted on the full X while final_estimator\_ is trained using cross-validated predictions of the base estimators using ``cross_val_predict``."* .. GENERATED FROM PYTHON SOURCE LINES 1200-1211 This is important to know. It is crucial that the predictions from the KNNs are calculated using ``cross_val_predict``, i.e *out of fold*, otherwise, we run into the risk of overfitting. To understand this point, let’s take an extreme example. Let’s say we use a KNN with ``n_neighbors=1``. If the predictions were calculated *in fold*, this KNN would completely overfit on the training data and the prediction would be perfect (okay, not quite, since we have duplicate coordinates in the data). Then the final model would only rely on the seemingly perfect KNN predictions, ignoring all other features. This would be really bad. That’s why the "cross-validated predictions" part is so crucial. .. GENERATED FROM PYTHON SOURCE LINES 1214-1218 .. code-block:: Python # Note: We may be tempted to simply add the new feature to the ``DataFrame`` # like this: ``df['pred_knn'] = knn.predict(df[['Longitude', 'Latitude']])``. .. GENERATED FROM PYTHON SOURCE LINES 1219-1223 However, this is problematic. First of all, we need to make out of fold predictions, as explained above. This code snippet would result in overfitting. We thus need to calculate the out of fold predictions manually, which adds more (errror prone) custom code. .. GENERATED FROM PYTHON SOURCE LINES 1225-1231 Second, let’s assume we want to deploy the final model. When we call it, we need to pass all the features, i.e. we would need to generate the KNN prediction before passing the features to the model. This requires even more custom code to be added, making the whole application even more error prone. By sticking with the tools sklearn gives us, we avoid the two issues. .. GENERATED FROM PYTHON SOURCE LINES 1233-1239 By default, ``StackingRegressor`` trains the final estimator only on the predictions of the individual estimators. However, we want it to be trained on all the other features too. This can be achieved by setting ``StackingRegressor(..., passthrough=True)``, which will *pass through* the original input and concatenate it with the predictions from our KNN. .. GENERATED FROM PYTHON SOURCE LINES 1241-1248 Another issue we need to solve is that we want the final estimator to be trained on all features, but the KNN is supposed to be only trained on longitude and latitude. When we pass all features as ``X``, the KNN would be trained on all these features, which, as we discussed, wouldn’t be a good idea. If we only pass longitude and latitude as ``X``, the final estimator cannot make use of the other features. What do we do? .. GENERATED FROM PYTHON SOURCE LINES 1250-1261 The solution here is to pass all the features, but to put the KNN into a ``Pipeline`` with one step *selecting* only the longitude and latitude, and the second step being the KNN itself. Unless we’re missing something, sklearn does not directly provide a transformer that is only used for selecting columns, but we can cobble one together ourselves. For this, we choose the ``FunctionTransformer``, which is a transformer that calls an arbitrary function. The function itself should simpy select the two columns ``["Longitude", "Latitude"]``. This can be done using the ``itemgetter`` function from the builtin ``operator`` library. The resulting ``Pipeline`` looks like this: .. GENERATED FROM PYTHON SOURCE LINES 1263-1270 .. code-block:: Python Pipeline( [ ("select_cols", FunctionTransformer(itemgetter(["Longitude", "Latitude"]))), ("knn", KNeighborsRegressor()), ] ) .. raw:: html 
Pipeline(steps=[('select_cols',
                     FunctionTransformer(func=operator.itemgetter(['Longitude', 'Latitude']))),
                    ('knn', KNeighborsRegressor())])
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.


.. GENERATED FROM PYTHON SOURCE LINES 1271-1277 Note: Alternatively, we could have made use of sklearn’s ``ColumnTransformer``, which does have a builtin way of selecting columns. It also wants to apply some transformer to these selected columns, even though we don’t need that. This can be circumvented by setting this transformer to ``"passthrough"``. The end result is: .. GENERATED FROM PYTHON SOURCE LINES 1279-1291 .. code-block:: Python Pipeline( [ ( "select_cols", ColumnTransformer( [("long_and_lat", "passthrough", ["Longitude", "Latitude"])] ), ), ("knn", KNeighborsRegressor()), ] ) .. raw:: html 
Pipeline(steps=[('select_cols',
                     ColumnTransformer(transformers=[('long_and_lat', 'passthrough',
                                                      ['Longitude', 'Latitude'])])),
                    ('knn', KNeighborsRegressor())])
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.


.. GENERATED FROM PYTHON SOURCE LINES 1292-1294 At the end of the day, it doesn’t really matter which variation we use. Let’s go with the 2nd approach. .. GENERATED FROM PYTHON SOURCE LINES 1296-1303 With this out of the way, what hyper-parameters do we want to use for our KNN? In our grid search of the KNN regressor, we found that small values of ``n_neighbors`` work best. As to the other hyper-parameters, we already saw that there is no point in changing the ``p`` parameter from the default value of 2 to 1. For the ``weights`` parameter, we saw that low ``n_neighbors`` work better with ``uniform``, the default, so let’s use that here. .. GENERATED FROM PYTHON SOURCE LINES 1305-1309 By the way, ``StackingRegressor`` can also take multiple estimators for the initial prediction. Therefore, we could pass a list of multiple KNNs with different hyper-parameters. This would be a good idea to test if we want to further improve the models. .. GENERATED FROM PYTHON SOURCE LINES 1311-1313 Linear regression ^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 1315-1324 Now let’s finally get started with training and evaluating our first ML model on the whole dataset. As a start, we try out a simple linear regression, which often provides a good benchmark for regression tasks. Usually, with linear regressions, we would like to normalize the data first, but sklearn’s ``LinearRegression`` already does this for us, so we don’t need to bother with that. (Given what we found out about the distributions of some of the features as shown in the earlier plot with the feature histograms, it would, however, be worth to spend some time thinking about whether we could improve upon the default preprocessing.) .. GENERATED FROM PYTHON SOURCE LINES 1326-1329 The ``StackingRegressor`` expects a list of tuples, where the first element of the tuple is a name and the second element is the estimator. Plugging our different parts together, we get: .. GENERATED FROM PYTHON SOURCE LINES 1332-1349 .. code-block:: Python knn_regressor = [ ( "knn@5", Pipeline( [ ( "select_cols", ColumnTransformer( [("long_and_lat", "passthrough", ["Longitude", "Latitude"])] ), ), ("knn", KNeighborsRegressor()), ] ), ), ] .. GENERATED FROM PYTHON SOURCE LINES 1350-1356 .. code-block:: Python lin = StackingRegressor( estimators=knn_regressor, final_estimator=LinearRegression(n_jobs=N_JOBS), passthrough=True, ) .. GENERATED FROM PYTHON SOURCE LINES 1357-1359 .. code-block:: Python lin.fit(df_train, y_train) .. raw:: html 
StackingRegressor(estimators=[('knn@5',
                                   Pipeline(steps=[('select_cols',
                                                    ColumnTransformer(transformers=[('long_and_lat',
                                                                                     'passthrough',
                                                                                     ['Longitude',
                                                                                      'Latitude'])])),
                                                   ('knn',
                                                    KNeighborsRegressor())]))],
                      final_estimator=LinearRegression(n_jobs=1), passthrough=True)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.


.. GENERATED FROM PYTHON SOURCE LINES 1360-1362 .. code-block:: Python -get_scorer("neg_root_mean_squared_error")(lin, df_test, y_test) .. rst-class:: sphx-glr-script-out .. code-block:: none 44316.26978449049 .. GENERATED FROM PYTHON SOURCE LINES 1363-1368 The final score is already quite an improvement over the results from the dummy model, so we can be happy about that. Moreover, if we compare this score the score we got when we trained a KNN purely on longitude and latitude, it’s also much better, which confirms our decision that using the other features is helpful for the model. .. GENERATED FROM PYTHON SOURCE LINES 1370-1373 When comparing the results from other people, it also doesn’t look too bad, but we have to keep in mind that the datasets and preprocessing steps are not identical, so differences should be expected. .. GENERATED FROM PYTHON SOURCE LINES 1375-1377 Just out of curiosity, let’s check the score without using the KNN predictions as features: .. GENERATED FROM PYTHON SOURCE LINES 1379-1381 .. code-block:: Python lin_raw = LinearRegression(n_jobs=N_JOBS) .. GENERATED FROM PYTHON SOURCE LINES 1382-1384 .. code-block:: Python lin_raw.fit(df_train, y_train) .. raw:: html 
LinearRegression(n_jobs=1)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.


.. GENERATED FROM PYTHON SOURCE LINES 1385-1387 .. code-block:: Python -get_scorer("neg_root_mean_squared_error")(lin_raw, df_test, y_test) .. rst-class:: sphx-glr-script-out .. code-block:: none 66647.75686249179 .. GENERATED FROM PYTHON SOURCE LINES 1388-1393 We see a quite substantial increase in the prediction error. This shouldn’t be too surprising. When the linear regressor tries to fit longitude and latitude, the only thing it can do is try to fit a plane on top of it, which is far too simple to fit the geospatial patterns we observed. .. GENERATED FROM PYTHON SOURCE LINES 1395-1397 Random forest ^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 1399-1402 Next let’s use our first decision tree-based model, the ``RandomForestRegressor``. For this, we use the same approach as above, we only need to swap the ``final_estimator``: .. GENERATED FROM PYTHON SOURCE LINES 1404-1412 .. code-block:: Python rf = StackingRegressor( estimators=knn_regressor, final_estimator=RandomForestRegressor( n_estimators=100, random_state=0, n_jobs=N_JOBS ), passthrough=True, ) .. GENERATED FROM PYTHON SOURCE LINES 1413-1415 .. code-block:: Python rf.fit(df_train, y_train) .. raw:: html 
StackingRegressor(estimators=[('knn@5',
                                   Pipeline(steps=[('select_cols',
                                                    ColumnTransformer(transformers=[('long_and_lat',
                                                                                     'passthrough',
                                                                                     ['Longitude',
                                                                                      'Latitude'])])),
                                                   ('knn',
                                                    KNeighborsRegressor())]))],
                      final_estimator=RandomForestRegressor(n_jobs=1,
                                                            random_state=0),
                      passthrough=True)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.


.. GENERATED FROM PYTHON SOURCE LINES 1416-1418 .. code-block:: Python -get_scorer("neg_root_mean_squared_error")(rf, df_test, y_test) .. rst-class:: sphx-glr-script-out .. code-block:: none 42136.84558918467 .. GENERATED FROM PYTHON SOURCE LINES 1419-1421 We can see a nice, but not huge, improvement over using the ``LinearRegressor``. .. GENERATED FROM PYTHON SOURCE LINES 1423-1425 Gradient boosted decision trees ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 1427-1428 Finally, let’s use sklearn’s ``GradientBoostingRegressor``: .. GENERATED FROM PYTHON SOURCE LINES 1430-1436 .. code-block:: Python gb = StackingRegressor( estimators=knn_regressor, final_estimator=GradientBoostingRegressor(n_estimators=100, random_state=0), passthrough=True, ) .. GENERATED FROM PYTHON SOURCE LINES 1437-1439 .. code-block:: Python gb.fit(df_train, y_train) .. raw:: html 
StackingRegressor(estimators=[('knn@5',
                                   Pipeline(steps=[('select_cols',
                                                    ColumnTransformer(transformers=[('long_and_lat',
                                                                                     'passthrough',
                                                                                     ['Longitude',
                                                                                      'Latitude'])])),
                                                   ('knn',
                                                    KNeighborsRegressor())]))],
                      final_estimator=GradientBoostingRegressor(random_state=0),
                      passthrough=True)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.


.. GENERATED FROM PYTHON SOURCE LINES 1440-1442 .. code-block:: Python -get_scorer("neg_root_mean_squared_error")(gb, df_test, y_test) .. rst-class:: sphx-glr-script-out .. code-block:: none 41826.816445595665 .. GENERATED FROM PYTHON SOURCE LINES 1443-1448 The score is almost identical to the ``RandomForestRegressor``, as is the training time. For other choices of hyper-parameters, this will certainly differ, but as is, it doesn’t really matter which model we choose. Here it would be a good exercise to perform a hyper-parameter search to determine what model is truly better. .. GENERATED FROM PYTHON SOURCE LINES 1450-1452 Aside: Checking the importance of longitude and latitude as predictive features ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 1454-1459 Remember that earlier on, we formed the hypothesis that tree-based models would have difficulties making use of the longitude and latitude features because they require a high amount of splits to really be useful? Even if it’s not strictly necessary, let’s take some time to check this hypothesis. .. GENERATED FROM PYTHON SOURCE LINES 1461-1467 To do this, what we would want to do is to check how much the tree-based model relies on said features, given that we change the depth of the tree. To determine the importance of the feature, we will use the :func:`sklearn.inspection.permutation_importance` function from sklearn. Then we would like to check if longitude and latitude become more important if we increase the depth of the decision trees. .. GENERATED FROM PYTHON SOURCE LINES 1469-1473 Be aware that we don’t want to use the ``StackingRegressor`` with the KNN predictions here, because for that model, the longitude and latitude influence the KNN prediction feature, obfuscating the result. So let’s use the pure ``GradientBoostingRegressor`` here. .. GENERATED FROM PYTHON SOURCE LINES 1475-1480 Our approach will be to train 3 models with different values for ``max_depth``. We choose relative small values here because gradient boosting typically uses very shallow trees (the default depth is 3). After training, we calculate the permutation importances of each feature and create a bar plot of the results: .. GENERATED FROM PYTHON SOURCE LINES 1482-1498 .. code-block:: Python max_depths = [2, 4, 6] fig, axes = plt.subplots(1, 3, figsize=(12, 8)) for md, ax in zip(max_depths, axes): gb = GradientBoostingRegressor(max_depth=md, random_state=0) gb.fit(df_train, y_train) pi = permutation_importance(gb, df_train, y_train, random_state=0) score = -get_scorer("neg_root_mean_squared_error")(gb, df_test, y_test) ax.barh(df_train.columns, pi["importances_mean"]) ax.set_xlim([0, 1.5]) title = f"permutation importances for max_depths={md}\n(test RMSE: {score:.0f})" ax.set_title(title) if md > max_depths[0]: ax.set_yticklabels([]) plt.tight_layout() .. image-sg:: /auto_examples/images/sphx_glr_plot_california_housing_011.png :alt: permutation importances for max_depths=2 (test RMSE: 54063), permutation importances for max_depths=4 (test RMSE: 47080), permutation importances for max_depths=6 (test RMSE: 44283) :srcset: /auto_examples/images/sphx_glr_plot_california_housing_011.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 1499-1507 As we can see, longitude and latitude are not very important for ``max_depth=2``, even being below the importance of "MedInc". In contrast, they are more important for ``max_depth=4``, with the importance growing even more for ``max_depth=6``. This confirms our hypothesis, though we should be aware that besides ``max_depth``, other hyper-parameters influence the effective number of splits we have, so in reality it’s not that simple. .. GENERATED FROM PYTHON SOURCE LINES 1509-1514 Out of curiosity, we also show the RMSE on the test set for the individual models. Interestingly, we find that its considerably worse than the scores we got earlier when we included the KNN predictions, not even beating the linear regression! This is a nice validation that our KNN feature really helps a lot. .. GENERATED FROM PYTHON SOURCE LINES 1516-1518 Final model ^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 1520-1526 Let’s settle on a final model for now. We will use gradient boosting again, only this time using more estimators. In general, with gradient boosting, more trees help more. The tradeoff is mostly that the resulting model will be bigger and slower. We go with 500 trees here (the default is 100), but ideally we should run a hyper-parameter search to get the best results. .. GENERATED FROM PYTHON SOURCE LINES 1528-1534 .. code-block:: Python gb_final = StackingRegressor( estimators=knn_regressor, final_estimator=GradientBoostingRegressor(n_estimators=500, random_state=0), passthrough=True, ) .. GENERATED FROM PYTHON SOURCE LINES 1535-1537 .. code-block:: Python gb_final.fit(df_train, y_train) .. raw:: html 
StackingRegressor(estimators=[('knn@5',
                                   Pipeline(steps=[('select_cols',
                                                    ColumnTransformer(transformers=[('long_and_lat',
                                                                                     'passthrough',
                                                                                     ['Longitude',
                                                                                      'Latitude'])])),
                                                   ('knn',
                                                    KNeighborsRegressor())]))],
                      final_estimator=GradientBoostingRegressor(n_estimators=500,
                                                                random_state=0),
                      passthrough=True)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.


.. GENERATED FROM PYTHON SOURCE LINES 1538-1540 .. code-block:: Python -get_scorer("neg_root_mean_squared_error")(gb_final, df_test, y_test) .. rst-class:: sphx-glr-script-out .. code-block:: none 40744.46986035733 .. GENERATED FROM PYTHON SOURCE LINES 1541-1544 The final RMSE is a bit better than we got earlier when using 100 trees, while the training time is still reasonably fast, so we can be happy with the outcome. .. GENERATED FROM PYTHON SOURCE LINES 1546-1548 Sharing the model ----------------- .. GENERATED FROM PYTHON SOURCE LINES 1550-1552 Saving the model artifact ~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 1554-1560 Now that we trained the final model, we should save it for later usage. We could use pickle (or joblib) to do this, and there’s nothing wrong with that. However, the resulting file is insecure because it can theoretically execute arbitrary code when loading. Thus, if we want to share this models with other people, they might be reluctant to unpickle the file if they don’t trust us completely. .. GENERATED FROM PYTHON SOURCE LINES 1562-1567 An alternative to the pickle format is the skops format. It is built with security in mind, therefore, other people can open your skops file without worries, even if they don’t trust us. More on the skops format can be found `here `__. .. GENERATED FROM PYTHON SOURCE LINES 1569-1573 For the purpose of this exercise, let’s use skops by calling ``skops.io.dump`` (``skops.io`` was imported as ``sio``) and store the model in a temporary directory, as shown below: .. GENERATED FROM PYTHON SOURCE LINES 1576-1579 .. code-block:: Python temp_dir = Path(mkdtemp()) file_name = temp_dir / "model.skops" .. GENERATED FROM PYTHON SOURCE LINES 1580-1582 .. code-block:: Python sio.dump(gb_final, file_name) .. GENERATED FROM PYTHON SOURCE LINES 1583-1585 Creating a model card ~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 1587-1589 When we want to share the model with others, it’s good practice to add a model card. That way, interested users can quickly learn what to expect. .. GENERATED FROM PYTHON SOURCE LINES 1591-1598 As a first approximation, a model card is a text document, often written in markdown, that contains sections talking about what kind of problem we’re dealing with, what kind of model to is used, what the intended purpose is, how to contact the authors, etc. It may also contain some metadata, which is targeted at machines and contains, say, tags that indicate what type of model or task being used. For now, we start without metadata. .. GENERATED FROM PYTHON SOURCE LINES 1600-1605 To help getting started, we can use the ``skops.card.Card`` class. It comes with a few default sections and provides some convenient methods for adding figures etc. The `documentation `__ goes into more details. .. GENERATED FROM PYTHON SOURCE LINES 1607-1611 For now, let’s start by creating a new model card and adding a few bits of information. We pass our final model as an argument to the ``Card`` class, which is used to create a table of hyper-parameters and a diagram of the model. .. GENERATED FROM PYTHON SOURCE LINES 1614-1617 .. code-block:: Python model_card = card.Card(model=gb_final) model_card .. rst-class:: sphx-glr-script-out .. code-block:: none Card( model=StackingRegressor(estimators=[('kn... random_state=0), passthrough=True), Model description/Training Procedure/Hyperparameters=TableSection(51x2), Model description/Training Procedure/Model Plot=