Selección y regularización de modelos lineales

10. Selección y regularización de modelos lineales#

A pesar de su simplicidad, los modelos lineales tienen ventajas de interpretabilidad y con frecuencia muestran buen rendimiento en predicción

Y = β_{0} + β_{1} X_{1} + \dots + β_{p} X_{p} + ϵ .

Esta sección aborda el cambiar el ajuste de mínimos cuadrados ordinarios (OLS) tradicional por métodos alternativos que se apalancan en OLS.

¿Por que usar alternativas a OLS?

Predictibilidad: Si la relación entre las variables es aproximadamente lineal y si $n ≫ p$ (el número de observaciones es más grande que el número de variables), entonces OLS tiene poco sesgo y rinde bien en datos test. Sin embargo, si $n$ no es tan grande respecto a $p$ , puede haber mucha variabilidad en OLS, lo que resulta en sobreajuste y mal rendimiento en test.

Algunos de los problemas que surgen cuando $n \approx p$ :

Sobreajuste (Overfitting)

El error cuadrático medio (MSE) en el conjunto de entrenamiento se define como:

M S E = \frac{1}{n} \sum_{i = 1}^{n} (y_{i} - {\hat{y}}_{i})^{2}

donde:

$y_{i}$ es el valor observado,
${\hat{y}}_{i}$ es la predicción del modelo,
$n$ es el número de observaciones.

Cuando el número de variables $p$ se aproxima a $n$ , el modelo puede ajustarse casi perfectamente, lo que minimiza el MSE en el conjunto de entrenamiento, pero aumenta el error en nuevos datos (error de generalización). Este es el problema del sobreajuste.

Multicolinealidad

El problema de multicolinealidad se detecta cuando las variables explicativas están altamente correlacionadas entre sí. Un indicador común es el Factor de Inflación de la Varianza (VIF) para cada variable $j$ , que se define como:

V I F_{j} = \frac{1}{1 - R_{j}^{2}}

donde $R_{j}^{2}$ es el coeficiente de determinación de la regresión de la variable $X_{j}$ sobre todas las demás variables.

Cuando $V I F_{j}$ es grande (mayor a 10, por ejemplo), indica una alta colinealidad. Esto causa inestabilidad en las estimaciones de los coeficientes ${\hat{β}}_{j}$ , lo que puede hacer que el modelo sea sensible a pequeños cambios en los datos.

Modelo No Identificable

En regresión lineal, los coeficientes de los parámetros $β = (β_{1}, β_{2}, \dots, β_{p})$ se obtienen resolviendo el sistema de ecuaciones lineales:

y = X β + ϵ

donde:

$y$ es el vector de observaciones (de dimensión $n \times 1$ ),
$X$ es la matriz de diseño (de dimensión $n \times p$ ),
$β$ es el vector de coeficientes,
$ϵ$ es el término de error.

Para resolver $β$ , necesitamos invertir la matriz $X^{T} X$ . Si $p \geq n$ , la matriz $X^{T} X$ no es invertible, lo que significa que el sistema es no identificable y no existe una única solución para $β$ .

Si det (X^{T} X) = 0, entonces no hay solución única.

Varianza Alta

Cuando hay más variables que observaciones, el error estándar de los coeficientes de regresión se incrementa. La varianza de los coeficientes ${\hat{β}}_{j}$ está dada por:

Var ({\hat{β}}_{j}) = σ^{2} {((X^{T} X)^{- 1})}_{j j}

donde $σ^{2}$ es la varianza del error residual. Si $X^{T} X$ es mal condicionada (debido a la colinealidad o alta dimensionalidad), los elementos de la matriz inversa $(X^{T} X)^{- 1}$ pueden ser muy grandes, lo que resulta en una varianza elevada para ${\hat{β}}_{j}$ . Esto significa que los coeficientes son altamente sensibles a los cambios en los datos.

Interpretabilidad: Incluir variables independientes sin asociación con la variable dependiente en un modelo de regresión resulta en incluir complejidades innecesarias en el modelo. Si son removidas (se fijan sus coeficientes en cero), se puede obtener un modelo que es más interpretable. En OLS hay muy poca probabilidad de que se tengan coeficientes iguales a cero. Al eliminar variables estamos haciendo selección de variables o feature selection: selección de variables, contracción o regularización y reducción de dimensionalidad.

10.1. Selección de variables#

10.1.1. Selección de las mejores variables#

Algoritmo: Selección de las mejores variables
1. Sea $M_{0}$ el modelo nulo (sin predictores). Este modelo solo predice la media muestral para cada observación.
2. Para $k = 1, 2, \dots, p$ :
a. Ajusta todos $(\binom{p}{k})$ modelos que contienen exactamente $k$ predictores.
b. Selecciona el mejor de los $(\binom{p}{k})$ modelos, nómbralo $M_{k}$ . El mejor modelo se define como aquel que tiene el mejor $R R S$ (suma de residuos al cuadrado), o de manera equivalente el más alto $R^{2}$ .

3. Selecciona el mejor modelo de los $M_{0}, \dots, M_{k}$ usando el error de predicción en un conjunto de validación, $C_{p}$ , AIC, BIC o $R^{2}$ ajustado. O usa el método de validación cruzada.

