Para este ejemplo consideramos los datos sobre desarrollo humano de las entidades del país en los años 2008, 2010 y 2014. El archivo se puede descargar desde aquí. La fuente original de los datos esta acá.
Primero hacermos una exploración del conjunto de datos
file_url <- 'http://sigma.iimas.unam.mx/jsantibanez/Cursos/Regresion/2017_2/dh.csv'
idh.ent <- read.csv(file_url)
str(idh.ent)
## 'data.frame': 96 obs. of 10 variables:
## $ Ent : Factor w/ 32 levels "Ags","BC","BCS",..: 1 2 3 4 8 9 7 6 5 10 ...
## $ Year: int 2008 2008 2008 2008 2008 2008 2008 2008 2008 2008 ...
## $ Esp : num 75.1 73.1 75.3 74.4 75.1 ...
## $ APE : num 8.7 8.89 8.99 7.78 8.97 ...
## $ AEE : num 11.8 12.3 13.4 11.9 12.4 ...
## $ ING : num 17507 21346 25716 16127 17266 ...
## $ IS : num 0.847 0.817 0.851 0.837 0.848 ...
## $ IE : num 0.618 0.638 0.672 0.588 0.643 ...
## $ II : num 0.78 0.81 0.838 0.768 0.778 ...
## $ IDH : num 0.742 0.75 0.782 0.723 0.751 ...
El conjunto de datos contiene 10 variables, una alfanumérica y nueve numéricas. La descripción de las variables es la siguiente:
Ent
: entidad federativaYear
: Año de referenciaEsp
: Esperanza de vida al nacerAPE
: Años promedio de escolaridadAEE
: Años esperados de escolaridadING
: Ingreso nacional bruto (per cápita y en USD PPC)IS
: Índice de saludIE
: Índice de educaciónII
: Índice de ingresoIDH
: Índice de desarrollo humanoNos interesa explorar las relaciones entre ingreso y eduación, en específico, nos interesa saber cómo el nivel educativo (APE
) es explicado por el ingreso (ING
). En este caso usaremos el lograritmo del ingreso.
library(dplyr)
idh10 <- filter(idh.ent, Year == 2010)
x <- log(idh10$ING)
y <- idh10$APE
Primero graficamos las variables de interés:
par(mar = c(4.1, 4.1, 1, 1))
plot(x, y, xlab = 'log(Ingreso)', ylab = 'Educación',
pch = 16, axes = F, xlim = c(8.75, 10.5), ylim = c(5.5, 10.75))
axis(1, seq(8.75, 10.5, by = 0.25))
axis(2, seq(5.5, 10.75, by = 1), las = 2)
box()
Recordamos que los estimadores de minímos cuadrados son:
\[\begin{align*} \hat{\beta}_0 &= \bar{y}_n - \frac{S_{xy}}{S_{xx}} \bar{x}_n \\ \hat{\beta}_1 &= \frac{S_{xy}}{S_{xx}} \\ \hat{\sigma}^2 &= \frac{1}{n-2}\sum_{i=1}^n (y_i - \hat{\beta}_0 - \hat{\beta}_1x_i)^2 \end{align*}\]sxx <- sum((x-mean(x))**2)
sxy <- sum((x-mean(x)) * (y - mean(y)))
b0.h <- mean(y) - sxy/sxx * mean(x); b0.h
## [1] -21.04291
b1.h <- sxy / sxx; b1.h
## [1] 3.018977
s2.h <- sum((y - b0.h - b1.h * x)**2) / (32 - 2); s2.h
## [1] 0.1861986
Entonces \(\hat{\beta}_0 =\) -21.043, \(\hat{\beta}_1 =\) 3.019 y \(\hat{\sigma}^2 =\) 0.186. Podemos graficar la recta ajustada con la nube de puntos.
par(mar = c(4.1, 4.1, 1, 1))
plot(x, y, xlab = 'log(Ingreso)', ylab = 'Educación',
pch = 16, axes = F, xlim = c(8.75, 10.5), ylim = c(5.5, 10.75))
axis(1, seq(8.75, 10.5, by = 0.25))
axis(2, seq(5.5, 10.75, by = 1), las = 2)
box()
abline(a = b0.h, b = b1.h, col = 'red2')
También podemos estimar las varianzas de nuestros estimadores. Recordemos:
\[\begin{align*} V(\hat{\beta}_0) &= \sigma^2 \left(\frac{1}{n} + \frac{\bar{x}_n^2}{S_{xx}} \right) \\ V(\hat{\beta}_1) &= \frac{\sigma^2}{S_{xx}} \end{align*}\]Podemos reemplazar \(\sigma^2\) por \(\hat{\sigma}^2\), entonces:
v_b0.h <- (1/96 + mean(x)**2/sxx)*s2.h; v_b0.h
## [1] 6.00916
v_b1.h <- s2.h / sxx; v_b1.h
## [1] 0.06423522
De lo anterior se sigue que \(V(\hat{\beta}_0) =\) 6.009 y \(V(\hat{\beta}_1) =\) 0.064.