Machine learning con Python y Scikit-learn
Scikit-learn¶
Python es uno de los lenguajes de programación que domina dentro del ámbito de la estadística, data mining y machine learning. Al tratarse de un software libre, innumerables usuarios han podido implementar sus algoritmos, dando lugar a un número muy elevado de librerías donde encontrar prácticamente todas las técnicas de machine learning existentes. Sin embargo, esto tiene un lado negativo, cada paquete tiene una sintaxis, estructura e implementación propia, lo que dificulta su aprendizaje. Scikit-learn, es una librería de código abierto que unifica bajo un único marco los principales algoritmos y funciones, facilitando en gran medida todas las etapas de preprocesado, entrenamiento, optimización y validación de modelos predictivos.
Scikit-learn ofrece tal cantidad de posibilidades que, difícilmente, pueden ser mostradas con un único ejemplo. En este documento, se emplean solo algunas de sus funcionalidades. Si en algún caso se requiere una explicación detallada, para que no interfiera con la narrativa del análisis, se añadirá un anexo. Aun así, para conocer bien todas las funcionalidades de Scikit-learn se recomienda leer su documentación.
Una de las características a destacar de esta librería es su elevado grado de madurez, lo que la hace adecuada para crear modelos predictivos que se quieren poner en producción.
Este documento está reproducido también con tidymodels, mlr3. Otros proyectos similares son caret y H2O, todos ellos basados en el lenguaje de programación R.
Introducción¶
Durante los últimos años, el interés y la aplicación de machine learning ha experimentado tal expansión que se ha convertido en una disciplina aplicada en prácticamente todos los ámbitos de investigación académica e industrial. El creciente número de personas dedicadas a esta disciplina ha dado como resultado todo un repertorio de herramientas con las que acceder a métodos predictivos potentes. El lenguaje de programación Python es un ejemplo de ello.
El término machine learning engloba al conjunto de algoritmos que permiten identificar patrones presentes en los datos y crear con ellos estructuras (modelos) que los representan. Una vez que los modelos han sido generados, se pueden emplear para predecir información sobre hechos o eventos que todavía no se han observado. Es importante recordar que, los sistemas de machine learning, "aprender" patrones que estén presentes en los datos con los que se entrenan, por lo tanto, solo pueden reconocer escenarios similares a lo que han visto antes. Al emplear sistemas entrenados con datos pasados para predecir futuros se está asumiendo que, en el futuro, el comportamiento será el similar, cosa que no siempre ocurre.
Aunque con frecuencia, términos como machine learning, data mining, inteligencia artificial, data science... son utilizados como sinónimos, es importante destacar que los métodos de machine learning son solo una parte de las muchas estrategias que se necesita combinar para extraer, entender y dar valor a los datos. El siguiente documento pretende ser un ejemplo del tipo de problema al que se suele enfrentar un analista: partiendo de un conjunto de datos más o menos procesado (la preparación de los datos es una etapa clave que precede al machine learning), se desea crear un modelo que permita predecir con éxito el comportamiento o valor que toman nuevas observaciones.
Etapas de un problema de machine learning
Definir el problema: ¿Qué se pretende predecir? ¿De qué datos se dispone? o ¿Qué datos es necesario conseguir?
Explorar y entender los datos que se van a emplear para crear el modelo.
Métrica de éxito: definir una forma apropiada de cuantificar cómo de buenos son los resultados obtenidos.
Preparar la estrategia para evaluar el modelo: separar las observaciones en un conjunto de entrenamiento, un conjunto de validación (o validación cruzada) y un conjunto de test. Es muy importante asegurar que ninguna información del conjunto de test participa en el proceso de entrenamiento del modelo.
Preprocesar los datos: aplicar las transformaciones necesarias para que los datos puedan ser interpretados por el algoritmo de machine learning seleccionado.
Ajustar un primer modelo capaz de superar unos resultados mínimos. Por ejemplo, en problemas de clasificación, el mínimo a superar es el porcentaje de la clase mayoritaria (la moda). En un modelo de regresión, la media de la variable respuesta.
Gradualmente, mejorar el modelo incorporando-creando nuevas variables u optimizando los hiperparámetros.
Evaluar la capacidad del modelo final con el conjunto de test para tener una estimación de la capacidad que tiene el modelo cuando predice nuevas observaciones.
Entrenar el modelo final con todos los datos disponibles.
A diferencia de otros documentos, este pretende ser un ejemplo práctico con menos desarrollo teórico. El lector podrá darse cuenta de lo sencillo que es aplicar un gran abanico de métodos predictivos con python y Scikit-learn. Sin embargo, es crucial que cualquier analista entienda los fundamentos teóricos en los que se basa cada uno de ellos para que un proyecto de este tipo tenga éxito. Aunque aquí solo se describen brevemente, estarán acompañados de links donde encontrar información detallada.
# Tratamiento de datos
# ==============================================================================
import numpy as np
import pandas as pd
from tabulate import tabulate
# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
from matplotlib import style
import matplotlib.ticker as ticker
import seaborn as sns
import statsmodels.api as sm
# Preprocesado y modelado
# ==============================================================================
from sklearn.model_selection import train_test_split
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import RepeatedKFold
from sklearn.neighbors import KNeighborsClassifier
from sklearn.datasets import make_blobs
from sklearn.metrics import euclidean_distances
from sklearn.neighbors import KNeighborsClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import Ridge
from skopt.space import Real, Integer
from skopt.utils import use_named_args
from skopt import gp_minimize
from skopt.plots import plot_convergence
import optuna
# Varios
# ==============================================================================
import multiprocessing
import random
from itertools import product
from fitter import Fitter, get_common_distributions
# Configuración matplotlib
# ==============================================================================
plt.rcParams['image.cmap'] = "bwr"
#plt.rcParams['figure.dpi'] = "100"
plt.rcParams['savefig.bbox'] = "tight"
style.use('ggplot') or plt.style.use('ggplot')
# Configuración warnings
# ==============================================================================
import warnings
warnings.filterwarnings('ignore')
Datos¶
El set de datos SaratogaHouses del paquete mosaicData de R contiene información sobre el precio de 1728 viviendas situadas en Saratoga County, New York, USA en el año 2006. Además del precio, incluye 15 variables adicionales:
price: precio de la vivienda.lotSize: metros cuadrados de la vivienda.age: antigüedad de la vivienda.landValue: valor del terreno.livingArea: metros cuadrados habitables.pctCollege: porcentaje del vecindario con título universitario.bedrooms: número de dormitorios.firplaces: número de chimeneas.bathrooms: número de cuartos de baño (el valor 0.5 hace referencia a cuartos de baño sin ducha).rooms: número de habitaciones.heating: tipo de calefacción.fuel: tipo de alimentación de la calefacción (gas, electricidad o diesel).sewer: tipo de desagüe.waterfront: si la vivienda tiene vistas al lago.newConstruction: si la vivienda es de nueva construcción.centralAir: si la vivienda tiene aire acondicionado.
Pueden descargarse los datos en formato csv de SaratogaHouses.csv
El objetivo es obtener un modelo capaz de predecir el precio del alquiler.
url = (
"https://raw.githubusercontent.com/JoaquinAmatRodrigo/Estadistica-machine-learning-python/"
"master/data/SaratogaHouses.csv"
)
datos = pd.read_csv(url, sep=",")
# Se renombran las columnas para que sean más descriptivas
datos.columns = ["precio", "metros_totales", "antiguedad", "precio_terreno", "metros_habitables",
"universitarios", "dormitorios", "chimenea", "banyos", "habitaciones",
"calefaccion","consumo_calefacion", "desague", "vistas_lago", "nueva_construccion",
"aire_acondicionado"]
Análisis exploratorio¶
Antes de entrenar un modelo predictivo, o incluso antes de realizar cualquier cálculo con un nuevo conjunto de datos, es muy importante realizar una exploración descriptiva de los mismos. Este proceso permite entender mejor qué información contiene cada variable, así como detectar posibles errores. Algunos ejemplos frecuentes son:
Que una columna se haya almacenado con el tipo incorrecto: una variable numérica está siendo reconocida como texto o viceversa.
Que una variable contenga valores que no tienen sentido: por ejemplo, para indicar que no se dispone del precio de una vivienda se introduce el valor 0 o un espacio en blanco.
Que en una variable de tipo numérico se haya introducido una palabra en lugar de un número.
Además, este análisis inicial puede dar pistas sobre qué variables son adecuadas como predictores en un modelo (más sobre esto en los siguientes apartados).
datos.head(4)
Tipo de cada columna¶
# Tipo de cada columna
# ==============================================================================
# En pandas, el tipo "object" hace referencia a strings
# datos.dtypes
datos.info()
Todas las columnas tienen el tipo adecuado.
Número de observaciones y valores ausentes¶
Junto con el estudio del tipo de variables, es básico conocer el número de observaciones disponibles y si todas ellas están completas. Los valores ausentes son muy importantes a la hora de crear modelos, la mayoría de algoritmos no aceptan observaciones incompletas o bien se ven muy influenciados por ellas. Aunque la imputación de valores ausentes es parte del preprocesado y, por lo tanto, debe de aprenderse únicamente con los datos de entrenamiento, su identificación se tiene que realizar antes de separar los datos para asegurar que se establecen todas las estrategias de imputación necesarias.
# Dimensiones del dataset
# ==============================================================================
datos.shape
# Número de datos ausentes por variable
# ==============================================================================
datos.isna().sum().sort_values()
Ninguna variable contiene valores ausente. En el apartado imputación de valores ausentes se muestra varias estrategias de imputación cuando el set de datos está incompleto.
Variable respuesta¶
Cuando se crea un modelo, es muy importante estudiar la distribución de la variable respuesta, ya que, a fin de cuentas, es lo que interesa predecir. La variable precio tiene una distribución asimétrica con una cola positiva debido a que, unas pocas viviendas, tienen un precio muy superior a la media. Este tipo de distribución suele visualizarse mejor tras aplicar el logarítmica o la raíz cuadrada.
fig, axes = plt.subplots(nrows=3, ncols=1, figsize=(6, 6))
sns.kdeplot(
datos.precio,
fill = True,
color = "blue",
ax = axes[0]
)
sns.rugplot(
datos.precio,
color = "blue",
ax = axes[0]
)
axes[0].set_title("Distribución original", fontsize = 'medium')
axes[0].set_xlabel('precio', fontsize='small')
axes[0].tick_params(labelsize = 6)
sns.kdeplot(
np.sqrt(datos.precio),
fill = True,
color = "blue",
ax = axes[1]
)
sns.rugplot(
np.sqrt(datos.precio),
color = "blue",
ax = axes[1]
)
axes[1].set_title("Transformación raíz cuadrada", fontsize = 'medium')
axes[1].set_xlabel('sqrt(precio)', fontsize='small')
axes[1].tick_params(labelsize = 6)
sns.kdeplot(
np.log(datos.precio),
fill = True,
color = "blue",
ax = axes[2]
)
sns.rugplot(
np.log(datos.precio),
color = "blue",
ax = axes[2]
)
axes[2].set_title("Transformación logarítmica", fontsize = 'medium')
axes[2].set_xlabel('log(precio)', fontsize='small')
axes[2].tick_params(labelsize = 6)
fig.tight_layout()
Algunos modelos de machine learning y aprendizaje estadístico requieren que la variable respuesta se distribuya de una forma determinada. Por ejemplo, para los modelos de regresión lineal (LM), la distribución tiene que ser de tipo normal. Para los modelos lineales generalizados (GLM), la distribución tiene que ser de la familia exponencial.
Existen varias librerías en python que permiten identificar a qué distribución se ajustan mejor los datos, una de ellas es fitter. Esta librería permite ajustar cualquiera de las 80 distribuciones implementadas en scipy.
distribuciones = ['cauchy', 'chi2', 'expon', 'exponpow', 'gamma',
'norm', 'powerlaw', 'beta', 'logistic']
fitter = Fitter(datos.precio, distributions=distribuciones)
fitter.fit()
fitter.summary(Nbest=10, plot=False)
Variables numéricas¶
# Variables numéricas
# ==============================================================================
datos.select_dtypes(include=['float64', 'int']).describe()
# Gráfico de distribución para cada variable numérica
# ==============================================================================
# Ajustar número de subplots en función del número de columnas
fig, axes = plt.subplots(nrows=3, ncols=3, figsize=(9, 5))
axes = axes.flat
columnas_numeric = datos.select_dtypes(include=['float64', 'int']).columns
columnas_numeric = columnas_numeric.drop('precio')
for i, colum in enumerate(columnas_numeric):
sns.histplot(
data = datos,
x = colum,
stat = "count",
kde = True,
color = (list(plt.rcParams['axes.prop_cycle'])*2)[i]["color"],
line_kws = {'linewidth': 2},
alpha = 0.3,
ax = axes[i]
)
axes[i].set_title(colum, fontsize = 7, fontweight = "bold")
axes[i].tick_params(labelsize = 6)
axes[i].set_xlabel("")
fig.tight_layout()
plt.subplots_adjust(top = 0.9)
fig.suptitle('Distribución variables numéricas', fontsize = 10, fontweight = "bold");
La variable chimenea, aunque es de tipo numérico, apenas toma unos pocos valores y la gran mayoría de observaciones pertenecen a solo dos de ellos. En casos como este, suele ser conveniente tratar la variable como cualitativa.
# Valores observados de chimenea
# ==============================================================================
datos.chimenea.value_counts()
# Se convierte la variable chimenea tipo string
# ==============================================================================
datos.chimenea = datos.chimenea.astype("str")
Como el objetivo del estudio es predecir el precio de las viviendas, el análisis de cada variable se hace también en relación a la variable respuesta precio. Analizando los datos de esta forma, se pueden empezar a extraer ideas sobre qué variables están más relacionadas con el precio y de qué forma.
# Gráfico de distribución para cada variable numérica
# ==============================================================================
# Ajustar número de subplots en función del número de columnas
fig, axes = plt.subplots(nrows=3, ncols=3, figsize=(9, 5))
axes = axes.flat
columnas_numeric = datos.select_dtypes(include=['float64', 'int']).columns
columnas_numeric = columnas_numeric.drop('precio')
for i, colum in enumerate(columnas_numeric):
sns.regplot(
x = datos[colum],
y = datos['precio'],
color = "gray",
marker = '.',
scatter_kws = {"alpha":0.4},
line_kws = {"color":"r","alpha":0.7},
ax = axes[i]
)
axes[i].set_title(f"precio vs {colum}", fontsize = 7, fontweight = "bold")
#axes[i].ticklabel_format(style='sci', scilimits=(-4,4), axis='both')
axes[i].yaxis.set_major_formatter(ticker.EngFormatter())
axes[i].xaxis.set_major_formatter(ticker.EngFormatter())
axes[i].tick_params(labelsize = 6)
axes[i].set_xlabel("")
axes[i].set_ylabel("")
# Se eliminan los axes vacíos
for i in [8]:
fig.delaxes(axes[i])
fig.tight_layout()
plt.subplots_adjust(top=0.9)
fig.suptitle('Correlación con precio', fontsize = 10, fontweight = "bold");
Correlación variables numéricas¶
Algunos modelos (LM, GLM, ...) se ven perjudicados si incorporan predictores altamente correlacionados. Por esta razón, es conveniente estudiar el grado de correlación entre las variables disponibles.
# Correlación entre columnas numéricas
# ==============================================================================
def tidy_corr_matrix(corr_mat):
'''
Función para convertir una matrix de correlación de pandas en formato tidy
'''
corr_mat = corr_mat.stack().reset_index()
corr_mat.columns = ['variable_1','variable_2','r']
corr_mat = corr_mat.loc[corr_mat['variable_1'] != corr_mat['variable_2'], :]
corr_mat['abs_r'] = np.abs(corr_mat['r'])
corr_mat = corr_mat.sort_values('abs_r', ascending=False)
return(corr_mat)
corr_matrix = datos.select_dtypes(include=['float64', 'int']).corr(method='pearson')
tidy_corr_matrix(corr_matrix).head(10)
# Heatmap matriz de correlaciones
# ==============================================================================
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(4, 4))
sns.heatmap(
corr_matrix,
annot = True,
cbar = False,
annot_kws = {"size": 6},
vmin = -1,
vmax = 1,
center = 0,
cmap = sns.diverging_palette(20, 220, n=200),
square = True,
ax = ax
)
ax.set_xticklabels(
ax.get_xticklabels(),
rotation = 45,
horizontalalignment = 'right',
)
ax.tick_params(labelsize = 8)
Variables cualitativas¶
# Variables cualitativas (tipo object)
# ==============================================================================
datos.select_dtypes(include=['object']).describe()
# Gráfico para cada variable cualitativa
# ==============================================================================
# Ajustar número de subplots en función del número de columnas
fig, axes = plt.subplots(nrows=3, ncols=3, figsize=(9, 5))
axes = axes.flat
columnas_object = datos.select_dtypes(include=['object']).columns
for i, colum in enumerate(columnas_object):
datos[colum].value_counts().plot.barh(ax = axes[i])
axes[i].set_title(colum, fontsize = 7, fontweight = "bold")
axes[i].tick_params(labelsize = 6)
axes[i].set_xlabel("")
# Se eliminan los axes vacíos
for i in [7, 8]:
fig.delaxes(axes[i])
fig.tight_layout()
plt.subplots_adjust(top=0.9)
fig.suptitle('Distribución variables cualitativas',
fontsize = 10, fontweight = "bold");
Si alguno de los niveles de una variable cualitativa tiene muy pocas observaciones en comparación a los otros niveles, puede ocurrir que, durante la validación cruzada o bootstrapping, algunas particiones no contengan ninguna observación de dicha clase (varianza cero), lo que puede dar lugar a errores. En estos casos, suele ser conveniente:
Eliminar las observaciones del grupo minoritario si es una variable multiclase.
Eliminar la variable si solo tiene dos niveles.
Agrupar los niveles minoritarios en un único grupo.
Asegurar que, en la creación de las particiones, todos los grupos estén representados en cada una de ellas.
Para este caso, hay que tener precaución con la variable chimenea. Se unifican los niveles de 2, 3 y 4 en un nuevo nivel llamado "2_mas".
datos.chimenea.value_counts().sort_index()
dic_replace = {'2': "2_mas",
'3': "2_mas",
'4': "2_mas"}
datos['chimenea'] = (
datos['chimenea']
.map(dic_replace)
.fillna(datos['chimenea'])
)
datos.chimenea.value_counts().sort_index()
# Gráfico relación entre el precio y cada cada variables cualitativas
# ==============================================================================
# Ajustar número de subplots en función del número de columnas
fig, axes = plt.subplots(nrows=3, ncols=3, figsize=(9, 5))
axes = axes.flat
columnas_object = datos.select_dtypes(include=['object']).columns
for i, colum in enumerate(columnas_object):
sns.violinplot(
x = colum,
y = 'precio',
data = datos,
color = "white",
ax = axes[i]
)
axes[i].set_title(f"precio vs {colum}", fontsize = 7, fontweight = "bold")
axes[i].yaxis.set_major_formatter(ticker.EngFormatter())
axes[i].tick_params(labelsize = 6)
axes[i].set_xlabel("")
axes[i].set_ylabel("")
# Se eliminan los axes vacíos
for i in [7, 8]:
fig.delaxes(axes[i])
fig.tight_layout()
plt.subplots_adjust(top=0.9)
fig.suptitle('Distribución del precio por grupo', fontsize = 10, fontweight = "bold");
División train y test¶
Evaluar la capacidad predictiva de un modelo consiste en comprobar cómo de próximas son sus predicciones a los verdaderos valores de la variable respuesta. Para poder cuantificarlo de forma correcta, se necesita disponer de un conjunto de observaciones, de las que se conozca la variable respuesta, pero que el modelo no haya "visto", es decir, que no hayan participado en su ajuste. Con esta finalidad, se dividen los datos disponibles en un conjunto de entrenamiento y un conjunto de test. El tamaño adecuado de las particiones depende en gran medida de la cantidad de datos disponibles y la seguridad que se necesite en la estimación del error, 80%-20% suele dar buenos resultados. El reparto debe hacerse de forma aleatoria o aleatoria-estratificada.
# Reparto de datos en train y test
# ==============================================================================
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(
datos.drop('precio', axis = 'columns'),
datos['precio'],
train_size = 0.8,
random_state = 1234,
shuffle = True
)
Es importante verificar que la distribución de la variable respuesta es similar en el conjunto de entrenamiento y en el de test. Para asegurar que esto se cumple, la función train_test_split() de scikit-learn permite, en problemas de clasificación, identificar con el argumento stratify la variable en base a la cual hacer el reparto.
Este tipo de reparto estratificado asegura que el conjunto de entrenamiento y el de test sean similares en cuanto a la variable respuesta, sin embargo, no garantiza que ocurra lo mismo con los predictores. Por ejemplo, en un set de datos con 100 observaciones, un predictor binario que tenga 90 observaciones de un grupo y solo 10 de otro, tiene un alto riesgo de que, en alguna de las particiones, el grupo minoritario no tenga representantes. Si esto ocurre en el conjunto de entrenamiento, algunos algoritmos darán error al aplicarlos al conjunto de test, ya que no entenderán el valor que se les está pasando. Este problema puede evitarse eliminando variables con varianza próxima a cero (ver más adelante).
print("Partición de entrenamento")
print("-----------------------")
print(y_train.describe())
print("Partición de test")
print("-----------------------")
print(y_test.describe())
Preprocesado¶
El preprocesado engloba todas aquellas transformaciones realizadas sobre los datos con el objetivo que puedan ser interpretados por el algoritmo de machine learning lo más eficientemente posible. Todo preprocesado de datos debe aprenderse con las observaciones de entrenamiento y luego aplicarse al conjunto de entrenamiento y al de test. Esto es muy importante para no violar la condición de que ninguna información procedente de las observaciones de test participe o influya en el ajuste del modelo. Este principio debe aplicarse también si se emplea validación cruzada (ver más adelante). En tal caso, el preprocesado debe realizarse dentro de cada iteración de validación, para que las estimaciones que se hacen con cada partición de validación no contengan información del resto de particiones. Aunque no es posible crear un único listado, a continuación se resumen algunos de los pasos de preprocesado que más se suelen necesitar.
Imputación de valores ausentes¶
La gran mayoría de algoritmos no aceptan observaciones incompletas, por lo que, cuando el set de datos contiene valores ausentes, se puede:
Eliminar aquellas observaciones que estén incompletas.
Eliminar aquellas variables que contengan valores ausentes.
Tratar de estimar los valores ausentes empleando el resto de información disponible (imputación).
Las primeras dos opciones, aunque sencillas, suponen perder información. La eliminación de observaciones solo puede aplicarse cuando se dispone de muchas y el porcentaje de registros incompletos es muy bajo. En el caso de eliminar variables, el impacto dependerá de cuánta información aporten dichas variables al modelo. Cuando se emplea la imputación, es muy importante tener en cuenta el riesgo que se corre al introducir valores en predictores que tengan mucha influencia en el modelo. Supóngase un estudio médico en el que, cuando uno de los predictores es positivo, el modelo predice casi siempre que el paciente está sano. Para un paciente cuyo valor de este predictor se desconoce, el riesgo de que la imputación sea errónea es muy alto, por lo que es preferible obtener una predicción basada únicamente en la información disponible. Esta es otra muestra de la importancia que tiene que el analista conozca el problema al que se enfrenta y pueda así tomar la mejor decisión.
El módulo sklearn.impute incorpora varios métodos de imputación distintos:
SimpleImputer: permite imputaciones empleando un valor constante o un estadístico (media, mediana, valor más frecuente) de la misma columna en la que se encuentra el valor ausente.IterativeImputer: permite imputar el valor de una columna teniendo en cuenta el resto de columnas. En concreto, se trata de un proceso iterativo en el que, en cada iteración, una de las variables se emplea como variable respuesta y el resto como predictores. Una vez obtenido el modelo, se emplea para predecir las posiciones vacías de esa variable. Este proceso se lleva a cabo con cada variable y se repite el ciclomax_iterveces para ganar estabilidad. La implementación desklearn.impute.IterativeImputerpermite que se emplee casi cualquiera de sus algoritmos para crear los modelos de imputación (KNN, RandomForest, GradientBoosting...).KNNImputer: es un caso concreto deIterativeImputeren el que se emplea k-Nearest Neighbors como algoritmo de imputación.
A pesar de ser un método muy utilizado, imputar utilizando KNN presenta dos problemas: su coste computacional elevado hace que solo sea aplicable en conjuntos de datos de tamaño pequeño o moderado. Si hay variables categóricas, debido a la dificultad de medir "distancias" en este contexto, puede dar lugar a resultados poco realistas. Por estas dos razones, es más recomendable utilizar un modelo tipo Random Forest IterativeImputer(predictor = RandomForestRegressor()).
Con el argumento add_indicator=True se crea automáticamente una nueva columna en la que se indica con el valor 1 qué valores han sido imputados. Esto puede ser útil tanto para identificar las observaciones en las que se ha realizado alguna imputación como para utilizarla como un predictor más en el modelo.
Exclusión de variables con varianza próxima a cero¶
No se deben incluir en el modelo predictores que contengan un único valor (cero-varianza) ya que no aportan información. Tampoco es conveniente incluir predictores que tengan una varianza próxima a cero, es decir, predictores que toman solo unos pocos valores, de los cuales, algunos aparecen con muy poca frecuencia. El problema con estos últimos es que pueden convertirse en predictores con varianza cero cuando se dividen las observaciones por validación cruzada o bootstrap.
La clase VarianceThreshold del módulo sklearn.feature_selection identifica y excluye todos aquellos predictores cuya varianza no supera un determinado threshold. En el caso de variables cualitativas, cabe recordar que scikitlearn requiere que se binaricen (one hot encoding o dummy) para poder entrenar los modelos. Una variable booleana sigue una distribución de Bernoulli, por lo que su varianza puede ser calculada como: $$\mathrm{Var}[X] = p(1 - p)$$
Si bien la eliminación de predictores no informativos podría considerarse un paso propio del proceso de selección de predictores, dado que consiste en un filtrado por varianza, tiene que realizarse antes de estandarizar los datos, ya que después, todos los predictores tienen varianza 1.
Estandarización y escalado de variables numéricas¶
Cuando los predictores son numéricos, la escala en la que se miden, así como la magnitud de su varianza pueden influir en gran medida en el modelo. Muchos algoritmos de machine learning (SVM, redes neuronales, lasso...) son sensibles a esto, de forma que, si no se igualan de alguna forma los predictores, aquellos que se midan en una escala mayor o que tengan más varianza dominarán el modelo aunque no sean los que más relación tienen con la variable respuesta. Existen principalmente 2 estrategias para evitarlo:
- Centrado: consiste en restarle a cada valor la media del predictor al que pertenece. Si los datos están almacenados en un dataframe, el centrado se consigue restándole a cada valor la media de la columna en la que se encuentra. Como resultado de esta transformación, todos los predictores pasan a tener una media de cero, es decir, los valores se centran en torno al origen.
StandardScaler(with_std=False)
Normalización (estandarización): consiste en transformar los datos de forma que todos los predictores estén aproximadamente en la misma escala. Hay dos formas de lograrlo:
Normalización Z-score (
StandardScaler): dividir cada predictor entre su desviación típica después de haber sido centrado, de esta forma, los datos pasan a tener una distribución normal.$$z = \frac{x - \mu}{\sigma}$$
Estandarización max-min (
MinMaxScaler): transformar los datos de forma que estén dentro del rango [0, 1].$$X_{norm} = \frac{X - X_{min}}{X_{max}-X_{min}}$$
Nunca se deben estandarizar las variables después de ser binarizadas (ver a continuación).
Binarización de las variables cualitativas¶
La binarización (one-hot-encoding) consiste en crear nuevas variables dummy con cada uno de los niveles de las variables cualitativas. Por ejemplo, una variable llamada color que contenga los niveles rojo, verde y azul, se convertirá en tres nuevas variables (color_rojo, color_verde, color_azul), todas con el valor 0 excepto la que coincide con la observación, que toma el valor 1.
Por defecto, la clase OneHotEncoder binariza todas las variables, por lo que hay que aplicarlo únicamente a las variables cualitativas (ver como hacerlo en el apartado ColumnTransformer). Con el argumento drop=‘first’ se elimina uno de los niveles para evitar redundancias. Volviendo al ejemplo anterior, no es necesario almacenar las tres variables, ya que, si color_rojo y color_verde toman el valor 0, la variable color_azul toma necesariamente el valor 1. Si color_rojo o color_verde toman el valor 1, entonces color_azul es necesariamente 0. Esto es importante en modelos que sufren problemas si los predictores están perfectamente correlacionados (modelos lineales sin regularización, redes neuronales...).
En ciertos escenarios puede ocurrir que, en los datos de test, aparezca un nuevo nivel que no estaba en los datos de entrenamiento. Si no se conoce de antemano cuáles son todos los posibles niveles, se puede evitar errores en las predicciones indicando OneHotEncoder(handle_unknown='ignore').
La forma de preprocesar los datos dentro del ecosistema scikit-learn es empleando los ColumnTransformer y pipeline. Además de las ya mencionadas, pueden encontrarse muchas más transformaciones de preprocesado en el módulo sklearn.preprocessing.
Pipeline y ColumnTransformer¶
Las clases ColumnTransformer y make_column_transformer del módulo sklearn.compose permiten combinar múltiples transformaciones de preprocesado, especificando a qué columnas se aplica cada una. Como todo transformer, tiene un método de entrenamiento (fit) y otro de transformación (transform) . Esto permite que el aprendizaje de las transformaciones se haga únicamente con observaciones de entrenamiento, y se puedan aplicar después a cualquier conjunto de datos. La idea detrás de este módulo es la siguiente:
Definir todas las transformaciones (escalado, selección, filtrado...) que se desea aplicar y a qué columnas
ColumnTransformer(). La selección de columnas puede hacerse por: nombre. índice, máscara booleana, slice, patrón regex, por tipo de columna o con las funciones de selecciónmake_column_selector.Aprender los parámetros necesarios para dichas transformaciones con las observaciones de entrenamiento
.fit().Aplicar las transformaciones aprendidas a cualquier conjunto de datos
.transform().
# Selección de las variables por tipo
# ==============================================================================
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OneHotEncoder
from sklearn.preprocessing import StandardScaler
from sklearn.compose import make_column_selector
# Se estandarizan las columnas numéricas y se hace one-hot-encoding de las
# columnas cualitativas. Para mantener las columnas a las que no se les aplica
# ninguna transformación se tiene que indicar remainder='passthrough'.
numeric_cols = X_train.select_dtypes(include=['float64', 'int']).columns.to_list()
cat_cols = X_train.select_dtypes(include=['object', 'category']).columns.to_list()
preprocessor = ColumnTransformer(
[('scale', StandardScaler(), numeric_cols),
('onehot', OneHotEncoder(handle_unknown='ignore'), cat_cols)],
remainder='passthrough')
Una vez que se ha definido el objeto ColumnTransformer, con el método fit() se aprenden las transformaciones con los datos de entrenamiento y se aplican a los dos conjuntos con transform().
X_train_prep = preprocessor.fit_transform(X_train)
X_test_prep = preprocessor.transform(X_test)
Por defecto, el resultado devuelto por ColumnTransformer es un numpy array, por lo que se pierden los nombres de las columnas. Suele ser interesante poder inspeccionar cómo queda el set de datos tras el preprocesado en formato dataframe. Por defecto, OneHotEncoder ordena las nuevas columnas de izquierda a derecha por orden alfabético.
# Convertir el output en dataframe y añadir el nombre de las columnas
# ==============================================================================
encoded_cat = preprocessor.named_transformers_['onehot'].get_feature_names_out(cat_cols)
nombre_columnas = np.concatenate([numeric_cols, encoded_cat])
X_train_prep = preprocessor.transform(X_train)
X_train_prep = pd.DataFrame(X_train_prep, columns=nombre_columnas)
X_train_prep.head(3)
A partir de la versión 1.2.1, scikit-learn incluye la posibilidad de que el resultado de los transformers sea un DataFrame de Pandas en lugar de un numpy array. Esto tiene la ventaja de que el resultado de toda transformación incorpora directamente el nombre de las columnas. Este comportamiento se puede indicar a nivel individual de transformer con el método .set_output(transform="pandas") o como una configuración global con set_config(transform_output="pandas").
Cabe destacar que, si se utiliza un transformer de tipo OneHotEncoder en combinación con .set_output(transform="pandas"), se tiene que indicar que el OneHotEncoder no devuelva los datos en forma de matriz sparse (sparse_output = False).
preprocessor = ColumnTransformer(
[('scale', StandardScaler(), numeric_cols),
('onehot', OneHotEncoder(handle_unknown='ignore', sparse_output=False), cat_cols)],
remainder = 'passthrough',
verbose_feature_names_out = False
).set_output(transform="pandas")
X_train_prep = preprocessor.fit_transform(X_train)
X_test_prep = preprocessor.transform(X_test)
X_train_prep.head(3)
ColumnTransformer aplica las operaciones de forma paralela, no de forma secuencial, esto significa que no permite aplicar más de una transformación a una misma columna. En el caso de que sea necesario hacerlo, hay que recurrir a los pipeline, que también agrupan operaciones pero las ejecutan de forma secuencial, de forma que la salida de una operación es la entrada de la siguiente. Si se quieren aplicar varias transformaciones de preprocesado sobre una misma columna, es necesario agruparlas primero en un pipeline
En el siguiente ejemplo se combinan las transormaciones:
Columnas numéricas: se imputan los valores ausentes con la mediana y a continuación se estandarizan.
Columnas categóricas (cualitativas): se imputan los valores ausentes con el valor más frecuente y a continuación se aplica one hot encoding.
# Selección de las variables por tipo
# ==============================================================================
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import OneHotEncoder
from sklearn.preprocessing import StandardScaler
from sklearn.compose import make_column_selector
numeric_cols = X_train.select_dtypes(include=['float64', 'int']).columns.to_list()
cat_cols = X_train.select_dtypes(include=['object', 'category']).columns.to_list()
# Transformaciones para las variables numéricas
numeric_transformer = Pipeline(
steps=[
('imputer', SimpleImputer(strategy='median')),
('scaler', StandardScaler())
]
)
# Transformaciones para las variables categóricas
categorical_transformer = Pipeline(
steps=[
('imputer', SimpleImputer(strategy='most_frequent')),
('onehot', OneHotEncoder(handle_unknown='ignore', sparse_output=False))
]
)
preprocessor = ColumnTransformer(
transformers=[
('numeric', numeric_transformer, numeric_cols),
('cat', categorical_transformer, cat_cols)
],
remainder='passthrough',
verbose_feature_names_out = False
).set_output(transform="pandas")
X_train_prep = preprocessor.fit_transform(X_train)
X_test_prep = preprocessor.transform(X_test)
X_train_prep = preprocessor.fit_transform(X_train)
X_test_prep = preprocessor.transform(X_test)
X_train_prep.head(3)
A partir de las versión scikit-learn 0.23 se puede crear una representación interactiva de un objeto pipeline.
from sklearn import set_config
set_config(display='diagram')
preprocessor
set_config(display='text')
Crear un modelo¶
El siguiente paso tras definir los datos de entrenamiento, es seleccionar el algoritmo que se va a emplear. En scikit-learn, esto se hace mediante la creación de un objeto estimator. En concreto, este objeto almacena el nombre del algoritmo, sus parámetros e hiperparámetros y contiene los métodos fit(X, y) y predict(T) que le permiten aprender de los datos y predecir nuevas observaciones. El siguiente listado contiene todos los algoritmos implementados en scikit-learn.
Entrenamiento¶
Se ajusta un primer modelo de regresión lineal con regularización ridge para predecir el precio de la vivienda en función de todos los predictores disponibles. Todos los argumentos de sklearn.linear_model.Ridge se dejan por defecto.
Es importante tener en cuenta que, cuando un modelo de regresión lineal incluye regularización en los coeficientes (ridge, lasso, elasticnet), deben estandarizarse los predictores. Para asegurar que el preprocesado se realiza únicamente con los datos de entrenamiento, se combinan las transformaciones y el entrenamiento en un mismo pipeline.
from sklearn.linear_model import Ridge
# Preprocedado
# ==============================================================================
# Identificación de columnas numéricas y categóricas
numeric_cols = X_train.select_dtypes(include=['float64', 'int']).columns.to_list()
cat_cols = X_train.select_dtypes(include=['object', 'category']).columns.to_list()
# Transformaciones para las variables numéricas
numeric_transformer = Pipeline(
steps=[('scaler', StandardScaler())]
)
# Transformaciones para las variables categóricas
categorical_transformer = Pipeline(
steps=[('onehot', OneHotEncoder(handle_unknown='ignore', sparse_output=False))]
)
preprocessor = ColumnTransformer(
transformers=[
('numeric', numeric_transformer, numeric_cols),
('cat', categorical_transformer, cat_cols)
],
remainder='passthrough',
verbose_feature_names_out = False
).set_output(transform="pandas")
# Pipeline
# ==============================================================================
# Se combinan los pasos de preprocesado y el modelo en un mismo pipeline
pipe = Pipeline([('preprocessing', preprocessor),
('modelo', Ridge())])
# Train
# ==============================================================================
# Se asigna el resultado a _ para que no se imprima por pantalla
_ = pipe.fit(X=X_train, y=y_train)
Validación¶
La finalidad última de un modelo es predecir la variable respuesta en observaciones futuras o en observaciones que el modelo no ha "visto" antes. El error mostrado por defecto tras entrenar un modelo suele ser el error de entrenamiento, el error que comete el modelo al predecir las observaciones que ya ha "visto". Si bien estos errores son útiles para entender cómo está aprendiendo el modelo (estudio de residuos), no es una estimación realista de cómo se comporta el modelo ante nuevas observaciones (el error de entrenamiento suele ser demasiado optimista). Para conseguir una estimación más certera, y antes de recurrir al conjunto de test, se pueden emplear estrategias de validación basadas en resampling. Scikit-learn incorpora en el módulo sklearn.model_selection varias estrategias de validación.
Todas ellas reciben como primer argumento un
estimatorque puede ser directamente un modelo o unpipeline.Las métricas de error de regresión se devuelven siempre en negativo de forma que, cuanto más próximo a 0 sea el valor, mejor el ajuste. Esto es así para que, los procesos de optimización siempre sean de maximización.
La forma más sencilla es emplear la función cross_val_score(), que utiliza por defecto KFold.
# Validación cruzada
# ==============================================================================
from sklearn.model_selection import cross_val_score
cv_scores = cross_val_score(
estimator = pipe,
X = X_train,
y = y_train,
scoring = 'neg_root_mean_squared_error',
cv = 5
)
print(f"Métricas validación cruzada: {cv_scores}")
print(f"Média métricas de validación cruzada: {cv_scores.mean()}")
También es posible utilizar otras estrategias de validación cruzada (KFold, RepeatedKFold, LeaveOneOut, LeavePOut, ShuffleSplit) realizando previamente el reparto con las funciones auxiliares de sklearn.model_selection y pasando los índices al argumento cv.
Cada método funciona internamente de forma distinta, pero todos ellos se basan en la idea: ajustar y evaluar el modelo de forma repetida, empleando cada vez distintos subconjuntos creados a partir de los datos de entrenamiento y obteniendo en cada repetición una estimación del error. El promedio de todas las estimaciones tiende a converger en el valor real del error de test. $^{\text{Anexo 1}}$
# Validación cruzada repetida
# ==============================================================================
from sklearn.model_selection import RepeatedKFold
cv = RepeatedKFold(n_splits=5, n_repeats=5, random_state=123)
cv_scores = cross_val_score(
estimator = pipe,
X = X_train,
y = y_train,
scoring = 'neg_root_mean_squared_error',
cv = cv
)
print(f"Métricas de validación cruzada: {cv_scores}")
print("")
print(f"Média métricas de validación cruzada: {cv_scores.mean()}")
La función cross_validate es similar a cross_val_score pero permite estimar varias métricas a la vez, tanto para test como para train, y devuelve los resultados en un diccionario.
# Validación cruzada repetida con múltiples métricas
# ==============================================================================
from sklearn.model_selection import cross_validate
cv = RepeatedKFold(n_splits=3, n_repeats=5, random_state=123)
cv_scores = cross_validate(
estimator = pipe,
X = X_train,
y = y_train,
scoring = ('r2', 'neg_root_mean_squared_error'),
cv = cv,
return_train_score = True
)
# Se convierte el diccionario a dataframe para facilitar la visualización
cv_scores = pd.DataFrame(cv_scores)
cv_scores
# Distribución del error de validación cruzada
# ==============================================================================
fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(6, 5), sharex=True)
sns.kdeplot(
cv_scores['test_neg_root_mean_squared_error'],
fill = True,
alpha = 0.3,
color = "firebrick",
ax = axes[0]
)
sns.rugplot(
cv_scores['test_neg_root_mean_squared_error'],
color = "firebrick",
ax = axes[0]
)
axes[0].set_title('neg_root_mean_squared_error', fontsize = 7, fontweight = "bold")
axes[0].tick_params(labelsize = 6)
axes[0].set_xlabel("")
sns.boxplot(
x = cv_scores['test_neg_root_mean_squared_error'],
ax = axes[1]
)
axes[1].set_title('neg_root_mean_squared_error', fontsize = 7, fontweight = "bold")
axes[1].tick_params(labelsize = 6)
axes[1].set_xlabel("")
fig.tight_layout()
plt.subplots_adjust(top=0.9)
fig.suptitle('Distribución error de validación cruzada', fontsize = 10,
fontweight = "bold");
La función cross_val_predict, en lugar de devolver la métrica de cada partición, devuelve las predicciones de cada partición. Esto es útil para poder evaluar los residuos del modelo y diagnosticar su comportamiento. Si se emplea validación cruzada repetida o bootstrapping, una misma observación puede formar parte de la partición de validación varias veces.
# Diagnóstico errores (residuos) de las predicciones de validación cruzada
# ==============================================================================
from sklearn.model_selection import cross_val_predict
from sklearn.model_selection import KFold
import statsmodels.api as sm
# Validación cruzada
# ==============================================================================
cv = KFold(n_splits=5, random_state=123, shuffle=True)
cv_prediccones = cross_val_predict(
estimator = pipe,
X = X_train,
y = y_train,
cv = cv
)
# Gráficos
# ==============================================================================
fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(9, 5))
axes[0, 0].scatter(y_train, cv_prediccones, edgecolors=(0, 0, 0), alpha = 0.4)
axes[0, 0].plot(
[y_train.min(), y_train.max()],
[y_train.min(), y_train.max()],
'k--', color = 'black', lw=2
)
axes[0, 0].set_title('Valor predicho vs valor real', fontsize = 10, fontweight = "bold")
axes[0, 0].set_xlabel('Real')
axes[0, 0].set_ylabel('Predicción')
axes[0, 0].tick_params(labelsize = 7)
axes[0, 1].scatter(list(range(len(y_train))), y_train - cv_prediccones,
edgecolors=(0, 0, 0), alpha = 0.4)
axes[0, 1].axhline(y = 0, linestyle = '--', color = 'black', lw=2)
axes[0, 1].set_title('Residuos del modelo', fontsize = 10, fontweight = "bold")
axes[0, 1].set_xlabel('id')
axes[0, 1].set_ylabel('Residuo')
axes[0, 1].tick_params(labelsize = 7)
sns.histplot(
data = y_train - cv_prediccones,
stat = "density",
kde = True,
line_kws= {'linewidth': 1},
color = "firebrick",
alpha = 0.3,
ax = axes[1, 0]
)
axes[1, 0].set_title('Distribución residuos del modelo', fontsize = 10,
fontweight = "bold")
axes[1, 0].set_xlabel("Residuo")
axes[1, 0].tick_params(labelsize = 7)
sm.qqplot(
y_train - cv_prediccones,
fit = True,
line = 'q',
ax = axes[1, 1],
color = 'firebrick',
alpha = 0.4,
lw = 2
)
axes[1, 1].set_title('Q-Q residuos del modelo', fontsize = 10, fontweight = "bold")
axes[1, 1].tick_params(labelsize = 7)
fig.tight_layout()
plt.subplots_adjust(top=0.9)
fig.suptitle('Diagnóstico residuos', fontsize = 12, fontweight = "bold");
A efectos prácticos, cuando se aplican métodos de resampling para validar un modelo hay que tener en cuenta dos cosas: el coste computacional que implica ajustar múltiples veces un modelo, cada vez con un subconjunto de datos distinto, y la reproducibilidad en la creación de las particiones. Las funciones cross_val_score() y cross_val_predict permiten paralelizar el proceso mediante el argumento n_jobs.
# Validación cruzada repetida paralelizada (multicore)
# ==============================================================================
from sklearn.model_selection import RepeatedKFold
cv = RepeatedKFold(n_splits=10, n_repeats=5, random_state=123)
cv_scores = cross_val_score(
estimator = pipe,
X = X_train,
y = y_train,
scoring = 'neg_root_mean_squared_error',
cv = cv,
n_jobs = -1 # todos los cores disponibles
)
print(f"Média métricas de validación cruzada: {cv_scores.mean()}")
El root_mean_squared_error promedio estimado mediante validación cruzada para el modelo ridge es de 56735. Este valor será contrastado más adelante cuando se calcule el error del modelo con el conjunto de test.
Predicción¶
Una vez que el modelo ha sido entrenado, bien empleando directamente un estimator o un pipeline, con el método .predict() se pueden predecir nuevas observaciones. Si se emplea un pipeline, se aplican automáticamente las transformaciones aprendidas durante el entrenamiento.
predicciones = pipe.predict(X_test)
# Se crea un dataframe con las predicciones y el valor real
df_predicciones = pd.DataFrame({'precio' : y_test, 'prediccion' : predicciones})
df_predicciones.head()
Error de test¶
Aunque mediante los métodos de validación (Kfold, LeaveOneOut) se consiguen buenas estimaciones del error que tiene un modelo al predecir nuevas observaciones, la mejor forma de evaluar un modelo final es prediciendo un conjunto test, es decir, un conjunto de observaciones que se ha mantenido al margen del proceso de entrenamiento y optimización. Dependiendo del problema en cuestión, pueden ser interesantes unas métricas$^{\text{Anexo 2}}$ u otras. El módulo sklearn.metrics incorpora una variedad considerable de métricas para evaluar la calidad de las predicciones.
# mean_squared_error de test
# ==============================================================================
from sklearn.metrics import mean_squared_error
rmse = mean_squared_error(
y_true = y_test,
y_pred = predicciones,
squared = False
)
rmse
En el apartado de validación, se estimó, mediante validación cruzada repetida, que el rmse del modelo era de 56735, un valor próximo al obtenido con el conjunto de test 65372.
Hiperparámetros (tuning)¶
Muchos modelos, entre ellos la regresión lineal con regularización Ridge, contienen parámetros que no pueden aprenderse a partir de los datos de entrenamiento y, por lo tanto, deben de ser establecidos por el analista. A estos se les conoce como hiperparámetros. Los resultados de un modelo pueden depender en gran medida del valor que tomen sus hiperparámetros, sin embargo, no se puede conocer de antemano cuál es el adecuado. Aunque con la práctica, los especialistas en machine learning ganan intuición sobre qué valores pueden funcionar mejor en cada problema, no hay reglas fijas. La forma más común de encontrar los valores óptimos es probando diferentes posibilidades.
Escoger un conjunto de valores para el o los hiperparámetros.
grid search: se hace una búsqueda exhaustiva sobre un conjunto de valores previamente definidos por el usuario.
random search: se evalúan valores aleatorios dentro de unos límites definidos por el usuario.
Para cada valor (combinación de valores si hay más de un hiperparámetro), entrenar el modelo y estimar su error mediante un método de validación $^{\text{Anexo 1}}$.
Finalmente, ajustar de nuevo el modelo, esta vez con todos los datos de entrenamiento y con los mejores hiperparámetros encontrados.
Scikilearn permite explorar diferentes valores de hiperparámetros mediante model_selection.GridSearchCV() y model_selection.RandomizedSearchCV()
Grid search¶
El modelo Ridge empleado hasta ahora tiene un hiperparámetro llamado alpha, que por defecto tiene el valor 1.0. Este hiperparámetro controla la penalización que se aplica a los coeficientes del modelo. Cuanto mayor es su valor, más restricción se impone sobre los coeficientes, reduciendo así varianza, atenuando el efecto de la correlación entre predictores y minimizando el riesgo de overfitting.
Se vuelve a ajustar un modelo Ridge con diferentes valores de alpha empleando validación cruzada repetida para identificar con cuál se obtienen mejores resultados.
from sklearn.model_selection import GridSearchCV, RepeatedKFold
from sklearn.linear_model import Ridge
# Pipe: preprocesado + modelo
# ==============================================================================
# Identificación de columnas numéricas y catégoricas
numeric_cols = X_train.select_dtypes(include=['float64', 'int']).columns.to_list()
cat_cols = X_train.select_dtypes(include=['object', 'category']).columns.to_list()
# Transformaciones para las variables numéricas
numeric_transformer = Pipeline(
steps=[('scaler', StandardScaler())]
)
# Transformaciones para las variables categóricas
categorical_transformer = Pipeline(
steps=[('onehot', OneHotEncoder(handle_unknown='ignore', sparse_output=False))]
)
preprocessor = ColumnTransformer(
transformers=[
('numeric', numeric_transformer, numeric_cols),
('cat', categorical_transformer, cat_cols)
],
remainder='passthrough',
verbose_feature_names_out = False
).set_output(transform="pandas")
# Se combinan los pasos de preprocesado y el modelo en un mismo pipeline
pipe = Pipeline([('preprocessing', preprocessor),
('modelo', Ridge())])
# Grid de hiperparámetros
# ==============================================================================
param_grid = {'modelo__alpha': np.logspace(-5, 3, 10)}
# Búsqueda por validación cruzada
# ==============================================================================
grid = GridSearchCV(
estimator = pipe,
param_grid = param_grid,
scoring = 'neg_root_mean_squared_error',
n_jobs = multiprocessing.cpu_count() - 1,
cv = RepeatedKFold(n_splits = 5, n_repeats = 5),
verbose = 0,
return_train_score = True
)
# Se asigna el resultado a _ para que no se imprima por pantalla
_ = grid.fit(X = X_train, y = y_train)
# Resultados del grid
# ==============================================================================
resultados = pd.DataFrame(grid.cv_results_)
resultados.filter(regex = '(param.*|mean_t|std_t)')\
.drop(columns = 'params')\
.sort_values('mean_test_score', ascending = False)
# Mejores hiperparámetros
# ==============================================================================
print("-----------------------------------")
print("Mejores hiperparámetros encontrados")
print("-----------------------------------")
print(f"{grid.best_params_} : {grid.best_score_} ({grid.scoring})")
# Gráfico resultados validación cruzada para cada hiperparámetro
# ==============================================================================
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10, 3.84), sharey=True)
# Gráfico 1
# ------------------------------------------------------------------------------
resultados.plot('param_modelo__alpha', 'mean_train_score', ax=axes[0])
resultados.plot('param_modelo__alpha', 'mean_test_score', ax=axes[0])
axes[0].fill_between(resultados.param_modelo__alpha.astype(float),
resultados['mean_train_score'] + resultados['std_train_score'],
resultados['mean_train_score'] - resultados['std_train_score'],
alpha=0.2)
axes[0].fill_between(resultados.param_modelo__alpha.astype(float),
resultados['mean_test_score'] + resultados['std_test_score'],
resultados['mean_test_score'] - resultados['std_test_score'],
alpha=0.2)
axes[0].legend()
axes[0].set_xscale('log')
axes[0].set_title('Evolución del error CV')
axes[0].set_ylabel('neg_root_mean_squared_error');
# Gráfico 2
# ------------------------------------------------------------------------------
numero_splits = grid.n_splits_
resultados.plot(
x = 'param_modelo__alpha',
y = [f'split{i}_train_score' for i in range(numero_splits)],
alpha = 0.3,
c = 'blue',
ax = axes[1]
)
resultados.plot(
x = 'param_modelo__alpha',
y = [f'split{i}_test_score' for i in range(numero_splits)],
alpha = 0.3,
c = 'red',
ax = axes[1]
)
axes[1].legend(
(axes[1].get_children()[0], axes[1].get_children()[numero_splits]),
('training scores', 'test_scores')
)
axes[1].set_xscale('log')
axes[1].set_title('Evolución del error CV')
axes[1].set_ylabel('neg_root_mean_squared_error');
fig.tight_layout()
fig, ax = plt.subplots(figsize=(6, 3.84))
ax.barh(
[str(d) for d in resultados['params']],
resultados['mean_test_score'],
xerr=resultados['std_test_score'],
align='center',
alpha=0
)
ax.plot(
resultados['mean_test_score'],
[str(d) for d in resultados['params']],
marker="D",
linestyle="",
alpha=0.8,
color="r"
)
ax.set_title('Comparación de Hiperparámetros')
ax.set_xlabel('Test neg_root_mean_squared_error');
Si en GridSearchCV() se indica refit=True, tras identificar los mejores hiperparámetros, se reentrena el modelo con ellos y se almacena en .best_estimator_.
Random grid search¶
GridSearchCV() hace una búsqueda exhaustiva evaluando todas las combinaciones de parámetros. Esta estrategia tiene el inconveniente de que se puede invertir mucho tiempo en regiones de poco interés antes de evaluar otras combinaciones. Una alternativa es hacer una búsqueda aleatoria, de esta forma, se consigue explorar el espacio de búsqueda de una forma más distribuida. RandomizedSearchCV() permite este tipo de estrategia, únicamente requiere que se le indique el espacio de búsqueda de cada hiperparámetro (lista de opciones o una distribución) y el número de combinaciones aleatorias a evaluar.
from itertools import product
import random
fig, axs = plt.subplots(nrows = 1, ncols = 2,figsize=(8, 4),
sharex = True, sharey = False)
# Grid exhaustivo
# ==============================================================================
hyperparametro_1 = np.linspace(start = 0, stop = 100, num=20)
hyperparametro_2 = np.linspace(start = 0, stop = 10, num=20)
# Lista con todas las combinaciones
combinaciones = [list(x) for x in product(hyperparametro_1, hyperparametro_2)]
combinaciones = pd.DataFrame.from_records(
combinaciones,
columns=['hyperparametro_1', 'hyperparametro_2']
)
combinaciones.plot(
x = 'hyperparametro_1',
y = 'hyperparametro_2',
kind = 'scatter',
ax = axs[0]
)
axs[0].set_title('Distribución uniforme')
# Grid aleatorio (random grid)
# ==============================================================================
hyperparametro_1 = np.random.uniform(low = 0, high = 100, size = 400)
hyperparametro_2 = np.random.uniform(low = 0, high = 10, size = 400)
combinaciones = pd.DataFrame(
{
'hyperparametro_1': hyperparametro_1,
'hyperparametro_2': hyperparametro_2,
}
)
combinaciones.plot(
x = 'hyperparametro_1',
y = 'hyperparametro_2',
kind = 'scatter',
ax = axs[1]
)
axs[1].set_title('Distribución aleatoria');
from sklearn.model_selection import RandomizedSearchCV, RepeatedKFold
# Espacio de búsqueda de cada hiperparámetro
# ==============================================================================
param_distributions = {'modelo__alpha': np.logspace(-5, 3, 100)}
# Búsqueda por validación cruzada
# ==============================================================================
grid = RandomizedSearchCV(
estimator = pipe,
param_distributions = param_distributions,
n_iter = 50,
scoring = 'neg_root_mean_squared_error',
n_jobs = multiprocessing.cpu_count() - 1,
cv = RepeatedKFold(n_splits = 5, n_repeats = 5),
verbose = 0,
random_state = 123,
return_train_score = True
)
grid.fit(X = X_train, y = y_train)
# Resultados del grid
# ==============================================================================
resultados = pd.DataFrame(grid.cv_results_)
resultados.filter(regex = '(param.*|mean_t|std_t)')\
.drop(columns = 'params')\
.sort_values('mean_test_score', ascending = False)\
.head(1)
Optimización bayesiana¶
Aunque estas tres estrategias son totalmente válidas y generan buenos resultados, sobretodo cuando se tiene criterio para acotar el rango de búsqueda, comparten una carencia común: ninguna tiene en cuenta los resultados obtenidos hasta el momento, lo que les impide focalizar la búsqueda en las regiones de mayor interés y evitar regiones innecesarias.
Una alternativa es la búsqueda de hiperparámetros con métodos de optimización bayesiana. En términos generales, la optimización bayesiana de hiperparámetros consiste en crear un modelo probabilístico en el que la función objetivo es la métrica de validación del modelo (rmse, auc, precisión..). Con esta estrategia, se consigue que la búsqueda se vaya redirigiendo en cada iteración hacia las regiones de mayor interés. El objetivo final es reducir el número de combinaciones de hiperparámetros con las que se evalúa el modelo, eligiendo únicamente los mejores candidatos. Esto significa que, la ventaja frente a las otras estrategias mencionadas, se maximiza cuando el espacio de búsqueda es muy amplio o la evaluación del modelo es muy lenta.
Dos librerías que permiten aplicar estrategias de optimización Bayesiana para encontrar hiperparámetros en modelos scikitlearn son: scikit-optimize y optuna.
Nota: scikit-optimize ha dejado de ser compatible con versiones de numpy igual o superiores a 1.24.0
# Búsqueda de hiperparámetros con scikit-optimize
# ==============================================================================
from skopt.space import Real, Integer
from skopt.utils import use_named_args
from skopt import gp_minimize
from skopt.plots import plot_convergence
espacio_busqueda = [Real(1e-6, 1e+3, "log-uniform", name='modelo__alpha')]
@use_named_args(espacio_busqueda)
def objective(**params):
pipe.set_params(**params)
return -np.mean(cross_val_score(pipe, X_train, y_train, cv=5, n_jobs=-1,
scoring="neg_root_mean_squared_error"))
resultados_opt = gp_minimize(
func = objective,
dimensions = espacio_busqueda,
n_calls = 50,
random_state = 0
)
print(f"Mejor score validación: {resultados_opt.fun}")
print(f"Mejores hiperparámetros: {list(zip([x.name for x in espacio_busqueda], resultados_opt.x))}")
# Evolución de la optimización
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 3.84))
plot_convergence(resultados_opt, ax = ax);
# Búsqueda de hiperparámetros con Optuna
# ==============================================================================
import optuna
optuna.logging.set_verbosity(optuna.logging.WARNING)
def objective(trial):
modelo__alpha = trial.suggest_float("modelo__alpha", 1e-6, 1e+3, log=True)
pipe.set_params(**{'modelo__alpha':modelo__alpha})
score = -np.mean(cross_val_score(pipe, X_train, y_train, cv=5, n_jobs=-1,
scoring="neg_root_mean_squared_error"))
return score
study = optuna.create_study(direction="minimize")
study.optimize(objective, n_trials=50, show_progress_bar=True)
print(study.best_trial)
print("")
print(f"Mejor score validación: {study.best_value}")
print(f"Mejores hiperparámetros: {study.best_params}")
Tuning del preprocesado¶
En la mayoría de casos, el proceso de optimización se centra en los hiperparámetros del modelo. Sin embargo, también puede ser muy interesante comparar diferentes transformaciones de preprocesado. Gracias a los pipeline, esto puede hacerse de la misma forma que con los hiperparametros. Véase el siguiente ejemplo en el que se compara el incorporar interacciones entre predictores.
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import OneHotEncoder
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import PolynomialFeatures
from sklearn.compose import make_column_selector
# Pipe: preprocesado + modelo
# ==============================================================================
# Transformaciones para las variables numéricas
numeric_transformer = Pipeline(
steps=[('scaler', StandardScaler())]
)
# Transformaciones para las variables categóricas
categorical_transformer = Pipeline(
steps=[('onehot', OneHotEncoder(handle_unknown='ignore', sparse_output=False))]
)
preprocessor = ColumnTransformer(
transformers=[
('numeric', numeric_transformer, numeric_cols),
('cat', categorical_transformer, cat_cols)
],
remainder='passthrough',
verbose_feature_names_out = False
).set_output(transform="pandas")
# Se combinan los pasos de preprocesado y el modelo en un mismo pipeline.
pipe = Pipeline(
[('preprocessing', preprocessor),
('interactions', PolynomialFeatures(degree=2)),
('modelo', Ridge())])
# Grid de hiperparámetros
# ==============================================================================
param_grid = {'interactions': [PolynomialFeatures(degree=2), 'passthrough'],
'modelo__alpha': np.logspace(-5, 3, 10)}
# Búsqueda por validación cruzada
# ==============================================================================
grid = GridSearchCV(
estimator = pipe,
param_grid = param_grid,
scoring = 'neg_root_mean_squared_error',
n_jobs = multiprocessing.cpu_count() - 1,
cv = RepeatedKFold(n_splits = 5, n_repeats = 5),
verbose = 0,
return_train_score = True
)
grid.fit(X = X_train, y = y_train)
# Resultados del grid
# ==============================================================================
resultados = pd.DataFrame(grid.cv_results_)
resultados.filter(regex = '(param.*|mean_t|std_t)')\
.drop(columns = 'params')\
.sort_values('mean_test_score', ascending = False)\
.head(1)
Algoritmos¶
En los siguientes apartados se entrenan diferentes modelos de machine learning con el objetivo de compararlos e identificar el que mejor resultado obtiene prediciendo el precio de las viviendas.
K-Nearest Neighbor (kNN)¶
K-Nearest Neighbor es uno de los algoritmos de machine learning más simples. Su funcionamiento es el siguiente: para predecir una observación se identifican las K observaciones del conjunto de entrenamiento que más se asemejan a ella (en base a sus predictores) y se emplea como valor predicho el promedio de la variable respuesta en dichas observaciones. Dada su sencillez, suele dar peores resultados que otros algoritmos, pero es un buen referente como baseline.
from sklearn.model_selection import RandomizedSearchCV, RepeatedKFold
from sklearn.neighbors import KNeighborsRegressor
# Pipeline: preprocesado + modelo
# ==============================================================================
# Identificación de columnas numéricas y catégoricas
numeric_cols = X_train.select_dtypes(include=['float64', 'int']).columns.to_list()
cat_cols = X_train.select_dtypes(include=['object', 'category']).columns.to_list()
# Transformaciones para las variables numéricas
numeric_transformer = Pipeline(
steps=[('scaler', StandardScaler())]
)
# Transformaciones para las variables categóricas
categorical_transformer = Pipeline(
steps=[('onehot', OneHotEncoder(handle_unknown='ignore', sparse_output=False))]
)
preprocessor = ColumnTransformer(
transformers=[
('numeric', numeric_transformer, numeric_cols),
('cat', categorical_transformer, cat_cols)
],
remainder='passthrough',
verbose_feature_names_out = False
).set_output(transform="pandas")
# Se combinan los pasos de preprocesado y el modelo en un mismo pipeline.
pipe = Pipeline([('preprocessing', preprocessor),
('modelo', KNeighborsRegressor())])
# Optimización de hiperparámetros
# ==============================================================================
# Espacio de búsqueda de cada hiperparámetro
param_distributions = {'modelo__n_neighbors': np.linspace(1, 100, 500, dtype=int)}
# Búsqueda random grid
grid = RandomizedSearchCV(
estimator = pipe,
param_distributions = param_distributions,
n_iter = 20,
scoring = 'neg_root_mean_squared_error',
n_jobs = multiprocessing.cpu_count() - 1,
cv = RepeatedKFold(n_splits = 5, n_repeats = 3),
refit = True,
verbose = 0,
random_state = 123,
return_train_score = True
)
grid.fit(X = X_train, y = y_train)
# Resultados del grid
# ==============================================================================
resultados = pd.DataFrame(grid.cv_results_)
resultados.filter(regex = '(param.*|mean_t|std_t)')\
.drop(columns = 'params')\
.sort_values('mean_test_score', ascending = False)\
.head(1)
# Gráfico resultados validación cruzada para cada hiperparámetro
# ------------------------------------------------------------------------------
fig, ax = plt.subplots(figsize=(6, 3.84))
hiperparametro = 'param_modelo__n_neighbors'
resultados = resultados.sort_values(hiperparametro, ascending = False)
metrica = grid.scoring
resultados.plot(hiperparametro, 'mean_train_score', ax=ax)
resultados.plot(hiperparametro, 'mean_test_score', ax=ax)
ax.fill_between(resultados[hiperparametro].astype(int),
resultados['mean_train_score'] + resultados['std_train_score'],
resultados['mean_train_score'] - resultados['std_train_score'],
alpha=0.2)
ax.fill_between(resultados[hiperparametro].astype(int),
resultados['mean_test_score'] + resultados['std_test_score'],
resultados['mean_test_score'] - resultados['std_test_score'],
alpha=0.2)
ax.legend()
ax.set_title('Evolución del error CV')
ax.set_ylabel(metrica);
fig, ax = plt.subplots(figsize=(8, 5))
resultados = resultados.sort_values('mean_test_score', ascending = True)
ax.barh(
[str(d) for d in resultados['params']],
resultados['mean_test_score'],
xerr=resultados['std_test_score'],
align='center',
alpha=0
)
ax.plot(
resultados['mean_test_score'],
[str(d) for d in resultados['params']],
marker="D",
linestyle="",
alpha=0.8,
color="r"
)
ax.set_title('Comparación de Hiperparámetros')
ax.set_ylabel(metrica);
Una vez identificados los mejores hiperparámetros, se reentrena el modelo indicando los valores óptimos en sus argumentos. Si en el GridSearchCV() se indica refit=True, este reentrenamiento se hace automáticamente y
el modelo resultante se encuentra almacenado en .best_estimator_.
# Error de test del modelo final
# ==============================================================================
modelo_final = grid.best_estimator_
predicciones = modelo_final.predict(X = X_test)
rmse_knn = mean_squared_error(
y_true = y_test,
y_pred = predicciones,
squared = False
)
print(f"El error (rmse) de test es: {rmse_knn}")
Regresión lineal (Ridge y Lasso)¶
La regresión lineal es un método estadístico que trata de modelar la relación lineal entre una variable continua (variable dependiente, respuesta o dependiente) y una o más variables independientes (regresores o predictores) mediante el ajuste de una ecuación lineal. Se llama regresión lineal simple cuando solo hay una variable independiente y regresión lineal múltiple cuando hay más de una.
El modelo de regresión lineal (Legendre, Gauss, Galton y Pearson) considera que, dado un conjunto de observaciones, la media $\mu$ de la variable respuesta $Y$ se relaciona de forma lineal con la o las variables regresoras $X$ acorde a la ecuación:
$$\mu_Y = \beta_0 + \beta_1 x_{1} + \beta_2 x_{2} + ... + \beta_p x_{p}$$o en notación matricial (incorporando $\beta_0$ en el vector $\mathbf{\beta}$):
$$\mu_Y = \mathbf{X}^T \mathbf{\beta}$$Además, el modelo lineal puede incluir regularización durante su ajuste, por lo que también incluye los modelos ridge regression, lasso y elastic net.
Scikit-Learn incorpora 3 tipos de regularización para los modelos lineales con el objetivo de evitar overfitting, reducir varianza y atenuar el efecto de la correlación entre predictores. Por lo general, aplicando regularización se consigue modelos con mayor poder predictivo (generalización).
El modelo lasso es un modelo lineal por mínimos cuadrados que incorpora una regularización que penaliza la suma del valor absolutos de los coeficientes de regresión $(||\beta||_1 = \sum_{k=1} ^p |\beta_k|)$. A esta penalización se le conoce como l1 y tiene el efecto de forzar a que los coeficientes de los predictores tiendan a cero. Dado que un predictor con coeficiente de regresión cero no influye en el modelo, lasso consigue seleccionar los predictores más influyentes. El grado de penalización está controlado por el hiperparámetro $\lambda$. Cuando $\lambda = 0$, el resultado es equivalente al de un modelo lineal por mínimos cuadrados ordinarios. A medida que $\lambda$ aumenta, mayor es la penalización y más predictores quedan excluidos.
El modelo ridge es un modelo lineal por mínimos cuadrados que incorpora una regularización que penaliza la suma de los coeficientes elevados al cuadrado $(||\beta||^2_2 = \sum_{k=1} ^p \beta^2_k)$. A esta penalización se le conoce como l2 y tiene el efecto de reducir de forma proporcional el valor de todos los coeficientes del modelo pero sin que estos lleguen a cero. Al igual que lasso, el grado de penalización está controlado por el hiperparámetro $\lambda$.
La principal diferencia práctica entre lasso y ridge es que el primero consigue que algunos coeficientes sean exactamente cero, por lo que realiza selección de predictores, mientras que el segundo no llega a excluir ninguno. Esto supone una ventaja notable de lasso en escenarios donde no todos los predictores son importantes para el modelo y se desea que los menos influyentes queden excluidos. Por otro lado, cuando existen predictores altamente correlacionados (linealmente), ridge reduce la influencia de todos ellos a la vez y de forma proporcional, mientras que lasso tiende a seleccionar uno de ellos, dándole todo el peso y excluyendo al resto. En presencia de correlaciones, esta selección varía mucho con pequeñas perturbaciones (cambios en los datos de entrenamiento), por lo que, las soluciones de lasso, son muy inestables si los predictores están altamente correlacionados.
Para conseguir un equilibrio óptimo entre estas dos propiedades, se puede emplear lo que se conoce como penalización elastic net, que combina ambas estrategias.
El modelo elastic net incluye una regularización que combina la penalización l1 y l2 $(\alpha \lambda ||\beta||_1 + \frac{1}{2}(1- \alpha)||\beta||^2_2)$. El grado en que influye cada una de las penalizaciones está controlado por el hiperparámetro $\alpha$. Su valor debe estar comprendido en el intervalo [0,1], cuando $\alpha = 0$, se aplica ridge regression y cuando $\alpha = 1$ se aplica lasso. La combinación de ambas penalizaciones suele dar lugar a buenos resultados. Una estrategia frecuentemente utilizada es asignarle casi todo el peso a la penalización l1 ($\alpha$ muy próximo a 1) para conseguir seleccionar predictores y un poco a la l2 para dar cierta estabilidad en el caso de que algunos predictores estén correlacionados.
Encontrar el mejor modelo implica identificar los valores óptimos de los hiperparámetros de regularización $\alpha$ y $\lambda$.
from sklearn.model_selection import RandomizedSearchCV, RepeatedKFold
from sklearn.linear_model import Ridge
# Pipeline: preprocesado + modelo
# ==============================================================================
# Identificación de columnas numéricas y catégoricas
numeric_cols = X_train.select_dtypes(include=['float64', 'int']).columns.to_list()
cat_cols = X_train.select_dtypes(include=['object', 'category']).columns.to_list()
# Transformaciones para las variables numéricas
numeric_transformer = Pipeline(
steps=[('scaler', StandardScaler())]
)
# Transformaciones para las variables categóricas
categorical_transformer = Pipeline(
steps=[('onehot', OneHotEncoder(handle_unknown='ignore', sparse_output=False))]
)
preprocessor = ColumnTransformer(
transformers=[
('numeric', numeric_transformer, numeric_cols),
('cat', categorical_transformer, cat_cols)
],
remainder='passthrough',
verbose_feature_names_out = False
).set_output(transform="pandas")
# Se combinan los pasos de preprocesado y el modelo en un mismo pipeline.
pipe = Pipeline([('preprocessing', preprocessor),
('modelo', Ridge())])
# Optimización de hiperparámetros
# ==============================================================================
# Espacio de búsqueda de cada hiperparámetro
param_distributions = {'modelo__alpha': np.logspace(-5, 5, 500)}
# Búsqueda random grid
grid = RandomizedSearchCV(
estimator = pipe,
param_distributions = param_distributions,
n_iter = 20,
scoring = 'neg_root_mean_squared_error',
n_jobs = multiprocessing.cpu_count() - 1,
cv = RepeatedKFold(n_splits = 5, n_repeats = 3),
refit = True,
verbose = 0,
random_state = 123,
return_train_score = True
)
grid.fit(X = X_train, y = y_train)
# Resultados del grid
# ==============================================================================
resultados = pd.DataFrame(grid.cv_results_)
resultados.filter(regex = '(param.*|mean_t|std_t)')\
.drop(columns = 'params')\
.sort_values('mean_test_score', ascending = False)\
.head(1)
# Gráfico resultados validación cruzada para cada hiperparámetro
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 3.84))
hiperparametro = 'param_modelo__alpha'
resultados = resultados.sort_values(hiperparametro, ascending = False)
metrica = grid.scoring
resultados.plot(hiperparametro, 'mean_train_score', ax=ax)
resultados.plot(hiperparametro, 'mean_test_score', ax=ax)
ax.fill_between(resultados[hiperparametro].astype(int),
resultados['mean_train_score'] + resultados['std_train_score'],
resultados['mean_train_score'] - resultados['std_train_score'],
alpha=0.2)
ax.fill_between(resultados[hiperparametro].astype(int),
resultados['mean_test_score'] + resultados['std_test_score'],
resultados['mean_test_score'] - resultados['std_test_score'],
alpha=0.2)
ax.legend()
ax.set_title('Evolución del error CV')
ax.set_ylabel(metrica);
# Error de test del modelo final
# ==============================================================================
modelo_final = grid.best_estimator_
predicciones = modelo_final.predict(X = X_test)
rmse_lm = mean_squared_error(
y_true = y_test,
y_pred = predicciones,
squared = False
)
print(f"El error (rmse) de test es: {rmse_lm}")
Random Forest¶
Un modelo Random Forest está formado por un conjunto de árboles de decisión individuales, cada uno entrenado con una muestra ligeramente distinta de los datos de entrenamiento generada mediante bootstrapping). La predicción de una nueva observación se obtiene agregando las predicciones de todos los árboles individuales que forman el modelo. Para conocer más sobre este tipo de modelo visitar Random Forest con Python.
from sklearn.model_selection import RandomizedSearchCV, RepeatedKFold
from sklearn.ensemble import RandomForestRegressor
# Pipeline: preprocesado + modelo
# ==============================================================================
# Identificación de columnas numéricas y catégoricas
numeric_cols = X_train.select_dtypes(include=['float64', 'int']).columns.to_list()
cat_cols = X_train.select_dtypes(include=['object', 'category']).columns.to_list()
# Transformaciones para las variables numéricas
numeric_transformer = Pipeline(
steps=[('scaler', StandardScaler())]
)
# Transformaciones para las variables categóricas
categorical_transformer = Pipeline(
steps=[('onehot', OneHotEncoder(handle_unknown='ignore', sparse_output=False))]
)
preprocessor = ColumnTransformer(
transformers=[
('numeric', numeric_transformer, numeric_cols),
('cat', categorical_transformer, cat_cols)
],
remainder='passthrough',
verbose_feature_names_out = False
).set_output(transform="pandas")
# Se combinan los pasos de preprocesado y el modelo en un mismo pipeline.
pipe = Pipeline([('preprocessing', preprocessor),
('modelo', RandomForestRegressor())])
# Optimización de hiperparámetros
# ==============================================================================
# Espacio de búsqueda de cada hiperparámetro
param_distributions = {
'modelo__n_estimators': [50, 100, 1000, 2000],
'modelo__max_features': [3, 5, 7, 1.0],
'modelo__max_depth' : [None, 3, 5, 10, 20]
}
# Búsqueda random grid
grid = RandomizedSearchCV(
estimator = pipe,
param_distributions = param_distributions,
n_iter = 20,
scoring = 'neg_root_mean_squared_error',
n_jobs = multiprocessing.cpu_count() - 1,
cv = RepeatedKFold(n_splits = 5, n_repeats = 3),
refit = True,
verbose = 0,
random_state = 123,
return_train_score = True
)
grid.fit(X = X_train, y = y_train)
# Resultados del grid
# ==============================================================================
resultados = pd.DataFrame(grid.cv_results_)
resultados.filter(regex = '(param.*|mean_t|std_t)')\
.drop(columns = 'params')\
.sort_values('mean_test_score', ascending = False)\
.head(1)
# Error de test del modelo final
# ==============================================================================
modelo_final = grid.best_estimator_
predicciones = modelo_final.predict(X = X_test)
rmse_rf = mean_squared_error(
y_true = y_test,
y_pred = predicciones,
squared = False
)
print(f"El error (rmse) de test es: {rmse_rf}")
Gradient Boosting Trees¶
Un modelo Gradient Boosting Trees está formado por un conjunto de árboles de decisión individuales, entrenados de forma secuencial, de forma que cada nuevo árbol trata de mejorar los errores de los árboles anteriores. La predicción de una nueva observación se obtiene agregando las predicciones de todos los árboles individuales que forman el modelo. Para conocer más sobre este tipo de modelo visitar Gradient Boosting con Python.
from sklearn.model_selection import RandomizedSearchCV, RepeatedKFold
from sklearn.ensemble import GradientBoostingRegressor
# Pipeline: preprocesado + modelo
# ==============================================================================
# Identificación de columnas numéricas y catégoricas
numeric_cols = X_train.select_dtypes(include=['float64', 'int']).columns.to_list()
cat_cols = X_train.select_dtypes(include=['object', 'category']).columns.to_list()
# Transformaciones para las variables numéricas
numeric_transformer = Pipeline(
steps=[('scaler', StandardScaler())]
)
# Transformaciones para las variables categóricas
categorical_transformer = Pipeline(
steps=[('onehot', OneHotEncoder(handle_unknown='ignore', sparse_output=False))]
)
preprocessor = ColumnTransformer(
transformers=[
('numeric', numeric_transformer, numeric_cols),
('cat', categorical_transformer, cat_cols)
],
remainder='passthrough',
verbose_feature_names_out = False
).set_output(transform="pandas")
# Se combinan los pasos de preprocesado y el modelo en un mismo pipeline.
pipe = Pipeline([('preprocessing', preprocessor),
('modelo', GradientBoostingRegressor())])
# Optimización de hiperparámetros
# ==============================================================================
# Espacio de búsqueda de cada hiperparámetro
param_distributions = {
'modelo__n_estimators': [50, 100, 1000, 2000],
'modelo__max_features': [3, 5, 7, 1.0],
'modelo__max_depth' : [None, 3, 5, 10, 20],
'modelo__subsample' : [0.5,0.7, 1]
}
# Búsqueda random grid
grid = RandomizedSearchCV(
estimator = pipe,
param_distributions = param_distributions,
n_iter = 20,
scoring = 'neg_root_mean_squared_error',
n_jobs = multiprocessing.cpu_count() - 1,
cv = RepeatedKFold(n_splits = 5, n_repeats = 3),
refit = True,
verbose = 0,
random_state = 123,
return_train_score = True
)
grid.fit(X = X_train, y = y_train)
# Resultados del grid
# ==============================================================================
resultados = pd.DataFrame(grid.cv_results_)
resultados.filter(regex = '(param.*|mean_t|std_t)')\
.drop(columns = 'params')\
.sort_values('mean_test_score', ascending = False)\
.head(1)
# Error de test del modelo final
# ==============================================================================
modelo_final = grid.best_estimator_
predicciones = modelo_final.predict(X = X_test)
rmse_gbm = mean_squared_error(
y_true = y_test,
y_pred = predicciones,
squared = False
)
print(f"El error (rmse) de test es: {rmse_gbm}")
Stacking¶
A modo general, el término model ensembling hace referencia a la combinación de las predicciones de dos o más modelos distintos con el objetivo de mejorar las predicciones finales. Esta estrategia se basa en la asunción de que, distintos modelos entrenados independientemente, emplean distintos aspectos de los datos para realizar las predicciones, es decir, cada uno es capaz de identificar parte de la “verdad” pero no toda ella. Combinando la perspectiva de cada uno de ellos, se obtiene una descripción más detallada de la verdadera estructura subyacente en los datos. A modo de analogía, imagínese un grupo de estudiantes que se enfrentan a un examen multidisciplinal. Aunque todos obtengan aproximadamente la misma nota, cada uno de ellos habrá conseguido más puntos con las preguntas que tratan sobre las disciplinas en las que destacan. Si en lugar de hacer el examen de forma individual, lo hacen en grupo, cada uno podrá contribuir en los aspectos que más domina, y el resultado final será probablemente superior a cualquiera de los resultados individuales.
La clave para que el ensembling consiga mejorar los resultados es la diversidad de los modelos. Si todos los modelos combinados son similares entre ellos, no podrán compensarse unos a otros. Por esta razón, se tiene que intentar combinar modelos que sean lo mejor posible a nivel individual y lo más diferentes entre ellos.
Las formas más simples de combinar las predicciones de varios modelos es emplear la media para problemas de regresión y la moda para problemas de clasificación. Sin embargo existen otras aproximaciones más complejas capaces de conseguir mejores resultados:
- Ponderar las agregaciones dando distinto peso a cada modelo, por ejemplo, en proporción al accuracy que han obtenido de forma individual.
- Super Learner (stacked regression): emplear un nuevo modelo para que decida la mejor forma de combinar las predicciones de los otros modelos.
Algoritmo Super Learner¶
La implementación de Super Learner disponible en scikit-learn en las clases StackingRegressor y StackingClassifier sigue el siguiente algoritmo:
Definición del ensemble
Definir un listado con los algoritmos base (cada uno con los hiperparámetros pertinentes).
Seleccionar el algoritmo de metalearning que defina cómo se entrena en modelo superior. Por defecto, se emplea
RidgeCVpara regresión yLogisticRegressionpara clasificación.
Entrenamiento del ensemble
Entrenar cada uno de los algoritmos base con el conjunto de entrenamiento.
Realizar k-fold cross-validation con cada uno de los algoritmos base y almacenar las predicciones hechas en cada una de las k particiones.
Combinar las predicciones del paso 2 en una única matriz NxL (N = número de observaciones en el conjunto de entrenamiento, L = número de modelos base).
Entrenar el metalearning con la variable respuesta y la matriz NxL como predictores.
El Super learner final está formado por los modelos base y el modelo metalearning.
Predecir
Predecir la nueva observación con cada uno de los modelos base.
Emplear las predicciones de los modelos base como input del metalearner para obtener la predicción final.
Se procede a crear un stacking con los modelos Ridge y RandomForest, empleando en cada caso lo mejores hiperparámetros encontrados en los apartados anteriores.
from sklearn.linear_model import Ridge
from sklearn.linear_model import RidgeCV
from sklearn.ensemble import RandomForestRegressor
from sklearn.ensemble import StackingRegressor
# Pipeline: preprocesado + modelos para el stacking
# ==============================================================================
# Identificación de columnas numéricas y catégoricas
numeric_cols = X_train.select_dtypes(include=['float64', 'int']).columns.to_list()
cat_cols = X_train.select_dtypes(include=['object', 'category']).columns.to_list()
# Transformaciones para las variables numéricas
numeric_transformer = Pipeline(
steps=[('scaler', StandardScaler())]
)
# Transformaciones para las variables categóricas
categorical_transformer = Pipeline(
steps=[('onehot', OneHotEncoder(handle_unknown='ignore', sparse_output=False))]
)
preprocessor = ColumnTransformer(
transformers=[
('numeric', numeric_transformer, numeric_cols),
('cat', categorical_transformer, cat_cols)
],
remainder='passthrough',
verbose_feature_names_out = False
).set_output(transform="pandas")
# Se combinan los pasos de preprocesado y los modelos creando varios pipeline.
pipe_ridge = Pipeline([('preprocessing', preprocessor),
('ridge', Ridge(alpha=3.4))])
pipe_rf = Pipeline([('preprocessing', preprocessor),
('random_forest', RandomForestRegressor(
n_estimators = 1000,
max_features = 7,
max_depth = 20
)
)])
# Definición y entrenamiento del StackingRegressor
# ==============================================================================
estimators = [('ridge', pipe_ridge),
('random_forest', pipe_rf)]
stacking_regressor = StackingRegressor(estimators=estimators,
final_estimator=RidgeCV())
# Se asigna el resultado a _ para que no se imprima por pantalla
_ = stacking_regressor.fit(X = X_train, y = y_train)
# Error de test del stacking
# ==============================================================================
modelo_final = stacking_regressor
predicciones = modelo_final.predict(X = X_test)
rmse_stacking = mean_squared_error(
y_true = y_test,
y_pred = predicciones,
squared = False
)
print(f"El error (rmse) de test es: {rmse_stacking}")
error_modelos = pd.DataFrame({
'modelo': ['knn', 'lm', 'random forest', 'gradient boosting',
'stacking'],
'rmse': [rmse_knn, rmse_lm, rmse_rf, rmse_gbm, rmse_stacking]
})
error_modelos = error_modelos.sort_values('rmse', ascending=False)
fig, ax = plt.subplots(figsize=(6, 3.84))
ax.hlines(error_modelos.modelo, xmin=0, xmax=error_modelos.rmse)
ax.plot(error_modelos.rmse, error_modelos.modelo, "o", color='black')
ax.tick_params(axis='y', which='major', labelsize=12)
ax.set_title('Comparación de error de test modelos'),
ax.set_xlabel('Test rmse');
Anexos¶
Anexo1: Métodos de validación¶
Los métodos de validación, también conocidos como resampling, son estrategias que permiten estimar la capacidad predictiva de los modelos cuando se aplican a nuevas observaciones, haciendo uso únicamente de los datos de entrenamiento. La idea en la que se basan todos ellos es la siguiente: el modelo se ajusta empleando un subconjunto de observaciones del conjunto de entrenamiento y se evalúa (calcular una métrica que mida cómo de bueno es el modelo, por ejemplo, accuracy) con las observaciones restantes. Este proceso se repite múltiples veces y los resultados se agregan y promedian. Gracias a las repeticiones, se compensan las posibles desviaciones que puedan surgir por el reparto aleatorio de las observaciones. La diferencia entre métodos suele ser la forma en la que se generan los subconjuntos de entrenamiento/validación.
k-Fold-Cross-Validation (CV)
Las observaciones de entrenamiento se reparten en k folds (conjuntos) del mismo tamaño. El modelo se ajusta con todas las observaciones excepto las del primer fold y se evalúa prediciendo las observaciones del fold que ha quedado excluido, obteniendo así la primera métrica. El proceso se repite k veces, excluyendo un fold distinto en cada iteración. Al final, se generan k valores de la métrica, que se agregan (normalmente con la media y la desviación típica) generando la estimación final de validación.
Leave-One-Out Cross-Validation (LOOCV)
LOOCV es un caso especial de k-Fold-Cross-Validation en el que el número k de folds es igual al número de observaciones disponibles en el conjunto de entrenamiento. El modelo se ajusta cada vez con todas las observaciones excepto una, que se emplea para evaluar el modelo. Este método supone un coste computacional muy elevado, el modelo se ajusta tantas veces como observaciones de entrenamiento, por lo que, en la práctica, no suele compensar emplearlo.
Repeated k-Fold-Cross-Validation (repeated CV)
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.
Leave-Group-Out Cross-Validation (LGOCV)
LGOCV, también conocido como repeated train/test splits o Monte Carlo Cross-Validation, consiste simplemente en generar múltiples divisiones aleatorias entrenamiento-test (solo dos conjuntos por repetición). La proporción de observaciones que va a cada conjunto se determina de antemano, 80%-20% suele dar buenos resultados. Este método, aunque más simple de implementar que CV, requiere de muchas repeticiones (>50) para generar estimaciones estables.
Bootstrapping
Una muestra bootstrap es una muestra obtenida a partir de la muestra original por muestreo aleatorio con reposición, y del mismo tamaño que la muestra original. Muestreo aleatorio con reposición (resampling with replacement) significa que, después de que una observación sea extraída, se vuelve a poner a disposición para las siguientes extracciones. Como resultado de este tipo de muestreo, algunas observaciones aparecerán múltiples veces en la muestra bootstrap y otras ninguna. Las observaciones no seleccionadas reciben el nombre de out-of-bag (OOB). Por cada iteración de bootstrapping se genera una nueva muestra bootstrap, se ajusta el modelo con ella y se evalúa con las observaciones out-of-bag.
Obtener una nueva muestra del mismo tamaño que la muestra original mediante muestreo 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.
La naturaleza del proceso de bootstrapping genera cierto bias en las estimaciones que puede ser problemático cuando el conjunto de entrenamiento es pequeño. Existen ciertas modificaciones del algoritmo original para corregir este problema, algunos de ellos son: 632 method y 632+ method.
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.
Puede encontrarse un estudio comparativo de los diferentes métodos en Comparing Different Species of Cross-Validation.
Anexo2: Métricas¶
Existe una gran variedad de métricas que permiten evaluar cómo de bueno es un algoritmo realizando predicciones. La idoneidad de cada una depende completamente del problema en cuestión, y su correcta elección dependerá de cómo de bien entienda el analista el problema al que se enfrenta. A continuación, se describen algunas de las más utilizadas.
Accuracy y Kappa
Estas dos métricas son las más empleadas en problemas de clasificación binaria y multiclase. Accuracy es el porcentaje de observaciones correctamente clasificadas respecto al total de predicciones. Kappa o Cohen’s Kappa es el valor de accuracy normalizado respecto del porcentaje de acierto esperado por azar. A diferencia de accuracy, cuyo rango de valores puede ser [0, 1], el de kappa es [-1, 1]. En problemas con clases desbalanceadas, donde el grupo mayoritario supera por mucho a los otros, Kappa es más útil porque evita caer en la ilusión de creer que un modelo es bueno cuando realmente solo supera por poco lo esperado por azar.
MSE, RMSE, MAE
Estas son las métricas más empleadas en problemas de regresión.
MSE (Mean Squared Error) es la media de los errores elevados al cuadrado. Suele ser muy buen indicativo de cómo funciona el modelo en general, pero tiene la desventaja de estar en unidades cuadradas. Para mejorar la interpretación, suele emplearse RMSE (Root Mean Squared Error), que es la raíz cuadrada del MSE y por lo tanto sus unidades son las mismas que la variable respuesta.
MAE (Mean Absolute Error) es la media de los errores en valor absoluto. La diferencia respecto a MSE es que, este último, eleva al cuadrado los errores, lo que significa que penaliza mucho más las desviaciones grandes. A modo general, MSE favorece modelos que se comportan aproximadamente igual de bien en todas las observaciones, mientras que MAE favorece modelos que predicen muy bien la gran mayoría de observaciones aunque en unas pocas se equivoque por mucho.
Información de sesión¶
import session_info
session_info.show(html=False)
Bibliografía¶
Scikit-learn: Machine Learning in Python, Pedregosa et al., JMLR 12, pp. 2825-2830, 2011.
API design for machine learning software: experiences from the scikit-learn project, Buitinck et al., 2013.
Uso de pipes y ColumnTransformer: https://www.dataschool.io/encoding-categorical-features-in-python
Introduction to Statistical Learning, Gareth James, Daniela Witten, Trevor Hastie and Robert Tibshirani
Applied Predictive Modeling Max Kuhn, Kjell Johnson
T.Hastie, R.Tibshirani, J.Friedman. The Elements of Statistical Learning. Springer, 2009
¿Cómo citar este documento?
Machine learning con Python y Scikit-learn por Joaquín Amat Rodrigo, disponible con licencia CC BY-NC-SA 4.0 en https://www.cienciadedatos.net/documentos/py06_machine_learning_python_scikitlearn.html 
¿Te ha gustado el artículo? Tu ayuda es importante
Mantener un sitio web tiene unos costes elevados, tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊
Este contenido, creado por Joaquín Amat Rodrigo, tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.
Se permite:
-
Compartir: copiar y redistribuir el material en cualquier medio o formato.
-
Adaptar: remezclar, transformar y crear a partir del material.
Bajo los siguientes términos:
-
Atribución: Debes otorgar el crédito adecuado, proporcionar un enlace a la licencia e indicar si se realizaron cambios. Puedes hacerlo de cualquier manera razonable, pero no de una forma que sugiera que el licenciante te respalda o respalda tu uso.
-
NoComercial: No puedes utilizar el material para fines comerciales.
-
CompartirIgual: Si remezclas, transformas o creas a partir del material, debes distribuir tus contribuciones bajo la misma licencia que el original.