1. Introducción

El presente trabajo se basa en el conjunto de datos publicado por el usuario Vipul L Rathod en Kaggle: Fish Market.

Este dataset recoge información de diferentes especies de peces a través de medidas morfológicas y su peso.

Las variables incluidas son:

• Species: especie del pez (string).
• Weight: peso en gramos (float).
• Length1, Length2, Length3: longitudes en cm medidas en diferentes ejes (vertical, diagonal y transversal) (float).
• Height: altura en cm (float).
• Width: ancho en cm (float).

En este trabajo se plantea como objetivo predecir el peso (Weight) a partir del resto de variables. Si bien es posible utilizar el dataset para clasificar la especie, en este caso se prioriza la estimación del peso porque puede aportar un valor práctico en la industria pesquera.

Por ejemplo, una empresa podría automatizar la estimación del peso de los peces mediante la captura de imágenes en una cinta transportadora, extrayendo las medidas morfológicas y especie con técnicas de procesamiento de imagen. De este modo, se evitaría el proceso de pesado, reduciendo costes y tiempo. De igual manera, se podría estimar el beneficio —conociendo el coste del pescado por kilo—, permitiendo replantear la logística y la distribución de venta.

2. Carga, Exploración y Preprocesamiento

a) Carga

import pandas as pd

file = 'Fish.csv'
df = pd.read_csv(file)
print("Filas:", df.shape[0], " |  Columnas:", df.shape[1])
df.head()

Filas: 159  |  Columnas: 7

	Species	Weight	Length1	Length2	Length3	Height	Width
0	Bream	242.0	23.2	25.4	30.0	11.5200	4.0200
1	Bream	290.0	24.0	26.3	31.2	12.4800	4.3056
2	Bream	340.0	23.9	26.5	31.1	12.3778	4.6961
3	Bream	363.0	26.3	29.0	33.5	12.7300	4.4555
4	Bream	430.0	26.5	29.0	34.0	12.4440	5.1340

b) Consistencia

Conociendo que se registran medidas y peso será sencillo comprobar su consistencia.

# ¿Hay algún valor negativo o 0? Es decir, sin sentido
numeric_cols = df.select_dtypes(include=['number'])
no_sense = (numeric_cols <= 0).any(axis=1)
print(df[no_sense])
print("\nHay 1 valor sin sentido (peso 0), se elimina:")
print(df.shape)
df = df[~no_sense]
print(df.shape)

   Species  Weight  Length1  Length2  Length3  Height   Width
40   Roach     0.0     19.0     20.5     22.8  6.4752  3.3516

Hay 1 valor sin sentido (peso 0), se elimina:
(159, 7)
(158, 7)

Con el mismo propósito, se comprueban las especies y nulos.

print("Valores únicos de Species:", df['Species'].unique())
print("\nNulos por columna:")
df.isnull().sum()

Valores únicos de Species: ['Bream' 'Roach' 'Whitefish' 'Parkki' 'Perch' 'Pike' 'Smelt']

Nulos por columna:

	0
Species	0
Weight	0
Length1	0
Length2	0
Length3	0
Height	0
Width	0

dtype: int64

c) Exploración

Véanse en primer lugar, las distribuciones de variables numéricas.

import matplotlib.pyplot as plt

print(df.describe())
numeric_cols = df.select_dtypes(include=['number'])
numeric_cols.hist(bins=15, layout=(2, 3), edgecolor='black', grid=False, color='skyblue', alpha=0.7)
plt.show()

            Weight     Length1     Length2     Length3      Height       Width
count   158.000000  158.000000  158.000000  158.000000  158.000000  158.000000
mean    400.847468   26.293038   28.465823   31.280380    8.986790    4.424232
std     357.697796   10.011427   10.731707   11.627605    4.295191    1.689010
min       5.900000    7.500000    8.400000    8.800000    1.728400    1.047600
25%     121.250000   19.150000   21.000000   23.200000    5.940600    3.398650
50%     281.500000   25.300000   27.400000   29.700000    7.789000    4.277050
75%     650.000000   32.700000   35.750000   39.675000   12.371850    5.586750
max    1650.000000   59.000000   63.400000   68.000000   18.957000    8.142000

Observando lo siguiente:

• Width es la única similar a una normal.
• El resto siguen una distribución asimétrica a la derecha, en especial Weight.
• Las variables Length_n tienen distribuciones muy similares.
• Como era de esperar, las escalas son muy distintas (ver ejes X o print).

Para conocer la relación del conjunto numérico con la variable objetivo se propone el siguiente gráfico.

import seaborn as sns

sns.pairplot(df, hue='Weight')
plt.show()

Gracias a este gráfico, se aprecia como las variables siguen una tendencia creciente en peso a mayor es su valor. Mientras, la variable Height también sigue esta tendencia pero de forma menos marcada y mucho menos lineal.

Otra opción recomendable es realizar el mismo gráfico empleando hue=Species, sin embargo, para no extender los outputs se evita. Con todo, este muestra la pertenencia de especies en pesos muy concretos y otras que se extienden en todo el rango, por ejemplo: "Smelt" tiene un peso siempre muy bajo, mientras que "Perch" tiene un rango muy extendido.

Por otro lado, ambos gráficos muestran una evidente multicolinealidad entre las variables referentes a Lenght. Para demostrar este hecho, véase la siguiente celda.

Nota: se evita este cálculo con los datos codificados en OneHot (de las especies) porque se estaría dividiendo entre 0, se trata simplemente de confrimar la hipótesis.

from statsmodels.stats.outliers_influence import variance_inflation_factor

# Importante eliminar variable objetivo y categórica
df_vif = df.drop(['Weight', 'Species'], axis=1)
vif = pd.DataFrame()
vif['Features'] = df_vif.columns
vif['VIF'] = [variance_inflation_factor(df_vif.values, i) for i in range(df_vif.shape[1])]
vif

	Features	VIF
