Regresión Múltiple¶

La regresión múltiple extiende el modelo simple a varios predictores simultáneos. Esto introduce nuevas preguntas: ¿el modelo en su conjunto es significativo? ¿cada variable aporta algo? ¿los predictores están correlacionados entre sí? ¿se cumplen los supuestos? Este notebook cubre la inferencia y el diagnóstico necesarios para responder estas preguntas con statsmodels.

---

Librerías¶

In [ ]:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.stats as stats
import statsmodels.formula.api as smf
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.stats.diagnostic import het_breuschpagan, het_white
from statsmodels.stats.stattools import durbin_watson

plt.rcParams['figure.figsize'] = (9, 4)
plt.rcParams['axes.spines.top'] = False
plt.rcParams['axes.spines.right'] = False

np.random.seed(42)

import numpy as np import pandas as pd import matplotlib.pyplot as plt import scipy.stats as stats import statsmodels.formula.api as smf import statsmodels.api as sm from statsmodels.stats.outliers_influence import variance_inflation_factor from statsmodels.stats.diagnostic import het_breuschpagan, het_white from statsmodels.stats.stattools import durbin_watson plt.rcParams['figure.figsize'] = (9, 4) plt.rcParams['axes.spines.top'] = False plt.rcParams['axes.spines.right'] = False np.random.seed(42)

---

1. El Modelo de Regresión Múltiple¶

El modelo poblacional con $k$ predictores es:

$$Y_i = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \cdots + \beta_k X_{ki} + \varepsilon_i$$

donde $\varepsilon_i \sim \mathcal{N}(0, \sigma^2)$ i.i.d.

Interpretación de $\beta_j$: cambio esperado en $Y$ por unidad de aumento en $X_j$, manteniendo constantes los demás predictores. Esta es la diferencia clave respecto al modelo simple.

Supuestos (mismos que el modelo simple, ahora más críticos):

Supuesto	Descripción	Cómo verificar
Linealidad	Relación lineal entre cada $X_j$ e $Y$	Residuos vs $\hat{y}$
Independencia	Errores no correlacionados	Durbin-Watson
Homocedasticidad	Varianza constante de los errores	Breusch-Pagan
Normalidad	Errores distribuidos Normal	Q-Q plot, Jarque-Bera
No multicolinealidad	Predictores no altamente correlacionados	VIF

In [ ]:

# Dataset: precio de viviendas según características
# Variables: superficie (m²), habitaciones, distancia al centro (km), antigüedad (años)
np.random.seed(42)
n = 200

superficie   = np.random.uniform(40, 200, n)
habitaciones = np.random.randint(1, 6, n).astype(float)
distancia    = np.random.uniform(1, 30, n)
antiguedad   = np.random.uniform(0, 40, n)

precio = (1_500_000
          + 35_000 * superficie
          + 800_000 * habitaciones
          - 120_000 * distancia
          - 50_000  * antiguedad
          + np.random.normal(0, 2_000_000, n))

df = pd.DataFrame({
    'precio':       precio,
    'superficie':   superficie,
    'habitaciones': habitaciones,
    'distancia':    distancia,
    'antiguedad':   antiguedad
})

print(df.describe().round(0))

# Dataset: precio de viviendas según características # Variables: superficie (m²), habitaciones, distancia al centro (km), antigüedad (años) np.random.seed(42) n = 200 superficie = np.random.uniform(40, 200, n) habitaciones = np.random.randint(1, 6, n).astype(float) distancia = np.random.uniform(1, 30, n) antiguedad = np.random.uniform(0, 40, n) precio = (1_500_000 + 35_000 * superficie + 800_000 * habitaciones - 120_000 * distancia - 50_000 * antiguedad + np.random.normal(0, 2_000_000, n)) df = pd.DataFrame({ 'precio': precio, 'superficie': superficie, 'habitaciones': habitaciones, 'distancia': distancia, 'antiguedad': antiguedad }) print(df.describe().round(0))

---

2. Ajuste y Lectura del summary()¶

In [ ]:

modelo = smf.ols(
    formula='precio ~ superficie + habitaciones + distancia + antiguedad',
    data=df
).fit()

print(modelo.summary())

modelo = smf.ols( formula='precio ~ superficie + habitaciones + distancia + antiguedad', data=df ).fit() print(modelo.summary())

In [ ]:

# Acceso programático a los resultados
print("=== Coeficientes estimados ===")
for var, coef in modelo.params.items():
    pval = modelo.pvalues[var]
    sig  = '***' if pval < 0.001 else ('**' if pval < 0.01 else ('*' if pval < 0.05 else ''))
    print(f"  {var:<15}: {coef:>14,.0f}   p={pval:.4f} {sig}")

print(f"\nR²           : {modelo.rsquared:.4f}")
print(f"R² ajustado  : {modelo.rsquared_adj:.4f}")
print(f"F-statistic  : {modelo.fvalue:.2f}")
print(f"Prob(F)      : {modelo.f_pvalue:.6f}")
print(f"AIC          : {modelo.aic:.2f}")
print(f"BIC          : {modelo.bic:.2f}")

