Vamos a profundizar en el uso de la regresión lineal (RL) con un ejemplo más sofisticado; también de la ENBIARE. Vamos a ver la asociación entre prestamós (VI/ Tr) y depresión (VD / outcome)

Limpieza de datos

Empezamos por crear el ID

library(tidyverse)
library(skimr)
library(ggdist)
library(sjPlot)
library(sjlabelled)
library(sjmisc)
library(broom)
library(effectsize)
library(performance)
library(see)



library(readr)
TSDEM <- read_csv("Datos/TSDEM.csv")

TSDEM1 <- TSDEM %>% 
   mutate(HOGAR= as.character(HOGAR)) %>% 
  mutate(id_personal= paste(UPM_DIS, VIV_SEL, HOGAR, N_REN, sep=".")) # ID


TENBIARE <- read_csv("Datos/TENBIARE.csv")

TENBIARE1 <- TENBIARE %>% 
   mutate(HOGAR= as.character(HOGAR)) %>% 
  mutate(id_personal= paste(UPM_DIS, VIV_SEL, HOGAR, N_REN, sep=".")) # ID

enbiare <- TENBIARE1 %>% # nuevo objeto bd
  left_join(TSDEM1, by= "id_personal")

Depresión

La variable dependiente es depresión medida con CESD-7. Es una vieja conocida por la encovid

enbiare <- enbiare %>% 
  mutate(PD2_6= case_when(  
    PD2_6== 0 ~ 3,
    PD2_6== 1 ~ 2,
    PD2_6== 2 ~ 1,
    PD2_6== 3 ~ 0)) 

enbiare <- enbiare %>% 
  mutate(d_total = PD2_1 + PD2_2 + PD2_3 + PD2_4 + PD2_5 + PD2_6 + PD2_7)

summary(enbiare$d_total)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   1.000   3.000   4.373   7.000  21.000

Y la gráfica como histograma; señalo la media en azul y el punto de corte en negro (9). Noten que no tiene una distribución normal porque hay un sesgo hacia los valores bajos. En otras palabras, la mayoría de las personas tienen cero o poca depresión y el valor medio (4.3) no es el más común. Más adelante esto será un problema porque un supuesta de la RL. También noten que el rango es restringido: números enteros que van de 0 a 21. Estos son los límites y la forma de la distribución que buscamos entender/ predecir con la RL.

enbiare %>%
  ggplot(aes(x= d_total))+
  geom_histogram(binwidth = 1, fill= "darkorange", color="white")+
  geom_vline(xintercept = mean(enbiare$d_total), color= "blue")+
  geom_vline(xintercept = 9, color= "black")+ # punto de corte
   theme_minimal()

Veamos de una vez la prevalencia: 16.4%

enbiare <- enbiare %>% 
  mutate(dep_level= if_else(d_total >= 9, 1, 0))

enbiare <- enbiare %>% 
  mutate(dep_level2= ordered(dep_level, levels=c("0", "1"), labels= c("Sin Síntomas", "Con Síntomas")))

table(enbiare$dep_level2, useNA = "ifany")

## 
## Sin Síntomas Con Síntomas 
##        26050         5116

prop.table(table(enbiare$dep_level2, useNA = "ifany"))*100

## 
## Sin Síntomas Con Síntomas 
##     83.58468     16.41532

Préstamos

Esta es nuestra principal variable independiente. Siguiendo el libro de texto, esta es la que consideraremos el “Tratamiento” aún cuando no estemos en un contexto experimental. El sentido de esto es que consideramos que cambios en “préstamos” se traducen en cambios en “depresión”. Asumiremos que hay causalidad aunque sabemos, a priori, que nuestro diseño no nos permite concluir causalidad. Sólo aspiramos a hacer comparaciones cuidadosas que nos muestren que hay asociación entre ciertas subpoblaciones o con ajustes estadísticos necesarios.

Veamos cómo se construyen “préstamos”.

PreguntaPrestamos

Tenemos varias maneras de usar estos datos. Cada uno es una variación téorica y cada uno tiene ciertos supuestos. Me quedo con alimentos que es la más importante (para mi)? Excluyo colegiaturas porque no todo mundo tiene hijos? y medicinas si no me enfermé? Renta si soy dueño?

Voy a sumarlas. Mi hipótesis es que mientras más necesidades tengo para pedir préstamos, más depresión tendré. Si tengo X y además tengo Y y luego me cayó Z, estoy más deprimido que quien solo tiene Z. Es decir, un puntaje de préstamos más alto se traduce en mayores puntos en depresión. El supuesto es que todos los tipos de préstamo valen lo mismo. colegiaturas = agua. Se vale? Estoy dispuesto a vivir con ello.

Vamos a combinarla. Vean las diferencias de medias en el skim. El préstamo más común es alimentos y medicinas. Los menos comunes son colegiaturas y rentas. Vean las proporciones: 62% no pide préstamos para nada. 38% tuvo que pedir prestado para pagar cosas básicas! De ahí vemos cómo desciende. 10% para 1; 8% para 2; 7% para 3… 1.5% para las 6.

table(enbiare$PF1_1, useNA = "ifany") # alimentos

## 
##     1     2 
##  9083 22083

table(enbiare$PF1_5, useNA = "ifany") # colegiatura

## 
##     1     2 
##  2002 29164

enbiare <- enbiare %>% 
  mutate_at(
            vars(starts_with("PF1_")),
            funs(case_when(
                           .== 1 ~ 1,
                           .== 2 ~ 0)))

enbiare %>%
  select(starts_with("PF1_")) %>% 
  skim()

Data summary
Name	Piped data
Number of rows	31166
Number of columns	6
_______________________
Column type frequency:
numeric	6
________________________
Group variables	None

Variable type: numeric

skim_variable	complete_rate	mean	sd	p75	p100	hist
PF1_1	1	0.29	0.45	1	1	▇▁▁▁▃
PF1_2	1	0.09	0.28	0	1	▇▁▁▁▁
PF1_3	1	0.14	0.35	0	1	▇▁▁▁▁
PF1_4	1	0.20	0.40	0	1	▇▁▁▁▂
PF1_5	1	0.06	0.25	0	1	▇▁▁▁▁
PF1_6	1	0.21	0.41	0	1	▇▁▁▁▂

enbiare <- enbiare %>% 
  mutate(prest_cont = PF1_1 + PF1_2 + PF1_3 + PF1_4 + PF1_5 + PF1_6)

