6.2 Diagnóstico del modelo

En el análisis de las encuestas de hogares cuando se ajusten modelo estadístico 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 (Tellez, 2016)

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

Estimación del \(R^{2}\) y \(R_{adj}^{2}\)

Una medida del ajuste del 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 0 y 1. 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 ambigüo puesto que por ejemplo, los físicos pueden obtienen \(R^{2}\) altísimos (mayores a 0.98–0.99) mientras que los químicos obtienen \(R^{2}\) superiores a 0.90. Sin embargo, los científicos sociales y demás que trabajen con poblaciones humanas encontrarán que su mejor modelo de regresión a menudo explicará solo entre el 20 % y el 40 % de la variación en la variable dependiente (Heringa). A continuación, se presenta como se calcula este coeficiente:

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

Donde, SST es la suma de cuadrados totales y SSE es la suma de cuadrados del error. El estimador de este parámetro usando muestras complejas está dado por:

\[ R_{\omega}^{2} = 1-\frac{WSSE}{WSST} \]

Donde, WSST es la variabilidad total estimada y WSSE es la suma de cuadrados estimada y se define como:

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

Por último, \(R_{adj}^{2}\) se estima:

\[ R_{adj}^{2} = 1-\frac{\left(n-1\right)}{\left(n-p\right)}R_{\omega}^{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:

fit_svy <- svyglm(Income ~ Expenditure, 
                  design = diseno)

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

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.509 19005

y, para estimar el \(R_{adj}^{2}\) se requiere definir el diseño muestral pero incluyendo los q-weigthed (Pffeferman, 2011). A continuación, se muestra los pasos para encontrar los q-weigthed:

  • Ajustar un modelo de regresión a los pesos finales de la encuesta utilizando las variables predictoras en el modelo de regresión de interés.
fit_Nul <- lm(wk ~ 1, data = encuesta)
  • Obtener las predicciones de los pesos de la encuesta para cada caso como una función de las variables predictoras en el conjunto de datos
qw <- predict(fit_Nul)
  • Dividir los pesos finales de la encuesta por los valores predichos en el paso anterior:
encuesta %<>% mutate(wk1 = wk/qw)
  • Usar los nuevos pesos obtenidos para el ajuste de los modelos de regresión:
diseno_qwgt <- encuesta %>%
  as_survey_design(
    strata = Stratum,
    ids = PSU,
    weights = wk1,
    nest = T)

Ahora bien, una vez definido el diseño muestral con los nuevos pesos q-weigthed, se procede a calcular el \(R_{adj}^{2}\) como sigue:

n = sum(diseno_qwgt$variables$wk)
p<- 2
R2Adj = 1-( ( (n-1)/(n-p) )*R2 )
R2Adj             
##      variance    SE
## [1,]    0.491 19005

Como se puede observar, el \(R_{adj}^{2}\) es un poco más bajo que el \(R^{2}\) y cercanos al 50% que como se comentó anteriormente, dependiendo del contexto del problema se podrá concluir si es grande o pequeño.

Después de realizar la comparación entre las diferentes formas de estimar los coeficientes del modelo se opta por la metodología consolidadas en svyglm:

diseno_qwgt %<>% mutate(Age2 = Age^2) 
mod_svy <- svyglm( Income ~ Expenditure + Zone + Sex + Age2 ,
                       design = diseno_qwgt)
s1 <- summary(mod_svy)
s0 <- summary(modNul)

mod_svy

Stratified 1 - level Cluster Sampling design (with replacement) With (238) clusters. Called via srvyr Sampling variables: - ids: PSU - strata: Stratum - weights: wk1

Call: svyglm(formula = Income ~ Expenditure + Zone + Sex + Age2, design = diseno_qwgt)

Coefficients: (Intercept) Expenditure ZoneUrban SexMale Age2
62.18419 1.22548 63.46000 21.73256 0.00852

Degrees of Freedom: 2604 Total (i.e. Null); 115 Residual Null Deviance: 6.35e+08 Residual Deviance: 3.08e+08 AIC: 38300

stargazer(mod_svy, header = FALSE,single.row = T,
           title = "Modelo propuesto", 
           style = "ajps",  omit.stat=c("bic", "ll"))

Diagnósticos de los residuales

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 como sigue (Heeringa)

\[ r_{p_{i}} = \left(y_{i}-\mu_{i}\left(\hat{\beta}_{\omega}\right)\right)\sqrt{\frac{\omega_{i}}{V\left(\hat{\mu}_{i}\right)}} \] 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.

Otra definición que se debe tener en consideración para el análisis de los residuales es el de la matriz hat, la cual se estima como:

\[ H = W^{1/2}X\left(X'WX\right)^{-1}X'W^{1/2} \] donde,

\[ W = diag\left\{ \frac{\omega_{1}}{V\left(\mu_{1}\right)\left[g'\left(\mu_{1}\right)\right]^{2}},...,\frac{\omega_{n}}{V\left(\mu_{n}\right)\left[g'\left(\mu_{n}\right)\right]^{2}}\right\} \] \(W\) es una matriz diagonal de \(n\times n\) y \(g()\) es la función de enlace del modelo lineal generalizado.

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 removerlo de la nube de puntos este causa un cambio grande en el ajuste del modelo. Una observación importante para resaltar es 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:

Distancia de Cook’s: 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. Se calcula de la siguiente manera:

\[ c_{i}=\frac{w_{i}^{*}w_{i}e_{i}^{2}}{p\phi V\left(\hat{\mu}_{i}\right)\left(1-h_{ii}\right)^{2}}\boldsymbol{x}_{i}^{t}\left[\widehat{Var}\left(U_{w}\left(\hat{\boldsymbol{B}}_{w}\right)\right)\right]^{-1}\boldsymbol{x}_{i} \]

donde,

  • \(w_i^* =\) Pesos de la encuesta.
  • \(w_i\) Elementos por fuera de la diagonal de la matriz hat
  • \(e_i=\) residuales
  • \(p=\) número de parámetros del Modelo de regresión.
  • \(\phi =\) parámetro de dispersión en el glm
  • \(\widehat{Var}\left(U_{w}\left(\hat{\boldsymbol{B}}_{w}\right)\right) =\) estimación de varianza linealizada de la ecuación de puntuación, que se utiliza para pseudo MLE en modelos lineales generalizados ajustados a datos de encuestas de muestras complejas.

Luego del cálculo de la distancia de Cook’s para las observaciones de la muestra, se procede a calcular el siguiente estadístico de prueba para evaluar la importancia de la estadística:

\[ \frac{\left(df-p+1\right)\times c_{i}}{df} \doteq F_{\left(p,df-p\right)} \]

donde \(df=\) grados de liberta basados en el diseño. Por otro lado, la literatura como Tellez (2016), Heeringa considera a las observaciones influyentes cuando \(c_{i}\) sean mayores a 2 o 3.

\(D_fBeta_{(i)}\): Este estadístico mide el cambio en la estimación del vector de coeficientes de regresión cuando la i-ésima observación es eliminada. Se evalúa con la siguiente expresión:

\[ D_fBeta_{(i)} = \hat{B}-\hat{B}_{\left(i\right)}=\frac{\boldsymbol{A}^{-1}\boldsymbol{X}_{\left(i\right)}^{t}\hat{e}_{i}w_{i}}{1-h_{ii}} \]

Donde \(\boldsymbol{A} =\boldsymbol{X}^{t}\boldsymbol{WX}\) \(\hat{B}_{(i)}\) es el vector de parámetros estimados una vez se ha eliminado la i-ésima observación, \(h_{ii}\) es el correspondiente elemento de la diagonal de H y \(\hat{e}_i\) es el residual de la i-ésima observación.

Otra forma de reescribir este estadístico en términos de la matriz \(H\) es:

\[ D_fBetas_{\left(i\right)}=\frac{{c_{ji}e_{i}}\big/{\left(1-h_{ii}\right)}}{\sqrt{v\left(\hat{B}_{j}\right)}} \]

donde:

  • \(c_{ji}=\) es el ji-estimo elemento de \(\boldsymbol{A}^{-1}w_{i}^{2}\boldsymbol{X}_{\left(i\right)}\boldsymbol{X}_{\left(i\right)}^{t}\boldsymbol{A}^{-1}\)

  • El estimador de \(v\left(\hat{B}_{j}\right)\) basado en el Modelo se obtiene como: \(v_{m}\left(\hat{B}_{j}\right)=\hat{\sigma}\sum_{i=1}^{n}c_{ji}^{2}\) con \(\hat{\sigma}=\sum_{i\in s}w_{i}e^2/ \left( \hat{N} - p \right)\) y \(\hat{N} = \sum_{i \in s}w_{i}\)

  • La i-ésima observación es influyente para \(B_j\) si \(\mid D_{f}Betas_{\left(i\right)j}\mid\geq\frac{z}{\sqrt{n}}\) con \(z=\) 2 o 3

  • Como alternativa puede usar \(t_{0.025,n-p}/\sqrt(n)\) donde \(t_{0.025,n-p}\) es el percentil \(97.5\)

Estadístico \(D_{f}Fits_{\left(i\right)}\): Este estadístico mide el cambio en el ajuste del modelo cuando se elimina el registro i-ésimo. Se calcula de la siguiente manera:

\[ D_{f}Fits_{\left(i\right)}= \frac{h_{ii}e_{i}\big/\left(1-h_{ii}\right)}{\sqrt{v\left(\hat{\beta}_{j}\right)}} \]

Donde, \(\sqrt{v\left(\hat{\beta}_{j}\right)}\) puede ser aproximada por el diseño o el Modelo. La i-ésima observación se considera influyente en el ajuste del Modelo si \(\mid DfFits\left(i\right)\mid\geq z\sqrt{\frac{p}{n}}\) con \(z =\) 2 o 3

Por otro lado, un análisis que es de vital importancia en el ajuste de modelos de regresión más específicamente en el análisis de residuales es el de varianza constante en los errores. La principal consecuencia de no tener en cuenta la violación de este supuesto es que los estimadores pierden eficiencia. 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. Como consecuencia de lo anterior, los intervalos de confianza serán más amplios y las pruebas t y F darán resultados imprecisos (Tellez, 2016).

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 un patrón (funciones cuadráticas, cúbicas, logarítmicas, etc), se puede decir que la varianza de los errores no es constante.

Otro supuesto que se debe revisar en los errores al momento de realizar ajustes es la normalidad en lo errores. Una forma muy común para hacer dicha evaluación 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.

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. Primero, se realizará el análisis de la normalidad en los errores de manera gráfica como se muestra a continuación:

par(mfrow = c(2,2))
plot(mod_svy)

Como se puedo observar en el QQplot, hay evidencia gráfica de que los errores no se distribuyen según una distribución normal.

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(mod_svy)$stdresids)
diseno_qwgt$variables %<>% mutate(stdresids = stdresids)

Podemos hacer el análisis de normalidad también por medio del histograma de los residuales estandarizados como sigue:

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 = "")

y como se puede observar gráficamente los errores no siguen una distribución normal.

Por otro lado, el otro análisis que se realiza de manera gráfica es el de varianzas constantes el cual se realizará a continuación:

Primero, agreguemos las predicciones a la base de datos para poder realizar las gráficas.

library(patchwork)
diseno_qwgt$variables %<>% 
  mutate(pred = predict(mod_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, ambas muestran tendencias y no un comportamiento aleatorio. Por lo anterior, se puede decir que las varianzas no son constantes.

Otros de os análisis a realizar es revisar si existen datos influyentes en la base de datos. Para ejemplificar los conceptos definidos, se seguirán 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’s usando la función svyCooksDdel paquete svydiags como sigue:

library(svydiags) 
d_cook = data.frame(
   cook = svyCooksD(mod_svy), 
     id = 1:length(svyCooksD(mod_svy)))

table(d_cook$cook>3)


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 lo que, podemos decir que no existen observaciones influyentes.

Ahora bien, se desea observar si hay observaciones influyentes pero utilizando \(D_{f}Betas_{\left(i\right)j}\) se realiza con la función svydfbetas como se muestra a continuación:

d_dfbetas = data.frame(t(svydfbetas(mod_svy)$Dfbetas))
colnames(d_dfbetas) <- paste0("Beta_", 1:5)
d_dfbetas %>% slice(1:10L)
Beta_1 Beta_2 Beta_3 Beta_4 Beta_5
0.0006 -2e-04 0.0021 -0.0045 -0.0077
-0.0006 -1e-04 0.0014 0.0026 -0.0031
-0.0009 -1e-04 0.0009 0.0022 0.0008
-0.0004 -1e-04 0.0012 -0.0031 0.0007
-0.0009 0e+00 0.0008 0.0021 0.0014
0.0009 6e-04 -0.0036 -0.0063 0.0098
0.0027 4e-04 -0.0031 -0.0076 -0.0028
0.0011 3e-04 -0.0028 0.0077 -0.0043
0.0030 4e-04 -0.0030 -0.0078 -0.0051
-0.0003 4e-04 0.0012 -0.0037 -0.0040

Una vez calculado los \(D_{f}Betas_{\left(i\right)j}\) se procede a acomodar la salida con 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(mod_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_1 0.2781 Si
890 Beta_2 -0.2593 Si
891 Beta_1 0.2559 Si
889 Beta_2 -0.2537 Si
891 Beta_2 -0.2491 Si
890 Beta_1 0.2456 Si
2311 Beta_5 0.2056 Si
889 Beta_5 -0.1993 Si
890 Beta_5 -0.1788 Si
890 Beta_4 0.1616 Si

Como se pudo observar en la salida anterior hay varias observaciones que resultan influyentes dado el criterio del \(D_{f}Betas_{\left(i\right)j}\). A continuación, y de manera ilustrativa, se grafican los \(D_{f}Betas_{\left(i\right)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 ahora es detectar observaciones influyentes pero considerando ahora la estadística \(D_{f}Fits_{\left(i\right)}\), se utiliza la función svydffits y se siguen los mismos pasos mostrados para el estadístico \(D_{f}Betas_{\left(i\right)j}\):

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

cutoff <- svydffits(mod_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_{\left(i\right)}\), las cuales se muestran en rojo en el gráfico.

Un último acercamiento que se trabajará en este texto para la detección de datos influyentes está encaminado al uso de la matriz H. En este sentido, la matriz asociada al Estimador de Pseudo Máxima Verosimilitud (PMLE) de \(\hat{\boldsymbol{B}}\) es \(\boldsymbol{H}=\boldsymbol{XA}^{-1}\boldsymbol{X}^{-t}\boldsymbol{W}\) cuya diagonal esta dado por \(h_{ii} = \boldsymbol{x_{i}^tA}^{-1}\boldsymbol{x_{i}}^{-t}w_{i}\). Utilizando la matriz H, una observación puede ser grande y, como resultado, influir en las predicciones, cuando un \(x_i\) es considerablemente diferente del promedio ponderado \(\bar{x}_w=\sum_{i\in s}w_{i}\boldsymbol{x_{i}}\big/\sum_{i\in s}w_i\). Según (Tellez, 2016) una observación es considerada grande si es mayor a tres veces el promedio de los \(h_{ii}\). A continuación, se muestra el procedimiento en R cuya función a utilizar es svyhat:

vec_hat <- svyhat(mod_svy, doplot = FALSE)
d_hat = data.frame(hat = vec_hat, id = 1:length(vec_hat))
d_hat %<>% mutate(C_cutoff = ifelse(hat > (3 * mean(hat)),"Si", "No"))

ggplot(d_hat, aes(y = hat, x = id)) +
  geom_point(aes(col = C_cutoff)) + 
  geom_hline(yintercept = (3 * mean(d_hat$hat))) +
  scale_color_manual(
    values = c("Si" = "red", "No" = "black"))+
  theme_cepal()

Dado que esta última técnica es empírica, se puede observar en el gráfico anterior que hay varias observaciones posiblemente influyentes en el conjunto de datos de la muestra de hogares.