<jupyter_start><jupyter_text>Notebook 2: Linear regression and random forests with scikit-learnNow that we have prepared our data we can start building machine learning models.We will start with the simplest model: Linear regression<jupyter_code>import pickle
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
import pandas as pd
import seaborn as sns
DATADIR = './dataset1/'
# Copy from previous notebook
def create_sub(preds, fn=None):
    df =  pd.DataFrame({'Id': range(len(preds)), 'Expected': preds})
    if fn is not None: df.to_csv(DATADIR + fn, index=False)
    return df
# Load the preprocessed data
with open('./tmp/preproc_data.pkl', 'rb') as f:
    X_train, y_train, X_valid, y_valid, X_test = pickle.load(f)
with open('./tmp/dfs.pkl', 'rb') as f:
    df_train, df_test = pickle.load(f)<jupyter_output><empty_output><jupyter_text>Linear regression$y = \sum_i a_i x_i + b$  Compute baseline mseBefore we start let's check which score we have to beat.<jupyter_code>def mse(y_true, y_pred): return ((y_true - y_pred)**2).mean()
mse(y_valid, df_train.t2m_fc_mean[df_train.time >= '2015-01-01'])<jupyter_output><empty_output><jupyter_text>sklearn APIAll sklearn models work in exactly the same way.<jupyter_code>from sklearn.linear_model import LinearRegression
# Create the model object
lr = LinearRegression()
# Fit the model parameters with the training data
lr.fit(X_train, y_train)<jupyter_output><empty_output><jupyter_text>Now the parameters $a_i$ and $b$ have been fit to our training data. We can not look at the individual parameters.<jupyter_code>X_train.columns
lr.coef_
lr.intercept_<jupyter_output><empty_output><jupyter_text>Wouldn't you expect the coefficient for `t2m_fc_mean` to be around 0 and the bias to be close to one? Why is this not the case?<jupyter_code># Make predictions for validation set
preds = lr.predict(X_valid)
mse(preds, y_valid)
# Submit to Kaggle
df_sub = create_sub(lr.predict(X_test), 'lr.csv'); df_sub.head()<jupyter_output><empty_output><jupyter_text>The R$^2$ scoreThe R$^2$ score is a common skill metric in machine learning: https://en.wikipedia.org/wiki/Coefficient_of_determination<jupyter_code>lr.score(X_valid, y_valid)
def print_scores(m, X_train=X_train, X_valid=X_valid):
    print('Train R2 = ', m.score(X_train, y_train), 
          ', Valid R2 = ', m.score(X_valid, y_valid), ', Valid MSE = ', 
          mse(m.predict(X_valid), y_valid))
print_scores(lr)<jupyter_output>Train R2 =  0.9319430866570234 , Valid R2 =  0.9187372901848856 , Valid MSE =  3.2422963614329707<jupyter_text>We can see that the training and validation score are roughly similar, indicating that the linear regression doesn't seem to be very prone to overfitting.  Decicion treesLinear regression improves the score but we can do better. One very common machine learning algorithm are decision trees. They are capable of representing non-linear relationships. In scikit-learn we can use the exact same API as before.<jupyter_code>from sklearn.tree import DecisionTreeRegressor
from sklearn.tree import export_graphviz
from IPython.display import display
import IPython, graphviz, re<jupyter_output><empty_output><jupyter_text>Visualize a tree<jupyter_code>def draw_tree(t, df, size=10, ratio=0.6, precision=0):
    """ From fast.ai machine learning course.
    Draws a representation of a random forest in IPython.
    Parameters:
    -----------
    t: The tree you wish to draw
    df: The data used to train the tree. This is used to get the names of the features.
    """
    s=export_graphviz(t, out_file=None, feature_names=df.columns, filled=True,
                      special_characters=True, rotate=True, precision=precision)
    IPython.display.display(graphviz.Source(re.sub('Tree {',
       f'Tree {{ size={size}; ratio={ratio}', s)))<jupyter_output><empty_output><jupyter_text>Let's start with a very shallow tree which we can visualize. Again we are using the exact same sklearn API.<jupyter_code>dt = DecisionTreeRegressor(max_depth=3)
