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)
```

# Import data

```{r}
serology <- read_csv("./01_data/raw/FFI.peru.serology.kmeans2.4000.csv") %>%
  rename(ffi_is_code = Codigo) %>%
  select(-1)

ffi_individual <- read_csv("./01_data/raw/individuals_result_share_072022.csv")
ffi_household <- read_csv("./01_data/raw/household_gps_share_072022.csv")
data_2000_2019 <- read_csv("./01_data/raw/Data_FFI_2000_2019_20212405.csv")
data_2020_2021 <- read_csv("./01_data/raw/Data_FFI_2020_2021_20210629.csv")
poblacion <- readRDS("./01_data/raw/poblacion_inei_2017.rds")
```

# Format data

## Individuals

```{r}
ffi_total <- ffi_individual %>%
  full_join(
    serology %>%
      mutate(
        ffi_is_code = paste0(0, ffi_is_code)
      ),
    by = "ffi_is_code"
  )

ffi_is_symp_top3 <- ffi_total %>%
  pivot_longer(
    cols = ffi_is_symp_month_head:ffi_is_symp_month_nosymptopm,
    names_to = "ffi_is_symp",
    values_to = "ffi_is_symp_ans"
  ) %>% 
  count(ffi_is_symp, ffi_is_symp_ans) %>% 
  filter(ffi_is_symp_ans != 0,
         ffi_is_symp != "ffi_is_symp_month_nosymptopm") %>%
  slice_max(n, n = 3)

ffi_total <- ffi_total %>%
  mutate(
    ffi_is_district = fct_recode(
      ffi_is_district,
      "Indiana" = "INDIANA",
      "Belen" = "BELEN"
    ),
    age_cat = case_when(
      ffi_is_age_fixed >= 70 ~ "[70+)",
      ffi_is_age_fixed >= 60 ~ "[60-70)",
      ffi_is_age_fixed >= 50 ~ "[50-60)",
      ffi_is_age_fixed >= 40 ~ "[40-50)",
      ffi_is_age_fixed >= 30 ~ "[30-40)",
      ffi_is_age_fixed >= 20 ~ "[20-30)",
      ffi_is_age_fixed >= 10 ~ "[10-20)",
      ffi_is_age_fixed >= 0 ~ "[0-10)"
    ),
    age_cat = factor(age_cat),
    age_code = case_when(
      age_cat == "[70+)" ~ 8,
      age_cat == "[60-70)" ~ 7,
      age_cat == "[50-60)" ~ 6,
      age_cat == "[40-50)" ~ 5,
      age_cat == "[30-40)" ~ 4,
      age_cat == "[20-30)" ~ 3,
      age_cat == "[10-20)" ~ 2,
      age_cat == "[0-10)" ~ 1,
    ),
    gender = factor(
      ffi_is_sex,
      labels = c("Male", "Female")
    ),
    ffi_is_malaria = factor(
      ffi_is_malaria,
      labels = c("No", "Yes", "Don't Know / No Answer")
    ),
    across(
      c(pf_recent:pv_historic),
      ~ factor(., labels = c("Negative", "Positive"))
    ),
    across(
      c(
        ffi_is_fever_month,
        ffi_is_antimal_drug_use, 
        ffi_is_mosq_net
      ),
      ~ factor(., labels = c("No", "Yes"))
    ),
    ffi_is_trip_month = factor(
      ffi_is_trip_month,
      labels = c("No", "Yes", "Don't Know/No Answer")
    ),
    ffi_is_mal_lifetime = factor(
      ffi_is_mal_lifetime,
      labels = c(
        "1 to 3 times",
        "3 to 7 times",
        "More than 7 times"
      )
    ),
    ffi_is_place_shower = factor(
      ffi_is_place_shower,
      labels = c(
        "Bathroom inside the dwelling",
        "Bathroom outside the dwelling",
        "In the countryside/river",
        "Other"
      )
    ),
    education_level = case_when(
      ffi_is_inst_level %in% 0 ~ "No schooling", 
      ffi_is_inst_level %in% 1:2 ~ "Primary school",
      ffi_is_inst_level %in% 3:4 ~ "Secondary school",
      ffi_is_inst_level %in% 5:6 ~ "Higher education"
    ),
    education_level = fct_relevel(education_level, "No schooling",
                                  "Primary school",
                                  "Secondary school"),
    economic_activities = factor(ffi_is_main_econ_act,
                                 labels = c("None",
                                            "Day labourer",
                                            "Wood extractor",
                                            "Fisherman",
                                            "Livestock farmer",
                                            "Farmer",
                                            "Trader",
                                            "Housewife",
                                            "Student",
                                            "Motorcycle taxi driver",
                                            "Other")),
    economic_activities2 = fct_collapse(
      economic_activities,
      "Forest Related" = c("Day labourer",
                           "Wood extractor",
                           "Fisherman",
                           "Livestock farmer",
                           "Farmer"),
      "Housewife/student" = c("Housewife",
                              "Student"),
      "Others" = c("Motorcycle taxi driver",
                   "Other")
    ),
    ffi_is_place_sleep = factor(ffi_is_place_sleep,
                                labels = c("Inside the dwelling",
                                           "Outside the dwelling",
                                           "Don't Know / Don't Answer")),
    ffi_is_mosq_net_freq = factor(ffi_is_mosq_net_freq,
                                  levels = c(0:2, 99),
                                  labels = c("Never",
                                             "Every day",
                                             "Sometimes",
                                             "Don't Know / Don't Answer")),
    across(
      c(ffi_is_symp_month_head,
        ffi_is_symp_month_musclejoint,
        ffi_is_symp_month_general),
      ~ factor(
        .,
        labels = c("No", "Yes")
      )
    ),
    ffi_is_knowl_mal_inf = factor(ffi_is_knowl_mal_inf,
                                  labels = c("I was diagnosed at the health facility",
                                             "Because of symptoms",
                                             "Other")),
    pv_exposure = case_when(
      pv_recent == "Negative" & pv_historic == "Negative" ~ "Negative",
      is.na(pv_recent) & is.na(pv_historic) ~ NA_character_,
      TRUE ~ "Positive"
    ),
    pf_exposure = case_when(
      pf_recent == "Negative" & pf_historic == "Negative" ~ "Negative",
      is.na(pf_recent) & is.na(pf_historic) ~ NA_character_,
      TRUE ~ "Positive"
    ),
    recent_exposure = case_when(
      pv_recent == "Negative" & pf_recent == "Negative" ~ "Negative",
      is.na(pv_recent) & is.na(pf_recent) ~ NA_character_,
      TRUE ~ "Positive"
    ),
    historical_exposure = case_when(
      pv_historic == "Negative" & pf_historic == "Negative" ~ "Negative",
      is.na(pv_historic) & is.na(pf_historic) ~ NA_character_,
      TRUE ~ "Positive"
    ),
    any_pv_pf_exposure = case_when(
      recent_exposure == "Negative" & historical_exposure == "Negative" ~ "Negative",
      is.na(recent_exposure) & is.na(historical_exposure) ~ NA_character_,
      TRUE ~ "Positive"
    ),
    only_pv_exposure = case_when(
      pv_exposure == "Positive" & pf_exposure == "Negative" ~ "Positive", 
      is.na(pv_exposure) & is.na(pf_exposure) ~ NA_character_,
      TRUE ~ "Negative"
    ),
    only_pf_exposure = case_when(
      pf_exposure == "Positive" & pv_exposure == "Negative" ~ "Positive", 
      is.na(pf_exposure) & is.na(pv_exposure) ~ NA_character_,
      TRUE ~ "Negative"
    ),
    pv_pf_exposure = case_when(
      pv_exposure == "Positive" & pf_exposure == "Positive" ~ "Positive",
      is.na(pv_exposure) & is.na(pf_exposure) ~ NA_character_,
      TRUE ~ "Negative"
    ),
    freedom_malaria = case_when(
      pv_exposure == "Negative" & pf_exposure == "Negative" ~ "Positive", 
      is.na(pv_exposure) & is.na(pf_exposure) ~ NA_character_,
      TRUE ~ "Negative"
    ),
    only_pv_recent = case_when(
      pv_recent == "Positive" & pf_recent == "Negative" ~ "Positive", 
      is.na(pv_recent) & is.na(pf_recent) ~ NA_character_,
      TRUE ~ "Negative"
    ),
    only_pf_recent = case_when(
      pf_recent == "Positive" & pv_recent == "Negative" ~ "Positive", 
      is.na(pf_recent) & is.na(pv_recent) ~ NA_character_,
      TRUE ~ "Negative"
    ),
    pv_pf_recent = case_when(
      pv_recent == "Positive" & pf_recent == "Positive" ~ "Positive",
      is.na(pv_recent) & is.na(pf_recent) ~ NA_character_,
      TRUE ~ "Negative"
    ),
    freedom_malaria_recent = case_when(
      pv_recent == "Negative" & pf_recent == "Negative" ~ "Positive", 
      is.na(pv_recent) & is.na(pf_recent) ~ NA_character_,
      TRUE ~ "Negative"
    ),
    only_pv_historic = case_when(
      pv_historic == "Positive" & pf_historic == "Negative" ~ "Positive", 
      is.na(pv_historic) & is.na(pf_historic) ~ NA_character_,
      TRUE ~ "Negative"
    ),
    only_pf_historic = case_when(
      pf_historic == "Positive" & pv_historic == "Negative" ~ "Positive", 
      is.na(pf_historic) & is.na(pv_historic) ~ NA_character_,
      TRUE ~ "Negative"
    ),
    pv_pf_historic = case_when(
      pv_historic == "Positive" & pf_historic == "Positive" ~ "Positive",
      is.na(pv_historic) & is.na(pf_historic) ~ NA_character_,
      TRUE ~ "Negative"
    ),
    freedom_malaria_historic = case_when(
      pv_historic == "Negative" & pf_historic == "Negative" ~ "Positive", 
      is.na(pv_historic) & is.na(pf_historic) ~ NA_character_,
      TRUE ~ "Negative"
    ),
    across(c(pv_exposure:freedom_malaria_historic), factor)
  )


labelled::var_label(ffi_total) <- list(
  ffi_is_district = "Districs",
  gender = "Gender",
  ffi_is_age = "Age",
  education_level = "Education Level",
  economic_activities = "Economic Activities",
  economic_activities2 = "Economic Activities",
  ffi_is_trip_month = "Trip in the last month",
  ffi_is_place_shower	= "Usual places to bathe",
  ffi_is_place_sleep = "Place to sleep last night",
  ffi_is_mosq_net = "Mosquito Net",
  ffi_is_mosq_net_freq = "Frequency of Mosquito Net",
  ffi_is_malaria = "Ever had malaria",
  ffi_is_fever_month = "History of fever in the last month",
  ffi_is_knowl_mal_inf = "Knowledge of whether it was malaria",
  ffi_is_symp_month_head = "Symptoms or signs in the last month: Headache",
  ffi_is_symp_month_musclejoint = "Symptoms or signs in the last month: Muscle or joint pain",
  ffi_is_symp_month_general = "Symptoms or signs in the last month: General malaise",
  pf_recent = "Recent P. Falciparum",
  pf_historic = "Historical P. Falciparum",
  pv_recent = "Recent P. Vivax",
  pv_historic = "Historical P. Vivax",
  pv_exposure = "P. Vivax Exposure",
  pf_exposure = "P. Falciparum Exposure"
)
```

## Household

```{r}
ffi_household <- ffi_household %>% 
  mutate(
    ffi_h_district = fct_recode(
      ffi_h_district,
      "Indiana" = "INDIANA",
      "Belen" = "BELEN"
    ),
    ffi_h_type = factor(ffi_h_type,
                        labels = c("HOUSE TYPE A: House without walls, with or without wooden floor",
                                   "HOUSE TYPE B: House that has one, two or three walls without completely enclosing a room, with or without wooden floor",
                                   "HOUSE TYPE C: House that has at least one room enclosed, with or without wooden floor",
                                   "HOUSE TYPE D: Completely enclosed house, with or without wooden floor")),
    ffi_h_sprayed = factor(ffi_h_sprayed,
                           labels = c("No", "Yes",
                                      "Don't know / No answer")),
    ffi_h_watersource = factor(ffi_h_watersource,
                               levels = c(1:6, 88, 99),
                               labels = c("Drinking water inside the house",
                                          "Drinking water outside the house",
                                          "Well",
                                          "Public well",
                                          "River or stream",
                                          "Rainwater",
                                          "Other",
                                          "Don't know/No answer"))
  )

labelled::var_label(ffi_household) <- list(
  ffi_h_district = "District",
  ffi_h_mti = "Insecticide-treated mosquito nets",
  ffi_h_type = "Type of house",
  ffi_h_sprayed = "Househould sprayed with insecticides in the last 12 months by the ministry of health",
  ffi_h_watersource = "Main source of water for human consumption (food preparation)"
)
```


# Distance Group about communities (Method Sampling)

```{r}
communities_sf <- ffi_household %>%
  select(
    ffi_h_district,
    ffi_h_health_facility_name,
    ffi_h_code_community,
    ffi_h_community,
    ffi_gps_long,
    ffi_gps_lat
  ) %>%
  st_as_sf(
    coords = c("ffi_gps_long", "ffi_gps_lat"),
    crs = 4326
  ) %>%
  group_by(ffi_h_district, ffi_h_health_facility_name, ffi_h_code_community, ffi_h_community) %>%
  summarise() %>%
  st_centroid()

names_hf <- ffi_household %>% 
  count(ffi_h_health_facility_name) %>% 
  pull(ffi_h_health_facility_name)

hf_shp <- readRDS("01_data/raw/HF_shp.rds")

hf_shp <- hf_shp %>% 
  filter(`Health Facility Name` %in% names_hf) %>% 
  select(District, 
         ffi_h_health_facility_name = `Health Facility Name`)
