01-synchrony_plots.Rmd

---
title: "Co-occurrence and spatiotemporal distribution of malaria and dengue in Loreto, Peru - Script for analysis and plots"
author: "Paloma M. Carcamo, Jesus M. Quispe, Gabriel Carrasco-Escobar"
output: html_document
---

```{r}
if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, RColorBrewer, biscale, cowplot, patchwork, sf, pals, broom, metR, colorspace, spdep, gt, reactable)
```

### Load data

```{r}
malaria_raw <- read_csv("data/malaria.csv")
dengue_raw <- read.csv("data/dengue.csv", sep = ";", encoding="UTF-8")
pop_raw <- read_rds("data/district_pop_1990-2022.rds")
data(Peru, package = "innovar")
```

### Format data

```{r}
# Aggregate malaria data by week
malaria <- malaria_raw |> 
  mutate(week_start =  floor_date(date_not, unit = "week"),
         ubigeo = as.character(ubigeo),
         year = epiyear(week_start),
         week = epiweek(week_start)) |> 
  group_by(year, week, week_start, ubigeo) |> 
  summarise(cases_m = n())

# Aggregate dengue data by week
dengue <- dengue_raw |> 
  mutate(date_not = as.Date(date_not, format = "%d/%m/%Y"),
         week_start =  floor_date(date_not, unit = "week"),
         year = epiyear(week_start),
         week = epiweek(week_start),
         ubigeo = as.character(ubigeo)) |> 
  group_by(year, week, week_start, ubigeo) |> 
  summarise(cases_d = n())

# Create grid with all weeks (2000-2021) and ubigeos (Loreto, per 2017 census)
all_weeks <- data.frame(date = seq(as.Date("2000-01-02"), as.Date("2021-12-31"), by = "days")) |> 
  mutate(week_start = floor_date(date, unit = "week")) |> 
  group_by(week_start) |> 
  summarise() |> 
  mutate(year = epiyear(week_start),
         week = epiweek(week_start))

ubi <- data.frame(ubigeo = unique(Peru[Peru$dep == "LORETO",]$ubigeo))

grid <- cross_join(all_weeks, ubi)

# Process population dataset
  # Using linear splines (na.approx) to estimate population in districts created during the study period
  # (cubic splines create negative population values)
pop <- pop_raw |> 
  filter(year %in% c(2000:2021),
         ubigeo %in% ubi$ubigeo) |> 
  group_by(year) |> 
  complete(ubigeo = as.character(ubi$ubigeo)) |> 
  group_by(ubigeo) |> 
  mutate(population = zoo::na.approx(population, rule = 2, na.rm = FALSE))

# Join dengue, malaria and population datasets
  # By week
data_week <- grid |> 
  left_join(malaria, by = c("week_start", "year", "week", "ubigeo")) |> 
  left_join(dengue, by = c("week_start", "year", "week", "ubigeo")) |> 
  mutate(cases_m = if_else(is.na(cases_m), 0, cases_m),
         cases_d = if_else(is.na(cases_d), 0, cases_d)) |> 
  left_join(pop, by = c("year", "ubigeo")) |> 
  mutate(inc_m = cases_m/population*1000,
         inc_d = cases_d/population*1000) |> 
  left_join(Peru[Peru$dep == "LORETO", c("ubigeo", "distr")] |> st_drop_geometry(), by = "ubigeo") |> 
  arrange(week_start,distr)

  # By year
data_year <- data_week |> 
  group_by(year, ubigeo) |> 
  summarise(cases_m = sum(cases_m),
            cases_d = sum(cases_d),
            population = mean(population),
            inc_m = cases_m/population*1000,
            inc_d = cases_d/population*1000) |> 
  left_join(Peru[Peru$dep == "LORETO", c("ubigeo", "distr")] |> st_drop_geometry(), by = "ubigeo") |> 
  ungroup() |> 
  mutate(cut_m = cut(log(inc_m + 1), breaks = cartography::getBreaks(log(inc_m + 1), nclass = 3, method = "fisher"), labels = c(1:3), include.lowest = TRUE),
         cut_d = cut(log(inc_d + 1), breaks = cartography::getBreaks(log(inc_d + 1), nclass = 3, method = "fisher"), labels = c(1:3), include.lowest = TRUE)) 
```


### Time-lagged cross-correlation