0	Length1	12749.616323
1	Length2	16580.478064
2	Length3	3382.291295
3	Height	76.053658
4	Width	92.653797

Demostrando una fuerte multicolinealidad entre variables explicativas, puesto que VIF apenas debería superar 10. Según el libro de Freddy Hernández, Olga Usuga y Mauricio Mazo, Modelos de Regresión con R, capítulo 18 (https://fhernanb.github.io/libro_regresion/multicoli.html), esto puede implicar varios aspectos negativos:

• Grandes varianzas y covarianzas de los estimadores.
• Estimaciones para los coeficientes demasiado grandes.
• Pequeños cambios en los datos provocan grandes cambios en las estimaciones de los coeficientes.
• Las estimaciones pueden darse en magnitudes poco razonables.

Además, por motivos obvios, se dificulta la interpretación de los coeficientes.

Entre las soluciones a este problema se plantean 2 opciones que más tarde se evaluarán en la práctica.

1. Creación de una sola variable para Lenght.
2. Eliminación de estas variables.

# Lenght fusionadas
df_vif_merged = df.copy()
df_vif_merged['Length'] = df_vif_merged[['Length1', 'Length2', 'Length3']].mean(axis=1)
df_vif_merged = df_vif_merged.drop(['Length1', 'Length2', 'Length3', 'Weight', 'Species'], axis=1)
vif = pd.DataFrame()
vif['Features'] = df_vif_merged.columns
vif['VIF'] = [variance_inflation_factor(df_vif_merged.values, i) for i in range(df_vif_merged.shape[1])]
vif

	Features	VIF
0	Height	14.700943
1	Width	50.745948
2	Length	32.533848

# Sin Lenght
df_vif_no_length = df.drop(['Weight', 'Species', 'Length1', 'Length2', 'Length3'], axis=1)
vif = pd.DataFrame()
vif['Features'] = df_vif_no_length.columns
vif['VIF'] = [variance_inflation_factor(df_vif_no_length.values, i) for i in range(df_vif_no_length.shape[1])]
vif

	Features	VIF
0	Height	14.525751
1	Width	14.525751

Claramente, los valores VIF han mejorado notoriamente, con todo, no dejan de ser relativamente altos dada la naturaleza de los datos. Pues sería extraño no mantener relaciones proporcionales en altura, anchura, etc.

Para terminar con la exploración, es necesario comprobar los valores atípicos.

# Escalado de las numéricas para facilitar la visualización
from sklearn.preprocessing import MinMaxScaler
scaler = MinMaxScaler()
df_scaled = df.copy()
df_scaled[numeric_cols.columns] = scaler.fit_transform(numeric_cols)

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
sns.histplot(df["Species"], ax=axes[0])
axes[0].set_title('Species Distribution')
sns.boxplot(data=df_scaled, orient="v", ax=axes[1])
axes[1].set_title('Outliers')
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()

En referencia a las especies, es claro que "Whitefish" es un outlier, sin embargo, para acercar el ejercicio a un uso real, se mantienen estos datos puesto que idealmente, también se querrá clasificar esta especie. Para el resto de valores numéricos, se eliminarán datos atípicos mediante el rango intercuantílico.

d) Preprocesamiento

Gracias a la exploración se han identificado estas tareas:

1. Eliminación de datos inconsitentes: hecha ✔️.
2. Creación de 3 conjuntos según la variable Lenght_n: todas incluidas, fusionadas, todas excluidas.
3. Eliminación de valores outliers.

Además, para desarrollar modelos predictivos, se proponen:

4. Codificación One-Hot para Species, puesto que es nominal.
5. Normalización, especialmente importante para algunos modelos y preferible por encima de Min-Max dadas las distribuciones asimétricas.

Puesto que las tareas relacionadas con outliers y One-Hot son comunes, se realizan como primer paso.

# Outliers
Q1 = df[numeric_cols.columns].quantile(0.25)
Q3 = df[numeric_cols.columns].quantile(0.75)
IQR = Q3 - Q1
outliers = (df[numeric_cols.columns] < (Q1 - 1.5 * IQR)) | (df[numeric_cols.columns] > (Q3 + 1.5 * IQR))
print("Size with outliers:", df.shape)
df_no_outliers = df[~outliers.any(axis=1)]
print("Size without outliers:", df_no_outliers.shape)

# One-hot
from sklearn.preprocessing import OneHotEncoder
encoder = OneHotEncoder(sparse_output=False)
species_encoded = encoder.fit_transform(df_no_outliers[['Species']])
species_df = pd.DataFrame(species_encoded, columns=encoder.get_feature_names_out(['Species']), index=df_no_outliers.index)
df_no_outliers = df_no_outliers.drop('Species', axis=1)
df_no_outliers = pd.concat([df_no_outliers, species_df], axis=1)
df_no_outliers.head()

Size with outliers: (158, 7)
Size without outliers: (155, 7)

	Weight	Length1	Length2	Length3	Height	Width	Species_Bream	Species_Parkki	Species_Perch	Species_Pike	Species_Roach	Species_Smelt	Species_Whitefish
0	242.0	23.2	25.4	30.0	11.5200	4.0200	1.0	0.0	0.0	0.0	0.0	0.0	0.0
1	290.0	24.0	26.3	31.2	12.4800	4.3056	1.0	0.0	0.0	0.0	0.0	0.0	0.0
2	340.0	23.9	26.5	31.1	12.3778	4.6961	1.0	0.0	0.0	0.0	0.0	0.0	0.0
3	363.0	26.3	29.0	33.5	12.7300	4.4555	1.0	0.0	0.0	0.0	0.0	0.0	0.0
4	430.0	26.5	29.0	34.0	12.4440	5.1340	1.0	0.0	0.0	0.0	0.0	0.0	0.0

Tal y como se aprecia, los 8 valores extremos eran de solo 3 instancias y las especies han quedado correctamente codificadas.

Los últimos 2 pasos son la creación de conjuntos y normalización.

