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 construir verificar el supuesto de linealidad del modelo RLM de fecundidad sobre el producto y la población en zonas urbanas.
R
file_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]
Primero expliramos la relación marginal de la respuesta con cada variable explicativa. Agregamos una curva de suavizamiento para entender mejor la relación entre cada par de variables.
logFertility
vs. logPPgdp
plot(datos$logPPgdp, datos$logFertility, pch = 16, col = 'steelblue',
xlab = 'logPPgdp (v. explicativa 1)', ylab = 'logFertility (v. respuesta)')
points(loess.smooth(datos$logPPgdp, datos$logFertility), type = 'l', lwd = 2, lty = 2)
La gráfica sugiere que la relación marginal entre logFertility
y logPPgdp
es no lineal.
logFertility
vs. Purban
plot(datos$Purban, datos$logFertility, pch = 16, col = 'steelblue',
xlab = 'Purban (v. explicativa 2)', ylab = 'logFertility (v. respuesta)')
points(loess.smooth(datos$Purban, datos$logFertility), type = 'l', lwd = 2, lty = 2)
La gráfica sugiere que la relación marginal entre logFertility
y Purban
es no lineal.
Ahora exploramos la relación entre los residuos y cada variable explicativa. En este caso también agregamos una curva de suavizamiento.
Primero ajustamos el modelo y extraemos los residuos.
modelo <- lm(logFertility ~ ., datos)
residuos <- residuals(modelo)
Residuos
vs. logPPgdp
plot(datos$logPPgdp, residuos, pch = 16, col = 'steelblue',
xlab = 'logPPgdp (v. explicativa 1)', ylab = 'Residuos')
points(loess.smooth(datos$logPPgdp, residuos), type = 'l', lwd = 2, lty = 2)
Esta gráfica también sugiere que la relación entre logFertility
y logPPgdp
es no lineal.
Residuos
vs. Purban
plot(datos$Purban, residuos, pch = 16, col = 'steelblue',
xlab = 'Purban (v. explicativa 2)', ylab = 'Residuos')
points(loess.smooth(datos$Purban, residuos), type = 'l', lwd = 2, lty = 2)
Esta gráfica también sugiere que la relación entre logFertility
y Purban
es no lineal.
Ahora procedemos a verificar nuestras sospechas de no linealidad con las pruebas de falta de ajuste.
logPPgdp
Ajustamos el modelo
\[ \texttt{logFertility} = \beta_0 + \beta_1\texttt{logPPgdp} + \beta_2\texttt{Purban} + \eta_1\texttt{logPPgdp}^2 \]
lack.of.fit1 <- lm(logFertility ~ logPPgdp + Purban + I(logPPgdp^2), datos)
summary(lack.of.fit1)
##
## Call:
## lm(formula = logFertility ~ logPPgdp + Purban + I(logPPgdp^2),
## data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.08171 -0.23352 0.05329 0.24583 0.94550
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.603987 0.641007 7.182 1.53e-11 ***
## logPPgdp -0.507502 0.120142 -4.224 3.73e-05 ***
## Purban -0.002711 0.001857 -1.460 0.14595
## I(logPPgdp^2) 0.016951 0.005266 3.219 0.00152 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3843 on 189 degrees of freedom
## Multiple R-squared: 0.4965, Adjusted R-squared: 0.4885
## F-statistic: 62.13 on 3 and 189 DF, p-value: < 2.2e-16
Los resultados anteriores indican que se rechaza la hipótesis \(H_0: \eta_1 = 0\), es decir, hay evidencia para sostener que la relación entre logFertility
y logPPgdp
es no lineal.
Purban
Ajustamos el modelo
\[ \texttt{logFertility} = \beta_0 + \beta_1\texttt{logPPgdp} + \beta_2\texttt{Purban} + \eta_2\texttt{Purban}^2 \]
lack.of.fit2 <- lm(logFertility ~ logPPgdp + Purban + I(Purban^2), datos)
summary(lack.of.fit2)
##
## Call:
## lm(formula = logFertility ~ logPPgdp + Purban + I(Purban^2),
## data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.99912 -0.23736 0.06394 0.25862 0.94249
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.010e+00 1.892e-01 15.912 < 2e-16 ***
## logPPgdp -1.291e-01 1.863e-02 -6.930 6.45e-11 ***
## Purban -2.082e-02 5.452e-03 -3.818 0.000182 ***
## I(Purban^2) 1.593e-04 4.728e-05 3.368 0.000916 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3833 on 189 degrees of freedom
## Multiple R-squared: 0.499, Adjusted R-squared: 0.491
## F-statistic: 62.75 on 3 and 189 DF, p-value: < 2.2e-16
Los resultados anteriores indican que se rechaza la hipótesis \(H_0: \eta_n = 0\), es decir, hay evidencia para sostener que la relación entre logFertility
y Purban
es no lineal.
Ahora aplicaremos la prueba de no aditividad d Tukey. Primero ajustamos un modelo RLM y calculamos \(Z\) a partir de los valores ajustados, posterioremente ajustamos un modelo RLM incluyendo a \(Z\).
Z <- fitted(modelo)^2/(2*mean(datos$logFertility))
summary(lm(logFertility ~ logPPgdp + Purban + Z, datos))
##
## Call:
## lm(formula = logFertility ~ logPPgdp + Purban + Z, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.06422 -0.23936 0.04417 0.24320 0.95933
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.834426 0.949077 -0.879 0.380411
## logPPgdp 0.070907 0.056860 1.247 0.213922
## Purban 0.002642 0.002486 1.063 0.289304
## Z 1.610369 0.440880 3.653 0.000336 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3814 on 189 degrees of freedom
## Multiple R-squared: 0.5039, Adjusted R-squared: 0.4961
## F-statistic: 64 on 3 and 189 DF, p-value: < 2.2e-16
Los resultados anteriores indican no adividad.
De los resultados anteriores obtenemos que \(\hat{\gamma} = 1.61\), entonces \(1-\hat{gamma} \approx -0.5\), entonces aplicamos la transformación
\[ \texttt{sqlogFert} = \frac{1}{\sqrt{1 + \texttt{logFertility}}} \]
se suma 1 para evitar indeterminaciones.
datos$sqlogFert <- 1/sqrt(1 + datos$logFertility)
Ahora ajustamos el modelo
\[ \texttt{sqlogFert} = \beta_0 + \beta_1 \texttt{logPPgdp} + \beta_2 \texttt{Purban} + \epsilon \]
los resultados se presentan a continuación.
modelo3 <- lm(sqlogFert ~ logPPgdp + Purban, datos)
summary(modelo3)
##
## Call:
## lm(formula = sqlogFert ~ logPPgdp + Purban, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.15922 -0.04969 -0.02266 0.03833 0.23703
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4526117 0.0305630 14.809 < 2e-16 ***
## logPPgdp 0.0210704 0.0039737 5.302 3.15e-07 ***
## Purban 0.0007235 0.0003922 1.845 0.0666 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08192 on 190 degrees of freedom
## Multiple R-squared: 0.3855, Adjusted R-squared: 0.379
## F-statistic: 59.6 on 2 and 190 DF, p-value: < 2.2e-16
Ahora exploramos las graficas de residuos para este nuevo modelo.
plot(datos$logPPgdp, residuals(modelo3), pch = 16, col = 'steelblue',
xlab = 'logPPgdp (v. explicativa 1)', ylab = 'Residuos')
points(loess.smooth(datos$logPPgdp, residuals(modelo3)), type = 'l', lwd = 2, lty = 2)
plot(datos$Purban, residuals(modelo3), pch = 16, col = 'steelblue',
xlab = 'Purban (v. explicativa 2)', ylab = 'Residuos')
points(loess.smooth(datos$Purban, residuals(modelo3)), type = 'l', lwd = 2, lty = 2)
Vemos como se ha corregido el probema de linealidad pero ahora parece haber problemas de varianza no constante en los errores.