<jupyter_start><jupyter_text>Predictive Models  "Housekeeping"<jupyter_code>import json

import numpy as np
import pandas as pd

from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import Lasso, LinearRegression, Ridge
from sklearn.metrics import (
    make_scorer,
    mean_absolute_error,
    mean_squared_error,
    r2_score,
)
from sklearn.model_selection import GridSearchCV, KFold
from sklearn.svm import SVR
from sklearn.utils import shuffle
from tqdm import tqdm as progress_bar

from utils import (
    CONTINUOUS_VARIABLES,
    DISCRETE_VARIABLES,
    NOMINAL_VARIABLES,
    ORDINAL_VARIABLES,
    TARGET_VARIABLES,
    bias_score,
    encode_ordinals,
    load_clean_data,
    max_deviation,
)
random_state = np.random.RandomState(42)
pd.set_option("display.max_columns", 250)<jupyter_output><empty_output><jupyter_text>Load the Data  1. Original DataThe DataFrame `df1` holds the cleaned data from notebook 1 with the all the nominal and ordinal features automatically translated to factor variables and ordered integer values.<jupyter_code>df1 = load_clean_data("data/data_clean.csv")<jupyter_output><empty_output><jupyter_text>This cell basically **replaces** all the manual work that went into generating new and identifying "interesting" features in notebooks 2 and 3.<jupyter_code>df1 = pd.concat([
    df1[CONTINUOUS_VARIABLES + DISCRETE_VARIABLES + ORDINAL_VARIABLES + TARGET_VARIABLES],
    pd.get_dummies(df1[NOMINAL_VARIABLES], dtype=int),
], axis=1)
# Re-order the columns for convenience.
df1 = df1[sorted(set(df1.columns) - set(TARGET_VARIABLES)) + TARGET_VARIABLES]
df1 = encode_ordinals(df1)
df1 = shuffle(df1, random_state=random_state)
df1.info()
df1.head()<jupyter_output><empty_output><jupyter_text>Obtain the raw numpy arrays:<jupyter_code>X1 = df1.drop(columns="SalePrice").values
y1 = df1["SalePrice"].values<jupyter_output><empty_output><jupyter_text>2. Improved DataThe DataFrame `df2` holds the data as manually processed in notebooks 2 and 3.New features like the *years_since_\** columns were generated or derived from other variables like *has 2nd Flr* (from the continuous *2nd Flr SF*). Further, factor variables were created taking into account patterns in the visualizations. For example, *Bldg Type*'s (from `df1`) five categories were condensed into just three. In summary, `df2` has less than half as many dimensions as `df1` to accomodate for a potential curse of dimensionality.<jupyter_code>df2 = load_clean_data("data/data_clean_with_transformations_and_factors.csv")
df2 = encode_ordinals(df2)
df2 = shuffle(df2, random_state=random_state)
df2.info()
df2.head()<jupyter_output><empty_output><jupyter_text>Obtain the raw numpy arrays:<jupyter_code>X2 = df2.drop(columns=["SalePrice", "SalePrice (box-cox-0)"]).values
y2 = df2["SalePrice"].values
y2l = df2["SalePrice (box-cox-0)"].values<jupyter_output><empty_output><jupyter_text>3. Manual Feature Pre-SelectionAlso, notebook 2 collects variables that correlate either weakly ($0.33  0.66$) with the *SalePrice* (or the logarithm thereof). These variables serve as a "naive" feature pre-selection.<jupyter_code>with open("data/correlated_variables.json", "r") as file:
    _ = json.loads(file.read())
    weakly_correlated = _["weakly_correlated"]
    strongly_correlated = _["strongly_correlated"]
