Regresión logística
Introducción
La regresión logística es un modelo estadístico ampliamente utilizado para predecir variables binarias, donde la respuesta toma valores de 0 o 1. Su desarrollo se basa en la función logística, que transforma el predictor lineal en probabilidades entre 0 y 1 mediante la función sigmoide. Este enfoque permite modelar la relación entre múltiples covariables y la probabilidad de ocurrencia de un evento. Las aplicaciones son diversas, incluyendo diagnóstico médico, análisis de crédito, marketing, predicción de fraude y riesgo de abandono de clientes. La estimación de los parámetros se realiza típicamente mediante máxima verosimilitud, resolviendo ecuaciones no lineales iterativamente con el método de Newton-Raphson. Cada iteración ajusta los coeficientes usando la Hessiana y el gradiente de la log-verosimilitud, asegurando convergencia hacia los estimadores óptimos. A diferencia de la regresión lineal, el modelo logit maneja correctamente la no linealidad y restricciones de probabilidad. Los resultados permiten interpretar los coeficientes en términos de odds ratios, facilitando su explicación práctica. La matriz de covarianzas de los estimadores, calculada como \((𝑋^{\top}R𝑋)^{-1}\), permite obtener errores estándar y evaluar la significancia estadística. En conclusión, la regresión logística combina rigor matemático, flexibilidad y amplio potencial de aplicación en problemas de clasificación binaria.
Definición
La forma general del modelo de regresión logística es \[\begin{equation} y_i = E[y_i] + \epsilon_i \end{equation}\]
donde las observaciones \(y_i\) son variables aleatorias independientes con distribución Bernoulli, cuyos valores esperados son \[\begin{eqnarray} E[y_i] & = & \pi_i\\ & = & \frac{exp(x_i^{\top}\beta)}{1+exp(x_i^{\top}\beta)}\\ & = & \frac{1}{1+exp(-x_i^{\top}\beta)} \end{eqnarray}\]
despejando, el predictor lineal queda definido de la siguiente manera \[\begin{eqnarray} log\left(\frac{\pi}{1-\pi} \right) & = & x_i^{\top}\beta \end{eqnarray}\]
y, expresado de la siguiente manera se le conoce como razones de momios (log-odds)
\[\begin{eqnarray} \frac{\pi}{1-\pi} & = & \exp(x_i^{\top}\beta ) \end{eqnarray}\]
o, de manera equivalente en términos logarítmicos,
\[ \log \left( \frac{\pi_i}{1 - \pi_i} \right) = \alpha + \beta_1 X_{i1} + \beta_2 X_{i2} + \cdots + \beta_k X_{ik}. \]
Además, las variables \(X\) pueden ser tan generales como en el modelo lineal general, incluyendo, por ejemplo:
- Variables explicativas cuantitativas.
- Transformaciones de variables explicativas cuantitativas.
- Regresores polinomiales (o splines de regresión) formados a partir de variables explicativas cuantitativas.
- Variables indicadoras (dummy) que representan variables explicativas cualitativas.
- Regresores de interacción.
Estimación
En esta sección desarrollaré los detalles de la estimación por máxima verosimilitud para el modelo logit lineal general. Es conveniente reescribir el modelo en forma vectorial como
\[ \pi_i = \frac{1}{1 + \exp(-x_i^{\top}\beta)} \]
donde \(x_i^{\top} = (1, X_{i1}, \dots, X_{ik})\) es la fila \(i\)-ésima de la matriz del modelo \(X\), y \(\beta = (\alpha, \beta_1, \dots, \beta_k)^{\top}\) es el vector de parámetros.
Verisimilitud
Cada observación de la muestra sigue una distribución Bernoulli, por lo que la distribución de probabilidad es
\[ \begin{eqnarray} f_i & = & \pi_i^{y_i} (1-\pi_i)^{1-y_i}, \ \ i= 1,2,..., n \\ & = & \left( \frac{\pi_i}{1- \pi_i} \right)^{y_i} (1-\pi_i), \ \ i= 1,2,..., n \end{eqnarray} \] donde cada observación \(y_i\) toma valores de \(0\) o \(1\), Cuya verosimilitud
\[ p(y_1, \dots, y_n \mid X) = \prod_{i=1}^{n} \left[ \left(exp(x_i^{\top}\beta) \right)^{y_i} \left( \frac{1}{1 + \exp(x_i^{\top}\beta)} \right) \right]. \]
La función de log-verosimilitud es
\[ \ell(\beta) = \sum_{i=1}^{n} y_i x_i^{\top}\beta - \sum_{i=1}^{n} \log\left(1 + \exp(x_i^{\top}\beta)\right). \]
Primera derivada parcial
Las derivadas parciales de la log-verosimilitud con respecto a \(\beta\) son
\[ \frac{\partial \ell(\beta)}{\partial \beta} = \sum_{i=1}^{n} y_i x_i - \sum_{i=1}^{n} \left[ \frac{\exp(x_i^{\top}\beta)}{1 + \exp(x_i^{\top}\beta)} \right] x_i \]
\[ = \sum_{i=1}^{n} y_i x_i - \sum_{i=1}^{n} \left[ \frac{1}{1 + \exp(-x_i^{\top}\beta)} \right] x_i. \]
\[ = \sum_{i=1}^{n} y_i x_i - \sum_{i=1}^{n} P_i x_i \]
donde
\[ P_i = \frac{1}{1 + \exp(-x_i^{\top} \hat{\beta})} \]
en términos matriciales
\[ \frac{\partial \ell(\beta)}{\partial \beta} = X^{\top} (y-p) \]
Igualando el vector de derivadas parciales a cero para maximizar la verosimilitud, se obtienen las ecuaciones de estimación
\[ \sum_{i=1}^{n} \frac{1}{1 + \exp(-x_i^{\top} \hat{\beta})} x_i = \sum_{i=1}^{n} y_i x_i, \]
donde \(\hat{\beta} = (\hat{\alpha}, \hat{\beta}_1, \dots, \hat{\beta}_k)^{\top}\) es el vector de estimadores de máxima verosimilitud. En forma matricial, las ecuaciones de estimación pueden escribirse como
\[ X^{\top} p = X^{\top} y, \]
donde \(p = (P_1, \dots, P_n)^{\top}\) es el vector de valores ajustados. Obsérvese la similitud con las ecuaciones normales de mínimos cuadrados
\[ X^{\top} X \hat{\beta} = X^{\top} y, \]
las cuales pueden escribirse como
\[ X^{\top} \hat{y} = X^{\top} y. \]
A diferencia de las ecuaciones normales del modelo lineal, las ecuaciones de estimación del modelo logit
\[ \sum_{i=1}^{n} P_i x_i = \sum_{i=1}^{n} y_i x_i \]
son funciones no lineales de \(\hat{\beta}\) y, por lo tanto, requieren un procedimiento iterativo para su solución. Un enfoque común para resolverlas es el método de Newton-Raphson, antes vamos a obtener la derivadas parciales y segundas derivadas parciales de la función de verodomilitud.
Las derivadas parciales de la log-verosimilitud con respecto a \(\beta\) son
\[ \frac{\partial \ell(\beta)}{\partial \beta} = \sum_{i=1}^{n} y_i x_i - \sum_{i=1}^{n} \left[ \frac{\exp(x_i^{\top}\beta)}{1 + \exp(x_i^{\top}\beta)} \right] x_i \]
\[ = \sum_{i=1}^{n} y_i x_i - \sum_{i=1}^{n} \left[ \frac{1}{1 + \exp(-x_i^{\top}\beta)} \right] x_i. \]
Segunda derivada parcial
Ahora vamos a obtener la segundas derivadas parciales con respecto a \(\beta\). Sea \(H = \frac{\partial^2 \ell(\beta)}{\partial^2 \beta}\). luego
\[ H = \sum_{i=1}^{n} \left[ \frac{\exp(-x_i^{\top}\beta)}{\left[1 + \exp(- x_i^{\top}\beta)\right]^2} x_i x^{\top}_i \right]^{-1} \]
\[ = \left( \sum_{i=1}^{n} P_i (1 - P_i)\, x_i x_i^\top \right)^{-1} \]
\[ = (X^\top R X)^{-1} \]
donde \(R = diag\{P_i(1-P_i)\}\)
Método de Newton-Raphson
El método de Newton-Raphson se define como
\[ \begin{equation} x^{(k+1)} = x^{(k)} - \frac{f(x^{(k)})}{f'(x^{(k)})} \end{equation} \]
en nuestro caso se expresaría como
\[ \begin{equation} \hat{\beta}^{(\ell + 1)} = \hat{\beta}^{(\ell)} - \frac{\frac{\partial \ell(\beta)}{\partial \beta}}{\frac{\partial^2 \ell(\beta)}{\partial^2 \beta}} \end{equation} \]
Para facilitar la escritura de varuables vamos a considerar \(b = \hat{\beta}\), entonces, seguimos el siguiente procedimiento:
1.- Seleccionar valores iniciales \(b^{(0)}\); una elección simple es
\[ b^{(0)} = 0. \]
2.- En cada iteración \(\ell + 1\), calcular
\[ \boxed{ b^{(\ell+1)} = b^{(\ell)} + H^{-1} X^{\top} \left( y - p^{(\ell)} \right) } \]
donde \(p^{(\ell)} = \left(\frac{1}{1 + \exp(-x_i^{\top} b^{(\ell)})} \right)\) es el vector de valores ajustados en la iteración previa, y \(R^{(\ell)} = \mathrm{diag}\{ P_i^{(\ell)} (1 - P_i^{(\ell)}) \}\).
3.- Las iteraciones continúan hasta que \(b^{(\ell+1)} \approx b^{(\ell)}\) con el grado de precisión deseado. Cuando ocurre convergencia,
\[ H^{-1} X^{\top} \left( y - p^{(\ell)} \right) \approx 0, \]
y por lo tanto las ecuaciones de estimación
\[ X^{\top} p = X^{\top} y \]
se satisfacen aproximadamente.
Implementación
import pandas as pd
import numpy as np
# matplotlib
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import seaborn as sns
from statsmodels.formula.api import logit
import warnings
import math
warnings.filterwarnings("ignore", category=UserWarning, module="pandas")
sns.set_theme()Para la implementación vamos a considerar los siguientes datos
El modelo general \[\begin{equation} y = \frac{e^{\beta_0 + \beta_1 x_1 + \beta_2 x_2}}{1+ e^{\beta_0 + \beta_1 x_1 + \beta_2 x_2}} \end{equation}\]
Predictor lineal
\[z_n = \omega_0 + \omega_1 x_1 + \omega_2 x_2\]
vamos a simular:
\[ z_n = 4.5 - 0.1 x_1 + 1.8 x_2\]
# muestras
n = 1000
# variables
k = 2
# coeficientes
coef= np.array([4.5, -0.1, 1.8])
# generamos 2 variables x_1 y x_2
X = np.random.randn(n, k)
# agregamos el término constante
X = np.hstack((np.ones((n, 1)), X))
# predictor lineal
zn = X.dot(coef)
# calculamos las probabilidades
prob = sigmoid(zn)
# generamos la variable respuesta con la distribución binomial
y = np.random.binomial(1, prob)## [[ 1. 0.24711162 -0.80119847]
## [ 1. -0.08328256 -0.27390787]
## [ 1. 0.95835831 -0.24877155]
## ...
## [ 1. 0.76265347 1.12866174]
## [ 1. -0.81737812 -0.52785491]
## [ 1. 1.26675885 0.25756795]]## [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1
## 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1
## 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 0 1 1 0 0 1 1 1 0 1 1 1 1 1 1
## 1]Implementación de Newton-Raphson
n, k = X.shape
# inicializamos coeficientes
b = np.zeros(k)
# iniciamos iteración
for i in range(10000):
# obtenemos probabilidades
Pi = sigmoid(X.dot(b))
# matriz diagonal
R = np.diag(Pi * (1 - Pi))
nabla = X.T.dot(Pi - y) #
# Calculamos H
H = X.T.dot(R).dot(X)
# actualización de coeficientes
b_new = b - np.linalg.inv(H).dot(nabla)
b = b_new## [4.35711132 0.11966899 1.55407206]Empleando statsmodels
names = ['x0', 'x1', 'x2', 'y']
df = pd.DataFrame(pd.concat([pd.DataFrame(X),pd.Series(y)], axis=1))
df.columns = names
df = df.drop('x0', axis=1)
df.head()## x1 x2 y
## 0 0.247112 -0.801198 1
## 1 -0.083283 -0.273908 1
## 2 0.958358 -0.248772 1
## 3 0.894632 0.287636 1
## 4 1.298548 -2.044638 1## Optimization terminated successfully.
## Current function value: 0.117421
## Iterations 9## Logit Regression Results
## ==============================================================================
## Dep. Variable: y No. Observations: 1000
## Model: Logit Df Residuals: 997
## Method: MLE Df Model: 2
## Date: mié., 04 mar. 2026 Pseudo R-squ.: 0.2260
## Time: 22:06:11 Log-Likelihood: -117.42
## converged: True LL-Null: -151.71
## Covariance Type: nonrobust LLR p-value: 1.278e-15
## ==============================================================================
## coef std err z P>|z| [0.025 0.975]
## ------------------------------------------------------------------------------
## Intercept 4.3571 0.313 13.910 0.000 3.743 4.971
## x1 0.1197 0.184 0.651 0.515 -0.240 0.480
## x2 1.5541 0.212 7.316 0.000 1.138 1.970
## ==============================================================================Conclusión
En el ejemplo implementado, se simuló un conjunto de datos con dos variables predictoras y un término constante, generando probabilidades mediante la función sigmoide y observaciones binarias a partir de una distribución binomial. El algoritmo de Newton-Raphson permitió estimar los coeficientes de manera iterativa, utilizando la Hessiana y el gradiente de la log-verosimilitud, mostrando convergencia hacia los valores verdaderos de los parámetros \(𝑏≈ [4.5,−0.1,1.8]\). Los resultados obtenidos con la implementación manual se validaron correctamente mediante el modelo logit de statsmodels, confirmando la consistencia entre ambos métodos. Este ejercicio demuestra cómo la regresión logística puede estimarse desde cero y cómo los conceptos teóricos de máxima verosimilitud se aplican en la práctica.
Referencias:
Agresti Alan (2015). Foundations of linear and generalized linear models. Wiley series in probability and statistics.
Lindsey James K. (1997). Applying Generalized Linear Models. Springer-Verlag New York, Inc.