7.4 Diagnóstico del modelo

En el análisis de las encuestas de hogares cuando se ajusten modelos estadísticos es importante realizar verificaciones de calidad y con esto tener certezas de las conclusiones que se obtienen. La mayoría de textos académicos dan un panorama bastante detallado de los supuestos y consideraciones que se deben tener en cuenta para tener un modelo correctamente definido. A continuación, se enlistan algunas de ellas:

  • Determinar si el modelo proporciona un adecuado ajuste a los datos.
  • Examinar si los errores están normalmente distribuidos.
  • Examinar si los errores tienen varianza constante.
  • Verificar si los errores se pueden asumir no correlacionados.
  • Determinar si alguno de los datos tiene valores con un efecto inusualmente grande sobre el modelo de regresión estimado, estos se conocen como datos influyentes.
  • Determinar si algún punto no sigue la tendencia de la mayoría de los datos cuando se toma en cuenta el modelo, estos puntos se conocen como outliers.

En este capítulo se abordarán alguno de los supuestos que se deben tener en cuenta al momento de ajustar un modelo de regresión lineal.

7.4.1 Coeficientes de determinación

Una medida de la bondad del ajuste en un modelo de regresión es el coeficiente de determinación o coeficiente de correlación múltiple (\(R^{2})\)). Dicho parámetro estima la proporción de la varianza de la población explicada por la regresión y oscila entre cero y uno. Entre más cercano esté de uno significa que mayor variabilidad explica y lo contrario ocurrirá si está cerca de cero. Lo anterior, en ocasiones es muy ambiguo puesto que, por ejemplo, en algunas disciplinas es posible obtener \(R^{2}\) con valores muy altos, mientras que en otras se obtienen \(R^{2}\) bajos. El cálculo de este parámetro a nivel poblacional se lleva a cabo de la siguiente manera:

\[ R^{2} = 1-\frac{SSE}{SST} \]

Donde, \(SST= \sum_{i=1}^N (y_i - \bar{y})^2\) es la suma de cuadrados totales, que representa la variabilidad total en la variable dependiente, y \(SSE= \sum_{i=1}^N (y_i - x_i \beta)^2\) es la suma de cuadrados del error, que representa la variabilidad no explicada por el modelo de regresión. El estimador de este parámetro usando muestras complejas está dado por:

\[ \widehat{R}_{\omega}^{2} = 1-\frac{\widehat{SSE}_{\omega}}{\widehat{SST}_{\omega}} \]

En donde \(\widehat{SSE}_{\omega}\) es la suma de cuadrados del error estimada, dada por:

\[ \widehat{SSE}_{\omega} = \sum_{h}^{H}\sum_{\alpha}^{a_{h}}\sum_{i=1}^{n_{h\alpha}}\omega_{h\alpha i}(y_{h\alpha i}-x_{h\alpha i}\hat{\beta})^{2} \]

Para continuar con los modelos ajustados en la sección anterior, se procede a estimar los \(R^{2}\) utilizando R. Inicialmente, se procede a estimar los parámetros del modelo utilizando la función svyglm de survey como se mostró anteriormente y también, se ajusta un modelo solo con el intercepto para obtener la estimación de la SST:

modNul <- svyglm(Income ~ 1, design = diseno_qwgt)

s1 <- summary(fit_svy)
s0 <- summary(modNul)

wSST <- s0$dispersion
wSSE <- s1$dispersion

Por tanto, la estimación del \(R^{2}\) es:

R2 = 1- wSSE/wSST
R2
##      variance    SE
## [1,]    0.513 20579

7.4.2 Residuales estandarizados

En el diagnóstico de los modelos, es crucial el análisis de los residuales. Estos análisis proporcionan, bajo el supuesto que el modelo ajustado es adecuado, una estimación de los errores. Por tanto, un estudio cuidadoso de los residuales deberá ayudar al investigador a concluir si el procedimiento de ajuste no ha violado los supuestos o si, por el contrario, uno o varios de los supuestos no se verifican y hay necesidad de revisar el procedimiento de ajuste.

