Ivermectin and the odds of hospitalization due to COVID-19: evidence from a quasi-experimental analysis based on a public intervention in Mexico City
code replication
Model 1
df_mod1 <- read_rds(here::here("", "base_publica.rds")) %>%
mutate(kit = kit_administrativo)
kit_hosp(df_mod1)A tibble: 4 x 4
Recibió kit Hospitalizado Total Porcentaje
1 No No 154560 66.09%
2 No Sí 1908 0.82%
3 Sí No 77067 32.96%
4 Sí Sí 314 0.13%
match_mod1 <- matchit(kit ~ hombre + comor_graves + sint_graves + sint_mod + grupo_edad,
data = df_mod1 %>%
select(hosp, kit, hombre, comor_graves, comor_leves, sint_graves, sint_mod,
sint_leves, grupo_edad, contacto_locatel) %>%
na.omit, method = "cem")
match_sum_mod1 <- summary(match_mod1)
match_data_mod1 <- match.data(match_mod1)
glm_mod1 <- glm(hosp ~ kit + hombre + comor_graves + sint_graves + sint_mod + grupo_edad,
data = match_data_mod1, weight = weights, family = "binomial")
match_mod1A matchit object - method: Coarsened exact matching - number of obs.: 233849 (original), 233338 (matched) - target estimand: ATT - covariates: hombre, comor_graves, sint_graves, sint_mod, grupo_edad
print(xtable(glm_mod1), type = "html")| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -5.5444 | 0.0643 | -86.28 | 0.0000 |
| kit | -1.1511 | 0.0613 | -18.79 | 0.0000 |
| hombre | 0.1956 | 0.0424 | 4.62 | 0.0000 |
| comor_graves | 0.2682 | 0.0390 | 6.87 | 0.0000 |
| sint_graves | 0.2194 | 0.0337 | 6.52 | 0.0000 |
| sint_mod | 0.1805 | 0.0186 | 9.71 | 0.0000 |
| grupo_edad30 a 40 | 0.5316 | 0.0720 | 7.39 | 0.0000 |
| grupo_edad40 a 50 | 0.7502 | 0.0694 | 10.80 | 0.0000 |
| grupo_edad50 a 60 | 0.9556 | 0.0710 | 13.46 | 0.0000 |
| grupo_edad60 a 70 | 1.1801 | 0.0803 | 14.70 | 0.0000 |
| grupo_edadMayores de 70 | 1.4632 | 0.0920 | 15.91 | 0.0000 |
Margins 1
mar_mod1 <- margins(glm_mod1, vcov = vcovHC(glm_mod1, type = "HC0")) %>%
broom::tidy()
mar_mod1 %>%
mutate(low = estimate - 1.96 * std.error, up = estimate + 1.96 * std.error, term = case_when(term ==
"kit" ~ "Medical kit", term == "hombre" ~ "Male", term == "comor_graves" ~
"Severe comorbidities", term == "sint_graves" ~ "Severe symptoms", term ==
"sint_mod" ~ "Moderate symptoms", term == "grupo_edad30 a 40" ~ "31-40 years",
term == "grupo_edad40 a 50" ~ "41-50 years", term == "grupo_edad50 a 60" ~
"51-60 years", term == "grupo_edad60 a 70" ~ "61-70 years", term == "grupo_edadMayores de 70" ~
"70+ years", term == "ocupacion" ~ "Occupancy", term == "contacto_locatel" ~
"Locatel tracking")) %>%
ggplot(aes(x = term)) + geom_point(aes(y = estimate)) + geom_errorbar(aes(ymin = low,
ymax = up)) + geom_hline(yintercept = 0, linetype = "dotdash", color = "red") +
theme_bw() + labs(x = "", y = "", title = "Marginal average effects") + coord_flip()
Model 2
df_mod2 <- read_rds(here::here("", "base_publica.rds")) %>%
mutate(kit = kit_administrativo) %>%
filter(contacto_locatel == 1)
kit_hosp(df_mod2)A tibble: 4 x 4
Recibió kit Hospitalizado Total Porcentaje
1 No No 72615 55.13%
2 No Sí 1029 0.78%
3 Sí No 57871 43.93%
4 Sí Sí 211 0.16%
match_mod2 <- matchit(kit ~ hombre + comor_graves + sint_graves + sint_mod + grupo_edad,
data = df_mod2 %>%
select(hosp, kit, hombre, comor_graves, comor_leves, sint_graves, sint_mod,
sint_leves, grupo_edad, contacto_locatel) %>%
na.omit, method = "cem")
match_sum_mod2 <- summary(match_mod2)
match_data_mod2 <- match.data(match_mod2)
glm_mod2 <- glm(hosp ~ kit + hombre + comor_graves + sint_graves + sint_mod + grupo_edad,
data = match_data_mod2, weight = weights, family = "binomial")
match_mod2A matchit object - method: Coarsened exact matching - number of obs.: 131726 (original), 131357 (matched) - target estimand: ATT - covariates: hombre, comor_graves, sint_graves, sint_mod, grupo_edad
print(xtable(glm_mod2), type = "html")| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -5.3369 | 0.0830 | -64.31 | 0.0000 |
| kit | -1.4511 | 0.0757 | -19.17 | 0.0000 |
| hombre | 0.1988 | 0.0557 | 3.57 | 0.0004 |
| comor_graves | 0.2764 | 0.0529 | 5.22 | 0.0000 |
| sint_graves | 0.2456 | 0.0442 | 5.56 | 0.0000 |
| sint_mod | 0.1937 | 0.0246 | 7.87 | 0.0000 |
| grupo_edad30 a 40 | 0.5379 | 0.0921 | 5.84 | 0.0000 |
| grupo_edad40 a 50 | 0.6838 | 0.0901 | 7.59 | 0.0000 |
| grupo_edad50 a 60 | 0.8919 | 0.0926 | 9.64 | 0.0000 |
| grupo_edad60 a 70 | 1.1222 | 0.1052 | 10.66 | 0.0000 |
| grupo_edadMayores de 70 | 1.3970 | 0.1216 | 11.49 | 0.0000 |
Margins 2
mar_mod2 <- margins(glm_mod2, vcov = vcovHC(glm_mod2, type = "HC0")) %>%
broom::tidy()
mar_mod2 %>%
mutate(low = estimate - 1.96 * std.error, up = estimate + 1.96 * std.error, term = case_when(term ==
"kit" ~ "Medical kit", term == "hombre" ~ "Male", term == "comor_graves" ~
"Severe comorbidities", term == "sint_graves" ~ "Severe symptoms", term ==
"sint_mod" ~ "Moderate symptoms", term == "grupo_edad30 a 40" ~ "31-40 years",
term == "grupo_edad40 a 50" ~ "41-50 years", term == "grupo_edad50 a 60" ~
"51-60 years", term == "grupo_edad60 a 70" ~ "61-70 years", term == "grupo_edadMayores de 70" ~
"70+ years", term == "ocupacion" ~ "Occupancy", term == "contacto_locatel" ~
"Locatel tracking")) %>%
ggplot(aes(x = term)) + geom_point(aes(y = estimate)) + geom_errorbar(aes(ymin = low,
ymax = up)) + geom_hline(yintercept = 0, linetype = "dotdash", color = "red") +
theme_bw() + labs(x = "", y = "", title = "Marginal average effects") + coord_flip()
Model 3
df_mod3 <- read_rds(here::here("", "base_publica.rds")) %>%
mutate(kit = kit_administrativo) %>%
filter(contacto_locatel != 1)
kit_hosp(df_mod3)A tibble: 4 x 4
Recibió kit Hospitalizado Total Porcentaje
1 No No 81945 80.24%
2 No Sí 879 0.86%
3 Sí No 19196 18.8%
4 Sí Sí 103 0.1%
match_mod3 <- matchit(kit ~ hombre + comor_graves + sint_graves + sint_mod + grupo_edad,
data = df_mod3 %>%
select(hosp, kit, hombre, comor_graves, comor_leves, sint_graves, sint_mod,
sint_leves, grupo_edad, contacto_locatel) %>%
na.omit, method = "cem")
match_sum_mod3 <- summary(match_mod3)
match_data_mod3 <- match.data(match_mod3)
glm_mod3 <- glm(hosp ~ kit + hombre + comor_graves + sint_graves + sint_mod + grupo_edad,
data = match_data_mod3, weight = weights, family = "binomial")
match_mod3A matchit object - method: Coarsened exact matching - number of obs.: 102123 (original), 101181 (matched) - target estimand: ATT - covariates: hombre, comor_graves, sint_graves, sint_mod, grupo_edad
print(xtable(glm_mod3), type = "html")| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -5.7259 | 0.1016 | -56.35 | 0.0000 |
| kit | -0.6926 | 0.1051 | -6.59 | 0.0000 |
| hombre | 0.1977 | 0.0651 | 3.04 | 0.0024 |
| comor_graves | 0.2038 | 0.0614 | 3.32 | 0.0009 |
| sint_graves | 0.2134 | 0.0516 | 4.13 | 0.0000 |
| sint_mod | 0.1423 | 0.0282 | 5.04 | 0.0000 |
| grupo_edad30 a 40 | 0.4989 | 0.1149 | 4.34 | 0.0000 |
| grupo_edad40 a 50 | 0.8591 | 0.1083 | 7.93 | 0.0000 |
| grupo_edad50 a 60 | 1.0596 | 0.1096 | 9.67 | 0.0000 |
| grupo_edad60 a 70 | 1.2505 | 0.1242 | 10.07 | 0.0000 |
| grupo_edadMayores de 70 | 1.5019 | 0.1428 | 10.52 | 0.0000 |
Margins 3
mar_mod3 <- margins(glm_mod3, vcov = vcovHC(glm_mod3, type = "HC0")) %>%
broom::tidy()
mar_mod3 %>%
mutate(low = estimate - 1.96 * std.error, up = estimate + 1.96 * std.error, term = case_when(term ==
"kit" ~ "Medical kit", term == "hombre" ~ "Male", term == "comor_graves" ~
"Severe comorbidities", term == "sint_graves" ~ "Severe symptoms", term ==
"sint_mod" ~ "Moderate symptoms", term == "grupo_edad30 a 40" ~ "31-40 years",
term == "grupo_edad40 a 50" ~ "41-50 years", term == "grupo_edad50 a 60" ~
"51-60 years", term == "grupo_edad60 a 70" ~ "61-70 years", term == "grupo_edadMayores de 70" ~
"70+ years", term == "ocupacion" ~ "Occupancy", term == "contacto_locatel" ~
"Locatel tracking")) %>%
ggplot(aes(x = term)) + geom_point(aes(y = estimate)) + geom_errorbar(aes(ymin = low,
ymax = up)) + geom_hline(yintercept = 0, linetype = "dotdash", color = "red") +
theme_bw() + labs(x = "", y = "", title = "Marginal average effects") + coord_flip()
Model 4
df_mod4 <- read_rds(here::here("", "base_publica.rds")) %>%
mutate(kit = kit_administrativo) %>%
filter(kit == 1)
kit_hosp(df_mod4)A tibble: 2 x 4
Recibió kit Hospitalizado Total Porcentaje
1 Sí No 77067 99.59%
2 Sí Sí 314 0.41%
match_mod4 <- matchit(contacto_locatel ~ hombre + comor_graves + sint_graves + sint_mod +
grupo_edad, data = df_mod4 %>%
select(contacto_locatel, hosp, hombre, comor_graves, comor_leves, sint_graves,
sint_mod, sint_leves, grupo_edad, contacto_locatel) %>%
na.omit, method = "cem")
match_sum_mod4 <- summary(match_mod4)
match_data_mod4 <- match.data(match_mod4)
glm_mod4 <- glm(hosp ~ contacto_locatel + hombre + comor_graves + sint_graves + sint_mod +
grupo_edad, data = match_data_mod4, weight = weights, family = "binomial")
match_mod4A matchit object - method: Coarsened exact matching - number of obs.: 77381 (original), 76889 (matched) - target estimand: ATT - covariates: hombre, comor_graves, sint_graves, sint_mod, grupo_edad
print(xtable(glm_mod4), type = "html")| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -6.9068 | 0.2160 | -31.98 | 0.0000 |
| contacto_locatel | -0.3499 | 0.1233 | -2.84 | 0.0045 |
| hombre | 0.3999 | 0.1158 | 3.45 | 0.0006 |
| comor_graves | 0.4963 | 0.1016 | 4.89 | 0.0000 |
| sint_graves | 0.1543 | 0.0978 | 1.58 | 0.1145 |
| sint_mod | 0.1819 | 0.0529 | 3.44 | 0.0006 |
| grupo_edad30 a 40 | 0.7053 | 0.2236 | 3.15 | 0.0016 |
| grupo_edad40 a 50 | 1.0984 | 0.2090 | 5.25 | 0.0000 |
| grupo_edad50 a 60 | 1.3178 | 0.2123 | 6.21 | 0.0000 |
| grupo_edad60 a 70 | 1.6491 | 0.2285 | 7.22 | 0.0000 |
| grupo_edadMayores de 70 | 2.0710 | 0.2461 | 8.42 | 0.0000 |
Margins 4
mar_mod4 <- margins(glm_mod4, vcov = vcovHC(glm_mod4, type = "HC0")) %>%
broom::tidy()
mar_mod4 %>%
mutate(low = estimate - 1.96 * std.error, up = estimate + 1.96 * std.error, term = case_when(term ==
"kit" ~ "Medical kit", term == "hombre" ~ "Male", term == "comor_graves" ~
"Severe comorbidities", term == "sint_graves" ~ "Severe symptoms", term ==
"sint_mod" ~ "Moderate symptoms", term == "grupo_edad30 a 40" ~ "31-40 years",
term == "grupo_edad40 a 50" ~ "41-50 years", term == "grupo_edad50 a 60" ~
"51-60 years", term == "grupo_edad60 a 70" ~ "61-70 years", term == "grupo_edadMayores de 70" ~
"70+ years", term == "ocupacion" ~ "Occupancy", term == "contacto_locatel" ~
"Locatel tracking")) %>%
ggplot(aes(x = term)) + geom_point(aes(y = estimate)) + geom_errorbar(aes(ymin = low,
ymax = up)) + geom_hline(yintercept = 0, linetype = "dotdash", color = "red") +
theme_bw() + labs(x = "", y = "", title = "Marginal average effects") + coord_flip()
Model 5
df_mod5 <- read_rds(here::here("", "base_publica.rds")) %>%
mutate(kit = kit_administrativo) %>%
filter(kit != 1)
kit_hosp(df_mod5)A tibble: 2 x 4
Recibió kit Hospitalizado Total Porcentaje
1 No No 154560 98.78%
2 No Sí 1908 1.22%
match_mod5 <- matchit(contacto_locatel ~ hombre + comor_graves + sint_graves + sint_mod +
grupo_edad, data = df_mod5 %>%
select(contacto_locatel, hosp, hombre, comor_graves, comor_leves, sint_graves,
sint_mod, sint_leves, grupo_edad, contacto_locatel) %>%
na.omit, method = "cem")
match_sum_mod5 <- summary(match_mod5)
match_data_mod5 <- match.data(match_mod5)
glm_mod5 <- glm(hosp ~ contacto_locatel + hombre + comor_graves + sint_graves + sint_mod +
grupo_edad, data = match_data_mod5, weight = weights, family = "binomial")
match_mod5A matchit object - method: Coarsened exact matching - number of obs.: 156468 (original), 156065 (matched) - target estimand: ATT - covariates: hombre, comor_graves, sint_graves, sint_mod, grupo_edad
print(xtable(glm_mod5), type = "html")| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -5.7874 | 0.0741 | -78.14 | 0.0000 |
| contacto_locatel | 0.3636 | 0.0476 | 7.64 | 0.0000 |
| hombre | 0.1653 | 0.0473 | 3.49 | 0.0005 |
| comor_graves | 0.2425 | 0.0437 | 5.55 | 0.0000 |
| sint_graves | 0.2378 | 0.0381 | 6.24 | 0.0000 |
| sint_mod | 0.2263 | 0.0190 | 11.90 | 0.0000 |
| grupo_edad30 a 40 | 0.5016 | 0.0785 | 6.39 | 0.0000 |
| grupo_edad40 a 50 | 0.6800 | 0.0768 | 8.86 | 0.0000 |
| grupo_edad50 a 60 | 0.8967 | 0.0785 | 11.42 | 0.0000 |
| grupo_edad60 a 70 | 1.0443 | 0.0919 | 11.37 | 0.0000 |
| grupo_edadMayores de 70 | 1.2974 | 0.1049 | 12.37 | 0.0000 |
Margins 5
mar_mod5 <- margins(glm_mod5, vcov = vcovHC(glm_mod5, type = "HC0")) %>%
broom::tidy()
mar_mod5 %>%
mutate(low = estimate - 1.96 * std.error, up = estimate + 1.96 * std.error, term = case_when(term ==
"kit" ~ "Medical kit", term == "hombre" ~ "Male", term == "comor_graves" ~
"Severe comorbidities", term == "sint_graves" ~ "Severe symptoms", term ==
"sint_mod" ~ "Moderate symptoms", term == "grupo_edad30 a 40" ~ "31-40 years",
term == "grupo_edad40 a 50" ~ "41-50 years", term == "grupo_edad50 a 60" ~
"51-60 years", term == "grupo_edad60 a 70" ~ "61-70 years", term == "grupo_edadMayores de 70" ~
"70+ years", term == "ocupacion" ~ "Occupancy", term == "contacto_locatel" ~
"Locatel tracking")) %>%
ggplot(aes(x = term)) + geom_point(aes(y = estimate)) + geom_errorbar(aes(ymin = low,
ymax = up)) + geom_hline(yintercept = 0, linetype = "dotdash", color = "red") +
theme_bw() + labs(x = "", y = "", title = "Marginal average effects") + coord_flip()
Model 6
df_mod6 <- read_rds(here::here("", "base_publica.rds")) %>%
filter(kit_administrativo == 1 | (fecha < "2020-12-28" & sint_graves + sint_leves +
sint_mod > 0)) %>%
filter(fecha >= "2020-12-15") %>%
filter(ocupacion >= 0.8 & ocupacion <= 0.85) %>%
mutate(kit = kit_administrativo)
kit_hosp(df_mod6)A tibble: 4 x 4
Recibió kit Hospitalizado Total Porcentaje
1 No No 11806 37.09%
2 No Sí 254 0.8%
3 Sí No 19639 61.7%
4 Sí Sí 133 0.42%
match_mod6 <- matchit(kit ~ hombre + comor_graves + grupo_edad, data = df_mod6 %>%
select(hosp, kit, hombre, comor_graves, comor_leves, grupo_edad) %>%
na.omit, method = "cem")
match_sum_mod6 <- summary(match_mod6)
match_data_mod6 <- match.data(match_mod6)
glm_mod6 <- glm(hosp ~ kit + hombre + comor_graves + grupo_edad, data = match_data_mod6,
weight = weights, family = "binomial")
match_mod6A matchit object - method: Coarsened exact matching - number of obs.: 31832 (original), 31817 (matched) - target estimand: ATT - covariates: hombre, comor_graves, grupo_edad
print(xtable(glm_mod6), type = "html")| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -4.7115 | 0.1401 | -33.63 | 0.0000 |
| kit | -1.1571 | 0.1083 | -10.69 | 0.0000 |
| hombre | 0.4766 | 0.1044 | 4.56 | 0.0000 |
| comor_graves | 0.2975 | 0.0989 | 3.01 | 0.0026 |
| grupo_edad30 a 40 | 0.3369 | 0.1726 | 1.95 | 0.0509 |
| grupo_edad40 a 50 | 0.7622 | 0.1589 | 4.80 | 0.0000 |
| grupo_edad50 a 60 | 0.8089 | 0.1688 | 4.79 | 0.0000 |
| grupo_edad60 a 70 | 1.0943 | 0.1913 | 5.72 | 0.0000 |
| grupo_edadMayores de 70 | 1.1420 | 0.2383 | 4.79 | 0.0000 |
Margins 6
mar_mod6 <- margins(glm_mod6, vcov = vcovHC(glm_mod6, type = "HC0")) %>%
broom::tidy()
mar_mod6 %>%
mutate(low = estimate - 1.96 * std.error, up = estimate + 1.96 * std.error, term = case_when(term ==
"kit" ~ "Medical kit", term == "hombre" ~ "Male", term == "comor_graves" ~
"Severe comorbidities", term == "sint_graves" ~ "Severe symptoms", term ==
"sint_mod" ~ "Moderate symptoms", term == "grupo_edad30 a 40" ~ "31-40 years",
term == "grupo_edad40 a 50" ~ "41-50 years", term == "grupo_edad50 a 60" ~
"51-60 years", term == "grupo_edad60 a 70" ~ "61-70 years", term == "grupo_edadMayores de 70" ~
"70+ years", term == "ocupacion" ~ "Occupancy", term == "contacto_locatel" ~
"Locatel tracking")) %>%
ggplot(aes(x = term)) + geom_point(aes(y = estimate)) + geom_errorbar(aes(ymin = low,
ymax = up)) + geom_hline(yintercept = 0, linetype = "dotdash", color = "red") +
theme_bw() + labs(x = "", y = "", title = "Marginal average effects") + coord_flip()
Model 7
df_mod7 <- read_rds(here::here("", "base_publica.rds")) %>%
filter(fecha >= "2020-12-15") %>%
filter(!is.na(kit))
kit_hosp(df_mod7)A tibble: 4 x 4
Recibió kit Hospitalizado Total Porcentaje
1 No No 56924 75.22%
2 No Sí 674 0.89%
3 Sí No 18018 23.81%
4 Sí Sí 56 0.07%
match_mod7 <- matchit(kit ~ hombre + comor_graves + grupo_edad, data = df_mod7 %>%
select(hosp, kit, hombre, comor_graves, comor_leves, grupo_edad) %>%
na.omit, method = "cem")
match_sum_mod7 <- summary(match_mod7)
match_data_mod7 <- match.data(match_mod7)
glm_mod7 <- glm(hosp ~ kit + hombre + comor_graves + grupo_edad, data = match_data_mod7,
weight = weights, family = "binomial")
match_mod7A matchit object - method: Coarsened exact matching - number of obs.: 75672 (original), 75655 (matched) - target estimand: ATT - covariates: hombre, comor_graves, grupo_edad
print(xtable(glm_mod7), type = "html")| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | -5.0738 | 0.0971 | -52.26 | 0.0000 |
| kit | -1.3640 | 0.1392 | -9.80 | 0.0000 |
| hombre | 0.1870 | 0.0738 | 2.53 | 0.0113 |
| comor_graves | 0.4482 | 0.0655 | 6.84 | 0.0000 |
| grupo_edad30 a 40 | 0.4554 | 0.1200 | 3.79 | 0.0001 |
| grupo_edad40 a 50 | 0.5867 | 0.1165 | 5.04 | 0.0000 |
| grupo_edad50 a 60 | 0.7112 | 0.1206 | 5.90 | 0.0000 |
| grupo_edad60 a 70 | 0.6614 | 0.1448 | 4.57 | 0.0000 |
| grupo_edadMayores de 70 | 0.7401 | 0.1833 | 4.04 | 0.0001 |
Margins7
mar_mod7 <- margins(glm_mod7, vcov = vcovHC(glm_mod7, type = "HC0")) %>%
broom::tidy()
mar_mod7 %>%
mutate(low = estimate - 1.96 * std.error, up = estimate + 1.96 * std.error, term = case_when(term ==
"kit" ~ "Medical kit", term == "hombre" ~ "Male", term == "comor_graves" ~
"Severe comorbidities", term == "sint_graves" ~ "Severe symptoms", term ==
"sint_mod" ~ "Moderate symptoms", term == "grupo_edad30 a 40" ~ "31-40 years",
term == "grupo_edad40 a 50" ~ "41-50 years", term == "grupo_edad50 a 60" ~
"51-60 years", term == "grupo_edad60 a 70" ~ "61-70 years", term == "grupo_edadMayores de 70" ~
"70+ years", term == "ocupacion" ~ "Occupancy", term == "contacto_locatel" ~
"Locatel tracking")) %>%
ggplot(aes(x = term)) + geom_point(aes(y = estimate)) + geom_errorbar(aes(ymin = low,
ymax = up)) + geom_hline(yintercept = 0, linetype = "dotdash", color = "red") +
theme_bw() + labs(x = "", y = "", title = "Marginal average effects") + coord_flip()