```{r}
# Format data for TLCC
data_week_var <- data_week |> 
  ungroup() |> 
  group_by(distr) |> 
  mutate(cases_m_var = (cases_m/lag(cases_m) - 1) * 100,
         cases_d_var = (cases_d/lag(cases_d) - 1) * 100) |> 
  mutate(cases_m_var = if_else(is.finite(cases_m_var), cases_m_var, 0),
         cases_d_var = if_else(is.finite(cases_d_var), cases_d_var, 0),
         cases_m_var = if_else(is.nan(cases_m_var), 0, cases_m_var),
         cases_d_var = if_else(is.nan(cases_d_var), 0, cases_d_var)) |> 
  mutate(cases_m_var_abs = (cases_m - lag(cases_m)),
         cases_d_var_abs = (cases_d - lag(cases_d))) |> 
  mutate(cases_m_var_abs = if_else(is.finite(cases_m_var_abs), cases_m_var_abs, 0),
         cases_d_var_abs = if_else(is.finite(cases_d_var_abs), cases_d_var_abs, 0),
         cases_m_var_abs = if_else(is.nan(cases_m_var_abs), 0, cases_m_var_abs),
         cases_d_var_abs = if_else(is.nan(cases_d_var_abs), 0, cases_d_var_abs)) |> 
  filter(!is.na(cases_m_var)) |> 
  filter(!is.na(cases_d_var)) |> 
  filter(!is.na(cases_m_var_abs)) |> 
  filter(!is.na(cases_d_var_abs))

# write_rds(data_week_var, "data/dengue-malaria.rds")

# Create TLCC function for lags of 2 years (104 weeks)
cross_corr <- function(df, lag = 104) {
    tidy(ccf(x = df$cases_d_var,
             y = df$cases_m_var,
             lag.max = lag,
             plot = FALSE))
}

by_distr <- data_week_var |> 
  ungroup() |> 
  select(distr, week_start, 
         cases_m, cases_d, 
         cases_m_var, cases_d_var,
         cases_m_var_abs, cases_d_var_abs) |> 
  group_by(distr) |> 
  nest()

# Calculate TLCC coefficients
by_distr2 <- by_distr |> 
  mutate(crosscorr = purrr::map(data, cross_corr)) 

ccfs <- unnest(by_distr2, crosscorr)

# Create db with maximum coefficients
lines_vert <- ccfs |> 
  select(distr, lag, acf)|> 
  group_by(distr) |> 
  mutate(lag = if_else(is.na(acf), NA, lag)) |> 
  slice_max(order_by = abs(acf)) |> 
  mutate(acf = round(acf,2)) |> 
  unique() |> 
  mutate(district = case_when(distr == "TENIENTE MANUEL CLAVERO" ~ "T. MANUEL CLAVERO",
                              distr == "TENIENTE CESAR LOPEZ ROJAS" ~ "T. CESAR LOPEZ ROJAS",
                              .default = distr))

# Summary measures for coefficients and lags

## Positive coefficients
mean(lines_vert$acf[lines_vert$acf > 0], na.rm = TRUE)
sd(lines_vert$acf[lines_vert$acf > 0], na.rm = TRUE)

mean(lines_vert$lag[lines_vert$acf > 0], na.rm = TRUE)
sd(lines_vert$lag[lines_vert$acf > 0], na.rm = TRUE)

## Negative coefficients
mean(lines_vert$acf[lines_vert$acf < 0], na.rm = TRUE)
sd(lines_vert$acf[lines_vert$acf < 0], na.rm = TRUE)

mean(lines_vert$lag[lines_vert$acf < 0], na.rm = TRUE)
sd(lines_vert$lag[lines_vert$acf < 0], na.rm = TRUE)
```

### Figure 1A: Districts in Loreto

```{r}
distr <- Peru |> 
  filter(dep.code == 16) |> 
  mutate(Province = str_to_title(prov),
         District = str_to_title(distr),
         index = 1:53)

fig_1a <- distr |> 
    ggplot() +
    geom_sf(aes(fill = Province), linewidth = 0.5, color = "black") +
    # ggrepel::geom_label_repel(aes(label = index, geometry = geometry),
    #                           label.size = NA,
    #                           label.padding = 0.05,
    #                           na.rm = TRUE,
    #                           fill = alpha("white", 0.6),
    #                           stat = "sf_coordinates", 
    #                           max.overlaps = 100,
    #                           size = 3,
    #                           segment.linetype = 3) +
  scale_fill_manual(values = as.vector(brewer.seqseq2(8))) +
  labs(title = "A") +
  theme_void() +
  theme(legend.position = c(0.99, 0.2),
        legend.title = element_text(hjust = 0.5, size = 9),
        legend.text = element_text(size=7)) +
  guides(fill = guide_legend(title.position = "top", 
                             title.hjust = 0.5,
                             ncol = 2))

```

### Figure 1B: Incidence line plots

```{r}
fig_1b <- data_week |> 
  mutate(month = month(week_start)) |>
  group_by(year, month) |> 
  summarise(week_start = min(week_start),
            Dengue = sum(cases_d),
            Malaria = sum(cases_m)) |> 
  ggplot(aes(x = week_start)) +
  geom_line(aes(y = Malaria, color = "Malaria"), size = 1) +
  geom_line(aes(y = Dengue, color = "Dengue"), size = 1) +
  geom_point(aes(y = Malaria, color = "Malaria", shape = "Malaria"), shape = 18, size = 1.2) +
  geom_point(aes(y = Dengue, color = "Dengue", shape = "Dengue"), shape = 18, size = 1.2) +
  scale_color_discrete(type = c(brewer.pal(3, "Blues")[3], brewer.pal(3, "Reds")[3])) +
  labs(x = "", y = "Monthly reported cases", color = "", title = "B") +
  theme_classic() +
  theme(legend.position = c(0.9, 0.85),
        legend.title = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        text = element_text(size = 16),
        axis.line=element_line(size=0.8),
        panel.grid.major.y = element_line(),
        panel.grid.minor.y = element_line()) +
  guides(color = guide_legend(title = NULL))
```