Nota que este algoritmo no puede ser aplicado cuando $p$ es grande.
También puede tener problemas de estadísticos cuando $p$ es grande. Cuanto mayor sea el espacio de búsqueda, mayor será la posibilidad de encontrar modelos que se vean bien en los datos de entrenamiento, aunque podrían no tener ningún poder predictivo sobre datos futuros.
Por lo tanto, un espacio de búsqueda enorme puede dar lugar a un sobreajuste y a una alta varianza de las estimaciones de los coeficientes.
Por los motivos antes listado, los métodos stepwise o paso a paso son más atractivos.

$C_{p}$ de Mallow

C_{p} = \frac{1}{n} (R S S + 2 d {\hat{σ}}^{2})

donde $d$ es el total de parámetros usados y ${\hat{σ}}^{2}$ es la estimación de la varianza del error $ϵ$ .

Criterio AIC

Se puede usar en modelos donde el ajuste se realiza con máxima verosimilitud:

A I C = - 2 log L + 2 d

Donde $L$ es el valor máximo de la función de máxima verosimilitud del modelo.

B I C = \frac{1}{n} (R S S + log (n) d {\hat{σ}}^{2})

donde $n$ es el número de observaciones.

$R^{2}$ ajustado

R_{a d j}^{2} = 1 - \frac{R S S / (n - d - 1)}{T S S / (n - 1)}

donde $T S S$ es la suma total de cuadrados.

10.1.2. Selección stepwise hacia adelante#

*Algoritmo: Selección stepwise* hacia adelante**
1. Sea $M_{0}$ el modelo nulo (sin predictores). Este modelo solo predice la media muestral para cada observación.
2. Para $k = 0, 2, \dots, p - 1$ :
a. Considera todos los $p - k$ modelos que aumentan predictores en $M_{k}$ con un predictor adicional.
b. Elige el mejor de estos $p - k$ modelos, y nómbralo $M_{k + 1}$ . El mejor modelo se define como aquel que tiene el mejor $R R S$ , o el más alto $R^{2}$ .

3. Selecciona el mejor modelo de los $M_{0}, \dots, M_{p}$ usando el error de predicción en un conjunto de validación, $C_{p}$ , AIC, BIC o $R^{2}$ ajustado. O usa el método de validación cruzada.

10.1.3. Selección stepwise hacia atrás#

*Algoritmo: Selección stepwise* hacia atrás**
1. Sea $M_{p}$ el modelo completo, que contiene todos los $p$ predictores.
2. Para $k = p, p - 1, \dots, 1$ :
a. Considera todos los $k$ modelos que contienen todos menos uno de los predictores en en $M_{k}$ , para un total de $k - 1$ predictores.
b. Elige el mejor de estos $k$ modelos, y nómbralo $M_{k - 1}$ . El mejor modelo se define como aquel que tiene el mejor $R R S$ , o el más alto $R^{2}$ .

3. Selecciona el mejor modelo de los $M_{0}, \dots, M_{p}$ usando el error de predicción en un conjunto de validación, $C_{p}$ , AIC, BIC o $R^{2}$ ajustado. O usa el método de validación cruzada.

10.1.3.1. Ejemplo#

La superintendencia de compañias del Ecuador tiene información de estados financieros de las empresas que regula. Vamos a predecir los ingresos por ventas de las empresas grandes del 2023. En la pestaña Recursos del siguiente enlace: https://appscvsmovil.supercias.gob.ec/ranking/reporte.html

import pandas as pd
import statsmodels.api as sm

uu = "https://raw.githubusercontent.com/vmoprojs/DataLectures/refs/heads/master/superciasGrandes2023.csv"
datos = pd.read_csv(uu)
datos.drop(columns = ['Unnamed: 0'],inplace = True)

# Supongamos que 'X' es tu conjunto de características y 'y' la variable dependiente.
X = datos.loc[:,datos.columns!='ingresos_ventas']  # tus datos de características
X = X.loc[:,X.columns!='costos_ventas_prod'] 
y = datos['ingresos_ventas']/1000000  # tu variable dependiente

#import sweetviz as sv

#my_report = sv.analyze(datos)
#my_report.show_html()# altas correlaciones

Notamos que existen variables con alta correlación, se eliminan correlaciones superiores a $0.9$ .

import numpy as np

# Calcular la matriz de correlación
corr_matrix = X.corr().abs()

# Seleccionar las columnas con alta correlación (por ejemplo, > 0.9)
upper_triangle = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(bool))

# Encontrar las columnas que tienen alta correlación
high_corr_columns = [column for column in upper_triangle.columns if any(upper_triangle[column] > 0.8)]

# Eliminar las columnas con alta correlación de X
X = X.drop(columns=high_corr_columns)

# Mostrar las columnas eliminadas
print("Columnas eliminadas por alta correlación:\n", high_corr_columns)

Columnas eliminadas por alta correlación:
 ['patrimonio', 'impuesto_renta', 'utilidad_ejercicio', 'utilidad_neta', 'prueba_acida', 'end_largo_plazo', 'apalancamiento', 'apalancamiento_c_l_plazo', 'rot_activo_fijo', 'rent_ope_activo', 'roa', 'total_gastos']

start_model = sm.OLS(y, X).fit()
print(start_model.summary())

                                 OLS Regression Results                                
