1.1 Linealidad

Para la regresión lineal múltiple, un supuesto importante es que existe una relación lineal entre la variable dependiente e independiente.

Podemos saber si se cumple este supuesto a partir de un gráfico de dispersión de datos, que relacione la variable dependiente y la independiente, y verificar de manera “intuitiva” si la tendencia de esta relación se puede describir por una línea recta. El paquete performace nos permite hacer esto con su función check_model indicando en el argumento check = "ncv

check_model(model1, check = c("ncv", "linearity"))

En este caso podemos notar claramente que no hay una relación lineal. En caso de que no sea claro, una forma numérica para chequear este supuesto es que el valor promedio de los residuos sea cero. Si esto no es así los residuos están sesgados sistemáticamente, entonces probable que el modelo no sea realmente lineal, y esta desviación de los residuos deba corregirse re-especificando (medir de otra manera la variable) algún término de la ecuación de regresión al cuadrado o al cubo. Un modo que se ocupa para testear la necesidad de esta re-especificación es el test RESET de Ramsey que indica que:

\(H_o\) cuando el modelo tiene algún término al cuadrado o cubo, la media de los residuos es cero. Es decir, que si no podemos rechazar \(H_0\), nuestro modelo está bien especificado (es decir, es lineal). Como en este caso es claro que se necesita respecificar el predictor elevamos edad al cuadrado. Imaginemos que es el caso

✔️ La línea de referencia es plana y horizontal

1.2 Homocedasticidad

Homoce ¿qué? Sí, homocedasticidad. Este concepto indica que los residuos se distribuyen de forma homogénea (por eso el sufijo homo). Como ya podrás haber notado este supuesto se vincula con el de linealidad

Al igual que la linealidad también puede comprobarse con un gráfico de dispersión entre la variable dependiente (\(Y\)) e independiente (\(X\)), donde podamos ver de manera clara la recta de regresión estimada y la distribución de los residuos. Aceptaremos el supuesto de homocedasticidad si la variación de los residuos es homogénea, es decir, no veremos un patrón claro y más bien se distribuirán de forma aleatoria. De manera gráfica veremos una nube de puntos que tiene una forma similar en todo el rango de las observaciones de la variable independiente.

Para comprobar el supuesto de homocedasticidad de manera más certera utilizaremos la prueba Breusch-Pagan Godfrey cuya hipótesis nula indica que

\(H_0\): La varianza de los residuos del modelo de regresión no es constante (heterocedasticidad)

En este caso, buscaremos que rechazar la \(H_0\). Esto implicaría que “en suma y resta” si bien hay residuos, estos tienen una variación homogénea en todos los tramos de la relación de la variable dependiente con la independiente.

A partir de la función check_heteroscedasticity verificaremos qué ocurre con la hipótesis nula

check_heteroscedasticity(model1)

OK: Error variance appears to be homoscedastic (p = 0.494).

✔️ La varianza es homocedástica

1.3 Normalidad de residuos

Además de la linealidad (media 0), la homocedasticidad (varianza mínima y constante), el método de estimación de la regresión lineal (OLS o MCO en español) requiere asegurar una distribución normal de los residuos pues en caso contrario el modelo no es consistente a través de las variables y observaciones (esto significa que los errores no son aleteatorios).

Al igual que con los otros supuestos, la normalidad de los residuos se puede evaluar con métodos numéricos con pruebas que ya conocemos de otros cursos como la prueba de Shapiro-Wilk y Kolmogorov-Smirnov

A partir de la función check_normality utilizaremos la prueba Shapiro-Wilk para ver qué ocurre con la hipótesis nula a

check_normality(model1)

Warning: Non-normality of residuals detected (p < .001).

⚠️ Los residuos no son normales

1.4 Independencia

Evidentemente si los residuos no siguen una distribución normal, es probable que estos no sean independientes entre sí. Esto significa que buscaremos que los errores asociados a nuestro modelo de regresión sean independientes. Para saber si se cumple ese criterio se utiliza la prueba de Durbin-Watson, donde la \(H_0\) supone que los residuos son independientes

A partir de la función check_autocorrelation utilizaremos la prueba Durbin-Watson para ver qué ocurre con la hipótesis nula a