### Figure 1C: Co-occurrence maps

```{r}
data <- bi_class(data_year, x = cut_m, y = cut_d, dim = 3)

fig_1c <- Peru |> 
  right_join(data, by = "ubigeo") |> 
  filter(year %in% c(2000, 2005, 2010, 2015, 2020)) |> 
  ggplot(lwd = 0.3) +
  geom_sf(aes(fill = bi_class), color = "black", linewidth = 0.5, size = 0.01, show.legend = FALSE) +
  bi_scale_fill(pal = "DkViolet", dim = 3) +
  labs(title = "C") +
  facet_wrap(~ year, ncol = 5) +
  theme_void() +
  theme(plot.title = element_text(margin = margin(t = 40, b = -12))) 

legend <- bi_legend(pal = "DkViolet",
                    dim = 3,
                    xlab = "Malaria",
                    ylab = "Dengue")
```

### Figure 1

```{r}
layout <- "AAAAAABBBBBB
          AAAAAABBBBBB
          AAAAAABBBBBB
          AAAAAABBBBBB
          CCCCCCCCCCCD
          CCCCCCCCCCCD"

(fig1 <- fig_1a + wrap_elements(full = fig_1b) + fig_1c + legend + plot_layout(design = layout) & theme(text = element_text(size = 18)))

# ggsave("plots/figure1.png", plot = fig1, width = 14, height = 9)
# ggsave("plots/figure1.eps", plot = fig1, device = "eps", width = 14, height = 9, family = "sans")
```

### Figure 2: Results of TLCC for eight selected districts

```{r}
lines_vert2 <- lines_vert |> 
  filter(distr %in% c("TAPICHE", "ALTO TAPICHE", "YAGUAS", "YAQUERANA",
                      "NAUTA", "TENIENTE MANUEL CLAVERO", "ANDOAS", "YURIMAGUAS"))

tlcc_8d_df <- ccfs |> 
  filter(distr %in% c("TAPICHE", "ALTO TAPICHE", "YAGUAS", "YAQUERANA",
                      "NAUTA", "TENIENTE MANUEL CLAVERO", "ANDOAS", "YURIMAGUAS")) |>
  mutate(group = case_when(distr  %in% c("TAPICHE", "ALTO TAPICHE") ~ "a",
                           distr  %in% c("YAGUAS", "YAQUERANA") ~ "b",
                           distr  %in% c("NAUTA", "TENIENTE MANUEL CLAVERO") ~ "c",
                           distr  %in% c("ANDOAS", "YURIMAGUAS") ~ "d"),
         group = factor(group, levels = c("a", "b", "c", "d"), labels = c("Negative TLCC, \nnegative lag\n",
                                                                          "Negative TLCC, \npositive lag\n",
                                                                          "Positive TLCC, \nnegative lag\n",
                                                                          "Positive TLCC, \npositive lag\n")))

ggplot(tlcc_8d_df, aes(lag, acf, color = group)) +
  geom_point(alpha = 2 / 3, size=0.5) +
  geom_line() +
  geom_vline(data = lines_vert2, aes(xintercept = lag),
           linetype = 2, color = "tan2") +
  geom_label(data = lines_vert2, aes(label = 
                                       paste("\u03C1", "=", acf, ",", "Lag", "=", lag), 
                                     x = 15, y = 0.5),
          fontface = "bold", size = 3.5, color = "tan2") +
  facet_wrap(~ factor(distr,
                      levels = c("TAPICHE", "ALTO TAPICHE", "YAGUAS", "YAQUERANA",
                      "NAUTA", "TENIENTE MANUEL CLAVERO", "ANDOAS", "YURIMAGUAS")), 
             ncol = 2, dir = "h") +
  scale_color_manual(values = brewer.pal(4, "Paired")) + 
  geom_segment(x = -10, xend = -50, y = -0.65, yend = -0.65,
               arrow = arrow(length = unit(5, "pt")),
               data = tlcc_8d_df |> filter(distr %in% c("ANDOAS", "YURIMAGUAS")),
               color = "#3182BD") +
  geom_segment(x = 10, xend = 50, y = -0.65, yend = -0.65,
               arrow = arrow(length = unit(5, "pt")),
               data = tlcc_8d_df |> filter(distr %in% c("ANDOAS", "YURIMAGUAS")),
               color = "#DE2D26") +
  geom_text(aes(label = "Dengue precedes\nmalaria"),
            hjust = 1,
            x = -55, y = -0.75,
            data = tlcc_8d_df |> filter(distr %in% c("ANDOAS")),
            color = "#3182BD",
            size = 3) +
  geom_text(aes(label = "Malaria precedes\ndengue"),
            hjust = 0,
            x = 55, y = -0.75,
            data = tlcc_8d_df |> filter(distr %in% c("YURIMAGUAS")),
            color = "#DE2D26",
            size = 3) +
  theme_bw() +
  labs(x = "Lag (weeks)", y = "TLCC coefficient", color = "") +
  theme(strip.background = element_rect(fill = "black"),
        strip.text = element_text(colour = 'white'),
        plot.margin = unit(c(0.2, 0.5, 0.5, 0.2), "inches")) +
  coord_cartesian(clip = "off")

# ggsave("plots/figure2.png", width = 7, height = 8)
# ggsave("plots/figure2.eps", device = "eps", width = 7, height = 8)
```