# Acceso programático a los resultados print("=== Coeficientes estimados ===") for var, coef in modelo.params.items(): pval = modelo.pvalues[var] sig = '***' if pval < 0.001 else ('**' if pval < 0.01 else ('*' if pval < 0.05 else '')) print(f" {var:<15}: {coef:>14,.0f} p={pval:.4f} {sig}") print(f"\nR² : {modelo.rsquared:.4f}") print(f"R² ajustado : {modelo.rsquared_adj:.4f}") print(f"F-statistic : {modelo.fvalue:.2f}") print(f"Prob(F) : {modelo.f_pvalue:.6f}") print(f"AIC : {modelo.aic:.2f}") print(f"BIC : {modelo.bic:.2f}")

---

3. Prueba F Global¶

La prueba F evalúa si el modelo en su conjunto tiene poder predictivo:

$$H_0: \beta_1 = \beta_2 = \cdots = \beta_k = 0 \qquad H_1: \text{al menos un } \beta_j \neq 0$$

El estadístico compara la varianza explicada por el modelo contra la varianza no explicada:

$$F = \frac{MSR}{MSE} = \frac{SSR/k}{SSE/(n-k-1)} \sim F_{k,\ n-k-1}$$

donde:

• $SSR$ = suma de cuadrados de la regresión (varianza explicada)
• $SSE$ = suma de cuadrados del error (varianza no explicada)
• $k$ = número de predictores

⚠️ Diferencia con la prueba t: la prueba F dice si el modelo en conjunto sirve. La prueba t dice si cada variable individualmente aporta algo, dado que las demás están en el modelo. Un modelo puede tener F significativa pero con algunas $t$ no significativas (variables redundantes).

In [ ]:

# Descomposición de la varianza (ANOVA de la regresión)
anova_tabla = sm.stats.anova_lm(modelo, typ=1)
print("=== Tabla ANOVA de la regresión ===")
print(anova_tabla.round(2))

# Descomposición de la varianza (ANOVA de la regresión) anova_tabla = sm.stats.anova_lm(modelo, typ=1) print("=== Tabla ANOVA de la regresión ===") print(anova_tabla.round(2))

In [ ]:

# Visualización de la prueba F
k  = modelo.df_model
df_resid = modelo.df_resid
F_obs = modelo.fvalue
p_F   = modelo.f_pvalue
F_crit = stats.f.ppf(0.95, dfn=k, dfd=df_resid)

x = np.linspace(0, F_obs * 1.4, 500)
fig, ax = plt.subplots()
ax.plot(x, stats.f.pdf(x, dfn=k, dfd=df_resid), 'steelblue', lw=2)

x_fill = x[x >= F_crit]
ax.fill_between(x_fill, stats.f.pdf(x_fill, dfn=k, dfd=df_resid),
                color='salmon', alpha=0.6, label='Región de rechazo (α=0.05)')
ax.axvline(F_obs,  color='darkred', lw=2.5, label=f'F observado = {F_obs:.2f}  (p = {p_F:.2e})')
ax.axvline(F_crit, color='black',   lw=1.5, linestyle='--', label=f'F crítico = {F_crit:.2f}')

ax.set_title(f'Prueba F global  —  F({int(k)}, {int(df_resid)})')
ax.set_xlabel('F')
ax.set_ylabel('Densidad')
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()

print(f"\nDecisión: {'Rechazar H₀ → el modelo es globalmente significativo' if p_F < 0.05 else 'No rechazar H₀'}")

# Visualización de la prueba F k = modelo.df_model df_resid = modelo.df_resid F_obs = modelo.fvalue p_F = modelo.f_pvalue F_crit = stats.f.ppf(0.95, dfn=k, dfd=df_resid) x = np.linspace(0, F_obs * 1.4, 500) fig, ax = plt.subplots() ax.plot(x, stats.f.pdf(x, dfn=k, dfd=df_resid), 'steelblue', lw=2) x_fill = x[x >= F_crit] ax.fill_between(x_fill, stats.f.pdf(x_fill, dfn=k, dfd=df_resid), color='salmon', alpha=0.6, label='Región de rechazo (α=0.05)') ax.axvline(F_obs, color='darkred', lw=2.5, label=f'F observado = {F_obs:.2f} (p = {p_F:.2e})') ax.axvline(F_crit, color='black', lw=1.5, linestyle='--', label=f'F crítico = {F_crit:.2f}') ax.set_title(f'Prueba F global — F({int(k)}, {int(df_resid)})') ax.set_xlabel('F') ax.set_ylabel('Densidad') ax.legend(fontsize=9) plt.tight_layout() plt.show() print(f"\nDecisión: {'Rechazar H₀ → el modelo es globalmente significativo' if p_F < 0.05 else 'No rechazar H₀'}")

---

4. Pruebas t por Coeficiente¶

Cada coeficiente tiene su propia prueba t, condicional a los demás predictores en el modelo:

$$H_0: \beta_j = 0 \quad \text{(dado que los demás } X \text{ están en el modelo)}$$

$$t_j = \frac{\hat{\beta}_j}{SE(\hat{\beta}_j)} \sim t_{n-k-1}$$

