## Characterizing semi-nomadic populations
## Functions behind figures and tables for manuscript
## Date: February 3, 2024
## Author: Hannah Meredith

library(dplyr)
library(tidyr)
library(ggplot2)
library(epiDisplay)
library(lme4)
library(extrafont)
library(leaflet)
library(leaflet.esri)
library(geosphere)


make.pop.tracks.figure <- function(df){
  Facility <- data.frame(name = c("Kangirisae", "Kerio", "Louwae", "Nakurio", "Lodwar"),
                         abrev = c("K", "R", "L", "N", "Lod"),
                         Latitude = c(2.613988, 2.998699, 2.825620, 2.875389, 3.120799),
                         Longitude = c(36.253085,36.050420, 36.203688, 36.161556 , 35.602175))
  
  
  lat = median(df$Hotspot_lat)
  lon = median(df$Hotspot_lon)
  
  #### plot of all coordinates by catchment area
  FigXA <- leaflet(options = leafletOptions(zoomSnap = 0.25, zoomDelta=0.25)) %>%  # facilitates fractional zooms
    setView(lng =lon,  lat = lat, zoom = 10) %>%
    addTiles() %>%
    addProviderTiles("Esri.WorldImagery", group = "Satellite")%>%
    addCircleMarkers(lng= ~ Hotspot_lon, lat= ~ Hotspot_lat, data = df,
                     color= df$facility_color, weight = 1, radius = 5) %>%
    addCircleMarkers(lng= ~ Longitude, lat= ~ Latitude, data = Facility,
                     color= "black", fillOpacity = 0, weight = 1, radius = 5)%>%
    addScaleBar(position = c("topright"), scaleBarOptions(imperial = FALSE))
  
  #### plot of all coordinates by number of unique households logged
  hotspot.hh.counts <- df %>% 
    group_by(HotspotID, Hotspot_lon, Hotspot_lat)%>%
    summarise(hotspot.hh.num = length(unique(GPS_carrier_ID)))
  hotspot.hh.counts$hotspot.hh.bin <- cut(hotspot.hh.counts$hotspot.hh.num,
                                          breaks = c(-Inf, 1, 2,  5, 20),
                                          labels = c("1", "2", "3-5", ">5")) 
  
  hotspot.hh.counts <- hotspot.hh.counts %>% 
    mutate(color = case_when(stringr::str_detect(hotspot.hh.bin, "1") ~ "#440154FF",
                             stringr::str_detect(hotspot.hh.bin, "2") ~ "#31688EFF",
                             stringr::str_detect(hotspot.hh.bin, "3-5") ~ "#35B779FF",
                             stringr::str_detect(hotspot.hh.bin, ">5") ~ "#FDE725FF",
                             TRUE ~ "a default"))
  
  hotspot.hh.counts <- hotspot.hh.counts[order(hotspot.hh.counts$hotspot.hh.bin),]
  
  FigXB <- leaflet(options = leafletOptions(zoomSnap = 0.25, zoomDelta=0.25)) %>%  # facilitates fractional zooms
    setView(lng =lon,  lat = lat, zoom = 9.75) %>%
    addTiles() %>%
    addProviderTiles("Esri.WorldImagery", group = "Satellite")%>%
    addCircleMarkers(lng= ~ Hotspot_lon, lat= ~ Hotspot_lat, data = hotspot.hh.counts, 
                     color= hotspot.hh.counts$color, fillOpacity = 1, weight = 1, radius = 5) %>%
    addCircleMarkers(lng= ~ Longitude, lat= ~ Latitude, data = Facility, color= "black", weight = 1, radius = 5)%>%
    addScaleBar(position = c("topright"), scaleBarOptions(imperial = FALSE))
  
  output = list(FigXA, FigXB)
  return(output)
  
}