```

## Distance communities and their HF's

```{r}
communities_distances_hf_sampling <- communities_sf %>% 
  ungroup() %>% 
  st_transform(crs = 32718) %>% 
  st_distance(
    hf_shp %>% 
      st_transform(crs = 32718)
  ) %>% 
  as_tibble() 

names(communities_distances_hf_sampling) <- hf_shp$ffi_h_health_facility_name

communities_distances_hf_sampling <- communities_distances_hf_sampling %>% 
  mutate(
    ffi_is_community = communities_sf$ffi_h_community,
    ffi_is_district = communities_sf$ffi_h_district,
    .before = 1
  ) %>% 
  pivot_longer(
    cols = -c(ffi_is_community, ffi_is_district),
    names_to = "ffi_is_health_facility_name",
    values_to = "distance_to_hf_sampl_num"
  ) %>% 
  mutate(distance_to_hf_sampl_num = as.numeric(distance_to_hf_sampl_num))

communities_distances_hf_sampling_cat <- communities_sf %>% 
  as_tibble() %>% 
  ungroup() %>% 
  select(
    ffi_is_district = ffi_h_district,
    ffi_is_health_facility_name = ffi_h_health_facility_name, 
    ffi_is_community = ffi_h_community
  ) %>% 
  left_join(
    communities_distances_hf_sampling
  )

communities_distances_hf_sampling_cat <- communities_distances_hf_sampling_cat %>% 
  group_by(ffi_is_health_facility_name) %>% 
  mutate(
    distance_to_hf_sampl = as.numeric(cut_number(distance_to_hf_sampl_num, 4)),
    distance_to_hf_sampl = factor(distance_to_hf_sampl, 
                                  labels = c("Proximate", "Moderate",
                                             "Distant", "Extra Distant"))
  ) %>% 
  ungroup()
```

Update info:

```{r}
ffi_total <- ffi_total %>% 
  left_join(
    communities_distances_hf_sampling_cat
  ) 
```



## Distance communities and Regional Hospital

```{r}
iquitos_centro <- tibble(
  "ffi_h_health_facility_name" = "Hospital Regional de Loreto",
  Longitude = -73.25385902080906,
  Latitude = -3.7264060164148716,
) %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326) %>%
  st_transform(crs = 32718)

communities_distances_rh_sampling <- communities_sf %>%
  st_transform(crs = 32718) %>%
  st_distance(iquitos_centro)

communities_distances_rh_sampling <- enframe(communities_distances_rh_sampling) %>%
  mutate(
    ffi_is_cod_com = communities_sf$ffi_h_code_community,
    ffi_is_community = communities_sf$ffi_h_community,
    distance_to_rh_sampl_num = as.numeric(value),
    ffi_is_district = communities_sf$ffi_h_district
  ) %>%
  select(-c(name, value)) %>% 
  group_by(ffi_is_district) %>% 
  mutate(
    distance_to_rh_sampl = as.numeric(cut_number(distance_to_rh_sampl_num, 4)),
    distance_to_rh_sampl = factor(distance_to_rh_sampl, 
                                  labels = c("Proximate", "Moderate",
                                             "Distant", "Extra Distant"))
  ) %>% 
  ungroup()
```

Las Comunidades mas cercanas y mas lejanas por distrito:

```{r eval=FALSE}
communities_distances_rh_sampling %>%
  group_by(ffi_is_district) %>%
  slice_min(distance_to_rh_sampl, n = 1) %>%
  bind_rows(
    communities_distances_rh_sampling %>%
      mutate(
        ffi_is_district = communities_sf$ffi_h_district
      ) %>%
      group_by(ffi_is_district) %>%
      slice_max(distance_to_rh_sampl, n = 1)
  ) %>%
  arrange(ffi_is_district, distance_to_rh_sampl)
```

Update info:

```{r}
ffi_total <- ffi_total %>% 
  left_join(
    communities_distances_rh_sampling 
  ) 
```

## Distance communities to Regional Hospital and Health Facilities (min distance)

```{r}
ffi_total <- ffi_total %>% 
  rowwise() %>% 
  mutate(
    distance_to_rh_hf_sampl_num = min(distance_to_rh_sampl_num, distance_to_hf_sampl_num)
  ) %>% 
  group_by(ffi_is_district) %>% 
  mutate(
    distance_to_rh_hf_sampl = as.numeric(cut_number(distance_to_rh_hf_sampl_num, 4)),
    distance_to_rh_hf_sampl = factor(distance_to_rh_hf_sampl, 
                                     labels = c("Proximate", "Moderate",
                                                "Distant", "Extra Distant"))
  ) %>% 
  ungroup()


labelled::var_label(ffi_total) <- list(
  distance_to_rh_hf_sampl = "Distance Category to Regional Hospital and Health Facilities"
)
```

# Distance Group about communities (Method Census)

```{r}
comunidades_info_total <- readRDS("01_data/raw/Información poblacional comunidades.RData") %>% 
  st_as_sf()

st_geometry(comunidades_info_total[comunidades_info_total$MNOMCP == "CABO PANTOJA", ]) <-
  st_sfc(st_point(c(-73.248636296, -3.777654358)))


comunidades_info_total <- comunidades_info_total %>% 
  mutate(
    geometry2 = case_when(
      MNOMCP == "8 DE DICIEMBRE" ~ st_cast(
        st_sfc(st_point(c(-73.248636296, 
                          -3.777654358))), 
        "GEOMETRY"
        ),
      MNOMCP == "CABO LOPEZ" ~ st_cast(
        st_sfc(st_point(c(-73.253457539, 
                          -3.781405681))), 
        "GEOMETRY"
        ),
      MNOMCP == "TIMICURO GRANDE" ~ st_cast(
        st_sfc(st_point(c(-73.023866667, 
                          -3.558536667))), 
        "GEOMETRY"
        ),
      TRUE ~ st_cast(geometry, "GEOMETRY")
      )
  )
  