=======================================================================================
Dep. Variable:        ingresos_ventas   R-squared (uncentered):                   0.818
Model:                            OLS   Adj. R-squared (uncentered):              0.816
Method:                 Least Squares   F-statistic:                              425.0
Date:                Mon, 14 Apr 2025   Prob (F-statistic):                        0.00
Time:                        22:21:52   Log-Likelihood:                         -15624.
No. Observations:                2969   AIC:                                  3.131e+04
Df Residuals:                    2938   BIC:                                  3.150e+04
Df Model:                          31                                                  
Covariance Type:            nonrobust                                                  
=============================================================================================
                                coef    std err          t      P>|t|      [0.025      0.975]
---------------------------------------------------------------------------------------------
activos                    4.267e-07   1.98e-08     21.527      0.000    3.88e-07    4.66e-07
utilidad_an_imp             2.59e-07   8.59e-08      3.017      0.003    9.07e-08    4.27e-07
n_empleados                   0.0252      0.002     13.497      0.000       0.022       0.029
liquidez_corriente            0.0171      0.014      1.227      0.220      -0.010       0.044
end_activo                   -0.1281      1.138     -0.113      0.910      -2.359       2.103
end_patrimonial              -0.0379      0.089     -0.426      0.670      -0.212       0.136
end_activo_fijo           -1.679e-05   5.77e-05     -0.291      0.771      -0.000    9.64e-05
end_corto_plazo              21.5596      1.986     10.855      0.000      17.665      25.454
cobertura_interes          4.176e-06   7.97e-06      0.524      0.600   -1.14e-05    1.98e-05
apalancamiento_financiero    -0.0020      0.005     -0.403      0.687      -0.011       0.008
end_patrimonial_ct           -0.0012      0.150     -0.008      0.994      -0.295       0.292
end_patrimonial_nct           0.6346      0.577      1.100      0.271      -0.496       1.765
rot_cartera                1.821e-05   1.95e-05      0.932      0.352   -2.01e-05    5.65e-05
rot_ventas                    0.0416      0.039      1.078      0.281      -0.034       0.117
per_med_cobranza          -8.883e-07   5.97e-07     -1.488      0.137   -2.06e-06    2.83e-07
per_med_pago                 1.7e-08   1.25e-08      1.361      0.174   -7.49e-09    4.15e-08
impac_gasto_a_v              -1.5197      0.725     -2.096      0.036      -2.941      -0.098
impac_carga_finan          -107.0103     22.392     -4.779      0.000    -150.917     -63.104
rent_neta_activo             14.2754      5.624      2.539      0.011       3.249      25.302
margen_bruto                -27.5007      3.178     -8.654      0.000     -33.732     -21.270
margen_operacional         6.395e-08    1.5e-08      4.254      0.000    3.45e-08    9.34e-08
rent_neta_ventas            -32.6691      7.673     -4.258      0.000     -47.713     -17.625
rent_ope_patrimonio          -0.0029      0.012     -0.236      0.814      -0.027       0.021
roe                           0.0333      0.055      0.604      0.546      -0.075       0.141
fortaleza_patrimonial        -0.4488      0.758     -0.592      0.554      -1.934       1.037
gastos_financieros         1.278e-06    4.1e-07      3.116      0.002    4.74e-07    2.08e-06
gastos_admin_ventas        1.142e-06   7.42e-08     15.384      0.000    9.96e-07    1.29e-06
depreciaciones            -7.185e-07      5e-07     -1.437      0.151    -1.7e-06    2.62e-07
amortizaciones            -3.584e-06   6.33e-07     -5.665      0.000   -4.82e-06   -2.34e-06
deuda_total                2.197e-06   1.38e-07     15.935      0.000    1.93e-06    2.47e-06
deuda_total_c_plazo        4.409e-06   4.75e-07      9.275      0.000    3.48e-06    5.34e-06
==============================================================================
Omnibus:                     4187.902   Durbin-Watson:                   2.054
Prob(Omnibus):                  0.000   Jarque-Bera (JB):          2749251.147
Skew:                           7.826   Prob(JB):                         0.00
Kurtosis:                     151.252   Cond. No.                     3.07e+09
==============================================================================

Notes:
[1] R² is computed without centering (uncentered) since the model does not contain a constant.
[2] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[3] The condition number is large, 3.07e+09. This might indicate that there are
strong multicollinearity or other numerical problems.

Evaluación del modelo

from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
def evaluate_model(y_true, y_pred):
    mae = mean_absolute_error(y_true, y_pred)
    mse = mean_squared_error(y_true, y_pred)
    rmse = np.sqrt(mse)
    r2 = r2_score(y_true, y_pred)
    return mae, mse, rmse, r2

# Predicciones y evaluación para OLS
y_pred_ols = start_model.predict(X)
mae_ols, mse_ols, rmse_ols, r2_ols = evaluate_model(y, y_pred_ols)

# Mostrar evaluación 
print("\nOLS Métricas de evaluación:")
print(f"Mean Absolute Error (MAE): {mae_ols:.2f}")
print(f"Mean Squared Error (MSE): {mse_ols:.2f}")
print(f"Root Mean Squared Error (RMSE): {rmse_ols:.2f}")
print(f"R-squared: {r2_ols:.2f}")

