Uno de los primeros pasos para construir un modelo es separar los datos en entrenamiento/test. Construimos el modelo con los datos de entrenamiento y comprobamos su rendimiento con los datos de test. El objetivo es aplicar nuestro modelo sobre datos que no ha visto (test) durante la fase de entrenamiento para ver si es capaz de generalizar correctamente.

Cuando calculamos el error durante la fase de entrenamiento, éste suele ser bastante optimista, ya que el modelo se ajusta con esos mismos datos. Este error es el llamado error de entrenamiento. Por contra, el error de test es el error que comete el modelo con los datos nuevos (test), que suele ser peor que el de entrenamiento. Una vez validado el modelo con los datos de test nuestro modelo está terminado.

El problema de este enfoque es que durante la fase de creación y entrenamiento del modelo la elección de parámetros, tipo de modelo (regresión lineal, random forest…)… lo hacemos a ciegas. No podemos comprobar, por ejemplo, si un modelo de regresión lineal generalizará mejor que un random forest. Lo único que tenemos como estimación son nuestros errores de entrenamiento.

Por este motivo, es habitual dividir a su vez los datos de entrenamiento en otros dos conjuntos: entrenamiento y validación. Según la forma en que se generan los subconjuntos de entrenamiento/validación tenemos diferentes estrategias de validación, entre las cuales están:

• Validación simple
  
  
• Leave One Out Cross-Validation (LOOCV)
  
  
• K-Fold Cross-Validation
  
  
• Repeated K-Fold Cross-Validation
  
  
• Bootstraping
  
  

Es la estrategia más sencilla. Consiste en separar aleatoriamente los datos en dos grupos, uno para el entrenamiento y otro para la validación. Aunque es la opción más simple, tiene dos problemas:

• El error depende mucho de la separación de los datos, con lo que puede tener mucha variabilidad según la separación.
  
  
• No aprovechamos todos los datos del conjunto de entrenamiento para construir nuestro modelo, ya que dejamos un grupo aparte (lo habitual es sobre el 20-30%).
  
  

Por ejemplo, vamos a usar el dataframe boston de Pandas para construir un modelo que prediga el precio de las casas en la ciudad de Boston.

import pandas as pd
import numpy as np

from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split

from sklearn.linear_model import LinearRegression
from sklearn import metrics

import matplotlib.pyplot as plt
import seaborn as sns

data = load_boston()
df = pd.DataFrame(data = data.data, columns = data.feature_names)
df['target'] = data.target

X = df.drop(columns = 'target', axis = 1)
y = df['target']

X_train, X_val, y_train, y_val = train_test_split(X, y, train_size   = 0.7)

model = LinearRegression()
model.fit(X_train, y_train)

y_predict = model.predict(X_val)

MSE = metrics.mean_squared_error(y_val, y_predict)

print("MSE:", MSE)

MSE: 24.931924983034474

Como vemos, nuestro modelo tiene un error (usando MSE como métrica) de 24.93. ¿Que pasaría si repitiéramos el proceso 100 veces?

errors = np.empty(shape = [100])

for i in range(0, 100):
    model = LinearRegression() 
    X_train, X_val, y_train, y_val = train_test_split(X, y, train_size   = 0.7)
    model.fit(X_train, y_train)
    y_predict = model.predict(X_val)
    MSE = metrics.mean_squared_error(y_val, y_predict)
    errors[i] = MSE

print("MSE: Máximo: %.2f. Mínimo: %.2f" % (errors.max(), errors.min()))

MSE: Media: Máximo: 37.85. Mínimo: 15.68

Podemos apreciar que, según la separación que se haga de los datos, nuestro error puede variar desde 15.68 a 37.85.

Esta estrategia es un método iterativo donde se empieza entrenando el modelo con todos los datos disponible excepto uno, que se usa como validación. Para evitar la variabilidad del error según la muestra escogida, el proceso se repite tantas veces como datos haya, dejando cada vez uno como validación. El error estimado será el promedio de todos los errores obtenidos.

La ventaja de este método es que reducimos la variabilidad del método anterior (validación simple), ya que usamos todos los datos tanto como entrenamiento como validación.

La principal desventaja es que puede ser muy costoso computacionalmente, ya que entrenamos y validamos el modelo tantas veces como datos tengamos.

Para nuestro ejemplo, vamos a usar dos métodos de sklearn: LeaveOneOut() y Cross_validate().

from sklearn.model_selection import LeaveOneOut