comunidades_info_total <- comunidades_info_total %>% 
  mutate(
    geometry = ifelse(
      MNOMCP == "8 DE DICIEMBRE", st_sfc(st_point(c(-73.248636296, 
                                                    -3.777654358))),
      ifelse(MNOMCP == "CABO LOPEZ", st_sfc(st_point(c(-73.253457539, 
                                                       -3.781405681))),
             ifelse(MNOMCP == "TIMICURO GRANDE",
                    st_sfc(st_point(c(-73.023866667, 
                          -3.558536667))), geometry)
  )


comunidades_info_adic <- comunidades_info_total %>% 
  filter(st_is_empty(geometry)) %>% 
  distinct(MNOMCP, .keep_all = TRUE) %>% 
  filter(`Health Facility Name` %in% hf_shp$ffi_h_health_facility_name) %>% 
  select(-geometry) %>% 
  left_join(
    communities_sf %>% 
      select(
        District = ffi_h_district,
        `Health Facility Name` = ffi_h_health_facility_name,
        MNOMCP = ffi_h_community,
        geometry
      ) %>% 
      as_tibble()
  )

comunidades_info_total <- comunidades_info_total %>% 
  filter(!st_is_empty(geometry)) %>% 
  distinct(MNOMCP, .keep_all = TRUE) %>% 
  bind_rows(
    comunidades_info_adic
  )
```

## Distance communities and their HF's

```{r}
communities_distances_hf <- comunidades_info_total %>% 
  filter(`Health Facility Name` %in% hf_shp$ffi_h_health_facility_name) %>% 
  st_as_sf(crs = 4326) %>% 
  st_transform(crs = 32718) %>% 
  st_distance(
    hf_shp %>% 
      st_transform(crs = 32718)
  ) %>% 
  as_tibble() 

names_comunities_9_hf <- comunidades_info_total %>% 
  filter(`Health Facility Name` %in% hf_shp$ffi_h_health_facility_name) %>% 
  pull(MNOMCP)

names_districts_9_hf <- comunidades_info_total %>% 
  filter(`Health Facility Name` %in% hf_shp$ffi_h_health_facility_name) %>% 
  pull(District)

names(communities_distances_hf) <- hf_shp$ffi_h_health_facility_name

communities_distances_hf <- communities_distances_hf %>% 
  mutate(
    ffi_is_community = names_comunities_9_hf,
    ffi_is_district = names_districts_9_hf,
    .before = 1
  ) %>% 
  pivot_longer(
    cols = -c(ffi_is_community, ffi_is_district),
    names_to = "ffi_is_health_facility_name",
    values_to = "distance_to_hf_census"
  ) %>% 
  mutate(distance_to_hf_census = as.numeric(distance_to_hf_census))

communities_distances_hf_total <- comunidades_info_total %>% 
  filter(`Health Facility Name` %in% hf_shp$ffi_h_health_facility_name) %>% 
  select(
    ffi_is_district = District,
    ffi_is_health_facility_name = `Health Facility Name`, 
    ffi_is_community = MNOMCP
  ) %>% 
  left_join(
    communities_distances_hf
  )
  
  
communities_distances_hf_cat <- communities_distances_hf_total %>%
  filter(ffi_is_health_facility_name != "SOLEDAD DE VILLA BELEN") %>%
  group_by(ffi_is_health_facility_name) %>%
  mutate(
    distance_to_hf_census = as.numeric(cut_number(distance_to_hf_census, 4)),
    distance_to_hf_census = factor(
      distance_to_hf_census,
      labels = c(
        "Proximate",
        "Moderate",
        "Distant",
        "Extra Distant"
      )
    )
  ) %>%
  ungroup() %>%
  bind_rows(
    communities_distances_hf_total %>% 
      filter(ffi_is_health_facility_name == "SOLEDAD DE VILLA BELEN") %>% 
      mutate(
        distance_to_hf_census = "Moderate"
      )
  )
```

Update info:

```{r}
ffi_total <- ffi_total %>% 
  left_join(
    communities_distances_hf_cat
  )
```


## Distance communities and Regional Hospital

```{r}
communities_distances_rh <- comunidades_info_total %>% 
  st_as_sf(crs = 4326) %>% 
  st_transform(crs = 32718) %>% 
  st_distance(
    iquitos_centro
  ) %>% 
  enframe() %>% 
  mutate(
    ffi_is_community = comunidades_info_total$MNOMCP,
    distance_to_rh_census = as.numeric(value),
    ffi_is_district = comunidades_info_total$District
  ) %>%
  select(-c(name, value))
```

```{r}
communities_distances_rh_cat <- communities_distances_rh %>%
  drop_na(distance_to_rh_census) %>%
  arrange(ffi_is_district, distance_to_rh_census) %>%
  group_by(ffi_is_district) %>%
  mutate(
    distance_to_rh_census = as.numeric(cut_number(distance_to_rh_census, 4)),
    distance_to_rh_census = factor(
      distance_to_rh_census,
      labels = c(
        "Proximate",
        "Moderate",
        "Distant",
        "Extra Distant"
      )
    )
  ) %>%
  ungroup() %>%
  distinct(ffi_is_community, ffi_is_district, .keep_all = TRUE)
```

Update info:

```{r}
ffi_total <- ffi_total %>% 
  left_join(
    communities_distances_rh_cat
  )
```


```{r eval = FALSE}
saveRDS(ffi_total, 
        file = "01_data/processed/ffi_total.rds")
```





# 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

### Individuals

```{r}
table_1_ind <- ffi_total %>%
  select(
    ffi_is_district,
    gender,
    ffi_is_age,
    education_level,
    economic_activities2,
    ffi_is_trip_month,
    ffi_is_place_shower,
    ffi_is_place_sleep,
    ffi_is_mosq_net,
    ffi_is_mosq_net_freq,
    ffi_is_fever_month,
    ffi_is_malaria,
    ffi_is_symp_month_head,
    ffi_is_symp_month_musclejoint,
    ffi_is_symp_month_general,
    pf_exposure,
    pv_exposure,
    distance_to_rh_hf_sampl
    # ffi_is_knowl_mal_inf,
    # pf_recent:pv_historic
  ) %>%
  #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/table1_ind.docx")
```

### Households


```{r}
table_1_household <- ffi_household %>%
  select(
    ffi_h_district,
    ffi_h_mti,
    ffi_h_type,
    ffi_h_sprayed,
    ffi_h_watersource
  ) %>%
  tbl_summary(
    by = "ffi_h_district",
    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_household
```


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

## Table for Plasmodium and Exposure Time

```{r}
table2_ind <- ffi_total %>%
  select(
    ffi_is_district,
    gender,
    ffi_is_age,
    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 %>%
  select(
    ffi_is_district,
    gender,
    ffi_is_age,
    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/table2_ind.docx")
```




# Plots

## Map

```{r}
library(ggmap)
```

```{r}
count_communities <- ffi_total %>% 
  count(ffi_is_community, ffi_is_cod_com, name = "sampling") %>% 
  mutate(
    sampling_size = case_when(
      sampling >= 180 ~ "180+",
      sampling >= 120 ~ "120-179",
      sampling >= 70 ~ "70-119",
      TRUE ~ "< 70"
    ),
    sampling_size = fct_relevel(sampling_size,
                                "120-179",
                                "180+", after = Inf)
  ) 

communities_sf2 <- communities_sf %>% 
  left_join(count_communities,
            by = c("ffi_h_code_community" = "ffi_is_cod_com",
                   "ffi_h_community" = "ffi_is_community"))


# od_hf_communities <- communities_sf %>% 
#   as_tibble() %>% 
#   select(O = ffi_h_health_facility_name,
#          D = ffi_h_community)
# 
# places_togueter <- hf_shp %>% 
#   rename(Place = ffi_h_health_facility_name) %>% 
#   select(-District) %>% 
#   mutate(Type = "Health Facility") %>% 
#   bind_rows(
#     communities_sf %>% 
#       select(Place = ffi_h_community) %>% 
#       mutate(Type = "Community")
#   )
# 
# places_sf_linestring <- od::od_to_sf(od_hf_communities,
#                                      places_togueter)
```


```{r}
ph_basemap <- readRDS("01_data/raw/gmap_shp.rds")
data(Peru, package = "innovar")
distr_shp <- Peru %>% 
  filter(prov == "MAYNAS",
         distr %in% c("INDIANA", "BELEN")) 

map_by_hf <- ggmap(ph_basemap) +
  coord_sf(crs = st_crs(4326)) +
  geom_sf(data = distr_shp %>% 
            mutate(distr = str_to_title(distr)),
          #aes(fill = distr),
          size = 1,
          alpha = 0.5,
          color = "black",
          inherit.aes = FALSE) +
  # scale_fill_manual(
  #   values = c("#64b2af", "#feee67")
  # ) + 
  geom_sf(data = communities_sf2 %>% 
            mutate(Type = "Community (38)",
                   ffi_h_health_facility_name = str_to_title(ffi_h_health_facility_name)),
          aes(size = sampling_size,
              shape = Type,
              #color = Type,
              color = ffi_h_health_facility_name),
          alpha = 0.8,
          inherit.aes = FALSE) +
  geom_sf(data = hf_shp %>% 
            mutate(Type = "Health Facility (9)",
                   ffi_h_health_facility_name = str_to_title(ffi_h_health_facility_name)),
          aes(shape = Type,
              #color = Type,
              color = ffi_h_health_facility_name),
          size = 5,
          inherit.aes = FALSE
          ) +
  # ggsflabel::geom_sf_text_repel(data = hf_shp,
  #                               aes(label = ffi_h_health_facility_name),
  #                               inherit.aes = FALSE) +
  # geom_sf(data = places_sf_linestring,
  #         linetype = "dashed",
  #         inherit.aes = FALSE) +
  scale_shape_manual(
    values = c(16, 18)
  ) +
  innovar::scale_color_innova("npr") +
  # scale_color_manual(
  #   values = c("black", "red")
  # ) + 
  # ggsflabel::geom_sf_text_repel(data = hf_shp,
  #                               aes(label = label2),
  #                               force = 10,
  #                               force_pull = -30,
  #                               direction = "y",
  #                               nudge_x = 0.35,
  #                               seed = 2021,
  #                               inherit.aes = FALSE) +
  coord_sf(xlim = c(-73.42, -72.56),
           ylim = c(-4.42, -3.35), 
           expand = FALSE) +
  labs(x = NULL,
       y = NULL,
       title = "Communities and Health Facilities for the Study") +
  guides(
    size = guide_legend("Sampling size"),
    fill = guide_legend("District"),
    color = guide_legend("Health Facilities")
  ) +
  theme_bw(base_size = 12) + 
  theme(
    plot.title = element_text(size = 14,
                              hjust = 0.5,
                              face = "bold")
  )
map_by_hf
```

```{r}
ggsave("02_output/plots/map_communities_hf_.png",
       map_by_hf,
       dpi = 400,
       device = grDevices::png,
       height = 5,
       width = 7,
       scale = 1.2)
```


## Malaria Annual Parasite Index


```{r}
data_malaria <- bind_rows(
  data_2000_2019,
  data_2020_2021
) %>%
  drop_na(District)

pob_loreto <- poblacion %>%
  filter(dep == "LORETO") %>%
  select(District = distr, Total)
```


```{r}
malaria_cases_pob <- data_malaria %>%
  rowwise() %>%
  mutate(
    cases_malaria = sum(
      c_across(c(
        `Confirmed P. Falciparum`,
        `Confirmed P. Vivax`
      )),
      na.rm = TRUE
    )
  ) %>%
  group_by(District) %>%
  summarise(cases_malaria = sum(cases_malaria))

malaria_cases_pob <- malaria_cases_pob %>%
  left_join(
    pob_loreto
  ) %>%
  mutate(
    api = cases_malaria * 1000 / Total
  )

data(Peru, package = "innovar")

malaria_cases_pob_sf <- Peru %>%
  filter(dep == "LORETO") %>%
  select(District = distr, geometry) %>% 
  right_join(
    malaria_cases_pob
  )
```


```{r}
figure1 <- ggplot(malaria_cases_pob_sf) +
  geom_sf(aes(fill = api, geometry = geometry)) +
  geom_sf(data = malaria_cases_pob_sf %>% dplyr::filter(District == "BELEN"), colour = "#aa6439", fill = NA, size = 1.5, aes(fill = api, geometry = geometry)) +
  geom_sf(data = malaria_cases_pob_sf %>% dplyr::filter(District == "INDIANA"), colour = "#256e5d", fill = NA, size = 1.5, aes(fill = api, geometry = geometry)) +
  scale_fill_distiller(
    name = "API",
    na.value = "black",
    # limits = c(0, 4000),
    palette = "RdYlGn",
    direction = -1,
    breaks = scales::pretty_breaks(n = 5),
    values = c(
      0,
      0.001,
      0.002,
      0.005,
      0.008,
      0.01,
      0.02,
      0.03,
      0.05,
      0.1,
      0.6,
      0.9,
      1
    )
  ) +
  labs(
    x = NULL,
    y = NULL,
    title = NULL
  ) +
  theme_classic() +
  geom_sf_label_repel(
    data = malaria_cases_pob_sf %>% dplyr::filter(District == "BELEN"),
    aes(label = District),
    min.segment.length = 0,
    force = 100,
    nudge_x = -1,
    seed = 10,
    colour = "#aa6439"
  ) +
  geom_sf_label_repel(
    data = malaria_cases_pob_sf %>% dplyr::filter(District == "INDIANA"),
    aes(label = District),
    min.segment.length = 0,
    force = 100,
    nudge_x = 1,
    seed = 10,
    colour = "#256e5d"
  )

figure1
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/figure1.png",
  figure1,
  device = grDevices::png,
  dpi = 300
)
```


## Serological Results

### By communities and type of plasmodium/time
```{r}
# ffi_total %>%
#   ggplot(aes(
#     x = ffi_is_community,
#     fill = pv_historic
#   )) +
#   coord_flip(ylim = c(0, 0.25)) +
#   geom_bar(position = "fill")

plot1 <- ffi_total %>%
  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),
    ffi_is_community = paste0(
      ffi_is_community,
      " (",
      str_to_title(ffi_is_district),
      ")"
    ),
    ffi_is_community = fct_reorder(
      ffi_is_community,
      distance_to_rh_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(
    x = "Communities",
    y = "Percentage"
  ) +
  guides(
    fill = guide_legend("Results")
  ) +  
  theme_minimal() +
  theme(
    strip.text = element_text(
      face = "bold",
      size = 11
    ),
    axis.title = element_text(
      face = "bold",
      size = 11
    ),
    legend.title = element_text(
      face = "bold",
      size = 11
    )
  )

plot1
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot1_typeofmalaria_communities.png",
  plot1,
  dpi = 300,
  bg = "white",
  width = 10,
  height = 9
)
```


### By communities, Indiana and type of plasmodium/time
```{r}
plot1_indiana <- ffi_total %>%
  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_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 = "Communities",
    y = "Percentage"
  ) +
  guides(
    fill = guide_legend("Results")
  ) +  
  theme_minimal() +
  theme(
    strip.text = element_text(
      face = "bold",
      size = 11
    ),
    axis.title = element_text(
      face = "bold",
      size = 11
    ),
    legend.title = element_text(
      face = "bold",
      size = 11
    )
  )

plot1_indiana
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot1_typeofmalaria_communities_indiana.png",
  plot1_indiana,
  dpi = 300,
  bg = "white",
  width = 10,
  height = 9
)
```


### By communities, Belen and type of plasmodium/time
```{r}
plot1_belen <- ffi_total %>%
  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_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 Belen District", 100),
    x = "Communities",
    y = "Percentage"
  ) +
  guides(
    fill = guide_legend("Results")
  ) +  
  theme_minimal() +
  theme(
    strip.text = element_text(
      face = "bold",
      size = 11
    ),
    axis.title = element_text(
      face = "bold",
      size = 11
    ),
    legend.title = element_text(
      face = "bold",
      size = 11
    )
  )

plot1_belen
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot1_typeofmalaria_communities_belen.png",
  plot1_belen,
  dpi = 300,
  bg = "white",
  width = 10,
  height = 9
)
```

### By communities, type of plasmodium exposure
```{r}
plot2 <- ffi_total %>%
  drop_na(pf_exposure, pv_exposure) %>%
  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"
    ),
    ffi_is_community = str_to_title(ffi_is_community),
    ffi_is_community = paste0(
      ffi_is_community,
      " (",
      str_to_title(ffi_is_district),
      ")"
    ),
    ffi_is_community = fct_reorder(
      ffi_is_community,
      distance_to_rh_sampl_num,
      .desc = TRUE
    )
  ) %>%
  ggplot(aes(
    x = ffi_is_community,
    fill = result_exposure
  )) +
  facet_wrap(vars(exposure)) +
  geom_bar(position = "fill", color = "black") +
  coord_flip(ylim = c(0, 0.40)) +
  scale_y_continuous(
    labels = scales::percent_format()
  ) +
  # ggsci::scale_fill_lancet() +
  innovar::scale_fill_innova("npr") +
  labs(
    x = "Communities",
    y = "Percentage"
  ) +
  guides(
    fill = guide_legend("Results")
  ) +
  theme_minimal() +
  theme(
    strip.text = element_text(
      face = "bold",
      size = 11
    ),
    axis.title = element_text(
      face = "bold",
      size = 11
    ),
    legend.title = element_text(
      face = "bold",
      size = 11
    )
  )

plot2
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot2_typeofmalaria_communities.png",
  plot2,
  dpi = 300,
  bg = "white",
  width = 10,
  height = 9
)
```

### By communities, time of plasmodium exposure
```{r}
plot3 <- ffi_total %>%
  drop_na(recent_exposure, historical_exposure) %>%
  pivot_longer(
    cols = recent_exposure:historical_exposure,
    names_to = "exposure",
    values_to = "result_exposure"
  ) %>%
  mutate(
    exposure = case_when(
      exposure == "recent_exposure" ~ "Recent Exposure",
      exposure == "historical_exposure" ~ "Historical Exposure"
    ),
    ffi_is_community = str_to_title(ffi_is_community),
    ffi_is_community = paste0(
      ffi_is_community,
      " (",
      str_to_title(ffi_is_district),
      ")"
    ),
    ffi_is_community = fct_reorder(
      ffi_is_community,
      distance_to_rh_sampl_num,
      .desc = TRUE
    )
  ) %>%
  ggplot(aes(
    x = ffi_is_community,
    fill = result_exposure
  )) +
  facet_wrap(vars(exposure)) +
  geom_bar(position = "fill", color = "black") +
  coord_flip(ylim = c(0, 0.30)) +
  scale_y_continuous(
    labels = scales::percent_format()
  ) +
  # ggsci::scale_fill_lancet() +
  innovar::scale_fill_innova("npr") +
  labs(
    x = "Communities",
    y = "Percentage"
  ) +
  guides(
    fill = guide_legend("Results")
  ) +
  theme_minimal() +
  theme(
    strip.text = element_text(
      face = "bold",
      size = 11
    ),
    axis.title = element_text(
      face = "bold",
      size = 11
    ),
    legend.title = element_text(
      face = "bold",
      size = 11
    )
  )

