01_descriptive_analysis.Rmd

---
title: "Descriptive analysis"
date: "`r Sys.Date()`"
output: 
  rmdformats::downcute:
    self_contained: true
    highlight: kate
    toc_depth: 3
    default_style: dark
    code_folding: hide
    code_download: true
    highlight_downlit: true
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
## Global options
knitr::opts_chunk$set(
  cache = TRUE,
  warnings = FALSE,
  messages = FALSE
)

library(tidyverse)
library(sf)
library(ggsflabel)
library(ggspatial)
library(patchwork)
source("03_functions/seropositivy_summarise.R")
conflicted::conflicts_prefer(dplyr::select())
conflicted::conflicts_prefer(dplyr::filter())
```

# Import data

```{r}
ffi_total_simply <- readRDS("01_data/processed/ffi_total.rds")
ffi_w.getis_simply <- readRDS("01_data/share/ffi_w.getis_simply.rds")

# ffi_total_simply <- ffi_total_simply %>% 
#   mutate(
#     distance_to_rh_hf_sampl = fct_recode(
#       distance_to_rh_hf_sampl,
#       "Proximate Distance" = "Proximate",
#       "Medium Distance" = "Moderate",
#       "Long Distance" = "Distant",
#       "Extreme Distance" = "Extra Distant"
#     )
#   )
```

```{r}
ffi_total_simply %>% 
  mutate(
    UR = case_match(
      ffi_is_community,
      c("URB. BAGAZAN", "PUEBLO LIBRE",
        "PP.JJ. SANTO CRISTO BAGAZAN") ~ "Urban",
      .default = "Rural"
    )
  ) %>% 
  count(distance_to_rh_hf_sampl, UR)
```


# Figures

## Figure 3

```{r}
sero_timeofmalaria_district_distance_rh_hf_sampling <- seropositivy_summarise(
    ffi_total_simply,
    "timeofmalaria",
    age_cat,
    distance_to_rh_hf_sampl
  ) %>%
  ungroup() %>% 
  mutate(
    distance_to_rh_hf_sampl = fct_relevel(distance_to_rh_hf_sampl, 
                                          "Proximate Distance",
                                          "Medium Distance",
                                          "Long Distance",
                                          "Extreme Distance")
  ) %>% 
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper, fill = exposure), alpha = 0.1) +
  facet_grid(
    vars(distance_to_rh_hf_sampl),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format()
  ) +
  labs(
    x = "Age",
    y = "Seropositivity Rate"
  ) +
  guides(
    color = guide_legend("Malaria"),
    fill = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  innovar::scale_fill_innova("npr") +
  theme_bw() +
  theme(
    plot.title = element_text(
      size = 14,
      hjust = 0.5,
      face = "bold"
    ),
    strip.text = element_text(
      size = 12,
      face = "bold"
    ),
    legend.title = element_text(
      size = 15,
      face = "bold"
    ),
    legend.text = element_text(size = 14),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_text(size = 12),
    axis.title = element_text(size = 13)
  )

sero_timeofmalaria_district_distance_rh_hf_sampling
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/Fig 3. Seropositivity according to time of exposure to malaria by district and distance category.png",
  sero_timeofmalaria_district_distance_rh_hf_sampling,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)

ggsave(
  "./02_output/plots/Fig 3. Seropositivity according to time of exposure to malaria by district and distance category.eps",
  sero_timeofmalaria_district_distance_rh_hf_sampling,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)