loocv = LeaveOneOut()
model = LinearRegression()

from sklearn.model_selection import cross_validate

scoring = ['neg_mean_squared_error']

scores = cross_validate(model, X, y, 
    cv=loocv, 
    scoring = scoring)

mse_mean = scores['test_neg_mean_squared_error'].mean()

print("MSE: mean = %.2f" % (mse_mean))

MSE: mean = -23.73

El método cross_validate acepta como parámetros el modelo, los datos (X, y), el método de validación cruzada (loocv, en nuestro caso) y un conjunto de métricas a evaluar (neg_mean_squared_error) y devuelve un array con los valores de cada métrica evaluada por iteración. Para acceder a los valores de cada métrica usamos test_score (en nuestro caso test_neg_mean_squared_error).

El valor que hemos obtenido (-23.73) es la media de todos los errores de cada iteración. Es similar al obtenido con validación simple, pero eliminamos el problema de la variabilidad.

Esta estrategia también es un método iterativo. En este caso, se divide el conjunto de entrenamiento en k subconjuntos del mismo tamaño aproximadamente. En cada iteración, utilizamos k-1 grupos para entrenar y 1 para validar. Este proceso se repite k veces, usando cada vez un grupo distinto como validación.

Una ventaja de este método con respecto a loocv es que se reduce el coste computacional. Otra ventaja es que recudimos el riesgo de overfitting. Con loocv al final el modelo está siendo entrenando varias veces con casi los mismos datos (sólo se diferencia uno cada vez), por los grupos están altamente correlacionados. Con k-fold los grupos son más independientes, con lo que reducimos la varianza al promediar el error.

Normalmente se emplea un valor de k = [5, 10].

El ejemplo en Python es casi igual que el anterior (loocv), usando el método KFold en lugar de loocv() y definiendo la k con el parámetro n_splits.

from sklearn.model_selection import KFold

kfcv = KFold(n_splits=5)
model = LinearRegression()

from sklearn.model_selection import cross_validate

scoring = ['neg_mean_squared_error']

scores = cross_validate(model, X, y, 
    cv=kfcv, 
    scoring = scoring)

mse_mean = scores['test_neg_mean_squared_error'].mean()

print("MSE: mean = %.2f" % (mse_mean))

MSE: mean = -37.13

Es exactamente igual al método k-Fold-Cross-Validation pero repitiendo el proceso completo n veces. Por ejemplo, 10-Fold-Cross-Validation con 5 repeticiones implica a un total de 50 iteraciones ajuste-validación, pero no equivale a un 50-Fold-Cross-Validation.

En este caso, usamos RepeatedKFold de sklearn. Además de definir la k tenemos que indicar el número de repeticiones con el parámetro n_repeats.

from sklearn.model_selection import RepeatedKFold

rkfcv = RepeatedKFold(n_splits = 5, n_repeats = 10)
model = LinearRegression()

from sklearn.model_selection import cross_validate

scoring = ['neg_mean_squared_error']

scores = cross_validate(model, X, y, 
    cv=rkfcv, 
    scoring = scoring)

mse_mean = scores['test_neg_mean_squared_error'].mean()

print("MSE: mean = %.2f" % (mse_mean))

MSE: mean = -23.85

Este método consiste en ir creando muestras boostrap, que son muestras del mismo tamaño que la muestra original seleccionando datos aleatorios con reposición (lo que quiere decir que un dato puede aparecer repetido en varias ocasiones de cada muestra bootstrap). Los datos que no han sido seleccionados reciben el nombre de out-of-bag (OOB) y son los utilizados para la validación en cada iteración.

Los pasos son los siguientes:

• Obtener una nueva muestra del mismo tamaño que la muestra original mediante muestro aleatorio con reposición.
  
  
• Ajustar el modelo empleando la nueva muestra generada en el paso 1.
  
  
• Calcular el error del modelo empleando aquellas observaciones de la muestra original que no se han incluido en la nueva muestra. A este error se le conoce como error de validación.
  
  
• Repetir el proceso n veces y calcular la media de los n errores de validación.
  
  
• Finalmente, y tras las n repeticiones, se ajusta el modelo final empleando todas las observaciones de entrenamiento originales.
  
  

En este caso, vamos a crear nuestro propio algoritmo de bootsraping utilizando el método choice de la librearía numpy. Este método selecciona un conjunto de muestras aleatorias de un array unidimensional (en nuestro caso será los índices de nuestro dataframe de Pandas). Además, le podemos indicar con el parámetro replace que queremos que la selección sea con reemplazo.