plot3
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot3_typeofmalaria_communities.png",
  plot3,
  dpi = 300,
  bg = "white",
  width = 10,
  height = 9
)
```


```{r}
# ffi_total %>%
#   drop_na(pf_recent:pv_historic, age_cat) %>%
#   pivot_longer(
#     cols = pf_recent:pv_historic,
#     names_to = "malaria",
#     values_to = "result_malaria"
#   ) %>%
#   mutate(
#     result_malaria = case_when(
#       result_malaria == "Positive" ~ 1,
#       TRUE ~ 0
#     ),
#     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)
#   ) %>%
#   group_by(ffi_is_community, age_cat, malaria) %>%
#   summarise(result_malaria = mean(result_malaria)) %>%
#   ggplot(
#     aes(
#       x = age_cat,
#       y = result_malaria,
#       color = malaria,
#       group = malaria
#     )
#   ) +
#   geom_point() +
#   geom_line() +
#   facet_wrap(vars(ffi_is_community)) +
#   theme_bw()
```

## Seropositivity 4 malaria

### By district
```{r}
sero_4malaria_district <- ffi_total %>%
  drop_na(pf_recent:pv_historic, age_cat) %>%
  pivot_longer(
    cols = pf_recent:pv_historic,
    names_to = "malaria",
    values_to = "result_malaria"
  ) %>%
  mutate(
    result_malaria = case_when(
      result_malaria == "Positive" ~ 1,
      TRUE ~ 0
    ),
    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"
    ),
    across(
      c(ffi_is_community:ffi_is_health_facility_name),
      str_to_title
    )
  ) %>%
  group_by(ffi_is_district, age_cat, malaria) %>%
  summarise(result_malaria = mean(result_malaria)) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_malaria,
      color = malaria,
      group = malaria
    )
  ) +
  geom_point() +
  geom_line() +
  facet_wrap(vars(ffi_is_district)) +
  scale_y_continuous(
    labels = scales::percent_format(),
    limits = c(0, 0.40)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_4malaria_district
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot4_sero_4malaria_district.png",
  sero_4malaria_district,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 7.5
)
```


### By Health Facility

#### 4 Malaria

```{r}
sero_4malaria_hf <- ffi_total %>%
  drop_na(pf_recent:pv_historic, age_cat) %>%
  pivot_longer(
    cols = pf_recent:pv_historic,
    names_to = "malaria",
    values_to = "result_malaria"
  ) %>%
  mutate(
    result_malaria = case_when(
      result_malaria == "Positive" ~ 1,
      TRUE ~ 0
    ),
    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"
    ),
    across(
      c(ffi_is_community:ffi_is_health_facility_name),
      str_to_title
    ),
    ffi_is_health_facility_name = paste(
      ffi_is_district,
      ffi_is_health_facility_name,
      sep = " - "
    )
  ) %>%
  group_by(ffi_is_health_facility_name, age_cat, malaria) %>%
  summarise(result_malaria = mean(result_malaria)) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_malaria,
      color = malaria,
      group = malaria
    )
  ) +
  geom_point() +
  geom_line() +
  facet_wrap(vars(ffi_is_health_facility_name)) +
  scale_y_continuous(
    labels = scales::percent_format()
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") + 
  theme_bw() +
  theme(
    axis.text.x = element_text(
      angle = 45,
      vjust = 1,
      hjust = 1
    )
  )
  
sero_4malaria_hf
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot5_sero_4malaria_hf.png",
  sero_4malaria_hf,
  dpi = 300,
  bg = "white",
  width = 10,
  height = 7
)
```


#### Type of Exposure

```{r}
sero_4malaria_hf_typof_exposure <- ffi_total %>%
  drop_na(pf_exposure, pv_exposure, age_cat) %>%
  pivot_longer(
    cols = pv_exposure:pf_exposure,
    names_to = "exposure",
    values_to = "result_exposure"
  ) %>%
  mutate(
    result_exposure = case_when(
      result_exposure == "Positive" ~ 1,
      TRUE ~ 0
    ),
    exposure = case_when(
         exposure == "pv_exposure" ~ "P. Vivax Exposure",
         exposure == "pf_exposure" ~ "P. Falciparum Exposure"
    ),
    across(
      c(ffi_is_community:ffi_is_health_facility_name),
      str_to_title
    ),
    ffi_is_health_facility_name = paste(
      ffi_is_district,
      ffi_is_health_facility_name,
      sep = " - "
    )
  ) %>%
  group_by(ffi_is_health_facility_name, age_cat, exposure) %>%
  summarise(result_exposure = mean(result_exposure)) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_wrap(vars(ffi_is_health_facility_name)) +
  scale_y_continuous(
    labels = scales::percent_format()
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") + 
  theme_bw() +
  theme(
    axis.text.x = element_text(
      angle = 45,
      vjust = 1,
      hjust = 1
    )
  )
  
sero_4malaria_hf_typof_exposure
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot5_sero_4malaria_hf_typeofexposure.png",
  sero_4malaria_hf_typof_exposure,
  dpi = 300,
  bg = "white",
  width = 10,
  height = 7
)
```

#### Time of Exposure

```{r}
sero_4malaria_hf_timeof_exposure <- ffi_total %>%
  drop_na(recent_exposure, historical_exposure, age_cat) %>%
  pivot_longer(
    cols = recent_exposure:historical_exposure,
    names_to = "exposure",
    values_to = "result_exposure"
  ) %>%
  mutate(
    result_exposure = case_when(
      result_exposure == "Positive" ~ 1,
      TRUE ~ 0
    ),
    exposure = case_when(
      exposure == "recent_exposure" ~ "Recent Exposure",
      exposure == "historical_exposure" ~ "Historical Exposure"
    ),
    across(
      c(ffi_is_community:ffi_is_health_facility_name),
      str_to_title
    ),
    ffi_is_health_facility_name = paste(
      ffi_is_district,
      ffi_is_health_facility_name,
      sep = " - "
    )
  ) %>%
  group_by(ffi_is_health_facility_name, age_cat, exposure) %>%
  summarise(result_exposure = mean(result_exposure)) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_wrap(vars(ffi_is_health_facility_name)) +
  scale_y_continuous(
    labels = scales::percent_format()
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") + 
  theme_bw() +
  theme(
    axis.text.x = element_text(
      angle = 45,
      vjust = 1,
      hjust = 1
    )
  )
  
sero_4malaria_hf_timeof_exposure
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot5_sero_4malaria_hf_timeofexposure.png",
  sero_4malaria_hf_timeof_exposure,
  dpi = 300,
  bg = "white",
  width = 10,
  height = 7
)
```

## Seropositivity by Type of Malaria

### By District

```{r}
ffi_typeof_malaria <- ffi_total %>%
  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
    )
  )

sero_typeofmalaria_district <- ffi_typeof_malaria %>%
  group_by(ffi_is_district, age_cat, exposure) %>%
  summarise(result_exposure = mean(result_exposure)) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_wrap(vars(ffi_is_district)) +
  scale_y_continuous(
    labels = scales::percent_format(),
    limits = c(0, 0.40)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_typeofmalaria_district
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot6_sero_typeofmalaria_district.png",
  sero_typeofmalaria_district,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 7.5
)
```

### By Health Facility

```{r}
sero_typeofmalaria_hf <- ffi_typeof_malaria %>%
  group_by(ffi_is_health_facility_name, age_cat, exposure) %>%
  summarise(result_exposure = mean(result_exposure)) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_wrap(vars(ffi_is_health_facility_name)) +
  scale_y_continuous(
    labels = scales::percent_format()
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_typeofmalaria_hf
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot7_sero_typeofmalaria_hf.png",
  sero_typeofmalaria_hf,
  dpi = 300,
  bg = "white",
  width = 13,
  height = 9
)
```

## Seropositivity by Time of Exposure 

### By District
```{r}
ffi_timeofmalaria <- ffi_total %>%
  drop_na(recent_exposure, historical_exposure, age_cat) %>%
  pivot_longer(
    cols = recent_exposure:historical_exposure,
    names_to = "exposure",
    values_to = "result_exposure"
  ) %>%
  mutate(
    exposure = case_when(
      exposure == "recent_exposure" ~ "Recent Exposure",
      exposure == "historical_exposure" ~ "Historical Exposure"
    ),
    result_exposure = case_when(
      result_exposure == "Positive" ~ 1,
      TRUE ~ 0
    ),
    across(
      c(ffi_is_community:ffi_is_health_facility_name),
      str_to_title
    )
  )

sero_timeofmalaria_district <- ffi_timeofmalaria %>%
  group_by(ffi_is_district, age_cat, exposure) %>%
  summarise(result_exposure = mean(result_exposure)) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_wrap(vars(ffi_is_district)) +
  scale_y_continuous(
    labels = scales::percent_format(),
    limits = c(0, 0.40)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot8_sero_timeofmalaria_district.png",
  sero_timeofmalaria_district,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 7.5
)
```

### By Health Facility

```{r}
sero_timeofmalaria_hf <- ffi_timeofmalaria %>%
  group_by(ffi_is_health_facility_name, age_cat, exposure) %>%
  summarise(result_exposure = mean(result_exposure)) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_wrap(vars(ffi_is_health_facility_name)) +
  scale_y_continuous(
    labels = scales::percent_format()
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") + 
  theme_bw()
  
sero_timeofmalaria_hf
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot9_sero_timeofmalaria_hf.png",
  sero_timeofmalaria_hf,
  dpi = 300,
  bg = "white",
  width = 13,
  height = 9
)
```

## Relation with sex

```{r}
library(rlang)

seropositivy_summarise <- function(data, type, variable, ...) {
  if (type == "4malaria") {
    ffi_total %>%
      drop_na(pf_recent:pv_historic, {{ variable }}, ...) %>%
      pivot_longer(
        cols = pf_recent:pv_historic,
        names_to = "malaria",
        values_to = "result_malaria"
      ) %>%
      mutate(
        result_malaria = case_when(
          result_malaria == "Positive" ~ 1,
          TRUE ~ 0
        ),
        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"
        ),
        across(
          c(ffi_is_community:ffi_is_health_facility_name),
          str_to_title
        )
      ) %>%
      group_by(ffi_is_district, {{ variable }}, ..., malaria) %>%
      summarise(result_malaria = mean(result_malaria))  %>%
      ungroup()
  } else if (type == "typeofmalaria" ) {
    ffi_total %>%
      drop_na(pf_exposure, pv_exposure, {{ variable }}, ...) %>%
      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(ffi_is_district, {{ variable }}, ..., exposure) %>%
      summarise(result_exposure = mean(result_exposure)) %>%
      ungroup()
  } else if (type == "timeofmalaria") {
    ffi_total %>%
      drop_na(recent_exposure, historical_exposure, 
              {{ variable }}, ...) %>%
      pivot_longer(
        cols = recent_exposure:historical_exposure,
        names_to = "exposure",
        values_to = "result_exposure"
      ) %>%
      mutate(
        exposure = case_when(
          exposure == "recent_exposure" ~ "Recent Exposure",
          exposure == "historical_exposure" ~ "Historical 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(ffi_is_district, {{ variable }}, ..., exposure) %>%
      summarise(result_exposure = mean(result_exposure)) %>%
      ungroup()
  }  
}
```

### By 4 Malaria

```{r}
sero_4malaria_district_gender <- seropositivy_summarise(
  ffi_total,
  "4malaria",
  age_cat,
  gender
) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_malaria,
      color = malaria,
      group = malaria
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(gender),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    limits = c(0, 0.40)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_4malaria_district_gender
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot10_sero_4malaria_district_gender.png",
  sero_4malaria_district_gender,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

### By Type of Malaria

```{r}
sero_typeofmalaria_district_gender <- seropositivy_summarise(
  ffi_total,
  "typeofmalaria",
  age_cat,
  gender
) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(gender),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    limits = c(0, 0.45)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_typeofmalaria_district_gender
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot11_sero_typeofmalaria_district_gender.png",
  sero_typeofmalaria_district_gender,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

### By Time of Malaria

```{r}
sero_timeofmalaria_district_gender <- seropositivy_summarise(
  ffi_total,
  "timeofmalaria",
  age_cat,
  gender
) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(gender),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    limits = c(0, 0.45)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_timeofmalaria_district_gender
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot12_sero_timeofmalaria_district_gender.png",
  sero_timeofmalaria_district_gender,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

## Relation with Education Level

### By 4 Malaria

```{r}
sero_4malaria_district_edulevel <- seropositivy_summarise(
  ffi_total,
  "4malaria",
  age_cat,
  education_level
) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_malaria,
      color = malaria,
      group = malaria
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(education_level),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    #limits = c(0, 0.40)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_4malaria_district_edulevel
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot13_sero_4malaria_district_edulevel.png",
  sero_4malaria_district_edulevel,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

### By Type of Malaria

```{r}
sero_typeofmalaria_district_edulevel <- seropositivy_summarise(
  ffi_total,
  "typeofmalaria",
  age_cat,
  education_level
) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(education_level),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    #limits = c(0, 0.45)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_typeofmalaria_district_edulevel
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot14_sero_typeofmalaria_district_edulevel.png",
  sero_typeofmalaria_district_edulevel,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

### By Time of Malaria

```{r}
sero_timeofmalaria_district_edulevel <- seropositivy_summarise(
  ffi_total,
  "timeofmalaria",
  age_cat,
  education_level
) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(education_level),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    #limits = c(0, 0.45)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_timeofmalaria_district_edulevel
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot15_sero_timeofmalaria_district_edulevel.png",
  sero_timeofmalaria_district_edulevel,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

## Relation with Fever Month

#### By 4 malaria
```{r}
ffi_fever_month <- ffi_total %>%
  drop_na(pf_recent:pv_historic, age_cat, ffi_is_fever_month) %>%
  pivot_longer(
    cols = pf_recent:pv_historic,
    names_to = "malaria",
    values_to = "result_malaria"
  ) %>%
  mutate(
    result_malaria = case_when(
      result_malaria == "Positive" ~ 1,
      TRUE ~ 0
    ),
    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"
    ),
    across(
      c(ffi_is_community:ffi_is_health_facility_name),
      str_to_title
    )
  ) %>%
  group_by(ffi_is_fever_month, age_cat, malaria) %>%
  summarise(result_malaria = mean(result_malaria))