table(enbiare$prest_cont, useNA = "ifany")

## 
##     0     1     2     3     4     5     6 
## 19609  3385  2563  2156  1927  1064   462

prop.table(table(enbiare$prest_cont, useNA = "ifany"))*100

## 
##         0         1         2         3         4         5         6 
## 62.917923 10.861195  8.223705  6.917795  6.183020  3.413977  1.482385

Vamos a ver sus asociaciones de forma descriptiva. ¿qué patrón esperamos ver según nuestra hipótesis? A mayor puntaje, mayor depresión en promedio. Usemos las medias y las medianas.

Las medias muestran un crecimiento muy claro. Y la mediana es más contundente. Se duplica al pasar de 0 a 1.

enbiare %>% 
  group_by(prest_cont) %>% 
  summarise(medias= mean(d_total), medianas= median(d_total))

## # A tibble: 7 × 3
##   prest_cont medias medianas
##        <dbl>  <dbl>    <dbl>
## 1          0   3.65        2
## 2          1   4.91        4
## 3          2   5.31        4
## 4          3   5.72        5
## 5          4   6.37        5
## 6          5   6.38        5
## 7          6   6.73        5

Veamos la relación gráficamente. El brinco de 0 a 2 es fuerte. Después es menos claro.

enbiare %>%
  ggplot( aes(x= as.character(prest_cont), y=d_total)) +
    geom_boxplot()

Otra versión para practicar las capas de ggplot. Fíjense en las capas que añadí, el stat (media de d_total) y que cambié el geom para ver las observaciones en cada categoría.

enbiare %>%
  ggplot( aes(x= as.character(prest_cont), y=d_total)) +
  geom_boxplot()+
    geom_jitter(color="gray", size=0.4, alpha=0.3) +
  stat_summary(fun.y="mean", color= "red")+
  theme_classic()

## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: Removed 7 rows containing missing values (geom_segment).

Ahora voy a hacer una versión un poco más sofisticada para que no haya dudas de los que estmamos haciendo. Noten que los box plots están en horizontal. Lean de abajo para arriba. La gráfica se va estirando y moviendo hacia la derecha (poquito).

enbiare %>%
  ggplot(aes(x= as.factor(prest_cont), y=d_total)) +
  stat_halfeye(adjust = 1,  justification= -.2, .width = 0, point_color= NA)+
  geom_boxplot(width=.2, alpha=.5)+
  stat_summary(fun.y="mean", color= "red")+
  coord_flip()+
  theme_classic()

## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: Removed 7 rows containing missing values (geom_segment).

Otra representación para evidenciar que el 0 es sustancialmente distinto a las demás. Vean cómo la base de esas botellas (distribuciones) cambia.

enbiare %>%
  ggplot( aes(x= as.character(prest_cont), y=d_total)) +
  geom_violin()+
  stat_summary(fun.y="mean", color= "red")+
  theme_classic()

## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: Removed 7 rows containing missing values (geom_segment).

Vamos a hacerla dicotómica, por si se ocupa. La diferencia en medias y medianas entre algún préstamo o no es de 2 puntotes.

enbiare <- enbiare %>% 
  mutate(prest_dico= if_else(prest_cont > 0, 1, 0))

enbiare <- enbiare %>% 
  mutate(prest_dico_ord= ordered(prest_dico, levels = c("0", "1"), 
                                labels=c("No", "Sí")))

table(enbiare$prest_dico_ord, useNA = "ifany")

## 
##    No    Sí 
## 19609 11557

enbiare %>% 
  group_by(prest_dico_ord) %>% 
  summarise(medias= mean(d_total), medianas= median(d_total))

## # A tibble: 2 × 3
##   prest_dico_ord medias medianas
##   <ord>           <dbl>    <dbl>
## 1 No               3.65        2
## 2 Sí               5.60        4

Ya para ver la diferencia gráfica en términos más absolutos

enbiare %>%
  ggplot( aes(x= prest_dico_ord, y=d_total)) +
    geom_boxplot()+
  stat_summary(fun.y="mean", color= "red")+
  theme_classic()

## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: Removed 2 rows containing missing values (geom_segment).

enbiare %>%
  ggplot( aes(x= prest_dico_ord, y=d_total)) +
  geom_violin()+
  stat_summary(fun.y="mean", color= "red")+
  theme_classic()

## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: Removed 2 rows containing missing values (geom_segment).

Controles

Vamos a construir las vars independientes con las que vamos a ajustar nuestro modelo.

¿Cómo elegimos los controles? Este diagrama es lo que tenemos en mente.

Son variables que preceden temporalmente al Tx y que se asocian, de forma independiente, tanto con Tx como con la VD; especialmente con la VD.

Sexo

Vean cómo voy nombrando las nuevas variables. Me interesa tener muy claro quién es el 0 y quién el 1.

table(enbiare$SEXO, useNA = "ifany")

## 
##     1     2 
## 14255 16911

enbiare <- enbiare %>% 
  mutate(mujer= case_when(
    SEXO== 1 ~ 0,
    SEXO== 2 ~ 1))

table(enbiare$mujer, useNA = "ifany")

## 
##     0     1 
## 14255 16911

enbiare <- enbiare %>% 
  mutate(mujer_ord= ordered(mujer, levels = c("0", "1" ), labels=c("Hombre", "Mujer")))

table(enbiare$mujer_ord, useNA = "ifany")

## 
## Hombre  Mujer 
##  14255  16911

Queremos conocer descriptivamente las relaciones entre todas mis vars independientes y las dependientes

La media de depresión de las mujeres es considerablemente más alta

t.test(d_total ~ mujer_ord, data = enbiare)

## 
##  Welch Two Sample t-test
## 
## data:  d_total by mujer_ord
## t = -26.426, df = 31088, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Hombre and group Mujer is not equal to 0
## 95 percent confidence interval:
##  -1.410288 -1.215530
## sample estimates:
## mean in group Hombre  mean in group Mujer 
##             3.660540             4.973449

enbiare %>%
  ggplot( aes(x=d_total, fill=mujer_ord)) +
    geom_bar(position = "dodge2")+ 
    scale_fill_manual(values=c("blue", "red"))

enbiare %>%
  ggplot( aes(x= mujer_ord, y=d_total)) +
    geom_boxplot()+
  stat_summary(fun.y="mean", color= "red")+
  theme_classic()

## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: Removed 2 rows containing missing values (geom_segment).

