<jupyter_start><jupyter_text>Kaggle Machine Learning courseThis is a notebook for the Kaggle machine learning course available att www.kaggle.com. It uses `scikit-learn` version 0.19.2 and `pandas` version 0.23.3. The Iowa Housing data was downloaded from https://ww2.amstat.org/publications/jse/v19n3/decock/AmesHousing.txt. We will begin by loading and exploring the data a little bit to see how it is structured.<jupyter_code>import pandas as pd

iowa_data = pd.read_csv("./data/AmesHousing.txt", delimiter = '\t')

print("The columns are:")
print(iowa_data.columns)

print("Description:")
print(iowa_data.describe())
# Let's extract all the sales prices of each house.

selection = ["Order", "SalePrice"]
sale_prices = iowa_data[selection]
print(sale_prices.head())<jupyter_output>Order  SalePrice
0      1     215000
1      2     105000
2      3     172000
3      4     244000
4      5     189900<jupyter_text>It's time to make a first, albeit silly for reasons we will see later, model. For this one we will use a so called *decision tree*.<jupyter_code>from sklearn.tree import DecisionTreeRegressor

first_model = DecisionTreeRegressor(random_state = 0)
y = sale_prices.SalePrice
predictors = ["Lot Area", "Year Built", "1st Flr SF", "2nd Flr SF", "Full Bath", "Bedroom AbvGr", "TotRms AbvGrd"]
X = iowa_data[predictors]

first_model.fit(X, y)

# With the model trained let's have it "predict" some prices.

print("Actual prices:")
print(y.head())

print("Predicted prices:")
print(first_model.predict(X.head()))<jupyter_output>Actual prices:
0    215000
1    105000
2    172000
3    244000
4    189900
Name: SalePrice, dtype: int64
Predicted prices:
[215000. 105000. 172000. 244000. 189900.]<jupyter_text>To measure the accuracy of out first model we use the *mean absolute error*.<jupyter_code>from sklearn.metrics import mean_absolute_error as mae

predicted_prices = first_model.predict(X)
error = mae(y, predicted_prices)
print("The mean absolute error of our model on our sample is:")
print(error)<jupyter_output>The mean absolute error of our model on our sample is:
139.54379977246873<jupyter_text>This score doesn't really tell us much however since it was calculated on data that the model was trained on. In order to get a better understanding of the accuracy of our model we will split our data into two sets: a *training set* and a *validation set*. We will then train a new model on the training set and evaluate it using the validation set.<jupyter_code>from sklearn.model_selection import train_test_split as split

train_X, val_X, train_y, val_y = split(X, y, random_state = 0)
snd_model = DecisionTreeRegressor()
snd_model.fit(train_X, train_y)<jupyter_output><empty_output><jupyter_text>With our new model trained let's evaluate it on the validation set.<jupyter_code># With our new model trained let's evaluate it on the validation set.

predictions = snd_model.predict(val_X)

print("Actual prices:")
print(val_y.head())
print("Predicted price:")
print(predictions[:5])

print("Mean absolute error is:")
print(mae(predictions, val_y))<jupyter_output>Actual prices:
2216    220000
836     143000
2396    281000
1962    135000
305     102776
Name: SalePrice, dtype: int64
Predicted price:
[257500. 155000. 233230. 129000.  96500.]
Mean absolute error is:
27732.587994542973<jupyter_text>Now that we have the basics down of how to train and evaluate a model let's explore some ways of improving a model by tuning hyperparameters, in this case the number of leaf nodes. First we create a helper function to train a model with a specific number of leaf nodes on a training set and return its mean absolute error on a validation set. We will then loop through different values for the maximum number of leaf nodes and pick the best one we find.<jupyter_code>def get_mae(leaves, train_X, val_X, train_y, val_y):
    model = DecisionTreeRegressor(max_leaf_nodes = leaves, random_state = 0)
    model.fit(train_X, train_y)
    predictions = model.predict(val_X)
    return(mae(predictions, val_y))

for leaves in [50, 100, 250, 500, 1000, 2500, 5000]:
    print("Max leaf nodes: {arg1} \t\t Mean absolute error: {arg2}"
          .format(arg1 = leaves, arg2 = get_mae(leaves, train_X, val_X, train_y, val_y)))<jupyter_output>Max leaf nodes: 50 		 Mean absolute error: 24989.409446993875