# Fusión
df_merged = df_no_outliers.copy()
df_merged['Length'] = df_merged[['Length1', 'Length2', 'Length3']].mean(axis=1)
df_merged = df_merged.drop(['Length1', 'Length2', 'Length3'], axis=1)
# Sin Length
df_reduced = df_merged.drop(['Length'], axis=1)

# Normalización
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
def split_scale(df, target_column='Weight', test_size=0.2, random_state=1):
    X = df.drop(target_column, axis=1)
    y = df[target_column]
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=test_size, random_state=random_state
    )

    scaler = StandardScaler()
    scaler.fit(X_train)
    X_train_scaled = pd.DataFrame(scaler.transform(X_train), columns=X_train.columns, index=X_train.index)
    X_test_scaled = pd.DataFrame(scaler.transform(X_test), columns=X_test.columns, index=X_test.index)

    return X_train_scaled, X_test_scaled, y_train, y_test


X_train_all, X_test_all, y_train_all, y_test_all = split_scale(df_no_outliers)
X_train_merged, X_test_merged, y_train_merged, y_test_merged = split_scale(df_merged)
X_train_reduced, X_test_reduced, y_train_reduced, y_test_reduced = split_scale(df_reduced)

print("nº Rows - train:", X_train_all.shape[0])
print("nº Rows - test:", X_test_all.shape[0])
X_train_all.head()

nº Rows - train: 124
nº Rows - test: 31

	Length1	Length2	Length3	Height	Width	Species_Bream	Species_Parkki	Species_Perch	Species_Pike	Species_Roach	Species_Smelt	Species_Whitefish
155	-1.545952	-1.586587	-1.613530	-1.541899	-1.858907	-0.552506	-0.296174	-0.728604	-0.342224	-0.342224	3.205110	-0.204980
67	-0.751337	-0.744400	-0.711776	0.087342	-0.580481	-0.552506	3.376389	-0.728604	-0.342224	-0.342224	-0.312002	-0.204980
33	1.251526	1.315405	1.398512	2.248297	1.411664	1.809934	-0.296174	-0.728604	-0.342224	-0.342224	-0.312002	-0.204980
157	-1.382675	-1.393797	-1.455490	-1.438335	-1.382378	-0.552506	-0.296174	-0.728604	-0.342224	-0.342224	3.205110	-0.204980
60	1.240641	1.213936	1.175398	0.779173	1.278949	-0.552506	-0.296174	-0.728604	-0.342224	-0.342224	-0.312002	4.878524

e) Resumen final datos limpios

A partir de los datos iniciales, se han obtenido 155 registros sin valores atípicos, organizados en 3 conjuntos. En todos ellos aparecen al menos 2 variables explicativas relacionadas con la medida (alto y ancho), junto con la especie (codificada en 7 variables mediante One-Hot Encoding).

La creación de estos 3 conjuntos se debe a la clara multicolinealidad entre las variables "length":

• Conjunto completo: 7 (especies) + 2 (alto y ancho) + 3 (length).
• Conjunto fusionado: 7 (especies) + 2 (alto y ancho) + 1 (length promedio).
• Conjunto sin length.

Además, la exploración muestra un notable desbalance en la frecuencia de especies. En cuanto a las variables numéricas, todas en cm salvo el peso, presentan rangos muy diferentes (p. ej. 70 vs 10). Este hecho, unido a distribuciones asimétricas y a la sensibilidad de algunos modelos a escalas y distancias, ha motivado la decisión de normalizar todos los valores.

3. Modelado, Predicción y Evaluación "simples"

Es fundamental reconocer las limitaciones del conjunto. Con tan pocos datos, una simple división en entrenamiento y test puede generar problemas:

• Con 80% (124 datos) para entrenar y 20% (31) para test, la evaluación depende de un conjunto muy pequeño.
• Con 70% (108) para entrenar y 30% (47) para test, el modelo pierde demasiada información o aprende poco.

Para mitigar esta situación, se propone:

1. Test reducido (20%) como evaluación final.
2. k-fold cross-validation sobre el 80% : los datos de entreno se dividen en k folds. En cada iteración, se entrena con k-1 folds y se valida con el restante. Así, todas las observaciones participan en entrenamiento y validación (pero nunca en la misma iteración), promediando las métricas obtenidas.

Esta estrategia reduce la dependencia de una única partición. Asimismo, ofrece métricas más veridicas al reducir la dependencia de un test tan pequeño mientras se incluye la validación con datos variados.

Se emplean los siguientes modelos: Linear Regression, Tree Regresion y kNN Regression.

from sklearn.model_selection import RepeatedKFold
from sklearn.linear_model import LinearRegression
from sklearn.neighbors import KNeighborsRegressor
from sklearn.tree import DecisionTreeRegressor

# Modelos
lr = LinearRegression()
tree = DecisionTreeRegressor(max_depth=5, random_state=1) # Para limitar el overfitting
knn = KNeighborsRegressor(n_neighbors=3) # Impar para evitar empates

# cross validation
cv = RepeatedKFold(n_splits=10, n_repeats=3, random_state=1)

Puesto que el objetivo de los modelos no es el de clasificación sino el de predicción de un valor continuo, las métricas de evauluación están principalmente basadas en el error.

Métrica	Unidades
Mean Absolute Error (MAE): media del valor absoluto de los errores.	Weight (gramos)
Mean Squared Error (MSE): media de los errores cuadrados.	(gramos)²
Root Mean Squared Error (RMSE): raíz cuadrada del MSE, refleja mejor los errores extremos.	Weight (gramos)

Por otro lado, se ha usado el coeficiente de determinación R^2, cuyo valor ideal es 1 e indica la proporción de varianza explicada por el modelo respecto a la variable dependiente. Simplemente una métrica diferente y útil.

from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
from sklearn.model_selection import cross_val_score
import numpy as np

