Introducción

Para este ejemplo simularemos dos conjuntos de pares de observaciones \((x, y)\) tales que uno de ellos cumpla con el supuesto de varianza constante y el otro no. En ambos casos la relación entre la media de \(X\) con \(X\) es la misma, la diferencia radica en la función de varianza, \(V(\epsilon_i) = \sigma^2 x_i^2\).

n <- 250  # tamaño de muestra

set.seed(2201)
x <- rnorm(76, 20, 6)                            # Se generan las X
y1 <- 28 + 5 * x + rnorm(76, 0, 25)           # Se generan las Y con errores homocedásticos
y2 <- 28 + 5 * x + x * rnorm(76, 0, 5)  # Se generan las Y con errores heterocedásticos


Graficamos ambos conjuntos de datos, en ambos casos se observa que la relación entre \(X\) y \(Y\) es lineal. El la primera gráfica vemos como la dispersión de los errores alrededor de la recta de regresión, es constante en todo el intervalo, mientras que en la segunda gráfica vemos como la dispersión de los errores es mayor a medida que \(X\) crece.

par(mfrow = c(1, 2))
plot(x, y1, pch = 16, col = 'darkcyan',
xlab = 'V. explicativa', ylab = 'V. respuesta', main = 'Errores homocedásticos')
lines(loess.smooth(x, y1), col = 'red3', lwd = 3, lty = 2)
plot(x, y2, pch = 16, col = 'darkcyan',
xlab = 'V. explicativa', ylab = 'V. respuesta', main = 'Errores heterocedásticos')
lines(loess.smooth(x, y2), col = 'red3', lwd = 3, lty = 2)


Efectos de la heterocedasticidad

Ajustamos un modelo RLS a cada conjunto de datos y observamos los resultados.

modelo_1 <- lm(y1 ~ x) # Modelo con errores homocedásticos
summary(modelo_1)
## 
## Call:
## lm(formula = y1 ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -52.330 -15.390   2.478  17.118  50.004 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   19.560     11.073   1.766   0.0815 .  
## x              5.236      0.529   9.898 3.42e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.99 on 74 degrees of freedom
## Multiple R-squared:  0.5697, Adjusted R-squared:  0.5639 
## F-statistic: 97.97 on 1 and 74 DF,  p-value: 3.42e-15
modelo_2 <- lm(y2 ~ x) # Modelo con errores heterocedásticos
summary(modelo_2)
## 
## Call:
## lm(formula = y2 ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -289.161  -47.505    3.595   49.695  259.228 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   76.130     43.009    1.77   0.0808 .
## x              2.959      2.054    1.44   0.1540  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 93.16 on 74 degrees of freedom
## Multiple R-squared:  0.02726,    Adjusted R-squared:  0.01412 
## F-statistic: 2.074 on 1 and 74 DF,  p-value: 0.154


Observamos varias diferencias entre los resultados. La significnacia de los parámetros no es la misma, ni tampoco son iguales los coeficientes \(R^2\).


Detección de heterocedasticidad

Métodos gráficos

Para detectar heterocedasticidad podemos utilizar gráficas de dispersión de los valores ajustados contra los residuos estandarizados. Utilizamos la función stdres del paquete MASS para obtener los residuos estandarizados

ajustados_1 <- fitted(modelo_1)
ajustados_2 <- fitted(modelo_2)

library(MASS)
residuos_1 <- stdres(modelo_1)
residuos_2 <- stdres(modelo_2)



Es díficil identificar patrones en las gráficas anteriores, por ello gráficamos ahora los valores absolutos de los residuos estandarizados.


par(mfrow = c(1, 2))
plot(ajustados_1, abs(residuos_1), pch = 16, col = 'darkcyan',
xlab = 'V. explicativa', ylab = 'abs(Residuos)', main = 'Errores homocedásticos')
lines(loess.smooth(ajustados_1, abs(residuos_1)), col = 'red3', lwd = 3, lty = 3)
plot(ajustados_2, abs(residuos_2), pch = 16, col = 'darkcyan',
xlab = 'V. explicativa', ylab = 'abs(Residuos)', main = 'Errores heterocedásticos')
lines(loess.smooth(ajustados_2, abs(residuos_2)), col = 'red3', lwd = 3, lty = 3)