sero_4malaria_fever_month <- ffi_fever_month %>%
  ggplot(
    aes(
      x = result_malaria,
      y = age_cat,
      group = age_cat
    )
  ) +
  geom_path(
    color = "#bdbdbd"
  ) +
  geom_point(
    aes(color = malaria),
    size = 3
  ) +
  facet_wrap(vars(ffi_is_fever_month)) +
  scale_x_continuous(
    labels = scales::percent_format(),
    #limits = c(0, 0.5)
  ) +
  labs(
    y = "Age",
    x = "Seropositivity",
    title = str_wrap("Malaria seropositivity by type of exposure, plasmodium and presence of fever in the last month", 100)
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()
  
sero_4malaria_fever_month
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot16_sero_4malaria_fevermonth.png",
  sero_4malaria_fever_month,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```


#### By Type of Malaria

```{r}
ffi_fever_month_typeofmalaria <- ffi_total %>%
  drop_na(pf_recent:pv_historic, age_cat, ffi_is_fever_month) %>%
  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(ffi_is_fever_month, age_cat, exposure) %>%
    summarise(result_exposure = mean(result_exposure))

sero_typeofmalaria_fever_month <- ffi_fever_month_typeofmalaria %>%
  ggplot(
    aes(
      x = result_exposure,
      y = age_cat,
      group = age_cat
    )
  ) +
  geom_path(
    color = "#bdbdbd"
  ) +
  geom_point(
    aes(color = exposure),
    size = 3
  ) +
  facet_wrap(vars(ffi_is_fever_month)) +
  scale_x_continuous(
    labels = scales::percent_format(),
    # limits = c(0, 0.5)
  ) +
  labs(
    y = "Age",
    x = "Seropositivity",
    title = str_wrap("Malaria seropositivity by type of plasmodium and presence of fever in the last month", 100)
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_typeofmalaria_fever_month
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot16_sero_typeofmalaria_fevermonth.png",
  sero_typeofmalaria_fever_month,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

#### By Time of Malaria

```{r}
ffi_fever_month_timeofmalaria <- ffi_total %>%
  drop_na(pf_recent:pv_historic, age_cat, ffi_is_fever_month) %>%
  pivot_longer(
    cols = recent_exposure:historical_exposure,
    names_to = "exposure",
    values_to = "result_exposure"
  ) %>%
  mutate(
    exposure = case_when(
      exposure == "recent_exposure" ~ "Recent Exposure",
      exposure == "historical_exposure" ~ "Historical 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(ffi_is_fever_month, age_cat, exposure) %>%
  summarise(result_exposure = mean(result_exposure))

sero_timeofmalaria_fever_month <- ffi_fever_month_typeofmalaria %>%
  ggplot(
    aes(
      x = result_exposure,
      y = age_cat,
      group = age_cat
    )
  ) +
  geom_path(
    color = "#bdbdbd"
  ) +
  geom_point(
    aes(color = exposure),
    size = 3
  ) +
  facet_wrap(vars(ffi_is_fever_month)) +
  scale_x_continuous(
    labels = scales::percent_format(),
    # limits = c(0, 0.5)
  ) +
  labs(
    y = "Age",
    x = "Seropositivity",
    title = str_wrap("Malaria seropositivity by time of plasmodium and presence of fever in the last month", 100)
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_timeofmalaria_fever_month
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot16_sero_timeofmalaria_fevermonth.png",
  sero_timeofmalaria_fever_month,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```


## Relation with Antimalarial Drugs

### By 4 Malaria

```{r}
sero_4malaria_district_antimaldrugs <- seropositivy_summarise(
  ffi_total,
  "4malaria",
  age_cat,
  ffi_is_antimal_drug_use
) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_malaria,
      color = malaria,
      group = malaria
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(ffi_is_antimal_drug_use),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    limits = c(0, 0.50)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_4malaria_district_antimaldrugs
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot19_sero_4malaria_district_antimaldrugs.png",
  sero_4malaria_district_antimaldrugs,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

### By Type of Malaria

```{r}
sero_typeofmalaria_district_antimaldrugs <- seropositivy_summarise(
  ffi_total,
  "typeofmalaria",
  age_cat,
  ffi_is_antimal_drug_use
) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(ffi_is_antimal_drug_use),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    limits = c(0, 0.60)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_typeofmalaria_district_antimaldrugs
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot20_sero_typeofmalaria_district_antimaldrugs.png",
  sero_typeofmalaria_district_antimaldrugs,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

### By Time of Malaria

```{r}
sero_timeofmalaria_district_antimaldrugs <- seropositivy_summarise(
  ffi_total,
  "timeofmalaria",
  age_cat,
  ffi_is_antimal_drug_use
) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(ffi_is_antimal_drug_use),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    limits = c(0, 0.5)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_timeofmalaria_district_antimaldrugs
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot21_sero_timeofmalaria_district_antimaldrugs.png",
  sero_timeofmalaria_district_antimaldrugs,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```


## Relation with time of someone had malaria

### By general malaria

```{r}
ffi_4malaria_mal_lifetime <- ffi_total %>% 
  drop_na(only_pv_exposure:freedom_malaria, 
          age_code, ffi_is_mal_lifetime) %>% 
  pivot_longer(
    cols = only_pv_exposure:freedom_malaria,
    names_to = "malaria",
    values_to = "result_malaria"
  ) %>% 
  mutate(
    result_malaria = case_when(
      result_malaria == "Positive" ~ 1,
      TRUE ~ 0
    ),
    malaria = case_when(
      malaria == "only_pv_exposure" ~ "Only P. vivax",
      malaria == "only_pf_exposure" ~ "Only P. falciparum",
      malaria == "pv_pf_exposure" ~ "P. vivax & falciparum Exposure",
      malaria == "freedom_malaria" ~ "Freedom from Malaria"
    ),
    across(
      c(ffi_is_community:ffi_is_health_facility_name),
      str_to_title
    ),
    age_code = as.numeric(age_code)
  ) %>%
  group_by(ffi_is_mal_lifetime, age_code, malaria) %>%
  summarise(result_malaria = mean(result_malaria)) %>%
  ungroup()
  

add_na_ffi_4malaria_mal_lifetime <- ffi_4malaria_mal_lifetime %>%
  count(ffi_is_mal_lifetime, age_code) %>% 
  pivot_wider(
    names_from = age_code,
    values_from = n
  ) %>% 
  hacksaw::keep_na(.logic = "OR") %>% 
  select(1, where(is.na)) %>% 
  pivot_longer(
    cols = -1,
    names_to = "age_code",
    values_to = "malaria"
  ) %>% 
  add_case(
    ffi_is_mal_lifetime = "More Than 7 Times",
    age_code = "3"
  ) %>% 
  mutate(
    malaria = NA_character_,
    age_code = as.numeric(age_code),
    result_malaria = 1
  )

sero_4malaria_mal_lifetime <- ffi_4malaria_mal_lifetime %>%
  ggplot(
    aes(
      x = age_code,
      y = result_malaria,
      fill = malaria
    )
  ) +
  geom_area() +
  geom_area(data = add_na_ffi_4malaria_mal_lifetime,
            aes(x = age_code,
                y = result_malaria,
                fill = malaria)) +
  scale_x_continuous(
    limits = c(1, 8),
    breaks = seq(1, 8, 1),
    labels = c(
      "[0-10)",
      "[10-20)",
      "[20-30)",
      "[30-40)",
      "[40-50)",
      "[50-60)",
      "[60-70)",
      "[70+)"
    ),
    expand = c(0, 0)
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    expand = c(0, 0)
  ) +
  facet_wrap(vars(ffi_is_mal_lifetime)) +
  labs(
    y = "Seropositivity",
    x = "Age",
    title = str_wrap("Malaria seropositivity by type of exposure, plasmodium and how many times they think they have had malaria", 100)
  ) +
  guides(
    fill = guide_legend("Malaria")
  ) +
  innovar::scale_fill_innova("npr",
                             na.value = "#a2a2a2") +
  theme_bw() +
  theme(
    axis.text.x = element_text(
      angle = 45,
      vjust = 1,
      hjust = 1
    ),
    panel.spacing = unit(1, "lines")
  )

sero_4malaria_mal_lifetime
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot22_sero_4ofmalaria_mal_lifetime.png",
  sero_4malaria_mal_lifetime,
  dpi = 300,
  bg = "white",
  width = 11,
  height = 5
)
```

### By Recent Malaria

```{r}
ffi_recentmalaria_mal_lifetime <- ffi_total %>% 
  drop_na(only_pv_recent:freedom_malaria_recent, 
          age_code, ffi_is_mal_lifetime) %>% 
  pivot_longer(
    cols = only_pv_recent:freedom_malaria_recent,
    names_to = "malaria",
    values_to = "result_malaria"
  ) %>% 
  mutate(
    result_malaria = case_when(
      result_malaria == "Positive" ~ 1,
      TRUE ~ 0
    ),
    malaria = case_when(
      malaria == "only_pv_recent" ~ "Only P. vivax Recent",
      malaria == "only_pf_recent" ~ "Only P. falciparum Recent",
      malaria == "pv_pf_recent" ~ "P. vivax & falciparum Recent",
      malaria == "freedom_malaria_recent" ~ "Freedom from Malaria Recent"
    ),
    across(
      c(ffi_is_community:ffi_is_health_facility_name),
      str_to_title
    ),
    age_code = as.numeric(age_code)
  ) %>%
  group_by(ffi_is_mal_lifetime, age_code, malaria) %>%
  summarise(result_malaria = mean(result_malaria)) %>%
  ungroup()

add_na_recentmalaria_mal_lifetime <- ffi_recentmalaria_mal_lifetime %>% 
  count(ffi_is_mal_lifetime, age_code) %>% 
  pivot_wider(
    names_from = age_code,
    values_from = n
  ) %>% 
  hacksaw::keep_na(.logic = "OR") %>% 
  select(1, where(is.na)) %>% 
  pivot_longer(
    cols = -1,
    names_to = "age_code",
    values_to = "malaria"
  ) %>% 
  add_case(
    ffi_is_mal_lifetime = "More Than 7 Times",
    age_code = "3"
  ) %>% 
  mutate(
    malaria = NA_character_,
    age_code = as.numeric(age_code),
    result_malaria = 1
  )

sero_recentmalaria_mal_lifetime <- ffi_recentmalaria_mal_lifetime %>%
  ggplot(
    aes(
      x = age_code,
      y = result_malaria,
      fill = malaria
    )
  ) +
  geom_area() +
  geom_area(data = add_na_recentmalaria_mal_lifetime,
            aes(x = age_code,
                y = result_malaria,
                fill = malaria)) +
  scale_x_continuous(
    limits = c(1, 8),
    breaks = seq(1, 8, 1),
    labels = c(
      "[0-10)",
      "[10-20)",
      "[20-30)",
      "[30-40)",
      "[40-50)",
      "[50-60)",
      "[60-70)",
      "[70+)"
    ),
    expand = c(0, 0)
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    expand = c(0, 0)
  ) +
  facet_wrap(vars(ffi_is_mal_lifetime)) +
  labs(
    y = "Seropositivity",
    x = "Age",
    title = str_wrap("Malaria seropositivity by recent exposure and how many times they think they have had malaria", 100)
  ) +
  guides(
    fill = guide_legend("Malaria")
  ) +
  innovar::scale_fill_innova("npr",
                             na.value = "#a2a2a2") +
  theme_bw() +
  theme(
    axis.text.x = element_text(
      angle = 45,
      vjust = 1,
      hjust = 1
    ),
    panel.spacing = unit(1, "lines")
  )

sero_recentmalaria_mal_lifetime
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot23_sero_recent_malaria_mal_lifetime.png",
  sero_recentmalaria_mal_lifetime,
  dpi = 300,
  bg = "white",
  width = 11,
  height = 5
)
```


### By Historic Malaria

```{r}
ffi_historicmalaria_mal_lifetime <- ffi_total %>% 
  drop_na(only_pv_historic:freedom_malaria_historic, 
          age_code, ffi_is_mal_lifetime) %>% 
  pivot_longer(
    cols = only_pv_historic:freedom_malaria_historic,
    names_to = "malaria",
    values_to = "result_malaria"
  ) %>% 
  mutate(
    result_malaria = case_when(
      result_malaria == "Positive" ~ 1,
      TRUE ~ 0
    ),
    malaria = case_when(
      malaria == "only_pv_historic" ~ "Only P. vivax Historic",
      malaria == "only_pf_historic" ~ "Only P. falciparum Historic",
      malaria == "pv_pf_historic" ~ "P. vivax & falciparum Historic",
      malaria == "freedom_malaria_historic" ~ "Freedom from Malaria Historic"
    ),
    across(
      c(ffi_is_community:ffi_is_health_facility_name),
      str_to_title
    ),
    age_code = as.numeric(age_code)
  ) %>%
  group_by(ffi_is_mal_lifetime, age_code, malaria) %>%
  summarise(result_malaria = mean(result_malaria)) %>%
  ungroup()

add_na_historicmalaria_mal_lifetime <- ffi_historicmalaria_mal_lifetime %>% 
  count(ffi_is_mal_lifetime, age_code) %>% 
  pivot_wider(
    names_from = age_code,
    values_from = n
  ) %>% 
  hacksaw::keep_na(.logic = "OR") %>% 
  select(1, where(is.na)) %>% 
  pivot_longer(
    cols = -1,
    names_to = "age_code",
    values_to = "malaria"
  ) %>% 
  add_case(
    ffi_is_mal_lifetime = "More Than 7 Times",
    age_code = "3"
  ) %>% 
  mutate(
    malaria = NA_character_,
    age_code = as.numeric(age_code),
    result_malaria = 1
  )

sero_historicmalaria_mal_lifetime <- ffi_historicmalaria_mal_lifetime %>%
  ggplot(
    aes(
      x = age_code,
      y = result_malaria,
      fill = malaria
    )
  ) +
  geom_area() +
  geom_area(data = add_na_historicmalaria_mal_lifetime,
            aes(x = age_code,
                y = result_malaria,
                fill = malaria)) +
  scale_x_continuous(
    limits = c(1, 8),
    breaks = seq(1, 8, 1),
    labels = c(
      "[0-10)",
      "[10-20)",
      "[20-30)",
      "[30-40)",
      "[40-50)",
      "[50-60)",
      "[60-70)",
      "[70+)"
    ),
    expand = c(0, 0)
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    expand = c(0, 0)
  ) +
  facet_wrap(vars(ffi_is_mal_lifetime)) +
  labs(
    y = "Seropositivity",
    x = "Age",
    title = str_wrap("Malaria seropositivity by historic exposure and how many times they think they have had malaria", 100)
  ) +
  guides(
    fill = guide_legend("Malaria")
  ) +
  innovar::scale_fill_innova("npr",
                             na.value = "#a2a2a2") +
  theme_bw() +
  theme(
    axis.text.x = element_text(
      angle = 45,
      vjust = 1,
      hjust = 1
    ),
    panel.spacing = unit(1, "lines")
  )

sero_historicmalaria_mal_lifetime
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot24_sero_historic_malaria_mal_lifetime.png",
  sero_historicmalaria_mal_lifetime,
  dpi = 300,
  bg = "white",
  width = 11,
  height = 5
)
```


## Relation where someone usually bathe

### By 4 Malaria

```{r}
sero_4malaria_district_ussualybathe <- seropositivy_summarise(
  ffi_total,
  "4malaria",
  age_cat,
  ffi_is_place_shower
) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_malaria,
      color = malaria,
      group = malaria
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(ffi_is_place_shower),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    #limits = c(0, 0.50)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_4malaria_district_ussualybathe
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot25_sero_4malaria_district_ussualybathe.png",
  sero_4malaria_district_ussualybathe,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

### By Type of Malaria

```{r}
sero_typeofmalaria_district_ussualybathe <- seropositivy_summarise(
  ffi_total,
  "typeofmalaria",
  age_cat,
  ffi_is_place_shower
) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(ffi_is_place_shower),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    #limits = c(0, 0.60)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_typeofmalaria_district_ussualybathe
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot26_sero_typeofmalaria_district_ussualybathe.png",
  sero_typeofmalaria_district_ussualybathe,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

### By Time of Malaria

```{r}
sero_timeofmalaria_district_ussualybathe <- seropositivy_summarise(
  ffi_total,
  "timeofmalaria",
  age_cat,
  ffi_is_place_shower
) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(ffi_is_place_shower),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    #limits = c(0, 0.5)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_timeofmalaria_district_ussualybathe
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot27_sero_timeofmalaria_district_ussualybathe.png",
  sero_timeofmalaria_district_ussualybathe,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```


## Relation with those who have a mosquito net

### By 4 Malaria

```{r}
sero_4malaria_district_mosquitnet <- seropositivy_summarise(
  ffi_total,
  "4malaria",
  age_cat,
  ffi_is_mosq_net
) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_malaria,
      color = malaria,
      group = malaria
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(ffi_is_mosq_net),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    #limits = c(0, 0.50)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_4malaria_district_mosquitnet
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot28_sero_4malaria_district_mosquitnet.png",
  sero_4malaria_district_mosquitnet,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

### By Type of Malaria

```{r}
sero_typeofmalaria_district_mosquitnet <- seropositivy_summarise(
  ffi_total,
  "typeofmalaria",
  age_cat,
  ffi_is_mosq_net
) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(ffi_is_mosq_net),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    #limits = c(0, 0.60)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_typeofmalaria_district_mosquitnet
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot29_sero_typeofmalaria_district_mosquitnet.png",
  sero_typeofmalaria_district_mosquitnet,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

### By Time of Malaria

```{r}
sero_timeofmalaria_district_mosquitnet <- seropositivy_summarise(
  ffi_total,
  "timeofmalaria",
  age_cat,
  ffi_is_mosq_net
) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(ffi_is_mosq_net),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    #limits = c(0, 0.5)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_timeofmalaria_district_mosquitnet
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot30_sero_timeofmalaria_district_mosquitnet.png",
  sero_timeofmalaria_district_mosquitnet,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

## Relation with those who have a trip in the last month

### By 4 Malaria

```{r}
sero_4malaria_district_tripmonth <- seropositivy_summarise(
  ffi_total,
  "4malaria",
  age_cat,
  ffi_is_trip_month
) %>%
  mutate(
    ffi_is_trip_month = fct_relevel(
      ffi_is_trip_month,
      "No", "Yes", "Don't Know/No Answer"
    ),
  ) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_malaria,
      color = malaria,
      group = malaria
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(ffi_is_trip_month),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    #limits = c(0, 0.50)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_4malaria_district_tripmonth
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot31_sero_4malaria_district_tripmonth.png",
  sero_4malaria_district_tripmonth,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

### By Type of Malaria

```{r}
sero_typeofmalaria_district_tripmonth <- seropositivy_summarise(
  ffi_total,
  "typeofmalaria",
  age_cat,
  ffi_is_trip_month
) %>%
  mutate(
    ffi_is_trip_month = fct_relevel(
      ffi_is_trip_month,
      "No",
      "Yes",
      "Don't Know/No Answer"
    ),
  ) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(ffi_is_trip_month),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    #limits = c(0, 0.60)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_typeofmalaria_district_tripmonth
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot32_sero_typeofmalaria_district_tripmonth.png",
  sero_typeofmalaria_district_tripmonth,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

### By Time of Malaria

```{r}
sero_timeofmalaria_district_tripmonth <- seropositivy_summarise(
  ffi_total,
  "timeofmalaria",
  age_cat,
  ffi_is_trip_month
) %>%
  mutate(
    ffi_is_trip_month = fct_relevel(
      ffi_is_trip_month,
      "No",
      "Yes",
      "Don't Know/No Answer"
    ),
  ) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(ffi_is_trip_month),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    #limits = c(0, 0.5)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_timeofmalaria_district_tripmonth
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot33_sero_timeofmalaria_district_tripmonth.png",
  sero_timeofmalaria_district_tripmonth,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

# Descriptive 01_data

```{r eval=FALSE}
plot_4_2 <- ffi_total %>%
  filter(ffi_is_malaria %in% c("No", "Yes")) %>%
  drop_na(age_cat, pv_historic) %>%
  mutate(
    malaria_gender = paste(gender, ffi_is_malaria)
  ) %>%
  count(pv_historic, age_cat, malaria_gender) %>%
  mutate(Percentage = n / sum(n)) %>%
  mutate(
    Percentage = case_when(
      str_detect(malaria_gender, "Female") ~ Percentage * -1,
      TRUE ~ Percentage
    )
  ) %>%
  ggplot(aes(
    y = age_cat,
    x = Percentage,
    fill = malaria_gender
  )) +
  geom_col(
    position = position_stack(reverse = TRUE),
    color = "black"
  ) +
  labs(
    x = "Age Group",
    y = NULL,
    title = "Population structure by age groups, gender and malaria transmission"
  ) +
  facet_wrap(vars(pv_historic)) + 
  scale_x_continuous(
    limits = c(-0.25, 0.25),
    breaks = seq(-0.25, 0.25, 0.05),
    labels = c(paste0(seq(25, 0, -5), "%"), paste0(seq(5, 25, 5), "%"))
  ) +
  # innovar::scale_fill_innova("npr") +
  scale_fill_discrete(
    type = innovar::innova_pal("npr")(4),
    limits = c("Male Yes", "Male No", "Female No", "Female Yes")
  ) +
  guides(
    fill = guide_legend("Have you ever \nhad malaria?")
  ) +
  geom_vline(xintercept = 0, size = 1) +
  theme_bw(base_size = 12) +
  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 = 12,
      face = "bold"
    ),
    plot.margin = margin(5, 15, 5, 15)
  )

plot_4_2
```

# Descriptive 01_data

```{r eval=FALSE}
plot_4_2 <- ffi_total %>%
  filter(ffi_is_malaria %in% c("No", "Yes")) %>%
  drop_na(age_cat, pv_historic) %>%
  mutate(
    malaria_gender = paste(gender, ffi_is_malaria)
  ) %>%
  count(pv_historic, age_cat, malaria_gender) %>%
  mutate(Percentage = n / sum(n)) %>%
  mutate(
    Percentage = case_when(
      str_detect(malaria_gender, "Female") ~ Percentage * -1,
      TRUE ~ Percentage
    )
  ) %>%
  ggplot(aes(
    y = age_cat,
    x = Percentage,
    fill = malaria_gender
  )) +
  geom_col(
    position = position_stack(reverse = TRUE),
    color = "black"
  ) +
  labs(
    x = "Age Group",
    y = NULL,
    title = "Population structure by age groups, gender and malaria transmission"
  ) +
  facet_wrap(vars(pv_historic)) + 
  scale_x_continuous(
    limits = c(-0.25, 0.25),
    breaks = seq(-0.25, 0.25, 0.05),
    labels = c(paste0(seq(25, 0, -5), "%"), paste0(seq(5, 25, 5), "%"))
  ) +
  # innovar::scale_fill_innova("npr") +
  scale_fill_discrete(
    type = innovar::innova_pal("npr")(4),
    limits = c("Male Yes", "Male No", "Female No", "Female Yes")
  ) +
  guides(
    fill = guide_legend("Have you ever \nhad malaria?")
  ) +
  geom_vline(xintercept = 0, size = 1) +
  theme_bw(base_size = 12) +
  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 = 12,
      face = "bold"
    ),
    plot.margin = margin(5, 15, 5, 15)
  )

plot_4_2
```

## Sankey_seropositive

```{r}
#devtools::install_github("davidsjoberg/ggsankey")
library(ggsankey)
```


### Recent exposure any type of plasmodium

```{r}
ffi_format_sankey <- ffi_total %>%
  mutate(
    ffi_is_access_malaria_yn_hf = factor(
      ffi_is_access_malaria_yn_hf,
      label = c("No", "Yes")
    ),
    ffi_is_mal_lifetime = as.character(ffi_is_mal_lifetime),
    ffi_is_mal_lifetime = replace_na(ffi_is_mal_lifetime, "0 times"),
    ffi_is_mal_lifetime = factor(
      ffi_is_mal_lifetime,
      labels = c(
        "0 times",
        "1 to 3 times",
        "3 to 7 times",
        "More than 7 times"
      )
    )
  )

sankey_seropositive_answ1 <- ffi_format_sankey %>%
  pivot_longer(
    cols = c(recent_exposure, ffi_is_mal_lifetime, ffi_is_access_malaria_yn_hf),
    names_to = "malaria_behavior1",
    values_to = "malaria_answ1"
  ) %>%
  count(malaria_behavior1, malaria_answ1) %>%
  group_by(malaria_behavior1) %>%
  mutate(
    percentage = scales::percent(
      n / sum(n),
      accuracy = 0.1
    )
  ) %>%
  ungroup() %>%
  mutate(
    malaria_behavior1 = fct_relevel(
      malaria_behavior1,
      "ffi_is_mal_lifetime",
      "ffi_is_access_malaria_yn_hf",
      "recent_exposure"
    ),
    malaria_percent1 = paste0(
      malaria_answ1,
      paste0("\n(", percentage, ")")
    ),
    malaria_percent1 = as_factor(malaria_percent1)
  ) %>%
  select(malaria_behavior1, malaria_answ1, malaria_percent1)


sankey_seropositive <- ffi_format_sankey %>%
  drop_na(ffi_is_mal_lifetime, 
  ffi_is_access_malaria_yn_hf,
  recent_exposure
  ) %>%
  make_long(
    ffi_is_mal_lifetime,
    ffi_is_access_malaria_yn_hf, 
    recent_exposure
  ) %>%
  left_join(
    sankey_seropositive_answ1,
      by = c(
        "x" = "malaria_behavior1",
        "node" = "malaria_answ1"
      )
  ) %>%
  mutate(
    node = malaria_percent1,
    node = fct_rev(node)
  ) %>%
  select(-malaria_percent1) %>%
  left_join(
    sankey_seropositive_answ1,
    by = c(
      "next_x" = "malaria_behavior1",
      "next_node" = "malaria_answ1"
    )
  ) %>%
  mutate(next_node = malaria_percent1) %>%
  select(-malaria_percent1) %>%
  ggplot(aes(
    x = x,
    next_x = next_x,
    node = node,
    next_node = next_node,
    fill = factor(node),
    label = node
  )) +
  geom_sankey(
    flow.alpha = .8,
    node.color = "gray30"
  ) +
  geom_sankey_label(size = 4, color = "white", fill = "gray30") +
  scale_x_discrete(
    labels = str_wrap(
      c(
        "How many times do you think you have had malaria in your life?",
        "If you had malaria, would you go to the health facility?",
        "Recent exposure to any type of plasmodium"
      ),
      30
    )
  ) +
  theme_sankey(base_size = 18) +
  innovar::scale_fill_innova("npr") +
  labs(
    x = NULL,
    title = "Behaviour associated with recent malaria exposure"
  ) +
  theme(
    legend.position = "none",
    plot.title = element_text(hjust = .5),
    plot.margin = margin(5, 5, 5, 5)
  )

sankey_seropositive
```

```{r eval=FALSE}
ggsave("./02_output/plots/plot34_sankey_seropositive.png",
       sankey_seropositive,
       height = 7,
       width = 14,
       dpi = 300)
```


### Recent exposure Falciparum

```{r}
ffi_format_sankey <- ffi_total %>%
  mutate(
    ffi_is_access_malaria_yn_hf = factor(
      ffi_is_access_malaria_yn_hf,
      label = c("No", "Yes")
    ),
    ffi_is_mal_lifetime = as.character(ffi_is_mal_lifetime),
    ffi_is_mal_lifetime = replace_na(ffi_is_mal_lifetime, "0 times"),
    ffi_is_mal_lifetime = factor(
      ffi_is_mal_lifetime,
      labels = c(
        "0 times",
        "1 to 3 times",
        "3 to 7 times",
        "More than 7 times"
      )
    )
  )

sankey_seropositive_answ1 <- ffi_format_sankey %>%
  pivot_longer(
    cols = c(pf_recent, ffi_is_mal_lifetime, ffi_is_access_malaria_yn_hf),
    names_to = "malaria_behavior1",
    values_to = "malaria_answ1"
  ) %>%
  count(malaria_behavior1, malaria_answ1) %>%
  group_by(malaria_behavior1) %>%
  mutate(
    percentage = scales::percent(
      n / sum(n),
      accuracy = 0.1
    )
  ) %>%
  ungroup() %>%
  mutate(
    malaria_behavior1 = fct_relevel(
      malaria_behavior1,
      "ffi_is_mal_lifetime",
      "ffi_is_access_malaria_yn_hf",
      "pf_recent"
    ),
    malaria_percent1 = paste0(
      malaria_answ1,
      paste0("\n(", percentage, ")")
    ),
    malaria_percent1 = as_factor(malaria_percent1)
  ) %>%
  select(malaria_behavior1, malaria_answ1, malaria_percent1)


sankey_seropositive_falciparum <- ffi_format_sankey %>%
  drop_na(ffi_is_mal_lifetime, 
          ffi_is_access_malaria_yn_hf,
          pf_recent) %>%
  make_long(
    ffi_is_mal_lifetime,
    ffi_is_access_malaria_yn_hf, 
    pf_recent
  ) %>%
  left_join(
    sankey_seropositive_answ1,
      by = c(
        "x" = "malaria_behavior1",
        "node" = "malaria_answ1"
      )
  ) %>%
  mutate(
    node = malaria_percent1,
    node = fct_rev(node)
  ) %>%
  select(-malaria_percent1) %>%
  left_join(
    sankey_seropositive_answ1,
    by = c(
      "next_x" = "malaria_behavior1",
      "next_node" = "malaria_answ1"
    )
  ) %>%
  mutate(next_node = malaria_percent1) %>%
  select(-malaria_percent1) %>%
  ggplot(aes(
    x = x,
    next_x = next_x,
    node = node,
    next_node = next_node,
    fill = factor(node),
    label = node
  )) +
  geom_sankey(
    flow.alpha = .8,
    node.color = "gray30"
  ) +
  geom_sankey_label(size = 4, color = "white", fill = "gray30") +
  scale_x_discrete(
    labels = str_wrap(
      c(
        "How many times do you think you have had malaria in your life?",
        "If you had malaria, would you go to the health facility?",
        "Recent exposure to Plasmodium Falciparum"
      ),
      30
    )
  ) +
  theme_sankey(base_size = 18) +
  innovar::scale_fill_innova("npr") +
  labs(
    x = NULL,
    title = "Behaviour associated with recent malaria exposure (Falciparum Plasmodium)"
  ) +
  theme(
    legend.position = "none",
    plot.title = element_text(hjust = .5),
    plot.margin = margin(5, 5, 5, 5)
  )

sankey_seropositive_falciparum
```

```{r eval=FALSE}
ggsave("./02_output/plots/plot34_sankey_seropositive_falciparum.png",
       sankey_seropositive_falciparum,
       height = 7,
       width = 14,
       dpi = 300)
```

### Recent exposure Vivax

```{r}
sankey_seropositive_answ1 <- ffi_format_sankey %>%
  pivot_longer(
    cols = c(pv_recent, ffi_is_mal_lifetime, ffi_is_access_malaria_yn_hf),
    names_to = "malaria_behavior1",
    values_to = "malaria_answ1"
  ) %>%
  count(malaria_behavior1, malaria_answ1) %>%
  group_by(malaria_behavior1) %>%
  mutate(
    percentage = scales::percent(
      n / sum(n),
      accuracy = 0.1
    )
  ) %>%
  ungroup() %>%
  mutate(
    malaria_behavior1 = fct_relevel(
      malaria_behavior1,
      "ffi_is_mal_lifetime",
      "ffi_is_access_malaria_yn_hf",
      "pv_recent"
    ),
    malaria_percent1 = paste0(
      malaria_answ1,
      paste0("\n(", percentage, ")")
    ),
    malaria_percent1 = as_factor(malaria_percent1)
  ) %>%
  select(malaria_behavior1, malaria_answ1, malaria_percent1)


sankey_seropositive_vivax <- ffi_format_sankey %>%
  drop_na(ffi_is_mal_lifetime, 
          ffi_is_access_malaria_yn_hf,
          pv_recent) %>%
  make_long(
    ffi_is_mal_lifetime,
    ffi_is_access_malaria_yn_hf, 
    pv_recent
  ) %>%
  left_join(
    sankey_seropositive_answ1,
      by = c(
        "x" = "malaria_behavior1",
        "node" = "malaria_answ1"
      )
  ) %>%
  mutate(
    node = malaria_percent1,
    node = fct_rev(node)
  ) %>%
  select(-malaria_percent1) %>%
  left_join(
    sankey_seropositive_answ1,
    by = c(
      "next_x" = "malaria_behavior1",
      "next_node" = "malaria_answ1"
    )
  ) %>%
  mutate(next_node = malaria_percent1) %>%
  select(-malaria_percent1) %>%
  ggplot(aes(
    x = x,
    next_x = next_x,
    node = node,
    next_node = next_node,
    fill = factor(node),
    label = node
  )) +
  geom_sankey(
    flow.alpha = .8,
    node.color = "gray30"
  ) +
  geom_sankey_label(size = 4, color = "white", fill = "gray30") +
  scale_x_discrete(
    labels = str_wrap(
      c(
        "How many times do you think you have had malaria in your life?",
        "If you had malaria, would you go to the health facility?",
        "Recent exposure to Plasmodium Vivax"
      ),
      30
    )
  ) +
  theme_sankey(base_size = 18) +
  innovar::scale_fill_innova("npr") +
  labs(
    x = NULL,
    title = "Behaviour associated with recent malaria exposure (Vivax Plasmodium)"
  ) +
  theme(
    legend.position = "none",
    plot.title = element_text(hjust = .5),
    plot.margin = margin(5, 5, 5, 5)
  )

sankey_seropositive_vivax
```

```{r eval=FALSE}
ggsave("./02_output/plots/plot34_sankey_seropositive_vivax.png",
       sankey_seropositive_vivax,
       height = 7,
       width = 14,
       dpi = 300)
```

## Relation gender and ocupation

### By 4 Malaria

```{r}
structure_df_district_gender_economic <- SimDesign::createDesign(
  ffi_is_district = c("Belen", "Indiana"),
  economic_activities = factor(c("Day labourer", 
                                 "Wood extractor", "Fisherman", "Livestock farmer", "Farmer", 
                                 "Trader", "Housewife", "Student", "Motorcycle taxi driver", "None", 
                                 "Other")),
  gender = factor(c("Male", "Female")),
  malaria = c("Historical P. Falciparum", "Historical P. Vivax", 
              "Recent P. Falciparum", 
              "Recent P. Vivax")
) 

add_cases_0_district_gender_economic <- structure_df_district_gender_economic %>% 
  anti_join(
    seropositivy_summarise(
      ffi_total,
      "4malaria",
      economic_activities,
      gender
    ) %>% 
      select(-result_malaria)
  ) %>% 
  mutate(
    result_malaria = 0
  )
```


```{r}
sero_4malaria_district_gender_economic <- seropositivy_summarise(
  ffi_total,
  "4malaria",
  economic_activities,
  gender
) %>% 
  bind_rows(add_cases_0_district_gender_economic) %>%
  ggplot(
    aes(
      x = economic_activities,
      y = result_malaria,
      color = malaria,
      group = malaria
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(gender),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    limits = c(0, 0.40)
  ) +
  labs(
    x = "Economic Activities",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw() +
  theme(
    axis.text.x = element_text(
      angle = 45,
      vjust = 1,
      hjust = 1
    )
  )

sero_4malaria_district_gender_economic
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot35_sero_4malaria_district_gender_economic.png",
  sero_4malaria_district_gender_economic,
  dpi = 300,
  bg = "white",
  width = 11,
  height = 7.5,
  scale = 0.9
)
```

### By Type of Malaria

```{r}
sero_typeofmalaria_district_gender <- seropositivy_summarise(
  ffi_total,
  "typeofmalaria",
  age_cat,
  gender
) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(gender),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    limits = c(0, 0.45)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_typeofmalaria_district_gender
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot11_sero_typeofmalaria_district_gender.png",
  sero_typeofmalaria_district_gender,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

### By Time of Malaria

```{r}
sero_timeofmalaria_district_gender <- seropositivy_summarise(
  ffi_total,
  "timeofmalaria",
  age_cat,
  gender
) %>%
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(gender),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    limits = c(0, 0.45)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity"
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_innova("npr") +
  theme_bw()

sero_timeofmalaria_district_gender
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot12_sero_timeofmalaria_district_gender.png",
  sero_timeofmalaria_district_gender,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```



## Relation with distance category towards HF and Regional Hospital

### Method Census for distance categories

#### Towards HF

```{r}
sero_timeofmalaria_district_distance_hf_census <- seropositivy_summarise(
    ffi_total,
    "timeofmalaria",
    age_cat,
    distance_to_hf_census
  ) %>%
  ungroup() %>% 
  mutate(
    distance_to_hf_census = fct_relevel(distance_to_hf_census, 
                                 "Proximate", "Moderate",
                                 "Distant", "Extra Distant")
  ) %>% 
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(distance_to_hf_census),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    #limits = c(0, 0.50)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity",
    title = str_wrap("Seropositivity according to Time of Exposure to malaria by District and distance category of the communities to their respective Health Facilities on a Population Basis.", 80)
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_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 = 12,
      face = "bold"
    )
  )

