library(HistData)
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)
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
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
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:
minimum
: valor mínimo de la función q2
,estimate
: los estimadores de MCO de \(\beta_0\) y \(\beta_1\),code
: el valor 1 nos indica la convergencia del algoritmo iterativo.Observamos nuevamente que los estimadores obtenidos son los mismos que en las secciones anteriores.
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)