Considerar el conjunto de datos de fecundidad de la ONU (disponible aquí), que contiene la 193 observaciones de tres variables, cada una corresponde a un país. Las variables son:
logPPgdp: logaritmo base 2 del producto interno bruto per capita.logFertility: logaritmo base 2 del promedio de hijos nacidos vivos de las mujeres de 15 a 40 años.Purban: poncentaje de población en zonas urbanas.El objetivo es ajustar un modelo RLM para explicar la fecundidad a partir del PIB y la urbanización.
Rfile_url <- 'http://users.stat.umn.edu/~sandy/alr3ed/website/data/UN2.txt'
datos <- read.table(file_url, T)
rownames(datos) <- datos$Locality
datos <- datos[,-4]You must enable Javascript to view this page properly.
Primero ajustamos el modelo a mano, es decir, con las expresiones que calculamos en clase
\[ \hat{\boldsymbol{\beta}} = \left(\mathbf{X}'\mathbf{X}\right)^{-1}\mathbf{X}'\mathbf{y} \qquad \text{y} \qquad \hat{\sigma}^2 = \frac{1}{n-(p+1)}(\mathbf{y}-\mathbf{X}\hat{\boldsymbol{\beta}})' (\mathbf{y}-\mathbf{X}\hat{\boldsymbol{\beta}}) \]
Declaramos el vector respuesta y la matriz de diseño:
y <- datos[,'logFertility']
X <- cbind(1, datos[,c('logPPgdp', 'Purban')])
X <- as.matrix(X)
Verificamos que la matriz de diseño está bien declarada:
head(X, 6)## 1 logPPgdp Purban
## Afghanistan 1 6.614710 22
## Albania 1 10.363040 43
## Algeria 1 10.800900 58
## Angola 1 9.529431 35
## Argentina 1 12.806348 88
## Armenia 1 9.424166 67
Ahora estimamos \(\hat{\boldsymbol{\beta}}\) y \(\sigma^2\):
bh <- solve(t(X) %*% X) %*% t(X) %*% y; round(bh, 5)## [,1]
## 1 2.59300
## logPPgdp -0.12548
## Purban -0.00352s2h <- t(y - X%*%bh) %*% (y - X%*%bh) / (193 - 3); round(s2h, 5)## [,1]
## [1,] 0.15494
También podemos estimar la varianza de \(\hat{\boldsymbol{\beta}}\):
vbh <- as.numeric(s2h) * solve(t(X) %*% X); round(vbh, 6)## 1 logPPgdp Purban
## 1 0.021569 -0.002450 0.000111
## logPPgdp -0.002450 0.000365 -0.000028
## Purban 0.000111 -0.000028 0.000004
lmLa función lm de R nos permite ajuster modelos RLM. Hay varias formas de utilizarla, todas dan el mismo resultado pero resultan prácticas según el número de variables explicativas.
ajuste1 <- lm(y ~ X[,-1])
summary(ajuste1)##
## Call:
## lm(formula = y ~ X[, -1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.05378 -0.16952 0.06835 0.25587 0.89290
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.592996 0.146864 17.656 < 2e-16 ***
## X[, -1]logPPgdp -0.125475 0.019095 -6.571 4.67e-10 ***
## X[, -1]Purban -0.003522 0.001884 -1.869 0.0631 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3936 on 190 degrees of freedom
## Multiple R-squared: 0.4689, Adjusted R-squared: 0.4633
## F-statistic: 83.88 on 2 and 190 DF, p-value: < 2.2e-16ajuste2 <- lm(logFertility ~ ., data = datos)
summary(ajuste2)##
## Call:
## lm(formula = logFertility ~ ., data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.05378 -0.16952 0.06835 0.25587 0.89290
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.592996 0.146864 17.656 < 2e-16 ***
## logPPgdp -0.125475 0.019095 -6.571 4.67e-10 ***
## Purban -0.003522 0.001884 -1.869 0.0631 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3936 on 190 degrees of freedom
## Multiple R-squared: 0.4689, Adjusted R-squared: 0.4633
## F-statistic: 83.88 on 2 and 190 DF, p-value: < 2.2e-16ajuste3 <- lm(logFertility ~ logPPgdp + Purban, data = datos)
summary(ajuste3)##
## Call:
## lm(formula = logFertility ~ logPPgdp + Purban, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.05378 -0.16952 0.06835 0.25587 0.89290
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.592996 0.146864 17.656 < 2e-16 ***
## logPPgdp -0.125475 0.019095 -6.571 4.67e-10 ***
## Purban -0.003522 0.001884 -1.869 0.0631 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3936 on 190 degrees of freedom
## Multiple R-squared: 0.4689, Adjusted R-squared: 0.4633
## F-statistic: 83.88 on 2 and 190 DF, p-value: < 2.2e-16
Los coeficientes estimados coinciden con los resultados anteriores
coef(ajuste2)## (Intercept) logPPgdp Purban
## 2.592995910 -0.125475126 -0.003522218También es posible obtener el estimador de \(\sigma^2\) y de las varianzas de las componetes de \(\boldsymbol{\beta}\), aunque no las covarianzas.
You must enable Javascript to view this page properly.