OLS Métricas de evaluación:
Mean Absolute Error (MAE): 16.09
Mean Squared Error (MSE): 2179.19
Root Mean Squared Error (RMSE): 46.68
R-squared: 0.80

Modelo hacia adelante:

from sklearn.feature_selection import SequentialFeatureSelector
from sklearn.linear_model import LinearRegression
# Crear el metodo para modelo de regresión lineal
model = LinearRegression()

# Crear un selector de características secuencial
sfs = SequentialFeatureSelector(model, n_features_to_select='auto', direction='forward')

# Ajustar el modelo
sfs.fit(X, y)

# Ver las características seleccionadas
selected_features = sfs.get_support()

# Ajustar el modelo final con las características seleccionadas
X_selected = X.loc[:, selected_features]
final_model = sm.OLS(y, X_selected).fit()

# Mostrar el resumen del modelo final
print(final_model.summary())

                                 OLS Regression Results                                
=======================================================================================
Dep. Variable:        ingresos_ventas   R-squared (uncentered):                   0.811
Model:                            OLS   Adj. R-squared (uncentered):              0.810
Method:                 Least Squares   F-statistic:                              842.7
Date:                Mon, 14 Apr 2025   Prob (F-statistic):                        0.00
Time:                        22:21:56   Log-Likelihood:                         -15680.
No. Observations:                2969   AIC:                                  3.139e+04
Df Residuals:                    2954   BIC:                                  3.148e+04
Df Model:                          15                                                  
Covariance Type:            nonrobust                                                  
=========================================================================================
                            coef    std err          t      P>|t|      [0.025      0.975]
-----------------------------------------------------------------------------------------
activos                4.878e-07    1.4e-08     34.946      0.000     4.6e-07    5.15e-07
n_empleados               0.0208      0.002     12.851      0.000       0.018       0.024
end_corto_plazo          22.2718      1.836     12.130      0.000      18.672      25.872
end_patrimonial_ct       -0.0346      0.086     -0.401      0.689      -0.204       0.135
rot_cartera             1.94e-05   1.98e-05      0.978      0.328   -1.95e-05    5.83e-05
per_med_pago           1.367e-08   1.26e-08      1.080      0.280   -1.11e-08    3.85e-08
impac_gasto_a_v          -1.5230      0.731     -2.082      0.037      -2.957      -0.089
impac_carga_finan       -98.7982     21.518     -4.591      0.000    -140.990     -56.606
rent_neta_activo         17.6911      5.258      3.364      0.001       7.381      28.001
margen_bruto            -25.7358      3.143     -8.189      0.000     -31.898     -19.573
rent_neta_ventas        -35.7664      7.223     -4.952      0.000     -49.928     -21.604
fortaleza_patrimonial    -0.8198      0.719     -1.140      0.254      -2.230       0.591
gastos_admin_ventas    9.447e-07   5.34e-08     17.702      0.000     8.4e-07    1.05e-06
deuda_total            2.124e-06   1.33e-07     16.029      0.000    1.86e-06    2.38e-06
deuda_total_c_plazo     4.96e-06   4.71e-07     10.530      0.000    4.04e-06    5.88e-06
==============================================================================
Omnibus:                     4174.333   Durbin-Watson:                   2.035
Prob(Omnibus):                  0.000   Jarque-Bera (JB):          2485003.184
Skew:                           7.817   Prob(JB):                         0.00
Kurtosis:                     143.866   Cond. No.                     2.81e+09
==============================================================================

Notes:
[1] R² is computed without centering (uncentered) since the model does not contain a constant.
[2] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[3] The condition number is large, 2.81e+09. This might indicate that there are
strong multicollinearity or other numerical problems.

# Predicciones y evaluación para Adelante
y_pred_adelante = final_model.predict(X_selected)
mae_adelante, mse_adelante, rmse_adelante, r2_adelante = evaluate_model(y, y_pred_adelante)

# Mostrar evaluación 
print("\n(Adelante) Métricas de evaluación:")
print(f"Mean Absolute Error (MAE): {mae_adelante:.2f}")
print(f"Mean Squared Error (MSE): {mse_adelante:.2f}")
print(f"Root Mean Squared Error (RMSE): {rmse_adelante:.2f}")
print(f"R-squared: {r2_adelante:.2f}")

(Adelante) Métricas de evaluación:
Mean Absolute Error (MAE): 16.51
Mean Squared Error (MSE): 2263.94
Root Mean Squared Error (RMSE): 47.58
R-squared: 0.79

Modelo hacia atrás

from sklearn.feature_selection import SequentialFeatureSelector
from sklearn.linear_model import LinearRegression
# Crear el metodo para modelo de regresión lineal
model = LinearRegression()

# Crear un selector de características secuencial
sfs = SequentialFeatureSelector(model, n_features_to_select='auto', direction='backward')

# Ajustar el modelo
sfs.fit(X, y)

# Ver las características seleccionadas
selected_features = sfs.get_support()

# Ajustar el modelo final con las características seleccionadas
X_selected = X.loc[:, selected_features]
final_model = sm.OLS(y, X_selected).fit()

# Mostrar el resumen del modelo final
print(final_model.summary())

                                 OLS Regression Results                                