ggsave(
  "./02_output/plots/Fig 3. Seropositivity according to time of exposure to malaria by district and distance category.tiff",
  sero_timeofmalaria_district_distance_rh_hf_sampling,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

## Figure 3 - Alt

```{r}
sero_typeofmalaria_district_distance_rh_hf_sampling <- seropositivy_summarise(
    ffi_total_simply,
    "typeofmalaria",
    age_cat,
    distance_to_rh_hf_sampl
  ) %>%
  ungroup() %>% 
  mutate(
    distance_to_rh_hf_sampl = fct_relevel(distance_to_rh_hf_sampl, 
                                          "Proximate Distance",
                                          "Medium Distance",
                                          "Long Distance",
                                          "Extreme Distance"),
    exposure = case_when(
      exposure == "P. Falciparum Exposure" ~ "*P. Falciparum* Exposure",
      exposure == "P. Vivax Exposure" ~ "*P. Vivax* Exposure"
    )
  ) %>% 
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper, fill = exposure), alpha = 0.1) +
  facet_grid(
    vars(distance_to_rh_hf_sampl),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format()
  ) +
  labs(
    x = "Age",
    y = "Seropositivity Rate"
  ) +
  guides(
    color = guide_legend("Malaria"),
    fill = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  innovar::scale_fill_innova("npr") +
  theme_bw() +
  theme(
    plot.title = element_text(
      size = 14,
      hjust = 0.5,
      face = "bold"
    ),
    strip.text = element_text(
      size = 12,
      face = "bold"
    ),
    legend.title = element_text(
      size = 15,
      face = "bold"
    ),
    legend.text = ggtext::element_markdown(
      size = 14
    ),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_text(size = 12),
    axis.title = element_text(size = 13)
  )

sero_typeofmalaria_district_distance_rh_hf_sampling
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/Fig 3. Seropositivity according to type of exposure to malaria by district and distance category.png",
  sero_typeofmalaria_district_distance_rh_hf_sampling,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)

ggsave(
  "./02_output/plots/Fig 3. Seropositivity according to type of exposure to malaria by district and distance category.eps",
  sero_typeofmalaria_district_distance_rh_hf_sampling,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)

ggsave(
  "./02_output/plots/Fig 3. Seropositivity according to type of exposure to malaria by district and distance category.tiff",
  sero_typeofmalaria_district_distance_rh_hf_sampling,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

## Figure 3 - Alt 2

```{r}
sero_4malaria_district_distance_rh_hf_sampling <- seropositivy_summarise(
    ffi_total_simply,
    "4malaria",
    age_cat,
    distance_to_rh_hf_sampl
  ) %>%
  ungroup() %>% 
  mutate(
    distance_to_rh_hf_sampl = fct_relevel(distance_to_rh_hf_sampl, 
                                          "Proximate Distance",
                                          "Medium Distance",
                                          "Long Distance",
                                          "Extreme Distance"),
    malaria = case_when(
      malaria == "Recent P. Falciparum" ~ "Recent *P. falciparum*",
      malaria == "Historical P. Falciparum" ~ "Historical *P. falciparum*",
      malaria == "Recent P. Vivax" ~ "Recent *P. vivax*",
      malaria == "Historical P. Vivax" ~ "Historical *P. vivax*"
    ),
  ) %>% 
  ggplot(
    aes(
      x = age_cat,
      y = result_malaria,
      color = malaria,
      group = malaria
    )
  ) +
  geom_point() +
  geom_line() +
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper, fill = malaria), alpha = 0.1) +
  facet_grid(
    vars(distance_to_rh_hf_sampl),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format()
  ) +
  labs(
    x = "Age",
    y = "Seropositivity Rate"
  ) +
  guides(
    color = guide_legend("Malaria"),
    fill = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  innovar::scale_fill_innova("npr") +
  theme_bw() +
  theme(
    plot.title = element_text(
      size = 14,
      hjust = 0.5,
      face = "bold"
    ),
    strip.text = element_text(
      size = 12,
      face = "bold"
    ),
    legend.title = element_text(
      size = 15,
      face = "bold"
    ),
    legend.text = ggtext::element_markdown(
      size = 14
    ),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_text(size = 12),
    axis.title = element_text(size = 13)
  )

sero_4malaria_district_distance_rh_hf_sampling
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/Fig 3. Seropositivity according to 4 malaria by district and distance category.png",
  sero_4malaria_district_distance_rh_hf_sampling,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)

ggsave(
  "./02_output/plots/Fig 3. Seropositivity according to 4 malaria by district and distance category.eps",
  sero_4malaria_district_distance_rh_hf_sampling,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)

