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:

El objetivo es construir verificar el supuesto de linealidad del modelo RLM de fecundidad sobre el producto y la población en zonas urbanas.

Cargar los datos en 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]

Respuesta vs. variables explicativas

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.

Residuos vs. variables explicativas

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.

Pruebas de falta de ajuste

Ahora procedemos a verificar nuestras sospechas de no linealidad con las pruebas de falta de ajuste.

Para 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.

Para 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.

Prueba de no aditividad de Tukey

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.

Corrección de la 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.