=======================================================================================
Dep. Variable:        ingresos_ventas   R-squared (uncentered):                   0.810
Model:                            OLS   Adj. R-squared (uncentered):              0.809
Method:                 Least Squares   F-statistic:                              787.0
Date:                Mon, 14 Apr 2025   Prob (F-statistic):                        0.00
Time:                        22:22:03   Log-Likelihood:                         -15685.
No. Observations:                2969   AIC:                                  3.140e+04
Df Residuals:                    2953   BIC:                                  3.150e+04
Df Model:                          16                                                  
Covariance Type:            nonrobust                                                  
=======================================================================================
                          coef    std err          t      P>|t|      [0.025      0.975]
---------------------------------------------------------------------------------------
activos              4.838e-07   1.39e-08     34.757      0.000    4.56e-07    5.11e-07
n_empleados             0.0209      0.002     12.916      0.000       0.018       0.024
end_patrimonial        -0.0832      0.088     -0.946      0.344      -0.256       0.089
end_corto_plazo        23.6055      1.800     13.114      0.000      20.076      27.135
end_patrimonial_ct      0.0500      0.142      0.353      0.724      -0.228       0.328
end_patrimonial_nct     0.6900      0.586      1.177      0.239      -0.460       1.840
rot_cartera          2.081e-05   1.99e-05      1.048      0.295   -1.81e-05    5.98e-05
per_med_pago         1.239e-08   1.27e-08      0.977      0.329   -1.25e-08    3.73e-08
impac_gasto_a_v        -1.5354      0.733     -2.093      0.036      -2.973      -0.097
impac_carga_finan    -102.0002     21.526     -4.738      0.000    -144.208     -59.793
margen_bruto          -25.6326      3.150     -8.138      0.000     -31.809     -19.457
rent_neta_ventas      -29.5289      6.929     -4.261      0.000     -43.116     -15.942
roe                     0.0675      0.042      1.608      0.108      -0.015       0.150
gastos_admin_ventas   9.57e-07   5.33e-08     17.964      0.000    8.53e-07    1.06e-06
deuda_total          2.121e-06   1.33e-07     15.981      0.000    1.86e-06    2.38e-06
deuda_total_c_plazo  4.953e-06   4.72e-07     10.498      0.000    4.03e-06    5.88e-06
==============================================================================
Omnibus:                     4163.035   Durbin-Watson:                   2.040
Prob(Omnibus):                  0.000   Jarque-Bera (JB):          2455536.160
Skew:                           7.777   Prob(JB):                         0.00
Kurtosis:                     143.027   Cond. No.                     2.81e+09
==============================================================================

Notes:
[1] R² is computed without centering (uncentered) since the model does not contain a constant.
[2] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[3] The condition number is large, 2.81e+09. This might indicate that there are
strong multicollinearity or other numerical problems.

# Predicciones y evaluación para Atras
y_pred_atras = final_model.predict(X_selected)
mae_atras, mse_atras, rmse_atras, r2_atras = evaluate_model(y, y_pred_atras)

# Mostrar evaluación 
print("\n(Atrás) Métricas de evaluación:")
print(f"Mean Absolute Error (MAE): {mae_atras:.2f}")
print(f"Mean Squared Error (MSE): {mse_atras:.2f}")
print(f"Root Mean Squared Error (RMSE): {rmse_atras:.2f}")
print(f"R-squared: {r2_atras:.2f}")

(Atrás) Métricas de evaluación:
Mean Absolute Error (MAE): 16.66
Mean Squared Error (MSE): 2270.38
Root Mean Squared Error (RMSE): 47.65
R-squared: 0.79

10.2. Regularización#

La regularización es un método que nos permite restringir el proceso de estimación, se usa para evitar un posible sobre ajuste del modelo (overfitting).

Una forma de hacer regularización es donde los coeficientes de variables en el modelo no sean muy grandes (regresión ridge). Otra forma es contraer los coeficientes, llegando a tener coeficientes iguales a cero (regresión lasso).

Como regla general, el segundo enfoque suele ser mejor que el primero.

Los coeficientes regularizados se obtiene usando una función de penalidad $p (α)$ para restringir el tamaño del vector de coeficientes $α = (α_{1}, \dots, α_{M})^{T}$ del predictor $f (x) = \sum_{j = 1}^{M} α_{j} C_{j} (x)$ . Los coeficientes penalizados son obtenidos como solución al problema de minimización:

{\hat{α}}_{λ} = \underset{α}{argmin} = {\sum_{i = 1}^{n} L (y_{i}, f (x_{i})) + λ p (α)},

donde $L$ es una función de pérdida y $λ > 0$ es un parámetro de regularización también conocido como ratio de aprendizaje.

Existen diferentes opciones para funciones de pérdida:

Exponencial:

L (y, f (x)) = e^{- y f (x)}, y \in {- 1, + 1} .

Logística:

L (y, f_{t} (x)) = log {1 + e^{- 2 y f_{t} (x)}}, y \in {- 1, + 1} .

Error cuadrático:

L (y, f_{t} (x)) = \frac{1}{2} (y - f_{t} (x))^{2}, y \in R .

Error absoluto:

L (y, f_{t} (x)) = | y - f_{t} (x) |, y \in R .

Huber