sero_timeofmalaria_district_distance_hf_census
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot36_sero_timeofmalaria_district_distance_hf_census.png",
  sero_timeofmalaria_district_distance_hf_census,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

### Towards Regional Hospital

```{r}
sero_timeofmalaria_district_distance_rh_census <- seropositivy_summarise(
    ffi_total,
    "timeofmalaria",
    age_cat,
    distance_to_rh_census
  ) %>%
  ungroup() %>% 
  mutate(
    distance_to_rh_census = fct_relevel(distance_to_rh_census, 
                                      "Proximate", "Moderate",
                                      "Distant", "Extra Distant")
  ) %>% 
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(distance_to_rh_census),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    #limits = c(0, 0.50)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity",
    title = str_wrap("Seropositivity according to Time of Exposure to malaria by District and distance category of the communities to Loreto Regional Hospital on a Population Basis", 80)
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_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 = 12,
      face = "bold"
    )
  )

sero_timeofmalaria_district_distance_rh_census
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot37_sero_timeofmalaria_district_distance_rh_census.png",
  sero_timeofmalaria_district_distance_rh_census,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

### Method Sampling System for distance categories

#### Towards HF

```{r}
sero_timeofmalaria_district_distance_hf_sampling <- seropositivy_summarise(
    ffi_total,
    "timeofmalaria",
    age_cat,
    distance_to_hf_sampl
  ) %>%
  ungroup() %>% 
  mutate(
    distance_to_hf_sampl = fct_relevel(distance_to_hf_sampl, 
                                 "Proximate", "Moderate",
                                 "Distant", "Extra Distant")
  ) %>% 
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(distance_to_hf_sampl),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    #limits = c(0, 0.50)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity",
    title = str_wrap("Seropositivity according to Time of Exposure to malaria by District and distance category of the communities to their respective Health Facilities on a Sampling Basis.", 80)
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_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 = 12,
      face = "bold"
    )
  )

