"""
Cross-validation classes for GAM
Author: Luca Puggini
"""
from statsmodels.compat.python import with_metaclass
from abc import ABCMeta, abstractmethod
import itertools
import numpy as np
from statsmodels.gam.smooth_basis import GenericSmoothers, UnivariateGenericSmoother
class BaseCV(with_metaclass(ABCMeta)):
"""
Base class for cross validation over a given model
Computes the cross validation error of a given model. All the cross
validation classes can be derived from this one (e.g. GamCV, LassoCV,...)
Parameters
----------
cv_iterator : instance of cross-validation iterator
Iterator used to split the data into train and test indices.
endog : ndarray
dependent (response) variable of the model
exog : ndarray
design matrix, used only to determine the number of observations
passed to ``cv_iterator.split``
Attributes
----------
cv_iterator : instance of cross-validation iterator
Iterator used to split the data into train and test indices.
endog : ndarray
dependent (response) variable of the model
exog : ndarray
design matrix, used only to determine the number of observations
passed to ``cv_iterator.split``
train_test_cv_indices : generator
generator of (train_index, test_index) pairs produced by
``cv_iterator.split``
"""
def __init__(self, cv_iterator, endog, exog):
self.cv_iterator = cv_iterator
self.exog = exog
self.endog = endog
# TODO: cv_iterator.split only needs nobs from endog or exog
self.train_test_cv_indices = self.cv_iterator.split(
self.exog, self.endog, label=None
)
def fit(self, **kwargs):
"""
Compute the cross validation error over all train/test splits
Parameters
----------
**kwargs
Keyword arguments passed to the ``fit`` method of the
cross-validated model.
Returns
-------
cv_err : ndarray
array of the errors computed on each test fold
"""
cv_err = []
for train_index, test_index in self.train_test_cv_indices:
cv_err.append(self._error(train_index, test_index, **kwargs))
return np.array(cv_err)
@abstractmethod
def _error(self, train_index, test_index, **kwargs):
"""
Train the model on the train set and return the error on the test set
Parameters
----------
train_index : ndarray
index or boolean mask selecting the training observations
test_index : ndarray
index or boolean mask selecting the test observations
**kwargs
Keyword arguments passed to the ``fit`` method of the
cross-validated model.
Returns
-------
error : float
prediction error evaluated on the test set
"""
def _split_train_test_smoothers(x, smoothers, train_index, test_index):
"""
Split smoothers in test and train sets and create GenericSmoothers
Parameters
----------
x : ndarray
values of the independent variables for all the smoothers
smoothers : instance
additive smoother instance whose components are split
train_index : ndarray
index or boolean mask selecting the training observations
test_index : ndarray
index or boolean mask selecting the test observations
Returns
-------
train_multivariate_smoothers : GenericSmoothers
smoothers evaluated on the training observations
test_multivariate_smoothers : GenericSmoothers
smoothers evaluated on the test observations
Notes
-----
This does not take exog_linear into account.
"""
train_smoothers = []
test_smoothers = []
for smoother in smoothers.smoothers:
train_basis = smoother.basis[train_index]
train_der_basis = smoother.der_basis[train_index]
train_der2_basis = smoother.der2_basis[train_index]
train_cov_der2 = smoother.cov_der2
# TODO: Double check this part. cov_der2 is calculated with all data
train_x = smoother.x[train_index]
train_smoothers.append(
UnivariateGenericSmoother(
train_x,
train_basis,
train_der_basis,
train_der2_basis,
train_cov_der2,
smoother.variable_name + " train",
)
)
test_basis = smoother.basis[test_index]
test_der_basis = smoother.der_basis[test_index]
test_cov_der2 = smoother.cov_der2
# TODO: Double check this part. cov_der2 is calculated with all data
test_x = smoother.x[test_index]
test_smoothers.append(
UnivariateGenericSmoother(
test_x,
test_basis,
test_der_basis,
train_der2_basis,
test_cov_der2,
smoother.variable_name + " test",
)
)
train_multivariate_smoothers = GenericSmoothers(x[train_index], train_smoothers)
test_multivariate_smoothers = GenericSmoothers(x[test_index], test_smoothers)
return train_multivariate_smoothers, test_multivariate_smoothers
class MultivariateGAMCV(BaseCV):
"""
Cross validation error of a multivariate additive model
Parameters
----------
smoother : additive smoother instance
smoother providing the basis for the smooth terms of the model
alphas : list of float
penalty weights, one for each smooth term
gam : model class
model class for creating a model with the training data
cost : function
cost function used to compute the prediction error
endog : ndarray
dependent (response) variable of the model
exog : ndarray or None
design matrix of the linear (non-smooth) terms of the model
cv_iterator : instance of cross-validation iterator
iterator used to split the data into train and test indices
Attributes
----------
cost : function
cost function used to compute the prediction error
gam : model class
model class for creating a model with the training data
smoother : additive smoother instance
smoother providing the basis for the smooth terms of the model
exog_linear : ndarray or None
design matrix of the linear (non-smooth) terms of the model
alphas : list of float
penalty weights, one for each smooth term
cv_iterator : instance of cross-validation iterator
iterator used to split the data into train and test indices
"""
def __init__(self, smoother, alphas, gam, cost, endog, exog, cv_iterator):
self.cost = cost
self.gam = gam
self.smoother = smoother
self.exog_linear = exog
self.alphas = alphas
self.cv_iterator = cv_iterator
# TODO: super does not do anything with endog, exog, except get nobs
# refactor to clean up what where `exog` and `exog_linear` is attached
# exog is not used in super
super().__init__(cv_iterator, endog, self.smoother.basis)
def _error(self, train_index, test_index, **kwargs):
"""
Train the model on the train set and return the error on the test set
Parameters
----------
train_index : ndarray
index or boolean mask selecting the training observations
test_index : ndarray
index or boolean mask selecting the test observations
**kwargs
Keyword arguments passed to the ``fit`` method of the gam model.
Returns
-------
error : float
prediction error, evaluated with ``cost``, on the test set
"""
train_smoother, test_smoother = _split_train_test_smoothers(
self.smoother.x, self.smoother, train_index, test_index
)
endog_train = self.endog[train_index]
endog_test = self.endog[test_index]
if self.exog_linear is not None:
exog_linear_train = self.exog_linear[train_index]
exog_linear_test = self.exog_linear[test_index]
else:
exog_linear_train = None
exog_linear_test = None
gam = self.gam(
endog_train,
exog=exog_linear_train,
smoother=train_smoother,
alpha=self.alphas,
)
gam_res = gam.fit(**kwargs)
# exog_linear_test and test_smoother.basis will be column_stacked
# but not transformed in predict
endog_est = gam_res.predict(
exog_linear_test, test_smoother.basis, transform=False
)
return self.cost(endog_test, endog_est)
class BasePenaltiesPathCV(with_metaclass(ABCMeta)):
"""
Base class for cross validation over a grid of parameters
The best parameter is saved in ``alpha_cv``.
Parameters
----------
alphas : list of float
grid of penalty weights to search over
Attributes
----------
alphas : list of float
grid of penalty weights to search over
alpha_cv : float or None
best penalty weight found by cross-validation
cv_error : ndarray or None
cross-validation error for each value in ``alphas``
cv_std : ndarray or None
standard deviation of the cross-validation error for each value in
``alphas``
Notes
-----
This class is currently not used.
"""
def __init__(self, alphas):
self.alphas = alphas
self.alpha_cv = None
self.cv_error = None
self.cv_std = None
def plot_path(self):
"""Plot the cross validation error and standard deviation over the alphas grid"""
from statsmodels.graphics.utils import _import_mpl
plt = _import_mpl()
plt.plot(self.alphas, self.cv_error, c="black")
plt.plot(self.alphas, self.cv_error + 1.96 * self.cv_std, c="blue")
plt.plot(self.alphas, self.cv_error - 1.96 * self.cv_std, c="blue")
plt.plot(self.alphas, self.cv_error, "o", c="black")
plt.plot(self.alphas, self.cv_error + 1.96 * self.cv_std, "o", c="blue")
plt.plot(self.alphas, self.cv_error - 1.96 * self.cv_std, "o", c="blue")
# TODO add return
[docs]
class MultivariateGAMCVPath:
"""
K-fold cross-validation for GAM
Parameters
----------
smoother : additive smoother instance
smoother providing the basis for the smooth terms of the model
alphas : list of iterables
list of alpha for smooths. The product space will be used as alpha
grid for cross-validation
gam : model class
model class for creating a model with k-fold training data
cost : function
cost function for the prediction error
endog : ndarray
dependent (response) variable of the model
exog : ndarray or None
design matrix of the linear (non-smooth) terms of the model
cv_iterator : instance of cross-validation iterator
iterator used to split the data into train and test indices
Attributes
----------
cost : function
cost function for the prediction error
smoother : additive smoother instance
smoother providing the basis for the smooth terms of the model
gam : model class
model class for creating a model with k-fold training data
alphas : list of iterables
list of alpha for smooths
alphas_grid : list of tuples
product space of ``alphas`` searched during cross-validation
endog : ndarray
dependent (response) variable of the model
exog : ndarray or None
design matrix of the linear (non-smooth) terms of the model
cv_iterator : instance of cross-validation iterator
iterator used to split the data into train and test indices
cv_error : ndarray
mean cross-validation error for each element of ``alphas_grid``
cv_std : ndarray
standard deviation of the cross-validation error for each element
of ``alphas_grid``
alpha_cv : tuple or None
element of ``alphas_grid`` with the smallest cross-validation error
Warnings
--------
The API of this class is preliminary and will change.
"""
def __init__(self, smoother, alphas, gam, cost, endog, exog, cv_iterator):
self.cost = cost
self.smoother = smoother
self.gam = gam
self.alphas = alphas
self.alphas_grid = list(itertools.product(*self.alphas))
self.endog = endog
self.exog = exog
self.cv_iterator = cv_iterator
self.cv_error = np.zeros(
shape=(
len(
self.alphas_grid,
)
)
)
self.cv_std = np.zeros(
shape=(
len(
self.alphas_grid,
)
)
)
self.alpha_cv = None
[docs]
def fit(self, **kwargs):
"""
Perform the cross-validation search over the alphas grid
Parameters
----------
**kwargs
Keyword arguments passed to the ``fit`` method of the gam model.
Returns
-------
self : MultivariateGAMCVPath
instance with ``cv_error``, ``cv_std`` and ``alpha_cv`` set
"""
for i, alphas_i in enumerate(self.alphas_grid):
gam_cv = MultivariateGAMCV(
smoother=self.smoother,
alphas=alphas_i,
gam=self.gam,
cost=self.cost,
endog=self.endog,
exog=self.exog,
cv_iterator=self.cv_iterator,
)
cv_err = gam_cv.fit(**kwargs)
self.cv_error[i] = cv_err.mean()
self.cv_std[i] = cv_err.std()
self.alpha_cv = self.alphas_grid[np.argmin(self.cv_error)]
return self