Ya hora vemos con prestamos. Las mujeres también piden más prestamos. Pero es menos claro.

t.test(prest_cont ~ mujer_ord, data = enbiare)

## 
##  Welch Two Sample t-test
## 
## data:  prest_cont by mujer_ord
## t = -7.5065, df = 30751, p-value = 6.236e-14
## alternative hypothesis: true difference in means between group Hombre and group Mujer is not equal to 0
## 95 percent confidence interval:
##  -0.16793117 -0.09839122
## sample estimates:
## mean in group Hombre  mean in group Mujer 
##             0.915328             1.048489

Aquí no se ve claramente porque hay más mujeres en la muestra

enbiare %>%
  ggplot( aes(x=prest_cont, fill=mujer_ord)) +
    geom_bar(position = "dodge2")+ 
    scale_fill_manual(values=c("blue", "red"))

Atención a la variable en la Y. Noten cómo son medias cajas. Alguna idea de la razón?

enbiare %>%
  ggplot( aes(x= mujer_ord, y=prest_cont)) +
    geom_boxplot()+
  stat_summary(fun.y="mean", color= "red")+
  theme_classic()

## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: Removed 2 rows containing missing values (geom_segment).

Educación

Vamos a reducir los niveles a 3 para que sea más manejable.

table(enbiare$NIVEL)

## 
##   00   01   02   03   04   05   06   07   08   09   10   99 
##  504   28 7126 8926  609   99 6264  729 6172   91  557   61

class(enbiare$NIVEL) # Tengo que usar comillas

## [1] "character"

enbiare <- enbiare %>% 
  mutate(educ= case_when(
    NIVEL== "00" ~ 0,
    NIVEL== "01" ~ 0, 
    NIVEL== "02" ~ 0, # Primaria
    NIVEL== "03" ~ 1,
    NIVEL== "04" ~ 1,
    NIVEL== "05" ~ 1,
    NIVEL== "06" ~ 1, # Prepa
    NIVEL== "07" ~ 2,
    NIVEL== "08" ~ 2,
    NIVEL== "09" ~ 2,
    NIVEL== "10" ~ 2, # Superior
    NIVEL== "99" ~ NA_real_))

enbiare <- enbiare %>% 
  mutate(educ_ord= ordered(educ, levels = c("0", "1", "2"), labels=c("Primaria", "Preparatoria", "Superior")))

table(enbiare$educ_ord, useNA = "ifany")

## 
##     Primaria Preparatoria     Superior         <NA> 
##         7658        15898         7549           61

Vamos a comparar. Primero medias por educación.

enbiare %>% 
  select(educ_ord, d_total) %>% 
  group_by(educ_ord) %>%
  summarize(media= mean(d_total), medianas= median(d_total))

## # A tibble: 4 × 3
##   educ_ord     media medianas
##   <ord>        <dbl>    <dbl>
## 1 Primaria      5.17        4
## 2 Preparatoria  4.32        3
## 3 Superior      3.68        3
## 4 <NA>          4.25        3

Y nuestra vieja amiga la anova. Sí tenemos diferencias

anova_educ <- aov(d_total ~ educ_ord, data = enbiare)

summary(anova_educ)

##                Df Sum Sq Mean Sq F value Pr(>F)    
## educ_ord        2   8545    4272     214 <2e-16 ***
## Residuals   31102 621031      20                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 61 observations deleted due to missingness

TukeyHSD(anova_educ)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = d_total ~ educ_ord, data = enbiare)
## 
## $educ_ord
##                             diff        lwr        upr p adj
## Preparatoria-Primaria -0.8464372 -0.9921127 -0.7007617     0
## Superior-Primaria     -1.4915499 -1.6614071 -1.3216926     0
## Superior-Preparatoria -0.6451127 -0.7914963 -0.4987291     0

La gráfica que ya trae el paquete:

plot(TukeyHSD(anova_educ))

Ahora comparemos con prestamos. Menos clara la diferencia.

enbiare %>% 
  select(educ_ord, prest_cont) %>% 
  group_by(educ_ord) %>%
  summarize(media= mean(prest_cont), medianas= median(prest_cont))

## # A tibble: 4 × 3
##   educ_ord     media medianas
##   <ord>        <dbl>    <dbl>
## 1 Primaria     1.10         0
## 2 Preparatoria 1.09         0
## 3 Superior     0.649        0
## 4 <NA>         0.984        0

Sólo hay dif con superior

anova_educ2 <- aov(prest_cont ~ educ_ord, data = enbiare)

summary(anova_educ2)

##                Df Sum Sq Mean Sq F value Pr(>F)    
## educ_ord        2   1145   572.7   236.2 <2e-16 ***
## Residuals   31102  75422     2.4                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 61 observations deleted due to missingness

TukeyHSD(anova_educ2)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = prest_cont ~ educ_ord, data = enbiare)
## 
## $educ_ord
##                               diff         lwr         upr     p adj
## Preparatoria-Primaria -0.009737777 -0.06050442  0.04102887 0.8946021
## Superior-Primaria     -0.454073155 -0.51326692 -0.39487939 0.0000000
## Superior-Preparatoria -0.444335378 -0.49534878 -0.39332198 0.0000000

Vean cómo cruza el 0 cuando se compara Prepa con Primaria

plot(TukeyHSD(anova_educ2))

Desempleo

Vamos a otras variables más sustantivas… La E.15. Desempleo. Dice ”

En los últimos 12 meses, ¿se quedó sin empleo o tuvo que cerrar un negocio propio?

+====================================================================================+ +————————————————————————————+

table(enbiare$PE15)

## 
##     1     2     3 
##  4300  2127 24739

enbiare <- enbiare %>% 
  mutate(desempleo= case_when(
    PE15== 1 ~ 1, # sí
    PE15== 2 ~ 0, # Sí, lo recuperó
    PE15== 3 ~ 0)) # No

table(enbiare$desempleo, useNA = "ifany")

## 
##     0     1 
## 26866  4300

enbiare <- enbiare %>% 
  mutate(desempleo_ord= ordered(desempleo, levels = c("0", "1"), 
                                labels=c("No", "Sí")))

table(enbiare$desempleo_ord, useNA = "ifany")

## 
##    No    Sí 
## 26866  4300

Vamos a comparar. Sí parece que desempleo afecta a la depresión

t.test(d_total ~ desempleo_ord, data = enbiare)