ggsave(
  "./02_output/plots/Fig 3. Seropositivity according to 4 malaria by district and distance category.tiff",
  sero_4malaria_district_distance_rh_hf_sampling,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

## Segmented regression

```{r}
summarise_exposure_sero <- ffi_total_simply %>%
  drop_na(pf_exposure, pv_exposure, age_cat) %>%
  pivot_longer(
    cols = pv_exposure:pf_exposure,
    names_to = "exposure",
    values_to = "result_exposure"
  ) %>%
  mutate(
    exposure = case_when(
      exposure == "pv_exposure" ~ "P. Vivax Exposure",
      exposure == "pf_exposure" ~ "P. Falciparum Exposure"
    ),
    result_exposure = case_when(
      result_exposure == "Positive" ~ 1,
      TRUE ~ 0
    ),
    across(
      c(ffi_is_community:ffi_is_health_facility_name),
      str_to_title
    )
  ) %>%
  group_by(age_cat, exposure) %>%
  summarise(result_exposure = mean(result_exposure)) %>%
  ungroup()
```

```{r}
summarise_exposure_vivax <- summarise_exposure_sero %>% 
  filter(exposure == "P. Vivax Exposure") %>% 
  mutate(age = 1:8)

fit_vivax <- lm(result_exposure ~ age,
                data = summarise_exposure_vivax)


fit_vivax_segmented <- segmented::segmented(fit_vivax, ~ age)
summary(fit_vivax_segmented)

segmented::slope(fit_vivax_segmented)

plot(summarise_exposure_vivax$age, summarise_exposure_vivax$result_exposure)
plot(fit_vivax_segmented, add = TRUE)

summarise_exposure_falc <- summarise_exposure_sero %>% 
  filter(exposure == "P. Falciparum Exposure") %>% 
  mutate(age = 1:8)

fit_falc <- lm(result_exposure ~ age,
                data = summarise_exposure_falc)


fit_falc_segmented <- segmented::segmented(fit_falc, ~ age) 
summary(fit_falc_segmented)

segmented::slope(fit_falc_segmented)

plot(summarise_exposure_vivax$age, summarise_exposure_vivax$result_exposure)
plot(fit_vivax_segmented, add = TRUE)
```

## Figure 4 

### Belen

```{r}
malaria_seropositivity_type_exposure_belen <- ffi_total_simply %>%
  drop_na(pf_recent:pv_historic) %>%
  pivot_longer(
    cols = pf_recent:pv_historic,
    names_to = "malaria",
    values_to = "result_malaria"
  ) %>%
  filter(ffi_is_district == "Belen") %>%
  mutate(
    malaria = case_when(
      malaria == "pf_recent" ~ "Recent *P. falciparum*",
      malaria == "pf_historic" ~ "Historical *P. falciparum*",
      malaria == "pv_recent" ~ "Recent *P. vivax*",
      malaria == "pv_historic" ~ "Historical *P. vivax*"
    ),
    ffi_is_community = str_to_title(ffi_is_community),
    ffi_is_community = fct_reorder(
      ffi_is_community,
      distance_to_rh_hf_sampl_num, # distance_to_rh_hf_sampl_num
      .desc = TRUE
    )
  ) %>%
  ggplot(aes(
    x = ffi_is_community,
    fill = result_malaria
  )) +
  facet_wrap(vars(malaria)) +
  geom_bar(position = "fill", color = "black") +
  coord_flip(ylim = c(0, 0.25)) +
  scale_y_continuous(
    labels = scales::percent_format()
  ) +
  innovar::scale_fill_innova("npr") +
  labs(
    #title = str_wrap("Serological results by plasmodium type of malaria in the Belen District", 100),
    x = "Villages",
    y = "Percentage"
  ) +
  guides(
    fill = guide_legend("Results")
  ) +  
  theme_minimal() +
  theme(
    axis.title = element_text(
      face = "bold",
      size = 13
    ),
    axis.text.y = element_text(
      size = 12
    ),
    axis.text.x = element_text(
      size = 11
    ),
    legend.title = element_text(
      face = "bold",
      size = 14
    ),
    legend.text = element_text(size = 13),
    strip.text = ggtext::element_markdown(
      face = "bold",
      size = 12
    ),
    plot.tag.position = c(0, 1),
    plot.tag = element_text(face = 'bold')
  )

malaria_seropositivity_type_exposure_belen
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/Fig 4. Malaria Seropositivity by type of exposure and plasmodium in Belen.png",
  malaria_seropositivity_type_exposure_belen,
  dpi = 300,
  bg = "white",
  width = 10,
  height = 9
)
```

Alternative color:

```{r}
malaria_seropositivity_type_exposure_belen2 <- 
  malaria_seropositivity_type_exposure_belen +
  innovar::scale_fill_innova("dark_green",
                             reverse = TRUE)
malaria_seropositivity_type_exposure_belen2
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/Fig 4. Malaria Seropositivity by type of exposure and plasmodium in Belen2.png",
  malaria_seropositivity_type_exposure_belen2,
  dpi = 300,
  bg = "white",
  width = 10,
  height = 9
)
```

### Indiana

```{r}
malaria_seropositivity_type_exposure_indiana <- ffi_total_simply %>%
  drop_na(pf_recent:pv_historic) %>%
  pivot_longer(
    cols = pf_recent:pv_historic,
    names_to = "malaria",
    values_to = "result_malaria"
  ) %>%
  mutate(
    malaria = case_when(
      malaria == "pf_recent" ~ "Recent *P. falciparum*",
      malaria == "pf_historic" ~ "Historical *P. falciparum*",
      malaria == "pv_recent" ~ "Recent *P. vivax*",
      malaria == "pv_historic" ~ "Historical *P. vivax*"
    ),
    ffi_is_community = str_to_title(ffi_is_community),
  ) %>%
  filter(ffi_is_district == "Indiana") %>%
  mutate(
    ffi_is_community = fct_reorder(
      ffi_is_community,
      distance_to_rh_hf_sampl_num,
      .desc = TRUE
    )
  ) %>%
  ggplot(aes(
    x = ffi_is_community,
    fill = result_malaria
  )) +
  facet_wrap(vars(malaria)) +
  geom_bar(position = "fill", color = "black") +
  coord_flip(ylim = c(0, 0.25)) +
  scale_y_continuous(
    labels = scales::percent_format()
  ) +
  #ggsci::scale_fill_lancet() +
  innovar::scale_fill_innova("npr") +
  labs(
    #title = str_wrap("Serological results by plasmodium type of malaria in the Indiana District", 100),
    x = "Villages",
    y = "Percentage"
  ) +
  guides(
    fill = guide_legend("Results")
  ) +  
  theme_minimal() +
  theme(
    strip.text = ggtext::element_markdown(
      face = "bold",
      size = 12
    ),
    axis.title = element_text(
      face = "bold",
      size = 13
    ),
    axis.text.y = element_text(
      size = 12
    ),
    axis.text.x = element_text(
      size = 11
    ),
    legend.title = element_text(
      face = "bold",
      size = 14
    ),
    legend.text = element_text(size = 13),
    plot.tag.position = c(0.015, 1),
    plot.margin = margin(0, 0, 0, 0.5, "cm"),
    plot.tag = element_text(face = 'bold')
  )

malaria_seropositivity_type_exposure_indiana
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/Fig 4. Malaria Seropositivity by type of exposure and plasmodium in Indiana.png",
  malaria_seropositivity_type_exposure_indiana,
  dpi = 300,
  bg = "white",
  width = 10,
  height = 9
)
```

Alternative color:

```{r}
malaria_seropositivity_type_exposure_indiana2 <- 
  malaria_seropositivity_type_exposure_indiana +
  innovar::scale_fill_innova("dark_green",
                             reverse = TRUE)
malaria_seropositivity_type_exposure_indiana2
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/Fig 4. Malaria Seropositivity by type of exposure and plasmodium in Indiana2.png",
  malaria_seropositivity_type_exposure_indiana2,
  dpi = 300,
  bg = "white",
  width = 10,
  height = 9
)
```

### Combined plots

```{r}
malaria_seropositivity_type_exposure_belen <- 
  malaria_seropositivity_type_exposure_belen + 
  labs(
    #title = "Belen",
    y = NULL
  )

malaria_seropositivity_type_exposure_indiana <- 
  malaria_seropositivity_type_exposure_indiana + 
  labs(
    #title = "Indiana",
    y = NULL,
    x = NULL
  ) 
```

```{r}
malaria_seropositivity_type_exposure <- 
  malaria_seropositivity_type_exposure_belen + malaria_seropositivity_type_exposure_indiana +
  plot_layout(guides = "collect") +
  plot_annotation(tag_levels = "a") 
```

```{r eval = FALSE}
ggsave("02_output/plots/Fig 4. Malaria Seropositivity by type of exposure and plasmodium.png",
       malaria_seropositivity_type_exposure,
       dpi = 300,
       width = 15,
       height = 8,
       device = grDevices::png)

ggsave("02_output/plots/Fig 4. Malaria Seropositivity by type of exposure and plasmodium.eps",
       malaria_seropositivity_type_exposure,
       dpi = 300,
       width = 15,
       height = 8)

ggsave("02_output/plots/Fig 4. Malaria Seropositivity by type of exposure and plasmodium.tiff",
       malaria_seropositivity_type_exposure,
       dpi = 300,
       width = 15,
       height = 8)
```

# Supplementary Figure 4

```{r}
plot_ratio_lisa_communities_pf_belen <- ffi_w.getis_simply %>% 
  filter(ffi_is_district == "Belen",
         type_plasmodium %in% c("prevalence_pf_historic", 
                                "prevalence_pf_recent")) %>% 
  mutate(
    ffi_is_community = str_to_title(ffi_is_community),
    ffi_is_community = abbreviate(ffi_is_community, 15)
  ) %>% 
  ggplot(
    aes(x = ffi_is_community,
        fill = LISA_clust)
  ) + 
  geom_bar(position = "fill") +
  ggh4x::facet_grid2(ffi_is_district ~ type_plasmodium,
                     labeller = labeller(
                       type_plasmodium = c(
                         prevalence_pf_historic = "Historical Exposure",
                         prevalence_pf_recent = "Recent Exposure"
                       )
                     )) +
  innovar::scale_fill_innova("npr",
                             drop = FALSE) + 
  scale_y_continuous(labels = scales::percent) +
  labs(
    y = NULL,
    x = NULL,
    fill = "Local Indicator of \nSpatial Association"
  ) + 
  theme_minimal(base_size = 12) +
  theme(
    strip.text = ggtext::element_markdown(size = 14),
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    plot.tag = element_text(face = 'bold'),
    plot.tag.position = c(0, 0.98),
    plot.margin = margin(0, 0, 0, 0, "cm"),
    legend.title = element_text(
        face = "bold",
        size = 15
      ),
    legend.text = element_text(
      size = 14
    ),
    axis.text.y = element_text(size = 12),
    legend.margin = margin(l = -0.5, unit = "cm")
  )

plot_ratio_lisa_communities_pf_belen
```


```{r}
plot_ratio_lisa_communities_pv_belen <- ffi_w.getis_simply %>% 
  filter(ffi_is_district == "Belen",
         type_plasmodium %in% c("prevalence_pv_historic",
                                "prevalence_pv_recent")) %>% 
  mutate(
    ffi_is_community = str_to_title(ffi_is_community),
    ffi_is_community = abbreviate(ffi_is_community, 15)
  ) %>% 
  ggplot(
    aes(x = ffi_is_community,
        fill = LISA_clust)
  ) + 
  geom_bar(position = "fill") +
  ggh4x::facet_grid2(ffi_is_district ~ type_plasmodium,
                     labeller = labeller(
                       type_plasmodium = c(
                         prevalence_pv_historic = "Historical Exposure",
                         prevalence_pv_recent = "Recent Exposure"
                       )
                     )) +
  innovar::scale_fill_innova("npr",
                             drop = FALSE) + 
  scale_y_continuous(labels = scales::percent) +
  labs(
    y = NULL,
    x = NULL,
    fill = "Local Indicator of \nSpatial Association"
  ) + 
  theme_minimal(base_size = 12) +
  theme(
    axis.text.x = element_text(
      angle = 60, 
      hjust = 1, 
      #vjust = 0.5, 
      size = 10.5,
      lineheight = 0.85
    ),
    strip.text = ggtext::element_markdown(size = 14),
    plot.tag = element_text(face = 'bold'),
    plot.tag.position = c(0, 0.98),
    plot.margin = margin(0, 0, 0, 0, "cm"),
    legend.title = element_text(
        face = "bold",
        size = 15
      ),
    legend.text = element_text(
      size = 14
    ),
    axis.text.y = element_text(size = 12),
    legend.margin = margin(l = -0.5, unit = "cm")
  )

plot_ratio_lisa_communities_pv_belen
```

```{r}
plot_ratio_lisa_communities_pf_indiana <- ffi_w.getis_simply %>% 
  filter(ffi_is_district == "Indiana",
         type_plasmodium %in% c("prevalence_pf_historic", 
                                "prevalence_pf_recent")) %>% 
  mutate(
    ffi_is_community = str_to_title(ffi_is_community),
    ffi_is_community = abbreviate(ffi_is_community, 15)
  ) %>% 
  ggplot(
    aes(x = ffi_is_community,
        fill = LISA_clust)
  ) + 
  geom_bar(position = "fill") +
  ggh4x::facet_grid2(ffi_is_district ~ type_plasmodium,
                     labeller = labeller(
                       type_plasmodium = c(
                         prevalence_pf_historic = "Historical Exposure",
                         prevalence_pf_recent = "Recent Exposure"
                       )
                     )) +
  innovar::scale_fill_innova("npr",
                             drop = FALSE) + 
  scale_y_continuous(labels = scales::percent) +
  labs(
    y = NULL,
    x = NULL,
    fill = "Local Indicator of \nSpatial Association"
  ) + 
  theme_minimal(base_size = 12) +
  theme(
    strip.text = ggtext::element_markdown(size = 14),
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    plot.tag = element_text(face = 'bold'),
    plot.tag.position = c(0, 0.98),
    plot.margin = margin(0, 0, 0, 0, "cm"),
    legend.title = element_text(
        face = "bold",
        size = 15
      ),
    legend.text = element_text(
      size = 14
    ),
    axis.text.y = element_text(size = 12),
    legend.margin = margin(l = -0.5, unit = "cm")
  )

plot_ratio_lisa_communities_pf_indiana
```


```{r}
plot_ratio_lisa_communities_pv_indiana <- ffi_w.getis_simply %>% 
  filter(ffi_is_district == "Indiana",
         type_plasmodium %in% c("prevalence_pv_historic",
                                "prevalence_pv_recent")) %>% 
  mutate(
    ffi_is_community = str_to_title(ffi_is_community),
    ffi_is_community = abbreviate(ffi_is_community, 15)
  ) %>% 
  ggplot(
    aes(x = ffi_is_community,
        fill = LISA_clust)
  ) + 
  geom_bar(position = "fill") +
  ggh4x::facet_grid2(ffi_is_district ~ type_plasmodium,
                     labeller = labeller(
                       type_plasmodium = c(
                         prevalence_pv_historic = "Historical Exposure",
                         prevalence_pv_recent = "Recent Exposure" 
                       )
                     )) +
  innovar::scale_fill_innova("npr",
                             drop = FALSE) + 
  scale_y_continuous(labels = scales::percent) +
  labs(
    y = NULL,
    x = NULL,
    fill = "Local Indicator of \nSpatial Association"
  ) + 
  theme_minimal(base_size = 12) +
  theme(
    axis.text.x = element_text(
      angle = 60, 
      hjust = 1, 
      #vjust = 0.5, 
      size = 10.5,
      lineheight = 0.85
    ),
    strip.text = ggtext::element_markdown(size = 14),
    plot.tag = element_text(face = 'bold'),
    plot.tag.position = c(0, 0.98),
    plot.margin = margin(0, 0, 0, 0, "cm"),
    legend.title = element_text(
        face = "bold",
        size = 15
      ),
    legend.text = element_text(
      size = 14
    ),
    axis.text.y = element_text(size = 12),
    legend.margin = margin(l = -0.5, unit = "cm")
  )

plot_ratio_lisa_communities_pv_indiana
```

```{r}
combined_plot_lisa_communities <- plot_ratio_lisa_communities_pf_belen / plot_ratio_lisa_communities_pv_belen / plot_spacer() /
  plot_ratio_lisa_communities_pf_indiana / plot_ratio_lisa_communities_pv_indiana +
  plot_layout(guides = "collect", heights = c(2, 2, -0.1, 2, 2)) + # une las guías de los dos gráficos
  plot_annotation(tag_levels = "a")
```

```{r}
ggsave(
  "02_output/plots/Sup Fig 4. Local Indicator of Spatial Association by Communities.png",
  combined_plot_lisa_communities,
  width = 10,
  height = 12,
  dpi = 300,
  bg = "white"
)

ggsave(
  "02_output/plots/Sup Fig 4. Local Indicator of Spatial Association by Communities.tiff",
  combined_plot_lisa_communities,
  width = 10,
  height = 12,
  dpi = 300,
  bg = "white"
)
```


# Table descriptive

```{r}
library(gtsummary)

fisher.test.simulate.p.values <- function(data, variable, by, ...) {
  result <- list()
  test_results <- stats::fisher.test(data[[variable]], data[[by]], simulate.p.value = TRUE)
  result$p <- test_results$p.value
  result$test <- test_results$method
  result
}
```

## Table for district - Table 1

```{r}
labelled::var_label(ffi_total_simply) <- list(ffi_is_district = "Districs")
```


```{r}
table_1_ind <- ffi_total_simply %>%
  select(
    ffi_is_district,
    gender,
    education_level,
    economic_activities2,
    pf_exposure,
    pv_exposure,
    distance_to_rh_hf_sampl
  ) %>%
  #drop_na(pf_recent:pv_historic) %>% 
  tbl_summary(
    by = "ffi_is_district",
    missing = "no",
    # missing_text = "Missing",
    #type = all_dichotomous() ~ "categorical",
    digits = everything() ~ c(0, 1)
  ) %>%
  add_n() %>%
  add_overall(
    last = TRUE,
    digits = everything() ~ c(0, 1)
  ) %>%
  add_p(
     test = list(all_categorical() ~ "fisher.test.simulate.p.values"),
     pvalue_fun = scales::pvalue
  ) %>% 
  bold_p() %>%
  modify_header(label = "**Variable**") %>%
  bold_labels()
  
table_1_ind
```

```{r eval=FALSE}
table_1_ind %>%
  as_flex_table() %>%
  flextable::save_as_docx(path = "./02_output/reports/Table 1. Characteristics of the study population.docx")
```


## Table for Plasmodium and Exposure Time - Table 2

```{r}
table2_ind <- ffi_total_simply %>%
  select(
    ffi_is_district,
    gender,
    age_cat,
    economic_activities2,
    ffi_is_trip_month,
    ffi_is_place_shower,
    ffi_is_mosq_net,
    ffi_is_fever_month,
    ffi_is_symp_month_head,
    ffi_is_symp_month_musclejoint,
    ffi_is_symp_month_general,
    ffi_is_malaria,
    distance_to_rh_hf_sampl,
    pf_recent:pv_historic
  ) %>%
  drop_na(pf_recent:pv_historic) %>%
  pivot_longer(
    cols = pf_recent:pv_historic,
    names_to = "exposure",
    values_to = "result_exposure"
  ) %>%
  mutate(
    exposure = case_when(
      exposure == "pf_recent" ~ "P. Falciparum Recent",
      exposure == "pf_historic" ~ "P. Falciparum Historical",
      exposure == "pv_recent" ~ "P. Vivax Recent",
      exposure == "pv_historic" ~ "P. Vivax Historical"
      ),
    exposure = factor(exposure,
                      levels = c("P. Falciparum Recent",
                                 "P. Falciparum Historical",
                                 "P. Vivax Recent",
                                 "P. Vivax Historical"))
  ) %>%
  tbl_strata(
  strata = exposure,
  .tbl_fun =
    ~ .x %>%
        tbl_summary(
          by = result_exposure,
          missing = "no",
          #missing_text = "Missing",
          #type = all_dichotomous() ~ "categorical",
          digits = everything() ~ c(0, 1)
        ) %>%
        add_p(
          test = list(all_categorical() ~ "fisher.test.simulate.p.values"),
          pvalue_fun = scales::pvalue
        ) %>%
        bold_p() %>%
        modify_header(label = "**Variable**") %>%
        bold_labels() 
  ) 

table2_ind_overall <- ffi_total_simply %>%
  select(
    ffi_is_district,
    gender,
    age_cat,
    economic_activities2,
    ffi_is_trip_month,
    ffi_is_place_shower,
    ffi_is_mosq_net,
    ffi_is_fever_month,
    ffi_is_symp_month_head,
    ffi_is_symp_month_musclejoint,
    ffi_is_symp_month_general,
    ffi_is_malaria,
    distance_to_rh_hf_sampl
  ) %>%
  tbl_summary(
    missing = "no",
    # missing_text = "Missing",
    # type = all_dichotomous() ~ "categorical",
    digits = everything() ~ c(0, 1)
  ) %>%
  modify_header(
    label = "**Variable**"
  ) %>%
  bold_labels() %>%
  modify_spanning_header(stat_0 ~ "**Overall**")

table2_ind <- tbl_merge(
  tbls = list(table2_ind, table2_ind_overall),
  tab_spanner = FALSE
)

table2_ind
```

```{r eval=FALSE}
table2_ind %>%
  as_flex_table() %>%
  flextable::save_as_docx(path = "./02_output/reports/Characteristics of study population by P.vivax and falciparum exposure.docx")
```


## Table for Plasmodium and Exposure Time - Supplementary Table 2

```{r}
ffi_total_simply <- ffi_total_simply %>%
  mutate(
    age_cat2 = fct_collapse(
      age_cat,
      `[0-20)` = c("[0-10)", "[10-20)"),
      `[20+)` = c("[20-30)", "[30-40)", "[40-50)", "[50-60)", "[60-70)", "[70+)")
    )
  )


table3_ind_1 <- ffi_total_simply %>% 
  filter(age_cat2 == "[0-20)") %>% 
  select(
    ffi_is_district,
    gender,
    education_level,
    economic_activities2,
    pf_recent:pv_historic
  ) %>%
  drop_na(pf_recent:pv_historic) %>%
  pivot_longer(
    cols = pf_recent:pv_historic,
    names_to = "exposure",
    values_to = "result_exposure"
  ) %>%
  mutate(
    exposure = case_when(
      exposure == "pf_recent" ~ "P. Falciparum Recent",
      exposure == "pf_historic" ~ "P. Falciparum Historical",
      exposure == "pv_recent" ~ "P. Vivax Recent",
      exposure == "pv_historic" ~ "P. Vivax Historical"
      ),
    exposure = factor(exposure,
                      levels = c("P. Vivax Recent",
                                 "P. Vivax Historical",
                                 "P. Falciparum Recent",
                                 "P. Falciparum Historical"))
  ) %>%
  tbl_strata(
  strata = exposure,
  .tbl_fun =
    ~ .x %>%
        tbl_summary(
          by = result_exposure,
          missing = "no",
          #missing_text = "Missing",
          #type = all_dichotomous() ~ "categorical",
          digits = everything() ~ c(0, 1)
        ) %>%
        add_p(
          test = list(all_categorical() ~ "fisher.test.simulate.p.values"),
          pvalue_fun = scales::pvalue
        ) %>%
        bold_p() %>%
        modify_header(label = "**Variable**") %>%
        bold_labels() 
  ) 

table3_ind_1_overall <- ffi_total_simply %>%
  filter(age_cat2 == "[0-20)") %>% 
  select(
    ffi_is_district,
    gender,
    education_level,
    economic_activities2
  ) %>%
  tbl_summary(
    missing = "no",
    # missing_text = "Missing",
    # type = all_dichotomous() ~ "categorical",
    digits = everything() ~ c(0, 1)
  ) %>%
  modify_header(
    label = "**Variable**"
  ) %>%
  bold_labels() %>%
  modify_spanning_header(stat_0 ~ "**Overall**")

table3_ind_1 <- tbl_merge(
  tbls = list(table3_ind_1, table3_ind_1_overall),
  tab_spanner = FALSE
)

table3_ind_1
```

```{r eval=FALSE}
table3_ind_1 %>%
  as_flex_table() %>%
  flextable::save_as_docx(path = "./02_output/reports/sup_table_3_1.docx")
```

```{r}
table3_ind_2 <- ffi_total_simply %>% 
  filter(age_cat2 == "[20+)") %>% 
  select(
    ffi_is_district,
    gender,
    education_level,
    economic_activities2,
    pf_recent:pv_historic
  ) %>%
  drop_na(pf_recent:pv_historic) %>%
  pivot_longer(
    cols = pf_recent:pv_historic,
    names_to = "exposure",
    values_to = "result_exposure"
  ) %>%
  mutate(
    exposure = case_when(
      exposure == "pf_recent" ~ "P. Falciparum Recent",
      exposure == "pf_historic" ~ "P. Falciparum Historical",
      exposure == "pv_recent" ~ "P. Vivax Recent",
      exposure == "pv_historic" ~ "P. Vivax Historical"
      ),
    exposure = factor(exposure,
                      levels = c("P. Vivax Recent",
                                 "P. Vivax Historical",
                                 "P. Falciparum Recent",
                                 "P. Falciparum Historical"))
  ) %>%
  tbl_strata(
  strata = exposure,
  .tbl_fun =
    ~ .x %>%
        tbl_summary(
          by = result_exposure,
          missing = "no",
          #missing_text = "Missing",
          #type = all_dichotomous() ~ "categorical",
          digits = everything() ~ c(0, 1)
        ) %>%
        add_p(
          test = list(all_categorical() ~ "fisher.test.simulate.p.values"),
          pvalue_fun = scales::pvalue
        ) %>%
        bold_p() %>%
        modify_header(label = "**Variable**") %>%
        bold_labels() 
  ) 

table3_ind_2_overall <- ffi_total_simply %>%
  filter(age_cat2 == "[20+)") %>% 
  select(
    ffi_is_district,
    gender,
    education_level,
    economic_activities2
  ) %>%
  tbl_summary(
    missing = "no",
    # missing_text = "Missing",
    # type = all_dichotomous() ~ "categorical",
    digits = everything() ~ c(0, 1)
  ) %>%
  modify_header(
    label = "**Variable**"
  ) %>%
  bold_labels() %>%
  modify_spanning_header(stat_0 ~ "**Overall**")

table3_ind_2 <- tbl_merge(
  tbls = list(table3_ind_2, table3_ind_2_overall),
  tab_spanner = FALSE
)

table3_ind_2
```

```{r eval=FALSE}
table3_ind_2 %>%
  as_flex_table() %>%
  flextable::save_as_docx(path = "./02_output/reports/sup_table_3_2.docx")
```





```{r}
drugs_behavior <- ffi_total %>% 
  pivot_longer(
    cols = pf_recent:pv_historic,
    names_to = "malaria",
    values_to = "seropositivity"
  ) %>% 
  select(malaria, seropositivity, ffi_is_access_malaria_yn_hf, ffi_is_antimal_drug_use) %>% 
  drop_na() %>% 
  mutate(
    ffi_is_access_malaria_yn_hf = factor(ffi_is_access_malaria_yn_hf,
                                         labels = c("No", "Yes"))
  )

drugs_behavior %>% 
  group_by(malaria) %>% 
  summarise(
    broom::tidy(
      chisq.test(seropositivity, ffi_is_access_malaria_yn_hf)
    )
  )
drugs_behavior %>% 
  group_by(malaria) %>% 
  summarise(
    broom::tidy(
      chisq.test(seropositivity, ffi_is_antimal_drug_use)
    )
  )

ffi_total %>% 
  select(pf_historic, ffi_is_access_malaria_yn_hf) %>% 
  drop_na() %>% 
  table()

chisq.test(drugs_behavior$ffi_is_access_malaria_yn_hf, drugs_behavior$ffi_is_antimal_drug_use)

# Crear la tabla de contingencia
tab <- table(drugs_behavior$ffi_is_access_malaria_yn_hf, drugs_behavior$ffi_is_antimal_drug_use)
# Calcular el tamaño total de la muestra
n <- sum(tab)

# Calcular el valor de Chi-cuadrada
X_squared <- sum(chisq.test(tab)$statistic)

# Calcular el coeficiente de contingencia
C <- sqrt(X_squared / (n + X_squared))
C

library(epitools)
OR_result <- oddsratio(tab)
OR_result
```