Regresión con Statsmodels

Introducción

La implementación manual construida en el módulo anterior estima los parámetros pero no provee información sobre su significancia estadística. Statsmodels aborda la regresión como un problema de inferencia: además de los estimadores, entrega toda la información necesaria para validar el modelo bajo los supuestos de OLS.

import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt

df = pd.read_csv("data/csvs/examenes.csv").sample(100, random_state=1)

OLS con statsmodels

El flujo básico es: definir las variables, crear el modelo y ajustarlo:

# Variable independiente (matriz X)
X = np.array(df[["study_hours"]])
# Para agregar intercepto explícito:
# X = sm.add_constant(X)

# Variable dependiente
y = np.array(df["exam_score"])

# Crear y ajustar el modelo
model  = sm.OLS(y, X)
result = model.fit()

# Ver el reporte completo
print(result.summary())

El output de result.summary() contiene toda la información del modelo. A continuación se explica cada sección.

Interpretación del summary

El summary se divide en tres bloques. Este es un ejemplo real con el dataset de exámenes:

                    OLS Regression Results
===========================================================
Dep. Variable:              y   R-squared:          0.886
Model:                    OLS   Adj. R-squared:     0.885
Method:          Least Squares   F-statistic:        768.8
Date:          Sat, 28 Feb 2026   Prob (F-statistic): 1.82e-48
===========================================================
                coef   std err       t    P>|t|
-----------------------------------------------------------
x1           13.5268     0.488   27.728    0.000
===========================================================
Omnibus:             1.114   Durbin-Watson:       1.549
Prob(Omnibus):       0.573   Jarque-Bera (JB):    1.094
===========================================================

R² y ajuste global

Coeficientes

P-valor y significancia estadística

El p-valor (P>|t|) responde: si el coeficiente fuera realmente cero, ¿qué tan probable sería observar este valor o uno más extremo por azar?

Análisis de residuos

El residuo de cada observación es la diferencia entre el valor real y el predicho por el modelo:

Los residuos son clave para diagnosticar si el modelo es adecuado.

Supuestos del modelo lineal

La regresión lineal tiene supuestos que deben cumplirse para que la inferencia sea válida:

1. Linealidad: la relación entre x e y es lineal
2. Independencia: los residuos son independientes entre sí
3. Homocedasticidad: la varianza de los residuos es constante para todos los valores de x
4. Normalidad: los residuos siguen una distribución aproximadamente normal

Gráficas de diagnóstico

X = np.array(df[["study_hours"]])
y = np.array(df["exam_score"])
result = sm.OLS(y, X).fit()

# Calcular residuos
y_hat    = result.predict(X)
residuos = y - y_hat

fig, axes = plt.subplots(1, 2, figsize=(11, 4))

# 1. Residuos vs valores predichos
axes[0].scatter(y_hat, residuos, alpha=0.6, s=20)
axes[0].axhline(0, color="#6c7fea", linewidth=1.2)
axes[0].set_xlabel("Valores predichos")
axes[0].set_ylabel("Residuos")
axes[0].set_title("Residuos vs Predichos")
axes[0].grid(True, alpha=0.2)

# 2. Distribución de residuos
axes[1].hist(residuos, bins=20, edgecolor="none", alpha=0.8)
axes[1].set_xlabel("Residuo")
axes[1].set_ylabel("Frecuencia")
axes[1].set_title("Distribución de residuos")

plt.tight_layout()
plt.show()

Izquierda: residuos vs valores ajustados — idealmente los puntos deben distribuirse de forma aleatoria alrededor de cero, sin patrón. Derecha: histograma de residuos — una distribución aproximadamente normal indica que el supuesto de normalidad se cumple.

Predicción

# Predecir para un nuevo valor
nuevo_x = np.array([[5.0]])   # estudiante que estudia 5 horas
prediccion = result.predict(nuevo_x)
print(f"Calificación estimada: {prediccion[0]:.1f}")

# Visualizar el modelo ajustado
x_line = np.linspace(df["study_hours"].min(), df["study_hours"].max(), 100).reshape(-1, 1)
y_pred  = result.predict(x_line)

fig, ax = plt.subplots(figsize=(7, 5))
ax.scatter(X, y, alpha=0.5, s=20, label="Datos reales")
ax.plot(x_line, y_pred, color="#6c7fea", linewidth=2, label="Modelo OLS")
ax.set_xlabel("Horas de estudio")
ax.set_ylabel("Calificación")
ax.legend()
plt.show()

Datos reales (puntos) y modelo OLS ajustado (línea). El título muestra la ecuación y el R² — en este caso el modelo explica más del 85% de la variabilidad en las calificaciones.

X = np.array(df[["study_hours"]])
X = sm.add_constant(X)  # agrega columna de unos para el término b

result = sm.OLS(y, X).fit()
print(result.summary())
# Ahora el summary muestra dos coeficientes: const (b) y x1 (m)