## 
##  Welch Two Sample t-test
## 
## data:  d_total by desempleo_ord
## t = -14.868, df = 5474.5, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group No and group Sí is not equal to 0
## 95 percent confidence interval:
##  -1.335112 -1.024041
## sample estimates:
## mean in group No mean in group Sí 
##         4.210191         5.389767

enbiare %>%
  ggplot( aes(x= desempleo_ord, y=d_total)) +
    geom_boxplot() +
  stat_summary(fun.y="mean", color= "red")+
  theme_classic()

## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: Removed 2 rows containing missing values (geom_segment).

Vamos ahora con préstamos. En efecto, los que piden préstamos también suelen estar desempleados.

t.test(prest_cont ~ desempleo_ord, data = enbiare)

## 
##  Welch Two Sample t-test
## 
## data:  prest_cont by desempleo_ord
## t = -30.131, df = 5141.6, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group No and group Sí is not equal to 0
## 95 percent confidence interval:
##  -0.9784810 -0.8589332
## sample estimates:
## mean in group No mean in group Sí 
##        0.8608278        1.7795349

enbiare %>%
  ggplot( aes(x= desempleo_ord, y=prest_cont)) +
    geom_boxplot() +
  stat_summary(fun.y="mean", color= "red")+
  theme_classic()

## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: Removed 2 rows containing missing values (geom_segment).

Aquí la pregunta difícil es: ¿Desempleo es un confounder o un mecanismo? Tenemos que hipotetizar el orden causal de las variables.

Si lo definimos como confounder: Estoy desempleado, pido préstamos, y me deprimo. Además, pueso estar desempleado y deprimirme incluso sin pedir préstamos.

Si lo definimos como mecanismo: Pido préstamos, lo que me lleva al desempleo, y por ello me deprimo.

Si definimos como causalidad inversa: estoy deprimido, por eso pierdo el trabajo y por tanto pido préstamos.

Esto es lo que debemos tener presente en nuestro diseño:

confounders

Regresión lineal

Empecemos por la relación más simple: prestamos predice depresión

Veamos la correlación: 0.217. Débil pero significativa

enbiare %>% 
  select(d_total, prest_cont) %>% 
  cor()

##              d_total prest_cont
## d_total    1.0000000  0.2177347
## prest_cont 0.2177347  1.0000000

mod1 <- lm(d_total ~ prest_cont, enbiare)

summary(mod1)

## 
## Call:
## lm(formula = d_total ~ prest_cont, data = enbiare)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.503 -3.381 -1.005  2.244 17.244 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.75632    0.02939  127.81   <2e-16 ***
## prest_cont   0.62437    0.01585   39.38   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.391 on 31164 degrees of freedom
## Multiple R-squared:  0.04741,    Adjusted R-squared:  0.04738 
## F-statistic:  1551 on 1 and 31164 DF,  p-value: < 2.2e-16

Cómo se interpreta esto?

Intercept: Una persona que no tiene préstamos, en promedio, tiene 3.7 en el puntaje total de depresión
Pendiente: Cada rubro adicional de préstamo que pide una persona, en promedio, su puntaje de depresión sube 0.62 puntos.

Breve recordatorio para el intercept de 3.75. No es exacto, pero de aquí viene:

enbiare %>% 
  group_by(prest_cont) %>% 
  summarise(medias= mean(d_total))

## # A tibble: 7 × 2
##   prest_cont medias
##        <dbl>  <dbl>
## 1          0   3.65
## 2          1   4.91
## 3          2   5.31
## 4          3   5.72
## 5          4   6.37
## 6          5   6.38
## 7          6   6.73

Veamos qué geom nos ayuda.. Vean la línea y piensen en ambos datos.

enbiare %>% 
  ggplot(aes(x=prest_cont, y= d_total)) +
  geom_point() +
  geom_smooth(method='lm')

## `geom_smooth()` using formula 'y ~ x'

No veo mucho, la verdad. aunque la línea sí es reveladora

enbiare %>% 
  ggplot(aes(x=prest_cont, y= d_total)) +
  geom_jitter(alpha= 0.1) +
  geom_smooth(method='lm')

## `geom_smooth()` using formula 'y ~ x'

Voy a tomar una muestra del 3%

enbiare %>% 
  sample_frac(0.03) %>% # 3% de mis datos aleatoriamente
  ggplot(aes(x=prest_cont, y= d_total)) +
  geom_jitter(alpha= 0.4) +
  geom_smooth(method='lm')

## `geom_smooth()` using formula 'y ~ x'

algo conceptual pero se ve dónde se concentran los datos

enbiare %>% 
  ggplot(aes(x=prest_cont, y= d_total)) +
  stat_bin2d(bins=6)+
  geom_smooth(method='lm')

## `geom_smooth()` using formula 'y ~ x'

Controles

Vamos a empezar a “controlar”

Primero un recordatorio de la reg simple. Vamos a fijarnos en cómo van cambiando los coeficientes del Tx

mod1 <- lm(d_total ~ prest_cont, enbiare)

summary(mod1)

## 
## Call:
## lm(formula = d_total ~ prest_cont, data = enbiare)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.503 -3.381 -1.005  2.244 17.244 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.75632    0.02939  127.81   <2e-16 ***
## prest_cont   0.62437    0.01585   39.38   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.391 on 31164 degrees of freedom
## Multiple R-squared:  0.04741,    Adjusted R-squared:  0.04738 
## F-statistic:  1551 on 1 and 31164 DF,  p-value: < 2.2e-16

Dos X

mod2 <- lm(d_total ~ prest_cont + mujer, enbiare)

summary(mod2)

## 
## Call:
## lm(formula = d_total ~ prest_cont + mujer, data = enbiare)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.983 -3.104 -1.104  1.896 17.896 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.10418    0.03915   79.28   <2e-16 ***
## prest_cont   0.60783    0.01571   38.68   <2e-16 ***
## mujer        1.23197    0.04948   24.90   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.348 on 31163 degrees of freedom
## Multiple R-squared:  0.06599,    Adjusted R-squared:  0.06593 
## F-statistic:  1101 on 2 and 31163 DF,  p-value: < 2.2e-16

Como ya cambiaron los coeficientes, vamos a interpretar primero teóricamente y luego prácticamente

Intercepto: Es el valor promedio de Y cuando todas X1 y X2 son cero
X1: es el cambio en Y al pasar de una unidad a otra de X1, cuando X2 es cero
X2: es el cambio en Y al pasar de una unidad a otra de X2, cuando X1 es cero