\begin{array}{r} L (y, f_{t} (x)) = {\begin{cases} \frac{1}{2} (y - f_{t} (x))^{2}, & if | y - f_{t} (x) | \leq δ, \\ δ (| y - f_{t} (x) | - δ / 2), & en otro caso. \end{cases} \end{array}

Hay dos tipos de funciones de penalidad:

$L_{2}$ : esta función de penalidad restringe la suma de cuadrados de los coeficientes,

p_{2} (α) = \sum_{j = 1}^{M} α_{j}^{2} .

Cuando $L$ combinado es una combinación convexa y usamos la pérdida de error cuadrático, el predictor de regresión penalizado óptimo es el estimador de regresión ridge.

$L_{1}$ : Los coeficientes se restringen tal que su suma de valores absolutos,

p_{1} (α) = \sum_{j = 1}^{M} | α_{j} | .

sea menor que un valor dado. La evidencia empírica sugiere que la penalización $L_{1}$ (lasso) funciona mejor cuando hay un número pequeño o mediano de coeficientes verdaderos de tamaño moderado.

10.2.1. Regresión Ridge#

En mínimos cuadrados ordinarios, las estimaciones de $β_{0}, β_{1}, \dots, β_{p}$ se obtienen minimizando

RSS = \sum_{i = 1}^{n} {(y_{i} - β_{0} - \sum_{j = 1}^{p} β_{j} x_{i j})}^{2} .

La regresión ridge es muy similar al enfoque de mínimos cuadrados, pero los coeficientes de la regresión ridge ${\hat{β}}^{R}$ son los valores que minimizan

\sum_{i = 1}^{n} {(y_{i} - β_{0} - \sum_{j = 1}^{p} β_{j} x_{i j})}^{2} + λ \sum_{j = 1}^{p} β_{j}^{2} = RSS + λ \sum_{j = 1}^{p} β_{j}^{2},

donde $λ$ es un parámetro ajustable que se determina de manera separada. La influencia de la regularización se controla con $λ$ . Valores altos de $λ$ significa más regularización y modelos más simples. Sin embargo, la regresión ridge siempre generará un modelo que incluya todos los predictores. Incrementar el valor de $λ$ tenderá a reducir las magnitudes de los coeficientes, pero no resultará en la exclusión de ninguna de las variables. Notemos que la penalidad no se aplica a $β_{0}$ .

10.2.2. Regresión Lasso#

\sum_{i = 1}^{n} {(y_{i} - β_{0} - \sum_{j = 1}^{p} β_{j} x_{i j})}^{2} + λ \sum_{j = 1}^{p} | β_{j} | = RSS + λ \sum_{j = 1}^{p} | β_{j} |,

La penalización $L_{1}$ tiene el efecto de forzar algunas de las estimaciones de coeficientes a ser exactamente iguales a cero cuando el parámetro $λ$ de ajuste es suficientemente grande. Por lo tanto, lasso realiza selección de variables. Como resultado, los modelos generados a partir de lasso son generalmente mucho más fáciles de interpretar que los producidos por regresión ridge. Decimos que lasso produce modelos dispersos (sparse), es decir, modelos que involucran solo un subconjunto de variables predictoras.

10.2.2.1. Ejemplo#

import matplotlib.pyplot as plt
from sklearn.linear_model import Ridge, Lasso
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import StandardScaler


# Dividir los datos en entrenamiento y prueba
#X = X_selected
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=85)


# Escalar los datos
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# Crear modelos de Ridge y Lasso
ridge = Ridge(alpha=1.0)  # Alpha controla la penalización
lasso = Lasso(alpha=5.0)
regOLS = LinearRegression()

# Entrenar los modelos
ridge.fit(X_train_scaled, y_train)
lasso.fit(X_train_scaled, y_train)
regOLS.fit(X_train_scaled, y_train)

# Predicciones
y_pred_ridge_tr = ridge.predict(X_train_scaled)
y_pred_lasso_tr = lasso.predict(X_train_scaled)
y_pred_ols_tr = regOLS.predict(X_train_scaled)


# Calcular el error cuadrático medio (MSE) TRAIN
num_factor = 1
mse_ridge_tr = mean_squared_error(y_train, y_pred_ridge_tr)/num_factor
mse_lasso_tr = mean_squared_error(y_train, y_pred_lasso_tr)/num_factor
mse_ols_tr = mean_squared_error(y_train, y_pred_ols_tr)/num_factor

# Resultados
print("Evaluación del modelo en TRAIN:\n")
print(f"MSE Ridge: {mse_ridge_tr:.2f}")
print(f"MSE Lasso: {mse_lasso_tr:.2f}")
print(f"MSE OLS: {mse_ols_tr:.2f}")


# Predicciones
y_pred_ridge = ridge.predict(X_test_scaled)
y_pred_lasso = lasso.predict(X_test_scaled)
y_pred_ols = regOLS.predict(X_test_scaled)

# Calcular el error cuadrático medio (MSE) TEST
mse_ridge = mean_squared_error(y_test, y_pred_ridge)/num_factor
mse_lasso = mean_squared_error(y_test, y_pred_lasso)/num_factor
mse_ols = mean_squared_error(y_test, y_pred_ols)/num_factor

# Resultados
print("\nEvaluación del modelo en TEST:\n")
print(f"MSE Ridge: {mse_ridge:.2f}")
print(f"MSE Lasso: {mse_lasso:.2f}")
print(f"MSE OLS: {mse_ols:.2f}")

