Implementation¶
This section shows how the linear regression extensions discussed in this chapter are typically fit in Python. First let’s import the Boston housing dataset.
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn import datasets
boston = datasets.load_boston()
X_train = boston['data']
y_train = boston['target']
Regularized Regression¶
Both Ridge and Lasso regression can be easily fit using scikit-learn. A bare-bones implementation is provided below. Note that the regularization parameter alpha (which we called \(\lambda\)) is chosen arbitrarily.
from sklearn.linear_model import Ridge, Lasso
alpha = 1
# Ridge
ridge_model = Ridge(alpha = alpha)
ridge_model.fit(X_train, y_train)
# Lasso
lasso_model = Lasso(alpha = alpha)
lasso_model.fit(X_train, y_train);
In practice, however, we want to choose alpha through cross validation. This is easily implemented in scikit-learn by designating a set of alpha values to try and fitting the model with RidgeCV or LassoCV.
from sklearn.linear_model import RidgeCV, LassoCV
alphas = [0.01, 1, 100]
# Ridge
ridgeCV_model = RidgeCV(alphas = alphas)
ridgeCV_model.fit(X_train, y_train)
# Lasso
lassoCV_model = LassoCV(alphas = alphas)
lassoCV_model.fit(X_train, y_train);
We can then see which values of alpha performed best with the following.
print('Ridge alpha:', ridgeCV.alpha_)
print('Lasso alpha:', lassoCV.alpha_)
Ridge alpha: 0.01
Lasso alpha: 1.0
Bayesian Regression¶
We can also fit Bayesian regression using scikit-learn (though another popular package is pymc3). A very straightforward implementation is provided below.
from sklearn.linear_model import BayesianRidge
bayes_model = BayesianRidge()
bayes_model.fit(X_train, y_train);
This is not, however, identical to our construction in the previous section since it infers the \(\sigma^2\) and \(\tau\) parameters, rather than taking those as fixed inputs. More information can be found here. The hidden chunk below demonstrates a hacky solution for running Bayesian regression in scikit-learn using known values for \(\sigma^2\) and \(\tau\), though it is hard to imagine a practical reason to do so
By default, Bayesian regression in scikit-learn treats \(\alpha = \frac{1}{\sigma^2}\) and \(\lambda = \frac{1}{\tau}\) as random variables and assigns them the following prior distributions
Note that \(E(\alpha) = \frac{\alpha_1}{\alpha_2}\) and \(E(\lambda) = \frac{\lambda_1}{\lambda_2}\). To fix \(\sigma^2\) and \(\tau\), we can provide an extremely strong prior on \(\alpha\) and \(\lambda\), guaranteeing that their estimates will be approximately equal to their expected value.
Suppose we want to use \(\sigma^2 = 11.8\) and \(\tau = 10\), or equivalently \(\alpha = \frac{1}{11.8}\), \(\lambda = \frac{1}{10}\). Then let
This guarantees that \(\sigma^2\) and \(\tau\) will be approximately equal to their pre-determined values. This can be implemented in scikit-learn as follows
big_number = 10**5
# alpha
alpha = 1/11.8
alpha_1 = big_number*alpha
alpha_2 = big_number
# lambda
lam = 1/10
lambda_1 = big_number*lam
lambda_2 = big_number
# fit
bayes_model = BayesianRidge(alpha_1 = alpha_1, alpha_2 = alpha_2, alpha_init = alpha,
lambda_1 = lambda_1, lambda_2 = lambda_2, lambda_init = lam)
bayes_model.fit(X_train, y_train);
Poisson Regression¶
GLMs are most commonly fit in Python through the GLM class from statsmodels. A simple Poisson regression example is given below.
As we saw in the GLM concept section, a GLM is comprised of a random distribution and a link function. We identify the random distribution through the family argument to GLM (e.g. below, we specify the Poisson family). The default link function depends on the random distribution. By default, the Poisson model uses the link function
which is what we use below. For more information on the possible distributions and link functions, check out the statsmodels GLM docs.
import statsmodels.api as sm
X_train_with_constant = sm.add_constant(X_train)
poisson_model = sm.GLM(y_train, X_train, family=sm.families.Poisson())
poisson_model.fit();