Ahora en contexto:

Intercepto: Un varón que no pide préstamos tiene, en promedio, 3.1 del puntaje de depresión
Prestamo (X1): Cada rubro adicional de préstamo que pide un varón, en promedio, su puntaje de depresión sube 0.61 puntos.

Sexo (X2): Una mujer que no pide préstamos tiene, en promedio, un puntaje de depresión 1.23 más alto que un varón (3.1 + 1.2 = 4.3). Esta es la idea:

enbiare %>% 
  select(mujer_ord, d_total, prest_cont) %>% 
  group_by(mujer_ord, prest_cont) %>%
  summarize(media= mean(d_total))

## `summarise()` has grouped output by 'mujer_ord'. You can override using the
## `.groups` argument.

## # A tibble: 14 × 3
## # Groups:   mujer_ord [2]
##    mujer_ord prest_cont media
##    <ord>          <dbl> <dbl>
##  1 Hombre             0  3.12
##  2 Hombre             1  4.20
##  3 Hombre             2  4.59
##  4 Hombre             3  4.51
##  5 Hombre             4  5.30
##  6 Hombre             5  5.06
##  7 Hombre             6  6.19
##  8 Mujer              0  4.13
##  9 Mujer              1  5.51
## 10 Mujer              2  5.86
## 11 Mujer              3  6.60
## 12 Mujer              4  7.11
## 13 Mujer              5  7.33
## 14 Mujer              6  7.11

Veamos gráficamente esto. OJO: sólo es una muestra y por eso no son los coeficientes exactos. Favor de ignorar, por ahora, que las pendientes no son iguales.

enbiare %>% 
  sample_frac(0.10) %>% # 10% de mis datos aleatoriamente
  ggplot(aes(x=prest_cont, y= d_total, group= mujer_ord)) +
  geom_jitter(alpha= 0.3) +
  geom_smooth(method='lm')

## `geom_smooth()` using formula 'y ~ x'

Tres X

mod3 <- lm(d_total ~ prest_cont + mujer + educ, enbiare)

summary(mod3)

## 
## Call:
## lm(formula = d_total ~ prest_cont + mujer + educ, data = enbiare)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.435 -3.137 -1.137  2.065 18.458 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.73213    0.05391   69.22   <2e-16 ***
## prest_cont   0.58110    0.01574   36.93   <2e-16 ***
## mujer        1.21644    0.04932   24.66   <2e-16 ***
## educ        -0.59509    0.03529  -16.86   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.329 on 31101 degrees of freedom
##   (61 observations deleted due to missingness)
## Multiple R-squared:  0.07444,    Adjusted R-squared:  0.07436 
## F-statistic: 833.8 on 3 and 31101 DF,  p-value: < 2.2e-16

cómo interpretamos esto?

plot_model(mod3, vline.color = "red")

¿A quién represnta el valor del intercepto de 3.73? Es el valor promedio del puntaje de depresión de un varon, con educación máxima primaria, que no pide préstamos. Veamos aritméticamente.

enbiare %>% 
  select(mujer_ord, educ_ord, d_total, prest_cont) %>% 
  group_by(mujer_ord, educ_ord, prest_cont) %>%
  summarize(media= mean(d_total))

## `summarise()` has grouped output by 'mujer_ord', 'educ_ord'. You can override
## using the `.groups` argument.

## # A tibble: 53 × 4
## # Groups:   mujer_ord, educ_ord [8]
##    mujer_ord educ_ord     prest_cont media
##    <ord>     <ord>             <dbl> <dbl>
##  1 Hombre    Primaria              0  3.70
##  2 Hombre    Primaria              1  4.62
##  3 Hombre    Primaria              2  5.06
##  4 Hombre    Primaria              3  5.16
##  5 Hombre    Primaria              4  6.10
##  6 Hombre    Primaria              5  5.57
##  7 Hombre    Primaria              6  8.21
##  8 Hombre    Preparatoria          0  2.97
##  9 Hombre    Preparatoria          1  4.01
## 10 Hombre    Preparatoria          2  4.41
## # … with 43 more rows

Vean todas las relaciones que estamos ya integrando. Son muchos números para manejar… Voy a quitar los puntitos.

enbiare %>%
  select(prest_cont, d_total, mujer_ord, educ_ord)%>%
  drop_na()%>%
  ggplot(aes(x=prest_cont, y= d_total, group= mujer_ord, color= mujer_ord)) +
  geom_smooth(method='lm')+
  facet_wrap(~ educ_ord)+
  theme_classic()

## `geom_smooth()` using formula 'y ~ x'

Interpretemos. un par de cosas en estas gráficas.

Lo primero es la diferencia entre las dos líneas por sexo en cada panel (Hagamos de cuenta, por favor, como que son líneas paralelas; en un momento más regreso a este punto, lo prometo). Nuestro modelo de regresión nos dice que, en promedio, la diferencia entre varones y mujeres, controlando por educación y préstamos, es de 1.21 puntos del puntaje de depresión. al pasar de varón a mujer, el puntaje sube en 1.21. Esto hace sentido porque en los 3 niveles educativos, la depresión es mayor en mujeres.

Sin embargo, fíjense que en los tres casos, tanto para hombres y mujeres, mayor educación parece vincularse con menores puntajes de pobreza. Veamos con más claridad esta relación.

enbiare %>%
  select(prest_cont, d_total, mujer_ord, educ_ord)%>%
  drop_na()%>%
  ggplot(aes(x=prest_cont, y= d_total)) +
  geom_smooth(method='lm')+
  facet_wrap(~ educ_ord)+
  theme_classic()

## `geom_smooth()` using formula 'y ~ x'

El modelo dice que, en personas que no piden préstamos, ajustando por sexo, por cada nivel adicional de educación, el puntaje de depresión desciende en 0.59 puntos. Comparen los interceptos.

Efectos marginales

Hasta ahora estoy mostrando aritméticamente y gráficamente las asociaciones que estamos modelando. Pero no estamos usando el modelo. Por eso las líneas no son estrictamente paralelas. como dice el libro, todavía no hemos visto los “weighted averages”. Vamos a verlos directamente:

summary(mod3)