# Puesto que hay 3 conjuntos y 4 métricas, se crea una función para sintetizar código
def evaluate_model(title, model, X_train, X_test, Y_train, Y_test):
    # Cross-validation scores para cada métrica
    scores_MAE = cross_val_score(model, X_train, Y_train, scoring='neg_mean_absolute_error', cv=cv, n_jobs=-1, error_score='raise')
    scores_MSE = cross_val_score(model, X_train, Y_train, scoring='neg_mean_squared_error', cv=cv, n_jobs=-1, error_score='raise')
    scores_RMSE = cross_val_score(model, X_train, Y_train, scoring='neg_root_mean_squared_error', cv=cv, n_jobs=-1, error_score='raise')
    scores_r2 = cross_val_score(model, X_train, Y_train, scoring='r2', cv=cv, n_jobs=-1, error_score='raise')

    # Se convierten los errores en absolutos excepto r2
    scores_MAE = np.abs(scores_MAE)
    scores_MSE = np.abs(scores_MSE)
    scores_RMSE = np.abs(scores_RMSE)

    # Entrenamiento y test
    model.fit(X_train, Y_train)
    Y_pred = model.predict(X_test)
    mae_test = mean_absolute_error(Y_test, Y_pred)
    mse_test = mean_squared_error(Y_test, Y_pred)
    rmse_test = np.sqrt(mse_test)
    r2_test = r2_score(Y_test, Y_pred)

    # Devolver los resultados para más tarde concatenarlos
    r = lambda x: np.round(x, 4)
    results = pd.DataFrame({
        'DATASET': [title],
        'MODEL': [str(model)],
        'MAE': [r(np.mean(scores_MAE))],
        'MAE_TEST': [r(mae_test)],
        'MSE': [r(np.mean(scores_MSE))],
        'MSE_TEST': [r(mse_test)],
        'RMSE': [r(np.mean(scores_RMSE))],
        'RMSE_TEST': [r(rmse_test)],
        'R2': [r(np.mean(scores_r2))],
        'R2_TEST': [r(r2_test)]
    })

    return results

results_full_lr = evaluate_model("FULL DATASET", lr, X_train_all, X_test_all, y_train_all, y_test_all)
results_merged_lr = evaluate_model("MERGED DATASET", lr, X_train_merged, X_test_merged, y_train_merged, y_test_merged)
results_reduced_lr = evaluate_model("REDUCED DATASET", lr, X_train_reduced, X_test_reduced, y_train_reduced, y_test_reduced)

results_full_tree = evaluate_model("FULL DATASET", tree, X_train_all, X_test_all, y_train_all, y_test_all)
results_merged_tree = evaluate_model("MERGED DATASET", tree, X_train_merged, X_test_merged, y_train_merged, y_test_merged)
results_reduced_tree = evaluate_model("REDUCED DATASET", tree, X_train_reduced, X_test_reduced, y_train_reduced, y_test_reduced)

results_full_knn = evaluate_model("FULL DATASET", knn, X_train_all, X_test_all, y_train_all, y_test_all)
results_merged_knn = evaluate_model("MERGED DATASET", knn, X_train_merged, X_test_merged, y_train_merged, y_test_merged)
results_reduced_knn = evaluate_model("REDUCED DATASET", knn, X_train_reduced, X_test_reduced, y_train_reduced, y_test_reduced)

results = pd.concat([
    results_full_lr, results_merged_lr, results_reduced_lr,
    results_full_tree, results_merged_tree, results_reduced_tree,
    results_full_knn, results_merged_knn, results_reduced_knn
], ignore_index=True)

results

	DATASET	MODEL	MAE	MAE_TEST	MSE	MSE_TEST	RMSE	RMSE_TEST	R2	R2_TEST