check_autocorrelation(model1)

Warning: Autocorrelated residuals detected (p = 0.010).

⚠️ Hay correlación entre los residuos

En síntesis sabemos la regresión lineal requiere de una relación lineal entre sus variables independientes y dependiente. Para ello no solo es importante chequear la distribución de los residuos, sino dos posibilidades que pueden tendenciar esa relación lineal: como casos influyentes en la muestra y predictores que están altamente relacionados. Revisaremos la última de estas

1.5 Multicolinealidad

La multicolinealidad es la relación de dependencia lineal fuerte entre más de dos predictores de un modelo.

El problema que produce es que será difícil cuantificar con exactitud el efecto de cada predictor sobre la variable dependiente, precisamente pues puede ocurrir que el efecto que una variable predictora tenga sobre el fenómeno que se busca estudiar dependa del valor de otra variable del modelo.

Para la regresión múltiple esto implica un problema pues suponemos que podemos “controlar” por el otro valor de la variable.

Esta relación endógena entre predictores la podemos examinar ante la existencia de altas correlaciones (lineales) entre variables. La aproximación numérica más utilizada es el VIF (factor de inflación de varianza) que indica hasta que punto la varianza de los coeficientes de regresión se debe a la colinealidad (o dependencia) entre otras variables del modelo.

Para evaluar esto ocuparemos el comando check_collinearity(). Como podemos ver en el gráfico, todos los valores son menores a 5 (como recomienda el paquete).

plot(check_collinearity(model1))

Ahora bien, dado que sabemos que las correlaciones en ciencias sociales nunca son tan altas, un criterio que se ocupa en nuestras disciplinas para evaluar multicolinealidad es es evitar valores del VIF mayores a 2.5.

check_collinearity(model1)

2.1 Transformar predictor al cuadrado o cubo

Si estamos ante problemas de linealidad, indicamos que el término cuadrático o al cubo de algún predictor produce que la media de los residuos sea 0. Por lo general, por su distribución, esta variable es edad.

movid_proc <- movid_proc %>% mutate(edad2 = (edad)^2)

2.2 Logaritmizar variable dependiente

Otro caso similar ocurre con ingresos, donde lo que se hace frecuentemente es transformar ingresos en un logaritmo de ingresos (log(ingresos)).

movid_proc %>% select(ingreso, tingreso) %$%
print(dfSummary(., headings = FALSE, method = "render", varnumbers = F, lang = "es"))

movid_proc <- movid_proc %>%
  mutate(tingreso = case_when(tingreso == "Menos de $200 mil pesos" ~ 200000,
                              tingreso == "Entre $200 y 350 mil pesos" ~ 275000,
                              tingreso == "Entre $351 y 500 mil pesos" ~ 425500,
                              tingreso == "Entre 501 y 800 mil pesos" ~ 650500,
                              tingreso == "Entre 801 mil y 1 millón 200 mil pesos" ~ 1000500,
                              tingreso == "Entre 1 millón 201 mil y 2 millones de pesos" ~ 1600500,
                              tingreso == "Entre 2 millones y 5 millones de pesos" ~ 3500000,
                              tingreso == "Más de 5 millones de pesos" ~ 5000000), #Paso 1
         ingreso = if_else(is.na(ingreso), tingreso, ingreso))

movid_proc %>% select(ingreso, tingreso) %$%
print(dfSummary(., headings = FALSE, method = "render", varnumbers = F, lang = "es")) #Paso 3

movid_proc <- movid_proc %>%
  mutate(log_ing = log(ingreso), #Log ingresos
         ing_per = ingreso/tamanohogar, #Ingreso percapita
    quintiles = dplyr::ntile(ing_per,
                              n = 5)) # n de categorias, para quintiles usamos 5

2.3 Dicotomizar variable dependiente

Es recurrente que en las encuestas sociales nos encontremos con preguntas con Escala Likert. Sin embargo, muchas veces estas no tienen una distribución normal, y más bien las respuestas están concentradas en algunas categorías de referencia.

Si bien no hay criterios canónicos para trabajar esas variables, tiene sentido indicar que más bien no representan un constructo con 5 o más niveles, sino que probablemente de 2. Re-especificar la variable como dicotómica no solo ayudará a trabajar de manera más realista el constructo, sino que facilitará las interpretaciones que queramos hacer de nuestro modelo.

