Ajustar por MCO un modelo de RLS a los datos de Galton

Lectura de los datos
library(HistData)

Gráfica

datos <- data.frame(table(Galton))
library(dplyr)
datos <- filter(datos, Freq > 0)
datos$Freq <- 1 + 3 * (datos$Freq - 1)/ 47
datos$parent <- as.numeric(as.character(datos$parent))
datos$child  <- as.numeric(as.character(datos$child))

plot(datos$parent, datos$child, cex = datos$Freq, col = 'steelblue3', pch = 16,
     main="Regresión a la media", xlab="Estatura del Parde", ylab="Estatura del Hijo")
grid(5, 5, lwd = 2)

Ajuste: solución analítica

Recordemos que los estimadores de MCO son

\[\begin{align*} \hat{\beta}_1 &= \frac{S_{xy}}{S_{xx}} = \frac{\sum(x_i - \bar{x}_n)(y_i - \bar{y}_n)}{\sum(x_i - \bar{x}_n)^2}\\ \hat{\beta}_0 &= \bar{y}_n - \hat{\beta}_1 \bar{x}_n \end{align*}\] 
x <- Galton$parent
y <- Galton$child

Sxx <- sum((x - mean(x))^2)
Sxy <- sum((x - mean(x)) * (y - mean(y)))

b1_h <- Sxy / Sxx
b0_h <- mean(y) - b1_h * mean(x)

b0_h
## [1] 23.94153
b1_h
## [1] 0.6462906

El valor de la suma de cuadrados envaluada en \((\hat{\beta}_0, \hat{\beta}_1)\) es

q2_min <- sum((y - b0_h - b1_h * x)**2)

Aunque el método de mínimos cuadrados ordinarios no proporciona un estimador para \(\sigma^2\), un estimador razonable está dado por

\[\begin{equation*} \hat{\sigma}^2 = \frac{1}{n-2} \sum_{i=1}^n (y_i - \hat{\beta}_0 - \hat{\beta}_1x_i)^2 \end{equation*}\]

En palabras, la suma de cuadrados evaluada en los estimadores MCO de \(\beta_0\) y \(\beta_1\) dividida por \(n-2\).

Entonces, la estimación de la varianza para este conjunto de datos es

s_h <- q2_min / (dim(Galton)[1] - 2)
s_h
## [1] 5.011094

Ajuste: con la función lm de R.

La función de R que podemos utilizar para ajustar modelos de regresión lineal (simple o múltiple) es lm. Para consultar la documentación relacionada con esta función podemos utilizar el siguiente código:

help('lm')

Una vez revisada la información de la función lm, ajustamos un modelo RLS a los datos de Galton con el siguiente código

ajuste_lm <- lm(y ~ x)
ajuste_lm
## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##     23.9415       0.6463

Podemos oservar que las estimaciones de \(\beta_0\) y \(\beta_1\) coinciden con la solución analítica de la sección anterior. El método que usa la función lm aproxima numéricamente a la solución de las ecuaciones normales.

Como obtuvimos las mismas estimaciones de \(\hat{\beta}_0\) y \(\hat{\beta}_1\), el valor mínimo de la suma de cuadrados de los residuos y \(\hat{\sigma}^2\) son los mismos que en la sección anterior.

Podemos utilizar el siguiente código para calcular la suma de cuadrados a partir del resultado del ajuste con la función lm. La función residuals nos permite obtener los residuos del modelo ajustado.

sum(residuals(ajuste_lm) ** 2)
## [1] 4640.273

Ajuste: aproximación numérica

Finalmente, podemos obtener los estimadores MCO de \(\beta_0\) y \(\beta_1\) optimizando numéricamente la función suma de cuadrados de los residuos. Primero debemos declarar la función que calcula la suma de cuadrados para valores de \(\beta_0\) y \(\beta_1\) dados

q2 <- function(beta)
  sum((y - beta[1] - beta[2]*x) ** 2)

El argumento de la función q2 es un vector con dos entradas, la primera corresponde al valor de \(\beta_0\) y la segunda a \(\beta_1\), esto para permitir la optimización numérica utilizando la funicón nlm de R. De igual manera, se puede consultar la documentación relacaionada con la función nlm con el código help('nlm').

Por el momento no nos interesan todos los argumentos de la función nlm nos basta con indicar la función a minimizar y un par de valores iniciales para \(\beta_0\) y \(\beta_1\).

ajuste_nlm <- nlm(q2, c(0.5, 0.5))
ajuste_nlm
## $minimum
## [1] 4640.273
## 
## $estimate
## [1] 23.9415301  0.6462906
## 
## $gradient
## [1] -5.904747e-05 -4.033109e-03
## 
## $code
## [1] 1
## 
## $iterations
## [1] 22

De todos los resultados que proporciona nlm nos interesan tres:

Observamos nuevamente que los estimadores obtenidos son los mismos que en las secciones anteriores.

Gráfica con recta ajustada

plot(datos$parent, datos$child, cex = datos$Freq, col = 'steelblue3', pch = 16,
     main="Regresión a la media", xlab="Estatura del Parde", ylab="Estatura del Hijo")
grid(5, 5, lwd = 2)
abline(a = b0_h, b = b1_h, col = 'red2', lwd = 2)

Aplicaciones del modelo

¿Cuál es la estatura esperada de un hijo con un padre de 70 in?

x_new <- 70 
y_new <- b0_h + b1_h * x_new
y_new
## [1] 69.18187