En la gráfica de la derecha parece no haber una tendencia en los residuos (creciente o decreciente), por lo que se podría consideraar que los errores tienen varianza constantes. En la gráfica de la derecha se identifica una tendencia decreciente en los residuos, esto podría indicar presencia de heterocedasticidad.


Pruebas estadísticas

Para confirmar las conclusiones obtenidas con los métodos gráficos, utilizamos las pruebas de Breush-Pagan y de White.


Breusch-Pagan

Primero la obtenemos con la función bptest del paquete lmtest. La hipótesis nula es establece homocedasticidad en los errores.

library(lmtest)
bptest(modelo_1)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_1
## BP = 0.051916, df = 1, p-value = 0.8198
bptest(modelo_2)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_2
## BP = 13.63, df = 1, p-value = 0.0002226

Vemos como no se rechaza la hipótesis de homocedasticidad en el conjunto de datos generado con errores constantes, mientras que sí se rechaza la hipótesis en el conjunto de datos generado con errores no constantes. Con esto confirmamos las conclusiones obtenidas con métodos gráficos.


Ahora hacemos una aproximación a la prueba de Breush-Pagan, utilizando el modelo

\[ u_i = \alpha_0 + \alpha_1x_i + \eta_i \] donde \(u_i = t^2_i\) y \(t_i\) son los residuos studentizados, \(i = 1, \ldots, n\). Rechazamos que la varianza es constante con la prueba de \(F\) del ANOVA para el modelo anterior. Para calcular los residuos studentizados utilizamos la función studres del paquete MASS.


residuos_stu_1 <- studres(modelo_1)
residuos_stu_2 <- studres(modelo_2)
bp_modelo1 <- lm(I(residuos_stu_1^2) ~ x)
summary(bp_modelo1)
## 
## Call:
## lm(formula = I(residuos_stu_1^2) ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0760 -0.8752 -0.5322  0.5335  4.2792 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.845746   0.568038   1.489    0.141
## x           0.008729   0.027135   0.322    0.749
## 
## Residual standard error: 1.23 on 74 degrees of freedom
## Multiple R-squared:  0.001397,   Adjusted R-squared:  -0.0121 
## F-statistic: 0.1035 on 1 and 74 DF,  p-value: 0.7486
bp_modelo2 <- lm(I(residuos_stu_2^2) ~ x)
summary(bp_modelo2)
## 
## Call:
## lm(formula = I(residuos_stu_2^2) ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2730 -0.9291 -0.4914  0.3188  9.1897 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.26348    0.85598  -2.644 0.009990 ** 
## x            0.16383    0.04089   4.007 0.000145 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.854 on 74 degrees of freedom
## Multiple R-squared:  0.1783, Adjusted R-squared:  0.1672 
## F-statistic: 16.05 on 1 and 74 DF,  p-value: 0.0001452


Las conclusiones obtenidas de las pruebas anteriores coinciden con las obtenidas a partir de la función bptest, sin embargo las significancias y estadísticos de prueba tienen valores diferentes. Se debe tener especial atención en las aproximaciones asintóticas cuando el tamaño de muestra es pequeño, en nuestro ejemplo \(n = 76\).


Correcciones por heterocedasticidad

Se puede corregir el problema de varianza no constante en los errores aplicando transformaciones a los datos. Sin embargo, si no se tiene cuidado, se puede afectar la linealidad o la normalidad. En este ejemplo, dado que generamos los errores con una distribución normal y sabemos sí se cumple el supuesto de linealidad, la mejor opción no es transformar los datos. En su lugar, consideraremos dos alternativas, ajustar por mínimos cuadrados ponderados y bootstrap. En el primer caso se trata con una extensión del modelo y en el segundo de una aproximación por simulación.

Mínimos cuadrados ponderados

Ahora ajustamos el modelo por minimos cuadrados ponderados. Dado que sabemos cómo se generaron los errores del modelo con heterocedasticidad, utilizamos como ponderadores a \(w_i = 1/x_i^2\), \(i=1,\ldots, n\). En un ejercicio real, el mayor problema es determinar cuáles son los pesos adecuados.