Evaluación del modelo en TRAIN:

MSE Ridge: 1852.46
MSE Lasso: 2083.61
MSE OLS: 1852.46

Evaluación del modelo en TEST:

MSE Ridge: 3527.83
MSE Lasso: 3663.03
MSE OLS: 3527.39

# Mostrar coeficientes con nombres
ridge_coefficients = pd.DataFrame(ridge.coef_, index=X.columns, columns=['Ridge Coefficients'])
lasso_coefficients = pd.DataFrame(lasso.coef_, index=X.columns, columns=['Lasso Coefficients'])
ols_coefficients = pd.DataFrame(regOLS.coef_, index=X.columns, columns=['OLS Coefficients'])

# Combinar coeficientes en un solo DataFrame
coefficients = pd.concat([ridge_coefficients, lasso_coefficients,ols_coefficients], axis=1)
print("\nCoeficientes:")
print(coefficients)

Coeficientes:
                           Ridge Coefficients  Lasso Coefficients  \
activos                             46.452401           51.339290   
utilidad_an_imp                      7.140067            2.802706   
n_empleados                         20.230335           15.067053   
liquidez_corriente                   1.097188           -0.000000   
end_activo                          -0.310678            0.000000   
end_patrimonial                     -2.044002            0.000000   
end_activo_fijo                     -0.244253           -0.000000   
end_corto_plazo                      3.527861            0.000000   
cobertura_interes                    0.281633           -0.000000   
apalancamiento_financiero           -0.502044           -0.000000   
end_patrimonial_ct                   0.527377            0.000000   
end_patrimonial_nct                  1.076781            0.000000   
rot_cartera                          0.157103            0.000000   
rot_ventas                           1.052480            0.000000   
per_med_cobranza                    -0.919068           -0.000000   
per_med_pago                         1.745660            0.000000   
impac_gasto_a_v                     -2.108877           -0.000000   
impac_carga_finan                   -4.171406           -0.037732   
rent_neta_activo                     2.607665            0.000000   
margen_bruto                        -8.536271           -4.644804   
margen_operacional                   5.138157            0.000000   
rent_neta_ventas                    -4.865955           -0.924355   
rent_ope_patrimonio                 -0.019458           -0.000000   
roe                                  2.066205            0.000000   
fortaleza_patrimonial               -0.921451           -0.000000   
gastos_financieros                   3.886436            0.000000   
gastos_admin_ventas                 31.292825           23.488379   
depreciaciones                      -3.459897           -0.000000   
amortizaciones                      -6.807988           -0.000000   
deuda_total                         17.489693           14.618748   
deuda_total_c_plazo                 11.799476           10.246291   

                           OLS Coefficients  
activos                           46.520173  
utilidad_an_imp                    7.103579  
n_empleados                       20.224461  
liquidez_corriente                 1.098485  
end_activo                        -0.310745  
end_patrimonial                   -2.056202  
end_activo_fijo                   -0.244556  
end_corto_plazo                    3.530621  
cobertura_interes                  0.281992  
apalancamiento_financiero         -0.502272  
end_patrimonial_ct                 0.530975  
end_patrimonial_nct                1.078609  
rot_cartera                        0.157290  
rot_ventas                         1.052782  
per_med_cobranza                  -0.924138  
per_med_pago                       1.746044  
impac_gasto_a_v                   -2.110272  
impac_carga_finan                 -4.169610  
rent_neta_activo                   2.610887  
margen_bruto                      -8.539168  
margen_operacional                 5.150837  
rent_neta_ventas                  -4.872438  
rent_ope_patrimonio               -0.020343  
roe                                2.076751  
fortaleza_patrimonial             -0.920677  
gastos_financieros                 3.868392  
gastos_admin_ventas               31.299005  
depreciaciones                    -3.459313  
amortizaciones                    -6.817272  
deuda_total                       17.486285  
deuda_total_c_plazo               11.794246

# Graficar los resultados
plt.figure(figsize=(12, 6))
plt.scatter(y_test, y_pred_ridge, color='blue', label="Predicciones Ridge", alpha=0.5)
plt.scatter(y_test, y_pred_lasso, color='red', label="Predicciones Lasso", alpha=0.5)
plt.plot([min(y_test), max(y_test)], [min(y_test), max(y_test)], color='black', linestyle='--', label="Línea de referencia")
plt.xlabel("Valores reales")
plt.ylabel("Valores predichos")
plt.title("Comparación de Predicciones: Ridge vs Lasso")
plt.legend()
plt.show()

_images/52d712b991314b6107628fb9a14b186e4d597fef20569544cdab5c55f44ed33d.png

Ajuste con el mejor $λ$

from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import KFold

# Definir el rango de valores de alpha (lambda) a evaluar
param_grid = {'alpha': [0.1, 0.5, 0.9, 1, 2,3,4,5,10,20,30,40,50,100,150,200,300,400,500]}

# Crear el modelo de Ridge y Lasso
ridge = Ridge()
lasso = Lasso()

# Crear un KFold con shuffle y random_state
kfold = KFold(n_splits=5, shuffle=True, random_state=42)

# Implementar GridSearchCV para Ridge
ridge_search = GridSearchCV(ridge, param_grid, scoring='neg_mean_squared_error', cv=kfold) #cv=5 sería determinista
ridge_search.fit(X_train_scaled, y_train)