characterize.trip.types <- function(all.tracks, night.tracks){
  all.change.status <- all.tracks %>% group_by(GPS_carrier_ID) %>%   # considering all locations visited, did they tend to stay at one for majority of time? 
    summarise(max.prop = max(Proportion_locations),
              all.change.status = ifelse(max.prop >= 0.90, "All.Stay", "All.Change"))
  night.change.status <- night.tracks %>% group_by(GPS_carrier_ID) %>% # considering night locations visited, did they tend to stay at one for majority of time? 
    summarise(max.night = max(Proportion_locations),
              night.all.change.status = ifelse(max.night >= 0.9, "Night.Stay", "Night.Change"))
  change.status <- left_join(all.change.status, night.change.status, by = "GPS_carrier_ID")
  
  ## Of those who changed day and night - what were the movements like? Long term moves (i.e. >= 7 consec nights spent) or more transient stays? 
  distinct.nights <- night.tracks %>% group_by(GPS_carrier_ID, ClusterID, Date_adj) %>% arrange(Date_time_adj) %>% slice(n()) %>% 
    distinct(GPS_carrier_ID, ClusterID, Date_adj, Trip_length, Total_distance_km, Total_night_distance_km) %>%ungroup() # One date per ID (as opposed to multiple Date-times) - take last point recorded of each night
  distinct.nights <- distinct.nights %>%
    group_by(GPS_carrier_ID)%>%
    arrange(Date_adj) %>%
    mutate(idx = 1:n())
  
  distinct.nights <- distinct.nights %>%              # Assign cluster ID for consecutive nights      
    group_by(GPS_carrier_ID, ClusterID)%>%
    mutate(group = c(0, cumsum(diff(idx) > 1)) + 1)
  
  distinct.nights <- distinct.nights %>%
    group_by(GPS_carrier_ID, ClusterID, group)%>%
    mutate(consec.nights = n())
  
  distinct.nights <- distinct.nights %>% mutate(trip.type = ifelse(consec.nights < 7, "Transient", "LongerTerm"))
  distinct.nights <- distinct.nights %>%
    group_by(GPS_carrier_ID) %>%
    mutate(nights = n(), 
           longterm = length(which(trip.type == "LongerTerm")),
           transient = length(which(trip.type == "Transient")),
           prop.longterm = longterm/nights, 
           prop.transient = transient/nights
    )
  
  distinct.nights.summary <- distinct.nights %>% ungroup() %>% distinct(GPS_carrier_ID, Trip_length, Total_distance_km, Total_night_distance_km, 
                                                                        ClusterID, consec.nights, group, trip.type, 
                                                                        prop.longterm, prop.transient)
  
  distinct.nights.summary <- distinct.nights.summary %>%
    group_by(GPS_carrier_ID, Trip_length, Total_distance_km, Total_night_distance_km, prop.longterm, prop.transient) %>%
    summarise(num.changes = n()-1,
              count.longterm = table(trip.type)[1],
              count.transient = table(trip.type)[2])  
  
  num.unique.camps <- distinct.nights %>%
    group_by(GPS_carrier_ID)%>%
    summarise(num.campsites = length(unique(ClusterID)))
  
  gps.point.stats <- left_join(distinct.nights.summary, change.status, by = "GPS_carrier_ID")
  gps.point.stats$combo.change.stat <- paste(gps.point.stats$all.change.status, gps.point.stats$night.all.change.status, sep = ".")
  gps.point.stats <- left_join(gps.point.stats, num.unique.camps, by = "GPS_carrier_ID")
  gps.point.stats$trip.type <- ifelse(gps.point.stats$combo.change.stat == "All.Change.Night.Stay", "Day_trip",
                                      ifelse(gps.point.stats$combo.change.stat == "All.Stay.Night.Stay", "Static",
                                             ifelse(gps.point.stats$prop.longterm <=0.5 & gps.point.stats$night.all.change.status != "Night.Stay", "Transient", 
                                                    ifelse(gps.point.stats$prop.longterm >0.5 & gps.point.stats$night.all.change.status != "Night.Stay", "Longterm", NA))))

  
  trip.type <- gps.point.stats %>% ungroup() %>%  distinct(GPS_carrier_ID, trip.type)
  
  all.tracks <- left_join(all.tracks, trip.type, by = "GPS_carrier_ID")
  
  output = list(all.tracks, gps.point.stats)
  return(output)
}

