---
jupytext:
  formats: content_md///md:myst,content///ipynb
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.16.4
kernelspec:
  display_name: Python 3 (ipykernel)
  language: python
  name: python3
---

(supervised:linear-regression:example1)=
# Worked Example: House Prices I

To test techniques for supervised learning discussed so far we train a model for predicting house prices in Germany. Inputs are properties of a house and of the plot of land it has been built on. Output is the selling price.

Training data exists in form of advertisements on specialized websites for finding a buyer for a house. In principle we could scrape data from such a website, but usually it's not allowed by the website operator and we would have to write lots of code. [Erdogan Seref](https://www.kaggle.com/scriptsultan) (unreachable since 2023) already did this job at [www.immobilienscout24.de](https://www.immobilienscout24.de) and published the data set at [www.kaggle.com](https://www.kaggle.com/scriptsultan/german-house-prices) (unreachable since 2023) under a [Attribution-NonCommercial-ShareAlike 4.0 International License](https://creativecommons.org/licenses/by-nc-sa/4.0).

```{code-cell} ipython3
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

import sklearn.linear_model as linear_model
import sklearn.metrics as metrics
import sklearn.model_selection as model_selection
import sklearn.preprocessing as preprocessing
import sklearn.pipeline as pipeline

data_path = 'german_housing.csv'
```

+++ {"user_expressions": []}

## The Data Set

At first we load the data set and try to get an overview of features and data quality.

```{code-cell} ipython3
data = pd.read_csv(data_path)
```

+++ {"user_expressions": []}