%time dt.fit(X_train, y_train)
draw_tree(dt, df_train.drop(['t2m_obs', 'time', 'station'], 1))<jupyter_output><empty_output><jupyter_text>Here we can see that the splits are entirely based on `t2m_fc_mean` which is not surpising.  Build a full treeNow let's actually build a complicated tree.<jupyter_code>dt = DecisionTreeRegressor()
%time dt.fit(X_train, y_train)
print_scores(dt)
dt.tree_.node_count<jupyter_output><empty_output><jupyter_text>So we have a perfect score on the training dataset, but a pretty awful score on the validation set. This is an extreme case of overfitting!But this is to be expected. We built our tree to perfectly fit the training set, so that only one sample remains in the final leafs.  Reduce the complexity of the treeTo combat this we can set the minimum number of samples in each final leaf to a higher value.<jupyter_code>dt = DecisionTreeRegressor(min_samples_leaf=200)
%time dt.fit(X_train, y_train)
print_scores(dt)<jupyter_output>Train R2 =  0.939008732905577 , Valid R2 =  0.9112729295089091 , Valid MSE =  3.5401164749291705<jupyter_text>The overfitting is less bad now but we still cannot get a better score than with linear regression. We need to do better!  Random forestsHere is where random forests are coming in. As the name suggests, a random forest is an ensemble of decision trees, where each tree is build with a random subsample of the features and training data. This should help to prevent overfitting.<jupyter_code>from sklearn.ensemble import RandomForestRegressor<jupyter_output><empty_output><jupyter_text>Fit a sample model Now we are fitting several decision trees. This can be done in parallel by several processors. `n_jobs = -1` uses all available processors.<jupyter_code>rf = RandomForestRegressor(n_estimators=10, n_jobs=-1)
%time rf.fit(X_train, y_train)
print_scores(rf)<jupyter_output>Train R2 =  0.9923784978318793 , Valid R2 =  0.9171492011905503 , Valid MSE =  3.3056594363252723<jupyter_text>Again, we see that we are overfitting quite badly. As before, we can tune the hyper-parameters to prevent overfitting  Hyperparameters  min_samples_leafThe minimum number of sample is each leaf.<jupyter_code>rf = RandomForestRegressor(n_estimators=10, n_jobs=-1, min_samples_leaf=100)
%time rf.fit(X_train, y_train)
print_scores(rf)<jupyter_output>Train R2 =  0.9468109911673573 , Valid R2 =  0.9192623595469146 , Valid MSE =  3.221346648017444<jupyter_text>n_estimatorsThe number of trees.<jupyter_code>rf = RandomForestRegressor(n_estimators=40, n_jobs=-1, min_samples_leaf=10)
%time rf.fit(X_train, y_train)
print_scores(rf)<jupyter_output>Train R2 =  0.9758792913984392 , Valid R2 =  0.9236987002474748 , Valid MSE =  3.044341335934703<jupyter_text>We are still getting a much better score on the training set. Partly we are still overfitting, but partly this is to be expected from our dataset and the way we chose the validation set. Can you explain why?Let's submit our new best model to kaggle.<jupyter_code># Submit to Kaggle
df_sub = create_sub(rf.predict(X_test), 'rf.csv'); df_sub.head()<jupyter_output><empty_output><jupyter_text>Feature importanceGetting a good prediction is nice, but ideally we would also understand why we are getting the prediction.For this we can look at feature importance. In feature importance, you randomly shuffle each column of the dataset and check by how much the prediction skill decreases.<jupyter_code>rf.feature_importances_
# Create a new dataframe for easy plotting
fi_df = pd.DataFrame(
    data={'Feature': X_train.columns, 'Feature importance': rf.feature_importances_}, 
    columns = ['Feature', 'Feature importance']
)
fi_df.sort_values('Feature importance', inplace=True, ascending=False)
sns.barplot(data=fi_df, x='Feature importance', y='Feature')
plt.xscale('log')<jupyter_output><empty_output><jupyter_text>As expected, by far the most important feature is `t2m_fc_mean`. Does the ordering of the remaining features make meteorological sense?There some other really cool ways of interpreting a random forest. For more information check out this notebook: https://github.com/fastai/fastai/blob/master/courses/ml1/lesson2-rf_interpretation.ipynb  Your turnAgain, here are some challenges for all of you that are bored ;)1. The different stations have very different characteristics. Try building a model (e.g. linear regression) for each station individually. See if you can improve the validation and test score.2. Try adding some time information (for example the month as a feature. Does this improve the forecast?3. Try to get as good a score as possible. Be creative and use Google!  Build a LR model for each station individually<jupyter_code>split_date = '2015-01-01'
stations_train = df_train.station[df_train.time < split_date]
stations_valid = df_train.station[df_train.time >= split_date]
stations_test = df_test.station
%%time
models = {}
for s in stations_train.unique():
    m = LinearRegression()
    m.fit(X_train[stations_train == s], y_train[stations_train == s])
    models[s] = m
preds = np.zeros(len(y_valid))
for s in stations_valid.unique():
    s_idxs = stations_valid == s
    if s in stations_train.unique():
        preds[s_idxs] = models[s].predict(X_valid[s_idxs])
    else:
        preds[s_idxs] = lr.predict(X_valid[s_idxs])
mse(preds, y_valid)
len(X_test), X_test.shape
test_preds = np.zeros(len(X_test))
for s in stations_test.unique():
    s_idxs = stations_test == s
    if s in stations_train.unique():
        test_preds[s_idxs] = models[s].predict(X_test[s_idxs])
    else:
        test_preds[s_idxs] = lr.predict(X_test[s_idxs])
# Submit to Kaggle
df_sub = create_sub(test_preds, 'lr_local.csv'); df_sub.head()<jupyter_output><empty_output><jupyter_text>Add month as information<jupyter_code>month = pd.DatetimeIndex(df_train['time']).month
X_train_month = X_train.copy()
X_valid_month = X_valid.copy()
X_train_month['month'] = month[df_train.time < split_date]
X_valid_month['month'] = month[df_train.time >= split_date]
rf = RandomForestRegressor(n_estimators=40, n_jobs=-1, min_samples_leaf=10)
%time rf.fit(X_train_month, y_train)
print_scores(rf, X_train_month, X_valid_month)
fi_df = pd.DataFrame(
    data={'Feature': X_train_month.columns, 'Feature importance': rf.feature_importances_}, 
    columns = ['Feature', 'Feature importance']
)
fi_df.sort_values('Feature importance', inplace=True, ascending=False)
sns.barplot(data=fi_df, x='Feature importance', y='Feature')
plt.xscale('log')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09409523809523809, 'fraction_numerical': 0.031047619047619046, 'mean_word_length': 5.629419191919192, 'pattern_counts': {'":': 0, '<': 74, '<?xml version=': 0, '>': 75, 'https://': 2, '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': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '664783', 'license': 'MIT', 'n_tokens_mistral': 3692, 'n_tokens_neox': 3436, 'n_words': 1057, 'path': 'nb2_sklearn-lr-rf_solution.ipynb', 'repo_name': 'waves2weather/W2W-ML-Workshop', 'size': '10500'}