Diagnóstico de regresión: leverage, outliers, observaciones influyentes y multicolinealidad
Autor/a
Alberto Dorantes
Fecha de publicación
30 de septiembre de 2026
Resumen
Solución detallada del Workshop 2 del Bloque 2. Aplicamos, con álgebra matricial, las técnicas de diagnóstico de la Lección 9 (Influential Points) de STAT 462 de PennState al modelo del Challenge 4 del Workshop 1: leverages, residuales estandarizados, residuales eliminados (studentized deleted residuals), DFFITS y distancia de Cook. Después proponemos dos correcciones (eliminar observaciones influyentes y winsorizar) y comparamos los resultados con el modelo original.
1 Preparación: el modelo del Challenge 4 del Workshop 1
Todas las medidas de diagnóstico de la lección salen de la matriz de diseño X y de la hat matrixH. Construimos Xa mano: una columna de unos (intercepto) y una columna por cada variable del modelo.
La columna de unos al inicio de la matriz X es necesaria para que se estime el intercepto (beta0).
Genero las variables de X como arreglos NumPy (d[“columna”].values) a partir de las columnas del dataset. Esto para facilitar las operaciones con álgebra matricial. y generamos X como una matriz a partir de estos arreglos.
med = (d["tamano"] =="Mediana").astype(float).valuesgde = (d["tamano"] =="Grande").astype(float).valuesupa = d["upa"].valuesda = d["deuda_activos"].valuesX = np.column_stack([np.ones(len(d)), med, gde, upa, med*upa, gde*upa, da, da**2])y = d["roa"].valuesnombres = ["Intercepto", "Mediana", "Grande", "upa", "Mediana:upa","Grande:upa", "deuda_activos", "deuda_activos^2"]n, p = X.shape # p = número de parámetros (k + 1)k = p -1# k = número de variables X (predictores)print(f"n = {n}, p = k + 1 = {p}, k = {k}")
n = 1069, p = k + 1 = 8, k = 7
Recordemos que los coeficientes beta OLS (mínimos cuadrados) de una regresión múltiple se estiman como la siguiente multiplicación de matrices:
b = (X'X)^{-1}X'y Como ya tengo la matriz X y la matriz y, puedo fácilmente estimar el vector de coeficientes de beta con operaciones de álgebra matricial:
# Primero calcula la inversa de X'X:XtX_inv = np.linalg.inv(X.T @ X)# Calculo el vector de coeficientes b: b = XtX_inv @ X.T @ y# X.T es la transpuesta de X # Genero un dataset con los valores de b y con los valores de coeficientes del modelo de regresión # estimado anteriormente con statsmodel de Python:comp = pd.DataFrame({"b (matricial)": b, "b (statsmodels)": m_orig.params.values}, index=nombres)print(comp.round(6))
b (matricial) b (statsmodels)
Intercepto -0.068406 -0.068406
Mediana 0.047039 0.047039
Grande 0.117978 0.117978
upa 0.842516 0.842516
Mediana:upa 0.248500 0.248500
Grande:upa 0.082133 0.082133
deuda_activos 0.020671 0.020671
deuda_activos^2 -0.114831 -0.114831
Vemos que las estimaciones de los coeficientes calculados con álgebra matricial son idénticos a los estimados con statsmodel.
3 Leverages: valores extremos en X (Lección 9.2)
3.1 Concepto y cálculo
Recordemos la ecuación de las predicciones de y en un modelo de regresión múltiple que tiene más de 1 variable x, representadas en una matriz X (que tiene 1’s en la primera columna):
Y=Xb+\epsilon
Esta representación matricial indica que el valor de Y_i es una combinación lineal de cada valor X_i y su coeficiente b_i respectivo:
Y = b_0 + b_1x_1+b_2x_2+...+b_nx_n+\epsilon
Entonces, el valor esperado de Y, llamado la predicción \hat{y} es:
\hat{Y} = Xb
Esto mismo expresado en los valores puntuales de cada observación de Y:
\hat{y_i} = b_0 + b_1x_1+b_2x_2+...+b_nx_n
Como el vector de coeficientes b se estima como:
b = (X'X)^{-1}X'y
Entonces, puedo representar el valor esperado (predecido) de Y como:
\hat{Y} = X(X'X)^{-1}X'y
Si defino una matriz llamada hat matrix H como:
H = X(X'X)^{-1}X'
Entonces, puedo representar \hat{Y} como:
\hat{Y} = Hy
Se llama así porque la hat matriz H “le pone el gorro” a y: \hat{y} = Hy. Es decir, cada valor predecido de y es una combinación lineal de todos los valores observados de y:
\hat{y}_i = h_{i1}y_1 + h_{i2}y_2 + \dots + h_{ii}y_i + \dots + h_{in}y_n El elemento de la diagonal h_{ii} es el leverage de la observación i: mide cuánto pesa y_i en su propio valor predecido. El leverage sólo depende de las X (no de y), y mide qué tan lejos está la combinación de valores X de la observación respecto del centro de los datos.
Propiedades:
0 \le h_{ii} \le 1
\sum_{i=1}^{n} h_{ii} = p = k+1, así que el leverage promedio es \bar{h} = (k+1)/n
Calculamos la hat matrix utilizando álgebra matricial:
H = X @ XtX_inv @ X.T # Hat matrix: matriz n x nh = np.diag(H).copy() # valores leverages que están en la diagonal de Hprint("Dimensión de H:", H.shape)
Dimensión de H: (1069, 1069)
print(f"Suma de leverages = {h.sum():.4f} (debe ser p = {p})")
print("¿H y-hat = valores ajustados de statsmodels?", np.allclose(H @ y, m_orig.fittedvalues.values))
¿H y-hat = valores ajustados de statsmodels? True
Vemos que el cálculo de H es correcto.
3.2 Regla para identificar leverages altos
La lección propone marcar como valores X inusuales las observaciones con:
h_{ii} > 3\left(\frac{k+1}{n}\right) Que quiere decir que una observación x inusual (valor extremo) es aquella que tiene leverage 3 veces más grande que el promedio.
(y una regla menos estricta es que sea 2 veces más grande que el promedio).
Ordeno las observaciones de acuerdo a su leverage (h) de mayor a menor:
lim_h3 =3* p / nlim_h2 =2* p / nd["h"] = hprint(f"Límite 3p/n = {lim_h3:.4f}: {(h > lim_h3).sum()} observaciones con leverage alto")
Límite 3p/n = 0.0225: 36 observaciones con leverage alto
Leverage de cada observación. Las líneas marcan 2p/n y 3p/n
Interpretación: 36 empresas (3.4%) superan el límite de 3p/n = 0.0225. El caso más extremo es BBIO, con un leverage de 0.79: esta sola empresa tiene casi 80% del peso en su propio valor predicho. La razón es su nivel de deuda: su deuda equivale a 3.18 veces sus activos, y en el modelo aparece además al cuadrado (10.1), muy lejos del resto de las empresas. Otras empresas con leverage alto (RUN, PBF, NVAX, BYND) combinan valores de upa y deuda poco comunes para su grupo de tamaño.
La siguiente gráfica muestra el “espacio X” de upa y deuda: los puntos con leverage alto están en la periferia de la nube de datos.
Espacio de las variables X: las observaciones con leverage alto están lejos del centro
Recordemos que un leverage alto sólo indica potencial de influencia. Si la observación sigue el mismo patrón que el resto, su residual será pequeño y no cambiará mucho la regresión.
4 Residuales y residuales estandarizados (Lección 9.3)
4.1 Residuales ordinarios
e = y - \hat{y} = (I - H)\,y
El problema de los residuales ordinarios es que dependen de las unidades de Y y que su varianza no es constante: Var(e_i) = \sigma^2(1 - h_{ii}). Observaciones con leverage alto “jalan” la recta hacia ellas, así que tienden a tener residuales más pequeños.
Como r_i es aproximadamente una variable estándar, una regla práctica es marcar como outlier cualquier observación con |r_i| > 3.
e = y - H @ ySSE = e @ eMSE = SSE / (n - p)r = e / np.sqrt(MSE * (1- h))d["e"], d["r"] = e, rprint(f"SSE = {SSE:.4f}, MSE = {MSE:.6f}, raíz(MSE) = {np.sqrt(MSE):.4f}")
SSE = 42.2040, MSE = 0.039778, raíz(MSE) = 0.1994
print(f"Observaciones con |r| > 3: {(np.abs(r) >3).sum()}")
Observaciones con |r| > 3: 18
5 Residuales eliminados (Lección 9.4)
5.1 Deleted residuals (errores de predicción PRESS)
Un outlier “jala” la regresión hacia él y esconde su propio residual. La solución es eliminar la observación i, estimar el modelo con las n-1 observaciones restantes y predecir y_i con ese modelo:
d_i = y_i - \hat{y}_{(i)}
Afortunadamente no hay que correr n regresiones: se puede demostrar que
Comparación de residuales estandarizados vs studentized deleted residuals
Interpretación: 18 empresas tienen |t_i| > 3. El outlier más extremo es AWH (t = -13.96): una empresa pequeña con ROA de -301% (pierde 3 dólares por cada dólar de activos), mucho peor de lo que el modelo predice dado su upa. Otros outliers como IRWD, BFRG y LUCD también son empresas con pérdidas enormes respecto a sus activos. En la gráfica de r vs t se ve que los residuales eliminados son más extremos que los estandarizados para los outliers, porque al quitar la observación el modelo ya no se “acomoda” a ella.
6 Observaciones influyentes (Lección 9.5)
Una observación es influyente si al eliminarla cambian de manera importante los coeficientes o las predicciones. Esto suele pasar cuando una observación es outlier y además tiene leverage alto.
6.1 DFFITS
Mide cuántos errores estándar cambia el valor predicho de la observación i cuando la eliminamos:
La segunda forma muestra que DFFITS combina el tamaño del residual (t_i) con el leverage (h_{ii}). La lección marca como influyente una observación con:
print(f"Observaciones con |DFFITS| > límite: {(np.abs(dffits) > lim_dffits).sum()}")
Observaciones con |DFFITS| > límite: 60
# Verificación con la definición (re-estimando sin la obs. i)print(f"\nVerificación para {d.loc[i, 'empresa']}: "f"fórmula = {dffits[i]:.5f} | definición = "f"{(H[i] @ y - yhat_i) / np.sqrt(MSE_sin * h[i]):.5f}")
Verificación para AWH: fórmula = -1.10924 | definición = -1.10924
6.2 Distancia de Cook
Mide cuánto cambian todos los valores predecidos cuando se elimina la observación i:
DFFITS (arriba) y distancia de Cook (abajo) para cada observación
Interpretación:
Ninguna observación supera D_i > 0.5, pero dos empresas sobresalen claramente del resto: NVAX (D = 0.31) e IRWD (D = 0.31), seguidas de AWH (D = 0.13). Según el criterio de “sobresale del resto”, son influyentes.
60 observaciones (5.6%) superan el límite de DFFITS.
NVAX e IRWD son influyentes porque combinan leverage alto y residual grande. AWH tiene un residual enorme pero leverage moderado. BBIO, en cambio, tiene el leverage más alto pero un residual pequeño: sigue el patrón del modelo y por eso no es influyente (D = 0.016).
6.3 Gráfica resumen: leverage vs residual
Esta gráfica junta los tres conceptos: el eje X es el leverage, el eje Y el residual eliminado, y el tamaño del círculo es la distancia de Cook. Las observaciones influyentes están en las esquinas (leverage alto y residual grande).
La mayoría de las observaciones influyentes son empresas pequeñas con pérdidas muy grandes respecto a sus activos o a su valor de mercado. Esto es consistente con la distribución del ROA y la upa: tienen colas izquierdas muy largas.
7 Multicolinealidad
Como el modelo tiene interacciones y un término cuadrático, conviene revisar la multicolinealidad. El Factor de Inflación de Varianza de cada variable es:
VIF_j = \frac{1}{1 - R^2_j}
donde R^2_j es el R² de la regresión de X_j contra las demás X. Con álgebra matricial, los VIF son la diagonal de la inversa de la matriz de correlación de las variables X:
Todos los VIF son menores a 3 (la regla usual de alerta es VIF > 5 o 10), así que no hay un problema de multicolinealidad. El VIF más alto es el de la deuda y su cuadrado (≈ 2.8), lo cual es natural porque DA y DA^2 están correlacionadas. Si fuera necesario, se puede reducir centrando la variable antes de elevarla al cuadrado:
Antes de eliminar datos hay que preguntarse por qué son extremos:
¿Es un error de datos? Si fuera así, se corrige o se elimina. Aquí no parece ser el caso: son empresas reales (por ejemplo, biotecnológicas en etapa de investigación como AWH, NVAX o IRWD) con pérdidas muy grandes, o empresas con deuda mayor que sus activos (patrimonio negativo, como BBIO).
¿Pertenecen a la población que queremos estudiar? Si queremos entender la productividad de empresas “típicas” en operación, las empresas con pérdidas extremas o con más deuda que activos pueden distorsionar las conclusiones.
Proponemos y comparamos dos correcciones:
Corrección A - Eliminar observaciones influyentes: quitamos las observaciones con |DFFITS| mayor al límite de la lección. Es la corrección más directa, pero se pierden datos (y podríamos introducir sesgo si eliminamos sistemáticamente a las empresas con pérdidas).
Corrección B - Winsorizar: reemplazamos los valores menores al percentil 1 por el valor del percentil 1, y los mayores al percentil 99 por el valor del percentil 99, en las variables continuas del modelo (roa, upa y deuda/activos). Así conservamos todas las empresas, pero limitamos cuánto pueden “jalar” la recta los valores extremos.
En ambos casos mantenemos la clasificación de tamaño original (es un atributo de cada empresa).
Observa que siempre habrá observaciones marcadas: los criterios son relativos al modelo y a los datos que se usan. Al quitar observaciones, otras empresas pasan a verse “extremas”. Por eso no se debe repetir la eliminación hasta que no quede ninguna marcada (acabaríamos eliminando una gran parte de la muestra). Lo importante es que:
Con la corrección B, el caso de leverage extremo desaparece: la deuda/activos de BBIO se recorta a 1.22 y el leverage máximo baja de 0.79 a 0.15. El |t| máximo baja de 13.96 a 5.78, porque el ROA de AWH se recorta a -1.24.
Ojo: la distancia de Cook máxima no baja con la winsorización (0.38, de NVAX). Además, BBIO pasa a ser influyente (D = 0.27): al recortar su deuda queda en una zona donde su ROA de -111% ya no es consistente con el modelo. Aun así, ninguna observación supera el criterio de D > 0.5 y ya no hay casos con leverage cercano a 1.
Con la corrección A, las observaciones más extremas ya no están en la muestra: el |t| máximo baja a 5.40, el leverage máximo a 0.18 y la distancia de Cook máxima a 0.08.
9 Comparación con el Challenge 4 del Workshop 1
def tabla(modelo):return (modelo.params.round(4).astype(str) +" ("+ modelo.pvalues.map(lambda v: f"{v:.3f}") +")")comp = pd.concat([tabla(m_orig), tabla(m_A), tabla(m_B)], axis=1)comp.columns = ["Original (W1-C4)", "A: Eliminar influyentes", "B: Winsorizar 1-99%"]comp.loc["R²"] = [f"{m.rsquared:.4f}"for m in [m_orig, m_A, m_B]]comp.loc["n"] = [int(m.nobs) for m in [m_orig, m_A, m_B]]print("Coeficiente (p-value)")
Coeficientes e intervalos de confianza al 95% en las tres versiones del modelo
9.1 Efecto sobre las pendientes de upa por tamaño
def pendientes(mm): bb = mm.paramsreturn {"Pequeña": bb["upa"],"Mediana": bb["upa"] + bb["C(tamano)[T.Mediana]:upa"],"Grande": bb["upa"] + bb["C(tamano)[T.Grande]:upa"]}print(pd.DataFrame({nom: pendientes(mm) for nom, (mm, _) in modelos.items()}).round(4))
Original A: Eliminar B: Winsorizar
Pequeña 0.8425 0.8465 0.8421
Mediana 1.0910 1.2884 1.0535
Grande 0.9246 0.9762 1.0815
9.2 Efecto sobre la relación cuadrática de la deuda
for nom, (mm, cc) in modelos.items(): b6 = mm.params["deuda_activos"]; b7 = mm.params["I(deuda_activos ** 2)"]print(f"{nom:14s}: b6 = {b6:7.4f}, b7 = {b7:7.4f}, máximo en DA = {-b6/(2*b7):.3f}")
Original : b6 = 0.0207, b7 = -0.1148, máximo en DA = 0.090
A: Eliminar : b6 = 0.0095, b7 = -0.0610, máximo en DA = 0.078
B: Winsorizar : b6 = 0.0961, b7 = -0.2029, máximo en DA = 0.237
grid = np.linspace(0, 1.3, 200)upa_med = d["upa"].median()fig, ax = plt.subplots(figsize=(8, 4.2))for nom, (mm, cc) in modelos.items(): nd = pd.DataFrame({"tamano": pd.Categorical(["Pequeña"]*len(grid), categories=d["tamano"].cat.categories),"upa": upa_med, "deuda_activos": grid}) ax.plot(grid, mm.predict(nd), color=cc, lw=2.5, label=nom)ax.set_xlabel("Deuda / Activos"); ax.set_ylabel("ROA esperado")ax.set_title("La curvatura de la deuda depende de las observaciones extremas")ax.legend(fontsize=8)plt.tight_layout(); plt.show()
Efecto de la deuda sobre el ROA (empresa pequeña con upa mediana) en las tres versiones
9.3 Resumen de la comparación
El efecto de upa es prácticamente idéntico en las tres versiones (0.842 - 0.846): es el resultado más estable.
La pendiente de upa es mayor para las medianas en las tres versiones, aunque la magnitud de la diferencia cambia (0.22 a 0.46).
El efecto no lineal de la deuda es frágil: con winsorización se mantiene (máximo en DA ≈ 0.24), pero al eliminar las observaciones influyentes ningún término de deuda es significativo.
El R² sube de 0.594 a 0.655 (winsorizado) y a 0.743 (eliminando influyentes). Esto no significa que el modelo sea “mejor” por sí mismo: al reducir la variabilidad extrema de Y es natural que el R² aumente.
Aun después de las correcciones, los residuales tienen colas pesadas (ver Kurtosis y Jarque-Bera en los resúmenes) y, como veremos, su varianza no es constante. Por eso el siguiente paso es usar errores estándar robustos y regresión robusta.
10 Errores estándar robustos y regresión robusta
Hay dos herramientas con nombres parecidos que atacan problemas distintos:
Errores estándar robustos
Regresión robusta
¿Qué problema atiende?
Heterocedasticidad (varianza del error no constante)
Outliers en Y que “jalan” los coeficientes
¿Cambia los coeficientes?
No, son los mismos de OLS
Sí
¿Qué cambia?
Errores estándar, valores t, p-values e intervalos
Coeficientes y errores estándar
¿Cómo?
Fórmula “sándwich” de la varianza
Mínimos cuadrados ponderados iterativos: las observaciones atípicas pesan menos
Usaremos ambas sobre el modelo winsorizado (corrección B), porque la winsorización ya controló los problemas de leverage (valores extremos en X), algo que, como veremos, la regresión robusta por sí sola no resuelve.
10.1 Errores estándar robustos (HC3)
10.1.1 ¿Para qué sirven?
La fórmula de OLS para la varianza de los coeficientes supone que todos los errores tienen la misma varianza\sigma^2 (homocedasticidad):
Var(b) = MSE\,(X'X)^{-1}
Si la varianza del error es diferente entre observaciones (heterocedasticidad), los coeficientes de OLS siguen siendo insesgados, pero los errores estándar están mal calculados y, por lo tanto, los valores t y los p-values no son confiables. Podríamos concluir que un coeficiente es significativo cuando no lo es.
10.1.2 ¿Hay heterocedasticidad en nuestro modelo?
Veamos la dispersión de los residuales por grupo de tamaño y la gráfica de residuales vs valores ajustados:
import statsmodels.api as smfrom statsmodels.stats.diagnostic import het_breuschpaganres_B = m_B.residprint("Desviación estándar de los residuales por tamaño (modelo winsorizado):")
Desviación estándar de los residuales por tamaño (modelo winsorizado):
tamano
Pequeña 0.2342
Mediana 0.1253
Grande 0.1078
dtype: float64
bp = het_breuschpagan(res_B, m_B.model.exog)print(f"\nPrueba de Breusch-Pagan: LM = {bp[0]:.2f}, p-value = {bp[1]:.2e}")
Prueba de Breusch-Pagan: LM = 198.31, p-value = 2.61e-39
fig, ax = plt.subplots(figsize=(8, 4.2))for g in ["Pequeña", "Mediana", "Grande"]: mk = (dB["tamano"] == g).values ax.scatter(m_B.fittedvalues[mk], res_B[mk], s=10, alpha=0.4, color=COL[g], marker=MRK[g], label=g)ax.axhline(0, color="black", lw=0.8)ax.set_xlabel("ROA ajustado"); ax.set_ylabel("Residual")ax.set_title("Heterocedasticidad: residuales más dispersos en pequeñas")ax.legend(fontsize=8)plt.tight_layout(); plt.show()
Residuales vs valores ajustados del modelo winsorizado: la dispersión no es constante
La desviación estándar de los residuales de las empresas pequeñas es más del doble que la de las grandes, y la prueba de Breusch-Pagan rechaza contundentemente la hipótesis nula de homocedasticidad. Los errores estándar de OLS no son confiables.
10.1.3 La fórmula sándwich
Los errores robustos estiman la varianza permitiendo que cada observación tenga su propia varianza, usando su residual al cuadrado como estimación:
Hay varias versiones según cómo se define la matriz diagonal \Omega:
HC0 (White): \Omega_{ii} = e_i^2
HC3: \Omega_{ii} = \dfrac{e_i^2}{(1-h_{ii})^2}
HC3 es la más recomendada en muestras con observaciones de leverage alto: divide entre (1-h_{ii})^2 porque, como vimos en la Lección 9.4, las observaciones con leverage alto “jalan” la recta y tienen residuales artificialmente pequeños (es el mismo ajuste que convierte un residual en deleted residual).
Calculamos HC3 con álgebra matricial y lo comparamos con statsmodels:
Errores estándar de OLS vs HC3 (modelo winsorizado)
Interpretación:
Los coeficientes son los mismos; sólo cambian los errores estándar.
El error estándar de Mediana × upa se duplica (0.065 → 0.130) y su p-value pasa de 0.001 a 0.096: con errores robustos, la diferencia en pendientes entre medianas y pequeñas sólo es significativa al 10%.
Grande × upa pasa de p = 0.065 a p = 0.185: no es significativa.
Los términos de deuda también pierden significancia: el cuadrático pasa de p = 0.001 a p = 0.057 y el lineal de 0.082 a 0.214.
upa y la dummy Grande siguen siendo muy significativas.
Es decir, parte de la significancia que veníamos reportando se debía a suponer varianza constante cuando no la hay. ¿Por qué crece tanto el error de las interacciones? Porque la pendiente de upa de las empresas pequeñas se estima con observaciones muy ruidosas (residuales con el doble de dispersión), y OLS lo ignoraba.
10.2 Regresión robusta (estimadores M)
10.2.1 ¿Qué es y para qué sirve?
OLS minimiza \sum e_i^2. Al elevar al cuadrado, un residual de 1.0 pesa 100 veces más que uno de 0.1: unos cuantos outliers pueden mover toda la recta. La regresión robusta minimiza una función \sum \rho(e_i/s) que crece más lento que el cuadrado para residuales grandes (s es una medida robusta de la escala de los residuales).
En la práctica, esto equivale a una regresión de mínimos cuadrados ponderados donde cada observación recibe un peso w_i según el tamaño de su residual estandarizado u_i = e_i/s:
Huber: w_i = 1 si |u_i| \le 1.345; si no, w_i = 1.345/|u_i|. Los outliers pesan menos, pero nunca cero.
Tukey (bicuadrada): w_i = \left[1 - (u_i/4.685)^2\right]^2 si |u_i| \le 4.685; si no, w_i = 0. Los outliers extremos se ignoran por completo.
(Las constantes 1.345 y 4.685 hacen que el estimador sea 95% tan eficiente como OLS cuando los errores sí son normales.)
u = np.linspace(-7, 7, 400)w_hub = np.where(np.abs(u) <=1.345, 1, 1.345/ np.abs(u))w_tuk = np.where(np.abs(u) <=4.685, (1- (u /4.685)**2)**2, 0)fig, ax = plt.subplots(figsize=(8, 3.6))ax.plot(u, np.ones_like(u), color=GRIS, lw=2.5, label="OLS (todas pesan 1)")ax.plot(u, w_hub, color=NARANJA, lw=2.5, label="Huber")ax.plot(u, w_tuk, color=AZUL, lw=2.5, label="Tukey (bicuadrada)")ax.set_xlabel("Residual estandarizado u = e / s"); ax.set_ylabel("Peso w")ax.set_title("Regresión robusta: los residuales grandes pesan menos")ax.legend(fontsize=8)plt.tight_layout(); plt.show()
Peso que recibe una observación según su residual estandarizado
10.2.2 El algoritmo: mínimos cuadrados ponderados iterativos (IRLS)
Como los pesos dependen de los residuales y los residuales dependen de los coeficientes, se estima de forma iterativa:
Empezar con los coeficientes de OLS: b^{(0)} = (X'X)^{-1}X'y
Calcular residuales e = y - Xb y su escala robusta s = \text{mediana}(|e|)/0.6745 (MAD)
Calcular los pesos w_i con la función de Tukey y formar W = \text{diag}(w_1, \dots, w_n)
Re-estimar con mínimos cuadrados ponderados: b^{(t+1)} = (X'WX)^{-1}X'Wy
Repetir 2-4 hasta que los coeficientes ya no cambien
Lo implementamos con álgebra matricial sobre los datos winsorizados:
yw = m_B.model.endogb_rob = XtXi_w @ Xw.T @ yw # paso 1: arrancamos con OLSc_tukey =4.685for it inrange(200): e_rob = yw - Xw @ b_rob # paso 2 s_rob = np.median(np.abs(e_rob)) /0.6745 u_rob = e_rob / s_rob w_rob = np.where(np.abs(u_rob) <= c_tukey, (1- (u_rob / c_tukey)**2)**2, 0) # paso 3 W = np.diag(w_rob) b_new = np.linalg.inv(Xw.T @ W @ Xw) @ Xw.T @ W @ yw # paso 4if np.max(np.abs(b_new - b_rob)) <1e-8: # paso 5 b_rob = b_newbreak b_rob = b_newprint(f"Convergió en {it} iteraciones")
fig, ax = plt.subplots(figsize=(8, 4.5))sc = ax.scatter(dB["upa"], dB["roa"], c=dB["peso_tukey"], cmap="Blues", vmin=-0.15, vmax=1, s=14, edgecolor="#bbbbbb", lw=0.3)cb = plt.colorbar(sc, ax=ax); cb.set_label("Peso w (Tukey)")ax.set_xlabel("upa (winsorizada)"); ax.set_ylabel("ROA (winsorizado)")ax.set_title("Pesos de la regresión robusta (Tukey)")plt.tight_layout(); plt.show()
Peso de cada empresa en la regresión robusta (Tukey). Las más claras pesan menos
70 empresas (6.5%) reciben peso cero y 154 pesan menos de 0.5; la mayoría son empresas pequeñas con pérdidas muy grandes. Es el mismo grupo que identificaron los diagnósticos de la Lección 9, pero ahora no tuvimos que elegir un umbral para eliminarlas: el algoritmo decide cuánto pesa cada una.
ImportanteLimitación: la regresión robusta no protege contra leverage
Los estimadores M (Huber, Tukey) asignan los pesos según el residual. Una observación con leverage muy alto “jala” la recta hacia ella, termina con un residual pequeño y conserva su peso completo. Por eso aplicamos la regresión robusta después de winsorizar: la winsorización controla los valores extremos en X (el leverage máximo bajó de 0.79 a 0.15) y la regresión robusta controla los outliers en Y. Existen estimadores que atacan ambos problemas a la vez (MM, LTS), pero están fuera del alcance de este curso.
ROA esperado vs deuda (empresa pequeña con upa mediana) según cada método. La zona sombreada contiene al 2.3% de las empresas con deuda/activos > 0.8
Lectura de la comparación:
upa es positiva y muy significativa en todas las versiones (entre 0.74 y 0.84).
Mediana × upa: OLS con errores HC3 sólo la detecta al 10% (p = 0.096). La regresión robusta, que no se deja llevar por las empresas pequeñas extremas, la detecta con mucha claridad y con una magnitud mayor (0.34 con Huber y 0.50 con Tukey).
Grande × upa sólo es significativa al 5% con Tukey (p = 0.039).
Deuda: la forma de la relación no es estable. OLS (original y winsorizado) estima una parábola cóncava (término cuadrático negativo) con un máximo. La regresión de Huber estima una relación prácticamente lineal y negativa, y la de Tukey estima un término cuadrático positivo y significativo (parábola convexa, con un mínimo en DA ≈ 0.8). Este último resultado contradice lo que espera la teoría financiera y, además, en la zona donde la curva vuelve a subir (DA > 0.8) sólo está el 2.3% de las empresas.
11 Modelo final e interpretación
11.1 ¿Cuál es el modelo más adecuado?
Elegimos como modelo final el modelo OLS winsorizado (corrección B) con errores estándar robustos HC3. Las razones son:
Consistencia con la teoría. La teoría financiera (teoría del trade-off de la estructura de capital) predice un efecto cuadrático negativo del apalancamiento: niveles moderados de deuda pueden ayudar a la empresa (financiar proyectos productivos, disciplina en el uso de recursos, beneficio fiscal), pero a partir de cierto punto los costos de la deuda (intereses, riesgo de quiebra, restricciones para invertir) reducen la productividad. El modelo winsorizado estima justamente una parábola cóncava. En cambio, la regresión robusta de Tukey estima un término cuadrático positivo y significativo, que implicaría que, a partir de cierto nivel, más deuda aumenta la productividad; esto no tiene sustento teórico y además se basa en muy pocas empresas.
Representa a la empresa promedio. OLS estima el valor esperado (promedio) del ROA dadas las variables X, usando a todas las empresas. La regresión robusta, en cambio, describe a la empresa típica que sigue el patrón general, porque ignora (peso cero) a 70 empresas y le da poco peso a otras 84. La winsorización conserva a todas las empresas y sólo limita el efecto de los valores más extremos.
Atiende los problemas del diagnóstico. La winsorización controla el leverage (el leverage máximo bajó de 0.79 a 0.15) y reduce la influencia de los outliers; los errores HC3 corrigen la inferencia ante la heterocedasticidad que confirmó la prueba de Breusch-Pagan.
Con este modelo, los coeficientes son los de OLS winsorizado y toda la inferencia (errores estándar, p-values, intervalos y pruebas F) se hace con HC3.
Como la interpretación de las interacciones y del efecto cuadrático depende de varios coeficientes a la vez, hacemos pruebas F conjuntas usando la matriz de varianza HC3:
print("1) ¿El efecto de upa depende del tamaño? H0: b4 = b5 = 0")
1) ¿El efecto de upa depende del tamaño? H0: b4 = b5 = 0
Efecto marginal: dROA/dDA = 0.0961 + 2(-0.2029) DA = 0.0961 -0.4058 DA
print(f"ROA esperado máximo en DA = {da_opt_f:.3f} (mediana de DA en la muestra = {da_med:.3f})\n")
ROA esperado máximo en DA = 0.237 (mediana de DA en la muestra = 0.231)
for da_ in [0.0, 0.1, 0.237, 0.4, 0.6, 0.8, 1.0]: em = b6f +2* b7f * da_print(f" DA = {da_:5.3f}: +10pp de deuda => {em*0.10*100:6.2f} pp de ROA")
DA = 0.000: +10pp de deuda => 0.96 pp de ROA
DA = 0.100: +10pp de deuda => 0.55 pp de ROA
DA = 0.237: +10pp de deuda => -0.00 pp de ROA
DA = 0.400: +10pp de deuda => -0.66 pp de ROA
DA = 0.600: +10pp de deuda => -1.47 pp de ROA
DA = 0.800: +10pp de deuda => -2.29 pp de ROA
DA = 1.000: +10pp de deuda => -3.10 pp de ROA
Modelo final: ROA esperado vs deuda/activos (empresa pequeña con upa mediana) con intervalo de confianza al 95% (HC3)
11.5 Interpretación final
El modelo explica el 65.5% de la variabilidad del ROA (R² ajustada = 0.653).
upa (b_3 = 0.842, p < 0.001): para una empresa pequeña, un aumento de 0.10 en upa se asocia con +8.4 pp de ROA en promedio, manteniendo constante la deuda. Es el resultado más sólido del análisis: es positivo y muy significativo en todas las versiones del modelo, también con errores HC3.
Tamaño (diferencias cuando upa = 0, con la misma deuda): las empresas grandes tienen un ROA esperado 9.6 pp mayor que las pequeñas (b_2 = 0.096, p < 0.001) y las medianas 3.8 pp mayor (b_1 = 0.038, p = 0.002).
Interacción tamaño × upa: las pendientes estimadas son mayores en medianas (1.059) y grandes (1.081) que en pequeñas (0.842). Sin embargo, con errores HC3 la diferencia sólo es marginalmente significativa para las medianas (b_4 = 0.217, p = 0.096), no es significativa para las grandes (b_5 = 0.240, p = 0.185), y la prueba conjunta de ambas interacciones no rechaza H_0 (p = 0.12). Además, las pendientes de medianas y grandes son prácticamente iguales (p = 0.92). Conclusión: hay indicios de que la upa tiene un efecto mayor en las empresas medianas y grandes, pero la evidencia no es concluyente una vez que corregimos por heterocedasticidad.
Deuda: el término lineal es positivo (b_6 = 0.096) y el cuadrático negativo (b_7 = -0.202): la relación tiene forma de U invertida, consistente con la teoría. El ROA esperado es máximo cuando la deuda representa alrededor del 24% de los activos, prácticamente la mediana de la muestra (23%). Por debajo de ese nivel, más deuda se asocia con un ROA ligeramente mayor; por encima, más deuda se asocia con un ROA cada vez menor: para una empresa con 60% de deuda, 10 pp más de deuda se asocian con 1.5 pp menos de ROA, y con 100% de deuda, con 3.1 pp menos. En cuanto a la significancia, con errores HC3 el término cuadrático es significativo al 10% (p = 0.057) y la prueba conjunta de ambos términos de deuda también (p = 0.058). Es decir, la evidencia apoya la forma cóncava, pero con un nivel de confianza del 90%, no del 95%.
11.6 Comparación con el Challenge 4 del Workshop 1
Conclusión
W1 - Challenge 4 (OLS original)
Modelo final (winsorizado + HC3)
¿Se mantiene?
Efecto de upa
Positivo, 0.84 (p < 0.001)
Positivo, 0.84 (p < 0.001)
Sí
Grandes vs pequeñas
ROA mayor (p < 0.001)
ROA mayor (p < 0.001)
Sí
Medianas vs pequeñas
ROA mayor (p = 0.004)
ROA mayor (p = 0.002)
Sí
Mediana × upa
Positiva, significativa (p = 0.001)
Positiva, sólo al 10% (p = 0.096)
Parcialmente
Grande × upa
No significativa
No significativa
Sí
Deuda
Cóncava, máximo en DA ≈ 0.09, significativa al 1%
Cóncava, máximo en DA ≈ 0.24, significativa al 10%
Sí la forma; menor certeza
R²
0.594
0.655
Conclusiones:
Las conclusiones sobre upa y sobre las diferencias de ROA por tamaño son robustas: no dependen de las observaciones extremas ni del supuesto de varianza constante.
La interacción entre tamaño y upa, que en el Workshop 1 parecía clara, se debilita: su significancia se apoyaba en suponer varianza constante, cuando en realidad las empresas pequeñas tienen residuales mucho más dispersos.
El efecto cuadrático negativo de la deuda se mantiene, pero el punto óptimo se mueve de 9% a 24% de deuda sobre activos al controlar los valores extremos, y la evidencia es más débil (significativa al 10%).
Los diagnósticos de la Lección 9, la winsorización, los errores robustos y la regresión robusta se complementan: los diagnósticos se calculan primero sobre OLS para entender los datos; los errores robustos se deben reportar siempre que haya heterocedasticidad; y la regresión robusta sirve como prueba de robustez. Al elegir el modelo final no basta con el ajuste estadístico: también importa que los resultados sean consistentes con la teoría y que el modelo responda a la pregunta de interés (en este caso, el efecto en la empresa promedio).