### Figure 3A: TLCC categories map

```{r}
fig_3a <- Peru |> 
  filter(dep == "LORETO") |> 
  left_join(lines_vert, by = "distr") |>
  mutate(category = if_else(acf > 0,
                            if_else(lag > 0,
                                   "Positive TLCC, \npositive lag",
                                   "Positive TLCC, \nnegative lag"),
                            if_else(lag > 0,
                                    "Negative TLCC, \npositive lag",
                                    "Negative TLCC, \nnegative lag"))) |> 
  ggplot() +
  geom_sf(aes(fill = category)) +
  labs(fill = "", title = "A") +
  scale_fill_discrete(type = brewer.pal(4, "Paired"),
                      labels = function(breaks) {breaks[is.na(breaks)] <- "Not calculated"; breaks}) +
  theme_void() +
  theme(legend.justification = "top",
        legend.text = element_text(size = 7))
```

### Figure 3A inset: Coefficient scatterplot

```{r}
fig_3a_inset <- lines_vert |> 
  mutate(category = if_else(acf > 0,
                            if_else(lag > 0,
                                   "Positive TLCC, \npositive lag",
                                   "Positive TLCC, \nnegative lag"),
                            if_else(lag > 0,
                                    "Negative TLCC, \npositive lag",
                                    "Negative TLCC, \nnegative lag"))) |>
  ggplot(aes(x = lag, y = acf, color = category)) +
  geom_point() +
  geom_hline(aes(yintercept = 0)) +
  geom_vline(aes(xintercept = 0)) +
  ggrepel::geom_label_repel(data = lines_vert |> 
                              mutate(district = str_to_title(district)) |> 
                              filter(abs(acf) > 0.2), 
                            aes(label = district),
                            label.size = NA,
                            min.segment.length = 0,
                            fill = alpha(c("white"), 0.3),
                            size = 3,
                            color = "black") +
  labs(x = "Lags (weeks)",
       y = "TLCC coefficient",
       color = "") +
  scale_color_discrete(type = brewer.pal(4, "Paired")) +
  theme_classic() +
  theme(legend.position = "none",
        axis.line = element_blank(),
        panel.background = element_blank())
```

### Figure 3B: Local Moran's test for TLCC

```{r}
clor_map_corr <- Peru |> 
  filter(dep == "LORETO") |> 
  left_join(lines_vert, by = "distr")|> 
  replace_na(list(lag = 0, acf = 0))

nb <- poly2nb(clor_map_corr, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

# Global Moran's test
moran.t <- moran.test(clor_map_corr$acf, lw)
glance_htest <- function(ht) c(ht$estimate,
                               "Std deviate" = unname(ht$statistic),
                               "p.value" = unname(ht$p.value))

moran.glance <- moran.t |> 
  glance_htest() |> 
  data.frame()

reactable(moran.glance)
```

```{r}
lcl_moran_acf <- localmoran(clor_map_corr$acf, lw)
lcl_moran_acf_full <- cbind(clor_map_corr, lcl_moran_acf) |> 
  mutate(stat = if_else(Pr.z....E.Ii.. < 0.05, Ii, NA))

fig_3b_prelim <- lcl_moran_acf_full |> 
  ggplot() +
  geom_sf(aes(fill = stat)) +
  scale_fill_gradient2(low = "#d7191c", mid = "#ffffbf", high = "#1a9641", na.value = "white", limits = c(-1.2, 1.2)) +
  theme_void() +
  labs(title = "B", fill = "Local Moran \nstatistic") +
  theme(legend.text = element_text(size = 8))

fig_3b <- fig_3b_prelim + theme(legend.position = "none")
```

### Figure 3C: Local Moran's test for lags

```{r}
lcl_moran_lag <- localmoran(clor_map_corr$lag, lw)
lcl_moran_lag_full <- cbind(clor_map_corr, lcl_moran_lag) |> 
  mutate(stat = if_else(Pr.z....E.Ii.. < 0.05, Ii, NA))

fig_3c <- lcl_moran_lag_full |> 
  ggplot() +
  geom_sf(aes(fill = stat)) +
  scale_fill_gradient2(low = "#d7191c", mid = "#ffffbf", high = "#1a9641", na.value = "white", limits = c(-1.2, 1.2)) +
  theme_void() +
  labs(title = "C") +
  theme(legend.position = "none")
```