pre_selection = sorted(set(weakly_correlated + strongly_correlated) & set(df2.columns))<jupyter_output><empty_output><jupyter_text>The `df3` DataFrame is just a subset of `df2` (71 columns).<jupyter_code>df3 = df2[pre_selection + TARGET_VARIABLES]
df3.info(verbose=False)
df3.head()<jupyter_output><empty_output><jupyter_text>Obtain the raw numpy arrays:<jupyter_code>X3 = df3.drop(columns=["SalePrice", "SalePrice (box-cox-0)"]).values
y3 = df3["SalePrice"].values
y3l = df3["SalePrice (box-cox-0)"].values<jupyter_output><empty_output><jupyter_text>Re-usable ComponentsDefine a function to run k-fold cross validation to obtain unbiased estimators for the following scores / errors:- Bias- Mean Absolute Error- Maximum Deviation (just to see the worst case prediction of a model)- R2 (coefficient of determination)- Root Mean Squared Error (default for comparison)<jupyter_code>def cross_validation(X, y, *, model, k=10, log=False, desc=None):
    """Perform a k-fold cross validation."""
    bias, mae, max_dev, r2, rmse = [], [], [], [], []
    # Iterate over the k folds.
    for train, test in progress_bar(KFold(n_splits=k).split(X), desc=desc, total=k):
        model.fit(X[train], y[train])
        y_pred = model.predict(X[test])
        # If the sales price is provided on a log scale,
        # take the exponent first so that scores and
        # errors are comparable to the non-logged counterparts.
        if log:
            y_true, y_pred = np.exp(y[test]), np.exp(y_pred) 
        else:
            y_true, y_pred = y[test], y_pred
        # Collect the scores/errors for each fold.
        bias.append(bias_score(y_true, y_pred))
        mae.append(mean_absolute_error(y_true, y_pred))
        max_dev.append(max_deviation(y_true, y_pred))
        r2.append(r2_score(y_true, y_pred))
        rmse.append(mean_squared_error(y_true, y_pred))
    # Round for convenience.
    return {
        "bias": np.round(np.mean(bias)),
        "mae": np.round(np.mean(mae)),
        "max_dev": np.round(np.mean(max_dev)),
        "r2": np.round(np.mean(r2), 3),
        "rmse": np.round(np.sqrt(np.mean(rmse))),
    }<jupyter_output><empty_output><jupyter_text>Use one dictionary to store all the results in a systematic way.<jupyter_code>results = {}<jupyter_output><empty_output><jupyter_text>Linear RegressionA plain OLS regression model serves as the base case for benchmarking.<jupyter_code>lm = LinearRegression()<jupyter_output><empty_output><jupyter_text>1. Original DataGiven the unprocessed data, the linear model is not able to make a good fit at all.<jupyter_code>results[('lm','o')] = cross_validation(X1, y1, model=lm)