## 
## Call:
## lm(formula = d_total ~ prest_cont + mujer + educ, data = enbiare)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.435 -3.137 -1.137  2.065 18.458 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.73213    0.05391   69.22   <2e-16 ***
## prest_cont   0.58110    0.01574   36.93   <2e-16 ***
## mujer        1.21644    0.04932   24.66   <2e-16 ***
## educ        -0.59509    0.03529  -16.86   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.329 on 31101 degrees of freedom
##   (61 observations deleted due to missingness)
## Multiple R-squared:  0.07444,    Adjusted R-squared:  0.07436 
## F-statistic: 833.8 on 3 and 31101 DF,  p-value: < 2.2e-16

Empezamos la línea en 3.7. De ahí, por cada préstamo adicional, la pendiente sube 0.58 puntos de depresión, controlando por educación y sexo.

plot_model(mod3, type = "pred", terms = "prest_cont")

Vamos a verlo más complejo: líneas paralelas!

plot_model(mod3, type = "pred", terms = c("prest_cont", "mujer"))

La idea es la siguiente: imaginen que llega una nueva persona a su consultorio. Sólo saben que tiene 4 préstamos, pero no saben si es mujer o varón. Tampoco saben qué nivel de depresión podría tener. La gráfica con la línea negra les dice que su mejor apuesta es decir que tendrá un puntaje de 6.2 en la escala de depresión.

Peeero… si ustedes saben que es varón (línea roja), su estimación cambia. Ahora apuestan a que para 4 préstamos su puntaje será de 5.5. Y si es mujer, les convendría apostar a que su puntaje es de 6.7.

¿Y si es una mujer con educación superior? en ese caso yo bajaría mi apuesta a 6.

plot_model(mod3, type = "pred", terms = c("prest_cont", "mujer", "educ"))

Dejenme profundizar en qué son estas gráficas que vemos. Recuerden que un modelo de RL es una ecuación que nos permite predecir. Los coeficientes son los parámetros a los que les metemos un valor nuevo.

d_total = 3.7 + 0.58 (prest_cont) + 1.21(mujer) - 0.59 (educ)

mod3$coefficients

## (Intercept)  prest_cont       mujer        educ 
##   3.7321349   0.5810950   1.2164391  -0.5950904

3.7 + 0.58*4 + 1.21*1 - 0.59*2

## [1] 6.05

Sólo que a mano esto no es práctico. Imaginen que tenemos un bd más grande y nueva

nueva_bd <- tibble::tribble(
              ~mujer, ~educ, ~prest_cont,
                 0L,    0L,          4L,
                 0L,    1L,          4L,
                 0L,    2L,          4L,
                 1L,    0L,          4L,
                 1L,    1L,          4L,
                 1L,    2L,          4L
              )

Ahora vamos a predecir la depresion de estos valores. Cuál es el caso que usamos en nuestra ecuación?

predict(mod3, nueva_bd)

##        1        2        3        4        5        6 
## 6.056515 5.461425 4.866334 7.272954 6.677864 6.082773

Bastante cerca! Y nos dio las predicciones de todas nuestras nuevas combinaciones.

Los efectos marginales lo que hacen es predecir los valores de prestamos a partir de esta ecuacion.

Cuatro X

¿Cambia mucho si ajustamos por empleo?

mod4 <- lm(d_total ~ prest_cont + mujer + educ + desempleo, enbiare)

summary(mod4)

## 
## Call:
## lm(formula = d_total ~ prest_cont + mujer + educ + desempleo, 
##     data = enbiare)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.353 -3.066 -1.066  2.062 18.542 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.67323    0.05417  67.813   <2e-16 ***
## prest_cont   0.54881    0.01605  34.184   <2e-16 ***
## mujer        1.22548    0.04925  24.881   <2e-16 ***
## educ        -0.60744    0.03526 -17.229   <2e-16 ***
## desempleo    0.71049    0.07256   9.791   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.322 on 31100 degrees of freedom
##   (61 observations deleted due to missingness)
## Multiple R-squared:  0.07729,    Adjusted R-squared:  0.07717 
## F-statistic: 651.3 on 4 and 31100 DF,  p-value: < 2.2e-16

comparemos!

 plot_models(mod3, mod4, show.values = TRUE)

O en tabla: cómo interpretamos

tab_model(mod1, mod3, mod4)

	d total			d total			d total
Predictors	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p
(Intercept)	3.76	3.70 – 3.81	<0.001	3.73	3.63 – 3.84	<0.001	3.67	3.57 – 3.78	<0.001
prest cont	0.62	0.59 – 0.66	<0.001	0.58	0.55 – 0.61	<0.001	0.55	0.52 – 0.58	<0.001
mujer				1.22	1.12 – 1.31	<0.001	1.23	1.13 – 1.32	<0.001
educ				-0.60	-0.66 – -0.53	<0.001	-0.61	-0.68 – -0.54	<0.001
desempleo							0.71	0.57 – 0.85	<0.001
Observations	31166			31105			31105
R² / R² adjusted	0.047 / 0.047			0.074 / 0.074			0.077 / 0.077

Estandarización de coeficientes

Durante la clase anterior preguntaron si los coeficientes son comprables o no. Cuando comparamos en dos modelos distintos coeficientes de las mismas variables, entonces estamos comparando la misma escala de la variable; lo que cambia son los ajustes en uno y otro modelo. Eso es lo que se hace arriba al comparar azul y rojo. Sin embargo, estas comparaciones son más difíciles entre variables de un mismo modelo. Nosotros podríamos querer saber si importa más educación o préstamos para explicar depresión. el problema en este segundo caso es que educación tiene una escala diferente a préstamos; la primera tiene 3 valores y la segunda 7. Recuerden que interpretamos que una unidad adicional en X, implica un cambio de tanto en Y. Si comparamos X1 y X2, en este caso, prestamos y educación, la “unidad adicional” es de diferentes tamaños. En otras palabras, estamos comparando los pasos de un gigante con los de un enano.

 table(enbiare$prest_cont)

## 
##     0     1     2     3     4     5     6 
## 19609  3385  2563  2156  1927  1064   462

 table(enbiare$educ)

## 
##     0     1     2 
##  7658 15898  7549

Una salida común a este problema es la estandarización. Este es un procedimiento común para unificar la escala de variables de diferentes tamaños. a veces también se llaman puntaciones Z. La fórmula es sencilla: (x- media(x)) / SD (x).

Vamos a verlo más despacio. Estandaricemos edad

class(enbiare$EDAD)

## [1] "character"

enbiare <- enbiare %>% 
  mutate(edad= as.numeric(EDAD))

summary(enbiare$edad)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.00   30.00   41.00   43.01   55.00   98.00

sd(enbiare$edad)