household.epi <- function(df, time.pt){
  
  household.epi_num <- df %>%
    summarise(People_per_HH = median(Members_count, na.rm = TRUE),
              Traveler_count = median(Traveler_count, na.rm = TRUE),
              Remainer_count = median(Remainer_count, na.rm = TRUE),
              Any_children = sum(Children_count > 0, na.rm = TRUE),
              Nets_per_HH = median(Bednet_count, na.rm = TRUE),
              Kangirisae = sum(Facility == "Kangirisae", na.rm = TRUE),
              Kerio = sum(Facility == "Kerio", na.rm = TRUE),
              Louwae = sum(Facility == "Louwae", na.rm = TRUE),
              Nakurio = sum(Facility == "Nakurio", na.rm = TRUE),
              Water_open = sum(Household_water %in% c("Open only", "Open or closed"), na.rm = T),
              Water_closed = sum(Household_water %in% c("Closed only", "Open or closed"), na.rm = T),
              Income_livestock = sum(Income_livestock == 1, na.rm = TRUE),
              Income_charcoal = sum(Income_charcoal == 1, na.rm = TRUE), 
              Income_weaving = sum(Income_weaving == 1, na.rm = TRUE), 
              Income_farming = sum(Income_farming == 1, na.rm = TRUE), 
              Income_informallyEmployed = sum(Income_informallyEmployed, na.rm = TRUE),
              Own_goats = length(which(Own_goats>0)),
              Own_sheep = length(which(Own_sheep>0)),
              Own_camels = length(which(Own_camels>0)),
              Own_cattle = length(which(Own_cattle>0)),
              
    )
  
  
  household.epi_denom <- df%>%
    summarise(People_per_HH = NA,
              Traveler_count = NA,
              Remainer_count = NA,
              Any_children = length(which(!is.na(Children_count))),
              Nets_per_HH = NA,
              Kangirisae = length(which(!is.na(Facility))),
              Kerio = length(which(!is.na(Facility))),
              Louwae = length(which(!is.na(Facility))),
              Nakurio = length(which(!is.na(Facility))),
              Water_open = length(which(!is.na(Household_water))),
              Water_closed = length(which(!is.na(Household_water))),
              Income_livestock = length(which(!is.na(Income_livestock))), 
              Income_charcoal = length(which(!is.na(Income_charcoal))), 
              Income_weaving = length(which(!is.na(Income_weaving))), 
              Income_farming = length(which(!is.na(Income_farming))), 
              Income_informallyEmployed = length(which(!is.na(Income_informallyEmployed))),
              Own_goats = length(which(!is.na(Own_goats))), 
              Own_sheep = length(which(!is.na(Own_sheep))), 
              Own_camels = length(which(!is.na(Own_camels))), 
              Own_cattle = length(which(!is.na(Own_cattle)))
    )
  
  compiled.household <- rbind.data.frame(household.epi_num, household.epi_denom)
  compiled.household <- t(compiled.household)
  colnames(compiled.household) = c("Numerator", "Denominator")
  
  compiled.household <- as.data.frame(compiled.household)
  compiled.household$percentage <- round(compiled.household$Numerator / compiled.household$Denominator *100, 1)
  compiled.household$percentage[rownames(compiled.household) == "People_per_HH"] <- paste(quantile(df$Household_size, 1/4, na.rm = TRUE),
                                                                                         quantile(df$Household_size, 3/4, na.rm = TRUE), sep = "-")
  compiled.household$percentage[rownames(compiled.household) == "Traveler_count"] <- paste(quantile(df$Traveler_count, 1/4, na.rm = TRUE),
                                                                                        quantile(df$Traveler_count, 3/4, na.rm = TRUE), sep = "-")
  compiled.household$percentage[rownames(compiled.household) == "Remainer_count"] <- paste(quantile(df$Remainer_count, 1/4, na.rm = TRUE),
                                                                                        quantile(df$Remainer_count, 3/4, na.rm = TRUE), sep = "-")
  compiled.household$percentage[rownames(compiled.household) == "Nets_per_HH"] <- paste(quantile(df$Bednet_count, 1/4, na.rm = TRUE),
                                                                                        quantile(df$Bednet_count, 3/4, na.rm = TRUE), sep = "-")
 
  compiled.household$household <- paste(compiled.household$percentage, " (", compiled.household$Numerator, "/", compiled.household$Denominator ,")", sep = "")
  
  return(compiled.household)
}