model = LinearRegression()

n_iterations = 100
n_items = df.shape[0]

errors = np.empty(shape = [n_iterations])

for i in range(n_iterations):
    #Seleccionar índices con reemplazo
    chosen_idx = np.random.choice(n_items, replace = True, size = n_items)
    #Crear los conjuntos de prueba y validación según los índices seleccionados
    df_train = df.iloc[chosen_idx]
    X_train = df_train.drop(columns = 'target', axis = 1)
    y_train = df_train['target']
    df_val = df.drop(chosen_idx, axis=0)
    X_val = df_train.drop(columns = 'target', axis = 1)
    y_val = df_train['target']
    #Entrenar el modelo
    model.fit(X_train, y_train)
    #Calcular el error
    y_predict = model.predict(X_val)
    MSE = metrics.mean_squared_error(y_val, y_predict)
    errors[i] = MSE

print("MSE: mean = %.2f" % (np.absolute(errors.mean())))

MSE: mean = 20.92

No existe un método de validación que supere al resto en todos los escenarios, la elección debe basarse en varios factores.

• Si el tamaño de la muestra es pequeño, se recomienda emplear repeated k-Fold-Cross-Validation, ya que consigue un buen equilibrio bias-varianza y, dado que no son muchas observaciones, el coste computacional no es excesivo.
  
  
• Si el objetivo principal es comparar modelos mas que obtener una estimación precisa de las métricas, se recomienda bootstrapping ya que tiene menos varianza.
  
  
• Si el tamaño muestral es muy grande, la diferencia entre métodos se reduce y toma más importancia la eficiencia computacional. En estos casos, 10-Fold-Cross-Validation simple es suficiente.
  
  

Estudio comparativo de los diferentes métodos

En temas anteriores vimos la regresión polinomial. Vamos a construir datos con un polinomio de grado 2:

import numpy as np

import matplotlib.pyplot as plt

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import PolynomialFeatures

m = 100
X = 6 * np.random.rand(m, 1) -3
y = 0.5 * X**2 + X + 2 +np.random.rand(m, 1)

Mostramos los datos en una gráfica:

Separamos los datos en entrenamiento y test y mostramos los datos de entrenamiento en una gráfica:

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)

Creamos un modelo de regresión lineal simple:

model = LinearRegression()

model.fit(X_train, y_train)

y_predict = model.predict(X_train)

Como vemos, la recta no se ajusta de forma óptima a nuestros datos. Vamos a ver el MSE de entrenamiento y de test:

rmse = mean_squared_error(y_true = y_train, y_pred = y_predict)
print("MSE entrenamiento: %.2f" % rmse)

y_predict = model.predict(X_test)
rmse = mean_squared_error(y_true = y_test, y_pred = y_predict)
print("MSE test: %.2f" % rmse)

MSE entrenamiento: 1.56
MSE test: 1.49

Nuestro modelo no es capaz de adecuarse correctamente a los datos. Lo que está haciendo es subajustando los datos (underfitting).

Vamos a crear ahora un modelo de regresión polinomial de grado 2:

poly_features = PolynomialFeatures(degree = 2)
X_poly = poly_features.fit_transform(X_train)

model.fit(X_poly, y_train)

y_predict = model.predict(X_poly)

Ahora sí que nuestro modelo se ajusta a los datos. Si miramos el MSE de entrenamiento y de test vemos que mejora mucho el rendimiento:

rmse = mean_squared_error(y_true = y_train, y_pred = y_predict)
print("MSE entrenamiento: %.2f" % rmse)

X_poly = poly_features.transform(X_test)
y_predict = model.predict(X_poly)
rmse = mean_squared_error(y_true = y_test, y_pred = y_predict)
print("MSE test: %.2f" % rmse)

MSE entrenamiento: 0.08
MSE test: 0.07

¿Que pasaría si aumentáramos el grado de polinomio? Por ejemplo, a 35:

poly_features = PolynomialFeatures(degree = 35)
X_poly = poly_features.fit_transform(X_train)

model.fit(X_poly, y_train)

y_predict = model.predict(X_poly)

Parece que no ajusta mal del todo en la mayoría de datos. Vamos a ver el MSE de entrenamiento y error:

rmse = mean_squared_error(y_true = y_train, y_pred = y_predict)
print("MSE entrenamiento: %.2f" % rmse)

