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.


Segunda parte

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