results[('lm','o')]<jupyter_output>100%|██████████| 10/10 [00:01<00:00,  9.86it/s]<jupyter_text>2. Improved Data  a) Normal Scale<jupyter_code>results[('lm','i')] = cross_validation(X2, y2, model=lm)
results[('lm','i')]<jupyter_output>100%|██████████| 10/10 [00:00<00:00, 22.99it/s]<jupyter_text>b) Log Scale<jupyter_code>results[('lm','il')] = cross_validation(X2, y2l, model=lm, log=True)
results[('lm','il')]<jupyter_output>100%|██████████| 10/10 [00:00<00:00, 19.21it/s]<jupyter_text>3. Improved Data with pre-selected Features  a) Normal Scale<jupyter_code>results[('lm','p')] = cross_validation(X3, y3, model=lm)
results[('lm','p')]<jupyter_output>100%|██████████| 10/10 [00:00<00:00, 278.33it/s]<jupyter_text>b) Log Scale<jupyter_code>results[('lm','pl')] = cross_validation(X3, y3l, model=lm, log=True)
results[('lm','pl')]<jupyter_output>100%|██████████| 10/10 [00:00<00:00, 271.17it/s]<jupyter_text>LASSO<jupyter_code>tol = 0.1
grid_search = GridSearchCV(
    estimator=Lasso(tol=tol, selection="random", random_state=random_state),
    param_grid={"alpha": [2 ** x for x in range(-8, 4)] + list(range(12, 65, 4))},
    cv=KFold(n_splits=4),
    n_jobs=-1,
)<jupyter_output><empty_output><jupyter_text>1. Original Data<jupyter_code>grid_search.fit(X1, y1)
alpha = grid_search.best_params_["alpha"]
alpha
results[("lasso", "o")] = cross_validation(X1, y1, model=Lasso(alpha=alpha, tol=tol))
results[("lasso", "o")]<jupyter_output>100%|██████████| 10/10 [00:00<00:00, 38.22it/s]<jupyter_text>2. Improved Data  a) Normal Scale<jupyter_code>grid_search.fit(X2, y2)
alpha = grid_search.best_params_["alpha"]
alpha
results[("lasso", "i")] = cross_validation(X2, y2, model=Lasso(alpha=alpha, tol=tol))
results[("lasso", "i")]<jupyter_output>100%|██████████| 10/10 [00:00<00:00, 54.59it/s]<jupyter_text>b) Log Scale<jupyter_code>grid_search.fit(X2, y2l)
alpha = grid_search.best_params_["alpha"]
alpha
results[("lasso", "il")] = cross_validation(X2, y2l, model=Lasso(alpha=alpha, tol=tol), log=True)
results[("lasso", "il")]<jupyter_output>100%|██████████| 10/10 [00:00<00:00, 46.34it/s]<jupyter_text>3. Improved Data with pre-selected Features  a) Normal Scale<jupyter_code>grid_search.fit(X3, y3)
alpha = grid_search.best_params_["alpha"]
alpha
results[("lasso", "p")] = cross_validation(X3, y3, model=Lasso(alpha=alpha, tol=tol))
results[("lasso", "p")]<jupyter_output>100%|██████████| 10/10 [00:00<00:00, 328.75it/s]<jupyter_text>b) Log Scale<jupyter_code>grid_search.fit(X3, y3l)
alpha = grid_search.best_params_["alpha"]
alpha
results[("lasso", "pl")] = cross_validation(X3, y3l, model=Lasso(alpha=alpha, tol=tol), log=True)
results[("lasso", "pl")]<jupyter_output>100%|██████████| 10/10 [00:00<00:00, 299.57it/s]<jupyter_text>Ridge Regression<jupyter_code>grid_search = GridSearchCV(
    estimator=Ridge(),
    param_grid={"alpha": [2 ** x for x in range(-8, 4)] + list(range(12, 65, 4))},
    cv=KFold(n_splits=4),
    n_jobs=-1,
)<jupyter_output><empty_output><jupyter_text>1. Original Data<jupyter_code>grid_search.fit(X1, y1)
alpha = grid_search.best_params_["alpha"]
alpha
results[("ridge", "o")] = cross_validation(X1, y1, model=Ridge(alpha=alpha))
results[("ridge", "o")]<jupyter_output>100%|██████████| 10/10 [00:00<00:00, 41.08it/s]<jupyter_text>2. Improved Data  a) Normal Scale<jupyter_code>grid_search.fit(X2, y2)
alpha = grid_search.best_params_["alpha"]
alpha
results[("ridge", "i")] = cross_validation(X2, y2, model=Ridge(alpha=alpha))
results[("ridge", "i")]<jupyter_output>100%|██████████| 10/10 [00:00<00:00, 69.82it/s]<jupyter_text>b) Log Scale<jupyter_code>grid_search.fit(X2, y2l)
alpha = grid_search.best_params_["alpha"]
alpha
results[("ridge", "il")] = cross_validation(X2, y2l, model=Ridge(alpha=alpha), log=True)
results[("ridge", "il")]<jupyter_output>100%|██████████| 10/10 [00:00<00:00, 74.39it/s]<jupyter_text>3. Improved Data with pre-selected Features  a) Normal Scale<jupyter_code>grid_search.fit(X3, y3)
alpha = grid_search.best_params_["alpha"]
alpha
results[("ridge", "p")] = cross_validation(X3, y3, model=Ridge(alpha=alpha))
results[("ridge", "p")]<jupyter_output>100%|██████████| 10/10 [00:00<00:00, 64.44it/s]<jupyter_text>b) Log Scale<jupyter_code>grid_search.fit(X3, y3l)
alpha = grid_search.best_params_["alpha"]
alpha
results[("ridge", "pl")] = cross_validation(X3, y3l, model=Ridge(alpha=alpha), log=True)
results[("ridge", "pl")]<jupyter_output>100%|██████████| 10/10 [00:00<00:00, 60.47it/s]<jupyter_text>Random Forest<jupyter_code>rf = RandomForestRegressor(
    n_estimators=500,
    n_jobs=-1, random_state=random_state
)<jupyter_output><empty_output><jupyter_text>1. Original Data<jupyter_code>results[("rf", "o")] = cross_validation(X1, y1, model=rf)
results[("rf", "o")]<jupyter_output>100%|██████████| 10/10 [00:19<00:00,  1.96s/it]<jupyter_text>2. Improved Data  a) Normal Scale<jupyter_code>results[("rf", "i")] = cross_validation(X2, y2, model=rf)
results[("rf", "i")]<jupyter_output>100%|██████████| 10/10 [00:16<00:00,  1.69s/it]<jupyter_text>b) Log Scale<jupyter_code>results[("rf", "il")] = cross_validation(X2, y2l, model=rf, log=True)
results[("rf", "il")]<jupyter_output>100%|██████████| 10/10 [00:19<00:00,  1.97s/it]<jupyter_text>3. Improved Data with pre-selected Features  a) Normal Scale<jupyter_code>results[("rf", "p")] = cross_validation(X3, y3, model=rf)
results[("rf", "p")]<jupyter_output>100%|██████████| 10/10 [00:12<00:00,  1.29s/it]<jupyter_text>b) Log Scale<jupyter_code>results[("rf", "pl")] = cross_validation(X3, y3l, model=rf, log=True)
results[("rf", "pl")]<jupyter_output>100%|██████████| 10/10 [00:11<00:00,  1.13s/it]<jupyter_text>Analysis of ResultsThis notebook did not focus on hyper-parameter optimization. Therefore, the predictions by Lasso, Ridge, and the Random Forest can potentially be improved by fine-graining the grid search even more.In general, the manually "improved" data clearly outperform the data that were only cleaned with the minimum effort. Also, the result suggests to allow the model to select its features. The manually pre-selected features perform well but not as good as the full feature set.<jupyter_code>def scores_by_source(source, score="rmse", *, ascending=True):
    rv = [
        (model, scores[score])
        for (model, data_source), scores in results.items()
        if data_source == source
    ]
    return sorted(rv, key=lambda x: x[1], reverse=(not ascending))<jupyter_output><empty_output><jupyter_text>Root Mean Squared Error<jupyter_code>scores_by_source("o", "rmse")
scores_by_source("i", "rmse")
scores_by_source("il", "rmse")
scores_by_source("p", "rmse")
scores_by_source("pl", "rmse")<jupyter_output><empty_output><jupyter_text>R2<jupyter_code>scores_by_source("o", "r2", ascending=False)
scores_by_source("i", "r2", ascending=False)
scores_by_source("il", "r2", ascending=False)
scores_by_source("p", "r2", ascending=False)
scores_by_source("pl", "r2", ascending=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16474413563520943, 'fraction_numerical': 0.04412789347371681, 'mean_word_length': 6.781927710843373, 'pattern_counts': {'":': 7, '<': 153, '<?xml version=': 0, '>': 133, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '38', 'dir': 'jupyter-structured-clean-dedup', 'id': '532365', 'license': 'MIT', 'n_tokens_mistral': 5365, 'n_tokens_neox': 4914, 'n_words': 976, 'path': '03_predictive_models.ipynb', 'repo_name': 'webartifex/ames-housing', 'size': '12917'}