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 \(y\) con \(x\) es la misma, la diferncia radica en la función de varianza.
n <- 250
set.seed(111)
x <- rnorm(250, 70, 20)
# summary(x)
y1 <- 80 + 4 * x + rnorm(250, 0, 40)
y2 <- 80 + 4 * x + x * rnorm(250, 0, 1)
Graficamos ambos conjuntos de datos.
par(mfrow = c(1, 2))
plot(x, y1, pch = 16, col = 'steelblue', main = 'Varianza constante',
xlab = 'V. explicativa', ylab = 'V. respuesta')
plot(x, y2, pch = 16, col = 'steelblue', main = 'Varianza no constante',
xlab = 'V. explicativa', ylab = 'V. respuesta')
Ajustamos un modelo RLS a cada conjunto de datos y observamos los resultados.
fit1 <- lm(y1 ~ x)
summary(fit1)
##
## Call:
## lm(formula = y1 ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -126.28 -26.97 0.85 24.83 119.21
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 78.2490 9.2023 8.503 1.75e-15 ***
## x 3.9964 0.1255 31.848 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 40.14 on 248 degrees of freedom
## Multiple R-squared: 0.8035, Adjusted R-squared: 0.8027
## F-statistic: 1014 on 1 and 248 DF, p-value: < 2.2e-16
fit2 <- lm(y2 ~ x)
summary(fit2)
##
## Call:
## lm(formula = y2 ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -206.930 -46.689 -3.606 42.066 236.796
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 75.7778 16.1998 4.678 4.77e-06 ***
## x 4.0978 0.2209 18.550 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 70.66 on 248 degrees of freedom
## Multiple R-squared: 0.5811, Adjusted R-squared: 0.5795
## F-statistic: 344.1 on 1 and 248 DF, p-value: < 2.2e-16
Vemos como las estimaciones puntuales son similares, no así la estimación de la varianza.
Ahora exploramos los residuos de cada conjunto de datos.
res1 <- residuals(fit1)
res2 <- residuals(fit2)
par(mfrow = c(1, 2))
plot(x, res1, pch = 16, col = 'steelblue', main = 'Varianza constante',
xlab = 'V. explicativa', ylab = 'Residuos')
plot(x, res2, pch = 16, col = 'steelblue', main = 'Varianza no constante',
xlab = 'V. explicativa', ylab = 'Residuos')
En la gráfica de la izquierda no se observan patrones en los residuos, mientras que en la gráfica de la derecha podemos apreciar la figura de megáfono.
Para confirmar nuestras sospechas hacemos la prueba de varianza constante. Primero estimamos la varianza en cada uno de los modelos y después ajustamos los modelos correspondientes.
sh.1 <- sum(res1^2) / (n - 2)
sh.2 <- sum(res2^2) / (n - 2)
u1 <- res1^2/sh.1
u2 <- res2^2/sh.2
fit.hom1 <- lm(u1 ~ x)
summary(fit.hom1)
##
## Call:
## lm(formula = u1 ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.3176 -0.8635 -0.5538 0.1848 8.9947
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.502949 0.354743 1.418 0.158
## x 0.006938 0.004837 1.434 0.153
##
## Residual standard error: 1.547 on 248 degrees of freedom
## Multiple R-squared: 0.008227, Adjusted R-squared: 0.004227
## F-statistic: 2.057 on 1 and 248 DF, p-value: 0.1528
fit.hom2 <- lm(u2 ~ x)
summary(fit.hom2)
##
## Call:
## lm(formula = u2 ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7631 -0.8413 -0.3808 0.2376 9.0523
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.63663 0.34611 -1.839 0.0671 .
## x 0.02311 0.00472 4.896 1.77e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.51 on 248 degrees of freedom
## Multiple R-squared: 0.08812, Adjusted R-squared: 0.08444
## F-statistic: 23.97 on 1 and 248 DF, p-value: 1.768e-06
Podemos ver que en el primer conjunto de datos el modelo no es significativo, lo que nos dice que hay evidencia para sostener que la varianza es constante.
En el otro conjunto de datos, vemos como el coeficiente de \(x\) es distinto de 0 con una significancia alta. Esto significa que hay evidencia para sostener que la varianza no es constante.
Ahora ajustamos el modelo por minimós cuadrados ponderados.
fit3 <- lm(y2 ~ x, weights = 1/x^2)
summary(fit3)
##
## Call:
## lm(formula = y2 ~ x, weights = 1/x^2)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -2.34963 -0.65609 -0.06243 0.64557 2.73609
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 72.76231 3.14440 23.14 <2e-16 ***
## x 4.14886 0.08085 51.31 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9574 on 248 degrees of freedom
## Multiple R-squared: 0.9139, Adjusted R-squared: 0.9136
## F-statistic: 2633 on 1 and 248 DF, p-value: < 2.2e-16
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 bastantes.
plot(x, y2, pch = 16, col = 'steelblue', main = 'Varianza no constante',
xlab = 'V. explicativa', ylab = 'V. respuesta')
abline(fit2, col = 'red2', lty = 2, lwd = 2)
abline(fit3, col = 'orange2', lty = 2, lwd = 2)
Ahora calculamos los residuos ponderados y gráficamos contra la variable explicativa.
res3 <- residuals(fit3)/x
par(mfrow = c(1, 2))
plot(x, res2, pch = 16, col = 'steelblue', main = 'MCO',
xlab = 'V. explicativa', ylab = 'Residuos')
plot(x, res3, pch = 16, col = 'steelblue', main = 'MCP',
xlab = 'V. explicativa', ylab = 'Residuos')
sh.3 <- sum(res3^2)/(n - 2)
u3 <- res3^2 / sh.3
fit.hom3 <- lm(u3 ~ x)
summary(fit.hom3)
##
## Call:
## lm(formula = u3 ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.0767 -0.8863 -0.5287 0.2573 7.1649
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.799199 0.311626 2.565 0.0109 *
## x 0.002735 0.004249 0.644 0.5204
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.359 on 248 degrees of freedom
## Multiple R-squared: 0.001668, Adjusted R-squared: -0.002358
## F-statistic: 0.4143 on 1 and 248 DF, p-value: 0.5204
Podemos ver que se ha corregido el problema de varianza no constante.
Supongamos que no tomamos en cuenta la heterocedasticidad y construimos intervalos de confianza 95% para los parámetros del modelo.
X <- matrix(c(rep(1, n), x), ncol = 2, byrow = F)
vh.b <- sh.2 * solve(t(X)%*%X)
b0_IC <- coef(fit2)[1] + c(-1, 1)*qt(0.95, n - 2)*sqrt(vh.b[1, 1]); b0_IC
## [1] 49.03156 102.52395
b1_IC <- coef(fit2)[2] + c(-1, 1)*qt(0.95, n - 2)*sqrt(vh.b[2, 2]); b1_IC
## [1] 3.73304 4.46248
Ahora calculamos los intervalos de confianza 95% con simulación
b0 <- c(); b1 <- c()
for (i in 1:5E3){
s <- sample(1:n, n, replace = T)
xs <- x[s]; y2s <- y2[s]
fit <- lm(y2s ~ xs)
b0[i] <- coef(fit)[1]
b1[i] <- coef(fit)[2]
}
quantile(b0, probs = c(0.05, 0.95))
## 5% 95%
## 53.11215 98.66261
quantile(b1, probs = c(0.05, 0.95))
## 5% 95%
## 3.740921 4.454484
Vemos como la longitud de los intervalos generados con bootstrap tiene menor longitud que el primer par de intervaos.
Finalmente calculamos los intervalos ajustando por mínimos cuadrados.
confint(fit3)
## 2.5 % 97.5 %
## (Intercept) 66.569167 78.955451
## x 3.989619 4.308103