Max leaf nodes: 100 		 Mean absolute error: 24401.58376237393
Max leaf nodes: 250 		 Mean absolute error: 24664.182995823816
Max leaf nodes: 500 		 Mean absolute error: 26177.366842277108
Max leaf nodes: 1000 		 Mean absolute error: 27688.295652880544
Max leaf nodes: 2500 		 Mean absolute error: 27549.8690313779
Max leaf nodes: 5000 		 Mean absolute error: 27549.8690313779<jupyter_text>In this case we can see that 100 maximum leaf nodes seems to be the optimal choice. With the decision tree model explored let's move on to a generalization of it called a *random decision forest*.<jupyter_code>from sklearn.ensemble import RandomForestRegressor

forest_model = RandomForestRegressor(random_state = 0)
forest_model.fit(train_X, train_y)
predictions = forest_model.predict(val_X)
print("The mean absolute error of the random forest model is:")
print(mae(predictions, val_y))<jupyter_output>The mean absolute error of the random forest model is:
21554.63075423894<jupyter_text>Advanced topicsThat covers the basics of creating and evaluating machine learning models. We will now move on to more advanced topics and begin by exploring the problem of handling missing values. First we need to find which columns actually contain some missing data. We will try three progressively more advanced methods: dropping columns with missing values, imputing missing values and finally imputing missing values and marking which fields were missing.<jupyter_code># Finding all columns with missing data. We will restrict ourselves to numeric data for convenience.

import numpy as np

numeric_data = iowa_data.select_dtypes(include = np.number)

columns_with_missing_data = [col for col in numeric_data.columns if iowa_data[col].isnull().any()]
print("The following columns are missing data:")
print(columns_with_missing_data)

# Extend predictors with some columns missing data and create new training and validation sets.

extra_columns = ["Lot Frontage", "Total Bsmt SF", "Garage Yr Blt", "Garage Area"]
predictors_with_missing = predictors + extra_columns
X = numeric_data[predictors_with_missing]
y = numeric_data.SalePrice
train_X, val_X, train_y, val_y = split(X, y, random_state = 0)
# Remove columns with missing data from training and validation sets.

reduced_train_X = train_X.drop(extra_columns, axis = 1)
reduced_val_X = val_X.drop(extra_columns, axis = 1)

# Creating and evaluating the model with dropped columns. This is the same as our previous model.

reduced_model = RandomForestRegressor(random_state = 0)
reduced_model.fit(reduced_train_X, train_y)
print("Mean absolute error of the reduced mdoel:")
print(mae(reduced_model.predict(reduced_val_X), val_y))
# Next we will impute values.

from sklearn.preprocessing import Imputer

imputer = Imputer(strategy = "mean", axis = 0)

impute_train_X = pd.DataFrame(imputer.fit_transform(train_X), index = train_X.index, columns = train_X.columns)
impute_val_X = pd.DataFrame(imputer.transform(val_X), index = val_X.index, columns = val_X.columns)
impute_model = RandomForestRegressor(random_state = 0)
impute_model.fit(impute_train_X, train_y)
print("Mean absolute error of the imputed model:")
print(mae(impute_model.predict(impute_val_X), val_y))
# Lastly we impute and mark which fields were imputed.

marked_train_X = train_X.copy()
marked_val_X = val_X.copy()

for col in extra_columns:
    marked_train_X[col + "_was_missing"] = marked_train_X[col].isnull()
    marked_val_X[col + "_was_missing"] = marked_val_X[col].isnull()
    
imputer = Imputer(strategy = "mean", axis = 0)
marked_train_X = imputer.fit_transform(marked_train_X)
marked_val_X = imputer.transform(marked_val_X)

marked_model = RandomForestRegressor(random_state = 0)
marked_model.fit(marked_train_X, train_y)

print("Mean absolute error for the marked model:")
print(mae(marked_model.predict(marked_val_X), val_y))<jupyter_output>Mean absolute error for the marked model:
21812.31574741766<jupyter_text>In our case imputing values did not do much to help our model, perhaps because we chose the wrong predictors. Let's move on with including categorical fields into our model. For now we will consider only a one-hot encoding.<jupyter_code># First let's select only the categorical data and see its columns.

categorical_data = iowa_data.select_dtypes(include = object)
categorical_data.columns
# Let's chose some categorical columns and add them to our predictors.

categorical_columns = ["MS Zoning", "Neighborhood", "Bldg Type", "Pool QC", "Sale Type", "Sale Condition"]
predictors_with_categories = predictors + categorical_columns
X = iowa_data[predictors_with_categories]
encoded_X = pd.get_dummies(X)

# Let's make new training and validation sets.

train_X, val_X, train_y, val_y = split(encoded_X, y, random_state = 0)

