Modelos Lineales

Modelos Lineales

Los modelos lineales son una clase de modelos estadísticos que asumen que la relación entre las variables dependientes e independientes es lineal. En este tutorial, aprenderemos cómo ajustar un modelo lineal a un conjunto de datos y cómo interpretar los resultados.

Entre los modelos lineales más comunes se encuentran la regresión lineal simple y la regresión lineal múltiple. Se pueden utilizar para predecir una variable dependiente a partir de una o más variables independientes.

Otros modelos lineales incluyen algunas pruebas de hipótesis, como el análisis de varianza (ANOVA) y la comparación de medias (prueba t).

Todos estos modelos caen en la categoría de modelos lineales generalizados (GLM), que son una extensión de los modelos lineales tradicionales que permiten una mayor flexibilidad en la especificación de la distribución de los errores y la función de enlace.

Regresión Lineal Simple

El modelo de regresión lineal simple es uno de los modelos lineales más simples y fáciles de entender. En este modelo, se asume que la relación entre la variable dependiente y la variable independiente es lineal. La ecuación de regresión lineal simple se puede expresar de la siguiente manera:

\[Y = \beta_0 + \beta_1X + \epsilon\]

Donde:

• \(Y\) es la variable dependiente
• \(X\) es la variable independiente
• \(\beta_0\) es la intersección
• \(\beta_1\) es la pendiente
• \(\epsilon\) es el error

El objetivo de la regresión lineal simple es encontrar los valores de \(\beta_0\) y \(\beta_1\) que minimizan la suma de los cuadrados de los errores (SSE), que se define como:

\[SSE = \sum_{i=1}^{n} (Y_i - \hat{Y_i})^2\]

Donde:

• \(Y_i\) es el valor observado de la variable dependiente
• \(\hat{Y_i}\) es el valor predicho de la variable dependiente
• \(n\) es el número de observaciones

Una vez que se han estimado los valores de \(\beta_0\) y \(\beta_1\), se pueden utilizar para predecir los valores de la variable dependiente para nuevas observaciones.

Regresión Lineal Simple por Mínimos Cuadrados

El método más común para ajustar un modelo de regresión lineal simple es el método de mínimos cuadrados. Este método consiste en encontrar los valores de \(\beta_0\) y \(\beta_1\) que minimizan la suma de los cuadrados de los errores. La fórmula para calcular los valores de \(\beta_0\) y \(\beta_1\) es la siguiente:

\[\beta_1 = \frac{\sum_{i=1}^{n} (X_i - \bar{X})(Y_i - \bar{Y})}{\sum_{i=1}^{n} (X_i - \bar{X})^2}\]

\[\beta_0 = \bar{Y} - \beta_1\bar{X}\]

Donde:

• \(\bar{X}\) es la media de la variable independiente
• \(\bar{Y}\) es la media de la variable dependiente
• \(n\) es el número de observaciones

Usemos los datos de la altura y el peso de un grupo de personas para ajustar un modelo de regresión lineal simple y predecir el peso de una persona en función de su altura.

import numpy as np
import matplotlib.pyplot as plt

# Datos
altura = np.array([150, 160, 170, 180, 190])
peso = np.array([50, 60, 70, 80, 90])

# Calcular medias
altura_media = np.mean(altura)
peso_media = np.mean(peso)

# Calcular beta_1
beta_1 = np.sum((altura - altura_media) * (peso - peso_media)) / np.sum((altura - altura_media) ** 2)

# Calcular beta_0
beta_0 = peso_media - beta_1 * altura_media

# Imprimir coeficientes
print(f'beta_0: {beta_0}')
print(f'beta_1: {beta_1}')

beta_0: -100.0
beta_1: 1.0

Podemos graficar los datos y la línea de regresión para visualizar la relación entre la altura y el peso.

# Graficar datos
plt.scatter(altura, peso, color='blue')
plt.plot(altura, beta_0 + beta_1 * altura, color='red')
plt.xlabel('Altura')
plt.ylabel('Peso')
plt.title('Regresión Lineal Simple')
plt.show()

Ahora podemos crear datos datos mas realistas y ajustar un modelo de regresión lineal simple.

# Crear datos
np.random.seed(0)

altura = np.random.normal(170, 10, 100)
peso = 50 + 0.05 * altura + np.random.normal(0, 5, 100)

# Calcular medias
altura_media = np.mean(altura)
peso_media = np.mean(peso)

# Calcular beta_1
beta_1 = np.sum((altura - altura_media) * (peso - peso_media)) / np.sum((altura - altura_media) ** 2)

# Calcular beta_0
beta_0 = peso_media - beta_1 * altura_media

# Imprimir coeficientes
print(f'beta_0: {beta_0}')
print(f'beta_1: {beta_1}')

# Graficar datos
plt.figure(figsize=(10, 6))
plt.scatter(altura, peso, color='blue')
plt.plot(altura, beta_0 + beta_1 * altura, color='red')
plt.xlabel('Altura')
plt.ylabel('Peso')
plt.title('Regresión Lineal Simple')
plt.show()

beta_0: 40.626398573820204
beta_1: 0.10734921677319048

Usando seaborn podemos ajustar un modelo de regresión lineal simple y visualizar los resultados.

import seaborn as sns
import pandas as pd

# Crear DataFrame
df = pd.DataFrame({'altura': altura, 'peso': peso})

# Ajustar modelo
plt.figure(figsize=(10, 6))
sns.lmplot(x='altura', y='peso', data=df)

plt.show()

<Figure size 960x576 with 0 Axes>

Nos muestra la gráfica de la regresión lineal simple y su intervalo de confianza.

La libreria statsmodels nos permite ajustar un modelo de regresión lineal simple y obtener un resumen de los resultados.

Puedes instalar la libreria usando el siguiente comando:

!pip install statsmodels

import statsmodels.api as sm

# Ajustar modelo
X = sm.add_constant(altura)
model = sm.OLS(peso, X).fit()