## [1] 16.46359

Con esta información, lo que haremos es que a cada valor de cada persona, le vamos a restar la media (43.01) y el resultado lo dividiremos entre la desviación estandard (16.4). Y ahora vean el resumen. La característica central de la estandarización es que toda variable la pone con una media de 0 y una sd de 1.

enbiare <- enbiare %>% 
  mutate(edad_z= ((edad - 43.01) / 16.4))

summary(enbiare$edad_z)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -1.525000 -0.793293 -0.122561  0.000261  0.731098  3.353049

sd(enbiare$edad_z)

## [1] 1.003878

No me creen? Hagamos lo mismo con depresión

summary(enbiare$d_total)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   1.000   3.000   4.373   7.000  21.000

enbiare <- enbiare %>% 
  mutate(depresion_z= ((d_total - mean(d_total)) / sd(d_total)))

summary(enbiare$depresion_z)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.9720 -0.7497 -0.3052  0.0000  0.5839  3.6958

sd(enbiare$depresion_z)

## [1] 1

Cuando estamos en el contexto de RL hay varias técnicas para hacerlo. Yo sólo les voy a enseñar una, la más robusta. Si ustedes quieren estandarizar una regresión, revisen bien las funcionalidades completas del paquete que usaremos: effectsize de easystats.

La idea es que primero estimamos nuestro modelo con variables no-estandarizadas. Siempre pongan los coeficientes no-estandarizados y después los estandarizados. Después, lo que hace effectsize es reestimar el modelo, sólo que primero convierte todas estandariza todas las variables independientes y luego las mete como insumo. Notn que no estandariza la Y.

modelo original

summary(mod4)

## 
## Call:
## lm(formula = d_total ~ prest_cont + mujer + educ + desempleo, 
##     data = enbiare)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.353 -3.066 -1.066  2.062 18.542 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.67323    0.05417  67.813   <2e-16 ***
## prest_cont   0.54881    0.01605  34.184   <2e-16 ***
## mujer        1.22548    0.04925  24.881   <2e-16 ***
## educ        -0.60744    0.03526 -17.229   <2e-16 ***
## desempleo    0.71049    0.07256   9.791   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.322 on 31100 degrees of freedom
##   (61 observations deleted due to missingness)
## Multiple R-squared:  0.07729,    Adjusted R-squared:  0.07717 
## F-statistic: 651.3 on 4 and 31100 DF,  p-value: < 2.2e-16

modelo con coeficientes estandarizados:

standardize_parameters(mod4, method = "refit")

## # Standardization method: refit
## 
## Parameter   | Coefficient (std.) |         95% CI
## -------------------------------------------------
## (Intercept) |          -4.98e-15 | [-0.01,  0.01]
## prest_cont  |               0.19 | [ 0.18,  0.20]
## mujer       |               0.14 | [ 0.13,  0.15]
## educ        |              -0.09 | [-0.11, -0.08]
## desempleo   |               0.05 | [ 0.04,  0.07]

Noten que cambiaron mucho los coeficientes. También cambió un poco la interpretación. antes decíamos que por una unidad adicional de X1 (un préstamo adicional), el puntaje de Y (depresión) aumentaba en 0.54. Ahora la unidad adicional ya no es un préstamo sino una desviación estandard. Ahora decimos:

por cada desviación estandard adicional en préstamos, depresión aumenta en 0.19.
De igual manera, por cada desviación estandard adicional en educación, depresión disminuye 0.09.

Ahora sí tenemos una mejor forma de comparar entre peras y manzanas. La influencia de préstamos es poco más del doble que en educación. Y en la dirección opuesta.

Es importante señalar que no todo mundo está de acuerdo en estandarizar como solución a los problemas de comparación, pero es mejor que hacerlo de forma cruda. Además, noten algo: ganamos en comparabilidad, pero perdemos en interpretabilidad. Es mucho más dificil interpretar en SD. Las unidades originales son mil veces más intuitivas: “un préstamos adicional…”

Un problema adicional es que la estandarización de variables dicotómicas no es tan sencilla. De hecho, no deberíamos de comparar variables dicotómicas con continuas, aun cuando estén estandarizadas. Si de todas maneras lo quieren hacer, o al menos comparar entre dicotómicas, esta es una mejor solución; se divide entre 2 SD.

standardize_parameters(mod4, method = "refit", two_sd = TRUE)

## # Standardization method: refit
## 
## Parameter   | Coefficient (std.) |         95% CI
## -------------------------------------------------
## (Intercept) |          -4.98e-15 | [-0.01,  0.01]
## prest_cont  |               0.38 | [ 0.36,  0.40]
## mujer       |               0.27 | [ 0.25,  0.29]
## educ        |              -0.19 | [-0.21, -0.17]
## desempleo   |               0.11 | [ 0.09,  0.13]
## 
## - Scaled by two SD(s) from the mean.

Noten que prestamos sigue siendo el doble de educación. Las relaciones no cambiaron, sólo los términos de la comparación. Aquí vemos que el sexo importa más que el desempleo.

¿Listos para publicar? Tenemos que revisar nuestros supuestos…

Supuestos

Al principio de este ejercicio les pedí que me permitieran una violación a la RL. lo que no dije con suficiente claridad es cuáles son las reglas de la RL. en otras palabras: ¿Cuándo se vale usar una regresión lineal y cuándo no? A estas reglas las llamamos supuestos

Normalidad de la Y.
- Esto lo verificamos con gráficas (historgramas) y pruebas de normalidad (Shapiro-Wilk, kolmogorov)
- si no es normal: 1) transformo la variable para hacerla normal; o 2) cambio de tipo de regresión
Linearidad de las asociaciones entre Y y X
- Veo relaciones con gráficas; scatterplots
- Si no es lineal, meto términos cuadráticos
Independencia de Y y X
- Asegurar que el muestreo es independiente. Y verificar que X no es parte de Y y viceversa (ej, predecir IMC con peso)
- Recodificar o eliminar Xs
Multicolinearidad: X1 y X2 miden cosas distintas
- Revisar la matriz de correlaciones entre Xs. Cuidar que no sean altas en algun par de predictores (ej., X1= depresión/ X2= ansiedad)
- Quitar Xs que sean muy parecidas entre sí
Homoscedasticidad: la línea de predicción es igualmente a lo largo de todos los valores. Varianza continua y homogénea.
- Comparar gráficamente residuales vs valores estimados y verificar si hay varianza constante
- si la hay, ver si faltan variables de control o si hay zonas donde la asociación es distinta.