### Figure 3 legend

```{r}
fig_3legend <- get_legend(fig_3b_prelim)
```

### Figure 3

```{r}
fig_3a_full <- fig_3a + inset_element(fig_3a_inset, left = 0.7, bottom = 0, right = 1.4, top = 0.45)

layout2 <- "AAAAAAA#BBBBB##
           AAAAAAA#BBBBB##
           AAAAAAA#BBBBBDD
           AAAAAAA#CCCCCDD
           AAAAAAA#CCCCC##
           AAAAAAA#CCCCC##"

(fig3 <- fig_3a_full + fig_3b + fig_3c + fig_3legend + plot_layout(design = layout2))

ggsave("plots/figure3.png", width = 10, height = 6)
ggsave("plots/figure3.eps", device = "eps", plot = fig3, width = 10, height = 6)
```

### Figure S1: Co-occurrence maps for all years

```{r}
fig_s1 <- Peru |> 
  right_join(data, by = "ubigeo") |> 
  ggplot() +
  geom_sf(aes(fill = bi_class), color = "black", size = 0.01, show.legend = FALSE) +
  bi_scale_fill(pal = "DkViolet", dim = 3) +
  facet_wrap(~ year) +
  theme_void() +
  theme(plot.title = element_text(margin = margin(t = 40, b = -12)))


fig_s1 + inset_element(legend, left = 0, bottom = 0.05, right = 1.7, top = 0.2)

# ggsave("plots/figureS1.png", height = 10, width = 8)
# ggsave("plots/figureS1.eps", height = 10,  device = "eps", width = 8, family = "sans")
```

### Figure S2: Dengue incidence

```{r}
(fig_s2 <- Peru |> 
  right_join(data_year, by = "ubigeo") |> 
  ggplot() +
  geom_sf(aes(fill = inc_d), lwd = 0.2, color = "black") +
  scale_fill_gradientn(name = "Dengue incidence \nper 1000",
                       colours = brewer.pal(9, "Blues"),
                       trans = "log1p", 
                       breaks = c(0, 3, 30), 
                       labels = c(0, 3, 30)) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        panel.background = element_blank(),
        strip.background = element_rect(fill = "black"),
        strip.text = element_text(colour = 'white')) +
  facet_wrap(~year, ncol = 6))

ggsave("plots/figureS2.png", width = 14, height = 10)
```

### Figure S3: Malaria incidence

```{r}
(fig_s3 <- Peru |> 
  right_join(data_year, by = "ubigeo") |> 
  ggplot() +
  geom_sf(aes(fill = inc_m), lwd = 0.2, color = "black") +
  scale_fill_gradientn(name = "Malaria incidence \nper 1000",
                       colours = brewer.pal(9, "Reds"),
                       trans = "log1p", 
                       breaks = c(0, 6, 60, 600), 
                       labels = c(0, 6, 60, 600)) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        panel.background = element_blank(),
        strip.background = element_rect(fill = "black"),
        strip.text = element_text(colour = 'white')) +
  facet_wrap(~year, ncol = 6))

ggsave("plots/figureS3.png", width = 14, height = 10)
```

### Figure S4: Results of TLCC for all districts

```{r}
bottom_row = c("VARGAS GUERRA", "YAGUAS", "YAQUERANA", "YAVARI", "YURIMAGUAS", "TORRES CAUSANA", "TROMPETEROS", "URARINAS")

ccfs2 <- ccfs |> 
  mutate(district = case_when(distr == "TENIENTE MANUEL CLAVERO" ~ "T. MANUEL CLAVERO",
                              distr == "TENIENTE CESAR LOPEZ ROJAS" ~ "T. CESAR LOPEZ ROJAS",
                              .default = distr))

ggplot(ccfs2, aes(lag, acf)) +
  geom_point(alpha = 2 / 3, size = 0.5) +
  geom_line(color = "#303030") +
  geom_vline(data = lines_vert, 
             aes(xintercept = lag),
             linetype = 2,
             color = "tan2") +
  geom_label(data = lines_vert |> filter(!is.na(lag)),
             aes(label = paste("\u03C1","=",acf,",","Lag","=",lag), x = 15, y = 0.5),
             size = 3.5,
             color = "tan2",
             fontface = "bold")  +
  facet_wrap(~ district) +
  geom_segment(x = -10, xend = -90, y = -0.65, yend = -0.65,
               arrow = arrow(length = unit(5, "pt")),
               data = ccfs2 |> filter(distr %in% bottom_row),
               color = "#3182BD") +
  geom_segment(x = 10, xend = 90, y = -0.65, yend = -0.65,
               arrow = arrow(length = unit(5, "pt")),
               data = ccfs2 |> filter(distr %in% bottom_row),
               color = "#DE2D26") +
  geom_text(data = ccfs2 |> filter(distr %in% bottom_row) |> filter(distr != "YURIMAGUAS"),
            aes(label = "Dengue\nprecedes\nmalaria"),
            hjust = 1,
            x = -10, y = -0.85,
            color = "#3182BD",
            size = 3) +
  geom_text(data = ccfs2 |> filter(distr %in% bottom_row) |> filter(distr != "YAVARI"),
            aes(label = "Malaria\nprecedes\ndengue"),
            hjust = 0,
            x = 10, y = -0.85,
            color = "#DE2D26",
            size = 3) +
  theme_bw() +
  theme(legend.position = "none") +
  labs(x = "\n\nLag (weeks)", y = "TLCC coefficient") +
  theme(strip.background = element_rect(fill = "black"),
        strip.text = element_text(colour = 'white',
                                  face = "bold"),
        plot.margin = unit(c(0.2, 0.5, 0.5, 0.2), "inches")) +
  coord_cartesian(clip = "off")

# ggsave("plots/figureS4.png", width = 15, height = 15)
# ggsave("plots/figureS4.eps", width = 15, height = 15)
```