encoded_model = RandomForestRegressor(random_state = 0)
encoded_model.fit(train_X, train_y)

print("Mean absolute error of the one-hot encoded model:")
print(mae(encoded_model.predict(val_X), val_y))<jupyter_output>Mean absolute error of the one-hot encoded model:
20189.36913532125<jupyter_text>Next up we'll look into another way of improving our models called *gradient boosting*. For this we will use the `XGBoost` library which includes a gradient boosted decision tree regressor. We will use the same test and validation sets as in the encoding example for comparison.<jupyter_code>from xgboost import XGBRegressor

boosted_model = XGBRegressor(silent = True, random_state = 0)
boosted_model.fit(train_X, train_y)
print("Mean absolute error of the gradient boosted decision tree model:")
print(mae(val_y, boosted_model.predict(val_X)))<jupyter_output>Mean absolute error of the gradient boosted decision tree model:
19189.39043315143<jupyter_text>For extra comparison let's also compare it to `scikit-learn`s built in implementation of a gradient boosted decision tree regressor.<jupyter_code>from sklearn.ensemble import GradientBoostingRegressor

boosted_model_2 = GradientBoostingRegressor(random_state = 0)
boosted_model_2.fit(train_X, train_y)
print("Mean absolute error of the built-in gradient boosted decision tree model:")
print(mae(val_y, boosted_model_2.predict(val_X)))<jupyter_output>Mean absolute error of the built-in gradient boosted decision tree model:
18990.403118871367<jupyter_text>Let's also explore some of the ways of tuning hyperparameters for the `XGBRegressor`. First we will try to find an optimal value for `n_estimators` using the `early_stopping_rounds` argument by stopping the iterations if we haven't improved the mean absolute error in `early_stopping_rounds` number of iterations.<jupyter_code>tuned_boosted_model = XGBRegressor(n_estimators = 2000, random_state = 0)
tuned_boosted_model.fit(train_X, train_y, early_stopping_rounds = 5,
                       eval_set = [(val_X, val_y)], eval_metric = "mae")
tuned_boosted_model.best_iteration<jupyter_output>[0]	validation_0-mae:163353
Will train until validation_0-mae hasn't improved in 5 rounds.
[1]	validation_0-mae:147121
[2]	validation_0-mae:132471
[3]	validation_0-mae:119235
[4]	validation_0-mae:107372
[5]	validation_0-mae:96811.6
[6]	validation_0-mae:87466.9
[7]	validation_0-mae:79103.9
[8]	validation_0-mae:71581.9
[9]	validation_0-mae:64935.2
[10]	validation_0-mae:58970.3
[11]	validation_0-mae:53687.4
[12]	validation_0-mae:49006.6
[13]	validation_0-mae:44858.3
[14]	validation_0-mae:41263.4
[15]	validation_0-mae:38225.3
[16]	validation_0-mae:35571.1
[17]	validation_0-mae:33324.4
[18]	validation_0-mae:31487.1
[19]	validation_0-mae:29970.5
[20]	validation_0-mae:28645.5
[21]	validation_0-mae:27371.6
[22]	validation_0-mae:26328.1
[23]	validation_0-mae:25502.1
[24]	validation_0-mae:24878
[25]	validation_0-mae:24281.9
[26]	validation_0-mae:23837.8
[27]	validation_0-mae:23370.3
[28]	validation_0-mae:22901.7
[29]	validation_0-mae:22654.6
[30]	validation_0-mae:22333.4
[31]	validation_0-mae:22[...]<jupyter_text>In this case we can see that 111 iterations give us the best model. We will continue trying to tune this model a little further, this time by adjusting the `learning_rate` parameter.<jupyter_code>import numpy as np

for rate in np.linspace(0.01, 0.1, 10):
    model = XGBRegressor(n_estimators = 2000, learning_rate = rate, random_state = 0)
    model.fit(train_X, train_y, early_stopping_rounds = 5,
             eval_set = [(val_X, val_y)], eval_metric = "mae",
             verbose = False)
    print("n iterations: {arg1}\tLearning rate: {arg2}\tMAE: {arg3}"
         .format(arg1 = model.best_iteration, arg2 = rate, arg3 = mae(val_y, model.predict(val_X))))<jupyter_output>n iterations: 901	Learning rate: 0.01	MAE: 19341.33800520123