Para realizar el análisis de los residuales, en primera instancia, se definen los residuales de Pearson (Heeringa, West, y Berglund 2017) como sigue:

\[ r_{p_{i}} = (t(y_{i}-\mu_{i}(t(\hat{\beta}_{\omega}))\sqrt{\frac{\omega_{i}}{V(t(\hat{\mu}_{i})}} \]

Donde, \(\mu_{i}\) es el valor esperado de \(y_{i}\), \(w_{i}\) es la ponderación de la encuesta para el i-ésimo individuo del diseño muestral complejo, Por último, \(V(\mu_{i})\) es la función de varianza del resultado. Sobre estos residuales se realizan los análisis de normalidad y varianza constante.

Si el supuesto de varianza constante no se cumple, los estimadores siguen siendo insesgados y consistentes, pero dejan de ser eficientes, es decir, dejan de ser los mejores en cuanto a que ya no tienen la menor varianza entre todos los estimadores insesgados. Una de las formas de analizar el supuesto de varianzas constantes en los errores es hacerlo de manera gráfica. Para ello, se grafica los residuos del modelo contra \(\hat{y}\) o los residuos del modelo contra \(X_{i}\). Si al realizar estos gráficos se logra evidenciar cualquier patrón que no sea una nube de puntos constante, se puede concluir que la varianza de los errores no es constante.

A manera de ejemplificar los conceptos vistos, se van a utilizar los modelos previamente ajustados. En primero instancia, el análisis del modelo se centrará en los supuestos de normalidad y varianza constante en los errores. La librería svydiags está pensada en ayudar en el diagnostico de modelos de regresión lineal, siendo una extensión más para complementar el paquete survey. Con las librerías svydiags se extraen los residuales estandarizados como sigue:

library(svydiags)
stdresids = as.numeric(svystdres(fit_svy)$stdresids)
diseno_qwgt$variables %<>% mutate(stdresids = stdresids)

Una forma muy común para hacer la evaluación de la normalidad de los residuales es realizar un gráfico cuantil-cuantil normal o QQplot. El QQplot es una gráfica de cuantiles para los residuos observados frente a los calculados a partir de una distribución normal teórica que tiene la misma media y varianza que la distribución de los residuos observados. Por lo tanto, una línea recta de 45° en este gráfico sugeriría que la normalidad es una suposición razonable para los errores aleatorios en el modelo.

qqnorm(stdresids)
qqline(stdresids, col = 2)

También podemos hacer el análisis de normalidad por medio del histograma de los residuales estandarizados con el siguiente código, en cuya salida se puede observar que el histograma de los residuales estandarizados (barras y línea azul) no necesariamente sigue el comportamiento de una distribución normal (línea roja).

ggplot(data = diseno_qwgt$variables,
       aes(x = stdresids)) +
  geom_histogram(
    aes(y = ..density..),
    colour = "black",
    fill = "blue",
    alpha = 0.3
  ) +
  geom_density(size = 2, colour = "blue") +
  geom_function(fun = dnorm,
                colour = "red",
                size = 2) +
  theme_cepal() + labs(y = "")

Como se puedo observar en el QQplot, hay evidencia gráfica de que los errores no se distribuyen según una distribución normal. Por otro lado, un análisis que se realiza de manera gráfica es la verificación de que las varianzas de los residuales sean constantes, condicionadas a los valores de las covariables. Esto se puede observar a través de diagramas de dispersión de los residuales estandarizados contra las covariables usadas en el modelo. Primero, se agregan las predicciones a la base de datos para poder realizar las gráficas.

library(patchwork)
diseno_qwgt$variables %<>%
  mutate(pred = predict(fit_svy))

g2 <- ggplot(data = diseno_qwgt$variables,
             aes(x = Expenditure, y = stdresids)) +
  geom_point() +
  geom_hline(yintercept = 0) + theme_cepal()

g3 <- ggplot(data = diseno_qwgt$variables,
             aes(x = Age2, y = stdresids)) +
  geom_point() +
  geom_hline(yintercept = 0) + theme_cepal()

g4 <- ggplot(data = diseno_qwgt$variables,
             aes(x = Zone, y = stdresids)) +
  geom_point() +
  geom_hline(yintercept = 0) + theme_cepal()
g5 <- ggplot(data = diseno_qwgt$variables,
             aes(x = Sex, y = stdresids)) +
  geom_point() +  geom_hline(yintercept = 0) +
  theme_cepal()

(g2 | g3) / (g4 | g5)

Como se puede observar en las gráficas de gastos y edad al cuadrado, ambas muestran que las varianzas no son constantes, puesto que la dispersión crece o decrece a medida que cambia el valor de estas covariables.

7.4.3 Observaciones influyentes

Otras técnicas utilizadas también para el análisis de los modelos consisten en el análisis de observaciones influyentes. Una observación se denomina influyente si al removerla de la nube de puntos esta causa un cambio grande en el ajuste del modelo. Es importante resaltar que un punto influyente podría o no ser un dato atípico. Para detectar observaciones influyentes es necesario tener claro qué tipo de influencia se quiere detectar. Lo anterior puesto que, por ejemplo, una observación puede ser influyente sobre la estimación de los parámetros, pero no para la estimación de la varianza del error. A continuación, se presentan las distintas técnicas estadísticas para la detección de datos influyentes:

  1. Distancia de Cook: diagnostica si la i-ésima observación es influyente en la estimación del modelo, por estar lejos del centro de masa de los datos. En general, distintos autores consideran que las observaciones son influyentes cuando esta cantidad es mayor que 2 o 3.
  2. Estadístico \(D_fBeta_{(i)}\): este estadístico mide el cambio en la estimación del vector de coeficientes de regresión cuando la observación es eliminada. Se determina que la i-ésima observación es influyente para \(B_j\) si \(\mid D_{f}Betas_{(t(i)j}\mid\geq\frac{z}{\sqrt{n}}\) con \(z = 2\). Como alternativa, se puede usar \(t_{0.025,n-p}/\sqrt(n)\) donde \(t_{0.025,n-p}\) es el percentil \(97.5\)
  3. Estadístico \(D_{f}Fits_{(t(i)}\)*: este estadístico mide el cambio en el ajuste del modelo cuando se elimina una observación particular. En esta instancia, la \(i\)-ésima observación se considera influyente en el ajuste del modelo si \(\mid DfFits(t(i)\mid\geq z\sqrt{\frac{p}{n}}\) con \(z = 2\).

Para ejemplificar los conceptos definidos, se seguirá con los modelos ajustados en la sección anterior. Una vez ajustados estos modelos y verificados los supuestos, se procede a hacer el cálculo de la distancia de Cook usando la función svyCooksD del paquete svydiags como sigue:

library(svydiags)

d_cook = data.frame(cook = svyCooksD(fit_svy),
                    id = 1:length(svyCooksD(fit_svy)))

ggplot(d_cook, aes(y = cook, x = id)) +
  geom_point() +
  theme_bw(20)

Como se puede observar, ninguna de las distancias de Cook’s es mayor a 3; por ende podemos afirmar que no existen observaciones influyentes. Ahora bien, si se desea corroborar que no hay observaciones influyentes utilizando el estadístico \(D_{f}Betas_{(t(i)j}\), esto se puede realizar con la función svydfbetas como se muestra a continuación:

d_dfbetas <- data.frame(t(svydfbetas(fit_svy)$Dfbetas))
colnames(d_dfbetas) <- paste0("Beta_", 0:4)
d_dfbetas %>% slice(1:10)
Beta_0 Beta_1 Beta_2 Beta_3 Beta_4
0.0005 -2e-04 0.0020 -0.0045 -0.0076
-0.0005 -1e-04 0.0013 0.0026 -0.0031
-0.0008 -1e-04 0.0008 0.0022 0.0008
-0.0004 -1e-04 0.0011 -0.0031 0.0007
-0.0008 0e+00 0.0008 0.0021 0.0014
0.0009 5e-04 -0.0037 -0.0066 0.0103
0.0027 4e-04 -0.0032 -0.0080 -0.0029
0.0011 3e-04 -0.0029 0.0081 -0.0045
0.0030 3e-04 -0.0031 -0.0082 -0.0054
-0.0003 4e-04 0.0012 -0.0039 -0.0042

Una vez calculado los \(D_{f}Betas_{(t(i)j}\) se procede a reacomodar la salida para verificar cuáles observaciones son influyentes. Para esto, de calcula el umbral (cutoff) para definir si es o no influyente la observación. Ese umbral es tomado de las salidas de la función svydfbetas. Por último, se genera una variable dicotómica que indique si la observación es o no influyente como se muestra a continuación:

d_dfbetas$id <- 1:nrow(d_dfbetas)
d_dfbetas <- reshape2::melt(d_dfbetas, id.vars = "id")
cutoff <- svydfbetas(fit_svy)$cutoff
d_dfbetas %<>% mutate(Criterio = ifelse(abs(value) > cutoff, "Si", "No"))

tex_label <- d_dfbetas %>%
  filter(Criterio == "Si") %>%
  arrange(desc(abs(value))) %>%
  slice(1:10L)
tex_label
id variable value Criterio
889 Beta_0 0.3284 Si
890 Beta_1 -0.3043 Si
891 Beta_0 0.3037 Si
889 Beta_1 -0.2973 Si
891 Beta_1 -0.2934 Si
890 Beta_0 0.2903 Si
889 Beta_4 -0.2526 Si
890 Beta_4 -0.2270 Si
2311 Beta_4 0.2122 Si
890 Beta_3 0.2029 Si

Como se pudo observar en la salida anterior hay varias observaciones que resultan influyentes dado el criterio del \(D_{f}Betas_{(t(i)j}\). A continuación, y de manera ilustrativa, se grafican los \(D_{f}Betas_{(t(i)j}\) y el umbral con el fin de ver de manera gráfica aquellas observaciones influyentes, teniendo en cuenta que, aquellos puntos rojos en la gráfica representan observaciones influyentes.

ggplot(d_dfbetas, aes(y = abs(value), x = id)) +
  geom_point(aes(col = Criterio)) +
  geom_text(data = tex_label,
            angle = 45,
            vjust = -1,
            aes(label = id)) +
  geom_hline(aes(yintercept = cutoff)) +
  facet_wrap(. ~ variable, nrow = 2) +
  scale_color_manual(values = c("Si" = "red", "No" = "black")) +
  theme_cepal()

Si el objetivo es detectar observaciones influyentes pero considerando ahora la estadística \(D_{f}Fits_{(t(i)}\), se utiliza la función svydffits y se siguen los mismos pasos mostrados para el estadístico \(D_{f}Betas_{(t(i)j}\).A continuación se muestra el código computacional apropiado para realizarlo:

d_dffits <- data.frame(dffits = svydffits(fit_svy)$Dffits,
                      id = 1:length(svydffits(fit_svy)$Dffits))

cutoff <- svydffits(fit_svy)$cutoff

d_dffits %<>% mutate(C_cutoff = ifelse(abs(dffits) > cutoff, "Si", "No"))
ggplot(d_dffits, aes(y = abs(dffits), x = id)) +
  geom_point(aes(col = C_cutoff)) +
  geom_hline(yintercept = cutoff) +
  scale_color_manual(values = c("Si" = "red", "No" = "black")) +
  theme_cepal()

Como se puede observar en el gráfico anterior, también hay observaciones influyentes utilizando \(D_{f}Fits_{(t(i)}\), las cuales se muestran en rojo en el gráfico.

References

Heeringa, Steven G., Brady T. West, y Patricia A. Berglund. 2017. Applied survey data analysis. Chapman y Hall CRC statistics en the social y behavioral sciences series. CRC Press.