sero_timeofmalaria_district_distance_hf_sampling
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot38_sero_timeofmalaria_district_distance_hf_sampling.png",
  sero_timeofmalaria_district_distance_hf_sampling,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

### Towards Regional Hospital

```{r}
sero_timeofmalaria_district_distance_rh_sampling <- seropositivy_summarise(
    ffi_total,
    "timeofmalaria",
    age_cat,
    distance_to_rh_sampl
  ) %>%
  ungroup() %>% 
  mutate(
    distance_to_rh_sampl = fct_relevel(distance_to_rh_sampl, 
                                      "Proximate", "Moderate",
                                      "Distant", "Extra Distant")
  ) %>% 
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(distance_to_rh_sampl),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    #limits = c(0, 0.50)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity",
    title = str_wrap("Seropositivity according to Time of Exposure to malaria by District and distance category of the communities to Loreto Regional Hospital on a Sampling Basis", 80)
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_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 = 12,
      face = "bold"
    )
  )

sero_timeofmalaria_district_distance_rh_sampling
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot39_sero_timeofmalaria_district_distance_rh_sampling.png",
  sero_timeofmalaria_district_distance_rh_sampling,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

### Towards Regional Hospital and Health Facilities

```{r}
sero_time_age_distance <- seropositivy_summarise(
    ffi_total,
    "timeofmalaria",
    age_cat,
    distance_to_rh_hf_sampl
  ) %>%
  ungroup() %>%
  mutate(
    distance_to_rh_hf_sampl = fct_relevel(
      distance_to_rh_hf_sampl,
      "Proximate",
      "Moderate",
      "Distant",
      "Extra Distant"
    )
  ) 
  

sero_time_age_distance %>%
  group_by(ffi_is_district, exposure) %>%
  summarise(
    broom::tidy(
      aov(result_exposure ~ age_cat + distance_to_rh_hf_sampl)
    )
  )
  # ggpubr::ggline(
  #   x = "age_cat",
  #   y = "result_exposure",
  #   fill = "distance_to_rh_hf_sampl",
  #   add = c("mean_se", "dotplot")
  # ) +
  # facet_grid(
  #   vars(ffi_is_district),
  #   vars(exposure)
  # )
```

```{r}
sero_timeofmalaria_district_distance_rh_hf_sampling <- seropositivy_summarise(
    ffi_total,
    "timeofmalaria",
    age_cat,
    distance_to_rh_hf_sampl
  ) %>%
  ungroup() %>% 
  mutate(
    distance_to_rh_hf_sampl = fct_relevel(distance_to_rh_hf_sampl, 
                                      "Proximate", "Moderate",
                                      "Distant", "Extra Distant")
  ) %>% 
  ggplot(
    aes(
      x = age_cat,
      y = result_exposure,
      color = exposure,
      group = exposure
    )
  ) +
  geom_point() +
  geom_line() +
  facet_grid(
    vars(distance_to_rh_hf_sampl),
    vars(ffi_is_district)    
  ) +
  scale_y_continuous(
    labels = scales::percent_format(),
    #limits = c(0, 0.50)
  ) +
  labs(
    x = "Age",
    y = "Seropositivity",
    title = str_wrap("Seropositivity according to Time of Exposure to malaria by District and distance category of the communities to Loreto Regional Hospital and Health Facilities on a Sampling Basis", 80)
  ) +
  guides(
    color = guide_legend("Malaria")
  ) +
  innovar::scale_color_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 = 12,
      face = "bold"
    )
  )