member.epi.follow <- function(df){
  
  member.epi_num <- df %>%
    summarise(female = sum(Gender == 0, na.rm = TRUE),   
              male = sum(Gender == 1, na.rm = TRUE), 
              age.under.15 = sum(Age <= 15, na.rm = TRUE),
              age.16.40 = sum(Age > 15 & Age <= 40, na.rm = TRUE),
              age.40.plus = sum(Age > 40, na.rm = TRUE),
              Kangirisae = sum(Facility == "Kangirisae", na.rm = TRUE),
              Kerio = sum(Facility == "Kerio", na.rm = TRUE),
              Louwae = sum(Facility == "Louwae", na.rm = TRUE),
              Nakurio = sum(Facility == "Nakurio", na.rm = TRUE),
              Water_open = sum(Household_water %in% c("Open only", "Open or closed"), na.rm = T),
              Water_closed = sum(Household_water %in% c("Closed only", "Open or closed"), na.rm = T),
              No.symptoms = sum(Followup_noSymptoms == 1, na.rm = TRUE), 
              recent.sick = sum(Sick_since_enrollment == 1 , na.rm = TRUE), 
              hospital = sum(Sick_action_HealthFacility == 1, na.rm = TRUE),
              self.med = sum(Sick_action_SelfMedicate == 1, na.rm = TRUE),
              mal.test = sum(Sick_action_MalariaTest == 1, na.rm = T),
              mal.test.pos = sum(Sick_malariaTest_Results == 1, na.rm = T),
              take.meds = sum(Sick_action_takeMedicine == 1, na.rm = T),
              antimals = sum(Medicine_antimalarial == 1, na.rm = T),
              antibiotics = sum(Medicine_antibiotic == 1, na.rm = T),
              pain.killers = sum(Medicine_painKiller == 1, na.rm = T)
    )
  
  
  member.epi_denom <- df%>%
    summarise(female = length(which(!is.na(Gender))),
              male = length(which(!is.na(Gender))), 
              age.under.15 = length(which(!is.na(Age))), 
              age.16.40 = length(which(!is.na(Age))), 
              age.40.plus = length(which(!is.na(Age))), 
              Kangirisae = length(which(!is.na(Facility))),
              Kerio = length(which(!is.na(Facility))),
              Louwae = length(which(!is.na(Facility))),
              Nakurio = length(which(!is.na(Facility))),
              Water_open = length(which(!is.na(Water_open))),
              Water_closed = length(which(!is.na(Water_closed))),
              No.symptoms = length(which(!is.na(Followup_noSymptoms))), 
              recent.sick = length(which(!is.na(Sick_since_enrollment))), 
              hospital = length(which(Sick_since_enrollment==1)),
              self.med = length(which(Sick_since_enrollment == 1)),
              mal.test = length(which(Sick_since_enrollment==1)),
              mal.test.pos = length(which(!is.na(Sick_malariaTest_Results))),
              take.meds = length(which(Sick_since_enrollment==1)),
              antimals = length(which(Sick_since_enrollment==1)),
              antibiotics = length(which(Sick_since_enrollment==1)),
              pain.killers = length(which(Sick_since_enrollment==1))
    )
  
  compiled.member <- rbind.data.frame(member.epi_num, member.epi_denom)
  compiled.member <- t(compiled.member)
  colnames(compiled.member) = c("Numerator", "Denominator")
  compiled.member <- as.data.frame(compiled.member)
  compiled.member$percentage <- round(compiled.member$Numerator / compiled.member$Denominator *100, 1)
  compiled.member$member <- paste(compiled.member$percentage, " (", compiled.member$Numerator, "/", compiled.member$Denominator ,")", sep = "")
  
  return(compiled.member)
}

