Considerar el conjunto de datos de desarrollo humano 2010 de las entidades del país. Construir intervalos de confianza \(95\%\) para los parámetros del modelo RLS de logaritmo del ingreso contra años promedio de escolaridad.
file_url <- 'http://sigma.iimas.unam.mx/jsantibanez/Cursos/Regresion/2017_2/dh.csv'
idh.ent <- read.csv(file_url)
library(dplyr)
idh10 <- filter(idh.ent, Year == 2010)
x <- log(idh10$ING, 2)
y <- idh10$APE
sxx <- sum((x - mean(x))**2)
sxy <- sum((x - mean(x)) * (y - mean(y)))
b0 <- mean(y) - (sxy / sxx) * mean(x); b0
## [1] -21.04291
b1 <- sxy / sxx; b1
## [1] 2.092595
n <- length(x)
s2 <- sum((y - b0 - b1*x)**2) / (n - 2); s2
## [1] 0.1861986
vb0 <- (1/n + mean(x)**2/sxx) * s2; vb0
## [1] 6.013039
vb1 <- s2 / sxx; vb1
## [1] 0.03086201
tt <- qt(0.025, n-2, lower.tail = F)
IC_b0 <- b0 + c(-1, 1) * tt * sqrt(vb0); IC_b0
## [1] -26.05087 -16.03495
IC_b1 <- b1 + c(-1, 1) * tt * sqrt(vb1); IC_b1
## [1] 1.733817 2.451373
cl <- qchisq(0.025, n-2)
cu <- qchisq(0.975, n-2)
IC_s2 <- (n-2) * s2 / c(cu, cl); IC_s2
## [1] 0.1189027 0.3326803
Cuando el modelo se ajusta con la función lm
de R
, es posible obenter los intervalos de confianza para \(\beta_0\) y \(\beta_1\) con la función confint
.
modelo <- lm(y ~ x)
confint(modelo)
## 2.5 % 97.5 %
## (Intercept) -26.050866 -16.034950
## x 1.733817 2.451373
Los resultados coinciden con los anteriores intervalos calculados anteriormente.
Lo que se pide es contrastar las hipótesis
\[\begin{equation*} H_0: \beta_1 = 0 \qquad \text{vs.} \qquad H_1: \beta_1 \neq 0 \end{equation*}\]Calculamos el estadístico de prueba y lo comparamos con el valor crítico
tp <- abs(b1/sqrt(vb1)); tp
## [1] 11.91168
tt
## [1] 2.042272
Como el estadístico de prueba es mayor que el valor crítico, se rechaza \(H_0\) con una significancia \(\alpha = 0.05\). Los datos indican que hay algún efecto del ingreso en la educación.
Lo que se pide es contrastar las hipótesis
\[\begin{equation*} H_0: \beta_0 = 0 \qquad \text{vs.} \qquad H_1: \beta_0 \neq 0 \end{equation*}\]Calculamos el estadístico de prueba y lo comparamos con el valor crítico
tp <- abs(b0/sqrt(vb0)); tp
## [1] 8.581412
tt
## [1] 2.042272
Como el estadístico de prueba es mayor que el valor crítico, se rechaza \(H_0\) con una significancia \(\alpha = 0.05\). Los datos indican que el modelo de regresión debe incluir a \(\beta_0\).
Los resultados de las pruebas para contrastar \(H_0: \beta_0 = 0\) y \(H_0: \beta_1\) se pueden consultar utilizando la función summary
en el modelo ajustado con la función lm
.
summary(modelo)
##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.87394 -0.23529 -0.03492 0.19134 1.27516
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -21.0429 2.4521 -8.581 1.42e-09 ***
## x 2.0926 0.1757 11.912 6.71e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4315 on 30 degrees of freedom
## Multiple R-squared: 0.8255, Adjusted R-squared: 0.8197
## F-statistic: 141.9 on 1 and 30 DF, p-value: 6.71e-13
Se puede observar como los resultados coinciden con los calculados previamente.
Lo que se pide es contrastar las hipótesis
\[\begin{equation*} H_0: \sigma^2 \leq 0.25 \text{vs.} \qquad H_1: \sigma^2 > 0.25 \end{equation*}\]Calculamos el estadístico de prueba y lo comparamos con el valor crítico
tp <- (n - 2)*s2 / 0.25; tp
## [1] 22.34384
qchisq(0.95, n - 2)
## [1] 43.77297
El estadístico de prueba es menor que el valor critico, no se rechaza \(H_0\) con una significancia \(\alpha = 0.05\). Los datos indican que la varianza no es mayor a 0.25.
Primero debemos obtener la estimación puntual y su varianza estimada
m40 <- b0 + b1 * log(40000, 2); m40
## [1] 10.94808
vm40 <- (1/n + (log(40000, 2) - mean(x))^2/sxx) * s2; vm40
## [1] 0.06091164
Ahora se calcula el intervalo de confianza
IC_m40 <- m40 + c(-1, 1) * tt * sqrt(vm40); IC_m40
## [1] 10.44405 11.45212
Este intervalo se puede calcular directamente utilizando la función predict
. Esta función se utiliza para hacer predicciones (calcular valores ajustados) en un nuevo conjunto de datos a partir de un modelo ajustado previamente. El parámetro interval
se utiliza para especificar si se requiere un intervalo de confianza o de predicción.
predict(modelo, data.frame(x = log(40000, 2)), interval = 'confidence')
## fit lwr upr
## 1 10.94808 10.44405 11.45212
Se puede observar que los resultados coinciden con los calculados previamente.
Ya tenemos la predicción puntual y hace falta la varianza de predicción
vp40 <- (1 + 1/n + (log(40000, 2) - mean(x))^2/sxx) * s2; vp40
## [1] 0.2471103
Ahora se calcula el intervalo de predicción
IC_p40 <- m40 + c(-1, 1) * tt * sqrt(vp40); IC_p40
## [1] 9.932866 11.963301
De nuevo, se puede utilizar la función predict
para obtener el intervalo de predicción, lo cual se debe especificar en el parámetro interval
.
predict(modelo, data.frame(x = log(40000, 2)), interval = 'prediction')
## fit lwr upr
## 1 10.94808 9.932866 11.9633
Las bandas de confianza son generalizaciones de los intervalos de confianza que se obtienen cuando se hacen inferencias sobre \(\beta_0 + \beta_1 x\) vista como función de \(x\). Se dice que estadísticos \(L(x, \mathbf{X})\) y \(U(x, \mathbf{X})\) conforman una banda de confianza \(100(1-\alpha)\%\) para \(\beta_0 + \beta_1 x\) si se cumple \[ P \left\lbrace L(x,\mathbf{X}) \leq \beta_0 + \beta_1 x \leq U(x, \mathbf{X}), \forall x \in \mathbb{R} \right\rbrace \geq 1 - \alpha \]
La interpretación de las bandas es la siguiente. Existe una verdadera función que relaciona \(E(Y \,\vert\, X = x)\) con \(x\) y tiene la forma \(\beta_0 + \beta_1 x\), por lo que es una recta. Si se observa un número grande de muestras y con cada una se calculan las funciones \(L(x, \mathbf{X})\) y \(U(x, \mathbf{X})\), resulta que en al menos el \(100(1-\alpha)\%\) de los casos la función \(\beta_0 + \beta_1 x\) está contenida entre \(L(x, \mathbf{X})\) y $U(x, ), para todo valor de \(x\).
El método de Hotelling-Scheffé para construir intervalos de confianza simultáneos nos permite calcular bandas de confinza para verdadera recta de regresión. utilizando la notación anterior, una banda de confianza \(100(1-\alpha)\%\) para \(\beta_0 + \beta_1 x\) está formada por
\[\begin{align*} \hat{\beta}_0 + \hat{\beta}_1 x &- \sqrt{2F^{(1-\alpha)}_{(2, n-2)}}\hat{\sigma}_{MCO} \sqrt{\left(\frac{1}{n} + \frac{(x-\bar{}_n)^2}{S_{xx}} \right)} \\ \hat{\beta}_0 + \hat{\beta}_1 x &+ \sqrt{2F^{(1-\alpha)}_{(2, n-2)}}\hat{\sigma}_{MCO} \sqrt{\left(\frac{1}{n} + \frac{(x-\bar{x}_n)^2}{S_{xx}} \right)} \\ \end{align*}\]Estas funciones se programan en el siguiente código
l_x <- function(x_n, alfa){
(b0 + b1*x_n) - sqrt(2*qf(1 - alfa, 2, n-2))*sqrt(1/n + (x_n - mean(x))^2/sxx)*sqrt(s2)
}
u_x <- function(x_n, alfa){
(b0 + b1*x_n) + sqrt(2*qf(1 - alfa, 2, n-2))*sqrt(1/n + (x_n - mean(x))^2/sxx)*sqrt(s2)
}
Finalmente, se agregan las gráficas de las bandas de confianza al gráfico de dispersión con la recta de regresión ajustada. Por construcción, la recta ajustada siempre queda al centro de las bandas, no así la verdadera recta de regresión, aunque esa es desconocida.
En la escala original de los datos las bandas de confianza 95% se ven así