sero_timeofmalaria_district_distance_rh_hf_sampling
```

```{r eval=FALSE}
ggsave(
  "./02_output/plots/plot39_5_sero_timeofmalaria_district_distance_rh_hf_sampling.png",
  sero_timeofmalaria_district_distance_rh_hf_sampling,
  dpi = 300,
  bg = "white",
  width = 12,
  height = 8.5
)
```

## Treemap by Communities

```{r}
library(treemapify)
```

```{r fig.height = 7, fig.width = 10, dpi = 300, fig.align = "center"}
district_label <- ffi_total %>% 
  count(ffi_is_district) %>% 
  mutate(label = paste0(ffi_is_district, 
                        " (n = ", n, ")")) %>% 
  pull(label)

names(district_label) <- unique(ffi_total$ffi_is_district)

plot_40 <- ffi_total %>% 
  count(ffi_is_district, 
        ffi_is_community,
        ffi_is_cod_com) %>% 
  mutate(
    label_com = paste(str_to_title(ffi_is_community),
                      paste0("(",  n, ")"), sep = "\n"),
    label_com = str_wrap(label_com, 10)
  ) %>% 
  ggplot(aes(area = n,
             fill = n,
             label = label_com,
             subgroup = ffi_is_district)) +
  geom_treemap() +
  geom_treemap_text(colour = "black",
                    place = "centre") +
  # geom_treemap_subgroup_text(place = "top", 
  #                            grow = TRUE, alpha = 1, 
  #                            colour = "black", fontface = "italic",
  #                            min.size = 0,
  #                            padding.y = grid::unit(-50, "mm")) +
  geom_treemap_subgroup_border() +
  innovar::scale_fill_innova("npr",
                             discrete = FALSE) +
  facet_wrap(vars(ffi_is_district),
             labeller = labeller(ffi_is_district = district_label)) +
  labs(
    title = "Sample size for each community by district"
  ) +
  theme_classic() +
  theme(
    plot.title = element_text(size = 14,
                              hjust = 0.5,
                              face = "bold"),
    strip.text = element_text(size = 14,
                              face = "bold")
  )

plot_40
```


```{r eval=FALSE}
ggsave("./02_output/plots/40. Sample size for each community by district.png",
       plot = plot_40,
       height = 7,
       width = 10,
       dpi = 300)
```

## General Population by age plot


```{r fig.height = 7, fig.width = 10, dpi = 300, fig.align = "center"}
plot_41 <- apyramid::age_pyramid(data = ffi_total,
                                 age_group = "age_cat",
                                 split_by = "gender",
                                 proportional = TRUE) +
  # scale_fill_brewer(palette = "Set2") +
  scale_fill_manual(
    values = c("#4cd477", "#dc8131")
  ) +
  labs(
    x = "Age Group",
    y = NULL,
    title = "Population structure by age groups and sex"
  ) +
  guides(
    fill = guide_legend("Gender")
  ) +
  innovar::scale_fill_innova("npr") +
  theme_bw(base_size = 12) +
  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 = 12,
                              face = "bold"),
    plot.margin = margin(5, 15, 5, 15)
  )
    
plot_41
```

```{r eval=FALSE}
ggsave("./02_output/plots/41. Population structure by age groups and sex.png",
       plot = plot_41,
       height = 7,
       width = 10,
       dpi = 300)
```

## Economic activity for each gender by education level and gender

```{r fig.height = 8, fig.width = 14, dpi = 300, fig.align = "center"}
na_education_activities <- ffi_total %>% 
  hacksaw::keep_na(education_level,
                   economic_activities,
                   .logic = "OR") %>% 
  count()

label_x_economics <- ffi_total %>% 
  drop_na(education_level, economic_activities) %>% 
  count(gender, economic_activities,
        .drop = FALSE) %>% 
  mutate(
    economic_activities = fct_reorder(economic_activities, n,
                                           .desc = TRUE)
  )

facet_label_gender <- ffi_total %>% 
  count(gender) %>% 
  mutate(label = paste0(gender, " (n = ", n, ")")) %>% 
  pull(label)

names(facet_label_gender) <- unique(ffi_total$gender)


add_economics_na <- label_x_economics %>% 
  filter(n == 0) %>% 
  mutate(education_level = NA,
         percentage = 1)


descriptive_economics_education <- ffi_total %>% 
  drop_na(education_level, economic_activities) %>% 
  count(gender, economic_activities, education_level,
        .drop = FALSE) %>% 
  group_by(gender, economic_activities) %>% 
  mutate(percentage = n/sum(n)) %>% 
  ungroup() %>% 
  mutate(
    economic_activities = fct_relevel(economic_activities,
                                      levels(label_x_economics$economic_activities))
  ) %>% 
  anti_join(
    add_economics_na %>% 
      select(gender, economic_activities)
  ) %>% 
  bind_rows(
    add_economics_na
  )


plot_42 <- descriptive_economics_education %>% 
  ggplot(aes(x = economic_activities,
             y = percentage,
             fill = education_level)) +
  geom_col() +
  facet_wrap(vars(gender),
             labeller = labeller(gender = facet_label_gender)) +
  # scale_y_continuous(labels = scales::percent_format()) +
  scale_y_continuous(
    breaks = c(0, 0.25, 0.5, 0.75, 1, 1.05),
    labels = c(scales::percent(seq(0, 1, 0.25)), "n")
  ) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 18)) +
  labs(
    x = NULL,
    y = NULL,
    title = str_wrap("Main economic activity for each gender by education level and gender",
                     40),
    caption = paste0("Note: ", na_education_activities, " data were excluded due to NA")
  ) +
  guides(
    fill = guide_legend("Education Level")
  ) +
  geom_text(data = label_x_economics,
            aes(x = economic_activities,
                y = 1.05,
                label = n),
            inherit.aes = FALSE) +
  innovar::scale_fill_innova("npr",
                             breaks = levels(ffi_total$education_level),
                             na.value = "#a2a2a2") +
  theme_void(base_size = 12) +
  theme(
    plot.title = element_text(size = 14,
                              hjust = 0.5,
                              face = "bold"),
    strip.text = element_text(size = 12,
                              face = "bold"),
    axis.text.y = element_text(face = c(rep("plain", 5), "bold"),
                               size = c(rep(10, 5), 14)),
    axis.text.x = element_text(size = 12,
                               angle = 45,
                               vjust = 1,
                               hjust = 1),
    plot.margin = margin(5, 15, 5, 15)
  )

plot_42
```


```{r eval=FALSE}
ggsave("./02_output/plots/42. Main economic activity for each gender by education level and gender.png",
       plot = plot_42,
       height = 8,
       width = 14,
       dpi = 300,
       bg = "white")
```

# Maps

```{r}
library(sf)
library(leaflet)
```

```{r}
ffi_total_gps <- ffi_total %>%
  select(
    ffi_is_cod_com:ffi_is_cod_ind,
    pf_recent:pv_historic,
    pv_exposure:historical_exposure
  ) %>%
  left_join(
    ffi_household %>%
      select(
        ffi_h_code_community:ffi_h_code_household,
        ffi_h_community
      ),
    by = c(
      "ffi_is_cod_com" = "ffi_h_code_community",
      "ffi_is_cod_household" = "ffi_h_code_household"
    )
  ) 

ffi_household_gps <- ffi_total_gps %>%
  select(
    ffi_is_cod_com:ffi_is_cod_household,
    ffi_h_community,
    recent_exposure
  ) %>%
  mutate(
    recent_exposure = case_when(
      recent_exposure == "Positive" ~ 1,
      TRUE ~ 0
    )
  ) %>%
  group_by(
    across(c(ffi_is_cod_com:ffi_h_community))
  ) %>%
  summarise(
    recent_exposure = mean(recent_exposure)
  ) %>%
  ungroup()
  # mutate(
  #   color_marker = case_when(
  #     recent_exposure == 0 ~ "red",
  #     TRUE ~ "black"
  #   )
  # )

ffi_household_gps <- ffi_household_gps %>%
  left_join(
    ffi_household %>%
      select(
        ffi_h_code_community:ffi_h_code_household,
        ffi_h_district,
        ffi_gps_long,
        ffi_gps_lat
      ),
    by = c(
      "ffi_is_cod_com" = "ffi_h_code_community",
      "ffi_is_cod_household" = "ffi_h_code_household"
    )
  ) %>%
  mutate(
    ffi_h_household_id = paste0(ffi_is_cod_com, ffi_is_cod_household),
    ffi_is_cod_com_label = paste0(
      "Cod: ",
      ffi_is_cod_household,
      " - ",
      str_to_title(ffi_h_community),
      " (",
      str_to_title(ffi_h_district),
      ")"
    )
  ) %>%
  st_as_sf(
    coords = c("ffi_gps_long", "ffi_gps_lat"),
    crs = 4326
  ) 
  

  
# colors_pal <- rev(innovar:::innova_palettes[["ecomst"]])
# pal <- colorNumeric(
#   colorRamp(colors_pal),
#   ffi_household_gps$recent_exposure
# )

pal <- colorNumeric(
  hcl.colors(17, palette = "zissou"),
  ffi_household_gps$recent_exposure
)

# icons_household <- awesomeIcons(
#   icon = "ios-home",
#   iconColor = "black",
#   library = "ion",
#   markerColor = innovar::innova_pal("ecomst", reverse = TRUE)(17)
# )

communities_sf <- communities_sf %>%
  mutate(
    ffi_h_community = str_to_title(ffi_h_community)
  ) %>%
  left_join(
    ffi_household %>%
      select(ffi_h_district, ffi_h_code_community) %>%
      distinct()
  )  %>%
  mutate(
    ffi_h_community_label = paste0(
      ffi_h_community,
      " (",
      str_to_title(ffi_h_district),
      ")"
    )
  )

  geom_sf(data = hf_shp %>% 
            mutate(Type = "Health Facility (9)",
                   ffi_h_health_facility_name = str_to_title(ffi_h_health_facility_name)),
          aes(shape = Type,
              #color = Type,
              color = ffi_h_health_facility_name),
          size = 5,
          inherit.aes = FALSE
          ) 


household_communities_leaft <- ffi_household_gps %>%
  leaflet() %>%
  addProviderTiles(
    providers$OpenStreetMap,
    group = "OpenStreetMap"
  ) %>%
  addCircleMarkers(
    layerId = ~ffi_h_household_id,
    label = ~ffi_is_cod_com_label,
    color = ~ pal(recent_exposure),
    group = "Households",
    radius = 7,
    weight = 5,
    opacity = 0.6,
    fillOpacity = 0.1,
    labelOptions = labelOptions(
      style = list(
        "font-weight" = "bold",
        padding = "3px 8px"
      ),
      textsize = "12px",
      direction = "auto"
    )
  ) %>%
  # addLegend(
  #   title = "Recent Exposure",
  #   pal = pal, 
  #   values = ~ recent_exposure, 
  #   group = "circles", 
  #   position = "bottomleft",
  #   transform = ~ scales::percent_format(),
  #   opacity = 1
  # ) %>%<font-awesome-icon icon="fa-solid fa-location-dot" />
  leaflegend::addLegendNumeric(
    title = "Recent Exposure",
    pal = pal,
    values = ~ recent_exposure,
    position = "bottomright",
    orientation = "horizontal",
    height = 20,
    width = 100,
    decreasing = FALSE,
    numberFormat = scales::percent_format(),
    group = "circles"
  )  %>%
  addAwesomeMarkers(
    data = communities_sf,
    layerId = ~ffi_h_code_community,
    label = ~ffi_h_community_label,
    group = "Communities"
  ) %>%
  addCircleMarkers(
    data = hf_shp %>% 
      mutate(
        ffi_h_health_facility_name = str_to_title(ffi_h_health_facility_name)
      ),
    layerId = ~ffi_h_health_facility_name,
    label = ~ffi_h_health_facility_name,
    group = "Health Facilities"
  ) %>% 
  addProviderTiles(
    providers$OpenStreetMap,
    group = "OpenStreetMap"
  ) %>%
  addTiles(
    urlTemplate = "http://mt0.google.com/vt/lyrs=m&hl=en&x={x}&y={y}&z={z}&s=Ga",
    attribution = "Google Maps",
    options = leaflet::tileOptions(
      maxNativeZoom = 19,
      maxZoom = 20
    ),
    group = "Google Maps"
  ) %>%
  addTiles(
    urlTemplate = "https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}",
    attribution = "Satellite View",
    options = leaflet::tileOptions(
      maxNativeZoom = 19,
      maxZoom = 20
    ),
    group = "Satellite View"
  ) %>%
  # Layers control*
  addLayersControl(
    overlayGroups = c("Households", "Communities", "Health Facilities"),
    baseGroups = c("OpenStreetMap", "Google Maps", "Satellite View"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  # addLegend(
  #   "topright",
  #   pal = innovar::innova_pal("ecomst", reverse = TRUE)(980)
  # )  %>%
  # addLegend(
  #   position = "bottomright",
  #   colors = c(
  #     "white; width:15px; height:15px;
  #                border:5px solid red; border-radius:50%;",
  #     "white; width:7.5px; height:7.5px; margin-top: 5px;
  #                border:5px solid black; border-radius:50%; margin-left:2px;"
  #   ),
  #   labels = c(
  #     "<div style='display: inline-block; height: 10px;
  #                margin-top: 8px;line-height: 10px;font-weight: bold;
  #                color: black; '>At least 1 person with recent malaria</div>",
  #     "<div style='display: inline-block;height: 10px;
  #                margin-top: 8px;line-height: 10px;font-weight: bold;
  #                color: black; margin-top: 10px;
  #                margin-left:5px'>No one with recent malaria at household</div>"
  #   ),
  #   opacity = 1
  # ) %>%
  leaflet.extras::addResetMapButton()

household_communities_leaft
```


```{r eval=FALSE}
htmlwidgets::saveWidget(
  household_communities_leaft,
  file = "./02_output/plots/household_communities_leaft.html"
)

webshot::webshot(
  "./02_output/plots/household_communities_leaft.html",
  "./02_output/plots/household_communities_leaft.png",
  cliprect = "viewport",
  zoom = 4
)
```