Especificación: se incluyeron todas las variables relevantes

Difícil de diagnósticar porque es parte del diagrama causal
Añadir variables de control faltantes; ojalá sean observadas.

Afortunadamente R, de la mano con easy stats, esta aquí para facilitarnos un poco la vida. Antes de mostrar quiero que que quede claro qué tipo de objeto es un modelo.

names(mod4) # cosas ahí adentro

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "na.action"     "xlevels"       "call"          "terms"        
## [13] "model"

mod4$coefficients

## (Intercept)  prest_cont       mujer        educ   desempleo 
##   3.6732252   0.5488090   1.2254756  -0.6074385   0.7104856

Ya vimos que un modelo es una fórmula. Así que con ella podemos estimar muchas más cosas con la ayuda de broom. Favor de darle a la cuadricula con la flecha, arriba a la derecha.

augment(mod4)

## # A tibble: 31,105 × 12
##    .rownames d_total prest_cont mujer  educ desempleo .fitted .resid      .hat
##    <chr>       <dbl>      <dbl> <dbl> <dbl>     <dbl>   <dbl>  <dbl>     <dbl>
##  1 1               0          0     0     2         0    2.46 -2.46  0.000145 
##  2 2               1          1     1     0         0    5.45 -4.45  0.000127 
##  3 3               0          0     1     2         0    3.68 -3.68  0.000139 
##  4 4               1          0     1     0         0    4.90 -3.90  0.000145 
##  5 5               1          0     1     1         0    4.29 -3.29  0.0000759
##  6 6               1          0     1     0         0    4.90 -3.90  0.000145 
##  7 7               8          0     0     2         0    2.46  5.54  0.000145 
##  8 8               2          0     0     2         0    2.46 -0.458 0.000145 
##  9 9               0          0     1     2         0    3.68 -3.68  0.000139 
## 10 10              0          0     1     2         0    3.68 -3.68  0.000139 
## # … with 31,095 more rows, and 3 more variables: .sigma <dbl>, .cooksd <dbl>,
## #   .std.resid <dbl>

Nos da estas cosas para cada observación!

broom augment

Esa es la información que necesitamos para probar nuestros supuestos

Y ahora sí nos apoyamos en easy stats, con su paquete performance

Primero la magia y luego las pruebas más formales

Empecemos por una prueba gobal (tarda un poco. Si tu compu es lenta toma precauciones). Y luego hazla grande!

check_model(mod4)

Normalidad

Aquí empezamos con el pie izquierdo porque ya sabemos que nuestra distribución de depresión tenía un sesgo claro hacia el cero. Vamos por pruebas más duras.

Nuestros residuales no se distribuyen con normalidad

check_normality(mod4)

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

No se ve bien… y esto ya lo sabíamos por las últimas dos gráficas…

#classify_distribution(mod4)

performance::check_distribution(mod4)

## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##   Distribution Probability
##         normal         47%
##        tweedie         47%
##  beta-binomial          3%
## 
## Predicted Distribution of Response
## 
##                Distribution Probability
##  neg. binomial (zero-infl.)         66%
##               beta-binomial         31%
##        poisson (zero-infl.)          3%

Cuál es esa de negative binomial… Tomada de los amigos de IRDRE de UCLA. Se parece a nuestro histograma de depresión?

negative binomial distribution

La verdad es que, antes de ver los datos, yo me hubiera por una poisson zero inflated:

Poisson

Multicolinearidad

Queremos los VIF (Variance inflation factor) abajo de 5. bien aquí!

performance::check_collinearity(mod4)

## # Check for Multicollinearity
## 
## Low Correlation
## 
##        Term  VIF Increased SE Tolerance
##  prest_cont 1.06         1.03      0.95
##       mujer 1.00         1.00      1.00
##        educ 1.01         1.01      0.99
##   desempleo 1.04         1.02      0.96

Heteroscedasticidad

Tampoco nos gusta mucho…

check_heteroscedasticity(mod4)

## Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).

Precisión

 performance_accuracy(mod4)

## # Accuracy of Model Predictions
## 
## Accuracy: 27.76%
##       SE: 0.74%-points
##   Method: Correlation between observed and predicted

performance_accuracy(mod1)

## # Accuracy of Model Predictions
## 
## Accuracy: 21.78%
##       SE: 1.17%-points
##   Method: Correlation between observed and predicted

compare_performance(mod1, mod3, mod4)

## Warning: When comparing models, please note that probably not all models were fit from
##   same data.

## # Comparison of Model Performance Indices
## 
## Name | Model |       AIC | AIC weights |       BIC | BIC weights |    R2 | R2 (adj.) |  RMSE | Sigma
## ----------------------------------------------------------------------------------------------------
## mod1 |    lm | 1.807e+05 |     < 0.001 | 1.807e+05 |     < 0.001 | 0.047 |     0.047 | 4.391 | 4.391
## mod3 |    lm | 1.794e+05 |     < 0.001 | 1.795e+05 |     < 0.001 | 0.074 |     0.074 | 4.328 | 4.329
## mod4 |    lm | 1.793e+05 |        1.00 | 1.794e+05 |        1.00 | 0.077 |     0.077 | 4.322 | 4.322

compare_performance(mod1, mod3, mod4, rank = TRUE)

## Warning: When comparing models, please note that probably not all models were fit from
##   same data.

## # Comparison of Model Performance Indices
## 
## Name | Model |    R2 | R2 (adj.) |  RMSE | Sigma | AIC weights | BIC weights | Performance-Score
## ------------------------------------------------------------------------------------------------
## mod4 |    lm | 0.077 |     0.077 | 4.322 | 4.322 |        1.00 |        1.00 |           100.00%
## mod3 |    lm | 0.074 |     0.074 | 4.328 | 4.329 |     < 0.001 |     < 0.001 |            60.31%
## mod1 |    lm | 0.047 |     0.047 | 4.391 | 4.391 |     < 0.001 |     < 0.001 |             0.00%

Regresión Lineal Múltiple 1

Pablo Gaitán Rossi

2/15/2022

Limpieza de datos

Depresión

Préstamos

Controles

Sexo

Educación

Desempleo

Regresión lineal

Controles

Dos X

Tres X

Efectos marginales

Cuatro X

Estandarización de coeficientes

Supuestos

Normalidad

Multicolinearidad

Heteroscedasticidad

Precisión