n iterations: 406	Learning rate: 0.020000000000000004	MAE: 19535.59997442019
n iterations: 229	Learning rate: 0.030000000000000006	MAE: 19701.711854109824
n iterations: 222	Learning rate: 0.04000000000000001	MAE: 19414.564269270122
n iterations: 149	Learning rate: 0.05000000000000001	MAE: 19579.37616174966
n iterations: 208	Learning rate: 0.06000000000000001	MAE: 19036.778398917122
n iterations: 167	Learning rate: 0.07	MAE: 19052.983650238744
n iterations: 79	Learning rate: 0.08	MAE: 19706.30301415416
n iterations: 66	Learning rate: 0.09000000000000001	MAE: 19640.327624062073
n iterations: 111	Learning rate: 0.1	MAE: 19054.502057042973<jupyter_text>It seems that we get the best model using `n_estimators =  208` and `learning_rate = 0.06`, at least by only comparing the mean absolute error. However none turned out better than `scikit-learn`'s built in version, hence we will use this one, though we note that `xgboost` has other advantages. Now that we have a pretty good model let's try to understand what this model is actually telling us by looking at some *partial dependence plots*. Let's pick a couple predictors and see what our model tells us about how the price depends on them.<jupyter_code>from sklearn.ensemble.partial_dependence import partial_dependence as pdep, plot_partial_dependence as plot_pdep

boosted_model = GradientBoostingRegressor(random_state = 0)
boosted_model.fit(train_X, train_y)

plots = plot_pdep(boosted_model, 
                  features = [0, 1, 5], 
                  X = train_X, 
                  feature_names = list(train_X.columns))<jupyter_output><empty_output><jupyter_text>From these plots it seems, for example, that our model predicts higher prices the more modern the house is. Let us now try to put all of this together into what is called a *pipeline*. Pipelines string several steps together which helps us keep our code clean and, hopefully, bug free. To really illustrate the power of pipelines we will go through the entire process again: loading the data, creating training and validation sets, taking care of categorical and missing data, training and evaluating the model. This time let's also include all the columns.<jupyter_code>from sklearn.pipeline import make_pipeline

# Loading the data
iowa_data = pd.read_csv("./data/AmesHousing.txt", delimiter = '\t')
y = iowa_data.SalePrice
X = iowa_data.drop(["Order", "SalePrice"], axis = 1)

# One-hot encode the categorical data
categorical = X.select_dtypes(include = "object")
X = X.drop(categorical.columns, axis = 1).join(pd.get_dummies(categorical))

# Splitting the data
train_X, val_X, train_y, val_y = split(X, y, random_state = 0)

# Creating the pipeline
pipe = make_pipeline(Imputer(strategy = "mean", axis = 0), GradientBoostingRegressor(random_state = 0))

# Fitting the model
pipe.fit(train_X, train_y)

# Evaluating the model
mae(val_y, pipe.predict(val_X))<jupyter_output><empty_output><jupyter_text>Next up we do some cross validation.<jupyter_code>from sklearn.model_selection import cross_val_score

scores = cross_val_score(pipe, X, y, scoring = "neg_mean_absolute_error")
print("Mean absolute error: {arg1}".format(arg1 = -1 * scores.mean()))<jupyter_output>Mean absolute error: 15416.831349617949<jupyter_text>When making machine learning models one should take care to avoid *data leakage*. This mainly occurs through *leaky predictors* or *leaky validation strategies*. To avoid such issues make sure that you have the correct predictors, that they are all available when the time comes to make a decision and that you only fit transformers and models on training sets. You should always keep an eye out for a model that seems unreasonably accurate on the validation set. It may be due to data leakage and it's predictions may be much more inaccurate on new data. In our case we can see, for example, the predictors `Mo Sold` and `Yr Sold`. If our goal is to predict the sale price of a particular house then this data will not be available to us since the house hasn't been sold yet. It therefore makes sense to remove it from consideration.<jupyter_code>X = X.drop(["Mo Sold", "Yr Sold"], axis = 1)

scores = cross_val_score(pipe, X, y, scoring = "neg_mean_absolute_error")
print("Mean absolute error: {arg1}".format(arg1 = -1 * scores.mean()))<jupyter_output>Mean absolute error: 15349.564622068012
{'fraction_non_alphanumeric': 0.08314515662092525, 'fraction_numerical': 0.06809101962827861, 'mean_word_length': 5.461653572764684, 'pattern_counts': {'":': 0, '<': 55, '<?xml version=': 0, '>': 55, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 6, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '340896', 'license': 'Apache-2.0', 'n_tokens_mistral': 6479, 'n_tokens_neox': 5519, 'n_words': 1903, 'path': 'machineLearning.ipynb', 'repo_name': 'dlilja/kaggleLearn', 'size': '17271'}