0	FULL DATASET	LinearRegression()	58.6017	67.0154	6943.2253	6880.1861	79.5771	82.9469	0.9039	0.9386
1	MERGED DATASET	LinearRegression()	58.0673	67.4560	6693.0946	6993.9075	78.1998	83.6296	0.9071	0.9376
2	REDUCED DATASET	LinearRegression()	59.8655	69.2021	7271.0819	7495.1281	81.4250	86.5744	0.9052	0.9331
3	FULL DATASET	DecisionTreeRegressor(max_depth=5, random_stat...	49.9134	40.0244	6188.6824	4314.9313	74.1643	65.6881	0.9289	0.9615
4	MERGED DATASET	DecisionTreeRegressor(max_depth=5, random_stat...	49.6292	56.1030	5907.6136	7121.2920	72.3075	84.3877	0.9337	0.9365
5	REDUCED DATASET	DecisionTreeRegressor(max_depth=5, random_stat...	60.8270	64.7768	11339.6729	8778.8098	96.3056	93.6953	0.8720	0.9217
6	FULL DATASET	KNeighborsRegressor(n_neighbors=3)	48.1915	39.8957	8703.7614	3816.0942	79.7525	61.7745	0.9100	0.9660
7	MERGED DATASET	KNeighborsRegressor(n_neighbors=3)	48.1594	41.3925	8591.5752	4750.7042	79.9808	68.9254	0.9113	0.9576
8	REDUCED DATASET	KNeighborsRegressor(n_neighbors=3)	52.9596	56.5409	9164.1179	7208.1439	83.7584	84.9008	0.9021	0.9357

Afortunadamente, se puede afirmar que no hay un overfitting, de hecho en muchos casos, los resultados en el test son incluso ligeramente mejores que en la validación cruzada.

Antes de abordar al análisis, véase la siguiente celda para comprender la magnitud de los errores.

print("Weight range:", y_test_all.min(), "to", y_test_all.max())
plt.figure(figsize=(7, 4))
sns.boxplot(x='Species', y='Weight', data=df, hue='Species', palette='Set3')
plt.title('Boxplot: Weight by Species')
plt.xticks(rotation=45)
plt.show()

Weight range: 7.5 to 1100.0

Ahora, es claro que todos los modelos tienen un buen desempeño. El mejor —que es el de árbol de decisión con el conjunto completo sin reducir—, tiene en el test un RMSE de 74g y un R^2 de 0.96. Por otro lado, el MAE que se registra en la unidad de Weight obtiene 40g en un rango absoluto de 1092.5g (1100 - 7.5).

Enfocándose en los modelos KNeighbors, el conjunto sin variables reducidas o completo es el segundo mejor y prácticamente el primero. Por ejemplo, tiene un MAE media décima más baja (mejor) que el primero en la validación cruzada, y es por otras pequeñísimas diferencias, ligeramente peor que el primero en otras métricas. Mientras, se puede considerar la peor opción la regresión lineal. Puesto que su mejor modelo, es por algunas décimas (o ~5-20g) peor en todas las métricas que los otros 2 algorimtos.

Respecto a los conjuntos de datos, aquel sin variables Lenght tiene un desempeño muy similar a los otros 2 en la regresión lineal, sin embargo, empeora entre un poco y bastante el resto de métricas en los otros modelos. Mientras, el conjunto con el promedio de las variables Lenght, le sucede lo mismo pero en menor medida. De forma más intuitiva se pueden ordenar por desempeño descendente:

1. Conjunto con todos los datos (mejor).
2. Conjunto Lenght promedidado.
3. Conjunto sin Lenght.

4. Alternativas de ensamblado

Con el propósito de mejorar la predicción, se investigan 2 estrategias de ensamblado. En primer lugar, un GradientBoostingRegressor (ensamblado secuencial de árboles débiles) y, en segundo lugar, un ensamblado personalizado.

from sklearn.ensemble import GradientBoostingRegressor

# Documentación: https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingRegressor.html
# De acuerdo con la fuente y después de algunas pruebas, se ha decidido mantener los valores por defecto de los parámetros
# excepto para el criterio en la calidad del split y la profunidad máxima

gbr = GradientBoostingRegressor(criterion='squared_error', max_depth=5, random_state=1)
cv = RepeatedKFold(n_splits=10, n_repeats=3, random_state=1)

gbr_full_dataset = evaluate_model("FULL DATASET", gbr, X_train_all, X_test_all, y_train_all, y_test_all)
gbr_merged_dataset = evaluate_model("MERGED DATASET", gbr, X_train_merged, X_test_merged, y_train_merged, y_test_merged)
gbr_reduced_dataset = evaluate_model("REDUCED DATASET", gbr, X_train_reduced, X_test_reduced, y_train_reduced, y_test_reduced)

gbr_results = pd.concat([
    gbr_full_dataset, gbr_merged_dataset, gbr_reduced_dataset
], ignore_index=True)

gbr_results

	DATASET	MODEL	MAE	MAE_TEST	MSE	MSE_TEST	RMSE	RMSE_TEST	R2	R2_TEST
0	FULL DATASET	GradientBoostingRegressor(criterion='squared_e...	39.2101	36.7589	4270.2565	2765.9324	60.0798	52.5921	0.9525	0.9753
1	MERGED DATASET	GradientBoostingRegressor(criterion='squared_e...	40.2388	34.6544	4524.5443	2990.2331	61.6007	54.6830	0.9496	0.9733
2	REDUCED DATASET	GradientBoostingRegressor(criterion='squared_e...	51.4215	53.6461	8600.4321	7080.1950	83.8593	84.1439	0.9069	0.9368

En estos modelos, las métricas entre el conjunto de prueba y la validación cruzada no son solo más parecidos, sino que en muchas ocasiones, ¡mejoran notablemente con datos nuevos! Esto indica con bastante contundencia que no se está sobreajustando y generaliza bien, por lo que los resultados son ahora mucho más fiables.

Los datasets que contienen las variables Lenght (las 3 sin modificar y también el promedio) tienen métricas prácticamente iguales. Siendo sorprendente la capacidad de evitar el daño de la multicolinealidad que se planteaba al principio. Por su parte, el peor modelo vuelve a ser aquel que no registra Lenght.

Comparativamente, el GradientBoostingRegressor supone una mejora considerable a los vistos anteriormente.

Métrica	Mejor Anterior	GradientBoostingRegressor
MAE_TEST	39.89	36.75
RMSE_TEST	72.30	52.59
R2_TEST	0.9660	0.9753

Véase ahora el ensamblaje personalizado con los 3 mejores modelos vistos hasta el momento.

import numpy as np
from sklearn.base import BaseEstimator, RegressorMixin

class CustomEnsembleRegressor(BaseEstimator, RegressorMixin):
    def __init__(self, models):
        self.models = models

    def fit(self, X, y):
        for model in self.models:
            model.fit(X, y)
        return self

    def predict(self, X):
        # Cada modelo aporta a la media de las predicciones
        predictions = []
        for model in self.models:
            predictions.append(model.predict(X))
        return np.mean(predictions, axis=0)

custom_ensemble = CustomEnsembleRegressor([
    GradientBoostingRegressor(criterion='squared_error', max_depth=5, random_state=1),
    DecisionTreeRegressor(max_depth=5, random_state=1),
    KNeighborsRegressor(n_neighbors=3),
])

ensemble_full_results = evaluate_model("Full Dataset", custom_ensemble, X_train_all, X_test_all, y_train_all, y_test_all)
ensemble_merged_results = evaluate_model("Merged Lenghts", custom_ensemble, X_train_merged, X_test_merged, y_train_merged, y_test_merged)
ensemble_reduced_results = evaluate_model("No Lenghts", custom_ensemble, X_train_reduced, X_test_reduced, y_train_reduced, y_test_reduced)

ensemble_results = pd.concat([
    ensemble_full_results, ensemble_merged_results, ensemble_reduced_results],
    ignore_index=True)

ensemble_results

	DATASET	MODEL	MAE	MAE_TEST	MSE	MSE_TEST	RMSE	RMSE_TEST	R2	R2_TEST
0	Full Dataset	CustomEnsembleRegressor(models=[GradientBoosti...	38.8697	33.6842	4329.3891	2599.8246	58.4152	50.9885	0.9536	0.9768
1	Merged Lenghts	CustomEnsembleRegressor(models=[GradientBoosti...	38.4730	36.9433	4303.5802	3462.4614	58.4038	58.8427	0.9542	0.9691
2	No Lenghts	CustomEnsembleRegressor(models=[GradientBoosti...	47.8149	52.1579	7860.4286	6605.2896	77.6991	81.2729	0.9168	0.9411

Complementariamente, se comprueba la evaluación de los 3 primeros modelos.

custom_ensemble2 = CustomEnsembleRegressor([LinearRegression(),DecisionTreeRegressor(max_depth=5, random_state=1),KNeighborsRegressor(n_neighbors=3)])
ensemble2_full_results = evaluate_model("Full Dataset", custom_ensemble2, X_train_all, X_test_all, y_train_all, y_test_all)
ensemble2_merged_results = evaluate_model("Merged Lenghts", custom_ensemble2, X_train_merged, X_test_merged, y_train_merged, y_test_merged)
ensemble2_reduced_results = evaluate_model("No Lenghts", custom_ensemble2, X_train_reduced, X_test_reduced, y_train_reduced, y_test_reduced)
ensemble2_results = pd.concat([ensemble2_full_results, ensemble2_merged_results, ensemble2_reduced_results])
ensemble2_results

	DATASET	MODEL	MAE	MAE_TEST	MSE	MSE_TEST	RMSE	RMSE_TEST	R2	R2_TEST
0	Full Dataset	CustomEnsembleRegressor(models=[LinearRegressi...	38.5943	36.8134	4332.8665	3045.2710	58.2879	55.1840	0.9533	0.9728
0	Merged Lenghts	CustomEnsembleRegressor(models=[LinearRegressi...	37.3253	41.8669	4174.3025	4008.2884	57.2897	63.3110	0.9554	0.9642
0	No Lenghts	CustomEnsembleRegressor(models=[LinearRegressi...	45.1098	50.8831	6433.9348	5810.3574	70.1776	76.2257	0.9323	0.9482

En referencia a estos métodos personalizados, se observa:

1. Modelo 1 (GB+Tree+KNN):

• No mejoran las métricas del mejor modelo (GB).
• Métricas de test notablemente mejores que la valdiación cruzada para el mejor modelo (buena generalización).
• Los mejores resultados siguen siendo del conjunto donde Lenght no se modifica.

2. Modelo 2 (LR+Tree+KNN):

• Tampoco mejoran el modelo GB
• Sí mejora —y notablemente— respecto al mejor modelo de esta combinación (Tree). Ver tabla.
• El mejor conjunto sigue siendo el mismo y sigue mejorando en el test.

Métrica	Modelo 2	Tree
MAE_TEST	36.8	40.0
RMSE_TEST	55.2	65.7
R2	0.95	0.93

En definitiva, promediar la predicción mediante 3 modelos ha mejorado los resultados individuales iniciales y los ha hecho más fiables aunque el rendimiento sea menor que el boosting en todos los casos. Con el fin de explotar esta técnica de modelaje, queda pendiente la ampliación a más modelos —como SVR o RF— con una búsqueda exhaustiva de hiperparámetros (grid search). De este modo, sería posible buscar combinaciones de n modelos iterativamente que mejorasen los resultados actuales.

5. Análisis de errores

El mejor modelo no es siempre aquel que maximiza la exactitud o minimiza el error, sino el que más se ajusta al propósito establecido. En este caso, se desconoce el precio de venta de cada especie de pescado, por lo que no es posible minimizar pérdidas por esta vía. En su defecto, se propone el siguiente gráfico que muestra el comportamiento del modelo según el peso.

Notar, que tan solo se valora el conjunto con todas las variables, claramente, el mejor en todos los casos.

# Puesto que se han usado diferentes características en el entreno (conjuntos),
# se vuelven a llamar los modelos con el dataset completo para poder predecir
# sino da error
evaluate_model("FULL DATASET", lr, X_train_all, X_test_all, y_train_all, y_test_all)
evaluate_model("FULL DATASET", tree, X_train_all, X_test_all, y_train_all, y_test_all)
evaluate_model("FULL DATASET", knn, X_train_all, X_test_all, y_train_all, y_test_all)
evaluate_model("FULL DATASET", gbr, X_train_all, X_test_all, y_train_all, y_test_all)
evaluate_model("FULL DATASET", custom_ensemble, X_train_all, X_test_all, y_train_all, y_test_all)
evaluate_model("FULL DATASET", custom_ensemble2, X_train_all, X_test_all, y_train_all, y_test_all)

pred_lr = lr.predict(X_test_all)
pred_tree = tree.predict(X_test_all)
pred_knn = knn.predict(X_test_all)
pred_gbr = gbr.predict(X_test_all)
pred_custom_ensemble = custom_ensemble.predict(X_test_all)
pred_custom_ensemble2 = custom_ensemble2.predict(X_test_all)

results = pd.DataFrame({
    'Real value': y_test_all,
    'LR': pred_lr,
    'Tree': pred_tree,
    'KNN': pred_knn,
    'GBR': pred_gbr,
    'Custom_Ensemble': pred_custom_ensemble,
    'Custom_Ensemble2': pred_custom_ensemble2
})

# Errores absolutos real - predicción
results['Abs Error LR'] = np.abs(results['Real value'] - results['LR'])
results['Abs Error Tree'] = np.abs(results['Real value'] - results['Tree'])
results['Abs Error KNN'] = np.abs(results['Real value'] - results['KNN'])
results['Abs Error GBR'] = np.abs(results['Real value'] - results['GBR'])
results['Abs Error Custom Ensemble'] = np.abs(results['Real value'] - results['Custom_Ensemble'])
results['Abs Error Custom Ensemble 2'] = np.abs(results['Real value'] - results['Custom_Ensemble2'])
results = results.round(2)
results.head(3)

	Real value	LR	Tree	KNN	GBR	Custom_Ensemble	Custom_Ensemble2	Abs Error LR	Abs Error Tree	Abs Error KNN	Abs Error GBR	Abs Error Custom Ensemble	Abs Error Custom Ensemble 2
118	820.0	873.89	897.50	921.67	918.30	912.49	897.69	53.89	77.50	101.67	98.30	92.49	77.69
76	70.0	0.67	83.57	76.50	75.92	78.66	53.58	69.33	13.57	6.50	5.92	8.66	16.42
52	290.0	339.65	273.38	217.33	276.70	255.80	276.79	49.65	16.62	72.67	13.30	34.20	13.21

plt.style.use("seaborn-v0_8-whitegrid")

models = ['LR', 'Tree', 'KNN', 'GBR', 'Custom_Ensemble', 'Custom_Ensemble2']
colors = sns.color_palette("Set2", len(models))

# Plot 3x2
fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(14, 12))
fig.suptitle("Models Error Comparison", fontsize=18, fontweight="bold")

for i, model in enumerate(models):
    row, col = i // 2, i % 2
    ax = axes[row, col]
    ax.set_title(f"{model}", fontsize=14, fontweight="semibold")
    # Puntos
    ax.scatter(results.index, results['Real value'],
               color='black', s=50, label='Real', alpha=0.8, marker="o")
    ax.scatter(results.index, results[model],
               color=colors[i], s=50, label='Prediction', alpha=0.7, marker="^")

    # Líneas de error
    for idx in results.index:
        ax.plot([idx, idx], [results['Real value'][idx], results[model][idx]],
                linestyle='dashed', color='red', alpha=0.4, linewidth=1)

    ax.legend(loc='upper center', frameon=True, framealpha=0.8)
    ax.set_xlabel("Index", fontsize=11)
    ax.set_ylabel("Weight", fontsize=11)

# Espacio entre plots
fig.subplots_adjust(hspace=0.4, wspace=0.3)
plt.show()

Lo más llamativo de estos gráficos son los errores negativos de los modelos que emplean regresión lineal (LR y Custom_Ensemble2). Si bien este problema parece ser puntual, el error escalado puede provocar pérdidas monetarias injustificadas. De igual manera, algunos factores predictivos tan grandes (ver error ~0 vs -200) suponen un riesgo inaceptable.

Por otro lado, todos los modelos pierden exactitud absoluta (no relativa) cuanto más aumenta el peso. Para comprender mejor su comportamiento, véase a escala relativa.

# Error relativo y categorización por peso
models = [model for model in models if model not in ['Custom_Ensemble2', 'LR']]
for model in models:
    results[f'Rel Error {model}'] = (
        abs(results['Real value'] - results[model]) / results['Real value']
) * 100
bins = [0, 200, 600, 1100]
labels = ['0-200g', '201-600g', '601-1100g']
results['Fish_Weight'] = pd.cut(results['Real value'], bins=bins, labels=labels)
size_map = {'0-200g': 50, '201-600g': 150, '601-1100g': 400}
line_colors = sns.color_palette("Set2", len(size_map))
cat_colors = dict(zip(size_map.keys(), sns.color_palette("Set2", len(size_map))))

# Plot
fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(14, 8))
fig.suptitle("Relative Error by Model", fontsize=18, fontweight="bold")
for i, model in enumerate(models):
    row, col = i // 2, i % 2
    ax = axes[row, col]
    ax.set_title(f"{model}", fontsize=14, fontweight="semibold")

    # Puntos
    ax.scatter(
        results.index,
        results[f'Rel Error {model}'],
        s=results['Fish_Weight'].map(size_map),
        c=results['Fish_Weight'].astype(str).map(cat_colors),
        alpha=0.7,
        edgecolor='black'
    )

    # Cálculo media por categoría
    for j, cat in enumerate(size_map.keys()):
        mean_val = results.loc[results['Fish_Weight'] == cat, f'Rel Error {model}'].mean()
        ax.axhline(mean_val, linestyle='dashed', color=line_colors[j],
                    linewidth=1.5, alpha=0.9,
                    label=f"Mean {cat}: {mean_val:.1f}%")
        ax.text(results.index[-1], mean_val, f"{mean_val:.1f}%", color=line_colors[j],
                ha='left', va='baseline', fontsize=12, fontweight='semibold')

    ax.set_xlabel("Index", fontsize=11)
    ax.set_ylabel("Relative Error (%)", fontsize=11)

    # Leyenda: tamaños y colores
    handles, labels_plot = [], []
    for cat in size_map.keys():
        handles.append(plt.scatter([], [], s=size_map[cat],
                                   color=cat_colors[cat], alpha=0.7, edgecolor='black'))
        labels_plot.append(cat)
        handles.append(plt.Line2D([0], [0], color=line_colors[list(size_map.keys()).index(cat)],
                                  linestyle='--', lw=1.5))
        labels_plot.append(f"Avg Error (%)")
    ax.legend(handles, labels_plot, loc='best', frameon=True, framealpha=0.8)

fig.subplots_adjust(hspace=0.4, wspace=0.3)
plt.show()

print("Avg Error (%) GBR:", results['Rel Error GBR'].mean().round(2))
print("Avg Error (%) Custom_Ensemble:", results['Rel Error Custom_Ensemble'].mean().round(2))

Avg Error (%) GBR: 9.89
Avg Error (%) Custom_Ensemble: 9.72

Los modelos kNN y Tree quedan descartados. Su desempeño en la predicción de pescado pesado (600g o más) es muy similar al de los modelos ensamblados, pero sus predicciones para peces menores es notablemente peor.

Tanto GBR como Custom_Ensemble cometen errores de forma muy similar, diferenciándose por algo más de 1 punto porcentual entre los peces medios y pequeños. Siendo ***GBR* el mejor predictor para los pequeños y Custom_Ensemble el mejor para los medianos**. Gracias a este gráfico es claro que el GBR ofrece mejores métricas en el test debido al buen desempeño en la estimación de los peces pequeños, que son los más numerosos en este conjunto tan reducido.

Para estas conclusiones, se ha empleado a consciencia la media, métrica en la que influyen fuertemente los valores extremos. De esta forma, se proporcionar una visión relalista de las ganancias y pérdidas producidas por modelo. Con todo, también se ha empleado un error relativo absoluto, por lo que del gráfico anterior no se concluye una equivocación acumulativa del ~10%.

Ver estimación:

for model in models:
    results[f'Signed Rel Error {model}'] = (
        (results[model] - results['Real value']) / results['Real value']
    ) * 100

for model in ['GBR', 'Custom_Ensemble']:
    col_name = f'Signed Rel Error {model}'
    # Media global
    mean_global = results[col_name].mean()
    print(f"---Avg Error Signed (%) {model}: {mean_global:.2f}---")
    # Media por Fish_Weight
    mean_by_weight = results.groupby('Fish_Weight', observed=False)[col_name].mean()
    for cat, val in mean_by_weight.items():
        print(f"Avg Error Signed (%) {model}, {cat}: {val:.2f}")
    print()

---Avg Error Signed (%) GBR: -0.70---
Avg Error Signed (%) GBR, 0-200g: -0.73
Avg Error Signed (%) GBR, 201-600g: 3.82
Avg Error Signed (%) GBR, 601-1100g: -5.18

---Avg Error Signed (%) Custom_Ensemble: -0.61---
Avg Error Signed (%) Custom_Ensemble, 0-200g: 2.68
Avg Error Signed (%) Custom_Ensemble, 201-600g: -1.99
Avg Error Signed (%) Custom_Ensemble, 601-1100g: -5.41

Por lo tanto se concluye:

1. Ambos modelos se equivocan en un rango [0%, 30%].
2. El error suele ser del 10%, ver distribución figura anterior.
3. Los errores se compensan obteniendo una tasa aproximada de -0.65%.
4. GBR es excelente para peces de peso inferior a 200g.
5. Custom_Ensemble es muy bueno (error compensado del 2%) en peces medianos.
6. Ambos son solamente buenos para peces más pesados. Probablemente causado por la falta de estos registros. Todos los modelos se comportan peor en estos casos.

6. Conclusiones y Trabajo Futuro

El objetivo fundamental de este proyecto ha sido desarrollar un modelo predictivo para estimar el peso del pescado con la máxima precisión posible, basándose en sus dimensiones y especie. A pesar de las limitaciones inherentes al conjunto de datos, los resultados obtenidos son notables.

Mediante la implementación de modelos de Gradient Boosting Regressor (GBR) y ensamblado personalizados, se han alcanzado tasas de errores relativos acumulados inferiores al 1% (e.g. unas veces predice +20g, otras -20g y se compensa en ~1%). El modelo GBR tiene un error absoulto promedio o MAE de tan solo 37g, mientras alcanza una bondad de ajuste R^2 de 0.975. La validación cruzada respalda la robustez de estos resultados, que demuestran una alta capacidad predictiva, especialmente en los rangos de peso más comunes.

El rendimiento de los modelos varía según el peso del pez, observándose una clara tendencia:

• Excelente rendimiento: En peces de peso inferior a 200g.
• Buen rendimiento: En peces con pesos entre 200g y 600g.
• Rendimiento moderado: En peces de más de 600g, donde la escasez de datos reduce la precisión.

A modo de ejemplo práctico, para una muestra agregada de 10 kg de pescado, se obtendría el siguiente error:

Tamaño	Error aproximado en 10kg
Pequeños (<200 g)	70 g
Medianos (200–600 g)	200 g
Grandes (>600 g)	520 g

Dada la variabilidad en el error de predicción individual, que puede oscilar entre un 0% y un 30% y que suele ser del 10%, se desaconseja firmemente el uso de este modelo para estimaciones con un número reducido de ejemplares (jamás en inferiores a 100).

Limitaciones del Estudio

El análisis y los modelos resultantes están condicionados por cuatro factores principales derivados del conjunto de datos:

1. Conjunto de datos reducido: el número limitado de registros reduce la capacidad de generalización de los modelos y su consistencia en datos no observados.
2. Distribuciones asimétricas: la mayor concentración de datos en peces de menor tamaño limita la capacidad de aprendizaje y evaluación para los ejemplares más grandes, que son menos frecuentes.
3. Especies infrarrepresentadas: la baja frecuencia de algunas especies (con tan solo 6 u 11 registros) compromete significativamente la utilidad y exactitud de las predicciones para dichos grupos.
4. Multicolinealidad: aunque es un factor de menor impacto para la predicción, la correlación natural entre las dimensiones de los peces es inevitable y limita, entre otros, la interpretabilidad de modelos más simples.

Líneas de Trabajo Futuro y Recomendaciones

Para iterar y mejorar sobre los resultados actuales, se proponen las siguientes líneas de acción:

1. Exploración de nuevos modelos: investigar y evaluar de forma iterativa otros algoritmos de regresión y combinaciones de ensamblado para identificar posibles mejoras en la precisión.
2. Especialización de modelos: desarrollar un sistema que dirija las predicciones a modelos especializados. Por ejemplo, entrenar un modelo específico para cada especie podría capturar mejor las relaciones únicas entre dimensiones y peso, solucionando parte del problema del desbalance.
3. Técnicas para datos desbalanceados: investigar los efectos de aplicar técnicas específicas para regresión en datos desbalanceados, como SMOTER, para generar datos sintéticos de las categorías infrarrepresentadas y mejorar el rendimiento.
4. Recolección de datos estratégica: si es posible ampliar la recolección, priorizar peces grandes (superiores a 600g) y especies infrarepresentadas como roach, whitefish, parkki, pike y smelt. Esta es la vía más directa y efectiva para mitigar las limitaciones principales del estudio.