### Figure S5a: District-level TLCC coefficient map

```{r}
fig_s5a <- clor_map_corr |> 
  ggplot() +
  geom_sf(aes(fill = acf)) +
  scale_fill_divergent(midpoint = 0,
                       name = "TLCC\ncoefficient") +
  theme_void()
```

### Figure S5b: District-level lag map

```{r}
fig_s5b <- clor_map_corr |> 
  ggplot() +
  geom_sf(aes(fill = lag)) +
  scale_fill_divergent(midpoint = 0,
                       name = "Lag (weeks)") +
  theme_void()
```

### Figure S5C: Density plot for TLCC coefficients

```{r}
(fig_s5c <- clor_map_corr |> 
  ggplot(aes(acf)) +
  geom_histogram(aes(y = ..density..)) +
  geom_density(color = "darkred", linewidth = 0.8) +
  theme_bw() +
  labs(x = "TLCC coefficient", y = "Density") +
  scale_y_continuous(expand = c(0,0)) +
  theme(panel.grid = element_blank()))
```

### Figure S5D: Density plot for lags

```{r}
(fig_s5d <- clor_map_corr |> 
  ggplot(aes(lag)) +
  geom_histogram(aes(y = ..density..)) +
   geom_density(color = "darkred", linewidth = 0.8) +
  theme_bw() +
  labs(x = "Lag (weeks)", y = "Density") +
  scale_y_continuous(expand = c(0,0)) +
  theme(panel.grid = element_blank()))
```

### Figure S5

```{r}
layout3 <- "AAAAAAAABBBBBBBB
            AAAAAAAABBBBBBBB
            AAAAAAAABBBBBBBB
            AAAAAAAABBBBBBBB
            AAAAAAAABBBBBBBB
            AAAAAAAABBBBBBBB
            CCCCCCCCDDDDDDDD
            CCCCCCCCDDDDDDDD"

fig_s5a + fig_s5b + fig_s5c + fig_s5d + 
  plot_annotation(tag_levels = "A") +
  plot_layout(design = layout3)

# ggsave("plots/figureS5.png", height = 8, width = 12)
# ggsave("plots/figureS5.eps", height = 8, width = 12)
```

## Sensitivity analysis:

```{r}
# # TLCC with max lag of 6 months (25 weeks)
# by_distr_6m <- by_distr |> 
#   mutate(crosscorr = purrr::map2(data, 25, cross_corr))
# 
# ccfs_6m <- unnest(by_distr_6m, crosscorr)
# 
# max_ccfs_6m <- ccfs_6m |> 
#   select(distr, lag, acf) |> 
#   group_by(distr) |> 
#   mutate(lag = if_else(is.na(acf), NA, lag)) |> 
#   slice_max(order_by = abs(acf)) |> 
#   mutate(acf = round(acf, 2)) |> 
#   unique() |> 
#   mutate(window = "6 months")

# TLCC with max lag of 1 year (52 weeks)
by_distr_1y <- by_distr |> 
  mutate(crosscorr = purrr::map2(data, 52, cross_corr))

ccfs_1y <- unnest(by_distr_1y, crosscorr)

max_ccfs_1y <- ccfs_1y |> 
  select(distr, lag, acf) |> 
  group_by(distr) |> 
  mutate(lag = if_else(is.na(acf), NA, lag)) |> 
  slice_max(order_by = abs(acf)) |> 
  mutate(acf = round(acf, 2)) |> 
  unique() |> 
  mutate(window = "1 year")

# TLCC with max lag of 3 years (156 weeks)
by_distr_3y <- by_distr |> 
  mutate(crosscorr = purrr::map2(data, 156, cross_corr))

ccfs_3y <- unnest(by_distr_3y, crosscorr)

max_ccfs_3y <- ccfs_3y |> 
  select(distr, lag, acf) |> 
  group_by(distr) |> 
  mutate(lag = if_else(is.na(acf), NA, lag)) |> 
  slice_max(order_by = abs(acf)) |> 
  mutate(acf = round(acf, 2)) |> 
  unique() |> 
  mutate(window = "3 years")

# TLCC with max lag of 4 years (208 weeks)
by_distr_4y <- by_distr |> 
  mutate(crosscorr = purrr::map2(data, 208, cross_corr))

ccfs_4y <- unnest(by_distr_4y, crosscorr)

max_ccfs_4y <- ccfs_4y |> 
  select(distr, lag, acf) |> 
  group_by(distr) |> 
  mutate(lag = if_else(is.na(acf), NA, lag)) |> 
  slice_max(order_by = abs(acf)) |> 
  mutate(acf = round(acf, 2)) |> 
  unique() |> 
  mutate(window = "4 years")

# TLCC with max lag of 5 years (260 weeks)
by_distr_5y <- by_distr |> 
  mutate(crosscorr = purrr::map2(data, 260, cross_corr))

ccfs_5y <- unnest(by_distr_5y, crosscorr)

max_ccfs_5y <- ccfs_5y |> 
  select(distr, lag, acf) |> 
  group_by(distr) |> 
  mutate(lag = if_else(is.na(acf), NA, lag)) |> 
  slice_max(order_by = abs(acf)) |> 
  mutate(acf = round(acf, 2)) |> 
  unique() |> 
  mutate(window = "5 years")
```