If a data frame has many columns Pandas by default does not show all columns. But we want to see all. Thus, we have to adjust the number of columns shown by [setting corresponding option](https://pandas.pydata.org/pandas-docs/stable/user_guide/options.html#frequently-used-options) to `None` (that is, unlimited).

```{code-cell} ipython3
pd.set_option('display.max_columns', None)
data.head(10)
```

```{code-cell} ipython3
data.info()
```

+++ {"user_expressions": []}

We should drop irrelevant columns and adjust data types.

* `Unnamed: 0`: Seems to be an integer index. We don't need it, so drop it.
* `Price`: This is our target variable.
* `Type`: An important column, because house prices are likely to depend on the type of house. We should convert this to categorical type.
* `Living_space` and `Lot`: Important features, keep them.
* `Usable_area`: Likely to have influence on the selling price, but available only for half the samples. If we want to use this for regression, we would have to drop half the training samples. Alternatively we could impute values, but it's very hard to guess usable area from other features. We should drop the column.
* `Free_of_Relation`: Not related to the selling price. Drop it.
* `Rooms`, `Bedrooms`, `Bathrooms`: Should have influence on prices, but not available for all samples. For the moment we keep all three columns. Later we should have a look on correlations between the three columns and possibly only keep the first one, which is available for all samples.
* `Floors`: Important feature, keep it.
* `Year_built`: Important feature, keep it.
* `Furnishing_quality`: Important, convert to categorical and keep.
* `Year_renovated`: Important, but half the data is missing. There is good chance that missing values indicate that there the house has not been renovated until today. Thus, a reasonable fill value is the year of construction.
* `Condition`: Important, convert to categorical and keep.
* `Heating` and `Energy_source`: Could be important, convert to categorical and keep.
* `Energy_certificate`, `Energy_certificate_type`, `Energy_consumption`: The first contains more or less only the value `'available'` (since energy certificates are required by law). The second is irrelevant and the third is missing for most samples. Drop them all.
* `Energy_efficiency_class`: Likely to have influence on the selling price, although classification procedure is very unreliable in practice. Keep and convert to categorical.
* `State`, `City`, `Place`: Geolocation surely influences selling prices. But it's hard to use location data for regression. For the moment we keep these columns.
* `Garages`: Could be important, keep.
* `Garagetype`: If we keep `Garages` then we also have to keep this column. Convert to categorical and rename to `Garage_type` to fit naming convention used for the other columns.

```{code-cell} ipython3
data = data.drop(columns=['Unnamed: 0', 'Usable_area', 'Free_of_Relation',
                          'Energy_certificate', 'Energy_certificate_type', 'Energy_consumption'])

data['Type'] = data['Type'].astype('category')
data['Furnishing_quality'] = data['Furnishing_quality'].astype('category')
data['Condition'] = data['Condition'].astype('category')
data['Heating'] = data['Heating'].astype('category')
data['Energy_source'] = data['Energy_source'].astype('category')
data['Energy_efficiency_class'] = data['Energy_efficiency_class'].astype('category')
data['Garagetype'] = data['Garagetype'].astype('category')

data = data.rename(columns={'Garagetype': 'Garage_type'})

nan_mask = data['Year_renovated'].isna()
data.loc[nan_mask, 'Year_renovated'] = data.loc[nan_mask, 'Year_built']
```

+++ {"user_expressions": []}

Categorical columns `Furnishing_quality`, `Condition` and `Energy_efficiency_class` should have a natural ordering, which should be represented by the data type.

```{code-cell} ipython3
print(data['Furnishing_quality'].cat.categories)
print(data['Condition'].cat.categories)
print(data['Energy_efficiency_class'].cat.categories)
```

+++ {"user_expressions": []}

We should rename same categories and sort them as good as possible.

```{code-cell} ipython3
data['Furnishing_quality'] = data['Furnishing_quality'] \
    .cat.rename_categories({'luxus': 'luxury'}) \
    .cat.reorder_categories(['basic', 'normal', 'refined', 'luxury'])

data['Condition'] = data['Condition'].cat.reorder_categories([
    'first occupation',
     'first occupation after refurbishment',
     'as new',
     'maintained',
     'renovated',
     'modernized',
     'refurbished',
     'by arrangement',
     'fixer-upper',
     'dilapidated'
])

data['Energy_efficiency_class'] = data['Energy_efficiency_class'] \
    .cat.rename_categories({
        ' A ': 'A',
        ' A+ ': 'A+',
        ' B ': 'B',
        ' C ': 'C',
        ' D ': 'D',
        ' E ': 'E',
        ' F ': 'F',
        ' G ': 'G',
        ' H ': 'H'
    }) \
    .cat.reorder_categories(['A+', 'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H'])
```

+++ {"user_expressions": []}

Now let's see how many complete samples we have.

```{code-cell} ipython3
len(data.dropna())
```

+++ {"user_expressions": []}

That's very few. So we should drop some columns with many missing values.

```{code-cell} ipython3
data.info()
```

+++ {"user_expressions": []}

`Energy_efficiency_class` is relatively unreliable and not too important for selling prices.

```{code-cell} ipython3
len(data.drop(columns=['Energy_efficiency_class']).dropna())
```

+++ {"user_expressions": []}

Better, but not good. The `Bedrooms` column has many missing values, too, and it's likely to be correlated to `Rooms`. So let's look at correlations between `Rooms`, `Bedrooms`, `Bathrooms`, `Floors`.

```{code-cell} ipython3
sns.pairplot(data[['Rooms', 'Bedrooms', 'Bathrooms', 'Floors']], plot_kws={"s": 5})
plt.show()
```

+++ {"user_expressions": []}

`Floors` is not correlated to the other columns, so keep it. `Bedrooms` show correlation to `Rooms` and `Bathrooms`, so drop `Bedrooms`. `Bathroom` shows some correlation to `Rooms`. Wether to drop `Bathrooms` should be decided by the increase in sample counts.

```{code-cell} ipython3
len(data.drop(columns=['Energy_efficiency_class', 'Bedrooms']).dropna())
```

```{code-cell} ipython3
len(data.drop(columns=['Energy_efficiency_class', 'Bedrooms', 'Bathrooms']).dropna())
```

+++ {"user_expressions": []}

We should keep `Bathrooms`, because dropping it only yields 300 more samples while neglecting possibly important information. Note that the number of bath rooms can be regarded as a measure for overall furnishing quality. Thus, there should be some correlation to `Furnishing_quality`.

```{code-cell} ipython3
data.groupby('Furnishing_quality')['Bathrooms'].mean()
```

+++ {"user_expressions": []}

In addition, judging about furnishing quality of a house is highly subjective. Thus, we should drop the column to get more samples without missing data.

```{code-cell} ipython3
len(data.drop(columns=['Energy_efficiency_class', 'Bedrooms', 'Furnishing_quality']).dropna())
```

+++ {"user_expressions": []}

The `Energy_source` is another candidate for dropping, because it has more than 1000 missing values and its influence on selling prices should be rather low.

```{code-cell} ipython3
for cat in data['Energy_source'].cat.categories:
    print(cat)
```

+++ {"user_expressions": []}

Values are very diverse and hard to preprocess for regression. We would have to convert the column to several boolean columns. In addition, some grouping would be necessary (`Holz` is a subcategory of `Bioenergie` and so on).

```{code-cell} ipython3
len(data.drop(columns=['Energy_efficiency_class', 'Bedrooms', 'Furnishing_quality', 'Energy_source']).dropna())
```

+++ {"user_expressions": []}

Now we have almost 5000 complete samples. Should be a good compromise between completeness and level of detail.

```{code-cell} ipython3
data = data.drop(columns=['Energy_efficiency_class', 'Bedrooms', 'Furnishing_quality', 'Energy_source'])
data = data.dropna()
```

+++ {"user_expressions": []}

## Outliers and Further Preprocessing

Now that we have a cleaned data set we should remove outliers. The simplest method of detecting outliers is to look at the ranges of all feature. With `describe` we get a first overview for numerical features.

```{code-cell} ipython3
data.describe()
```

+++ {"user_expressions": []}

### `Price` Column

```{code-cell} ipython3
sns.histplot(data['Price'])
plt.show()
```

+++ {"user_expressions": []}

There are only very few high prices and price distribution concentrates on low prices. If the target variable has wide range, but most samples concentrate on a small portion of the range, then 'learning' the target is much more difficult than for more uniformly distributed data.

A common trick is to use nonlinear scaling. Especially for market prices it is known from experience that they follow a [log-normal distribution](https://en.wikipedia.org/wiki/Log-normal_distribution), that is, after applying the logarithm we see a normal distribution. Before applying the logarithm we should drop samples with zeros in the `Price` column to avoid undefined results. A price of zero indicates that the seller did not provide a price in the advertisement. Thus, dropping such sample even is a good idea if we wouldn't want to apply the logarithm.

```{code-cell} ipython3
data = data.loc[data['Price'] > 0, :]
sns.histplot(np.log(data['Price'].to_numpy()))
plt.show()
```

+++ {"user_expressions": []}

There seem to be same very small values.

```{code-cell} ipython3
data.loc[data['Price'] <= np.exp(7), :]
```

+++ {"user_expressions": []}

Those samples should be dropped because house prices below $\mathrm{e}^7\approx 1000$ EUR are very uncommon.

```{code-cell} ipython3
data['Price'] = np.log(data['Price'].to_numpy())
data = data.loc[data['Price'] > 7, :]
```

+++ {"user_expressions": []}

### `Living_space` Column

```{code-cell} ipython3
sns.histplot(data['Living_space'])
plt.show()
```

+++ {"user_expressions": []}

Same here as for `Price`.

```{code-cell} ipython3
data = data.loc[data['Living_space'] > 0, :]
sns.histplot(np.log(data['Living_space'].to_numpy()))
plt.show()
```

```{code-cell} ipython3
data['Living_space'] = np.log(data['Living_space'].to_numpy())
```

+++ {"user_expressions": []}

### `Lot` Column

```{code-cell} ipython3
sns.histplot(data['Lot'])
plt.show()
```

+++ {"user_expressions": []}

Same here as for `Price` again.

```{code-cell} ipython3
data = data.loc[data['Lot'] > 0, :]
sns.histplot(np.log(data['Lot'].to_numpy()))
plt.show()
```

```{code-cell} ipython3
data.loc[data['Lot'] <= np.exp(2), :]
```

+++ {"user_expressions": []}

Lot size below $\mathrm{e}^2<8$ m² is very unlikely.

```{code-cell} ipython3
data['Lot'] = np.log(data['Lot'].to_numpy())
data = data.loc[data['Lot'] > 2, :]
```

+++ {"user_expressions": []}

### `Rooms` Column

```{code-cell} ipython3
sns.histplot(data['Rooms'])
plt.show()
```

```{code-cell} ipython3
data.loc[data['Rooms'] >= 30, :]
```

+++ {"user_expressions": []}

There are only very few sample with high number of rooms. There is no chance to get good predictions from those few samples.

```{code-cell} ipython3
data = data.loc[data['Rooms'] < 30, :]
```

+++ {"user_expressions": []}

### `Bathrooms` Column

```{code-cell} ipython3
sns.histplot(data['Bathrooms'])
plt.show()
```

```{code-cell} ipython3
data.loc[data['Bathrooms'] >= 15, :]
```

```{code-cell} ipython3
data = data.loc[data['Bathrooms'] < 15, :]
```

+++ {"user_expressions": []}

### `Floors` Column

```{code-cell} ipython3
sns.histplot(data['Floors'])
plt.show()
```

+++ {"user_expressions": []}

Nothing to do here.

### `Year_built` Column

```{code-cell} ipython3
sns.histplot(data['Year_built'])
plt.show()
```

```{code-cell} ipython3
data.loc[data['Year_built'] <= 1500, :]
```

```{code-cell} ipython3
data = data.loc[data['Year_built'] > 1500, :]
```

+++ {"user_expressions": []}

Values above 2020 obviously are wrong (data set is from 2020).

```{code-cell} ipython3
data.loc[data['Year_built'] > 2020, :]
```

```{code-cell} ipython3
data = data.loc[data['Year_built'] <= 2020, :]
```

+++ {"user_expressions": []}

To get a better distribution of the samples over the range, we again apply a logarithmic transform.

```{code-cell} ipython3
data.loc[:, 'Year_built'] = np.log(2021 - data['Year_built'].to_numpy())
sns.histplot(data['Year_built'])
plt.show()
```

+++ {"user_expressions": []}

### `Year_renovated` Column

```{code-cell} ipython3
sns.histplot(data['Year_renovated'])
plt.show()
```

+++ {"user_expressions": []}

There seem to be renovations before 1900, which seems somewhat strange. But remember that we filled missing values with values from `Year_built`. Values above 2020 obviously are wrong.

```{code-cell} ipython3
data.loc[data['Year_renovated'] > 2020, :]
```

```{code-cell} ipython3
data = data.loc[data['Year_renovated'] <= 2020, :]
```

```{code-cell} ipython3
data.loc[:, 'Year_renovated'] = np.log(2021 - data['Year_renovated'].to_numpy())
sns.histplot(data['Year_renovated'])
plt.show()
```

+++ {"user_expressions": []}

### `Garages` Column

```{code-cell} ipython3
sns.histplot(data['Garages'])
plt.show()
```

```{code-cell} ipython3
data.loc[data['Garages'] >= 20, :]
```

```{code-cell} ipython3
data = data.loc[data['Garages'] < 20, :]
```

+++ {"user_expressions": []}

### `Type` Column

```{code-cell} ipython3
data['Type'].value_counts()
```

```{code-cell} ipython3
data = data.loc[data['Type'] != 'Castle', :]
data['Type'] = data['Type'].cat.remove_categories('Castle')
```

+++ {"user_expressions": []}

### `Condition` Column

```{code-cell} ipython3
data['Condition'].value_counts()
```

+++ {"user_expressions": []}

We should remove `'as new'` and `'by arrangement'` because only few samples use these categories and both are somewhat dubious.

```{code-cell} ipython3
data = data.loc[~data['Condition'].isin(['as new', 'by arrangement']), :]
data['Condition'] = data['Condition'].cat.remove_categories(['as new', 'by arrangement'])
```

+++ {"user_expressions": []}

### `Heating` Column

```{code-cell} ipython3
data['Heating'].value_counts()
```

+++ {"user_expressions": []}

Something is wrong here! More than every second house sold in 2020 has stove heating? And what about `'floor heating'`? Is it gas powered or oil powered or what else? What's the difference between `'floor heating'` and `'underfloor heating'`. It's better to drop this column.

```{code-cell} ipython3
data = data.drop(columns=['Heating'])
```

+++ {"user_expressions": []}

### `Garage_type` Column

```{code-cell} ipython3
data['Garage_type'].value_counts()
```

+++ {"user_expressions": []}

There are many similar categories. We should join some.

```{code-cell} ipython3
data.loc[data['Garage_type'] == 'Car park lot', 'Garage_type'] = 'Outside parking lot'
data.loc[data['Garage_type'] == 'Duplex lot', 'Garage_type'] = 'Outside parking lot'
data.loc[data['Garage_type'] == 'Parking lot', 'Garage_type'] = 'Outside parking lot'
data['Garage_type'] = data['Garage_type'].cat.remove_categories(['Car park lot', 'Duplex lot', 'Parking lot'])

data['Garage_type'].value_counts()
```

+++ {"user_expressions": []}

## Save Cleaned Data

We save cleaned data for future use.

```{code-cell} ipython3
data.to_csv(data_path.replace('.csv', '_preprocessed.csv'))
```

+++ {"user_expressions": []}

## Linear Regression

Now data is almost ready for training a model. It remains to convert categorical data to numerical data. `Condition` is ordered and numeric representation is accessible with [`Series.cat.codes`](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.Series.cat.codes.html). Columns `Type` and `Garage_type` should be one-hot encoded.

```{code-cell} ipython3
data['Condition_codes'] = data['Condition'].cat.codes
data = pd.get_dummies(data, columns=['Type', 'Garage_type'], drop_first=True)
```

```{code-cell} ipython3
data.head()
```

+++ {"user_expressions": []}

We drop columns not used for regression and convert the data frame to NumPy arrays suitable for Scikit-Learn.

```{code-cell} ipython3
y = data['Price'].to_numpy()
X = data.drop(columns=['Price', 'Condition', 'State', 'City', 'Place']).to_numpy()

print(X.shape, y.shape)
```

+++ {"user_expressions": []}

We have relatively few data. Thus, test set should be small to have more training samples.

```{code-cell} ipython3
X_train, X_test, y_train, y_test = model_selection.train_test_split(X, y, test_size=0.2)
print(y_train.size, y_test.size)
```

+++ {"user_expressions": []}

We use polynomial regression with regularization.

```{code-cell} ipython3
steps = [('poly', preprocessing.PolynomialFeatures()),
         ('ridge', linear_model.Ridge())]

pipe = pipeline.Pipeline(steps)

param_grid = {'poly__degree': [1, 2, 3],
              'ridge__alpha': [0] + [2 ** k for k in range(5, 15)]}

gs = model_selection.GridSearchCV(pipe, param_grid=param_grid,
                                  scoring='neg_mean_squared_error', n_jobs=-1, cv=5)

gs.fit(X_train, y_train)
best_params = gs.best_params_
```

+++ {"user_expressions": []}

## Evaluating the Model

Now we use the test set to evaluate prediction quality of the model.

```{code-cell} ipython3
print(best_params)

pipe.set_params(**best_params)
pipe.fit(X_train, y_train)

y_test_pred = pipe.predict(X_test)
```

Root mean squared error between predicted and exact targets on its own does not tell much about fitting quality. We have to compare the value to standard deviation of the targets. Standard deviation is the root mean squared error of the exact targets and their mean. In other words, standard deviation tells us the prediction error if we would use constant predictions for all inputs. Obviously the constant should be the mean of the training (!) targets, but the mean of the training targets should be very close the mean of the test targets if test samples have been selected randomly.

```{code-cell} ipython3
rmse = metrics.root_mean_squared_error(y_test, y_test_pred)
sigma = np.std(y_test)
print('RMSE:', rmse)
print('standard deviation:', sigma)
print('ratio:', rmse / sigma)
```

+++ {"user_expressions": []}

We see that the model's prediction is better than constant prediction, but not so much.

We should have a closer look at the predictions. Since there is no natural ordering in the set of samples plotting `y_test` and `y_test_pred` with `plot` does not help much.

```{code-cell} ipython3
fig, ax = plt.subplots()
ax.plot(y_test, '-b', label='true targets')
ax.plot(y_test_pred, '-r', label='predictions')
ax.legend()
plt.show()
```

+++ {"user_expressions": []}

A better idea is to plot `y_test` versus `y_test_pred`. If true and predicted labels are close, then points should concentrate along the diagonal. Else they are far away from the diagonal.

```{code-cell} ipython3
fig, ax = plt.subplots()
ax.plot(y_test, y_test_pred, 'or', markersize=3)
ax.plot([9, 17], [9, 17], '-b')
ax.set_xlabel('true targets')
ax.set_ylabel('predictions')
ax.set_aspect('equal')
plt.show()
```

The red cloud shows some rotation compared to the blue line. Small target values get too high predictions and high target values get too low predictions. In other words, predictions tend to be too close to the targets' mean. Such behavior is typically observed if there are many similar samples with different targets in the training data. Then there is no clear functional dependence of the targets on the inputs and models tend to predict the mean targets.

To further investigate this issue we should look at the predictions on the training set. If we are right, then predictions on the training set should show similar behavior (predictions close to mean).

```{code-cell} ipython3
y_train_pred = pipe.predict(X_train)

fig, ax = plt.subplots()
ax.plot(y_train, y_train_pred, 'or', markersize=3)
ax.plot([9, 17], [9, 17], '-b')
ax.set_xlabel('true targets')
ax.set_ylabel('predictions')
ax.set_aspect('equal')
plt.show()
```

+++ {"user_expressions": []}

Again we see slight rotation. To summarize: our input data has too few details to explain the targets. There are simlar inputs with different targets leading to underestimation of high values and over estimation of low values. The only way out is gathering more data, either by dropping less columns or by getting relevant data from additional sources. We will come back to this issue soon.

## Predictions

When using our model for predicting house prices we have to keep in mind that we transformed some of the input data. All those transforms have to be applied to new inputs, too.

```{code-cell} ipython3
living_space = 80
lot = 3600
rooms = 5
bathrooms = 0
floors = 2
year_built = 1948
year_renovated = 1948
garages = 2
condition_codes = 7    # 0 = 'first occupation', 7 = 'dilapidated'
type_corner_house = 0
type_duplex = 0
type_farmhouse = 1
type_midterrace_house = 0
type_multiple_dwelling = 0
type_residential_property = 0
type_single_dwelling = 0
type_special_property = 0
type_villa = 0
garage_type_garage = 0
garage_type_outside_parking_lot = 1
garage_type_underground_parking_lot = 0

X = np.asarray([np.log(living_space), np.log(lot), rooms, bathrooms, floors,
                np.log(2021 - year_built), np.log(2021 - year_renovated),
                garages, condition_codes, type_corner_house, type_duplex, type_farmhouse,
                type_midterrace_house, type_multiple_dwelling, type_residential_property,
                type_single_dwelling, type_special_property, type_villa,
                garage_type_garage, garage_type_outside_parking_lot,
                garage_type_underground_parking_lot]).reshape(1, -1)

y = np.exp(pipe.predict(X))

print('predicted price: {:.0f} EUR'.format(y[0]))
```
