Both regularization methods, ridge regression and LASSO, are implemented in the package scikit-learn.
These methods fit a generalized linear model via penalized maximum likelihood. The degree of regularization depends on the regularization parameter $\lambda$, in scikit-learn given by the alpha
argument. Thus, it is useful to evaluate the regression function for a sequence of $\lambda$s. We evaluate a logarithmic sequence of 100 $\lambda$ values.
Let us load and prepare the data first. We assign the response vector and the model matrix of the training set and the test set to the variables y_train
, y_test
, X_train
, and X_test
.
from sklearn import linear_model
from sklearn.model_selection import train_test_split
from sklearn.datasets import make_regression
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import folium
dwd = pd.read_table(
"https://userpage.fu-berlin.de/soga/data/raw-data/DWD.csv",
index_col=0,
sep=',',
)
df = dwd.drop([
'DWD_ID',
'RECORD_LENGTH',
'STATION_NAME',
'FEDERAL_STATE',
'PERIOD',
# 'LON',
# 'LAT',
], axis=1).dropna()
X = df.drop('MEAN_ANNUAL_RAINFALL', axis=1)
y = df['MEAN_ANNUAL_RAINFALL']
scaler = StandardScaler()
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)
rmses = pd.read_feather('https://userpage.fu-berlin.de/soga/data/py-data/30222_rmses.feather').set_index('index')
rmses
name | train_RMSE | test_RMSE | |
---|---|---|---|
index | |||
0 | baseline model | 243.882152 | 180.877011 |
0 | simple alt model | 154.992815 | 138.854544 |
0 | max rainfall model | 119.953630 | 117.437897 |
0 | multi alt rain model | 118.095746 | 113.746363 |
4 | forward model | 87.392377 | 117.884639 |
5 | backward model | 87.236023 | 118.367974 |
6 | mlxtend SFS model | 91.976905 | 91.232370 |
We calculate coefficients for $alpha \in [10^{-1}, 10^5]$ with the ridge regression and plot the results.
n_alphas = 100
alphas = np.logspace(-1, 5, n_alphas)
coefs = []
for a in alphas:
ridge = linear_model.Ridge(alpha=a)
results = ridge.fit(X_train, y_train)
coefs.append(ridge.coef_)
ax = plt.gca()
ax.plot(alphas, coefs)
plt.xlabel("$\lambda$")
ax.set_xscale("log")
plt.ylabel("Coefficients")
plt.title("Ridge coefficients as a function of the regularization")
plt.axis("tight")
plt.grid()
plt.show()
Let us have a look at the coefficient difference to simple linear regressions for each parameter. To this end, we calculate simple linear regression coeffs reg_coeff
and compute the mean squared error and plot it.
reg_coeff = linear_model.LinearRegression().fit(X_train, y_train)
errors = []
for a in alphas:
ridge = linear_model.Ridge(alpha=a)
results = ridge.fit(X_train, y_train)
errors.append(mean_squared_error(ridge.coef_, reg_coeff.coef_))
ax = plt.gca()
ax.plot(alphas, errors)
ax.set_xscale("log")
plt.xlabel("log($\lambda$)")
plt.ylabel("error")
plt.title("Coefficient deviation as a function of the regularization")
plt.axis("tight")
plt.grid()
plt.show()
The plot shows for a sequence of values of $\alpha$ the values for the Ridge model coefficients. For large values of $\alpha$ the coefficients correspond to the least square solution. By decreasing $\alpha$ the coefficients are penalized and~ converge to $0$.
In the next step we want to find a model configuration, and thus the value of $\alpha$ that minimizes the residuals (difference between observations and prediction). The scikit-learn package ships with a cross-validation){target="_blank"} function. The class RidgeCV()
performs a leave-one-out cross-validation on the data set and stores the MSE values in .cv_values_
if the argument store_cv_values
is active. We can evaluate the model performance by calculating the mean-squared error (MSE){target="_blank"}.
By visualizing, we can demonstrate the relation between MSE and $\alpha$.
reg = linear_model.RidgeCV(alphas=alphas, store_cv_values=True)
reg = reg.fit(X_train, y_train)
ax = plt.gca()
ax.plot(alphas, np.mean(reg.cv_values_, axis=0))
ax.set_xscale("log")
plt.xlabel("alpha")
plt.ylabel("mean-squared error")
plt.title("MSE as a function of the regularization")
plt.axis("tight")
plt.grid()
plt.show()
The plot shows that the MSE increases with $\alpha$. The dotted vertical lines indicate the value of $\alpha$ that gives the minimum mean cross-validated error and the largest value of $\alpha$ such that the error is within 1 standard error of the minimum. We can assess these values of $\alpha$ by the object attributes .alpha_
.
print("Best score:", reg.best_score_)
print("Estimated regularization parameter alpha:", reg.alpha_)
print("Training score:", reg.score(X_train, y_train))
print("Test score:", reg.score(X_test, y_test))
Best score: -6440.036134991497 Estimated regularization parameter alpha: 1.6297508346206444 Training score: 0.9081942915442142 Test score: 0.8625788494567455
Finally we may use the best $\alpha$ to predict the response vector y_train
and y_test
, respectively. Once again we calculate the RMSE, for both the training data and the test data, and store the model together with the corresponding RMSE in our list
object called model_outcome
.
# prediction for the training set
m_ridge_cv_pred_train = reg.predict(X_train)
# calculate RMSE on training set
rmse_test = mean_squared_error(m_ridge_cv_pred_train, y_train, squared=False)
print("RMSE on training set:", rmse_test)
RMSE on training set: 71.84011342174531
# prediction for the test set
m_ridge_cv_pred_test = reg.predict(X_test)
# calculate RMSE for the test data set
rmse_train = mean_squared_error(m_ridge_cv_pred_test, y_test, squared=False)
print("RMSE on training set:", rmse_train)
RMSE on training set: 77.40882425619203
rmses.loc[len(rmses)] = ['Ridge regresion', rmse_train, rmse_test]
rmses
name | train_RMSE | test_RMSE | |
---|---|---|---|
index | |||
0 | baseline model | 243.882152 | 180.877011 |
0 | simple alt model | 154.992815 | 138.854544 |
0 | max rainfall model | 119.953630 | 117.437897 |
0 | multi alt rain model | 118.095746 | 113.746363 |
4 | forward model | 87.392377 | 117.884639 |
5 | backward model | 87.236023 | 118.367974 |
6 | mlxtend SFS model | 91.976905 | 91.232370 |
7 | Ridge regresion | 77.408824 | 71.840113 |
The mechanics for the LASSO regression in Python are very similar to those for the ridge regression. Let us build the model and plot it.
n_alphas = 100
alphas = np.logspace(-2, 5, n_alphas)
coefs = []
for a in alphas:
ridge = linear_model.Lasso(alpha=a)
results = ridge.fit(X_train, y_train)
coefs.append(ridge.coef_)
ax = plt.gca()
ax.plot(alphas, coefs, )
plt.xlabel("$\lambda$")
ax.set_xscale("log")
plt.ylabel("Coefficients")
plt.title("Ridge coefficients as a function of the regularization")
plt.axis("tight")
plt.grid()
plt.show()
The plot shows for a sequence (default is 100) of $\alpha$ the values for the model coefficients. By increasing $\alpha$ the coefficients are penalized and eventually become $0$. The upper x-axis shows the decreasing number of coefficients of the model with increasing values of $\alpha$. This plot nicely visualizes variable selection via the LASSO method.
In the next step we need to find a model configuration, and thus the value of $\alpha$, which minimizes the residuals (difference between observations and prediction). The glmnet
package ships with a build-in cross-validation) function. The function cv.glmnet()
performs by default 10-fold cross-validation on the data set and evaluates the model performance by calculating the mean-squared error (MSE). Calling the plot()
function on the cv.glmnet
object provides a nice graphical representation of the relation between MSE and $\lambda$.
reg = linear_model.LassoCV(alphas=alphas, cv=10)
reg = reg.fit(X_train, y_train)
_, ax = plt.subplots(figsize=(10,4))
ax.scatter(alphas, reg.alphas)
ax.set_xscale("log")
plt.xlabel("$\lambda$")
plt.ylabel("mean-squared error")
plt.title("MSE as a function of the regularization")
plt.axis("tight")
plt.grid()
plt.show()
The plot shows that the MSE increases with $\lambda$. The dotted vertical lines indicate the value of $\lambda$ that gives the minimum mean cross-validated error and the largest value of $\lambda$ such that the error is within 1 standard error of the minimum. Sometimes it is more recommendable to choose the simpler model, and thus the model with fewer model parameters. We can assess these values of $\lambda$ by the object attributes lambda.min and lambda.1se.
print("Best score:", reg.alpha_)
print("Estimated regularization parameter alpha:", reg.alpha_)
print("Training score:", reg.score(X_train, y_train))
print("Test score:", reg.score(X_test, y_test))
Best score: 1.1233240329780276 Estimated regularization parameter alpha: 1.1233240329780276 Training score: 0.9064845447705463 Test score: 0.8748521084827312
reg.score(X_test, y_test)
0.8748521084827312
# prediction for the training set
lasso_pred_train = reg.predict(X_train)
# calculate RMSE on training set
rmse_test = mean_squared_error(m_lasso_pred_train, y_train, squared=False)
print("RMSE on training set:", rmse_test)
RMSE on training set: 72.50598591177226
# prediction for the test set
lasso_pred_test = reg.predict(X_test)
# calculate RMSE for the test data set
rmse_train = mean_squared_error(lasso_pred_test, y_test, squared=False)
print("RMSE on training set:", rmse_train)
RMSE on training set: 73.87125026527004
rmses.loc[len(rmses)] = ['Lasso regresion', rmse_train, rmse_test]
rmses
name | train_RMSE | test_RMSE | |
---|---|---|---|
index | |||
0 | baseline model | 243.882152 | 180.877011 |
0 | simple alt model | 154.992815 | 138.854544 |
0 | max rainfall model | 119.953630 | 117.437897 |
0 | multi alt rain model | 118.095746 | 113.746363 |
4 | forward model | 87.392377 | 117.884639 |
5 | backward model | 87.236023 | 118.367974 |
6 | mlxtend SFS model | 91.976905 | 91.232370 |
7 | Ridge regresion | 77.408824 | 71.840113 |
8 | Lasso regresion | 73.871250 | 72.505986 |
Let us visualize the final results and check which models are leading our "model selection comparison contest".
_, ax = plt.subplots(figsize=(10,6))
rmses.plot(kind="bar", x="name", ax=ax)
ax.set_xlabel("")
ax.set_ylabel("RMSE in mm")
plt.show()
The plot shows the RMSE of each model on the training data set and on the test data set. The best models with respect to RMSE on the test data set are those models, which used a stepwise parameter selection procedures or regularization terms. For the best performing models the RMSE on the test data set is comparable low.
At the end we draw a nice interactive map of the weather stations in our data set and color and scale them according to the RMSE. The benefit of using the RMSE is that the RMSE is on the same scale as the data, hence the RMSE can be interpreted as deviance of the predicted value from the observed value as rainfall in mm.
Feel free to interact with the map. You may zoom in and zoom out. You may change the underlying basemap. Further, you can hoover over the data points, and by clicking you get additional information about the data point. Let us look at predictions of the Ridge regression model.
# m_ridge_cv_pred_train
train_set = dwd.loc[y_train.index]
DWD_ID | STATION_NAME | FEDERAL_STATE | LAT | LON | ALTITUDE | PERIOD | RECORD_LENGTH | MEAN_ANNUAL_AIR_TEMP | MEAN_MONTHLY_MAX_TEMP | ... | MEAN_ANNUAL_WIND_SPEED | MEAN_CLOUD_COVER | MEAN_ANNUAL_SUNSHINE | MEAN_ANNUAL_RAINFALL | MAX_MONTHLY_WIND_SPEED | MAX_AIR_TEMP | MAX_WIND_SPEED | MAX_RAINFALL | MIN_AIR_TEMP | MEAN_RANGE_AIR_TEMP | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
ID | |||||||||||||||||||||
950 | 5825 | Berge | Brandenburg | 52.6198 | 12.7867 | 40.0 | 1953-2016 | 63 | 9.1 | 13.8 | ... | 2.0 | 64.0 | 1677.0 | 514.0 | 3.0 | 34.5 | 25.6 | 31.0 | -13.8 | 8.4 |
431 | 2503 | Kaltennordheim | Thüringen | 50.6266 | 10.1455 | 487.0 | 1951-2003 | 52 | 6.8 | 10.6 | ... | 2.0 | 72.0 | 1427.0 | 749.0 | 3.0 | 29.3 | 27.5 | 35.0 | -17.9 | 7.7 |
170 | 891 | Cuxhaven | Niedersachsen | 53.8713 | 8.7058 | 5.0 | 1922-2016 | 94 | 9.1 | 12.0 | ... | 3.0 | 68.0 | 1667.0 | 810.0 | 4.0 | 30.0 | 32.7 | 33.0 | -9.7 | 5.3 |
112 | 599 | Bonn-Friesdorf | Nordrhein-Westfalen | 50.7055 | 7.1467 | 62.0 | 1931-1999 | 68 | 9.9 | 14.1 | ... | 2.0 | 71.0 | 1437.0 | 659.0 | 3.0 | 32.1 | 24.0 | 36.0 | -11.9 | 8.1 |
440 | 2559 | Kempten | Bayern | 47.7233 | 10.3348 | 705.0 | 1951-2016 | 65 | 7.2 | 12.5 | ... | 2.0 | 66.0 | 1751.0 | 1238.0 | 2.0 | 31.0 | 24.2 | 52.0 | -20.1 | 9.8 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
494 | 2907 | Leck | Schleswig-Holstein | 54.7903 | 8.9514 | 7.0 | 1974-2016 | 42 | 8.4 | 11.8 | ... | 3.0 | 66.0 | 1503.0 | 778.0 | 4.0 | 29.4 | 30.3 | 32.0 | -13.6 | 7.2 |
74 | 403 | Berlin-Dahlem (FU) | Berlin | 52.4537 | 13.3017 | 51.0 | 1719-2016 | 297 | 9.0 | 13.3 | ... | 2.0 | 65.0 | 1700.0 | 592.0 | 3.0 | 33.4 | 27.7 | 35.0 | -13.9 | 8.1 |
429 | 2497 | Kall-Sistig | Nordrhein-Westfalen | 50.5014 | 6.5264 | 505.0 | 1947-2016 | 69 | 7.6 | 11.6 | ... | 2.0 | 70.0 | 1530.0 | 812.0 | 3.0 | 30.2 | 27.0 | 38.0 | -13.9 | 7.5 |
838 | 5078 | Travemünde | Schleswig-Holstein | 53.9619 | 10.8892 | 2.0 | 1936-2001 | 65 | 8.5 | 11.8 | ... | 2.0 | 65.0 | 1626.0 | 629.0 | 3.0 | 30.3 | 31.1 | 32.0 | -12.6 | 6.4 |
458 | 2667 | Köln-Bonn | Nordrhein-Westfalen | 50.8646 | 7.1575 | 92.0 | 1957-2016 | 59 | 10.1 | 14.5 | ... | 2.0 | 68.0 | 1515.0 | 797.0 | 3.0 | 33.4 | 26.3 | 36.0 | -13.3 | 9.0 |
163 rows × 21 columns
# extract response vector
y = dwd['MEAN_ANNUAL_RAINFALL']
train_set = dwd.loc[y_train.index]
test_set = dwd.loc[y_test.index]
rmse_train_vector = (m_ridge_cv_pred_train - y.loc[y_train.index])**2 / len(y_train)
rmse_test_vector = (m_ridge_cv_pred_test - y.loc[y_test.index])**2 / len(y_test)
loc_center = [df['LAT'].mean(), df['LON'].mean()]
germany_map = folium.Map(location = loc_center, tiles='Openstreetmap', zoom_start = 5, control_scale=True)
for index, loc in train_set.iterrows():
folium.CircleMarker(
[loc['LAT'], loc['LON']],
radius = np.log(rmse_train_vector.loc[index]),
color = 'red',
opacity = .7,
fill = True,
popup = f'training data, rmse: {rmse_train_vector.loc[index]}',
labels = ['training data'],
).add_to(germany_map)
for index, loc in test_set.iterrows():
folium.CircleMarker(
[loc['LAT'], loc['LON']],
radius = np.log(rmse_test_vector.loc[index]),
color = 'blue',
opacity = .7,
fill = True,
popup = f'test data, rmse: {rmse_test_vector.loc[index]}',
labels = 'test data',
).add_to(germany_map)
folium.LayerControl().add_to(germany_map)
germany_map
The map shows the RMSE for each data point (weather station). The spatial distribution of the RMSE shows that the prediction of mean annual rainfall is quite good in the North German Lowlands. In the area near Hamburg and Cuxhaven the model predictions are slightly worse. In contrast, the predictions of annual mean rainfall for the Central German Uplands, the Alpine foreland and the Alps are much worse. Maybe the input data for those mountainous regions is too sparse for capturing the spatial variability in the rainfall pattern.
Citation
The E-Learning project SOGA-Py was developed at the Department of Earth Sciences by Annette Rudolph, Joachim Krois and Kai Hartmann. You can reach us via mail by soga[at]zedat.fu-berlin.de.
Please cite as follow: Rudolph, A., Krois, J., Hartmann, K. (2023): Statistics and Geodata Analysis using Python (SOGA-Py). Department of Earth Sciences, Freie Universitaet Berlin.