modelo_mcp <- lm(y2 ~ x, weights = 1/x^2)
summary(modelo_mcp)
## 
## Call:
## lm(formula = y2 ~ x, weights = 1/x^2)
## 
## Weighted Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.6306  -2.5880   0.1597   2.3220   9.9395 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   67.032     31.569   2.123   0.0371 *
## x              3.415      1.741   1.961   0.0536 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.221 on 74 degrees of freedom
## Multiple R-squared:  0.04942,    Adjusted R-squared:  0.03658 
## F-statistic: 3.847 on 1 and 74 DF,  p-value: 0.05359


Vemos en los resultados de este nuevo ajuste como los coeficientes estimados no difieren tanto de los obtenidos sin poderar. Sin embargo, la varianza estimada sí se redujo considerablemente. En la siguiente gráfica se representan los datos del conjunto generado con errores heterocedásticos y se agregan las rectas ajustadas con MCO y con MCP, se puede apreciar como las diferencias son mínimas.



Para verificar que el problema se haya resuelto, calculamos y gráficamos los valores ajustados contra los residuos estandarizados (y ponderados). Podemos observar como se ha eliminado la tendencia en la varianza de los errores.

ajustados_mcp <- fitted(modelo_mcp)
residuos_mcp <- stdres(modelo_mcp)



Ahora verificamos las conclusiones obtenidas a partir de las gráficas, para ello realizamos una prueba de Breusch-Pagan. En este caso procedemos sólo con la prueba hecha a mano.


bp_modelo_mcp <- lm(I(residuos_mcp^2) ~ x)
summary(bp_modelo_mcp)
## 
## Call:
## lm(formula = I(residuos_mcp^2) ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.4174 -0.8571 -0.5355  0.3327  5.2829 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.32945    0.66549   0.495    0.622
## x            0.03334    0.03179   1.049    0.298
## 
## Residual standard error: 1.441 on 74 degrees of freedom
## Multiple R-squared:  0.01465,    Adjusted R-squared:  0.001331 
## F-statistic:   1.1 on 1 and 74 DF,  p-value: 0.2977


Dado que el p-value anterior es mayor a 0.05, se corrobora que el problema de varianza no constante se ha corregido.


Inferencias usando bootstrap

Para tener una idea del efecto que tiene la heterocedasticidad en las inferencias sobre los \(\beta\)s, se compararán los intervalos de confianza para \(\beta_0\) y \(\beta_1\) que se obtienen con \(MCO\), con \(MCP\) y con bootstrap, con el conjunto de datos generados con errores heterocedásticos. En los tres casos se considera 90% de confianza.

  • Primero calculamos los intervalos con el modelo ajustado por \(MCO\).
confint(modelo_2, level = 0.9)
##                    5 %       95 %
## (Intercept)  4.4904922 147.769767
## x           -0.4634172   6.380927


  • Después calculamos los intervalos con el modelo ajustado por \(MCP\).
confint(modelo_mcp, level = 0.9)
##                    5 %       95 %
## (Intercept) 14.4479849 119.616299
## x            0.5148873   6.314812
  • Finalmente calculamos los intervalos con bootstrap. Como sabemos que hay problemas de heterocedasticidad, no se procede como cuando los errores no siguen una distribución normal, en este caso debemos seleccionar muestras de pares \((x, y)\) y ajustar un modelo con cada muestra.
beta <- data.frame(beta0 = rep(0, 1E4), beta1 = rep(0, 1E4))
for (i in 1:1E4){
  s <- sample(1:n, n, replace = T)
  xs <- x[s]; ys <- y2[s]
  fit <- lm(ys ~ xs, weights = 1/xs^2)
  beta[i,] <- coef(fit)
}

res <- apply(beta, 2, function(x) quantile(x, probs = c(0.05, 0.95)))
t(res)
##                 5%        95%
## beta0  5.534706179 129.614962
## beta1 -0.008398445   6.831847


Si comparamos los resultados anteriores podemos notar como los intervalos más amplios son los construidos con el ajuste basado en \(MCO\), depués están un poco más estrechos los intervalos generados con bootstrap y los más estrechos son los intervalos basados en \(MCP\). También se debe notar que con el par de intervalos más amplios, no se descarta que \(\beta_1\) sea significativamente diferente de cero, \(\alpha = 0.1\), mientras que con los intervalos basados en \(MCP\) la conclusión es que la variable explicativa sí es significativa.

Comentarios finales