```{r}
ccf_by_window <- lines_vert |> 
  mutate(window = "2 years") |> 
  select(-district) |> 
#  rbind(max_ccfs_6m) |> 
  rbind(max_ccfs_1y) |> 
  rbind(max_ccfs_3y) |> 
  rbind(max_ccfs_4y) |> 
  rbind(max_ccfs_5y)

map_ccf_by_window <- Peru |> 
  filter(dep == "LORETO") |> 
  left_join(ccf_by_window, by = "distr")|> 
  replace_na(list(lag = 0, acf = 0)) |> 
  mutate(window = factor(window, levels = c("1 year", "2 years", "3 years", "4 years", "5 years")))

coeff_windows <- map_ccf_by_window |> 
  ggplot() +
  geom_sf(aes(fill = acf)) +
  scale_fill_divergent(midpoint = 0,
                       name = "TLCC\ncoefficient") +
  theme_void() +
  facet_wrap(~window, ncol = 6)

lag_windows <- map_ccf_by_window |> 
  ggplot() +
  geom_sf(aes(fill = lag)) +
  scale_fill_divergent(midpoint = 0,
                       name = "Lag (weeks)") +
  theme_void() +
  facet_wrap(~window, ncol = 6)

coeff_windows / lag_windows

# ggsave("plots/sensitivity1.png", height = 8, width = 12)
# ggsave("plots/sensitivity1.eps", height = 8, width = 12)
```

## Figure S6: Sensitivity analysis of time windows for TLCC