# Implementar GridSearchCV para Lasso
lasso_search = GridSearchCV(lasso, param_grid, scoring='neg_mean_squared_error', cv=kfold)
lasso_search.fit(X_train_scaled, y_train)

# Obtener el mejor valor de alpha y el MSE para Ridge
best_ridge_alpha = ridge_search.best_params_['alpha']
best_ridge_mse = -ridge_search.best_score_/num_factor  # Negativo para convertirlo en MSE normal

# Obtener el mejor valor de alpha y el MSE para Lasso
best_lasso_alpha = lasso_search.best_params_['alpha']
best_lasso_mse = -lasso_search.best_score_/num_factor  # Negativo para convertirlo en MSE normal

import matplotlib.pyplot as plt
import numpy as np

# Obtener los resultados de Ridge y Lasso
ridge_results = ridge_search.cv_results_
lasso_results = lasso_search.cv_results_

# Valores de alpha
alphas = ridge_results['param_alpha'].data

# MSE para Ridge
ridge_mse = -ridge_results['mean_test_score']

# MSE para Lasso
lasso_mse = -lasso_results['mean_test_score']

# Crear la figura
plt.figure(figsize=(14, 7))

# Graficar Ridge
plt.subplot(1, 2, 1)
plt.plot(alphas, ridge_mse, marker='o', label='MSE Ridge', color='blue')
plt.axvline(x=best_ridge_alpha, color='red', linestyle='--', label=f'Best alpha: {best_ridge_alpha}')
plt.xscale('log')  # Escala logarítmica para los ejes
plt.title('Ridge Regression: MSE vs. Alpha')
plt.xlabel('Alpha (log scale)')
plt.ylabel('Mean Squared Error')
plt.legend()
plt.grid(True)

# Graficar Lasso
plt.subplot(1, 2, 2)
plt.plot(alphas, lasso_mse, marker='o', label='MSE Lasso', color='green')
plt.axvline(x=best_lasso_alpha, color='red', linestyle='--', label=f'Best alpha: {best_lasso_alpha}')
plt.xscale('log')  # Escala logarítmica para los ejes
plt.title('Lasso Regression: MSE vs. Alpha')
plt.xlabel('Alpha (log scale)')
plt.ylabel('Mean Squared Error')
plt.legend()
plt.grid(True)

# Mostrar la gráfica
plt.tight_layout()
plt.show()

_images/2b9951457e0ac79b69d4d395091d79e581857330075beaaebfd7f68ed8a67329.png

# Imprimir los mejores resultados
print("\nEvaluación del modelo en CV:\n")
print(f"Mejor alpha para Ridge: {best_ridge_alpha}")
print(f"Mejor MSE para Ridge: {best_ridge_mse:.2f} ")
print(f"Mejor alpha para Lasso: {best_lasso_alpha}")
print(f"Mejor MSE para Lasso: {best_lasso_mse:.2f} ")

# Ajustar los modelos nuevamente con los mejores alpha
ridge_best = Ridge(alpha=best_ridge_alpha)
lasso_best = Lasso(alpha=best_lasso_alpha)

ridge_best.fit(X_train_scaled, y_train)
lasso_best.fit(X_train_scaled, y_train)

# Predicciones con los mejores modelos TRAIN
y_pred_ridge_best_tr = ridge_best.predict(X_train_scaled)
y_pred_lasso_best_tr = lasso_best.predict(X_train_scaled)

# Calcular el MSE en el conjunto TRAIN
mse_ridge_best_tr = mean_squared_error(y_train, y_pred_ridge_best_tr)/num_factor
mse_lasso_best_tr = mean_squared_error(y_train, y_pred_lasso_best_tr)/num_factor

# Mostrar los resultados
print("\nEvaluación del modelo en TRAIN:\n")
print(f"\nMSE Ridge con mejor alpha: {mse_ridge_best_tr:.2f}")
print(f"MSE Lasso con mejor alpha: {mse_lasso_best_tr:.2f} ")



# Predicciones con los mejores modelos
y_pred_ridge_best = ridge_best.predict(X_test_scaled)
y_pred_lasso_best = lasso_best.predict(X_test_scaled)

# Calcular el MSE en el conjunto de prueba
mse_ridge_best = mean_squared_error(y_test, y_pred_ridge_best)/num_factor
mse_lasso_best = mean_squared_error(y_test, y_pred_lasso_best)/num_factor

# Mostrar los resultados
print("\nEvaluación del modelo en TEST:\n")
print(f"\nMSE Ridge con mejor alpha: {mse_ridge_best:.2f}")
print(f"MSE Lasso con mejor alpha: {mse_lasso_best:.2f} ")

Evaluación del modelo en CV:

Mejor alpha para Ridge: 300
Mejor MSE para Ridge: 2282.32 
Mejor alpha para Lasso: 0.9
Mejor MSE para Lasso: 2374.17 

Evaluación del modelo en TRAIN:


MSE Ridge con mejor alpha: 1903.47
MSE Lasso con mejor alpha: 1871.67 

Evaluación del modelo en TEST:


MSE Ridge con mejor alpha: 3617.36
MSE Lasso con mejor alpha: 3534.08