Un coeficiente no significativo puede indicar:

1. La variable genuinamente no aporta información adicional
2. La variable está correlacionada con otra (multicolinealidad)
3. El tamaño de muestra es insuficiente para detectar el efecto

In [ ]:

# Tabla completa de inferencia por coeficiente
ic = modelo.conf_int(alpha=0.05)
ic.columns = ['IC 2.5%', 'IC 97.5%']

tabla_inf = pd.DataFrame({
    'Coef':    modelo.params,
    'SE':      modelo.bse,
    't':       modelo.tvalues,
    'p-value': modelo.pvalues,
    'IC 2.5%': ic['IC 2.5%'],
    'IC 97.5%':ic['IC 97.5%'],
    'Sig.':    modelo.pvalues.apply(
        lambda p: '***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else '')))
})
print(tabla_inf.round(4))

# Tabla completa de inferencia por coeficiente ic = modelo.conf_int(alpha=0.05) ic.columns = ['IC 2.5%', 'IC 97.5%'] tabla_inf = pd.DataFrame({ 'Coef': modelo.params, 'SE': modelo.bse, 't': modelo.tvalues, 'p-value': modelo.pvalues, 'IC 2.5%': ic['IC 2.5%'], 'IC 97.5%':ic['IC 97.5%'], 'Sig.': modelo.pvalues.apply( lambda p: '***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else ''))) }) print(tabla_inf.round(4))

In [ ]:

# Visualización: coeficientes con IC (forest plot)
vars_plot = [v for v in modelo.params.index if v != 'Intercept']
coefs  = modelo.params[vars_plot]
ic_low = ic.loc[vars_plot, 'IC 2.5%']
ic_up  = ic.loc[vars_plot, 'IC 97.5%']

# Estandarizamos visualmente dividiendo por el rango del IC
fig, ax = plt.subplots(figsize=(9, 4))
colores_coef = ['steelblue' if modelo.pvalues[v] < 0.05 else 'gray' for v in vars_plot]

for i, (var, color) in enumerate(zip(vars_plot, colores_coef)):
    ax.plot([ic_low[var], ic_up[var]], [i, i], color=color, lw=4, alpha=0.7)
    ax.plot(coefs[var], i, 'o', color=color, ms=8, zorder=5)

ax.axvline(0, color='red', lw=1.5, linestyle='--', label='β = 0')
ax.set_yticks(range(len(vars_plot)))
ax.set_yticklabels(vars_plot)
ax.set_title('Coeficientes con IC 95%  (azul = significativo, gris = no significativo)')
ax.set_xlabel('Valor del coeficiente (CLP)')
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()

# Visualización: coeficientes con IC (forest plot) vars_plot = [v for v in modelo.params.index if v != 'Intercept'] coefs = modelo.params[vars_plot] ic_low = ic.loc[vars_plot, 'IC 2.5%'] ic_up = ic.loc[vars_plot, 'IC 97.5%'] # Estandarizamos visualmente dividiendo por el rango del IC fig, ax = plt.subplots(figsize=(9, 4)) colores_coef = ['steelblue' if modelo.pvalues[v] < 0.05 else 'gray' for v in vars_plot] for i, (var, color) in enumerate(zip(vars_plot, colores_coef)): ax.plot([ic_low[var], ic_up[var]], [i, i], color=color, lw=4, alpha=0.7) ax.plot(coefs[var], i, 'o', color=color, ms=8, zorder=5) ax.axvline(0, color='red', lw=1.5, linestyle='--', label='β = 0') ax.set_yticks(range(len(vars_plot))) ax.set_yticklabels(vars_plot) ax.set_title('Coeficientes con IC 95% (azul = significativo, gris = no significativo)') ax.set_xlabel('Valor del coeficiente (CLP)') ax.legend(fontsize=9) plt.tight_layout() plt.show()

---

5. Multicolinealidad y VIF¶

La multicolinealidad ocurre cuando dos o más predictores están altamente correlacionados entre sí. Esto no afecta las predicciones del modelo, pero infla los errores estándar de los coeficientes, haciendo que pruebas t sean no significativas aunque el efecto sea real.

El Factor de Inflación de la Varianza (VIF) cuantifica cuánto se infla la varianza de $\hat{\beta}_j$ debido a la correlación con los demás predictores:

$$VIF_j = \frac{1}{1 - R^2_j}$$

donde $R^2_j$ es el $R^2$ de la regresión de $X_j$ sobre todos los demás predictores.

Reglas prácticas:

VIF	Interpretación
$< 5$	Sin problemas
$5 - 10$	Multicolinealidad moderada, monitorear
$> 10$	Multicolinealidad severa, tomar acción

In [ ]:

# Cálculo del VIF para cada predictor
X = df[['superficie', 'habitaciones', 'distancia', 'antiguedad']]
X_const = sm.add_constant(X)

vif_data = pd.DataFrame({
    'Variable': X.columns,
    'VIF': [variance_inflation_factor(X_const.values, i+1)
            for i in range(X.shape[1])]
})
vif_data['Interpretación'] = vif_data['VIF'].apply(
    lambda v: 'Sin problema' if v < 5 else ('Moderada' if v < 10 else 'Severa')
)

print("=== Factor de Inflación de la Varianza (VIF) ===")
print(vif_data.round(3).to_string(index=False))

# Cálculo del VIF para cada predictor X = df[['superficie', 'habitaciones', 'distancia', 'antiguedad']] X_const = sm.add_constant(X) vif_data = pd.DataFrame({ 'Variable': X.columns, 'VIF': [variance_inflation_factor(X_const.values, i+1) for i in range(X.shape[1])] }) vif_data['Interpretación'] = vif_data['VIF'].apply( lambda v: 'Sin problema' if v < 5 else ('Moderada' if v < 10 else 'Severa') ) print("=== Factor de Inflación de la Varianza (VIF) ===") print(vif_data.round(3).to_string(index=False))

In [ ]:

# Demostración del problema: modelo con multicolinealidad artificial
np.random.seed(7)
# Creamos superficie_2 altamente correlacionada con superficie
df['superficie_2'] = df['superficie'] * 0.95 + np.random.normal(0, 3, n)

modelo_mc = smf.ols(
    'precio ~ superficie + superficie_2 + habitaciones + distancia + antiguedad',
    data=df
).fit()

X_mc = df[['superficie', 'superficie_2', 'habitaciones', 'distancia', 'antiguedad']]
X_mc_const = sm.add_constant(X_mc)

vif_mc = pd.DataFrame({
    'Variable': X_mc.columns,
    'VIF': [variance_inflation_factor(X_mc_const.values, i+1)
            for i in range(X_mc.shape[1])]
})

print("=== VIF con multicolinealidad severa ===")
print(vif_mc.round(1).to_string(index=False))
print()
print("Coeficientes y p-values del modelo con multicolinealidad:")
print(pd.DataFrame({
    'coef':    modelo_mc.params[['superficie', 'superficie_2']],
    'p-value': modelo_mc.pvalues[['superficie', 'superficie_2']]
}).round(4))

# Eliminar columna auxiliar
df.drop(columns=['superficie_2'], inplace=True)

# Demostración del problema: modelo con multicolinealidad artificial np.random.seed(7) # Creamos superficie_2 altamente correlacionada con superficie df['superficie_2'] = df['superficie'] * 0.95 + np.random.normal(0, 3, n) modelo_mc = smf.ols( 'precio ~ superficie + superficie_2 + habitaciones + distancia + antiguedad', data=df ).fit() X_mc = df[['superficie', 'superficie_2', 'habitaciones', 'distancia', 'antiguedad']] X_mc_const = sm.add_constant(X_mc) vif_mc = pd.DataFrame({ 'Variable': X_mc.columns, 'VIF': [variance_inflation_factor(X_mc_const.values, i+1) for i in range(X_mc.shape[1])] }) print("=== VIF con multicolinealidad severa ===") print(vif_mc.round(1).to_string(index=False)) print() print("Coeficientes y p-values del modelo con multicolinealidad:") print(pd.DataFrame({ 'coef': modelo_mc.params[['superficie', 'superficie_2']], 'p-value': modelo_mc.pvalues[['superficie', 'superficie_2']] }).round(4)) # Eliminar columna auxiliar df.drop(columns=['superficie_2'], inplace=True)

In [ ]:

# Visualización: correlación entre predictores (heatmap)
corr = df[['superficie', 'habitaciones', 'distancia', 'antiguedad']].corr()

fig, ax = plt.subplots(figsize=(6, 5))
im = ax.imshow(corr, cmap='RdBu_r', vmin=-1, vmax=1)
plt.colorbar(im, ax=ax)

vars_names = corr.columns.tolist()
ax.set_xticks(range(len(vars_names)))
ax.set_yticks(range(len(vars_names)))
ax.set_xticklabels(vars_names, rotation=30, ha='right')
ax.set_yticklabels(vars_names)

for i in range(len(vars_names)):
    for j in range(len(vars_names)):
        ax.text(j, i, f'{corr.iloc[i, j]:.2f}', ha='center', va='center',
                fontsize=11, color='white' if abs(corr.iloc[i, j]) > 0.5 else 'black')

ax.set_title('Correlación entre predictores')
plt.tight_layout()
plt.show()

# Visualización: correlación entre predictores (heatmap) corr = df[['superficie', 'habitaciones', 'distancia', 'antiguedad']].corr() fig, ax = plt.subplots(figsize=(6, 5)) im = ax.imshow(corr, cmap='RdBu_r', vmin=-1, vmax=1) plt.colorbar(im, ax=ax) vars_names = corr.columns.tolist() ax.set_xticks(range(len(vars_names))) ax.set_yticks(range(len(vars_names))) ax.set_xticklabels(vars_names, rotation=30, ha='right') ax.set_yticklabels(vars_names) for i in range(len(vars_names)): for j in range(len(vars_names)): ax.text(j, i, f'{corr.iloc[i, j]:.2f}', ha='center', va='center', fontsize=11, color='white' if abs(corr.iloc[i, j]) > 0.5 else 'black') ax.set_title('Correlación entre predictores') plt.tight_layout() plt.show()

---

6. Supuestos del Modelo¶

6.1 Homocedasticidad — Prueba de Breusch-Pagan¶

La heterocedasticidad ocurre cuando la varianza de los errores no es constante a lo largo de los valores ajustados. Esto no sesga los coeficientes pero invalida los errores estándar y, por tanto, toda la inferencia.

Prueba de Breusch-Pagan: $$H_0: \text{homocedasticidad (varianza constante)}$$ $$H_1: \text{heterocedasticidad (varianza depende de } X \text{)}$$

In [ ]:

residuos = modelo.resid
y_hat    = modelo.fittedvalues
X_exog   = modelo.model.exog  # matriz de diseño incluyendo intercepto

# Prueba de Breusch-Pagan
bp_lm, bp_pvalue, bp_fvalue, bp_fp = het_breuschpagan(residuos, X_exog)

print("=== Prueba de Breusch-Pagan ===")
print(f"Estadístico LM:  {bp_lm:.4f}")
print(f"p-value:         {bp_pvalue:.4f}")
if bp_pvalue < 0.05:
    print("→ Rechazar H₀: hay evidencia de heterocedasticidad")
else:
    print("→ No rechazar H₀: sin evidencia de heterocedasticidad")

residuos = modelo.resid y_hat = modelo.fittedvalues X_exog = modelo.model.exog # matriz de diseño incluyendo intercepto # Prueba de Breusch-Pagan bp_lm, bp_pvalue, bp_fvalue, bp_fp = het_breuschpagan(residuos, X_exog) print("=== Prueba de Breusch-Pagan ===") print(f"Estadístico LM: {bp_lm:.4f}") print(f"p-value: {bp_pvalue:.4f}") if bp_pvalue < 0.05: print("→ Rechazar H₀: hay evidencia de heterocedasticidad") else: print("→ No rechazar H₀: sin evidencia de heterocedasticidad")

In [ ]:

# Visualización: residuos vs valores ajustados
fig, axes = plt.subplots(1, 2, figsize=(13, 4))

# Residuos vs ajustados
ax = axes[0]
ax.scatter(y_hat/1e6, residuos/1e6, alpha=0.5, color='steelblue', edgecolor='white', s=50)
ax.axhline(0, color='red', lw=1.5, linestyle='--')
ax.set_title('Residuos vs Valores ajustados')
ax.set_xlabel('Valores ajustados (millones CLP)')
ax.set_ylabel('Residuos (millones CLP)')

# Escala-Ubicación (√|residuos estandarizados| vs ajustados)
ax = axes[1]
resid_std = residuos / residuos.std()
ax.scatter(y_hat/1e6, np.sqrt(np.abs(resid_std)),
           alpha=0.5, color='seagreen', edgecolor='white', s=50)
ax.set_title('Scale-Location (homocedasticidad)')
ax.set_xlabel('Valores ajustados (millones CLP)')
ax.set_ylabel('√|Residuos estandarizados|')

plt.suptitle('Diagnóstico de homocedasticidad', y=1.02)
plt.tight_layout()
plt.show()

# Visualización: residuos vs valores ajustados fig, axes = plt.subplots(1, 2, figsize=(13, 4)) # Residuos vs ajustados ax = axes[0] ax.scatter(y_hat/1e6, residuos/1e6, alpha=0.5, color='steelblue', edgecolor='white', s=50) ax.axhline(0, color='red', lw=1.5, linestyle='--') ax.set_title('Residuos vs Valores ajustados') ax.set_xlabel('Valores ajustados (millones CLP)') ax.set_ylabel('Residuos (millones CLP)') # Escala-Ubicación (√|residuos estandarizados| vs ajustados) ax = axes[1] resid_std = residuos / residuos.std() ax.scatter(y_hat/1e6, np.sqrt(np.abs(resid_std)), alpha=0.5, color='seagreen', edgecolor='white', s=50) ax.set_title('Scale-Location (homocedasticidad)') ax.set_xlabel('Valores ajustados (millones CLP)') ax.set_ylabel('√|Residuos estandarizados|') plt.suptitle('Diagnóstico de homocedasticidad', y=1.02) plt.tight_layout() plt.show()

6.2 Autocorrelación — Durbin-Watson¶

La autocorrelación ocurre cuando los errores están correlacionados entre sí, lo que viola el supuesto de independencia. Es especialmente relevante en datos temporales o espaciales.

El estadístico de Durbin-Watson varía entre 0 y 4:

Valor DW	Interpretación
$\approx 2$	Sin autocorrelación
$< 1.5$	Autocorrelación positiva
$> 2.5$	Autocorrelación negativa

In [ ]:

dw = durbin_watson(residuos)
print("=== Estadístico de Durbin-Watson ===")
print(f"DW = {dw:.4f}")
if dw < 1.5:
    print("→ Posible autocorrelación positiva")
elif dw > 2.5:
    print("→ Posible autocorrelación negativa")
else:
    print("→ Sin evidencia de autocorrelación (valor cercano a 2)")

dw = durbin_watson(residuos) print("=== Estadístico de Durbin-Watson ===") print(f"DW = {dw:.4f}") if dw < 1.5: print("→ Posible autocorrelación positiva") elif dw > 2.5: print("→ Posible autocorrelación negativa") else: print("→ Sin evidencia de autocorrelación (valor cercano a 2)")

In [ ]:

# Visualización: residuos en orden de observación
fig, axes = plt.subplots(1, 2, figsize=(13, 4))

# Residuos vs índice
ax = axes[0]
ax.plot(residuos.values/1e6, color='steelblue', lw=1, alpha=0.8)
ax.axhline(0, color='red', lw=1.5, linestyle='--')
ax.set_title(f'Residuos en orden de observación  (DW = {dw:.3f})')
ax.set_xlabel('Observación')
ax.set_ylabel('Residuo (millones CLP)')

# Autocorrelación: residuo(t) vs residuo(t-1)
ax = axes[1]
ax.scatter(residuos.values[:-1]/1e6, residuos.values[1:]/1e6,
           alpha=0.5, color='seagreen', edgecolor='white', s=40)
ax.set_title('Residuo(t) vs Residuo(t-1)')
ax.set_xlabel('Residuo(t-1) (millones CLP)')
ax.set_ylabel('Residuo(t) (millones CLP)')

plt.suptitle('Diagnóstico de autocorrelación', y=1.02)
plt.tight_layout()
plt.show()

# Visualización: residuos en orden de observación fig, axes = plt.subplots(1, 2, figsize=(13, 4)) # Residuos vs índice ax = axes[0] ax.plot(residuos.values/1e6, color='steelblue', lw=1, alpha=0.8) ax.axhline(0, color='red', lw=1.5, linestyle='--') ax.set_title(f'Residuos en orden de observación (DW = {dw:.3f})') ax.set_xlabel('Observación') ax.set_ylabel('Residuo (millones CLP)') # Autocorrelación: residuo(t) vs residuo(t-1) ax = axes[1] ax.scatter(residuos.values[:-1]/1e6, residuos.values[1:]/1e6, alpha=0.5, color='seagreen', edgecolor='white', s=40) ax.set_title('Residuo(t) vs Residuo(t-1)') ax.set_xlabel('Residuo(t-1) (millones CLP)') ax.set_ylabel('Residuo(t) (millones CLP)') plt.suptitle('Diagnóstico de autocorrelación', y=1.02) plt.tight_layout() plt.show()

6.3 Normalidad de Residuos¶

Los errores deben distribuirse normalmente para que los p-values e intervalos de confianza sean exactos. Con muestras grandes, el TCL hace que la inferencia sea robusta a desviaciones moderadas.

In [ ]:

# Pruebas formales de normalidad
jb_stat, jb_pvalue = stats.jarque_bera(residuos)
sw_stat, sw_pvalue = stats.shapiro(residuos[:50])  # Shapiro-Wilk (mejor para n pequeño)

print("=== Pruebas de normalidad de residuos ===")
print(f"Jarque-Bera:   stat={jb_stat:.4f}  p={jb_pvalue:.4f}  "
      f"→ {'Normal' if jb_pvalue > 0.05 else 'No normal'} (α=0.05)")
print(f"Shapiro-Wilk:  stat={sw_stat:.4f}  p={sw_pvalue:.4f}  "
      f"→ {'Normal' if sw_pvalue > 0.05 else 'No normal'} (α=0.05, n=50)")

# Pruebas formales de normalidad jb_stat, jb_pvalue = stats.jarque_bera(residuos) sw_stat, sw_pvalue = stats.shapiro(residuos[:50]) # Shapiro-Wilk (mejor para n pequeño) print("=== Pruebas de normalidad de residuos ===") print(f"Jarque-Bera: stat={jb_stat:.4f} p={jb_pvalue:.4f} " f"→ {'Normal' if jb_pvalue > 0.05 else 'No normal'} (α=0.05)") print(f"Shapiro-Wilk: stat={sw_stat:.4f} p={sw_pvalue:.4f} " f"→ {'Normal' if sw_pvalue > 0.05 else 'No normal'} (α=0.05, n=50)")

In [ ]:

# Q-Q plot e histograma
fig, axes = plt.subplots(1, 2, figsize=(12, 4))

ax = axes[0]
sm.qqplot(residuos, line='s', ax=ax, alpha=0.6, color='steelblue')
ax.set_title('Q-Q Plot de residuos')

ax = axes[1]
ax.hist(residuos/1e6, bins=25, color='steelblue', edgecolor='white',
        density=True, alpha=0.8)
x_range = np.linspace(residuos.min(), residuos.max(), 200)
ax.plot(x_range/1e6,
        stats.norm.pdf(x_range, residuos.mean(), residuos.std()) * 1e6,
        'red', lw=2, label='Normal teórica')
ax.set_title('Distribución de residuos')
ax.set_xlabel('Residuos (millones CLP)')
ax.set_ylabel('Densidad')
ax.legend()

plt.tight_layout()
plt.show()

# Q-Q plot e histograma fig, axes = plt.subplots(1, 2, figsize=(12, 4)) ax = axes[0] sm.qqplot(residuos, line='s', ax=ax, alpha=0.6, color='steelblue') ax.set_title('Q-Q Plot de residuos') ax = axes[1] ax.hist(residuos/1e6, bins=25, color='steelblue', edgecolor='white', density=True, alpha=0.8) x_range = np.linspace(residuos.min(), residuos.max(), 200) ax.plot(x_range/1e6, stats.norm.pdf(x_range, residuos.mean(), residuos.std()) * 1e6, 'red', lw=2, label='Normal teórica') ax.set_title('Distribución de residuos') ax.set_xlabel('Residuos (millones CLP)') ax.set_ylabel('Densidad') ax.legend() plt.tight_layout() plt.show()

---

7. Panel de Diagnóstico Completo¶

Un resumen de todos los diagnósticos en un solo bloque para uso rápido:

In [ ]:

def diagnostico_regresion(modelo, df_features):
    """
    Imprime un resumen de diagnósticos para un modelo OLS de statsmodels.
    """
    residuos = modelo.resid
    X_exog   = modelo.model.exog
    dw       = durbin_watson(residuos)
    bp_lm, bp_p, _, _ = het_breuschpagan(residuos, X_exog)
    jb_stat, jb_p = stats.jarque_bera(residuos)

    X_const = sm.add_constant(df_features)
    vifs = [variance_inflation_factor(X_const.values, i+1)
            for i in range(df_features.shape[1])]

    print("=" * 55)
    print("         DIAGNÓSTICO DEL MODELO")
    print("=" * 55)
    print(f"  R²           : {modelo.rsquared:.4f}")
    print(f"  R² ajustado  : {modelo.rsquared_adj:.4f}")
    print(f"  F-stat       : {modelo.fvalue:.2f}  (p={modelo.f_pvalue:.2e})")
    print()
    print("  [Homocedasticidad — Breusch-Pagan]")
    print(f"  p-value: {bp_p:.4f}  "
          f"{'✅ Sin evidencia de heterocedasticidad' if bp_p > 0.05 else '⚠️  Heterocedasticidad detectada'}")
    print()
    print("  [Autocorrelación — Durbin-Watson]")
    estado_dw = '✅ Sin autocorrelación' if 1.5 <= dw <= 2.5 else '⚠️  Posible autocorrelación'
    print(f"  DW: {dw:.4f}  {estado_dw}")
    print()
    print("  [Normalidad — Jarque-Bera]")
    print(f"  p-value: {jb_p:.4f}  "
          f"{'✅ Sin evidencia de no normalidad' if jb_p > 0.05 else '⚠️  No normalidad detectada'}")
    print()
    print("  [Multicolinealidad — VIF]")
    for var, vif in zip(df_features.columns, vifs):
        estado = '✅' if vif < 5 else ('⚠️ ' if vif < 10 else '❌')
        print(f"  {var:<15}: {vif:.2f}  {estado}")
    print("=" * 55)


X_feat = df[['superficie', 'habitaciones', 'distancia', 'antiguedad']]
diagnostico_regresion(modelo, X_feat)

def diagnostico_regresion(modelo, df_features): """ Imprime un resumen de diagnósticos para un modelo OLS de statsmodels. """ residuos = modelo.resid X_exog = modelo.model.exog dw = durbin_watson(residuos) bp_lm, bp_p, _, _ = het_breuschpagan(residuos, X_exog) jb_stat, jb_p = stats.jarque_bera(residuos) X_const = sm.add_constant(df_features) vifs = [variance_inflation_factor(X_const.values, i+1) for i in range(df_features.shape[1])] print("=" * 55) print(" DIAGNÓSTICO DEL MODELO") print("=" * 55) print(f" R² : {modelo.rsquared:.4f}") print(f" R² ajustado : {modelo.rsquared_adj:.4f}") print(f" F-stat : {modelo.fvalue:.2f} (p={modelo.f_pvalue:.2e})") print() print(" [Homocedasticidad — Breusch-Pagan]") print(f" p-value: {bp_p:.4f} " f"{'✅ Sin evidencia de heterocedasticidad' if bp_p > 0.05 else '⚠️ Heterocedasticidad detectada'}") print() print(" [Autocorrelación — Durbin-Watson]") estado_dw = '✅ Sin autocorrelación' if 1.5 <= dw <= 2.5 else '⚠️ Posible autocorrelación' print(f" DW: {dw:.4f} {estado_dw}") print() print(" [Normalidad — Jarque-Bera]") print(f" p-value: {jb_p:.4f} " f"{'✅ Sin evidencia de no normalidad' if jb_p > 0.05 else '⚠️ No normalidad detectada'}") print() print(" [Multicolinealidad — VIF]") for var, vif in zip(df_features.columns, vifs): estado = '✅' if vif < 5 else ('⚠️ ' if vif < 10 else '❌') print(f" {var:<15}: {vif:.2f} {estado}") print("=" * 55) X_feat = df[['superficie', 'habitaciones', 'distancia', 'antiguedad']] diagnostico_regresion(modelo, X_feat)

---

8. Ejemplo Integrador¶

Contexto: Un banco quiere modelar el monto de crédito aprobado en función de: ingreso mensual, antigüedad laboral, número de productos financieros que tiene el cliente, y su score crediticio.

Ajustamos el modelo completo, interpretamos los resultados, y verificamos los supuestos.

In [ ]:

np.random.seed(88)
n = 250

ingreso        = np.random.normal(1_200_000, 400_000, n).clip(400_000, 4_000_000)
antiguedad     = np.random.uniform(0, 30, n)
n_productos    = np.random.randint(1, 8, n).astype(float)
score          = np.random.normal(650, 80, n).clip(400, 850)

credito = (2_000_000
           + 4.5  * ingreso
           + 200_000 * antiguedad
           + 500_000 * n_productos
           + 30_000  * score
           + np.random.normal(0, 5_000_000, n))

df_banco = pd.DataFrame({
    'credito':     credito,
    'ingreso':     ingreso,
    'antiguedad':  antiguedad,
    'n_productos': n_productos,
    'score':       score
})

modelo_banco = smf.ols(
    'credito ~ ingreso + antiguedad + n_productos + score',
    data=df_banco
).fit()

print(modelo_banco.summary())

np.random.seed(88) n = 250 ingreso = np.random.normal(1_200_000, 400_000, n).clip(400_000, 4_000_000) antiguedad = np.random.uniform(0, 30, n) n_productos = np.random.randint(1, 8, n).astype(float) score = np.random.normal(650, 80, n).clip(400, 850) credito = (2_000_000 + 4.5 * ingreso + 200_000 * antiguedad + 500_000 * n_productos + 30_000 * score + np.random.normal(0, 5_000_000, n)) df_banco = pd.DataFrame({ 'credito': credito, 'ingreso': ingreso, 'antiguedad': antiguedad, 'n_productos': n_productos, 'score': score }) modelo_banco = smf.ols( 'credito ~ ingreso + antiguedad + n_productos + score', data=df_banco ).fit() print(modelo_banco.summary())

In [ ]:

# Diagnóstico completo
X_banco = df_banco[['ingreso', 'antiguedad', 'n_productos', 'score']]
diagnostico_regresion(modelo_banco, X_banco)

# Diagnóstico completo X_banco = df_banco[['ingreso', 'antiguedad', 'n_productos', 'score']] diagnostico_regresion(modelo_banco, X_banco)

In [ ]:

# Interpretación de negocio
print("=== Interpretación de los coeficientes ===")
vars_inter = ['ingreso', 'antiguedad', 'n_productos', 'score']
unidades   = ['CLP de ingreso', 'año de antigüedad',
              'producto adicional', 'punto de score']

for var, unidad in zip(vars_inter, unidades):
    coef  = modelo_banco.params[var]
    pval  = modelo_banco.pvalues[var]
    sig   = '(significativo)' if pval < 0.05 else '(no significativo)'
    print(f"  {var:<15}: {coef:>12,.0f} CLP por {unidad}  {sig}")

# Interpretación de negocio print("=== Interpretación de los coeficientes ===") vars_inter = ['ingreso', 'antiguedad', 'n_productos', 'score'] unidades = ['CLP de ingreso', 'año de antigüedad', 'producto adicional', 'punto de score'] for var, unidad in zip(vars_inter, unidades): coef = modelo_banco.params[var] pval = modelo_banco.pvalues[var] sig = '(significativo)' if pval < 0.05 else '(no significativo)' print(f" {var:<15}: {coef:>12,.0f} CLP por {unidad} {sig}")

---

Resumen¶

Concepto	¿Qué evalúa?	Herramienta en statsmodels
Prueba F global	Significancia del modelo completo	modelo.fvalue, modelo.f_pvalue
Prueba t por coeficiente	Aporte individual de cada predictor	modelo.tvalues, modelo.pvalues
IC para coeficientes	Rango plausible del efecto real	modelo.conf_int()
VIF	Multicolinealidad entre predictores	variance_inflation_factor
Breusch-Pagan	Homocedasticidad de los errores	het_breuschpagan
Durbin-Watson	Autocorrelación de los errores	durbin_watson
Jarque-Bera / Shapiro	Normalidad de los errores	stats.jarque_bera, modelo.summary()
Q-Q plot	Normalidad visual	sm.qqplot

Ideas clave:

• La prueba F y las pruebas t responden preguntas distintas: no confundirlas
• VIF > 10 es señal de alerta: considerar eliminar variables redundantes o usar regularización
• Heterocedasticidad no sesga los coeficientes pero invalida los errores estándar → usar errores robustos (cov_type='HC3')
• La normalidad de residuos es el supuesto menos crítico con $n$ grande (TCL)
• Siempre correr el diagnóstico completo antes de reportar resultados

---

Referencias¶

1. James, G., Witten, D., Hastie, T., Tibshirani, R. (2021). An Introduction to Statistical Learning. Springer. Cap. 3.
2. Greene, W. H. (2018). Econometric Analysis. Pearson. Cap. 4-5.
3. Statsmodels — Diagnósticos: https://www.statsmodels.org/stable/stats.html#residual-diagnostics-and-specification-tests
4. Statsmodels — VIF: https://www.statsmodels.org/stable/generated/statsmodels.stats.outliers_influence.variance_inflation_factor.html