X_poly = poly_features.transform(X_test)
y_predict = model.predict(X_poly)
rmse = mean_squared_error(y_true = y_test, y_pred = y_predict)
print("MSE test: %.2f" % rmse)

MSE entrenamiento: 0.17
MSE test: 46545.50

¿Qué está pasando? Nuestro modelo funciona bastante bien con los datos de entrenamiento (mucho mejor que con una regresión lineal simple), pero no generaliza nada bien. Ésto ocurre porque el modelo está sobreajustando mucho los datos (overfitting).

Obviamente, en nuestro caso está claro que el modelo que mejor se ajusta a nuestros datos es el segundo (el de grado 2), ya que hemos creado nosotros los datos con un polinomio de grado 2. Pero, en general, no sabemos que función ha generado los datos, con lo que no tenemos modo de saber, a priori, como de complejo tiene que ser nuestro modelo. ¿Cómo puedo saber si el modelo está sobreajustando o subajustando los datos?

Una forma de saberlo es utilizando la validación cruzada. Si un modelo tiene un buen rendimiento con los datos de entrenamiento pero malo con los de validación, está sobreajustando. Si se tiene un mal rendimiento en ambos casos, está subajustando.

Otra forma de saberlo es fijarse en las curvas de aprendizaje. Por ejemplo, vamos a crear una función que muestre el RMSE en función del tamaño del conjunto de entrenamiento:

def plot_learning_curves(model, X, y):
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
    train_errors, test_errors = [], []
    for m in range(1, len(X_train)):
        model.fit(X_train[:m], y_train[:m])
        y_train_predict = model.predict(X_train[:m])
        y_test_predict = model.predict(X_test)
        train_errors.append(mean_squared_error(y_train[:m], y_train_predict, squared = False))
        test_errors.append(mean_squared_error(y_test, y_test_predict, squared = False))

    figure = plt.figure(figsize = (15, 10))
    axes = figure.add_subplot()
    _ = axes.plot(train_errors, "r-+", linewidth = 1, label = "train")
    _ = axes.plot(test_errors, "b-", linewidth = 1, label = "test")
    axes.set_ylim(ymin=0,ymax=3)
    axes.legend(fontsize=15,facecolor='#CDCDCD',labelcolor="#000000")
    axes.set_xlabel('Tamaño conjunto entrenamiento', fontsize=25,labelpad=20)
    axes.set_ylabel('RMSE', fontsize=25,labelpad=20)

Nuestra función recibe un modelo y unos datos (X, y). Separa en entrenamiento y test y va entrenando el modelo con un conjunto de entrenamiento cada vez más grande. Por último, muestra el gráfico.

Aplicamos la función al primer modelo (regresión lineal simple):

model  = LinearRegression()
plot_learning_curves(model, X, y)

Si nos fijamos en el rendimiento de los datos de entrenamiento, vemos que cuando hay una o dos instancias, el modelo ajusta perfectamente. A medida que vamos aumentando el tamaño del conjunto de entrenamiento, el error va aumentando, hasta llegar a una especie de meseta. En los datos de test, nuestro modelo es incapaz de ajustar cuando hay pocos datos. Si aumentamos el tamaño de los datos de entrenamiento, va mejorando el rendimiento hasta llegar a otra meseta cerca de la otra curva. Esta gráfica es típica de modelo subajustados. Ambas curvas están cerca, han llegado a una meseta y el error es bastante elevado.

Vamos ahora con el segundo modelo (polinomio de grado 2):

poly_features = PolynomialFeatures(degree = 2)
X_poly = poly_features.fit_transform(X)

plot_learning_curves(model, X_poly, y)

En este caso, las curvas están más cerca y el error es mucho menor. Nuestro modelo está funcionando de manera correcta.

Último modelo (polinomio de grado 35):

poly_features = PolynomialFeatures(degree = 35)
X_poly = poly_features.fit_transform(X)

plot_learning_curves(model, X_poly, y)

Aquí vemos que el error de entrenamiento va aumentando poco a poco, pero se mantiene bastante bajo. Sin embargo, el error de test se mantiene muy alto y las curvas están muy separadas. Obviamente, si aumentáramos el tamaño de los datos de entrenamiento las dos curvas se acercarían cada vez más. Esta gráfica es típica de modelo sobreajustados.

Otra forma de intentar evitar el sobreajuste es, como veremos en el punto siguiente, la regularización.

Validación de modelos predictivos: Cross-validation, OneLeaveOut, Bootstraping by Joaquín Amat Rodrigo, available under a Attribution 4.0 International (CC BY 4.0)