Primera parte

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.

Primero cargamos los datos

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

Estimamos los parámetros del modelo

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

Graficamos los datos y el modelo ajustado

Estimamos las varianzas de los estimadores

vb0 <- (1/n + mean(x)**2/sxx) * s2; vb0
## [1] 6.013039
vb1 <- s2 / sxx; vb1
## [1] 0.03086201

Calculamos los intervalos de confianza 95%

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.

Segunda parte

Efecto del ingreso en la educación

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.

Modelo con o sin intercepto

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.

Hipótesis sobre la varianza del modelo

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.

Intervalo de confianza de la respuesta media

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.

Intervalo de predicción

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

Tercera parte

Bandas de confianza para la respuesta media

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í