Ocuparemos dos criterios para la dicotomización:

1. Medias: se ocupará como criterio discriminante la media de la variable (donde 1 puede ser los valores mayores a la media, y 0 los inferiores).

2. Mediana: la más frecuente en medidas ordinales como las escalas Likert es cuando el 50% de los casos se concentra en unas pocas categoría de respuesta (eg, “Muy de acuerdo” y “De acuerdo” serán 1 y el resto 0).

En el caso de la variable fatiga que indica “A medida que ha avanzado la crisis sanitaria, me siento cada vez más desmotivado para seguir las medidas de protección recomendadas”, recodificaremos a aquellos como 1 a quiénes asienten a esta frase (“Muy de acuerdo” y “De acuerdo”)

movid_proc <- movid_proc %>%
  mutate(fatigadummy = case_when(fatiga %in% c(5,4) ~ 1,
                                 fatiga %in% c(3,2,1) ~ 0, TRUE ~ NA_real_))

2.4 Errores estándares robustos

En caso de que estemos ante problemas de heterocedasticidad debemos re-estimar nuestro modelo considerando errores estándares robustos

model_robust<- lmtest::coeftest(model1, vcov=sandwich::vcovHC(model1))

2.5 Creación de índices

En el tutorial de dplyr ya habíamos abordado este punto. También está desarrollado de manera extensa en el práctico N°1. Por ello, solo lo cacularemos para aquellos ítems que mostraron posibilidad de ser colineales. Los pasos son

1. Correlacionar para verificar que estamos ante la presencia de ítems que podrían estar midiendo un constructo común. Como vemos en la siguiente tabla, las correlaciones son altas

movid_proc %>% select(starts_with("c2")) %>%
  mutate_all(~as.numeric(.)) %>%
sjPlot::tab_corr(., triangle = "lower")

2. Construcción de índice: este puede ser sumativo o promedio. Esto dependerá de la escala de los ítems y del sentido del constructo final que queremos utilizar. En nuestro caso crearemos una índice sumativa de salud mental

movid_proc <- movid_proc %>%
  mutate_at(vars(starts_with("c2")),~as.numeric(.)) %>%
  rowwise() %>%
  mutate(salud_mental = sum(c2_1,c2_2,c2_3,c2_4, na.rm = T))

2.6 Eliminar casos influyentes

En caso de que estemos ante la prescencia de casos influyentes debemos seguir los siguientes pasos

n<- nobs(model1) #n de observaciones
k<- length(coef(model1)) # n de parametros
dcook<- 4/(n-k-1) #Punto de corte

# Datos donde se filtran los valores sobre el punto de corte

movid_proc_so <- broom::augment_columns(model1,data = movid_proc) %>% filter(.cooksd<dcook)

Una vez finalizada las transformaciones calcular dos nuevos modelos

model1_fit <- lm(as.numeric(fatiga) ~ salud_mental +
               trabaja + sexo + edad + ing_per,
             data = movid_proc)

model2 <- lm(as.numeric(fatiga) ~ salud_mental +
               trabaja + sexo + edad2 + ing_per,
             data = movid_proc)

model3 <- glm(fatigadummy ~ salud_mental +
               trabaja + sexo + edad2 + ing_per, family = "binomial",
             data = movid_proc)

Luego, podemos hacer un chequeo general de diagnósticos de robustez con check_model, pero ahora indicando que queremos evaluar todos los indicadores posibles

check_model(model1_fit, check = c("vif","normality", "linearity", "ncv", "homogeneity"))

check_model(model2, check = c("vif",  "normality", "linearity", "ncv", "homogeneity"))

check_model(model3, check = c("vif",  "homogeneity"))

Existen otros diagnósticos posibles para abordar en nuestros modelos, todo con el propósito de mejorar la calidad de estos. Uno de ellos, y que no abordaremos por su extensión, es la ausencia de posibles interacciones entre las variables que no han sido modeladas (Fox & Weisberg 2018). En caso de su interés pueden revisar esto y ver su aplicación simple en R en el siguiente link.