# Purpose: Produce supplementary figures and tables (Figs. S1-S20 and Tables S1-S12) # Project: McMahon, Baylis, Sweeney, and Funk (2025) # Author: Katie McMahon #-------------------Base Setup-------------------------------------------------- #--clear environment rm(list=ls()) #--load libraries library(tidyverse) library(janitor) library(sf) library(fixest) library(modelsummary) dodge <- position_dodge(width = 0.5) # for plotting with modelsummary #--define directories inputdir <- "YOUR_FOLDER/" # this inputdir should now be your `outdir` from the 03 script outdir <- "YOUR_FOLDER/" # this is where you will store all tables and figures #-------------------Read data--------------------------------------------------- #--load clean linked heat-DHS dataset kid <- readRDS(paste0(inputdir, "merged_sa_dhs_clean.Rds")) #--load daily tmax and wbgtmax datasets tmax <- readRDS(paste0(inputdir, "tmax_daily_sa_dhs.Rds")) wbgtmax <- readRDS(paste0(inputdir, "wbgtmax_daily_sa_dhs.Rds")) #-------------------Create births dataset--------------------------------------- #--for births # need to aggregate to state-by-survey wave level admin1_mo <- kid %>% filter(kidid %in% kid_mod$kidid | kidid %in% kid_died$kidid) %>% mutate(birth = 1) %>% mutate(admin1_sample = paste0(admin1_lbl, sample_lbl), # we want to combine with sample_lbl because some states appear more than others b/c of nonstandard DHS survey frequency admin1_pop = case_when(kidbirthdate <= ymd("2002-06-15") ~ admin1_f_pop_2000, kidbirthdate > ymd("2002-06-15") & kidbirthdate <= ymd("2007-06-15") ~ admin1_f_pop_2005, kidbirthdate > ymd("2007-06-15") & kidbirthdate <= ymd("2012-06-15") ~ admin1_f_pop_2010, kidbirthdate > ymd("2012-06-15") ~ admin1_f_pop_2015)) %>% group_by(admin1_sample, kidbirthdate) %>% summarize( admin1_lbl = first(admin1_lbl), admin1_pop = first(admin1_pop), country_lbl = first(country_lbl), births = sum(birth), # total number of kids born in a given year-month in a given admin 1 in a given DHS wave neodeaths = sum(neodeath, na.rm = T), infdeaths = sum(infdeath, na.rm = T), tri0_t35_tmax = mean(tri0_t35_tmax), # mean no. of 35tmax days during tri 0 among all kids born in admin1-year-mo (so mean across DHS clusters within a state) tri1_t35_tmax = mean(tri1_t35_tmax), tri2_t35_tmax = mean(tri2_t35_tmax), tri3_t35_tmax = mean(tri3_t35_tmax), tri0_t29_wbgtmax = mean(tri0_t29_wbgtmax), tri1_t29_wbgtmax = mean(tri1_t29_wbgtmax), tri2_t29_wbgtmax = mean(tri2_t29_wbgtmax), tri3_t29_wbgtmax = mean(tri3_t29_wbgtmax), ) %>% ungroup() %>% mutate( month = as.factor(month(kidbirthdate)), year = as.factor(year(kidbirthdate)) ) #-------------------Figure S1--------------------------------------------------- ##----Panel A---- #--trim record # we trim to only the study period to best represent the exposure levels of our sample tmax_mindate <- tmax %>% dplyr::select(-latnum, -longnum) %>% # drop to speed up computations filter(date >= min(kid$triminus1_start), # minimum date dhsid %in% kid$dhsid) # drop any excess clusters #--density plot tmax_mindate %>% ggplot(aes(x = tmax)) + geom_density(color = 'black', fill = 'gray', alpha = 0.3) + scale_x_continuous(limits = c(0, 50), breaks = scales::pretty_breaks(n = 10)) + scale_y_continuous(limits = c(0, 0.0875), breaks = scales::pretty_breaks(n = 4)) + geom_vline(xintercept = 35, color = 'orange', lty = 2) + labs(title = NULL, x = NULL, y = "Density") + theme_classic() + theme(axis.text = element_text(size = 15)) + theme(axis.title = element_text(size = 20)) #--save output ggsave( "FigureS1A.jpeg", path = outdir, width = 8, height = 6, units = "in", dpi = 300 ) ##----Panel B---- #--trim record # we trim to only the study period to best represent the exposure levels of our sample wbgtmax_mindate <- wbgtmax %>% dplyr::select(-latnum, -longnum) %>% # drop to speed up computations filter(date >= min(kid$triminus1_start), # minimum date dhsid %in% kid$dhsid) # drop any excess clusters #--density plot wbgtmax_mindate %>% ggplot(aes(x = wbgtmax)) + geom_density(color = 'black', fill = 'gray', alpha = 0.3) + scale_x_continuous(limits = c(0, 35), breaks = scales::pretty_breaks(n = 10)) + scale_y_continuous(limits = c(0, 0.0875), breaks = scales::pretty_breaks(n = 4)) + geom_vline(xintercept = 29, color = 'orange', lty = 2) + labs(title = NULL, x = NULL, y = "Density") + theme_classic() + theme(axis.text = element_text(size = 15)) + theme(axis.title = element_text(size = 20)) #--save output ggsave( "FigureS1B.jpeg", path = outdir, width = 8, height = 6, units = "in", dpi = 300 ) #-------------------Figure S2--------------------------------------------------- #--density plot kid %>% ggplot(aes(x = haz)) + geom_density(color = 'black', fill = 'gray', alpha = 0.3) + scale_x_continuous(breaks = scales::pretty_breaks(n = 12)) + scale_y_continuous(limits = c(0, 0.35)) + geom_vline(xintercept = -2, color = 'orange', lty = 2) + geom_vline(xintercept = -3, color = 'red', lty = 2) + geom_vline(xintercept = median(kid$haz, na.rm = T), color = 'black', lty = 3, alpha = 0.5) + labs(title = NULL, x = "HAZ", y = "Density") + theme_classic() + theme(axis.text = element_text(size = 15)) + theme(axis.title = element_text(size = 20)) #--save output ggsave( "FigureS2.jpeg", path = outdir, width = 8, height = 6, units = "in", dpi = 300 ) #-------------------Figure S3--------------------------------------------------- # Materials for Figure S3 are in separate files. # ExploreChangesInExtremeWBGT_S_Asia.pro contains the main routine. #-------------------Figure S4--------------------------------------------------- #--set bins kid$bin0 <- kid$tri3_b20_wbgtmax kid$bin2 <- kid$tri3_b22_24_wbgtmax kid$bin3 <- kid$tri3_b24_26_wbgtmax kid$bin4 <- kid$tri3_b26_28_wbgtmax kid$bin5 <- kid$tri3_b28_30_wbgtmax kid$bin6 <- kid$tri3_b30_32_wbgtmax kid$bin7 <- kid$tri3_b32_wbgtmax #--run regression mod <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + bin0 + bin2 + bin3 + bin4 + bin5 + bin6 + bin7 | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") #--collect model output a <- as.data.frame(coef(mod, keep = "bin")) %>% rownames_to_column("var") b <- as.data.frame(se(mod, keep = "bin")) %>% rownames_to_column("var") c <- full_join(a, b) colnames(c)[2:3] <- c("est", "se") c[nrow(c) + 1,] = c("bin1", 0, 0) c <- c %>% mutate(est = as.numeric(est), se = as.numeric(se)) %>% mutate( conf.low = est - 1.96*se, conf.high = est + 1.96*se ) rm(a, b) #--coefficient plot ggplot(c, aes(y = var, x = est, xmin = conf.low, xmax = conf.high, )) + geom_point(size = 3, color = "purple") + geom_linerange(aes(xmin=conf.low, xmax=conf.high), color = "purple") + geom_vline(xintercept = 0, color = 'black', lty = 2) + theme(axis.line = element_line()) + labs( x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL ) + scale_x_continuous(limits = c(-0.0095, 0.0045)) + scale_y_discrete(labels=c("<20C", "20-22C", "22-24C", "24-26C", "26-28C", "28-30C", "30-32C", ">32C")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.3))) + theme(legend.text = element_text(size = rel(1))) + theme(axis.title.y = element_text(size = rel(1.25))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() #--save output ggsave( "FigureS4.jpeg", path = outdir, width = 8, height = 6, units = "in", dpi = 300 ) #-------------------Figure S5--------------------------------------------------- #--set bins kid$bin0 <- kid$tri3_b20_tmax kid$bin2 <- kid$tri3_b25_30_tmax kid$bin3 <- kid$tri3_b30_35_tmax kid$bin4 <- kid$tri3_b35_40_tmax kid$bin5 <- kid$tri3_b40_tmax #--run regression mod <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + bin0 + bin2 + bin3 + bin4 + bin5 | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") #--collect model output a <- as.data.frame(coef(mod, keep = "bin")) %>% rownames_to_column("var") b <- as.data.frame(se(mod, keep = "bin")) %>% rownames_to_column("var") c <- full_join(a, b) colnames(c)[2:3] <- c("est", "se") c[nrow(c) + 1,] = c("bin1", 0, 0) c <- c %>% mutate(est = as.numeric(est), se = as.numeric(se)) %>% mutate( conf.low = est - 1.96*se, conf.high = est + 1.96*se ) rm(a, b) #--coefficient plot ggplot(c, aes(y = var, x = est, xmin = conf.low, xmax = conf.high)) + geom_point(size = 3, color = "red") + geom_linerange(aes(xmin=conf.low, xmax=conf.high), color = "red") + geom_vline(xintercept = 0, color = 'black', lty = 2) + theme(axis.line = element_line()) + labs( x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL ) + scale_x_continuous(limits = c(-0.0065, 0.003)) + scale_y_discrete(labels=c("<20C", "20-25C", "25-30C", "30-35C", "35-40C", ">40C")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.4))) + theme(legend.text = element_text(size = rel(1))) + theme(axis.title.y = element_text(size = rel(1.25))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() #--save output ggsave( "FigureS5.jpeg", path = outdir, width = 8, height = 6, units = "in", dpi = 300 ) #-------------------Figure S6--------------------------------------------------- #--Run separate heat models mod_tmax <- feols(stunt ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t35_tmax + tri1_t35_tmax + tri2_t35_tmax + tri3_t35_tmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") mod_wbgtmax <- feols(stunt ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t29_wbgtmax + tri1_t29_wbgtmax + tri2_t29_wbgtmax + tri3_t29_wbgtmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") # save output using modelsummary modplot <- modelplot(list("Tmax>35C" = mod_tmax, "WBGTmax>29C" = mod_wbgtmax), keep = "tri", draw = FALSE) %>% mutate(term = str_sub(term, end = 4)) modplot_tmax <- modplot %>% filter(model == "Tmax>35C") modplot_wbgtmax <- modplot %>% filter(model == "WBGTmax>29C") #--Run combined heat model mod_both <- feols(stunt ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t35_tmax + tri1_t35_tmax + tri2_t35_tmax + tri3_t35_tmax + tri0_t29_wbgtmax + tri1_t29_wbgtmax + tri2_t29_wbgtmax + tri3_t29_wbgtmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") # save output a <- as.data.frame(coef(mod_both, keep = "tri")) %>% rownames_to_column("var") b <- as.data.frame(se(mod_both, keep = "tri")) %>% rownames_to_column("var") c <- full_join(a, b) colnames(c)[2:3] <- c("est", "se") c <- c %>% mutate( conf.low = est - 1.96*se, conf.high = est + 1.96*se, metric = str_sub(var, start = 6, end = -3) ) %>% mutate(var = str_sub(var, end = 4)) %>% mutate(metric = factor(metric, levels=c("t35_tmax", "t29_wbgtmax"))) rm(a, b) #--Plot panel A ggplot(modplot_tmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL) + scale_color_manual(labels = c("Tmax>35C"), values = c("red")) + scale_shape_manual(labels = c("Tmax>35C"), values = c("circle")) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_x_continuous(limits = c(-0.00075, 0.0011)) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.5))) + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + theme(legend.position = "bottom") + coord_flip() # save output ggsave( "FigureS6A.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #--Plot panel B ggplot(modplot_wbgtmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL) + scale_color_manual(labels = c("WBGTmax>29C"), values = c("purple")) + scale_shape_manual(labels = c("WBGTmax>29C"), values = c("square")) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_x_continuous(limits = c(-0.00075, 0.0011)) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.5))) + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + theme(legend.position = "bottom") + coord_flip() # save output ggsave( "FigureS6B.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #--Plot panel C ggplot(c, aes(y = var, x = est, xmin = conf.low, xmax = conf.high, color = metric, shape = metric)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(c$var)) == "tri0") + which(levels(as.factor(c$var)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL) + scale_x_continuous(limits = c(-0.00075, 0.0011)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("Tmax>35C", "WBGTmax>29C"), values = c("red", "purple")) + scale_shape_manual(labels = c("Tmax>35C", "WBGTmax>29C"), values = c("circle", "square")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.4))) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() # save output ggsave( "FigureS6C.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #-------------------Figure S7--------------------------------------------------- #--run separate heat models mod_births_tmax <- fepois(births ~ tri3_t35_tmax + tri2_t35_tmax + tri1_t35_tmax + tri0_t35_tmax + log(admin1_pop) | admin1_sample + month + year, data = admin1_mo, cluster = "admin1_sample", fixef.rm = "infinite_coef") mod_births_wbgtmax <- fepois(births ~ tri3_t29_wbgtmax + tri2_t29_wbgtmax + tri1_t29_wbgtmax + tri0_t29_wbgtmax + log(admin1_pop) | admin1_sample + month + year, data = admin1_mo, cluster = "admin1_sample", fixef.rm = "infinite_coef") #--collect model output modplot <- modelplot(list("Tmax>35C" = mod_births_tmax, "WBGTmax>29C" = mod_births_wbgtmax), keep = "tri", draw = FALSE) %>% mutate(term = as.factor(str_sub(term, end = 4))) modplot_tmax <- modplot %>% filter(model == "Tmax>35C") modplot_wbgtmax <- modplot %>% filter(model == "WBGTmax>29C") #--plot panel A ggplot(modplot_tmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL, lty = NULL) + scale_x_continuous(limits = c(-0.0075, 0.0085)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("Tmax>35C"), values = c("red")) + scale_shape_manual(labels = c("Tmax>35C"), values = c("circle")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.4))) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() # save output ggsave( "FigureS7A.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #--plot panel B ggplot(modplot_wbgtmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL, lty = NULL) + scale_x_continuous(limits = c(-0.0075, 0.0085)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("WBGTmax>29C"), values = c("purple")) + scale_shape_manual(labels = c("WBGTmax>29C"), values = c("square")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.4))) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() # save output ggsave( "FigureS7B.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #-------------------Figure S8--------------------------------------------------- #--run separate heat models mod_births_tmax <- fenegbin(births ~ tri3_t35_tmax + tri2_t35_tmax + tri1_t35_tmax + tri0_t35_tmax + log(admin1_pop) | admin1_sample + month + year, data = admin1_mo, cluster = "admin1_sample", fixef.rm = "infinite_coef") mod_births_wbgtmax <- fenegbin(births ~ tri3_t29_wbgtmax + tri2_t29_wbgtmax + tri1_t29_wbgtmax + tri0_t29_wbgtmax + log(admin1_pop) | admin1_sample + month + year, data = admin1_mo, cluster = "admin1_sample", fixef.rm = "infinite_coef") #--collect model output modplot <- modelplot(list("Tmax>35C" = mod_births_tmax, "WBGTmax>29C" = mod_births_wbgtmax), keep = "tri", draw = FALSE) %>% mutate(term = as.factor(str_sub(term, end = 4))) modplot_tmax <- modplot %>% filter(model == "Tmax>35C") modplot_wbgtmax <- modplot %>% filter(model == "WBGTmax>29C") #--plot panel A ggplot(modplot_tmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL, lty = NULL) + scale_x_continuous(limits = c(-0.0075, 0.0085)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("Tmax>35C"), values = c("red")) + scale_shape_manual(labels = c("Tmax>35C"), values = c("circle")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.4))) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() # save output ggsave( "FigureS8A.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #--plot panel B ggplot(modplot_wbgtmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL, lty = NULL) + scale_x_continuous(limits = c(-0.0075, 0.0085)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("WBGTmax>29C"), values = c("purple")) + scale_shape_manual(labels = c("WBGTmax>29C"), values = c("square")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.4))) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() # save output ggsave( "FigureS8B.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #-------------------Figure S9--------------------------------------------------- #--run separate heat models mod_infdeath_tmax <- feols(infdeath ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t35_tmax + tri1_t35_tmax + tri2_t35_tmax + tri3_t35_tmax | kidagemo_new + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid_died_mod$perweight, data = kid_died_mod, cluster = "dhsid", fixef.rm = "infinite_coef") # new final argument to return to default state before recent fixest update mod_infdeath_wbgtmax <- feols(infdeath ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t29_wbgtmax + tri1_t29_wbgtmax + tri2_t29_wbgtmax + tri3_t29_wbgtmax | kidagemo_new + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid_died_mod$perweight, data = kid_died_mod, cluster = "dhsid", fixef.rm = "infinite_coef") #--collect model output modplot <- modelplot(list("Tmax>35C" = mod_infdeath_tmax, "WBGTmax>29C" = mod_infdeath_wbgtmax), keep = "tri", draw = FALSE) %>% mutate(term = as.factor(str_sub(term, end = 4))) modplot_tmax <- modplot %>% filter(model == "Tmax>35C") modplot_wbgtmax <- modplot %>% filter(model == "WBGTmax>29C") #--plot panel A ggplot(modplot_tmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL, lty = NULL) + scale_x_continuous(limits = c(-0.0003, 0.0003)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("Tmax>35C"), values = c("red")) + scale_shape_manual(labels = c("Tmax>35C"), values = c("circle")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.4))) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() # save output ggsave( "FigureS9A.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #--plot panel B ggplot(modplot_wbgtmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL, lty = NULL) + scale_x_continuous(limits = c(-0.0003, 0.0003)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("WBGTmax>29C"), values = c("purple")) + scale_shape_manual(labels = c("WBGTmax>29C"), values = c("square")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.4))) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() # save output ggsave( "FigureS9B.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #-------------------Figure S10-------------------------------------------------- #--run separate heat models: 80th percentile, 95th percentile, main thresholds mod_q80_tmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_q80_tmax + tri1_q80_tmax + tri2_q80_tmax + tri3_q80_tmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") mod_q80_wbgtmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_q80_wbgtmax + tri1_q80_wbgtmax + tri2_q80_wbgtmax + tri3_q80_wbgtmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") mod_q95_tmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_q95_tmax + tri1_q95_tmax + tri2_q95_tmax + tri3_q95_tmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") mod_q95_wbgtmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_q95_wbgtmax + tri1_q95_wbgtmax + tri2_q95_wbgtmax + tri3_q95_wbgtmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") mod_t35_tmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t35_tmax + tri1_t35_tmax + tri2_t35_tmax + tri3_t35_tmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") mod_t29_wbgtmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t29_wbgtmax + tri1_t29_wbgtmax + tri2_t29_wbgtmax + tri3_t29_wbgtmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") #--collect model output modplot <- modelplot(list("Tmax>80th percentile" = mod_q80_tmax, "WBGTmax>80th percentile" = mod_q80_wbgtmax, "Tmax>95th percentile" = mod_q95_tmax, "WBGTmax>95th percentile" = mod_q95_wbgtmax, "Tmax>35C" = mod_t35_tmax, "WBGTmax>29C" = mod_t29_wbgtmax), keep = "tri", draw = FALSE) %>% mutate(term = str_sub(term, end = 4)) modplot_tmax <- modplot %>% filter(grepl("^Tmax", model)) modplot_wbgtmax <- modplot %>% filter(grepl("WBGTmax", model)) #--plot panel A ggplot(modplot_tmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, lty = model)) + geom_point(position = dodge, size = 3, color = "red", shape = "circle") + geom_linerange(aes(xmin=conf.low, xmax=conf.high), color = "red", position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL) + scale_x_continuous(limits = c(-0.006, 0.005)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_linetype_manual(labels = c("Tmax>80th %ile", "Tmax>95th %ile", "Tmax>35C"), values = c("dotted", "dashed", "solid")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.3))) + theme(legend.title = element_blank()) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() #--save output ggsave( "FigureS10A", path = outdir, width = 8, height = 6, units = "in", dpi = 300 ) #--plot panel B ggplot(modplot_wbgtmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, lty = model)) + geom_point(position = dodge, size = 3, color = "purple", shape = "square") + geom_linerange(aes(xmin=conf.low, xmax=conf.high), color = "purple", position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL) + scale_x_continuous(limits = c(-0.006, 0.005)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_linetype_manual(labels = c("WBGTmax>80th %ile", "WBGTmax>95th %ile", "WBGTmax>29C"), values = c("dotted", "dashed", "solid")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.3))) + theme(legend.title = element_blank()) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() #--save output ggsave( "FigureS10B.jpeg", path = outdir, width = 8, height = 6, units = "in", dpi = 300 ) #-------------------Figure S11-------------------------------------------------- #--run separate heat models mod_hist_mean_tmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t35_tmax + tri1_t35_tmax + tri2_t35_tmax + tri3_t35_tmax + tri0_hist_mean_tmax + tri1_hist_mean_tmax + tri2_hist_mean_tmax + tri3_hist_mean_tmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") mod_hist_mean_wbgtmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t29_wbgtmax + tri1_t29_wbgtmax + tri2_t29_wbgtmax + tri3_t29_wbgtmax + tri0_hist_mean_wbgtmax + tri1_hist_mean_wbgtmax + tri2_hist_mean_wbgtmax + tri3_hist_mean_wbgtmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") #--collect model output modplot <- modelplot(list("Tmax>35C" = mod_hist_mean_tmax, "WBGTmax>29C" = mod_hist_mean_wbgtmax), keep = "tri", draw = FALSE) %>% filter(!grepl("hist", term)) %>% mutate(term = str_sub(term, end = 4)) modplot_tmax <- modplot %>% filter(model == "Tmax>35C") modplot_wbgtmax <- modplot %>% filter(model == "WBGTmax>29C") #--plot panel A ggplot(modplot_tmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL) + scale_x_continuous(limits = c(-0.005, 0.004)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("Tmax>35C"), values = c("red")) + scale_shape_manual(labels = c("Tmax>35C"), values = c("circle")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.4))) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() #--save output ggsave( "FigureS11A.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #--plot panel B ggplot(modplot_wbgtmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL) + scale_x_continuous(limits = c(-0.005, 0.004)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("WBGTmax>29C"), values = c("purple")) + scale_shape_manual(labels = c("WBGTmax>29C"), values = c("square")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.4))) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() #--save output ggsave( "FigureS11B.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #-------------------Figure S12-------------------------------------------------- #--run separate heat models mod_mean_tmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_mean_tmax + tri1_mean_tmax + tri2_mean_tmax + tri3_mean_tmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") mod_mean_wbgtmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_mean_wbgtmax + tri1_mean_wbgtmax + tri2_mean_wbgtmax + tri3_mean_wbgtmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") #--collect model output modplot <- modelplot(list("Mean Tmax" = mod_mean_tmax, "Mean WBGTmax" = mod_mean_wbgtmax), keep = "tri", draw = FALSE) %>% mutate(term = str_sub(term, end = 4)) modplot_tmax <- modplot %>% filter(model == "Mean Tmax") modplot_wbgtmax <- modplot %>% filter(model == "Mean WBGTmax") #--plot panel A ggplot(modplot_tmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL) + scale_x_continuous(limits = c(-0.07, 0.03)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("Mean Tmax"), values = c("red")) + scale_shape_manual(labels = c("Mean Tmax"), values = c("circle")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.4))) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() #--save output ggsave( "FigureS12A.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #--plot panel B ggplot(modplot_wbgtmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL) + scale_x_continuous(limits = c(-0.07, 0.03)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("Mean WBGTmax"), values = c("purple")) + scale_shape_manual(labels = c("Mean WBGTmax"), values = c("square")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.4))) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() #--save output ggsave( "FigureS12B.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #-------------------Figure S13-------------------------------------------------- #--run separate heat models: main cluster fixed effects vs household fixed effects mod_tmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t35_tmax + tri1_t35_tmax + tri2_t35_tmax + tri3_t35_tmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") mod_wbgtmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t29_wbgtmax + tri1_t29_wbgtmax + tri2_t29_wbgtmax + tri3_t29_wbgtmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") mod_tmax_hhid <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t35_tmax + tri1_t35_tmax + tri2_t35_tmax + tri3_t35_tmax | kidagemo + kidbirthmo + momage + surveymo + hhid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "singleton") mod_wbgtmax_hhid <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t29_wbgtmax + tri1_t29_wbgtmax + tri2_t29_wbgtmax + tri3_t29_wbgtmax | kidagemo + kidbirthmo + momage + surveymo + hhid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "singleton") #--collect model output modplot_dhsid <- modelplot(list("Tmax>35C" = mod_tmax, "WBGTmax>29C" = mod_wbgtmax), keep = "tri", draw = FALSE) %>% mutate(term = str_sub(term, end = 4), fe = "cluster") modplot_hhid <- modelplot(list("Tmax>35C" = mod_tmax_hhid, "WBGTmax>29C" = mod_wbgtmax_hhid), keep = "tri", draw = FALSE) %>% mutate(term = str_sub(term, end = 4), fe = "household") modplot <- rbind(modplot_dhsid, modplot_hhid) modplot_tmax <- modplot %>% filter(grepl("^Tmax", model)) modplot_wbgtmax <- modplot %>% filter(grepl("WBGTmax", model)) #--plot panel A ggplot(modplot_tmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model, lty = fe)) + geom_point(position = dodge, size = 3, color = "red", shape = "circle") + geom_linerange(aes(xmin=conf.low, xmax=conf.high), color = "red", position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL, lty = NULL) + scale_x_continuous(limits = c(-0.0055, 0.003)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_linetype_manual(labels = c("Cluster FE", "Household FE"), values = c("solid", "dashed")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.4))) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() #--save output ggsave( "FigureS13A.jpeg", path = outdir, width = 8, height = 6, units = "in", dpi = 300 ) #--plot panel B ggplot(modplot_wbgtmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model, lty = fe)) + geom_point(position = dodge, size = 3, color = "purple", shape = "square") + geom_linerange(aes(xmin=conf.low, xmax=conf.high), color = "purple", position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(modplot$term) == "tri0") + which(levels(modplot$term) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL, lty = NULL) + scale_x_continuous(limits = c(-0.0055, 0.003)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_linetype_manual(labels = c("Cluster FE", "Household FE"), values = c("solid", "dashed")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.4))) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() #--save output ggsave( "FigureS13B.jpeg", path = outdir, width = 8, height = 6, units = "in", dpi = 300 ) #-------------------Figure S14-------------------------------------------------- # drop missings that fall out of main regression kid_mod <- kid %>% drop_na(haz, kidsex, kidtwin, birthorder, birthplace, momedu, parity, marstat, imptoilet, religion, kidagemo, kidbirthmo, momage, surveymo, dhsid, admin1, surveyyr, perweight) ##----Without BD---- # create dataset without observations from Bangladesh kid_nobd <- kid_mod %>% filter(!country_lbl == "Bangladesh") mod_tmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t35_tmax + tri1_t35_tmax + tri2_t35_tmax + tri3_t35_tmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid_nobd$perweight, data = kid_nobd, cluster = "dhsid", fixef.rm = "infinite_coef") mod_wbgtmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t29_wbgtmax + tri1_t29_wbgtmax + tri2_t29_wbgtmax + tri3_t29_wbgtmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid_nobd$perweight, data = kid_nobd, cluster = "dhsid", fixef.rm = "infinite_coef") # save output using modelsummary modplot <- modelplot(list("Tmax>35C" = mod_tmax, "WBGTmax>29C" = mod_wbgtmax), keep = "tri", draw = FALSE) %>% mutate(term = str_sub(term, end = 4)) modplot_tmax <- modplot %>% filter(model == "Tmax>35C") modplot_wbgtmax <- modplot %>% filter(model == "WBGTmax>29C") #--Run combined heat model mod_both <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t35_tmax + tri1_t35_tmax + tri2_t35_tmax + tri3_t35_tmax + tri0_t29_wbgtmax + tri1_t29_wbgtmax + tri2_t29_wbgtmax + tri3_t29_wbgtmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid_nobd$perweight, data = kid_nobd, cluster = "dhsid", fixef.rm = "infinite_coef") # save output a <- as.data.frame(coef(mod_both, keep = "tri")) %>% rownames_to_column("var") b <- as.data.frame(se(mod_both, keep = "tri")) %>% rownames_to_column("var") c <- full_join(a, b) colnames(c)[2:3] <- c("est", "se") c <- c %>% mutate( conf.low = est - 1.96*se, conf.high = est + 1.96*se, metric = str_sub(var, start = 6, end = -3) ) %>% mutate(var = str_sub(var, end = 4)) %>% mutate(metric = factor(metric, levels=c("t35_tmax", "t29_wbgtmax"))) rm(a, b) #--Plot panel A ggplot(modplot_tmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL, lty = NULL) + scale_x_continuous(limits = c(-0.125, 0.055)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("Tmax>35C"), values = c("red")) + scale_shape_manual(labels = c("Tmax>35C"), values = c("circle")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.position = "none") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() # save output ggsave( "FigureS14A.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #--Plot panel B ggplot(modplot_wbgtmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL, lty = NULL) + scale_x_continuous(limits = c(-0.125, 0.055)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("WBGTmax>29C"), values = c("purple")) + scale_shape_manual(labels = c("WBGTmax>29C"), values = c("square")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.position = "none") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() # save output ggsave( "FigureS14B.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #--Plot panel C ggplot(c, aes(y = var, x = est, xmin = conf.low, xmax = conf.high, color = metric, shape = metric)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(c$var)) == "tri0") + which(levels(as.factor(c$var)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL) + scale_x_continuous(limits = c(-0.125, 0.055)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("Tmax>35C", "WBGTmax>29C"), values = c("red", "purple")) + scale_shape_manual(labels = c("Tmax>35C", "WBGTmax>29C"), values = c("circle", "square")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.4))) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() # save output ggsave( "FigureS14C.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) ##----Without IA---- # create dataset without observations from India kid_noia <- kid_mod %>% filter(!country_lbl == "India") mod_tmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t35_tmax + tri1_t35_tmax + tri2_t35_tmax + tri3_t35_tmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid_noia$perweight, data = kid_noia, cluster = "dhsid", fixef.rm = "infinite_coef") mod_wbgtmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t29_wbgtmax + tri1_t29_wbgtmax + tri2_t29_wbgtmax + tri3_t29_wbgtmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid_noia$perweight, data = kid_noia, cluster = "dhsid", fixef.rm = "infinite_coef") # save output using modelsummary modplot <- modelplot(list("Tmax>35C" = mod_tmax, "WBGTmax>29C" = mod_wbgtmax), keep = "tri", draw = FALSE) %>% mutate(term = str_sub(term, end = 4)) modplot_tmax <- modplot %>% filter(model == "Tmax>35C") modplot_wbgtmax <- modplot %>% filter(model == "WBGTmax>29C") #--Run combined heat model mod_both <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t35_tmax + tri1_t35_tmax + tri2_t35_tmax + tri3_t35_tmax + tri0_t29_wbgtmax + tri1_t29_wbgtmax + tri2_t29_wbgtmax + tri3_t29_wbgtmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid_noia$perweight, data = kid_noia, cluster = "dhsid", fixef.rm = "infinite_coef") # save output a <- as.data.frame(coef(mod_both, keep = "tri")) %>% rownames_to_column("var") b <- as.data.frame(se(mod_both, keep = "tri")) %>% rownames_to_column("var") c <- full_join(a, b) colnames(c)[2:3] <- c("est", "se") c <- c %>% mutate( conf.low = est - 1.96*se, conf.high = est + 1.96*se, metric = str_sub(var, start = 6, end = -3) ) %>% mutate(var = str_sub(var, end = 4)) %>% mutate(metric = factor(metric, levels=c("t35_tmax", "t29_wbgtmax"))) rm(a, b) #--Plot panel A ggplot(modplot_tmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL, lty = NULL) + scale_x_continuous(limits = c(-0.125, 0.055)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("Tmax>35C"), values = c("red")) + scale_shape_manual(labels = c("Tmax>35C"), values = c("circle")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.position = "none") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() # save output ggsave( "FigureS14D.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #--Plot panel B ggplot(modplot_wbgtmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL, lty = NULL) + scale_x_continuous(limits = c(-0.125, 0.055)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("WBGTmax>29C"), values = c("purple")) + scale_shape_manual(labels = c("WBGTmax>29C"), values = c("square")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.position = "none") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() # save output ggsave( "FigureS14E.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #--Plot panel C ggplot(c, aes(y = var, x = est, xmin = conf.low, xmax = conf.high, color = metric, shape = metric)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(c$var)) == "tri0") + which(levels(as.factor(c$var)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL) + scale_x_continuous(limits = c(-0.125, 0.055)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("Tmax>35C", "WBGTmax>29C"), values = c("red", "purple")) + scale_shape_manual(labels = c("Tmax>35C", "WBGTmax>29C"), values = c("circle", "square")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.4))) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() # save output ggsave( "FigureS14F.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) ##----Without NP---- # create dataset without observations from Nepal kid_nonp <- kid_mod %>% filter(!country_lbl == "Nepal") mod_tmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t35_tmax + tri1_t35_tmax + tri2_t35_tmax + tri3_t35_tmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid_nonp$perweight, data = kid_nonp, cluster = "dhsid", fixef.rm = "infinite_coef") mod_wbgtmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t29_wbgtmax + tri1_t29_wbgtmax + tri2_t29_wbgtmax + tri3_t29_wbgtmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid_nonp$perweight, data = kid_nonp, cluster = "dhsid", fixef.rm = "infinite_coef") # save output using modelsummary modplot <- modelplot(list("Tmax>35C" = mod_tmax, "WBGTmax>29C" = mod_wbgtmax), keep = "tri", draw = FALSE) %>% mutate(term = str_sub(term, end = 4)) modplot_tmax <- modplot %>% filter(model == "Tmax>35C") modplot_wbgtmax <- modplot %>% filter(model == "WBGTmax>29C") #--Run combined heat model mod_both <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t35_tmax + tri1_t35_tmax + tri2_t35_tmax + tri3_t35_tmax + tri0_t29_wbgtmax + tri1_t29_wbgtmax + tri2_t29_wbgtmax + tri3_t29_wbgtmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid_nonp$perweight, data = kid_nonp, cluster = "dhsid", fixef.rm = "infinite_coef") # save output a <- as.data.frame(coef(mod_both, keep = "tri")) %>% rownames_to_column("var") b <- as.data.frame(se(mod_both, keep = "tri")) %>% rownames_to_column("var") c <- full_join(a, b) colnames(c)[2:3] <- c("est", "se") c <- c %>% mutate( conf.low = est - 1.96*se, conf.high = est + 1.96*se, metric = str_sub(var, start = 6, end = -3) ) %>% mutate(var = str_sub(var, end = 4)) %>% mutate(metric = factor(metric, levels=c("t35_tmax", "t29_wbgtmax"))) rm(a, b) #--Plot panel A ggplot(modplot_tmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL, lty = NULL) + scale_x_continuous(limits = c(-0.125, 0.055)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("Tmax>35C"), values = c("red")) + scale_shape_manual(labels = c("Tmax>35C"), values = c("circle")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.position = "none") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() # save output ggsave( "FigureS14G.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #--Plot panel B ggplot(modplot_wbgtmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL, lty = NULL) + scale_x_continuous(limits = c(-0.125, 0.055)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("WBGTmax>29C"), values = c("purple")) + scale_shape_manual(labels = c("WBGTmax>29C"), values = c("square")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.position = "none") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() # save output ggsave( "FigureS14H.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #--Plot panel C ggplot(c, aes(y = var, x = est, xmin = conf.low, xmax = conf.high, color = metric, shape = metric)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(c$var)) == "tri0") + which(levels(as.factor(c$var)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL) + scale_x_continuous(limits = c(-0.125, 0.055)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("Tmax>35C", "WBGTmax>29C"), values = c("red", "purple")) + scale_shape_manual(labels = c("Tmax>35C", "WBGTmax>29C"), values = c("circle", "square")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.4))) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() # save output ggsave( "FigureS14I.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #-------------------Figure S15-------------------------------------------------- #--run separate heat models mod_tmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t40_tmax + tri1_t40_tmax + tri2_t40_tmax + tri3_t40_tmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") mod_wbgtmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t31_wbgtmax + tri1_t31_wbgtmax + tri2_t31_wbgtmax + tri3_t31_wbgtmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") #--collect model output modplot <- modelplot(list("Tmax>35C" = mod_tmax, "WBGTmax>29C" = mod_wbgtmax), keep = "tri", draw = FALSE) %>% mutate(term = str_sub(term, end = 4)) modplot_tmax <- modplot %>% filter(grepl("^Tmax", model)) modplot_wbgtmax <- modplot %>% filter(grepl("WBGTmax", model)) #--Plot panel A ggplot(modplot_tmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL, lty = NULL) + scale_x_continuous(limits = c(-0.002, 0.0045)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("Tmax>40C"), values = c("brown")) + scale_shape_manual(labels = c("Tmax>40C"), values = c("triangle")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.position = "none") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() ggsave( "FigureS15A.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #--Plot panel B ggplot(modplot_wbgtmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL, lty = NULL) + scale_x_continuous(limits = c(-0.002, 0.0045)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("WBGTmax>31C"), values = c("darkblue")) + scale_shape_manual(labels = c("WBGTmax>31C"), values = c("triangle")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.position = "none") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() ggsave( "FigureS15B.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #-------------------Figure S16-------------------------------------------------- #--run separate heat models mod_precip_tmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t35_tmax + tri1_t35_tmax + tri2_t35_tmax + tri3_t35_tmax + tri0_tot_precip + tri1_tot_precip + tri2_tot_precip + tri3_tot_precip | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") mod_precip_wbgtmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + tri0_t29_wbgtmax + tri1_t29_wbgtmax + tri2_t29_wbgtmax + tri3_t29_wbgtmax + tri0_tot_precip + tri1_tot_precip + tri2_tot_precip + tri3_tot_precip | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") #--collect model output modplot <- modelplot(list("Tmax>35C" = mod_precip_tmax, "WBGTmax>29C" = mod_precip_wbgtmax), keep = "tri", draw = FALSE) %>% filter(!grepl("precip", term)) %>% mutate(term = str_sub(term, end = 4)) modplot_tmax <- modplot %>% filter(model == "Tmax>35C") modplot_wbgtmax <- modplot %>% filter(model == "WBGTmax>29C") #--plot panel A ggplot(modplot_tmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL) + scale_x_continuous(limits = c(-0.005, 0.003)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("Tmax>35C"), values = c("red")) + scale_shape_manual(labels = c("Tmax>35C"), values = c("circle")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.4))) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() #--save output ggsave( "FigureS16A.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #--plot panel B ggplot(modplot_wbgtmax, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL) + scale_x_continuous(limits = c(-0.005, 0.003)) + scale_y_discrete(labels=c("Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3")) + scale_color_manual(labels = c("WBGTmax>29C"), values = c("purple")) + scale_shape_manual(labels = c("WBGTmax>29C"), values = c("square")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.4))) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() #--save output ggsave( "FigureS16B.jpeg", path = outdir, width = 6, height = 6, units = "in", dpi = 300 ) #-------------------Figure S17-------------------------------------------------- #--run WBGTmax model mod <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + triminus1_t29_wbgtmax + tri0_t29_wbgtmax + tri1_t29_wbgtmax + tri2_t29_wbgtmax + tri3_t29_wbgtmax + tri4_t29_wbgtmax | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") #--collect model output modplot <- modelplot(list("WBGTmax>29C" = mod), keep = "tri", draw = FALSE) %>% mutate(term = str_sub(term, end = 4)) %>% mutate(term = case_when(term == "trim" ~ "tri-1", T ~ term)) #--coefficient plot ggplot(modplot, aes(y = term, x = estimate, xmin = conf.low, xmax = conf.high, color = model, shape = model)) + geom_point(position = dodge, size = 3) + geom_linerange(aes(xmin=conf.low, xmax=conf.high), position = dodge) + geom_vline(xintercept = 0, color = 'black', lty = 2) + geom_hline(yintercept = (which(levels(as.factor(modplot$term)) == "tri0") + which(levels(as.factor(modplot$term)) == "tri1")) / 2, color = 'orange', alpha = 0.2, size = 1.5, lty = 1) + theme(axis.line = element_line()) + labs(title = NULL, x = "Coefficient estimates and 95% CI", y = NULL, color = NULL, shape = NULL) + scale_y_discrete(labels=c("Tri. -1", "Tri. 0", "Tri. 1", "Tri. 2", "Tri. 3", "Tri. 4")) + scale_x_continuous(limits = c(-0.005, 0.0021)) + scale_color_manual(labels = c("WBGTmax>29C"), values = c("purple")) + scale_shape_manual(labels = c("WBGTmax>29C"), values = c("square")) + theme_classic() + theme(plot.title = element_text(size = rel(1.5))) + theme(axis.text = element_text(size = rel(1.5))) + theme(legend.text = element_text(size = rel(1.5))) + theme(legend.position = "bottom") + theme(axis.title.y = element_text(size = rel(1.5))) + theme(panel.grid.major = element_line(size = 0.3)) + theme(panel.grid.minor = element_line(size = 0.1)) + coord_flip() #--save output ggsave( "FigureS17.jpeg", path = outdir, width = 8, height = 6, units = "in", dpi = 300 ) #-------------------Figure S18-------------------------------------------------- #--trim record # we trim to only the study period to best represent the exposure levels of our sample tmax_mindate <- tmax %>% dplyr::select(-latnum, -longnum) %>% # drop to speed up computations filter(date >= min(kid$triminus1_start), # minimum date dhsid %in% kid$dhsid) # drop any excess clusters #--panel A # create subset of tmax dataset for just Bangladesh tmax_bd <- tmax_mindate %>% filter(grepl("^BD", dhsid)) # density plot tmax_bd %>% ggplot(aes(x = tmax)) + geom_density(color = 'black', fill = 'gray', alpha = 0.3) + scale_x_continuous(limits = c(0, 50), breaks = scales::pretty_breaks(n = 10)) + scale_y_continuous(limits = c(0, 0.15), breaks = scales::pretty_breaks(n = 4)) + geom_vline(xintercept = 35, color = 'orange', lty = 2) + labs(title = NULL, x = NULL, y = "Density") + theme_classic() + theme(axis.text = element_text(size = 15)) + theme(axis.title = element_text(size = 20)) # save output ggsave( "FigureS18A.jpeg", path = outdir, width = 8, height = 6, units = "in", dpi = 300 ) #--panel B # create subset of tmax dataset for just India tmax_ia <- tmax_mindate %>% filter(grepl("^IA", dhsid)) # density plot tmax_ia %>% ggplot(aes(x = tmax)) + geom_density(color = 'black', fill = 'gray', alpha = 0.3) + scale_x_continuous(limits = c(0, 50), breaks = scales::pretty_breaks(n = 10)) + scale_y_continuous(limits = c(0, 0.15), breaks = scales::pretty_breaks(n = 4)) + geom_vline(xintercept = 35, color = 'orange', lty = 2) + labs(title = NULL, x = NULL, y = "Density") + theme_classic() + theme(axis.text = element_text(size = 15)) + theme(axis.title = element_text(size = 20)) # save output ggsave( "FigureS18B.jpeg", path = outdir, width = 8, height = 6, units = "in", dpi = 300 ) #--panel C # create subset of tmax dataset for just Bangladesh tmax_np <- tmax_mindate %>% filter(grepl("^NP", dhsid)) # density plot tmax_np %>% ggplot(aes(x = tmax)) + geom_density(color = 'black', fill = 'gray', alpha = 0.3) + scale_x_continuous(limits = c(0, 50), breaks = scales::pretty_breaks(n = 10)) + scale_y_continuous(limits = c(0, 0.15), breaks = scales::pretty_breaks(n = 4)) + geom_vline(xintercept = 35, color = 'orange', lty = 2) + labs(title = NULL, x = NULL, y = "Density") + theme_classic() + theme(axis.text = element_text(size = 15)) + theme(axis.title = element_text(size = 20)) # save output ggsave( "FigureS18C.jpeg", path = outdir, width = 8, height = 6, units = "in", dpi = 300 ) #-------------------Figure S19-------------------------------------------------- #--trim record # we trim to only the study period to best represent the exposure levels of our sample wbgtmax_mindate <- wbgtmax %>% dplyr::select(-latnum, -longnum) %>% # drop to speed up computations filter(date >= min(kid$triminus1_start), # minimum date dhsid %in% kid$dhsid) # drop any excess clusters #--panel A # create subset of wbgtmax dataset for just Bangladesh wbgtmax_bd <- wbgtmax_mindate %>% filter(grepl("^BD", dhsid)) # density plot wbgtmax_bd %>% ggplot(aes(x = wbgtmax)) + geom_density(color = 'black', fill = 'gray', alpha = 0.3) + scale_x_continuous(limits = c(0, 35), breaks = scales::pretty_breaks(n = 10)) + scale_y_continuous(limits = c(0, 0.155), breaks = scales::pretty_breaks(n = 4)) + geom_vline(xintercept = 29, color = 'orange', lty = 2) + labs(title = NULL, x = NULL, y = "Density") + theme_classic() + theme(axis.text = element_text(size = 15)) + theme(axis.title = element_text(size = 20)) # save output ggsave( "FigureS19A.jpeg", path = outdir, width = 8, height = 6, units = "in", dpi = 300 ) #--panel B # create subset of wbgtmax dataset for just India wbgtmax_ia <- wbgtmax_mindate %>% filter(grepl("^IA", dhsid)) # density plot wbgtmax_ia %>% ggplot(aes(x = wbgtmax)) + geom_density(color = 'black', fill = 'gray', alpha = 0.3) + scale_x_continuous(limits = c(0, 35), breaks = scales::pretty_breaks(n = 10)) + scale_y_continuous(limits = c(0, 0.155), breaks = scales::pretty_breaks(n = 4)) + geom_vline(xintercept = 29, color = 'orange', lty = 2) + labs(title = NULL, x = NULL, y = "Density") + theme_classic() + theme(axis.text = element_text(size = 15)) + theme(axis.title = element_text(size = 20)) # save output ggsave( "FigureS19B.jpeg", path = outdir, width = 8, height = 6, units = "in", dpi = 300 ) #--panel C # create subset of wbgtmax dataset for just Bangladesh wbgtmax_np <- wbgtmax_mindate %>% filter(grepl("^NP", dhsid)) # density plot wbgtmax_np %>% ggplot(aes(x = wbgtmax)) + geom_density(color = 'black', fill = 'gray', alpha = 0.3) + scale_x_continuous(limits = c(0, 35), breaks = scales::pretty_breaks(n = 10)) + scale_y_continuous(limits = c(0, 0.155), breaks = scales::pretty_breaks(n = 4)) + geom_vline(xintercept = 29, color = 'orange', lty = 2) + labs(title = NULL, x = NULL, y = "Density") + theme_classic() + theme(axis.text = element_text(size = 15)) + theme(axis.title = element_text(size = 20)) # save output ggsave( "FigureS19C.jpeg", path = outdir, width = 8, height = 6, units = "in", dpi = 300 ) #-------------------Figure S20-------------------------------------------------- #--summarize daily heat data to monthly # tmax tmax_mo <- tmax %>% drop_na(tmax) %>% mutate( date = case_when( day(date) >= 15 ~ paste0(str_sub(date, end = 8), 15), day(date) < 15 ~ paste0(str_sub((date %m-% months(1)), end = 8), 15), T ~ NA_character_), date = lubridate::ymd(date)) %>% group_by(dhsid, date) %>% summarize( dhsid = first(dhsid), mean_tmax = mean(tmax) ) %>% ungroup() %>% filter(date >= min(kid$triminus1_start), # minimum date dhsid %in% kid$dhsid) # drop any excess clusters # wbgtmax wbgtmax_mo <- wbgtmax %>% drop_na(wbgtmax) %>% mutate( date = case_when( day(date) >= 15 ~ paste0(str_sub(date, end = 8), 15), day(date) < 15 ~ paste0(str_sub((date %m-% months(1)), end = 8), 15), T ~ NA_character_), date = lubridate::ymd(date)) %>% group_by(dhsid, date) %>% summarize( dhsid = first(dhsid), mean_wbgtmax = mean(wbgtmax) ) %>% ungroup() %>% filter(date >= min(kid$triminus1_start), # minimum date dhsid %in% kid$dhsid) # drop any excess clusters # merge monthly data frames mo <- full_join(tmax_mo, wbgtmax_mo) #--scatter plot # for text annotation on plot cor_coef <- as.character(round( cor(mo$mean_wbgtmax, mo$mean_tmax, method = "pearson", use = "complete.obs"), digits = 3)) mo %>% ggplot(aes(x = mean_wbgtmax, y = mean_tmax, )) + geom_point(size = 0.75, color = "lightgray" ) + geom_smooth(size = 1, color = "red", ) + annotate("text", x=30, y=10, label= paste0("r = ", cor_coef), size = 5) + scale_x_continuous(breaks = scales::pretty_breaks()) + theme_minimal() + labs(x = "Mean WBGTmax (°C)", y = "Mean Tmax (°C)", ) # save output ggsave( "FigureS20.jpeg", path = outdir, width = 6.5, height = 5, units = "in", dpi = 300 ) #-------------------Table S1---------------------------------------------------- # generate data frame for summary table sum_table <- kid %>% filter(is.na(haz) == F) %>% filter(kidalive == 1) %>% filter(is.na(imptoilet) == F) %>% dplyr::select(haz, stunt, sevstunt, sex_lbl, twin_lbl, kidagemo, birthorder, age, birthplace, edu_lbl, momstunt, parity, marstat_lbl, age5yr_lbl, religion_lbl, imptoilet, urban, triminus1_t35_tmax:tri4_t35_tmax, triminus1_t29_tmin:tri4_t29_tmin, triminus1_t29_wbgtmax:tri4_t29_wbgtmax, triminus1_t29_wbgtmax_2050_SSP585:tri4_t29_wbgtmax_2050_SSP585, triminus1_t29_wbgtmax_2050_SSP245:tri4_t29_wbgtmax_2050_SSP245, triminus1_t35_tmax_2050_SSP585:tri4_t35_tmax_2050_SSP585, triminus1_t35_tmax_2050_SSP245:tri4_t35_tmax_2050_SSP245) %>% mutate( twin = case_when( twin_lbl == "Single birth/not a twin" ~ 0, twin_lbl == NA_character_ ~ NA_real_, T ~ 1), girl = case_when( sex_lbl == "Female" ~ 1, sex_lbl == "Male" ~ 0, T ~ NA_real_), bornclinic = case_when( birthplace == "Clinic" ~ 1, birthplace == NA_character_ ~ NA_real_, T ~ 0), bornhome = case_when( birthplace == "Home" ~ 1, birthplace == NA_character_ ~ NA_real_, T ~ 0), bornother = case_when( birthplace == "Other" ~ 1, birthplace == NA_character_ ~ NA_real_, T ~ 0), married = case_when( marstat_lbl == "In union" ~ 1, marstat_lbl == NA_character_ ~ NA_real_, T ~ 0), edu0 = case_when( edu_lbl == "No education" ~ 1, edu_lbl == NA_character_ ~ NA_real_, T ~ 0), edu1 = case_when( edu_lbl == "Primary" ~ 1, edu_lbl == NA_character_ ~ NA_real_, T ~ 0), edu2 = case_when( edu_lbl == "Secondary" ~ 1, edu_lbl == NA_character_ ~ NA_real_, T ~ 0), edu3 = case_when( edu_lbl == "Higher" ~ 1, edu_lbl == NA_character_ ~ NA_real_, T ~ 0), christian = case_when( religion_lbl == "Christian" ~ 1, religion_lbl == NA_character_ ~ NA_real_, T ~ 0), hindu = case_when( religion_lbl == "Hindu" ~ 1, religion_lbl == NA_character_ ~ NA_real_, T ~ 0), muslim = case_when( religion_lbl == "Muslim" ~ 1, religion_lbl == NA_character_ ~ NA_real_, T ~ 0), buddhist = case_when( religion_lbl == "Buddhist" ~ 1, religion_lbl == NA_character_ ~ NA_real_, T ~ 0), sikh = case_when( religion_lbl == "Sikh" ~ 1, religion_lbl == NA_character_ ~ NA_real_, T ~ 0), jain = case_when( religion_lbl == "Jain" ~ 1, religion_lbl == NA_character_ ~ NA_real_, T ~ 0), imptoilet = case_when( imptoilet == "NA" ~ NA_real_, T ~ imptoilet), birthorder = as.numeric(birthorder), kidagemo = as.numeric(kidagemo), parity = as.numeric(parity), ) # output table -- then cleaned up in latex modelsummary::datasummary_skim(sum_table, fmt = 2, output = "latex") #-------------------Table S2---------------------------------------------------- # fit regression mod_tmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + csw(tri3_t35_tmax_z, tri2_t35_tmax_z, tri1_t35_tmax_z, tri0_t35_tmax_z) | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") # generate table - clean up in latex etable(mod_tmax, digits = 2, tex = T, drop = c("birthorder*", "parity*", "religion*", "marstat*", "kidtwin*")) #-------------------Table S3---------------------------------------------------- # fit regression mod_tmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + csw(tri3_t35_tmax, tri2_t35_tmax, tri1_t35_tmax, tri0_t35_tmax) | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") # generate table - clean up in latex etable(mod_tmax, digits = 2, tex = T, drop = c("birthorder*", "parity*", "religion*", "marstat*", "kidtwin*")) #-------------------Table S4---------------------------------------------------- # fit regression mod_wbgtmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + csw(tri3_t29_wbgtmax_z, tri2_t29_wbgtmax_z, tri1_t29_wbgtmax_z, tri0_t29_wbgtmax_z) | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") # generate table - clean up in latex etable(mod_wbgtmax, digits = 2, tex = T, drop = c("birthorder*", "parity*", "religion*", "marstat*", "kidtwin*")) #-------------------Table S5---------------------------------------------------- # fit regression mod_wbgtmax <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + csw(tri3_t29_wbgtmax, tri2_t29_wbgtmax, tri1_t29_wbgtmax, tri0_t29_wbgtmax) | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") # generate table - clean up in latex etable(mod_wbgtmax, digits = 2, tex = T, drop = c("birthorder*", "parity*", "religion*", "marstat*", "kidtwin*")) #-------------------Table S6---------------------------------------------------- # fit regression mod_both <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + csw(tri3_t35_tmax_z, tri2_t35_tmax_z, tri1_t35_tmax_z, tri0_t35_tmax_z, tri3_t29_wbgtmax_z, tri2_t29_wbgtmax_z, tri1_t29_wbgtmax_z, tri0_t29_wbgtmax_z) | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") # generate table - clean up in latex etable(mod_both, digits = 2, tex = T, drop = c("birthorder*", "parity*", "religion*", "marstat*", "kidtwin*")) #-------------------Table S7---------------------------------------------------- # fit regression mod_both <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + csw(tri3_t35_tmax, tri2_t35_tmax, tri1_t35_tmax, tri0_t35_tmax, tri3_t29_wbgtmax, tri2_t29_wbgtmax, tri1_t29_wbgtmax, tri0_t29_wbgtmax) | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") # generate table - clean up in latex etable(mod_both, digits = 2, tex = T, drop = c("birthorder*", "parity*", "religion*", "marstat*", "kidtwin*")) #-------------------Table S8---------------------------------------------------- # fit regression mod_tmax_hhid <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + csw(tri3_t35_tmax, tri2_t35_tmax, tri1_t35_tmax, tri0_t35_tmax) | kidagemo + kidbirthmo + momage + surveymo + hhid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "singleton") # generate table - clean up in latex etable(mod_tmax_hhid, digits = 2, tex = T, drop = c("birthorder*", "parity*", "religion*", "marstat*", "kidtwin*")) #-------------------Table S9---------------------------------------------------- # fit regression mod_wbgtmax_hhid <- feols(haz ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + csw(tri3_t29_wbgtmax, tri2_t29_wbgtmax, tri1_t29_wbgtmax, tri0_t29_wbgtmax) | kidagemo + kidbirthmo + momage + surveymo + hhid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "singleton") # generate table - clean up in latex etable(mod_wbgtmax_hhid, digits = 2, tex = T, drop = c("birthorder*", "parity*", "religion*", "marstat*", "kidtwin*")) #-------------------Table S10--------------------------------------------------- # fit regression mod_tmax <- feols(stunt ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + csw(tri3_t35_tmax, tri2_t35_tmax, tri1_t35_tmax, tri0_t35_tmax) | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") # generate table - clean up in latex etable(mod_tmax, digits = 2, tex = T, drop = c("birthorder*", "parity*", "religion*", "marstat*", "kidtwin*")) #-------------------Table S11--------------------------------------------------- # fit regression mod_wbgtmax <- feols(stunt ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + csw(tri3_t29_wbgtmax, tri2_t29_wbgtmax, tri1_t29_wbgtmax, tri0_t29_wbgtmax) | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") # generate table - clean up in latex etable(mod_wbgtmax, digits = 2, tex = T, drop = c("birthorder*", "parity*", "religion*", "marstat*", "kidtwin*")) #-------------------Table S12--------------------------------------------------- # fit regression mod_both <- feols(stunt ~ kidsex + kidtwin + birthorder + birthplace + momedu + parity + marstat + imptoilet + religion + csw(tri3_t35_tmax, tri2_t35_tmax, tri1_t35_tmax, tri0_t35_tmax, tri3_t29_wbgtmax, tri2_t29_wbgtmax, tri1_t29_wbgtmax, tri0_t29_wbgtmax) | kidagemo + kidbirthmo + momage + surveymo + dhsid + admin1^surveyyr, weights = kid$perweight, data = kid, cluster = "dhsid", fixef.rm = "infinite_coef") # generate table - clean up in latex etable(mod_both, digits = 2, tex = T, drop = c("birthorder*", "parity*", "religion*", "marstat*", "kidtwin*"))