member.epi.follow.diff.denom <- function(df){
  
  member.epi_num <- df %>%
    summarise(female = sum(Gender == 0 & PCR_result_Falciparum == 1, na.rm = TRUE), 
              male = sum(Gender == 1 & PCR_result_Falciparum == 1, na.rm = TRUE),
              age.under.15 = sum(Age <= 15 & PCR_result_Falciparum == 1, na.rm = TRUE),
              age.16.40 = sum(Age > 15 & Age <= 40 & PCR_result_Falciparum == 1, na.rm = TRUE),
              age.40.plus = sum(Age > 40 & PCR_result_Falciparum == 1, na.rm = TRUE),
              Kangirisae = sum(Facility == "Kangirisae" & PCR_result_Falciparum == 1, na.rm = TRUE),
              Kerio = sum(Facility == "Kerio" & PCR_result_Falciparum == 1, na.rm = TRUE),
              Louwae = sum(Facility == "Louwae" & PCR_result_Falciparum == 1, na.rm = TRUE),
              Nakurio = sum(Facility == "Nakurio" & PCR_result_Falciparum == 1, na.rm = TRUE),
              Water_open = sum(Water_open == 1 & PCR_result_Falciparum == 1, na.rm = TRUE),
              Water_closed = sum(Water_closed == 1 & PCR_result_Falciparum == 1, na.rm = TRUE),
              No.symptoms = sum(Followup_noSymptoms == 1 & PCR_result_Falciparum == 1, na.rm = TRUE), 
              recent.sick = sum(Sick_since_enrollment == 1  & PCR_result_Falciparum == 1, na.rm = TRUE), 
              hospital = sum(Sick_action_HealthFacility == 1 & PCR_result_Falciparum == 1, na.rm = TRUE),
              self.med = sum(Sick_action_SelfMedicate == 1 & PCR_result_Falciparum == 1, na.rm = TRUE),
              mal.test = sum(Sick_action_MalariaTest == 1 & PCR_result_Falciparum == 1, na.rm = T),
              mal.test.pos = sum(Sick_malariaTest_Results == 1 & PCR_result_Falciparum == 1, na.rm = T),
              take.meds = sum(Sick_action_takeMedicine == 1 & PCR_result_Falciparum == 1, na.rm = T),
              antimals = sum(Medicine_antimalarial == 1 & PCR_result_Falciparum == 1, na.rm = T),
              antibiotics = sum(Medicine_antibiotic == 1 & PCR_result_Falciparum == 1, na.rm = T),
              pain.killers = sum(Medicine_painKiller == 1 & PCR_result_Falciparum == 1, na.rm = T)
    )
  
  
  member.epi_denom <- df %>%
    summarise(female = sum(Gender == 0, na.rm = TRUE), 
              male = sum(Gender == 1, na.rm = TRUE), 
              age.under.15 = sum(Age <= 15, na.rm = TRUE),
              age.16.40 = sum(Age > 15 & Age <= 40, na.rm = TRUE),
              age.40.plus = sum(Age > 40, na.rm = TRUE),
              Kangirisae = sum(Facility == "Kangirisae", na.rm = TRUE),
              Kerio = sum(Facility == "Kerio", na.rm = TRUE),
              Louwae = sum(Facility == "Louwae", na.rm = TRUE),
              Nakurio = sum(Facility == "Nakurio", na.rm = TRUE),
              Water_open = sum(Water_open == 1, na.rm = TRUE),
              Water_closed = sum(Water_closed == 1, na.rm = TRUE),
              No.symptoms = sum(Followup_noSymptoms == 1, na.rm = TRUE), 
              recent.sick = sum(Sick_since_enrollment == 1 , na.rm = TRUE), 
              hospital = sum(Sick_action_HealthFacility == 1, na.rm = TRUE),
              self.med = sum(Sick_action_SelfMedicate == 1, na.rm = TRUE),
              mal.test = sum(Sick_action_MalariaTest == 1, na.rm = T),
              mal.test.pos = sum(Sick_malariaTest_Results == 1, na.rm = T),
              take.meds = sum(Sick_action_takeMedicine == 1, na.rm = T),
              antimals = sum(Medicine_antimalarial == 1, na.rm = T),
              antibiotics = sum(Medicine_antibiotic == 1, na.rm = T),
              pain.killers = sum(Medicine_painKiller == 1, na.rm = T)
              
    )
  
  compiled.member <- rbind.data.frame(member.epi_num, member.epi_denom)
  compiled.member <- t(compiled.member)
  colnames(compiled.member) = c("Numerator", "Denominator")
  compiled.member <- as.data.frame(compiled.member)
  compiled.member$percentage <- round(compiled.member$Numerator / compiled.member$Denominator *100, 1)
  compiled.member$member <- paste(compiled.member$percentage, " (", compiled.member$Numerator, "/", compiled.member$Denominator ,")", sep = "")
  
  return(compiled.member)
}