Taller 7: Análisis de Residuos en Modelos Lineales#
Objetivos de Aprendizaje#
Implementar pruebas de diagnóstico para residuos
Detectar violaciones de supuestos del modelo
Aplicar transformaciones correctivas
Comprender la importancia del análisis de residuos para la validez del modelo
Contenido#
Introducción
Marco Teórico
Supuestos del modelo lineal
Tipos de residuos
Pruebas formales para verificar supuestos
Implementación Computacional
Cálculo de diferentes tipos de residuos
Visualizaciones diagnósticas
Pruebas estadísticas para verificar supuestos
Transformaciones correctivas
Transformaciones para corregir heterocedasticidad
Transformaciones para corregir no normalidad
Transformaciones para corregir no linealidad
Ejemplos Prácticos
Ejercicios
# Configuración inicial e importación de librerías
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.stats as stats
from scipy import optimize
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error
import statsmodels.api as sm
from statsmodels.stats.diagnostic import het_breuschpagan, het_white, normal_ad
from statsmodels.stats.stattools import durbin_watson
from statsmodels.graphics.gofplots import ProbPlot
import warnings
# Configuración de visualización
plt.style.use('seaborn-v0_8-whitegrid')
sns.set_context("notebook", font_scale=1.2)
warnings.filterwarnings('ignore')
# Configuración para reproducibilidad
np.random.seed(123)
# Definimos una función para formatear mejor las tablas en Pandas
pd.set_option('display.float_format', '{:.4f}'.format)
pd.set_option('display.max_columns', None)
1. Introducción#
El análisis de residuos es una etapa crítica en el proceso de modelado estadístico, especialmente en los modelos lineales. Los residuos, definidos como la diferencia entre los valores observados y los valores predichos por el modelo, contienen información valiosa sobre la adecuación del modelo a los datos.
En los talleres anteriores hemos explorado aspectos clave de los modelos lineales como la estimación de parámetros, pruebas de hipótesis, multicolinealidad y métodos de selección de variables. Sin embargo, incluso un modelo con buenos estadísticos de ajuste (como \(R^2\) alto) puede ser inadecuado si no satisface los supuestos básicos del modelado lineal.
El análisis de residuos nos permite:
Verificar si se cumplen los supuestos del modelo lineal
Identificar patrones no explicados por el modelo
Detectar valores atípicos e influyentes
Guiar la mejora del modelo a través de transformaciones o especificaciones alternativas
En este taller, aprenderemos a realizar un análisis completo de residuos, implementando herramientas computacionales para el diagnóstico y aplicando transformaciones correctivas cuando sea necesario.
2. Marco Teórico#
2.1 Supuestos del Modelo Lineal Clásico#
Recordemos que el modelo lineal clásico puede expresarse como:
\[Y = X\beta + \varepsilon\]
Donde:
\(Y\) es el vector de respuestas (variable dependiente) de dimensión \(n \times 1\)
\(X\) es la matriz de diseño de dimensión \(n \times p\)
\(\beta\) es el vector de coeficientes de dimensión \(p \times 1\)
\(\varepsilon\) es el vector de errores aleatorios de dimensión \(n \times 1\)
Los supuestos básicos del modelo lineal clásico son:
Linealidad: \(E[\varepsilon_i] = 0\) o equivalentemente \(E[Y|X] = X\beta\)
La relación entre predictores y respuesta es lineal
No hay error sistemático en el modelo
Homocedasticidad: \(Var(\varepsilon_i) = \sigma^2\) para todo \(i\)
La varianza del error es constante para todas las observaciones
No depende de los valores de las variables predictoras o del valor esperado de Y
Independencia: \(Cov(\varepsilon_i, \varepsilon_j) = 0\) para todo \(i \neq j\)
Los errores no están correlacionados entre sí
No hay estructura temporal o espacial en los errores
Normalidad: \(\varepsilon \sim N(0, \sigma^2 I)\)
Los errores siguen una distribución normal
Este supuesto es necesario para la inferencia estadística (pruebas t, F, intervalos de confianza)
Cuando estos supuestos no se cumplen, las consecuencias pueden ser:
Estimadores sesgados o ineficientes
Errores estándar incorrectos
Inferencias estadísticas inválidas
Predicciones poco confiables
El análisis de residuos es nuestra principal herramienta para verificar si estos supuestos se cumplen en la práctica.
2.2 Tipos de Residuos#
Existen diferentes tipos de residuos que podemos analizar, cada uno con propósitos específicos:
2.2.1 Residuos Ordinarios o Crudos#
Son simplemente la diferencia entre los valores observados y ajustados:
\[r_i = y_i - \hat{y}_i = y_i - x_i^T\hat{\beta}\]
Estos residuos son útiles pero tienen limitaciones:
Su varianza no es constante (depende de la posición de \(x_i\) en el espacio de variables predictoras)
No siguen exactamente una distribución normal incluso cuando los errores sí lo hacen
Son sensibles a valores atípicos
2.2.2 Residuos Estandarizados#
Son los residuos ordinarios divididos por una estimación de su desviación estándar:
\[r_i^{std} = \frac{r_i}{\hat{\sigma}}\]
donde \(\hat{\sigma} = \sqrt{\frac{\sum_{i=1}^n r_i^2}{n-p}}\) es una estimación de la desviación estándar de los errores.
Estos residuos tienen varianza aproximadamente constante cuando se cumplen los supuestos del modelo.
2.2.3 Residuos Estudentizados Internos#
Son los residuos ordinarios divididos por su error estándar estimado, que tiene en cuenta la influencia de cada observación:
\[r_i^{stud} = \frac{r_i}{\hat{\sigma}\sqrt{1-h_{ii}}}\]
donde \(h_{ii}\) es el i-ésimo elemento diagonal de la matriz de proyección o “hat matrix” \(H = X(X^TX)^{-1}X^T\).
Estos residuos tienen en cuenta que las observaciones con valores extremos en las variables predictoras tienen mayor varianza en sus residuos.
2.2.4 Residuos Estudentizados Externos (o Eliminados)#
Se calculan ajustando el modelo sin la i-ésima observación y comparando la predicción para esa observación con su valor real:
\[t_i = \frac{r_i}{\hat{\sigma}_{(i)}\sqrt{1-h_{ii}}}\]
donde \(\hat{\sigma}_{(i)}\) es la estimación de la desviación estándar calculada sin la i-ésima observación.
Estos residuos siguen una distribución t de Student con n-p-1 grados de libertad bajo los supuestos del modelo lineal, lo que los hace útiles para detectar valores atípicos.
2.2.5 Residuos Parciales#
Se utilizan para examinar la relación entre una variable predictora específica y la variable respuesta, después de ajustar por los efectos de las otras variables:
\[r_{ij}^{partial} = r_i + \hat{\beta}_j x_{ij}\]
Estos residuos son útiles para detectar no linealidades en la relación entre predictores individuales y la respuesta.
2.3 Pruebas Formales para Verificar Supuestos#
Además de las visualizaciones, existen pruebas estadísticas para verificar formalmente los supuestos del modelo:
2.3.1 Pruebas de Normalidad#
Test de Shapiro-Wilk:
\(H_0\): Los residuos siguen una distribución normal
\(H_1\): Los residuos no siguen una distribución normal
Estadístico \(W\) basado en la correlación entre los residuos y los cuantiles normales
Más potente para muestras pequeñas (n < 50)
Test de Anderson-Darling:
Similar al Shapiro-Wilk pero con mayor énfasis en las colas de la distribución
Más sensible a desviaciones en las colas que otras pruebas
Test de Jarque-Bera:
Basado en la asimetría y curtosis de la distribución de residuos
\(JB = \frac{n}{6}\left(S^2 + \frac{(K-3)^2}{4}\right)\)
Donde \(S\) es la asimetría y \(K\) es la curtosis
Bajo \(H_0\), sigue una distribución chi-cuadrado con 2 grados de libertad
2.3.2 Pruebas de Homocedasticidad#
Test de Breusch-Pagan:
\(H_0\): La varianza de los errores es constante (homocedasticidad)
\(H_1\): La varianza de los errores no es constante (heterocedasticidad)
Procedimiento:
Ajustar el modelo y obtener los residuos \(r_i\)
Calcular \(z_i = r_i^2 / \hat{\sigma}^2\)
Regresar \(z_i\) sobre las variables predictoras
El estadístico de prueba es \(BP = n \cdot R^2\) de esta regresión auxiliar
Bajo \(H_0\), \(BP \sim \chi^2_{p-1}\) donde \(p\) es el número de predictores
Test de White:
Similar al Breusch-Pagan pero incluye términos cuadráticos e interacciones
Más general, no asume una forma específica de heterocedasticidad
Más potente cuando la heterocedasticidad está relacionada con interacciones o términos no lineales de los predictores
2.3.3 Pruebas de Independencia#
Test de Durbin-Watson:
Diseñado para detectar autocorrelación de primer orden en los residuos
\(H_0\): No hay autocorrelación de primer orden
\(H_1\): Existe autocorrelación de primer orden
Estadístico \(DW = \frac{\sum_{i=2}^n (r_i - r_{i-1})^2}{\sum_{i=1}^n r_i^2}\)
\(DW \approx 2\) sugiere no autocorrelación
\(DW < 2\) sugiere autocorrelación positiva
\(DW > 2\) sugiere autocorrelación negativa
Test de Breusch-Godfrey:
Más general que Durbin-Watson, puede detectar autocorrelación de orden superior
\(H_0\): No hay autocorrelación hasta el orden \(p\)
\(H_1\): Existe autocorrelación de algún orden hasta \(p\)
2.3.4 Pruebas de Linealidad#
Test RESET de Ramsey:
Prueba para detectar incorrecta especificación del modelo
\(H_0\): El modelo está correctamente especificado
\(H_1\): El modelo está incorrectamente especificado
Procedimiento:
Ajustar el modelo original y obtener \(\hat{y}\)
Ajustar un modelo aumentado incluyendo potencias de \(\hat{y}\) como predictores adicionales
Realizar una prueba F comparando los dos modelos
3. Implementación Computacional#
3.1 Cálculo de Diferentes Tipos de Residuos#
A continuación, implementaremos funciones para calcular los diferentes tipos de residuos que hemos discutido:
def calcular_residuos(model, X, y):
"""
Calcula diferentes tipos de residuos para un modelo lineal.
Parámetros:
-----------
model : objeto modelo de statsmodels con método 'fit'
El modelo lineal ajustado
X : array-like
Variables predictoras
y : array-like
Variable respuesta
Retorna:
--------
dict : Diccionario con diferentes tipos de residuos
"""
# Aseguramos que X tenga un intercept
X_with_const = sm.add_constant(X) if 'const' not in X.columns else X
# Ajustamos el modelo
fitted_model = model.fit()
# Obtenemos predicciones y residuos ordinarios
y_pred = fitted_model.predict(X_with_const)
residuos_ordinarios = y - y_pred
# Calcular la matriz hat (matriz de proyección)
hat_matrix = X_with_const @ np.linalg.inv(X_with_const.T @ X_with_const) @ X_with_const.T
leverage = np.diag(hat_matrix)
# Grados de libertad y sigma estimado
n, p = X_with_const.shape
sigma_estimado = np.sqrt(np.sum(residuos_ordinarios**2) / (n - p))
# Residuos estandarizados
residuos_estandarizados = residuos_ordinarios / sigma_estimado
# Residuos estudentizados internos
residuos_estudentizados = residuos_ordinarios / (sigma_estimado * np.sqrt(1 - leverage))
# Residuos estudentizados externos (aproximación)
# Para calcularlos correctamente, necesitaríamos ajustar el modelo n veces
# Esta es una aproximación utilizando la fórmula de actualización
sigma_i = np.sqrt(((n - p) * sigma_estimado**2 - residuos_ordinarios**2 / (1 - leverage)) / (n - p - 1))
residuos_estudentizados_externos = residuos_ordinarios / (sigma_i * np.sqrt(1 - leverage))
# Residuos parciales para cada predictor (excluyendo la constante)
residuos_parciales = {}
if len(X.columns) > 0: # Si hay predictores (además de la constante)
for j, predictor in enumerate(X.columns):
beta_j = fitted_model.params[predictor]
residuos_parciales[predictor] = residuos_ordinarios + beta_j * X[predictor]
return {
'ordinarios': residuos_ordinarios,
'estandarizados': residuos_estandarizados,
'estudentizados_internos': residuos_estudentizados,
'estudentizados_externos': residuos_estudentizados_externos,
'parciales': residuos_parciales,
'leverage': leverage,
'fitted_values': y_pred,
'model': fitted_model
}
Implementemos también una función para generar datos de prueba con diferentes violaciones de supuestos:
def generar_datos_prueba(n=100, p=3, problema='normal', intensidad=1.0):
"""
Genera datos de prueba con diferentes problemas para el análisis de residuos.
Parámetros:
-----------
n : int
Número de observaciones
p : int
Número de predictores
problema : str
Tipo de problema a simular: 'normal', 'heterocedasticidad', 'no_linealidad',
'autocorrelacion', 'no_normalidad', 'outliers'
intensidad : float
Controla la intensidad del problema (0.0 a 2.0, donde 1.0 es moderado)
Retorna:
--------
X : DataFrame
Variables predictoras
y : Series
Variable respuesta
"""
# Generamos predictores
X = pd.DataFrame(np.random.normal(0, 1, size=(n, p)),
columns=[f'X{i+1}' for i in range(p)])
# Coeficientes verdaderos
beta = np.array([0.5, 1.0, -0.7][:p])
# Término de error base (normal estándar)
epsilon = np.random.normal(0, 1, size=n)
# Aplicamos diferentes problemas
if problema == 'normal':
# Caso normal (sin problemas)
y = X @ beta + epsilon
elif problema == 'heterocedasticidad':
# Varianza creciente con X1
factor_varianza = (1 + intensidad * np.abs(X['X1']))
epsilon_hetero = epsilon * factor_varianza
y = X @ beta + epsilon_hetero
elif problema == 'no_linealidad':
# Relación cuadrática con X1
termino_cuadratico = intensidad * (X['X1']**2)
y = X @ beta + termino_cuadratico + epsilon
elif problema == 'autocorrelacion':
# Autocorrelación de primer orden
epsilon_auto = np.zeros(n)
epsilon_auto[0] = epsilon[0]
for i in range(1, n):
epsilon_auto[i] = intensidad * 0.8 * epsilon_auto[i-1] + epsilon[i]
y = X @ beta + epsilon_auto
elif problema == 'no_normalidad':
# Error con distribución asimétrica (chi-cuadrado)
df = 3 # Grados de libertad para chi-cuadrado
epsilon_no_normal = np.random.chisquare(df, size=n) - df
epsilon_no_normal = intensidad * epsilon_no_normal + (1 - intensidad) * epsilon
y = X @ beta + epsilon_no_normal / np.std(epsilon_no_normal)
elif problema == 'outliers':
# Añadir algunos valores atípicos
n_outliers = int(0.05 * n * intensidad)
outlier_idx = np.random.choice(range(n), n_outliers, replace=False)
epsilon_out = epsilon.copy()
epsilon_out[outlier_idx] = epsilon_out[outlier_idx] + np.random.choice([-1, 1], n_outliers) * np.random.uniform(4, 7, n_outliers)
y = X @ beta + epsilon_out
return X, pd.Series(y, name='y')
3.2 Visualizaciones Diagnósticas#
Las visualizaciones son herramientas esenciales para el análisis de residuos. Implementaremos varias funciones para crear gráficos diagnósticos comunes:
def plot_diagnostico_residuos(resultados_residuos, figsize=(18, 12)):
"""
Crea un conjunto de gráficos diagnósticos para analizar residuos.
Parámetros:
-----------
resultados_residuos : dict
Diccionario con residuos y valores ajustados (output de calcular_residuos)
figsize : tuple
Tamaño de la figura
"""
residuos = resultados_residuos['estudentizados_internos']
y_pred = resultados_residuos['fitted_values']
leverage = resultados_residuos['leverage']
modelo = resultados_residuos['model']
fig, axes = plt.subplots(2, 2, figsize=figsize)
# 1. Residuos vs. Valores Ajustados
axes[0, 0].scatter(y_pred, residuos, alpha=0.7)
axes[0, 0].axhline(y=0, color='r', linestyle='--')
# Añadir línea de tendencia suavizada
try:
sns.regplot(x=y_pred, y=residuos, lowess=True, scatter=False,
color='g', line_kws={'linewidth': 1}, ax=axes[0, 0])
except:
pass # Si lowess falla
axes[0, 0].set_xlabel('Valores Ajustados')
axes[0, 0].set_ylabel('Residuos Estudentizados')
axes[0, 0].set_title('Residuos vs. Valores Ajustados')
# Agregar límites de referencia en ±2 y ±3
axes[0, 0].axhline(y=2, color='orange', linestyle=':', alpha=0.7)
axes[0, 0].axhline(y=-2, color='orange', linestyle=':', alpha=0.7)
axes[0, 0].axhline(y=3, color='red', linestyle=':', alpha=0.7)
axes[0, 0].axhline(y=-3, color='red', linestyle=':', alpha=0.7)
# 2. QQ Plot
QQ = ProbPlot(residuos)
QQ.qqplot(line='45', ax=axes[0, 1])
axes[0, 1].set_title('Q-Q Plot')
# 3. Scale-Location Plot (raíz cuadrada de residuos abs. estandarizados vs. valores ajustados)
sqrt_abs_residuos = np.sqrt(np.abs(residuos))
axes[1, 0].scatter(y_pred, sqrt_abs_residuos, alpha=0.7)
try:
sns.regplot(x=y_pred, y=sqrt_abs_residuos, lowess=True, scatter=False,
color='g', line_kws={'linewidth': 1}, ax=axes[1, 0])
except:
pass
axes[1, 0].set_xlabel('Valores Ajustados')
axes[1, 0].set_ylabel('√|Residuos Estudentizados|')
axes[1, 0].set_title('Scale-Location')
# 4. Leverage vs. Residuos Estudentizados al Cuadrado
# Calculamos la distancia de Cook
p = len(modelo.params)
n = len(residuos)
residuos_est_sq = residuos ** 2
cook_dist = residuos_est_sq * leverage / (p * (1 - leverage))
axes[1, 1].scatter(leverage, residuos_est_sq, alpha=0.7)
axes[1, 1].set_xlabel('Leverage')
axes[1, 1].set_ylabel('Residuos Estudentizados²')
axes[1, 1].set_title('Leverage vs. Residuos²')
# Añadir contornos para la distancia de Cook
max_leverage = max(leverage) * 1.1
max_residuos_sq = max(residuos_est_sq) * 1.1
leverage_grid = np.linspace(0.001, max_leverage, 100)
for d in [0.5, 1.0]:
residuos_sq_grid = d * p * (1 - leverage_grid) / leverage_grid
valid_idx = residuos_sq_grid > 0
axes[1, 1].plot(leverage_grid[valid_idx], residuos_sq_grid[valid_idx],
'r--', alpha=0.6, label=f'Cook\'s D = {d}')
# Línea vertical en 2p/n y 3p/n (reglas prácticas para leverage alto)
axes[1, 1].axvline(x=2*p/n, color='orange', linestyle=':', alpha=0.7, label='Leverage = 2p/n')
axes[1, 1].axvline(x=3*p/n, color='red', linestyle=':', alpha=0.7, label='Leverage = 3p/n')
axes[1, 1].legend()
plt.tight_layout()
plt.show()
# Gráficos adicionales: Residuos Parciales
residuos_parciales = resultados_residuos['parciales']
if residuos_parciales and len(residuos_parciales) > 0:
n_predictores = len(residuos_parciales)
fig, axes = plt.subplots(1, n_predictores, figsize=(5*n_predictores, 5))
if n_predictores == 1:
axes = [axes] # Convertir a lista si solo hay un predictor
for i, (predictor, residuos_parc) in enumerate(residuos_parciales.items()):
X_predictor = modelo.model.exog[:, modelo.model.data.param_names.index(predictor)]
axes[i].scatter(X_predictor, residuos_parc, alpha=0.7)
try:
sns.regplot(x=X_predictor, y=residuos_parc, lowess=True, scatter=False,
color='g', line_kws={'linewidth': 1}, ax=axes[i])
except:
pass
axes[i].set_xlabel(f'{predictor}')
axes[i].set_ylabel(f'Residuos Parciales')
axes[i].set_title(f'Residuos Parciales para {predictor}')
plt.tight_layout()
plt.show()
def plot_residuos_vs_predictores(resultados_residuos, X,figsize=(15,16)):
figsize=(15, 4*len(X.columns))
"""
Crea gráficos de residuos versus cada predictor.
Parámetros:
-----------
resultados_residuos : dict
Diccionario con residuos y valores ajustados
X : DataFrame
Variables predictoras
figsize : tuple
Tamaño de la figura
"""
residuos = resultados_residuos['estudentizados_internos']
n_predictores = len(X.columns)
fig, axes = plt.subplots(n_predictores, 1, figsize=figsize)
if n_predictores == 1:
axes = [axes] # Convertir a lista si solo hay un predictor
for i, predictor in enumerate(X.columns):
axes[i].scatter(X[predictor], residuos, alpha=0.7)
axes[i].axhline(y=0, color='r', linestyle='--')
# Añadir línea de tendencia suavizada
try:
sns.regplot(x=X[predictor], y=residuos, lowess=True, scatter=False,
color='g', line_kws={'linewidth': 1}, ax=axes[i])
except:
pass
axes[i].set_xlabel(predictor)
axes[i].set_ylabel('Residuos Estudentizados')
axes[i].set_title(f'Residuos vs. {predictor}')
# Agregar límites de referencia en ±2 y ±3
axes[i].axhline(y=2, color='orange', linestyle=':', alpha=0.7)
axes[i].axhline(y=-2, color='orange', linestyle=':', alpha=0.7)
axes[i].axhline(y=3, color='red', linestyle=':', alpha=0.7)
axes[i].axhline(y=-3, color='red', linestyle=':', alpha=0.7)
plt.tight_layout()
plt.show()
def plot_secuencia_temporal(resultados_residuos, figsize=(12, 6)):
"""
Crea un gráfico de secuencia temporal de los residuos (útil para detectar autocorrelación).
Parámetros:
-----------
resultados_residuos : dict
Diccionario con residuos y valores ajustados
figsize : tuple
Tamaño de la figura
"""
residuos = resultados_residuos['estudentizados_internos']
plt.figure(figsize=figsize)
plt.plot(residuos, marker='o', linestyle='-', alpha=0.7)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel('Índice de Observación')
plt.ylabel('Residuos Estudentizados')
plt.title('Secuencia Temporal de Residuos')
# Agregar límites de referencia en ±2 y ±3
# Completando la función plot_secuencia_temporal que quedó incompleta
plt.axhline(y=2, color='orange', linestyle=':', alpha=0.7)
plt.axhline(y=-2, color='orange', linestyle=':', alpha=0.7)
plt.axhline(y=3, color='red', linestyle=':', alpha=0.7)
plt.axhline(y=-3, color='red', linestyle=':', alpha=0.7)
plt.tight_layout()
plt.show()
# Podemos agregar también una función para visualizar la autocorrelación de los residuos:
def plot_autocorrelacion(resultados_residuos, lags=20, figsize=(12, 8)):
"""
Crea gráficos de autocorrelación y autocorrelación parcial de los residuos.
Parámetros:
-----------
resultados_residuos : dict
Diccionario con residuos y valores ajustados
lags : int
Número de rezagos a mostrar
figsize : tuple
Tamaño de la figura
"""
residuos = resultados_residuos['estudentizados_internos']
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=figsize)
# Autocorrelación
sm.graphics.tsa.plot_acf(residuos, lags=lags, ax=ax1)
ax1.set_title('Función de Autocorrelación (ACF)')
# Autocorrelación parcial
sm.graphics.tsa.plot_pacf(residuos, lags=lags, ax=ax2)
ax2.set_title('Función de Autocorrelación Parcial (PACF)')
plt.tight_layout()
plt.show()
3.3 Pruebas Estadísticas para Verificar Supuestos#
Implementaremos funciones para realizar las pruebas estadísticas formales que verifican los supuestos del modelo lineal:
def test_normalidad(residuos):
"""
Realiza pruebas de normalidad sobre los residuos.
Parámetros:
-----------
residuos : array-like
Residuos del modelo
Retorna:
--------
dict : Diccionario con resultados de las pruebas
"""
# Test de Shapiro-Wilk
if len(residuos) < 5000: # Limitación del test
stat_sw, p_sw = stats.shapiro(residuos)
else:
stat_sw, p_sw = np.nan, np.nan
# Test de Anderson-Darling
result_ad = stats.anderson(residuos, dist='norm')
stat_ad = result_ad.statistic
# Verificamos si el estadístico supera los valores críticos
p_ad = "< 0.01" # Por defecto
for i, critical_val in enumerate(result_ad.critical_values):
if stat_ad < critical_val:
significancias = [15, 10, 5, 2.5, 1]
p_ad = f"> {significancias[i]}%"
break
# Test de Jarque-Bera
stat_jb, p_jb = stats.jarque_bera(residuos)
# Test de Normalidad de D'Agostino y Pearson
stat_dp, p_dp = stats.normaltest(residuos)
# Calcular asimetría y curtosis
asimetria = stats.skew(residuos)
curtosis = stats.kurtosis(residuos)
return {
'shapiro': {'estadistico': stat_sw, 'p_valor': p_sw},
'anderson_darling': {'estadistico': stat_ad, 'significancia': p_ad},
'jarque_bera': {'estadistico': stat_jb, 'p_valor': p_jb},
'dagostino_pearson': {'estadistico': stat_dp, 'p_valor': p_dp},
'asimetria': asimetria,
'curtosis': curtosis
}
def test_homocedasticidad(modelo, X):
"""
Realiza pruebas de homocedasticidad sobre los residuos.
Parámetros:
-----------
modelo : objeto modelo ajustado de statsmodels
El modelo lineal ajustado
X : array-like
Variables predictoras
Retorna:
--------
dict : Diccionario con resultados de las pruebas
"""
# Test de Breusch-Pagan
bp_test = het_breuschpagan(modelo.resid, X)
labels_bp = ['Estadístico LM', 'p-valor LM', 'Estadístico F', 'p-valor F']
# Test de White
# Necesitamos generar términos cuadráticos e interacciones
X_white = X.copy()
# Añadir términos cuadráticos
for col in X.columns:
X_white[f'{col}^2'] = X[col]**2
# Añadir interacciones (productos cruzados)
for i, col1 in enumerate(X.columns):
for col2 in X.columns[i+1:]:
X_white[f'{col1}*{col2}'] = X[col1] * X[col2]
white_test = het_white(modelo.resid, X_white)
labels_white = ['Estadístico LM', 'p-valor LM', 'Estadístico F', 'p-valor F']
# Calcular residuos cuadrados y correlación con valores ajustados
residuos_cuadrados = modelo.resid**2
valores_ajustados = modelo.fittedvalues
corr_res2_yhat, p_corr = stats.pearsonr(residuos_cuadrados, valores_ajustados)
return {
'breusch_pagan': {labels_bp[i]: bp_test[i] for i in range(len(labels_bp))},
'white': {labels_white[i]: white_test[i] for i in range(len(labels_white))},
'corr_res2_yhat': {'correlacion': corr_res2_yhat, 'p_valor': p_corr}
}
def test_independencia(residuos, max_lag=10):
"""
Realiza pruebas de independencia (no autocorrelación) sobre los residuos.
Parámetros:
-----------
residuos : array-like
Residuos del modelo
max_lag : int
Número máximo de rezagos a considerar
Retorna:
--------
dict : Diccionario con resultados de las pruebas
"""
# Test de Durbin-Watson
dw_stat = durbin_watson(residuos)
# Autocorrelaciones y p-valores
autocorrelaciones = {}
p_valores = {}
for lag in range(1, min(max_lag+1, len(residuos)//5)): # Limitamos a n/5
# Autocorrelación de orden lag
residuos_lag = residuos[lag:]
residuos_orig = residuos[:-lag]
if len(residuos_lag) > 1: # Aseguramos suficientes observaciones
corr, p_valor = stats.pearsonr(residuos_orig, residuos_lag)
autocorrelaciones[lag] = corr
p_valores[lag] = p_valor
# Test Box-Pierce y Ljung-Box
bp_stats = {}
bp_pvals = {}
lb_stats = {}
lb_pvals = {}
for lag in range(1, min(max_lag+1, len(residuos)//5)):
bp_stat, bp_pval = sm.stats.acorr_ljungbox(residuos, lags=[lag], boxpierce=True, return_df=False)
lb_stat, lb_pval = sm.stats.acorr_ljungbox(residuos, lags=[lag], boxpierce=False, return_df=False)
bp_stats[lag] = bp_stat[0]
bp_pvals[lag] = bp_pval[0]
lb_stats[lag] = lb_stat[0]
lb_pvals[lag] = lb_pval[0]
return {
'durbin_watson': dw_stat,
'autocorrelaciones': autocorrelaciones,
'p_valores': p_valores,
'box_pierce': {'estadisticos': bp_stats, 'p_valores': bp_pvals},
'ljung_box': {'estadisticos': lb_stats, 'p_valores': lb_pvals}
}
def test_linealidad(modelo, X, y):
"""
Realiza pruebas de linealidad (especificación correcta) sobre el modelo.
Parámetros:
-----------
modelo : objeto modelo ajustado de statsmodels
El modelo lineal ajustado
X : array-like
Variables predictoras
y : array-like
Variable respuesta
Retorna:
--------
dict : Diccionario con resultados de las pruebas
"""
# Test RESET de Ramsey
# Paso 1: Obtener predicciones del modelo original
y_pred = modelo.predict()
# Paso 2: Crear nuevas variables con potencias de las predicciones
X_reset = X.copy()
# Añadimos potencias de y_pred como predictores adicionales
X_reset['y_pred^2'] = y_pred**2
X_reset['y_pred^3'] = y_pred**3
# Paso 3: Ajustar un nuevo modelo con las variables originales más las potencias
X_reset_const = sm.add_constant(X_reset)
modelo_reset = sm.OLS(y, X_reset_const).fit()
# Paso 4: Realizar prueba F comparando los dos modelos
# Calculamos el estadístico F manualmente
n = len(y)
p_original = len(modelo.params)
p_reset = len(modelo_reset.params)
rss_original = np.sum(modelo.resid**2)
rss_reset = np.sum(modelo_reset.resid**2)
df1 = p_reset - p_original # Grados de libertad del numerador
df2 = n - p_reset # Grados de libertad del denominador
f_stat = ((rss_original - rss_reset) / df1) / (rss_reset / df2)
p_valor = 1 - stats.f.cdf(f_stat, df1, df2)
# Correlación entre residuos y valores ajustados al cuadrado
corr_resid_pred2, p_corr = stats.pearsonr(modelo.resid, y_pred**2)
return {
'ramsey_reset': {'estadistico_F': f_stat, 'p_valor': p_valor, 'df1': df1, 'df2': df2},
'corr_resid_pred2': {'correlacion': corr_resid_pred2, 'p_valor': p_corr}
}
def diagnostico_completo(modelo, X, y):
"""
Realiza un diagnóstico completo del modelo verificando todos los supuestos.
Parámetros:
-----------
modelo : objeto modelo ajustado de statsmodels
El modelo lineal ajustado
X : array-like
Variables predictoras
y : array-like
Variable respuesta
Retorna:
--------
dict : Diccionario con resultados de todas las pruebas
"""
# Calcular residuos
resultados_residuos = {}
resultados_residuos['ordinarios'] = modelo.resid
resultados_residuos['estudentizados_internos'] = modelo.get_influence().resid_studentized_internal
resultados_residuos['fitted_values'] = modelo.fittedvalues
resultados_residuos['leverage'] = modelo.get_influence().hat_matrix_diag
resultados_residuos['model'] = modelo
# Pruebas de supuestos
resultados_normalidad = test_normalidad(resultados_residuos['estudentizados_internos'])
resultados_homocedasticidad = test_homocedasticidad(modelo, X)
resultados_independencia = test_independencia(resultados_residuos['ordinarios'])
resultados_linealidad = test_linealidad(modelo, X, y)
return {
'residuos': resultados_residuos,
'normalidad': resultados_normalidad,
'homocedasticidad': resultados_homocedasticidad,
'independencia': resultados_independencia,
'linealidad': resultados_linealidad
}
def imprimir_resultados_diagnostico(resultados_diagnostico):
"""
Imprime los resultados del diagnóstico de forma legible.
Parámetros:
-----------
resultados_diagnostico : dict
Diccionario con resultados del diagnóstico
"""
print("="*80)
print("DIAGNÓSTICO DE RESIDUOS DEL MODELO LINEAL")
print("="*80)
# Normalidad
print("\n" + "-"*30 + " PRUEBAS DE NORMALIDAD " + "-"*30)
norm = resultados_diagnostico['normalidad']
print(f"Asimetría: {norm['asimetria']:.4f}")
print(f"Curtosis: {norm['curtosis']:.4f}")
print("\nTest de Shapiro-Wilk:")
print(f" Estadístico: {norm['shapiro']['estadistico']:.4f}")
print(f" p-valor: {norm['shapiro']['p_valor']:.4f}")
print(f" Conclusión: {'Rechazar H0 (no normal)' if norm['shapiro']['p_valor'] < 0.05 else 'No rechazar H0 (normal)'}")
print("\nTest de Jarque-Bera:")
print(f" Estadístico: {norm['jarque_bera']['estadistico']:.4f}")
print(f" p-valor: {norm['jarque_bera']['p_valor']:.4f}")
print(f" Conclusión: {'Rechazar H0 (no normal)' if norm['jarque_bera']['p_valor'] < 0.05 else 'No rechazar H0 (normal)'}")
print("\nTest de D'Agostino-Pearson:")
print(f" Estadístico: {norm['dagostino_pearson']['estadistico']:.4f}")
print(f" p-valor: {norm['dagostino_pearson']['p_valor']:.4f}")
print(f" Conclusión: {'Rechazar H0 (no normal)' if norm['dagostino_pearson']['p_valor'] < 0.05 else 'No rechazar H0 (normal)'}")
# Homocedasticidad
print("\n" + "-"*30 + " PRUEBAS DE HOMOCEDASTICIDAD " + "-"*30)
homo = resultados_diagnostico['homocedasticidad']
print("\nTest de Breusch-Pagan:")
print(f" Estadístico LM: {homo['breusch_pagan']['Estadístico LM']:.4f}")
print(f" p-valor: {homo['breusch_pagan']['p-valor LM']:.4f}")
print(f" Conclusión: {'Rechazar H0 (heterocedasticidad)' if homo['breusch_pagan']['p-valor LM'] < 0.05 else 'No rechazar H0 (homocedasticidad)'}")
print("\nTest de White:")
print(f" Estadístico LM: {homo['white']['Estadístico LM']:.4f}")
print(f" p-valor: {homo['white']['p-valor LM']:.4f}")
print(f" Conclusión: {'Rechazar H0 (heterocedasticidad)' if homo['white']['p-valor LM'] < 0.05 else 'No rechazar H0 (homocedasticidad)'}")
# Independencia
print("\n" + "-"*30 + " PRUEBAS DE INDEPENDENCIA " + "-"*30)
indep = resultados_diagnostico['independencia']
print(f"\nEstadístico Durbin-Watson: {indep['durbin_watson']:.4f}")
if indep['durbin_watson'] < 1.5:
conclusion_dw = "Posible autocorrelación positiva"
elif indep['durbin_watson'] > 2.5:
conclusion_dw = "Posible autocorrelación negativa"
else:
conclusion_dw = "No hay evidencia de autocorrelación de primer orden"
print(f" Conclusión: {conclusion_dw}")
# Mostrar autocorrelaciones significativas
print("\nAutocorrelaciones significativas (p < 0.05):")
hay_significativas = False
for lag, p_valor in indep['p_valores'].items():
if p_valor < 0.05:
hay_significativas = True
print(f" Lag {lag}: {indep['autocorrelaciones'][lag]:.4f} (p-valor: {p_valor:.4f})")
if not hay_significativas:
print(" Ninguna autocorrelación significativa")
# Linealidad
print("\n" + "-"*30 + " PRUEBAS DE LINEALIDAD " + "-"*30)
lineal = resultados_diagnostico['linealidad']
print("\nTest RESET de Ramsey:")
print(f" Estadístico F: {lineal['ramsey_reset']['estadistico_F']:.4f}")
print(f" p-valor: {lineal['ramsey_reset']['p_valor']:.4f}")
print(f" Conclusión: {'Rechazar H0 (especificación incorrecta)' if lineal['ramsey_reset']['p_valor'] < 0.05 else 'No rechazar H0 (especificación correcta)'}")
# Resumen general
print("\n" + "="*30 + " RESUMEN DEL DIAGNÓSTICO " + "="*30)
problemas = []
# Verificar normalidad
if (norm['shapiro']['p_valor'] < 0.05 and not np.isnan(norm['shapiro']['p_valor'])) or \
norm['jarque_bera']['p_valor'] < 0.05 or \
norm['dagostino_pearson']['p_valor'] < 0.05:
problemas.append("No normalidad de los residuos")
# Verificar homocedasticidad
if homo['breusch_pagan']['p-valor LM'] < 0.05 or homo['white']['p-valor LM'] < 0.05:
problemas.append("Heterocedasticidad")
# Verificar independencia
if indep['durbin_watson'] < 1.5 or indep['durbin_watson'] > 2.5 or hay_significativas:
problemas.append("Autocorrelación en los residuos")
# Verificar linealidad
if lineal['ramsey_reset']['p_valor'] < 0.05:
problemas.append("Especificación incorrecta del modelo (no linealidad)")
if problemas:
print("\nSe detectaron los siguientes problemas:")
for i, problema in enumerate(problemas):
print(f" {i+1}. {problema}")
else:
print("\nNo se detectaron problemas graves en los supuestos del modelo.")
4. Transformaciones Correctivas#
Cuando se detectan violaciones a los supuestos del modelo, es necesario aplicar transformaciones correctivas. En esta sección, implementaremos y analizaremos diferentes transformaciones para corregir problemas específicos.
4.1 Transformaciones para Corregir Heterocedasticidad#
La heterocedasticidad ocurre cuando la varianza de los errores no es constante. Las transformaciones más comunes para corregir este problema son:
def transformaciones_heterocedasticidad(X, y, metodo='log'):
"""
Aplica transformaciones para corregir heterocedasticidad.
Parámetros:
-----------
X : DataFrame
Variables predictoras
y : Series
Variable respuesta
metodo : str
Tipo de transformación: 'log', 'sqrt', 'reciproco', 'box-cox'
Retorna:
--------
tuple : (X_transformado, y_transformado, lambda_y)
"""
# Hacemos copias para no modificar los originales
X_trans = X.copy()
# Para Box-Cox necesitamos datos positivos
min_y = min(y)
if min_y <= 0:
y_adj = y - min_y + 0.01 # Ajustamos para asegurar positividad
else:
y_adj = y.copy()
if metodo == 'log':
# Transformación logarítmica (útil cuando la varianza aumenta con la media)
y_trans = np.log(y_adj)
lambda_y = 0 # Log equivale a Box-Cox con lambda = 0
elif metodo == 'sqrt':
# Transformación raíz cuadrada (menos drástica que log)
y_trans = np.sqrt(y_adj)
lambda_y = 0.5 # Sqrt equivale a Box-Cox con lambda = 0.5
elif metodo == 'reciproco':
# Transformación recíproca (1/y)
y_trans = 1 / y_adj
lambda_y = -1 # Recíproco equivale a Box-Cox con lambda = -1
elif metodo == 'box-cox':
# Transformación Box-Cox (busca el lambda óptimo)
y_trans, lambda_y = stats.boxcox(y_adj)
print(f"Lambda óptimo para Box-Cox: {lambda_y:.4f}")
else:
raise ValueError("Método no reconocido. Use 'log', 'sqrt', 'reciproco' o 'box-cox'")
return X_trans, pd.Series(y_trans, index=y.index, name=y.name), lambda_y
def comparar_modelos_heterocedasticidad(X, y):
"""
Compara diferentes transformaciones para corregir heterocedasticidad.
Parámetros:
-----------
X : DataFrame
Variables predictoras
y : Series
Variable respuesta
Retorna:
--------
DataFrame : Resultados comparativos
"""
# Modelo original
X_const = sm.add_constant(X)
modelo_original = sm.OLS(y, X_const).fit()
resultados_original = diagnostico_completo(modelo_original, X, y)
# Modelos con transformaciones
transformaciones = ['log', 'sqrt', 'reciproco', 'box-cox']
resultados = []
for metodo in transformaciones:
try:
X_trans, y_trans, lambda_y = transformaciones_heterocedasticidad(X, y, metodo)
X_trans_const = sm.add_constant(X_trans)
modelo_trans = sm.OLS(y_trans, X_trans_const).fit()
# Diagnóstico del modelo transformado
diag_trans = diagnostico_completo(modelo_trans, X_trans, y_trans)
# Extraer p-valores de las pruebas de heterocedasticidad
p_bp = diag_trans['homocedasticidad']['breusch_pagan']['p-valor LM']
p_white = diag_trans['homocedasticidad']['white']['p-valor LM']
# Añadir resultados
resultados.append({
'Transformación': metodo,
'R² Ajustado': modelo_trans.rsquared_adj,
'AIC': modelo_trans.aic,
'BIC': modelo_trans.bic,
'p-valor BP': p_bp,
'p-valor White': p_white,
'DW': diag_trans['independencia']['durbin_watson'],
'Lambda': lambda_y if metodo == 'box-cox' else None
})
except Exception as e:
print(f"Error en transformación {metodo}: {e}")
# Añadir resultados del modelo original
p_bp_orig = resultados_original['homocedasticidad']['breusch_pagan']['p-valor LM']
p_white_orig = resultados_original['homocedasticidad']['white']['p-valor LM']
resultados.append({
'Transformación': 'Original',
'R² Ajustado': modelo_original.rsquared_adj,
'AIC': modelo_original.aic,
'BIC': modelo_original.bic,
'p-valor BP': p_bp_orig,
'p-valor White': p_white_orig,
'DW': resultados_original['independencia']['durbin_watson'],
'Lambda': None
})
# Crear DataFrame con resultados
df_resultados = pd.DataFrame(resultados)
# Ordenar por AIC (menor es mejor)
df_resultados = df_resultados.sort_values('AIC')
return df_resultados
4.2 Transformaciones para Corregir No Normalidad#
La no normalidad de los residuos puede afectar la validez de las inferencias. Las transformaciones para corregir este problema son similares a las de heterocedasticidad, pero el enfoque es diferente:
def transformaciones_normalidad(X, y, metodo='yeo-johnson'):
"""
Aplica transformaciones para corregir no normalidad de los residuos.
Parámetros:
-----------
X : DataFrame
Variables predictoras
y : Series
Variable respuesta
metodo : str
Tipo de transformación: 'box-cox', 'yeo-johnson', 'quantile'
Retorna:
--------
tuple : (X_transformado, y_transformado, parametros)
"""
# Hacemos copias para no modificar los originales
X_trans = X.copy()
# Para transformar y
min_y = min(y)
if min_y <= 0 and metodo == 'box-cox':
y_adj = y - min_y + 0.01 # Ajustamos para asegurar positividad
else:
y_adj = y.copy()
if metodo == 'box-cox':
# Transformación Box-Cox
y_trans, lambda_y = stats.boxcox(y_adj)
parametros = {'lambda': lambda_y}
elif metodo == 'yeo-johnson':
# Transformación Yeo-Johnson (similar a Box-Cox pero funciona con valores negativos)
from sklearn.preprocessing import PowerTransformer
pt = PowerTransformer(method='yeo-johnson')
y_array = np.array(y_adj).reshape(-1, 1)
y_trans = pt.fit_transform(y_array).flatten()
parametros = {'lambda': pt.lambdas_[0]}
elif metodo == 'quantile':
# Transformación de quantiles a normal
from sklearn.preprocessing import QuantileTransformer
qt = QuantileTransformer(output_distribution='normal')
y_array = np.array(y_adj).reshape(-1, 1)
y_trans = qt.fit_transform(y_array).flatten()
parametros = {'distribution': 'normal'}
else:
raise ValueError("Método no reconocido. Use 'box-cox', 'yeo-johnson', o 'quantile'")
return X_trans, pd.Series(y_trans, index=y.index, name=y.name), parametros
def comparar_modelos_normalidad(X, y):
"""
Compara diferentes transformaciones para corregir no normalidad.
Parámetros:
-----------
X : DataFrame
Variables predictoras
y : Series
Variable respuesta
Retorna:
--------
DataFrame : Resultados comparativos
"""
# Modelo original
X_const = sm.add_constant(X)
modelo_original = sm.OLS(y, X_const).fit()
resultados_original = diagnostico_completo(modelo_original, X, y)
# Modelos con transformaciones
transformaciones = ['box-cox', 'yeo-johnson', 'quantile']
resultados = []
for metodo in transformaciones:
try:
X_trans, y_trans, parametros = transformaciones_normalidad(X, y, metodo)
X_trans_const = sm.add_constant(X_trans)
modelo_trans = sm.OLS(y_trans, X_trans_const).fit()
# Diagnóstico del modelo transformado
diag_trans = diagnostico_completo(modelo_trans, X_trans, y_trans)
# Extraer p-valores de las pruebas de normalidad
p_jb = diag_trans['normalidad']['jarque_bera']['p_valor']
p_sw = diag_trans['normalidad']['shapiro']['p_valor']
# Añadir resultados
resultados.append({
'Transformación': metodo,
'R² Ajustado': modelo_trans.rsquared_adj,
'AIC': modelo_trans.aic,
'BIC': modelo_trans.bic,
'p-valor JB': p_jb,
'p-valor SW': p_sw,
'Asimetría': diag_trans['normalidad']['asimetria'],
'Curtosis': diag_trans['normalidad']['curtosis'],
'Parámetros': str(parametros)
})
except Exception as e:
print(f"Error en transformación {metodo}: {e}")
# Añadir resultados del modelo original
p_jb_orig = resultados_original['normalidad']['jarque_bera']['p_valor']
p_sw_orig = resultados_original['normalidad']['shapiro']['p_valor']
resultados.append({
'Transformación': 'Original',
'R² Ajustado': modelo_original.rsquared_adj,
'AIC': modelo_original.aic,
'BIC': modelo_original.bic,
'p-valor JB': p_jb_orig,
'p-valor SW': p_sw_orig,
'Asimetría': resultados_original['normalidad']['asimetria'],
'Curtosis': resultados_original['normalidad']['curtosis'],
'Parámetros': None
})
# Crear DataFrame con resultados
df_resultados = pd.DataFrame(resultados)
# Ordenar por AIC (menor es mejor)
df_resultados = df_resultados.sort_values('AIC')
return df_resultados
4.3 Transformaciones para Corregir No Linealidad#
Cuando la relación entre predictores y respuesta no es lineal, podemos aplicar transformaciones a las variables predictoras o incluir términos no lineales:
def transformaciones_linealidad(X, y, metodo='polinomico', grado=2):
"""
Aplica transformaciones para corregir no linealidad.
Parámetros:
-----------
X : DataFrame
Variables predictoras
y : Series
Variable respuesta
metodo : str
Tipo de transformación: 'polinomico', 'log_predictores', 'splines'
grado : int
Grado del polinomio o número de nodos para splines
Retorna:
--------
tuple : (X_transformado, y_transformado, detalles)
"""
# Hacemos copias para no modificar los originales
y_trans = y.copy()
if metodo == 'polinomico':
# Añadir términos polinómicos
X_trans = X.copy()
for col in X.columns:
for g in range(2, grado + 1):
X_trans[f'{col}^{g}'] = X[col] ** g
detalles = {'grado': grado}
elif metodo == 'log_predictores':
# Aplicar log a los predictores
X_trans = X.copy()
for col in X.columns:
if (X[col] > 0).all(): # Solo aplicar log a columnas positivas
X_trans[f'log({col})'] = np.log(X[col])
else:
print(f"Advertencia: No se puede aplicar log a {col} - contiene valores no positivos")
detalles = {'transformaciones': 'logarítmicas'}
elif metodo == 'splines':
# Aplicar splines naturales
try:
from patsy import dmatrix
X_splines = {}
for col in X.columns:
# Crear una matriz de diseño con splines naturales
spline_matrix = dmatrix(f"bs({col}, df={grado})", data=X, return_type='dataframe')
# Renombrar columnas
new_cols = {col_name: f"{col}_spline_{i}" for i, col_name in enumerate(spline_matrix.columns)}
spline_matrix = spline_matrix.rename(columns=new_cols)
X_splines.update({new_name: spline_matrix[new_name] for new_name in new_cols.values()})
X_trans = pd.DataFrame(X_splines)
detalles = {'grados_libertad': grado}
except ImportError:
print("La biblioteca 'patsy' no está instalada. Utilizando polinomios en su lugar.")
# Caer en el caso polinómico si no está disponible patsy
X_trans = X.copy()
for col in X.columns:
for g in range(2, grado + 1):
X_trans[f'{col}^{g}'] = X[col] ** g
detalles = {'grado': grado}
else:
raise ValueError("Método no reconocido. Use 'polinomico', 'log_predictores', o 'splines'")
return X_trans, y_trans, detalles
def comparar_modelos_linealidad(X, y):
"""
Compara diferentes transformaciones para corregir no linealidad.
Parámetros:
-----------
X : DataFrame
Variables predictoras
y : Series
Variable respuesta
Retorna:
--------
DataFrame : Resultados comparativos
"""
# Modelo original
X_const = sm.add_constant(X)
modelo_original = sm.OLS(y, X_const).fit()
resultados_original = diagnostico_completo(modelo_original, X, y)
# Modelos con transformaciones
transformaciones = [
('polinomico', 2),
('polinomico', 3),
('log_predictores', None),
('splines', 3)
]
resultados = []
for metodo, grado in transformaciones:
try:
if grado is None:
X_trans, y_trans, detalles = transformaciones_linealidad(X, y, metodo)
else:
X_trans, y_trans, detalles = transformaciones_linealidad(X, y, metodo, grado)
X_trans_const = sm.add_constant(X_trans)
modelo_trans = sm.OLS(y_trans, X_trans_const).fit()
# Diagnóstico del modelo transformado
diag_trans = diagnostico_completo(modelo_trans, X_trans, y_trans)
# Extraer p-valores de la prueba de linealidad
p_reset = diag_trans['linealidad']['ramsey_reset']['p_valor']
# Añadir resultados
nombre_transf = f"{metodo}"
if grado is not None:
nombre_transf += f" (grado {grado})"
resultados.append({
'Transformación': nombre_transf,
'R² Ajustado': modelo_trans.rsquared_adj,
'AIC': modelo_trans.aic,
'BIC': modelo_trans.bic,
'p-valor RESET': p_reset,
'Num. Predictores': len(X_trans.columns),
'Detalles': str(detalles)
})
except Exception as e:
print(f"Error en transformación {metodo}: {e}")
# Añadir resultados del modelo original
p_reset_orig = resultados_original['linealidad']['ramsey_reset']['p_valor']
resultados.append({
'Transformación': 'Original',
'R² Ajustado': modelo_original.rsquared_adj,
'AIC': modelo_original.aic,
'BIC': modelo_original.bic,
'p-valor RESET': p_reset_orig,
'Num. Predictores': len(X.columns),
'Detalles': None
})
# Crear DataFrame con resultados
df_resultados = pd.DataFrame(resultados)
# Ordenar por AIC (menor es mejor)
df_resultados = df_resultados.sort_values('AIC')
return df_resultados
4.4 Corrección para Autocorrelación#
La autocorrelación en los residuos es común en datos de series temporales. Podemos corregirla mediante:
Método Cochrane-Orcutt: Ajusta el modelo iterativamente para eliminar la autocorrelación
Método Prais-Winsten: Similar al anterior pero ajusta el primer valor de la serie
def corregir_autocorrelacion(X, y, metodo='cochrane-orcutt', max_iter=20):
"""
Corrige la autocorrelación en los residuos.
Parámetros:
-----------
X : DataFrame
Variables predictoras
y : Series
Variable respuesta
metodo : str
Método de corrección: 'cochrane-orcutt', 'prais-winsten'
max_iter : int
Número máximo de iteraciones
Retorna:
--------
tuple : (modelo_corregido, rho_estimado)
"""
# Ajustar modelo inicial
X_const = sm.add_constant(X)
modelo_inicial = sm.OLS(y, X_const).fit()
# Obtener residuos
residuos = modelo_inicial.resid
if metodo == 'cochrane-orcutt':
# Implementación del método Cochrane-Orcutt
rho = 0.0
convergencia = False
for _ in range(max_iter):
# Estimar rho utilizando los residuos actuales
rho_nuevo = np.corrcoef(residuos[:-1], residuos[1:])[0, 1]
# Verificar convergencia
if abs(rho_nuevo - rho) < 1e-6:
convergencia = True
break
rho = rho_nuevo
# Transformar variables
y_trans = y.iloc[1:] - rho * y.iloc[:-1]
X_trans = X_const.iloc[1:] - rho * X_const.iloc[:-1]
# Reajustar el modelo con variables transformadas
modelo_trans = sm.OLS(y_trans, X_trans).fit()
# Actualizar residuos
residuos = modelo_trans.resid
# Mensaje según convergencia
if convergencia:
print(f"Convergencia alcanzada en {_ + 1} iteraciones. Rho estimado: {rho:.4f}")
else:
print(f"No se alcanzó convergencia en {max_iter} iteraciones. Último rho: {rho:.4f}")
# Ajustar modelo final con el rho estimado
y_final = y.iloc[1:] - rho * y.iloc[:-1]
X_final = X_const.iloc[1:] - rho * X_const.iloc[:-1]
modelo_final = sm.OLS(y_final, X_final).fit()
elif metodo == 'prais-winsten':
# Método Prais-Winsten (similar a Cochrane-Orcutt pero conserva la primera observación)
from statsmodels.regression.linear_model import GLSAR
# GLSAR implementa Prais-Winsten por defecto
modelo_final = GLSAR(y, X, rho=1).iterative_fit(maxiter=max_iter)
rho = modelo_final.rho
else:
raise ValueError("Método no reconocido. Use 'cochrane-orcutt' o 'prais-winsten'")
return modelo_final, rho
5. Ejemplos Prácticos#
En esta sección, aplicaremos las herramientas desarrolladas a conjuntos de datos con diferentes problemas de violaciones de supuestos. Para cada caso, mostraremos cómo identificar el problema mediante diagnósticos y cómo aplicar transformaciones correctivas.
5.1 Ejemplo 1: Detección y Corrección de Heterocedasticidad#
# Generamos datos con heterocedasticidad
np.random.seed(123)
X_heter, y_heter = generar_datos_prueba(n=200, p=3, problema='heterocedasticidad', intensidad=1.5)
# Ajustamos modelo a los datos originales
X_heter_const = sm.add_constant(X_heter)
modelo_heter = sm.OLS(y_heter, X_heter_const).fit()
# Diagnóstico completo
resultados_heter = diagnostico_completo(modelo_heter, X_heter, y_heter)
# Imprimir resultados del diagnóstico
imprimir_resultados_diagnostico(resultados_heter)
# Visualizaciones diagnósticas
plot_diagnostico_residuos(resultados_heter['residuos'])
plot_residuos_vs_predictores(resultados_heter['residuos'], X_heter)
# Aplicar transformaciones y comparar resultados
comparativa_heter = comparar_modelos_heterocedasticidad(X_heter, y_heter)
print("\nComparación de modelos para corregir heterocedasticidad:")
print(comparativa_heter)
# Seleccionamos la mejor transformación (por ejemplo, Box-Cox)
X_trans, y_trans, lambda_y = transformaciones_heterocedasticidad(X_heter, y_heter, metodo='box-cox')
# Ajustamos nuevo modelo con datos transformados
X_trans_const = sm.add_constant(X_trans)
modelo_trans = sm.OLS(y_trans, X_trans_const).fit()
# Diagnóstico del modelo transformado
resultados_trans = diagnostico_completo(modelo_trans, X_trans, y_trans)
# Visualizaciones comparativas
print("\nDiagnóstico del modelo transformado:")
imprimir_resultados_diagnostico(resultados_trans)
plot_diagnostico_residuos(resultados_trans['residuos'])
# Interpretación de coeficientes con variable transformada
print("\nCoeficientes del modelo original:")
print(modelo_heter.summary2().tables[1][['Coef.', 'Std.Err.', 'P>|t|']])
print("\nCoeficientes del modelo transformado:")
print(modelo_trans.summary2().tables[1][['Coef.', 'Std.Err.', 'P>|t|']])
print("\nNota sobre interpretación de coeficientes transformados:")
if lambda_y == 0: # Transformación logarítmica
print("- Con transformación logarítmica, los coeficientes representan cambios porcentuales aproximados")
print("- Un aumento de 1 unidad en Xj corresponde a un cambio de (100 * βj)% en la variable respuesta")
elif lambda_y == 0.5: # Raíz cuadrada
print("- Con transformación de raíz cuadrada, la interpretación no es directa")
print("- Para interpretar el efecto marginal: dY/dXj = βj / (2 * √Y)")
else:
print(f"- Con transformación Box-Cox (λ = {lambda_y:.4f}), la interpretación requiere cálculos adicionales")
print("- El efecto marginal es: dY/dXj = βj * Y^(1-λ)")
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[10], line 10
7 modelo_heter = sm.OLS(y_heter, X_heter_const).fit()
9 # Diagnóstico completo
---> 10 resultados_heter = diagnostico_completo(modelo_heter, X_heter, y_heter)
12 # Imprimir resultados del diagnóstico
13 imprimir_resultados_diagnostico(resultados_heter)
Cell In[5], line 231, in diagnostico_completo(modelo, X, y)
229 # Pruebas de supuestos
230 resultados_normalidad = test_normalidad(resultados_residuos['estudentizados_internos'])
--> 231 resultados_homocedasticidad = test_homocedasticidad(modelo, X)
232 resultados_independencia = test_independencia(resultados_residuos['ordinarios'])
233 resultados_linealidad = test_linealidad(modelo, X, y)
Cell In[5], line 66, in test_homocedasticidad(modelo, X)
51 """
52 Realiza pruebas de homocedasticidad sobre los residuos.
53
(...)
63 dict : Diccionario con resultados de las pruebas
64 """
65 # Test de Breusch-Pagan
---> 66 bp_test = het_breuschpagan(modelo.resid, X)
67 labels_bp = ['Estadístico LM', 'p-valor LM', 'Estadístico F', 'p-valor F']
69 # Test de White
70 # Necesitamos generar términos cuadráticos e interacciones
File ~\miniconda3\envs\USTA\Lib\site-packages\statsmodels\stats\diagnostic.py:801, in het_breuschpagan(resid, exog_het, robust)
733 r"""
734 Breusch-Pagan Lagrange Multiplier test for heteroscedasticity
735
(...)
798 heteroskedasticity". Journal of Econometrics 17 (1): 107–112.
799 """
800 x = array_like(exog_het, "exog_het", ndim=2)
--> 801 _check_het_test(x, "The Breusch-Pagan")
802 y = array_like(resid, "resid", ndim=1) ** 2
803 if not robust:
File ~\miniconda3\envs\USTA\Lib\site-packages\statsmodels\stats\diagnostic.py:726, in _check_het_test(x, test_name)
721 x_max = x.max(axis=0)
722 if (
723 not np.any(((x_max - x.min(axis=0)) == 0) & (x_max != 0))
724 or x.shape[1] < 2
725 ):
--> 726 raise ValueError(
727 f"{test_name} test requires exog to have at least "
728 "two columns where one is a constant."
729 )
ValueError: The Breusch-Pagan test requires exog to have at least two columns where one is a constant.
5.2 Ejemplo 2: Detección y Corrección de No Normalidad#
# Generamos datos con errores no normales
np.random.seed(456)
X_no_norm, y_no_norm = generar_datos_prueba(n=200, p=3, problema='no_normalidad', intensidad=1.2)
# Ajustamos modelo a los datos originales
X_no_norm_const = sm.add_constant(X_no_norm)
modelo_no_norm = sm.OLS(y_no_norm, X_no_norm_const).fit()
# Diagnóstico completo
resultados_no_norm = diagnostico_completo(modelo_no_norm, X_no_norm, y_no_norm)
# Imprimir resultados del diagnóstico
imprimir_resultados_diagnostico(resultados_no_norm)
# Visualizaciones diagnósticas
plot_diagnostico_residuos(resultados_no_norm['residuos'])
# Histograma de residuos
plt.figure(figsize=(10, 6))
sns.histplot(resultados_no_norm['residuos']['estudentizados_internos'], kde=True, bins=30)
plt.title('Histograma de Residuos Estudentizados')
plt.xlabel('Residuos Estudentizados')
plt.axvline(x=0, color='r', linestyle='--')
plt.show()
# Aplicar transformaciones y comparar resultados
comparativa_norm = comparar_modelos_normalidad(X_no_norm, y_no_norm)
print("\nComparación de modelos para corregir no normalidad:")
print(comparativa_norm)
# Seleccionamos la mejor transformación
X_trans, y_trans, parametros = transformaciones_normalidad(X_no_norm, y_no_norm, metodo='yeo-johnson')
# Ajustamos nuevo modelo con datos transformados
X_trans_const = sm.add_constant(X_trans)
modelo_trans = sm.OLS(y_trans, X_trans_const).fit()
# Diagnóstico del modelo transformado
resultados_trans = diagnostico_completo(modelo_trans, X_trans, y_trans)
# Visualizaciones comparativas
print("\nDiagnóstico del modelo transformado:")
imprimir_resultados_diagnostico(resultados_trans)
plot_diagnostico_residuos(resultados_trans['residuos'])
# Histograma de residuos transformados
plt.figure(figsize=(10, 6))
sns.histplot(resultados_trans['residuos']['estudentizados_internos'], kde=True, bins=30)
plt.title('Histograma de Residuos Estudentizados (Modelo Transformado)')
plt.xlabel('Residuos Estudentizados')
plt.axvline(x=0, color='r', linestyle='--')
plt.show()
5.3 Ejemplo 3: Detección y Corrección de No Linealidad#
# Generamos datos con relación no lineal
np.random.seed(789)
X_no_lin, y_no_lin = generar_datos_prueba(n=200, p=3, problema='no_linealidad', intensidad=1.0)
# Ajustamos modelo a los datos originales
X_no_lin_const = sm.add_constant(X_no_lin)
modelo_no_lin = sm.OLS(y_no_lin, X_no_lin_const).fit()
# Diagnóstico completo
resultados_no_lin = diagnostico_completo(modelo_no_lin, X_no_lin, y_no_lin)
# Imprimir resultados del diagnóstico
imprimir_resultados_diagnostico(resultados_no_lin)
# Visualizaciones diagnósticas
plot_diagnostico_residuos(resultados_no_lin['residuos'])
plot_residuos_vs_predictores(resultados_no_lin['residuos'], X_no_lin)
# Aplicar transformaciones y comparar resultados
comparativa_lin = comparar_modelos_linealidad(X_no_lin, y_no_lin)
print("\nComparación de modelos para corregir no linealidad:")
print(comparativa_lin)
# Seleccionamos la mejor transformación (polinomio de grado 2)
X_trans, y_trans, detalles = transformaciones_linealidad(X_no_lin, y_no_lin, metodo='polinomico', grado=2)
# Ajustamos nuevo modelo con datos transformados
X_trans_const = sm.add_constant(X_trans)
modelo_trans = sm.OLS(y_trans, X_trans_const).fit()
# Diagnóstico del modelo transformado
resultados_trans = diagnostico_completo(modelo_trans, X_trans, y_trans)
# Visualizaciones comparativas
print("\nDiagnóstico del modelo transformado:")
imprimir_resultados_diagnostico(resultados_trans)
plot_diagnostico_residuos(resultados_trans['residuos'])
# Comparación de modelos
print("\nResumen del modelo original:")
print(modelo_no_lin.summary().tables[0])
print("\nResumen del modelo con términos cuadráticos:")
print(modelo_trans.summary().tables[0])
5.4 Ejemplo 4: Detección y Corrección de Autocorrelación#
# Generamos datos con autocorrelación
np.random.seed(101112)
X_auto, y_auto = generar_datos_prueba(n=200, p=3, problema='autocorrelacion', intensidad=1.0)
# Ajustamos modelo a los datos originales
X_auto_const = sm.add_constant(X_auto)
modelo_auto = sm.OLS(y_auto, X_auto_const).fit()
# Diagnóstico completo
resultados_auto = diagnostico_completo(modelo_auto, X_auto, y_auto)
# Imprimir resultados del diagnóstico
imprimir_resultados_diagnostico(resultados_auto)
# Visualizaciones diagnósticas
plot_diagnostico_residuos(resultados_auto['residuos'])
plot_secuencia_temporal(resultados_auto['residuos'])
plot_autocorrelacion(resultados_auto['residuos'])
# Aplicar corrección y comparar resultados
print("\nCorrección mediante Cochrane-Orcutt:")
modelo_co, rho_co = corregir_autocorrelacion(X_auto, y_auto, metodo='cochrane-orcutt')
print(f"Rho estimado: {rho_co:.4f}")
print("\nCorrección mediante Prais-Winsten:")
modelo_pw, rho_pw = corregir_autocorrelacion(X_auto, y_auto, metodo='prais-winsten')
print(f"Rho estimado: {rho_pw:.4f}")
# Comparamos los errores estándar y significancia
print("\nComparación de coeficientes y errores estándar:")
print("\nModelo Original:")
print(modelo_auto.summary2().tables[1][['Coef.', 'Std.Err.', 'P>|t|']])
print("\nModelo Cochrane-Orcutt:")
print(modelo_co.summary2().tables[1][['Coef.', 'Std.Err.', 'P>|t|']])
print("\nModelo Prais-Winsten:")
print(modelo_pw.summary2().tables[1][['Coef.', 'Std.Err.', 'P>|t|']])
# Verificar si se corrigió la autocorrelación
# Calculamos estadístico Durbin-Watson para el modelo corregido (Prais-Winsten)
dw_orig = resultados_auto['independencia']['durbin_watson']
dw_corr = durbin_watson(modelo_pw.resid)
print(f"\nEstadístico Durbin-Watson (Original): {dw_orig:.4f}")
print(f"Estadístico Durbin-Watson (Corregido): {dw_corr:.4f}")
if abs(dw_corr - 2.0) < abs(dw_orig - 2.0):
print("La corrección mejoró la autocorrelación.")
else:
print("La corrección no mejoró significativamente la autocorrelación.")
6. Ejercicios#
A continuación se presentan ejercicios para practicar y profundizar en el análisis de residuos. Se recomienda realizarlos en orden, ya que cada uno construye sobre los conocimientos adquiridos en los anteriores.
Ejercicio 1: Análisis de Residuos en Dataset Real#
Utiliza el dataset de Boston Housing disponible en scikit-learn para:
Ajustar un modelo lineal utilizando todas las variables disponibles.
Realizar un diagnóstico completo de los residuos.
Identificar qué supuestos se violan.
Aplicar al menos dos transformaciones correctivas para el problema principal detectado.
Comparar los modelos original y corregido mediante criterios estadísticos (AIC, BIC, R² ajustado).
Código inicial:
from sklearn.datasets import load_boston
# Cargar datos
boston = load_boston()
X_boston = pd.DataFrame(boston.data, columns=boston.feature_names)
y_boston = pd.Series(boston.target, name='MEDV')
# TODO: Realizar el análisis de residuos
Ejercicio 2: Identificación de Valores Atípicos#
Utilizando el mismo dataset de Boston:
Calcula los residuos estudentizados externos.
Identifica observaciones con valores absolutos de residuos estudentizados externos mayores a 3.
Calcula la distancia de Cook para cada observación.
Identifica observaciones con distancia de Cook mayor a 4/n.
Compara el modelo con todas las observaciones y el modelo sin las observaciones influyentes.
Sugerencia: Puedes utilizar las funciones de influencia de statsmodels:
# Ejemplo
modelo_boston = sm.OLS(y_boston, sm.add_constant(X_boston)).fit()
influence = modelo_boston.get_influence()
student_resid = influence.resid_studentized_external
cook_dist = influence.cooks_distance[0]
Ejercicio 3: Transformaciones para Problemas Múltiples#
A menudo, los datos presentan múltiples violaciones de supuestos simultáneamente. Por ejemplo, puede haber heterocedasticidad y no normalidad al mismo tiempo.
Genera un conjunto de datos sintético con dos problemas simultáneos (por ejemplo, heterocedasticidad y no normalidad).
Implementa una estrategia para corregir ambos problemas.
Compara el modelo original con el modelo corregido.
Discute las ventajas y desventajas de corregir múltiples problemas simultáneamente versus abordarlos uno por uno.
Código inicial:
# Generar datos con múltiples problemas
np.random.seed(9876)
X_multi = pd.DataFrame(np.random.normal(0, 1, size=(200, 3)), columns=['X1', 'X2', 'X3'])
beta = np.array([0.5, 1.0, -0.7])
# Error con heterocedasticidad y no normalidad
epsilon_base = np.random.normal(0, 1, size=200)
factor_varianza = (1 + 1.2 * np.abs(X_multi['X1']))
epsilon_chi2 = np.random.chisquare(3, size=200) - 3
epsilon = factor_varianza * (0.7 * epsilon_chi2 + 0.3 * epsilon_base)
y_multi = X_multi @ beta + epsilon
y_multi = pd.Series(y_multi, name='y')
# TODO: Implementar estrategia para corregir ambos problemas
Ejercicio 4: Cross-Validation para Selección de Transformaciones#
Implementa un enfoque de validación cruzada para seleccionar la mejor transformación correctiva:
Divide el dataset en conjuntos de entrenamiento y prueba.
Ajusta diferentes transformaciones correctivas en el conjunto de entrenamiento.
Evalúa cada modelo transformado en el conjunto de prueba.
Selecciona la transformación que produce el menor error cuadrático medio en el conjunto de prueba.
Código inicial:
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
# Usar dataset de Boston
X_train, X_test, y_train, y_test = train_test_split(X_boston, y_boston, test_size=0.3, random_state=42)
# Transformaciones a evaluar
transformaciones_y = ['log', 'sqrt', 'reciproco', 'box-cox']
# TODO: Implementar validación cruzada para selección de transformaciones
Ejercicio 5: Investigación de Supuestos en la Literatura#
Investiga y presenta un breve resumen sobre uno de los siguientes temas:
¿Qué tan robustos son los modelos lineales ante la violación del supuesto de normalidad?
¿En qué situaciones prácticas es más crítica la violación del supuesto de homocedasticidad?
¿Qué alternativas existen al método de mínimos cuadrados ordinarios cuando hay heterocedasticidad?
¿Cómo afecta la multicolinealidad al análisis de residuos?
Incluye al menos tres referencias académicas y presenta tus conclusiones en forma de ensayo breve (máximo 500 palabras).