# Imprimir resumen
print(model.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.042
Model:                            OLS   Adj. R-squared:                  0.033
Method:                 Least Squares   F-statistic:                     4.341
Date:                Sat, 22 Jun 2024   Prob (F-statistic):             0.0398
Time:                        13:41:19   Log-Likelihood:                -305.62
No. Observations:                 100   AIC:                             615.2
Df Residuals:                      98   BIC:                             620.4
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         40.6264      8.805      4.614      0.000      23.152      58.100
x1             0.1073      0.052      2.083      0.040       0.005       0.210
==============================================================================
Omnibus:                        5.184   Durbin-Watson:                   1.995
Prob(Omnibus):                  0.075   Jarque-Bera (JB):                3.000
Skew:                           0.210   Prob(JB):                        0.223
Kurtosis:                       2.262   Cond. No.                     2.90e+03
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.9e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

Para ajustar el modelo primero agregamos una columna de unos a la matriz de características X y luego ajustamos el modelo usando la función OLS de statsmodels. Finalmente, imprimimos un resumen de los resultados del modelo.

Imprimamos el valor de X

X

array([[  1.        , 187.64052346],
       [  1.        , 174.00157208],
       [  1.        , 179.78737984],
       [  1.        , 192.40893199],
       [  1.        , 188.6755799 ],
       [  1.        , 160.2272212 ],
       [  1.        , 179.50088418],
       [  1.        , 168.48642792],
       [  1.        , 168.96781148],
       [  1.        , 174.10598502],
       [  1.        , 171.44043571],
       [  1.        , 184.54273507],
       [  1.        , 177.61037725],
       [  1.        , 171.21675016],
       [  1.        , 174.43863233],
       [  1.        , 173.33674327],
       [  1.        , 184.94079073],
       [  1.        , 167.94841736],
       [  1.        , 173.13067702],
       [  1.        , 161.45904261],
       [  1.        , 144.47010184],
       [  1.        , 176.53618595],
       [  1.        , 178.64436199],
       [  1.        , 162.5783498 ],
       [  1.        , 192.69754624],
       [  1.        , 155.45634325],
       [  1.        , 170.45758517],
       [  1.        , 168.1281615 ],
       [  1.        , 185.32779214],
       [  1.        , 184.6935877 ],
       [  1.        , 171.54947426],
       [  1.        , 173.7816252 ],
       [  1.        , 161.12214252],
       [  1.        , 150.19203532],
       [  1.        , 166.52087851],
       [  1.        , 171.56348969],
       [  1.        , 182.30290681],
       [  1.        , 182.02379849],
       [  1.        , 166.12673183],
       [  1.        , 166.97697249],
       [  1.        , 159.51447035],
       [  1.        , 155.79982063],
       [  1.        , 152.93729809],
       [  1.        , 189.50775395],
       [  1.        , 164.90347818],
       [  1.        , 165.61925698],
       [  1.        , 157.4720464 ],
       [  1.        , 177.77490356],
       [  1.        , 153.86102152],
       [  1.        , 167.8725972 ],
       [  1.        , 161.04533439],
       [  1.        , 173.86902498],
       [  1.        , 164.89194862],
       [  1.        , 158.19367816],
       [  1.        , 169.71817772],
       [  1.        , 174.28331871],
       [  1.        , 170.66517222],
       [  1.        , 173.02471898],
       [  1.        , 163.65677906],
       [  1.        , 166.37258834],
       [  1.        , 163.27539552],
       [  1.        , 166.40446838],
       [  1.        , 161.86853718],
       [  1.        , 152.73717398],
       [  1.        , 171.77426142],
       [  1.        , 165.98219064],
       [  1.        , 153.69801653],
       [  1.        , 174.62782256],
       [  1.        , 160.92701636],
       [  1.        , 170.51945396],
       [  1.        , 177.29090562],
       [  1.        , 171.28982911],
       [  1.        , 181.39400685],
       [  1.        , 157.6517418 ],
       [  1.        , 174.02341641],
       [  1.        , 163.15189909],
       [  1.        , 161.29202851],
       [  1.        , 164.21150335],
       [  1.        , 166.88447468],
       [  1.        , 170.56165342],
       [  1.        , 158.34850159],
       [  1.        , 179.00826487],
       [  1.        , 174.6566244 ],
       [  1.        , 154.63756314],
       [  1.        , 184.88252194],
       [  1.        , 188.95889176],
       [  1.        , 181.78779571],
       [  1.        , 168.20075164],
       [  1.        , 159.29247378],
       [  1.        , 180.54451727],
       [  1.        , 165.96823053],
       [  1.        , 182.2244507 ],
       [  1.        , 172.08274978],
       [  1.        , 179.76639036],
       [  1.        , 173.56366397],
       [  1.        , 177.06573168],
       [  1.        , 170.10500021],
       [  1.        , 187.85870494],
       [  1.        , 171.26912093],
       [  1.        , 174.01989363]])

Que corresponde a lo que se denomina como matriz de diseño. Que es una matriz que contiene las variables independientes y una columna de unos que representa la intersección. Podemos usar esta matriz para ajustar el modelo de regresión lineal simple por medio de la siguiente fórmula:

\[Y = X\beta + \epsilon\]

Y usando el método de mínimos cuadrados para encontrar los valores de \(\beta\) que minimizan la suma de los cuadrados de los errores, nos da la siguiente fórmula:

\[\beta = (X^TX)^{-1}X^TY\]

Donde:

• \(Y\) es el vector de la variable dependiente
• \(X\) es la matriz de diseño
• \(\beta\) es el vector de coeficientes

Hagamos esto en Python.

# Ajustar modelo
X = np.column_stack((np.ones(len(altura)), altura))

print(X)

beta = np.linalg.inv(X.T @ X) @ X.T @ peso

# Imprimir coeficientes
print(f'beta_0: {beta[0]}')
print(f'beta_1: {beta[1]}')

[[  1.         187.64052346]
 [  1.         174.00157208]
 [  1.         179.78737984]
 [  1.         192.40893199]
 [  1.         188.6755799 ]
 [  1.         160.2272212 ]
 [  1.         179.50088418]
 [  1.         168.48642792]
 [  1.         168.96781148]
 [  1.         174.10598502]
 [  1.         171.44043571]
 [  1.         184.54273507]
 [  1.         177.61037725]
 [  1.         171.21675016]
 [  1.         174.43863233]
 [  1.         173.33674327]
 [  1.         184.94079073]
 [  1.         167.94841736]
 [  1.         173.13067702]
 [  1.         161.45904261]
 [  1.         144.47010184]
 [  1.         176.53618595]
 [  1.         178.64436199]
 [  1.         162.5783498 ]
 [  1.         192.69754624]
 [  1.         155.45634325]
 [  1.         170.45758517]
 [  1.         168.1281615 ]
 [  1.         185.32779214]
 [  1.         184.6935877 ]
 [  1.         171.54947426]
 [  1.         173.7816252 ]
 [  1.         161.12214252]
 [  1.         150.19203532]
 [  1.         166.52087851]
 [  1.         171.56348969]
 [  1.         182.30290681]
 [  1.         182.02379849]
 [  1.         166.12673183]
 [  1.         166.97697249]
 [  1.         159.51447035]
 [  1.         155.79982063]
 [  1.         152.93729809]
 [  1.         189.50775395]
 [  1.         164.90347818]
 [  1.         165.61925698]
 [  1.         157.4720464 ]
 [  1.         177.77490356]
 [  1.         153.86102152]
 [  1.         167.8725972 ]
 [  1.         161.04533439]
 [  1.         173.86902498]
 [  1.         164.89194862]
 [  1.         158.19367816]
 [  1.         169.71817772]
 [  1.         174.28331871]
 [  1.         170.66517222]
 [  1.         173.02471898]
 [  1.         163.65677906]
 [  1.         166.37258834]
 [  1.         163.27539552]
 [  1.         166.40446838]
 [  1.         161.86853718]
 [  1.         152.73717398]
 [  1.         171.77426142]
 [  1.         165.98219064]
 [  1.         153.69801653]
 [  1.         174.62782256]
 [  1.         160.92701636]
 [  1.         170.51945396]
 [  1.         177.29090562]
 [  1.         171.28982911]
 [  1.         181.39400685]
 [  1.         157.6517418 ]
 [  1.         174.02341641]
 [  1.         163.15189909]
 [  1.         161.29202851]
 [  1.         164.21150335]
 [  1.         166.88447468]
 [  1.         170.56165342]
 [  1.         158.34850159]
 [  1.         179.00826487]
 [  1.         174.6566244 ]
 [  1.         154.63756314]
 [  1.         184.88252194]
 [  1.         188.95889176]
 [  1.         181.78779571]
 [  1.         168.20075164]
 [  1.         159.29247378]
 [  1.         180.54451727]
 [  1.         165.96823053]
 [  1.         182.2244507 ]
 [  1.         172.08274978]
 [  1.         179.76639036]
 [  1.         173.56366397]
 [  1.         177.06573168]
 [  1.         170.10500021]
 [  1.         187.85870494]
 [  1.         171.26912093]
 [  1.         174.01989363]]
beta_0: 40.62639857381858
beta_1: 0.10734921677318345

Demostración de la Fórmula de Mínimos Cuadrados*

Tenemos el modelo de regresión lineal general:

\[Y = \beta_0 + \beta_1X + ... + \beta_nX_n + \epsilon\]

La función de costo que queremos minimizar es el error cuadrático residual (RSS), que se define como:

\[RSS = \sum_{i=1}^{n} (Y_i - \hat{Y_i})^2\]

En forma matricial, la función de costo se puede expresar como:

\[RSS = (Y - X\beta)^T(Y - X\beta)\]

Que se conoce como la forma cuadrática de la función de costo. Para minimizar la función de costo, tomamos la derivada de la función de costo con respecto a \(\beta\):

\[\begin{align*} \frac{\partial RSS}{\partial \beta} & = \frac{\partial}{\partial \beta} (Y - X\beta)^T(Y - X\beta) \\ & = \frac{\partial}{\partial \beta} (Y^TY - Y^TX\beta - (X\beta)^TY + (X\beta)^TX\beta) \\ & = \frac{\partial}{\partial \beta} (Y^TY - Y^TX\beta - \beta^TX^TY + \beta^TX^TX\beta) \\ & = \frac{\partial}{\partial \beta} Y^TY - \frac{\partial}{\partial \beta} Y^TX\beta - \frac{\partial}{\partial \beta} \beta^TX^TY + \frac{\partial}{\partial \beta} \beta^TX^TX\beta \\ & = 0 - X^TY - X^TY + 2X^TX\beta \\ & = -2X^TY + 2X^TX\beta \\ & = 2X^T(X\beta - Y) \\ \end{align*}\]

Igualamos la derivada de la función de costo a cero para encontrar el mínimo:

\[\frac{\partial RSS}{\partial \beta} = 0\]

\[2X^T(X\beta - Y) = 0\]

Resolviendo para \(\beta\):

\[\begin{align*} 2X^T(X\beta - Y) & = 0 \\ X^T(X\beta - Y) & = 0 \\ X^TX\beta - X^TY & = 0 \\ X^TX\beta & = X^TY \\ \beta & = (X^TX)^{-1}X^TY \\ \end{align*}\]

Regresión Lineal Simple por Descenso del Gradiente

Otra forma de ajustar un modelo de regresión lineal simple es mediante el descenso del gradiente. Este método consiste en ajustar los valores de \(\beta_0\) y \(\beta_1\) iterativamente para minimizar una función de costo. La función de costo que se utiliza en la regresión lineal simple es el error cuadrático medio (MSE), que se define como:

\[MSE = \frac{1}{n} \sum_{i=1}^{n} (Y_i - \hat{Y_i})^2\]

Donde:

• \(Y_i\) es el valor observado de la variable dependiente
• \(\hat{Y_i}\) es el valor predicho de la variable dependiente
• \(n\) es el número de observaciones

El descenso de gradiente utiliza la derivada de la función de costo con respecto a los parámetros \(\beta_0\) y \(\beta_1\) para ajustar los valores de los parámetros en la dirección que minimiza la función de costo. La regla de actualización de los parámetros es la siguiente:

\[\beta_0 = \beta_0 - \alpha \frac{1}{n} \sum_{i=1}^{n} (Y_i - \hat{Y_i})\]

\[\beta_1 = \beta_1 - \alpha \frac{1}{n} \sum_{i=1}^{n} (Y_i -\hat{Y_i})X_i\]

Donde:

• \(\alpha\) es la tasa de aprendizaje
• \(\hat{Y_i}\) es el valor predicho de la variable dependiente
• \(Y_i\) es el valor observado de la variable dependiente
• \(X_i\) es el valor de la variable independiente
• \(n\) es el número de observaciones
• \(\alpha\) es la tasa de aprendizaje que controla el tamaño de los pasos de actualización de los parámetros

Derivación de las Ecuaciones de Descenso del Gradiente*

Aquí está la demostración de la fórmula de actualización de los parámetros en el descenso del gradiente para la regresión lineal simple.

La función de costo para la regresión lineal simple es el error cuadrático medio (MSE), que se define como:

\[MSE = \frac{1}{n} \sum_{i=1}^{n} (Y_i - \hat{Y_i})^2\]

Reemplazando \(\hat{Y_i}\) con la ecuación de regresión lineal simple, obtenemos:

\[MSE = \frac{1}{n} \sum_{i=1}^{n} (Y_i - \beta_0 - \beta_1X_i)^2\]

Derivamos respecto a \(\beta_0\) y \(\beta_1\) para obtener las derivadas parciales de la función de costo:

\[\frac{\partial MSE}{\partial \beta_0} = \frac{-2}{n} \sum_{i=1}^{n} (Y_i - \beta_0 - \beta_1X_i)\]

\[\frac{\partial MSE}{\partial \beta_1} = \frac{-2}{n} \sum_{i=1}^{n} X_i (Y_i - \beta_0 - \beta_1X_i)\]

Usamos las derivadas parciales para actualizar los parámetros en la dirección que minimiza la función de costo:

\[\beta_0 = \beta_0 - \alpha \frac{\partial MSE}{\partial \beta_0}\]

\[\beta_1 = \beta_1 - \alpha \frac{\partial MSE}{\partial \beta_1}\]

Reemplazando las derivadas parciales, obtenemos la fórmula de actualización de los parámetros en el descenso del gradiente para la regresión lineal simple:

\[\beta_0 = \beta_0 + \alpha \frac{2}{n} \sum_{i=1}^{n} (Y_i - \hat{Y_i})\]

\[\beta_1 = \beta_1 + \alpha \frac{2}{n} \sum_{i=1}^{n} X_i (Y _i - \hat{Y_i})\]

El vector de coeficientes \(\beta\) se puede actualizar iterativamente hasta que se alcance un criterio de convergencia, como un número máximo de iteraciones o una tolerancia en el cambio de los parámetros.

Implementación del Descenso del Gradiente

Vamos a implementar el descenso del gradiente para ajustar un modelo de regresión lineal simple.

# Normalización de los datos
X = (altura - np.mean(altura)) / np.std(altura)
Y = (peso - np.mean(peso)) / np.std(peso)
n = len(X)

# Valores iniciales de los parámetros
b0 = 0
b1 = 0

# Hiperparámetros
alpha = 1e-3
epochs = 10000

# Descenso del gradiente
for i in range(epochs):
    y_pred = b0 + b1 * X # Predicción

    db0 = (-2/n) * sum(Y - y_pred) # Derivada parcial de b0
    db1 = (-2/n) * sum(X * (Y - y_pred)) # Derivada parcial de b1

    b0 = b0 - alpha * db0 # Actualizar b0
    b1 = b1 - alpha * db1 # Actualizar b1

    if i % 500 == 0: # Imprimir resultados cada 500 iteraciones
        print(f"Epoch {i}: b0 = {b0 * np.std(peso) + np.mean(peso) - b1 * np.mean(altura)}, b1 = {b1 * (np.std(peso) / np.std(altura))}")

# Desnormalizar los parámetros
b1 = b1 * (np.std(peso) / np.std(altura))
b0 = b0 * np.std(peso) + np.mean(peso) - b1 * np.mean(altura)

# Imprimir coeficientes
print(f"Los valores óptimos son: b0 = {b0}, b1 = {b1}")

Epoch 0: b0 = 58.86970065046024, b1 = 0.00021469843354638102
Epoch 500: b0 = 36.69221818906457, b1 = 0.06797607549967058
Epoch 1000: b0 = 28.541743780104426, b1 = 0.09287914421609095
Epoch 1500: b0 = 25.546352702386315, b1 = 0.10203130224985668
Epoch 2000: b0 = 24.44551276862618, b1 = 0.10539482333349699
Epoch 2500: b0 = 24.040941703173104, b1 = 0.10663095518768421
Epoch 3000: b0 = 23.89225728322957, b1 = 0.1070852475565854
Epoch 3500: b0 = 23.83761408547845, b1 = 0.10725220511515218
Epoch 4000: b0 = 23.817532095303484, b1 = 0.10731356389700908
Epoch 4500: b0 = 23.81015173789428, b1 = 0.10733611393992575
Epoch 5000: b0 = 23.8074393734815, b1 = 0.10734440133449522
Epoch 5500: b0 = 23.80644254903254, b1 = 0.1073474470452729
Epoch 6000: b0 = 23.806076204828443, b1 = 0.10734856637826272
Epoch 6500: b0 = 23.80594156921029, b1 = 0.10734897774573438
Epoch 7000: b0 = 23.80589208910532, b1 = 0.10734912892791008
Epoch 7500: b0 = 23.805873904609804, b1 = 0.10734918448906133
Epoch 8000: b0 = 23.805867221603073, b1 = 0.10734920490840964
Epoch 8500: b0 = 23.805864765522912, b1 = 0.10734921241274988
Epoch 9000: b0 = 23.80586386288578, b1 = 0.10734921517067947
Epoch 9500: b0 = 23.805863531156447, b1 = 0.10734921618424971
Los valores óptimos son: b0 = 40.62639861081885, b1 = 0.1073492165563143

La normalización de los datos es importante para que el descenso del gradiente converja más rápido. En este caso, normalizamos las variables independientes y dependientes restando la media y dividiendo por la desviación estándar. Después de ajustar el modelo, desnormalizamos los parámetros para obtener los valores en la escala original.

Otro factor importante es la tasa de aprendizaje \(\alpha\), que controla el tamaño de los pasos de actualización de los parámetros. Si la tasa de aprendizaje es demasiado pequeña, el algoritmo puede converger lentamente. Si es demasiado grande, el algoritmo puede divergir, encontrar la tasa de aprendizaje adecuada es crucial para el éxito del descenso del gradiente y muchos algoritmos adaptativos ajustan automáticamente la tasa de aprendizaje durante el entrenamiento para mejorar la convergencia.

Veremos de nuevo el método de descenso del gradiente cuando hablemos de redes neuronales.

Regresión Lineal Múltiple

La regresión lineal múltiple es una extensión de la regresión lineal simple que permite modelar la relación entre una variable dependiente y dos o más variables independientes. La ecuación de regresión lineal múltiple se puede expresar de la siguiente manera:

\[Y = \beta_0 + \beta_1X_1 + \beta_2X_2 + ... + \beta_nX_n + \epsilon\]

Donde:

• \(Y\) es la variable dependiente
• \(X_1, X_2, ..., X_n\) son las variables independientes
• \(\beta_0\) es la intersección
• \(\epsilon\) es el error

Al igual que la regresión lineal simple, el objetivo de la regresión lineal múltiple es encontrar los valores de \(\beta_0, \beta_1, \beta_2, ..., \beta_n\) que minimizan la suma de los cuadrados de los errores (SSE).

Regresión Lineal Múltiple por Mínimos Cuadrados

Debido a que ahora tratamos con múltiples variables independientes, la fórmula para calcular los valores de \(\beta_0, \beta_1, \beta_2, ..., \beta_n\) esta en forma matricial:

\[\beta = (X^TX)^{-1}X^TY\]

Qué es idénctico a la fórmula de la regresión lineal simple, solo que ahora \(X\) es una matriz que contiene las variables independientes y una columna de unos que representa la intersección.

Usemos datos de ejemplo para ajustar un modelo de regresión lineal múltiple y predecir la variable dependiente. Generaremos datos de 5 variables independientes y una variable dependiente.

# Crear datos
np.random.seed(1014)

# Variables independientes
X1 = np.random.normal(0, 1, 100)
X2 = np.random.normal(10, 5, 100)
X3 = np.random.normal(-5, 2, 100)
X4 = np.random.normal(3, 1, 100)
X5 = np.random.normal(2, 0.5, 100)

# Efecto de cada variable independiente
betas = np.array([5, 10, 3, -2, -1, 4])

# Error
epsilon = np.random.normal(0, 2.5, 100)

# Variable dependiente
Y = betas[0] + betas[1] * X1 + betas[2] * X2 + betas[3] * X3 + betas[4] * X4 + betas[5] * X5 + epsilon

df = pd.DataFrame({'X1': X1, 'X2': X2, 'X3': X3, 'X4': X4, 'X5': X5, 'Y': Y})

df.head()

	X1	X2	X3	X4	X5	Y
0	0.759433	10.731036	-6.699854	2.615701	1.670699	61.504776
1	-1.007260	16.029547	-3.969604	1.957694	2.076248	55.376552
2	-0.644990	10.329352	-7.744539	2.040599	2.300200	52.803668
3	-0.266741	14.039842	-7.208872	3.020797	2.237179	62.968653
4	0.291256	13.169584	-5.827498	3.217659	1.341268	61.591722

Ahora ajustaremos un modelo de regresión lineal múltiple a los datos y obtendremos un resumen de los resultados.

# Ajustar modelo
X = sm.add_constant(df[['X1', 'X2', 'X3', 'X4', 'X5']])
model = sm.OLS(df['Y'], X).fit()

# Imprimir resumen
print(model.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      Y   R-squared:                       0.976
Model:                            OLS   Adj. R-squared:                  0.975
Method:                 Least Squares   F-statistic:                     758.3
Date:                Sat, 22 Jun 2024   Prob (F-statistic):           2.63e-74
Time:                        13:41:20   Log-Likelihood:                -243.14
No. Observations:                 100   AIC:                             498.3
Df Residuals:                      94   BIC:                             513.9
Df Model:                           5                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.4953      1.737      3.165      0.002       2.047       8.943
X1             9.9913      0.310     32.255      0.000       9.376      10.606
X2             2.8923      0.058     49.999      0.000       2.777       3.007
X3            -1.9544      0.142    -13.789      0.000      -2.236      -1.673
X4            -0.7991      0.281     -2.843      0.005      -1.357      -0.241
X5             4.0200      0.613      6.553      0.000       2.802       5.238
==============================================================================
Omnibus:                        2.713   Durbin-Watson:                   2.227
Prob(Omnibus):                  0.258   Jarque-Bera (JB):                1.890
Skew:                           0.134   Prob(JB):                        0.389
Kurtosis:                       2.383   Cond. No.                         73.8
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Podemos imprimir los coeficientes del modelo para ver cómo se relacionan con los valores reales.

# Imprimir coeficientes
print(f'Intercepto: {model.params[0]}')
print(f'Coeficientes: {model.params[1:]}')

Intercepto: 5.495327623636475
Coeficientes: X1    9.991312
X2    2.892250
X3   -1.954408
X4   -0.799148
X5    4.020035
dtype: float64

/tmp/ipykernel_11165/3979803380.py:2: FutureWarning:

Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`

El resumen del modelo nos proporciona información sobre la calidad del ajuste, los coeficientes de las variables independientes, los errores estándar de los coeficientes, los valores p, el coeficiente de determinación \(R^2\), entre otros.

Gráfiquemos los parámetros estimados y los valores reales.

# Graficar parámetros estimados y valores reales
plt.figure(figsize=(10, 6))
plt.plot(betas, label='Real', marker='o', markersize=10)
plt.plot(model.params, label='Estimado', marker='x', markersize=10)
plt.xlabel('Variables', fontsize=14)
plt.ylabel('Coeficientes', fontsize=14)
plt.title('Coeficientes Estimados vs. Coeficientes Reales', fontsize=16, fontweight='bold')
plt.legend()
plt.show()

La paquetería statsmodels es muy completa y nos permite ajustar modelos de regresión lineal múltiple con facilidad. También podemos usar la función ols de statsmodels.formula.api para ajustar modelos de regresión lineal múltiple con fórmulas de estilo R.

import statsmodels.formula.api as smf

# Ajustar modelo
model = smf.ols('Y ~ X1 + X2 + X3 + X4 + X5', data=df).fit()

# Imprimir resumen
print(model.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      Y   R-squared:                       0.976
Model:                            OLS   Adj. R-squared:                  0.975
Method:                 Least Squares   F-statistic:                     758.3
Date:                Sat, 22 Jun 2024   Prob (F-statistic):           2.63e-74
Time:                        13:41:20   Log-Likelihood:                -243.14
No. Observations:                 100   AIC:                             498.3
Df Residuals:                      94   BIC:                             513.9
Df Model:                           5                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      5.4953      1.737      3.165      0.002       2.047       8.943
X1             9.9913      0.310     32.255      0.000       9.376      10.606
X2             2.8923      0.058     49.999      0.000       2.777       3.007
X3            -1.9544      0.142    -13.789      0.000      -2.236      -1.673
X4            -0.7991      0.281     -2.843      0.005      -1.357      -0.241
X5             4.0200      0.613      6.553      0.000       2.802       5.238
==============================================================================
Omnibus:                        2.713   Durbin-Watson:                   2.227
Prob(Omnibus):                  0.258   Jarque-Bera (JB):                1.890
Skew:                           0.134   Prob(JB):                        0.389
Kurtosis:                       2.383   Cond. No.                         73.8
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Dos resultados importantes del resumen es la prueba Durbin-Watson y la prueba de Jarque-Bera. La prueba Durbin-Watson se utiliza para detectar la presencia de autocorrelación en los residuos del modelo. Un valor de Durbin-Watson cercano a 2 indica que no hay autocorrelación. La prueba de Jarque-Bera se utiliza para detectar la normalidad de los residuos del modelo. Un valor de Jarque-Bera cercano a 0 indica que los residuos son normales.

Detectar la presencia de autocorrelación y no normalidad en los residuos es importante para evaluar la calidad del ajuste del modelo y tomar decisiones informadas sobre su uso. Si se obtienen resultados significativos en estas pruebas, es posible que sea necesario realizar ajustes adicionales al modelo para mejorar su rendimiento.

Regresión Lineal con Scikit-Learn

La librería scikit-learn es una de las librerías más populares para el aprendizaje automático en Python. Proporciona una amplia gama de algoritmos de aprendizaje automático, incluidos los modelos de regresión lineal.

Vamos a ajustar un modelo de regresión lineal simple y un modelo de regresión lineal múltiple con scikit-learn.

from sklearn.linear_model import LinearRegression

# Regresión Lineal Simple
X = altura.reshape(-1, 1)
Y = peso

model = LinearRegression()
model.fit(X, Y)

print(f'Intercepto: {model.intercept_}')
print(f'Coeficiente: {model.coef_[0]}')

Intercepto: 40.6263985738202
Coeficiente: 0.1073492167731905

La función reshape(-1, 1) se utiliza para convertir el vector de altura en una matriz de una sola columna. Esto es necesario porque scikit-learn espera que las variables independientes sean matrices bidimensionales. Imprimamos esta matriz.

X

array([[187.64052346],
       [174.00157208],
       [179.78737984],
       [192.40893199],
       [188.6755799 ],
       [160.2272212 ],
       [179.50088418],
       [168.48642792],
       [168.96781148],
       [174.10598502],
       [171.44043571],
       [184.54273507],
       [177.61037725],
       [171.21675016],
       [174.43863233],
       [173.33674327],
       [184.94079073],
       [167.94841736],
       [173.13067702],
       [161.45904261],
       [144.47010184],
       [176.53618595],
       [178.64436199],
       [162.5783498 ],
       [192.69754624],
       [155.45634325],
       [170.45758517],
       [168.1281615 ],
       [185.32779214],
       [184.6935877 ],
       [171.54947426],
       [173.7816252 ],
       [161.12214252],
       [150.19203532],
       [166.52087851],
       [171.56348969],
       [182.30290681],
       [182.02379849],
       [166.12673183],
       [166.97697249],
       [159.51447035],
       [155.79982063],
       [152.93729809],
       [189.50775395],
       [164.90347818],
       [165.61925698],
       [157.4720464 ],
       [177.77490356],
       [153.86102152],
       [167.8725972 ],
       [161.04533439],
       [173.86902498],
       [164.89194862],
       [158.19367816],
       [169.71817772],
       [174.28331871],
       [170.66517222],
       [173.02471898],
       [163.65677906],
       [166.37258834],
       [163.27539552],
       [166.40446838],
       [161.86853718],
       [152.73717398],
       [171.77426142],
       [165.98219064],
       [153.69801653],
       [174.62782256],
       [160.92701636],
       [170.51945396],
       [177.29090562],
       [171.28982911],
       [181.39400685],
       [157.6517418 ],
       [174.02341641],
       [163.15189909],
       [161.29202851],
       [164.21150335],
       [166.88447468],
       [170.56165342],
       [158.34850159],
       [179.00826487],
       [174.6566244 ],
       [154.63756314],
       [184.88252194],
       [188.95889176],
       [181.78779571],
       [168.20075164],
       [159.29247378],
       [180.54451727],
       [165.96823053],
       [182.2244507 ],
       [172.08274978],
       [179.76639036],
       [173.56366397],
       [177.06573168],
       [170.10500021],
       [187.85870494],
       [171.26912093],
       [174.01989363]])

Podemos utilizar el método predict para hacer predicciones con el modelo.

# Predicción
altura_nueva = np.array([[160], [170], [180]])
peso_pred = model.predict(altura_nueva)

print(peso_pred)

[57.80227326 58.87576543 59.94925759]

Podemos usar otras funciones de scikit-learn para realizar una regresion Rodge o Lasso.

from sklearn.linear_model import Ridge, Lasso

# Regresión Ridge
model = Ridge(alpha=1.0)
model.fit(X, Y)

print(f'Intercepto: {model.intercept_}')
print(f'Coeficiente: {model.coef_[0]}')

# Regresión Lasso
model = Lasso(alpha=1.0)
model.fit(X, Y)

print(f'Intercepto: {model.intercept_}')
print(f'Coeficiente: {model.coef_[0]}')

Intercepto: 40.62820122077731
Coeficiente: 0.10733865014222098
Intercepto: 42.30580012649057
Coeficiente: 0.09750501717175815

Pruebas de Normalidad

Las pruebas de normalidad son pruebas estadísticas que se utilizan para determinar si una muestra de datos proviene de una distribución normal. Existen varias pruebas de normalidad, como la prueba de Shapiro-Wilk, la prueba de Kolmogorov-Smirnov y la prueba de Anderson-Darling. O métodos gráficos como los histogramas y el gráfico Q-Q.

Gráfico Q-Q

El gráfico Q-Q (quantile-quantile) es una forma visual de evaluar si una muestra de datos proviene de una distribución normal. En un gráfico Q-Q, los cuantiles de la muestra se comparan con los cuantiles de una distribución normal teórica. Si los puntos en el gráfico Q-Q siguen una línea recta, entonces los datos se ajustan a una distribución normal. Evidentemente, si los puntos se alejan de la línea recta, entonces los datos no se ajustan a una distribución normal.

Vamos a generar datos de una distribución normal y una distribución no normal para comparar los gráficos Q-Q.

# Datos de una distribución normal
np.random.seed(0)
normal_data = np.random.normal(0, 1, 1000)

# Datos de una distribución no normal
non_normal_data = np.random.exponential(1, 1000)

# Gráfico Q-Q
fig, ax = plt.subplots(1, 2, figsize=(10, 6))

# Distribución normal
ax[0].set_title('Distribución Normal')
sm.qqplot(normal_data, line='s', ax=ax[0])
ax[0].set_aspect('equal')
ax[0].set_xlabel('Cuantiles Teóricos')
ax[0].set_ylabel('Cuantiles de los Datos')

# Distribución no normal
ax[1].set_title('Distribución No Normal')
sm.qqplot(non_normal_data, line='s', ax=ax[1])
ax[1].set_aspect('equal')
ax[1].set_xlabel('Cuantiles Teóricos')
ax[1].set_ylabel('Cuantiles de los Datos')

plt.show()

En el gráfico Q-Q de la distribución normal, los puntos siguen una línea recta, lo que indica que los datos se ajustan a una distribución normal. En el gráfico Q-Q de la distribución no normal, los puntos se alejan de la línea recta, lo que indica que los datos no se ajustan a una distribución normal.

Prueba de Shapiro-Wilk

La prueba de Shapiro-Wilk es una prueba estadística que se utiliza para determinar si una muestra de datos proviene de una distribución normal.

Hipótesis:

• \(H_0\): Los datos provienen de una distribución normal
• \(H_1\): Los datos no provienen de una distribución normal

Vamos a aplicar la prueba de Shapiro-Wilk a los datos de las distribuciones normal y no normal.

from scipy.stats import shapiro

# Prueba de Shapiro-Wilk
stat, p = shapiro(normal_data)
print(f'Distribución Normal: Estadístico = {stat}, p-valor = {p}')

stat, p = shapiro(non_normal_data)
print(f'Distribución No Normal: Estadístico = {stat}, p-valor = {p}')

Distribución Normal: Estadístico = 0.9985554728235057, p-valor = 0.5912267898687746
Distribución No Normal: Estadístico = 0.8336406754912049, p-valor = 5.034540538267324e-31

El valor p de los datos normales fue de \(0.5912267898687746\) y el valor p de los datos no normales fue de \(5.034540538267324e-31\). Por lo que podemos rechazar la hipótesis nula en el segundo caso.

Ejercicios

1. Carga los datos desde https://github.com/Christian-F-Badillo/Ciencia-de-datos-con-Python-de-estadistica-descriptiva-a-redes-neuronales/blob/main/data/cars.csv.

Columnas:

• Car_ID: A unique identifier for each car listing.
• Brand: The brand or manufacturer of the car (e.g., Toyota, Honda, Ford, etc.).
• Model: The model of the car (e.g., Camry, Civic, Mustang, etc.).
• Year: The manufacturing year of the car.
• Kilometers_Driven: The total kilometers driven by the car.
• Fuel_Type: The type of fuel used by the car (e.g., Petrol, Diesel, Electric, etc.).
• Transmission: The transmission type of the car (e.g., Manual, Automatic).
• Owner_Type: The number of previous owners of the car (e.g., First, Second, Third).
• Mileage: The fuel efficiency of the car in kilometers per liter.
• Engine: The engine capacity of the car in CC (Cubic Centimeters).
• Power: The maximum power output of the car in bhp (Brake Horsepower).
• Seats: The number of seats available in the car.
• Price: The selling price of the car in INR (Indian Rupees), which is the target variable to predict.

2. Revisa si existen valores nulos en el conjunto de datos y elimina las filas que los contengan.

3. Usando la función de pandas get_dummies convierte las variables categóricas en variables dummy. O en su defecto, usa la función LabelEncoder de sklearn.

4. Realiza un análisis exploratorio de los datos y describe las características de las variables.

5. Visualiza los datos utilizando gráficos de dispersión, histogramas y diagramas de caja para identificar posibles relaciones entre las variables.

6. Realiza un análisis de correlación entre las variables y visualiza los resultados.

7. Comprueba la normalidad de la variable dependiente (precio de los autos) y de las variables independientes (kilómetros recorridos, eficiencia de combustible, capacidad del motor, potencia máxima, número de asientos) utilizando gráficos Q-Q y la prueba de Shapiro-Wilk. ¿Qué puedes concluir?

8. Ajusta un modelo de regresión lineal múltiple para predecir el precio de los autos en función de las variables independientes (usando statsmodels y scikit-learn). ¿Cuál libreria prefieres y por qué? ¿Qué variables tienen un mayor impacto en el precio de los autos? ¿Qué otras variables podrían ser importantes para predecir el precio de los autos?

9. Evalúa la calidad del ajuste del modelo y realiza predicciones para un conjunto de datos de prueba de tu elección.

Modelos Lineales Generalizados.

Recordemos que el modelo de regresión lineal se puede representar de la siguiente forma:

\[ y \sim N(\mu, \sigma) \]

\[ \mu = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_p x_p \]

Donde \(y\) es la variable dependiente, \(\mu\) es la media de la distribución normal, \(\beta_0\) es el intercepto, \(\beta_1, \beta_2, \ldots, \beta_p\) son los coeficientes de las variables independientes \(x_1, x_2, \ldots, x_p\) respectivamente y \(\sigma\) es la desviación estándar de la distribución normal.

En el caso de los modelos lineales generalizados, la distribución de la variable dependiente \(y\) no necesariamente tiene que ser normal. En general, la distribución de \(y\) puede ser cualquier distribución de la familia exponencial. La forma general de un modelo lineal generalizado es la siguiente:

\[ y \sim f(g(\mu), \phi) \]

\[ g(\mu) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_p x_p \]

Donde \(y\) es la variable dependiente, \(f\) es la función de densidad de probabilidad de la distribución de \(y\), \(g\) es la función de enlace, \(\mu\) es el valor esperado de la distribución de \(y\), \(\beta_0\) es el intercepto, \(\beta_1, \beta_2, \ldots, \beta_p\) son los coeficientes de las variables independientes \(x_1, x_2, \ldots, x_p\) respectivamente y \(\phi\) es un parámetro de dispersión.

Dentro de las distribuciones de la familia exponencial nos encontramos con:

• Distribución Normal
• Distribución Binomial
• Distribución Poisson
• Distribución Gamma

En el caso de la distribución normal, la función de enlace es la identidad, es decir, \(g(\mu) = \mu\). En el caso de la distribución binomial, la función de enlace es la función logit, es decir, \(g(\mu) = \log\left(\frac{\mu}{1-\mu}\right)\). En el caso de la distribución Poisson, la función de enlace es la función log, es decir, \(g(\mu) = \log(\mu)\). En el caso de la distribución Gamma, la función de enlace es la inversa, es decir, \(g(\mu) = \frac{1}{\mu}\).

Las funciones de enlace tienen como objetivo enlazar la linealidad de los coeficientes \(\beta_0, \beta_1, \ldots, \beta_p\) y las variables independientes con el valor esperado de la distribución de \(y\) \(\mu\) de forma no lineal.

Regresión Logística

Podemos derivar el modelo de regresión logística a partir del modelo lineal que ya vimos. En el caso de la regresión logística, la variable dependiente \(y\) es binaria, es decir, \(y \in \{0, 1\}\). Por tanto podemos pensar en modelar la probabilidad de que \(y = 1\) en función de las variables independientes \(x_1, x_2, \ldots, x_p\). El complemento de la probabilidad de que \(y = 1\) es la probabilidad de que \(y = 0\).

Ahora calculemos los “odds” de que \(y = 1\):

\[ \text{odds} = \frac{P(y = 1)}{P(y = 0)} \]

Los “odds” son la razón de la probabilidad de que \(y = 1\) entre la probabilidad de que \(y = 0\). Si los “odds” son mayores a 1, entonces la probabilidad de que \(y = 1\) es mayor que la probabilidad de que $y = 0. Si los “odds” son menores a 1, entonces la probabilidad de que \(y = 0\) es mayor que la probabilidad de que \(y = 1\).

Ahora vamos a suponer que los “odds” de que \(y = 1\) son una función lineal de las variables independientes \(x_1, x_2, \ldots, x_p\):

\[\text{odds} = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_p x_p\]

Un problema que tienen los “odds” es que están desbalanceados. Por ejemplo, si los “odds” son 10, entonces la probabilidad de que \(y = 1\) es 10 veces mayor que la probabilidad de que \(y = 0\). Pero si los “odds” son 0.1, entonces la probabilidad de que \(y = 1\) es 10 veces menor que la probabilidad de que \(y = 0\). Para resolver este problema, vamos a tomar el logaritmo de los “odds” y vamos a modelar el logaritmo de los “odds” en función de las variables independientes \(x_1, x_2, \ldots, x_p\):

\[\log\left(\frac{P(y = 1)}{P(y = 0)}\right) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_p x_p\]

\[\log\left(\frac{P(y = 1)}{1 - P(y = 1)}\right) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_p x_p\]

Ahora solo tenemos que despejar \(P(y = 1)\) paso a paso:

\[\log\left(\frac{P(y = 1)}{1 - P(y = 1)}\right) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_p x_p\]

\[\frac{P(y = 1)}{1 - P(y = 1)} = e^{\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_p x_p}\]

\[P(y = 1) = (1 - P(y = 1))e^{\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_p x_p}\]

\[P(y = 1) = e^{\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_p x_p} - P(y = 1)e^{\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_p x_p}\]

\[P(y = 1) + P(y = 1)e^{\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_p x_p} = e^{\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_p x_p}\]

\[P(y = 1)(1 + e^{\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_p x_p}) = e^{\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_p x_p}\]

\[P(y = 1) = \frac{e^{\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_p x_p}}{1 + e^{\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_p x_p}}\]

\[P(y = 1) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_p x_p)}}\]

Esta última ecuación es la función logística. La función logística es una función sigmoide que toma valores entre 0 y 1. Gráficamente, la función logística se ve de la siguiente forma:

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-10, 10, 1000)
y = 1 / (1 + np.exp(-x))

plt.plot(x, y)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Función Logística')
plt.show()

Ahora veamos como se ve la función logística dandole valor a dos coeficientes \(\beta_0 = 0\) y \(\beta_1 = 1\):

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-10, 10, 1000)
b0 = 5
b1 = 10

y = 1 / (1 + np.exp(-(b0 + b1 * x)))

plt.plot(x, y, label=r'$y = \frac{1}{1 + e^{-(\beta_0 + \beta_1 x)}}$')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Función Logística')
plt.legend()
plt.show()

Algunas propiedades de la función logística son:

• La función logística es simétrica respecto al punto \((0, 0.5)\).
• La función logística es monótona creciente.
• La función logística es acotada entre 0 y 1.

Estimar los coeficientes de la regresión logística

Para estimar los coeficientes de la regresión logística, podemos utilizar el método de máxima verosimilitud. La función de verosimilitud es la probabilidad de observar los datos dados los coeficientes del modelo. La función de verosimilitud de la regresión logística es la siguiente:

\[L(\beta_0, \beta_1, \ldots, \beta_p) = \prod_{i=1}^{n} P(y_i = 1)^{y_i} (1 - P(y_i = 1))^{1 - y_i}\]

Donde \(n\) es el número de observaciones, \(y_i\) es la variable dependiente de la observación \(i\) y \(P(y_i = 1)\) es la probabilidad de que la observación \(i\) tome el valor 1. La función de verosimilitud es el producto de las probabilidades de que cada observación tome el valor 1 o 0.

Dado que los datos son independientes e idénticamente distribuidos, podemos tomar el logaritmo de la función de verosimilitud para simplificar los cálculos:

\[\log L(\beta_0, \beta_1, \ldots, \beta_p) = \sum_{i=1}^{n} y_i \log P(y_i = 1) + (1 - y_i) \log (1 - P(y_i = 1))\]

Ahora solo tenemos que maximizar la función de verosimilitud logarítmica para encontrar los coeficientes del modelo. Para maximizar la función de verosimilitud logarítmica, podemos utilizar el método de Newton-Raphson o con descenso de gradiente. En la práctica, la mayoría de los paquetes de software utilizan el método de Newton-Raphson.

Dado que la explicación de la estimación de los coeficientes de la regresión logística es un poco extensa, vamos a ver un ejemplo en Python utilizando el paquete statsmodels.

import pandas as pd
import numpy as np
import statsmodels.api as sm

# Generar datos
np.random.seed(0)
n = 1000
x1 = np.random.normal(0, 1, n)
x2 = np.random.normal(0, 1, n)

# Definir coeficientes
b0 = 5
b1 = 3
b2 = 2

# Generar variable dependiente
p = 1 / (1 + np.exp(-(b0 + b1 * x1 + b2 * x2)))

y = np.random.binomial(1, p)

# Crear DataFrame
df = pd.DataFrame({'y': y, 'x1': x1, 'x2': x2})

# Muestra de los datos
df.head(10)

	y	x1	x2
0	1	1.764052	0.555963
1	1	0.400157	0.892474
2	1	0.978738	-0.422315
3	1	2.240893	0.104714
4	1	1.867558	0.228053
5	1	-0.977278	0.201480
6	1	0.950088	0.540774
7	1	-0.151357	-1.818078
8	1	-0.103219	-0.049324
9	1	0.410599	0.239034

Ahora estimemos los coeficientes del modelo de regresión logística:

# Estimar modelo
data = sm.add_constant(df[['x1', 'x2']])
model = sm.GLM(df['y'], data, family=sm.families.Binomial())
result = model.fit()

# Mostrar resultados
print(result.summary())

                 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                      y   No. Observations:                 1000
Model:                            GLM   Df Residuals:                      997
Model Family:                Binomial   Df Model:                            2
Link Function:                  Logit   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -143.21
Date:                Sat, 22 Jun 2024   Deviance:                       286.43
Time:                        13:41:21   Pearson chi2:                 1.41e+03
No. Iterations:                     8   Pseudo R-squ. (CS):             0.2988
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.4077      0.442     12.228      0.000       4.541       6.274
x1             3.2868      0.317     10.383      0.000       2.666       3.907
x2             2.0070      0.248      8.093      0.000       1.521       2.493
==============================================================================

Como antes tenemos que añaadir una columna de unos al DataFrame para poder estimar el intercepto del modelo. Después utilizamos la función GLM de statsmodels que estima distribuciones de la familia exponencial. En este caso, estamos utilizando la distribución binomial, con el parámetro family=sm.families.Binomial(). Finalmente, utilizamos el método fit para estimar los coeficientes del modelo.

Podemos estimar la regresión logística con scikit-learn:

from sklearn.linear_model import LogisticRegression

# Estimar modelo
model = LogisticRegression()
model.fit(df[['x1', 'x2']], df['y'])

# Imprimir coeficientes
print(f'Intercepto: {model.intercept_}')
print(f'Coeficientes: {model.coef_}')

Intercepto: [4.94144067]
Coeficientes: [[2.94011   1.7773892]]

Podemos probar que tan bueno es el modelo al verificar que predice con los datos que ya tenemos y compararlo con los datos reales:

# Predicción
y_pred = model.predict(df[['x1', 'x2']])

# Comparar predicciones con datos reales
df['y_pred'] = y_pred
df.head(10)

	y	x1	x2	y_pred
0	1	1.764052	0.555963	1
1	1	0.400157	0.892474	1
2	1	0.978738	-0.422315	1
3	1	2.240893	0.104714	1
4	1	1.867558	0.228053	1
5	1	-0.977278	0.201480	1
6	1	0.950088	0.540774	1
7	1	-0.151357	-1.818078	1
8	1	-0.103219	-0.049324	1
9	1	0.410599	0.239034	1

Para comparar modelos de regresión logística, podemos utilizar métricas como la precisión, la sensibilidad, la especificidad, el valor predictivo positivo y el valor predictivo negativo. Estas métricas nos permiten evaluar la capacidad del modelo para predecir correctamente los valores positivos y negativos y se basan en la matriz de confusión.

from sklearn.metrics import confusion_matrix

# Matriz de confusión
cm = confusion_matrix(df['y'], df['y_pred'])
print(cm)

[[ 52  46]
 [ 11 891]]

La matriz de confusión es una tabla que muestra el número de verdaderos positivos, falsos positivos, verdaderos negativos y falsos negativos del modelo. A partir de la matriz de confusión, podemos calcular métricas como la precisión, la sensibilidad, la especificidad, el valor predictivo positivo y el valor predictivo negativo.

La diagonal de la matriz de confusión contiene los valores correctos, es decir, los verdaderos positivos y los verdaderos negativos. Los valores fuera de la diagonal contienen los valores incorrectos, es decir, los falsos positivos y los falsos negativos.

\[\begin{array}{|c|c|} \hline \text{Verdaderos Positivos (TP)} & \text{Falsos Positivos (FP)} \\ \hline \text{Falsos Negativos (FN)} & \text{Verdaderos Negativos (TN)} \\ \hline \end{array}\]

La precisión es la proporción de predicciones correctas entre todas las predicciones realizadas por el modelo:

\[\text{Precisión} = \frac{TP + TN}{TP + FP + FN + TN}\]

La sensibilidad es la proporción de verdaderos positivos entre todos los valores positivos reales:

\[\text{Sensibilidad} = \frac{TP}{TP + FN}\]

La especificidad es la proporción de verdaderos negativos entre todos los valores negativos reales:

\[\text{Especificidad} = \frac{TN}{TN + FP}\]

El valor predictivo positivo es la proporción de verdaderos positivos entre todas las predicciones positivas realizadas por el modelo:

\[\text{Valor Predictivo Positivo} = \frac{TP}{TP + FP}\]

El valor predictivo negativo es la proporción de verdaderos negativos entre todas las predicciones negativas realizadas por el modelo:

\[\text{Valor Predictivo Negativo} = \frac{TN}{TN + FN}\]

En nuestro modelo de regresión logística tenemos los siguientes valores:

# Precisión
precision = (cm[0, 0] + cm[1, 1]) / np.sum(cm)
print(f'Precisión: {precision}')

# Especificidad
specificity = cm[1, 1] / (cm[1, 1] + cm[0, 1])
print(f'Especificidad: {specificity}')

# Sensibilidad
sensitivity = cm[0, 0] / (cm[0, 0] + cm[1, 0])
print(f'Sensibilidad: {sensitivity}')

# Valor Predictivo Positivo
ppv = cm[0, 0] / (cm[0, 0] + cm[0, 1])
print(f'Valor Predictivo Positivo: {ppv}')

# Valor Predictivo Negativo
npv = cm[1, 1] / (cm[1, 1] + cm[1, 0])
print(f'Valor Predictivo Negativo: {npv}')

Precisión: 0.943
Especificidad: 0.9509071504802561
Sensibilidad: 0.8253968253968254
Valor Predictivo Positivo: 0.5306122448979592
Valor Predictivo Negativo: 0.9878048780487805

Estas métricas son usadas también en la evaluación de modelos de clasificación como lo son las redes neuronales.

Otras métricas comunes para evaluar modelos de clasificación son el área bajo la curva ROC (AUC-ROC) y el área bajo la curva PR (AUC-PR). El AUC-ROC mide la capacidad del modelo para distinguir entre las clases positiva y negativa, mientras que el AUC-PR mide la precisión del modelo en la clase positiva.

from sklearn.metrics import roc_auc_score, average_precision_score

# Área bajo la curva ROC
roc_auc = roc_auc_score(df['y'], df['y_pred'])

# Área bajo la curva PR
pr_auc = average_precision_score(df['y'], df['y_pred'])

print(f'Área bajo la curva ROC: {roc_auc}')
print(f'Área bajo la curva PR: {pr_auc}')

Área bajo la curva ROC: 0.7592085614733698
Área bajo la curva PR: 0.9503107218158627

El AUC-ROC y el AUC-PR son valores entre 0 y 1, donde un valor de 1 indica un modelo perfecto y un valor de 0.5 indica un modelo aleatorio.

Podemos graficar la curva ROC y la curva PR para visualizar la calidad del modelo de regresión logística.

from sklearn.metrics import roc_curve, precision_recall_curve

# Curva ROC
fpr, tpr, _ = roc_curve(df['y'], df['y_pred'])

# Curva PR
precision, recall, _ = precision_recall_curve(df['y'], df['y_pred'])

# Graficar curva ROC
fig, ax = plt.subplots(1, 2, figsize=(10, 6))

ax[0].plot(fpr, tpr)
ax[0].plot([0, 1], [0, 1], linestyle='--', color='gray')
ax[0].set_xlabel('Tasa de Falsos Positivos')
ax[0].set_ylabel('Tasa de Verdaderos Positivos')
ax[0].set_title('Curva ROC')

# Graficar curva PR
ax[1].plot(recall, precision)
ax[1].set_xlabel('Recall')
ax[1].set_ylabel('Precision')
ax[1].set_title('Curva PR')

plt.show()

La forma de la curva ROC y la curva PR nos permite evaluar la calidad del modelo de regresión logística. Una curva ROC que se acerca al punto (0, 1) y una curva PR que se acerca al punto (1, 1) indican un modelo de alta calidad. Cada punto en la curva ROC representa una combinación de la tasa de verdaderos positivos y la tasa de falsos positivos del modelo, mientras que cada punto en la curva PR representa una combinación de la precisión y el recall del modelo. El recall es la sensibilidad del modelo, es decir, la proporción de verdaderos positivos entre todos los valores positivos reales.