---
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:example2)=
# Worked Example: House Prices II

We try to improve prediction of house prices based on [Erdogan Seref's](https://www.kaggle.com/scriptsultan) (unreachable since 2023) German housing dataset from [www.immobilienscout24.de](https://www.immobilienscout24.de) published 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).

We load preprocessed data and adjust data types.

```{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_preprocessed.csv'
regions_path = 'regions.csv'
```

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

data['Type'] = data['Type'].astype('category')
data['Condition'] = data['Condition'].astype('category')
data['Garage_type'] = data['Garage_type'].astype('category')

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

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

## More Data

Results obtained from linear regression showed that input variables do not suffice to explain the targets. Thus, we should add more input variables. When preprocessing the data we dropped several columns. Keeping them could increase prediction quality slightly, but there were several good reasons to drop those columns. The main reason were lots of missing values in those columns.

A far better idea is to collect additional data. What features of a house influence the selling price? Of course its location! Up to now we did not use location information at all, but we have location information available. There are columns `State`, `City`, `Place`. But city names do not help. We need something like proximity to big cities or nice landscape. Adding a layer of abstraction we might ask for the demand for houses and the whealth of potential buyers. So we should head out for statistical information about local real estate markets and about economic power of different regions in Germany.

Everything we need is publicly available at [www.regionalstatistik.de](https://www.regionalstatistik.de) provided by *Statistische Ämter des Bundes und der Länder* under the license [Datenlizenz Deutschland – Namensnennung – Version 2.0](https://www.govdata.de/dl-de/by-2-0). Clicking here and there we find two interesting tables:
* [annual income per inhabitant](https://www.regionalstatistik.de/genesis//online?operation=table&code=AI-S-01&bypass=true&levelindex=0&levelid=1744088175420)
* [prices for construction ground](https://www.regionalstatistik.de/genesis//online?operation=table&code=61511-01-03-4&bypass=true&levelindex=0&levelid=1744088230007)

From those tables we may compile a table with 4 columns for region id, region name, income, ground prices.

The difficult part is matching region names in German housing data set with region names in the region table. Here is some code doing the job:

```{code-cell} ipython3
# remove rows with missing location information
data = data.dropna(subset=('State', 'City'))

# reindex to remove gaps in the index
data.index = pd.RangeIndex(0, len(data))
```

```{code-cell} ipython3
regions = pd.read_csv(regions_path)

regions.head(5)
```

```{code-cell} ipython3
data['city_short'] = data['City'].str.replace(' (Kreis)', '', regex=False)
data['region_idx'] = 0

for (idx, city_short) in enumerate(data['city_short']):
    find_results = regions['region'].str.find(str(city_short))
    if not (find_results > -1).any():
        #print(city_short)
        if data.loc[idx, 'State'] == 'Hamburg':
            find_results = regions['region'].str.find('Hamburg')
            data.loc[idx, 'region_idx'] = regions.index[find_results > -1][-1]
        elif data.loc[idx, 'State'] == 'Bremen':
            find_results = regions['region'].str.find('Bremen')
            data.loc[idx, 'region_idx'] = regions.index[find_results > -1][-1]
        elif data.loc[idx, 'State'] == 'Berlin':
            district = city_short.split('(')[-1][0:-1]
            if district == 'Weißensee':
                district = 'Pankow'
            if district == 'Prenzlauer Berg':
                district = 'Pankow'
            if district == 'Hohenschönhausen':
                district = 'Lichtenberg'
            if district == 'Wedding':
                district = 'Berlin-Mitte'
            find_results = regions['region'].str.find(district)
            #print('***', city_short, ':', district, '-->', regions['Region'][find_results > -1])
            data.loc[idx, 'region_idx'] = regions.index[find_results > -1][-1]
        elif city_short == 'Neuss (Rhein-Kreis)':
            find_results = regions['region'].str.find('Neuss')
            data.loc[idx, 'region_idx'] = regions.index[find_results > -1][-1]
        elif city_short == 'Sankt Wendel':
            find_results = regions['region'].str.find('Wendel')
            data.loc[idx, 'region_idx'] = regions.index[find_results > -1][-1]
        elif city_short == 'Stadtverband Saarbrücken':
            find_results = regions['region'].str.find('Saarbrücken, Regionalverband')
            data.loc[idx, 'region_idx'] = regions.index[find_results > -1][-1]
        elif city_short.split()[1] in ('in', 'im', 'an', 'am'):
            if city_short.split()[0] == 'Neustadt':
                find_results = regions['region'].str.find(city_short.split()[-1])
            else:
                find_results = regions['region'].str.find(city_short.split()[0])
            #print('***', city_short, '-->', regions.loc[find_results > -1, 'Region'])
            data.loc[idx, 'region_idx'] = regions.index[find_results > -1][-1]
        else:
            print('NOT FOUND:', city_short)
            
    else:
        # take last match (smallest region)
        data.loc[idx, 'region_idx'] = regions.index[find_results > -1][-1]
```

```{code-cell} ipython3
data['Region_id'] = regions.loc[data['region_idx'], 'id'].values
data['Income'] = regions.loc[data['region_idx'], 'income'].values
data['Land_prices'] = regions.loc[data['region_idx'], 'prices'].values
```

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

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

+++ {"user_expressions": []}

## Preprocessing New Columns

Now that we have two now columns (`Income` and `Land_prices`) we should look at their values.

+++ {"user_expressions": []}

### `Income` Column

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

+++ {"user_expressions": []}

Log-scaling is not mandadory here, but it brings values to a similar range like the other columns. Without log-scaling we would have a column with very large range, which may result in problems when using regularization (see [](supervised:linear-regression:regularization)). Alternatively we could standardize all columns before doing linear regression.

```{code-cell} ipython3
sns.histplot(np.log(data['Income'].to_numpy()))
plt.show()
```

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

+++ {"user_expressions": []}

### `Land_prices` Column

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

```{code-cell} ipython3
sns.histplot(np.log(data['Land_prices'].to_numpy()))
plt.show()
```

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

+++ {"user_expressions": []}

## Save Data

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

+++ {"user_expressions": []}

## Linear Regression

Now we do linear regression as before, but with two additional columns.

```{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
y = data['Price'].to_numpy()
X = data.drop(columns=['Price', 'Condition', 'State', 'City', 'Place', 'Region_id']).to_numpy()

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

```{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)
```

```{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": []}

## Evaluation

Now the interesting part. Do we see an increase in prediction quality?

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

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

y_test_pred = pipe.predict(X_test)
```

```{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)
```

```{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()
```

+++ {"user_expressions": []}

Looks much better!

## Feature Importance

With a trained model we may look at feature importances to see which features have high influence on the selling price.

```{code-cell} ipython3
import sklearn.inspection as inspection
```

```{code-cell} ipython3
result = inspection.permutation_importance(pipe, X, y, n_jobs=-1)
```

```{code-cell} ipython3
cols = data.drop(columns=['Price', 'Condition', 'State', 'City', 'Place', 'Region_id']).columns
imp = pd.Series(result.importances_mean, index=cols)
imp = imp.sort_values(ascending=False)
imp
```