```{r}
pos_tlcc <- c("ANDOAS", "CAPELO", "LAGUNAS", "MANSERICHE", "MAZAN", "PUNCHANA", "RAMON CASTILLA", "SAN JUAN BAUTISTA", "SAQUENA", "TENIENTE MANUEL CLAVERO", "TIGRE", "YAVARI", "YURIMAGUAS")
neg_tlcc <- c("ALTO TAPICHE", "CAHUAPANAS", "EMILIO SAN MARTIN", "JENARO HERRERA", "PARINARI", "PASTAZA", "SAN PABLO", "SANTA CRUZ", "SARAYACU", "SOPLIN", "TAPICHE", "YAGUAS", "YAQUERANA")

ccf_plot <- map_ccf_by_window |> 
  mutate(level = case_when(window == "1 year" ~ 1,
                           window == "2 years" ~ 2,
                           window == "3 years" ~ 3,
                           window == "4 years" ~ 4,
                           window == "5 years" ~ 5),
         bi_class = case_when(acf < 0 & lag < 0 ~ "1-1",
                              acf > 0 & lag < 0 ~ "2-1",
                              acf < 0 & lag > 0 ~ "1-2",
                              acf > 0 & lag > 0 ~ "2-2"),
         category = case_when(distr %in% pos_tlcc ~ "Consistently positive TLCC",
                              distr %in% neg_tlcc ~ "Consistently negative TLCC",
                              distr == "ROSA PANDURO" ~ "Not calculated",
                              .default = "Inconsistent TLCC or lags"))

ccf_plot |> 
  ggplot(aes(x = level, y = acf, color = distr)) +
  geom_hline(yintercept = 0) +
  geom_line(show.legend = FALSE) +
  facet_wrap(~prov) +
  labs(x = "Years in window", y = "TLCC coefficient") +
  theme_bw()

ccf_plot |> 
  ggplot(aes(x = level, y = lag, color = distr)) +
  geom_hline(yintercept = 0) +
  geom_line(show.legend = FALSE) +
  facet_wrap(~prov) +
  labs(x = "Years in window", y = "lag") +
  theme_bw()

pos_plot <- ccf_plot |> 
  mutate(`ACF (*1000)` = acf*1000,
         `Lag (weeks)` = lag,
         distr = if_else(distr == "TENIENTE MANUEL CLAVERO", "T MANUEL CLAVERO", distr)) |> 
  pivot_longer(cols = c("ACF (*1000)", "Lag (weeks)"), names_to = "cat", values_to = "value") |> 
  filter(category == "Consistently positive TLCC") |> 
  ggplot(aes(x = level, y = value, color = cat)) +
  geom_hline(yintercept = 0) +
  geom_line() +
  facet_wrap(~ distr, ncol = 5) +
  labs(x = "Years in window", y = "", color = "", title = "Consistently positive TLCC coefficient") +
  theme_bw() +
  theme(strip.background = element_rect(fill = "black"),
        strip.text = element_text(colour = 'white'))

neg_plot <- ccf_plot |> 
  mutate(`ACF (*1000)` = acf*1000,
         `Lag (weeks)` = lag) |> 
  pivot_longer(cols = c("ACF (*1000)", "Lag (weeks)"), names_to = "cat", values_to = "value") |> 
  filter(category == "Consistently negative TLCC") |> 
  ggplot(aes(x = level, y = value, color = cat)) +
  geom_hline(yintercept = 0) +
  geom_line() +
  facet_wrap(~ distr, ncol = 5) +
  labs(x = "Years in window", y = "", color = "", title = "Consistently negative TLCC coefficient") +
  theme_bw() +
  theme(strip.background = element_rect(fill = "black"),
        strip.text = element_text(colour = 'white'))

others_plot <- ccf_plot |> 
  mutate(`ACF (*1000)` = acf*1000,
         `Lag (weeks)` = lag,
         distr = if_else(distr == "TENIENTE CESAR LOPEZ ROJAS", "T LOPEZ ROJAS", distr)) |> 
  pivot_longer(cols = c("ACF (*1000)", "Lag (weeks)"), names_to = "cat", values_to = "value") |> 
  filter(category == "Inconsistent TLCC or lags") |> 
  ggplot(aes(x = level, y = value, color = cat)) +
  geom_hline(yintercept = 0) +
  geom_line() +
  facet_wrap(~ distr, ncol = 5) +
  labs(x = "Years in window", y = "", color = "", title = "Inconsistent TLCC coefficient or lags") +
  theme_bw() +
  theme(strip.background = element_rect(fill = "black"),
        strip.text = element_text(colour = 'white'))

layout4 <- "
AAAAAAAAAAAAAAAAAAA
AAAAAAAAAAAAAAAAAAA
AAAAAAAAAAAAAAAAAAA
BBBBBBBBBBBBBBBBBBB
BBBBBBBBBBBBBBBBBBB
BBBBBBBBBBBBBBBBBBB
CCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCC"

fig_s6 <- pos_plot / neg_plot / others_plot + plot_annotation(tag_levels = "A") + plot_layout(design = layout4)

# ggsave("plots/figureS6.png", plot = fig_s6, width = 10, height = 15)
# ggsave("plots/figureS6.eps", plot = fig_s6, width = 10, height = 15)
```

## Figure S7
```{r}
fig_s7 <- ccf_plot |> 
  group_by(geometry, category) |> 
  summarise() |> 
  ggplot(aes(fill = category)) +
  geom_sf() +
  labs(fill = "") +
  scale_fill_manual(values = c("#FC8D59", "#FFFFBF", "#91BFDB", "gray")) +
  theme_void() +
  theme(plot.background = element_rect(fill = "white", color = NA))

# ggsave("plots/figureS7.png", plot = fig_s7, width = 5, height = 5)
# ggsave("plots/figureS7.eps", plot = fig_s7, width = 5, height = 5)
```

## Figure S8
```{r}
lvs <- unique(ccf_plot$distr[order(ccf_plot$acf)])

plot_s8 <- ccf_plot |> 
  mutate(distr.ord = factor(distr, levels = lvs)) |> 
  ggplot(aes(x = level, y = distr.ord, fill = bi_class)) +
  geom_tile(show.legend = FALSE) +
  bi_scale_fill(pal = "GrPink2", dim = 2) +
  theme_bw() +
  labs(x = "Years in window", y = "", fill = "")

leg_s8 <- bi_legend(pal = "GrPink2",
                    dim = 2,
                    xlab = "TLCC-    TLCC+",
                    ylab = "Lag-     Lag+",
                    size = 10, arrows = FALSE)
layout5 <- "
AAAAAAA###
AAAAAAA#BB
AAAAAAA#BB
AAAAAAA###"

plot_s8 + leg_s8 +  plot_layout(design = layout5)

# ggsave("plots/figureS8.png", width = 8, height = 10)
# ggsave("plots/figureS8.eps", width = 8, height = 10, family = "sans")
```