---
title: "Arboviral Vectors at Zoonotic Interfaces Rmd Script for P&V publication"
author: "Rebecca Brown"
date: "2025-10-31"
output: html_document
---

###### This rmd script analyses the mosquito data from the 2024 PANDASIA dry and wet season invertebrate sampling to publish on arboviral vectors in zoonotic interfaces in Thailand. Mosquitoes were collected from 7 interfaces hypothesised to be high risk for zoonotic pathogen spillover due to the presence of animals harbouring potential zoonotic viruses, and direct overlap or sharing the environment with humans. These were Houses (HO), Temples (TE), Dumpsites (DU), Orchards (OR), Forest edge (FE), Forest interior (FI) and Bat Caves (BC). Four districts were selected for the location of the study; Pong Nam Ron/Soi Dao in Chanthaburi, Wiang Kaen and Mae Fah Luang in Chiang Rai.

###### Four trapping techniques were used to sample the mosquito population: PK - prokopack aspiration, a timed 10 minute search of shelters or vegetation in the morning, PRO - Biogents PRO Traps, baited with BG lures mimicking human scent and dry ice releasing CO2 to simulate mammalian hosts, LT - Light Traps, baited with dry ice and light to attract host seeking mosquitoes, GAT - Gravid Aedes traps, black buckets baited with hay infused water to lure ovipositing mosquitoes. These four trapping techniques were clustered into groups, where 2 to 3 groups of traps were positioned at each interface with a minimum distance of 100 metres between groups, 20m between traps within a group. Two repeats of collections were performed on subsequent days for each location in both wet and dry seasons.

```{r setup, include=FALSE}
# knitr::opts_chunk$set(echo = FALSE, message = FALSE,
#                       root.dir = 'C:/Users/PUT-FILE-LOCATION-HERE/')
#install packages
# install.packages(c("knitr", "readxl", "kableExtra", "lubridate", "tidyverse","ggplot2", "multcomp", "cowplot", "openxlsx", "reshape", "glmmTMB", "emmeans", "ggeffects", "DHARMa", "wesanderson", "lsmeans", "scales", "viridis", ggpubr))


#load packages
library(knitr)
library(readxl)
library(kableExtra)
library(lubridate)
library(tidyverse)
library(ggplot2)
library(multcomp)
library(cowplot)
library(openxlsx)
library(reshape)
library(glmmTMB)
library(emmeans)
library(ggeffects)
library(DHARMa)
library(wesanderson)
library(lsmeans)
library(scales)
library(ggpubr)
library(ggsignif)
library(iNEXT)
library(patchwork)
library(lemon)

# set theme
my_theme <- theme(strip.background = element_rect(fill = "white", colour = "black"),
        panel.background = element_rect(fill = "aliceblue", color = "black"),
        panel.grid.major = element_line(color = "gray", linetype = "dotted"),
        panel.grid.minor = element_line(color = "gray", linetype = "dotted"),
        axis.text = element_text(size=12, color="black"),
        axis.title = element_text(size=14, color="black"),
        legend.title = element_text(size = 14),
        legend.text=element_text(size=12),
        plot.title = element_text(size = 16),
        plot.caption = element_text(size=12, hjust = 0),
        strip.text.x = element_text(size=14, color="black",face="bold"),
        strip.text.y = element_text(size=14, color="black",face="bold"))

```

```{r import data}

#import cleaned mosquito data

mq_data <- read.xlsx('C:/Users/XXXX.xlsx')
```

```{r other data cleaning}
#other formatting to do before publishing

mq_data <- mq_data %>% 
  mutate(closest_id = ifelse(closest_id == "Culex   spp.", "Culex  spp.",
                             ifelse(closest_id == "Coquillettidia  spp.", "Coquillettidia (Coq.) spp.",
                                    closest_id))) %>%
  mutate(district = ifelse(district == "Pong Nam Ron", "Pong Nam Ron & Soi Dao", district))
```

```{r, overall counts}
#data restructuring for Table 2

#group both seasons for 2024
#converting males to 0 for the summary for ease in reporting analyses later
mq_females <- mq_data %>% mutate(females = ifelse(sex == "F", count, 0))
mq_males <- mq_data %>% mutate(males = ifelse(sex == "M", count, 0))

#sum females and sum males each season
f_sum_d <- mq_females %>% filter(season == "dry") %>% summarise(sum(females))
m_sum_d <- mq_males %>% filter(season == "dry") %>% summarise(sum(males))
f_sum_w <- mq_females %>% filter(season == "wet") %>% summarise(sum(females))
m_sum_w <- mq_males %>% filter(season == "wet") %>% summarise(sum(males))

#make summary of replicates
replicates_2024 <- mq_data %>%
  dplyr::select(start_date, start_time, season, district, interface, trap, group) %>%
  distinct() %>%
  group_by(district, interface, trap) %>% #can group on season here if want to break into wet/dry
  mutate(replicates = n()) %>%
  dplyr::select(!c(start_date, start_time, season, group)) %>%
  distinct() #%>%
  #mutate(repeats = paste0("(", repeats, ")"))

#summary of female abundances  
f_summary_2024 <- mq_females %>%
  group_by(district, interface, trap) %>%
  summarise(total = sum(females)) %>% 
  mutate(total = ifelse(is.na(total), 0, total)) %>%
  # mutate(sex =c("F +")) %>%
  # unite(females, c(total, sex), sep = " ", remove = TRUE)
  dplyr::rename(females = total)

#summary of male abundances  
m_summary_2024 <- mq_males %>%
  group_by(district, interface, trap) %>%
  summarise(total = sum(males)) %>% 
  mutate(total = ifelse(is.na(total), 0, total)) %>%
  # mutate(sex =c("M")) %>%
  # unite(males, c(total, sex), sep = " ", remove = TRUE)
  dplyr::rename(males = total)

#create overall sum column

f_temp <- mq_females %>%
  group_by(district, interface) %>%
  summarise(f_total = sum(females)) %>% 
  mutate(f_total = ifelse(is.na(f_total), 0, f_total))

m_temp <- mq_males %>%
  group_by(district, interface) %>%
  summarise(m_total = sum(males)) %>% 
  mutate(m_total = ifelse(is.na(m_total), 0, m_total)) %>%
  merge(f_temp) %>% 
  mutate(Total = f_total + m_total) %>% 
  dplyr::select(district, interface, Total)

#bring both sexes together with replicates
full_summary_2024 <- f_summary_2024 %>%
  merge(m_summary_2024) %>%
  merge(replicates_2024) %>%
  group_by(district, interface, females, males, replicates) %>%
  pivot_wider(names_from = c(trap),
              values_from = c(females, males, replicates)) %>%
  merge(m_temp) %>%
  mutate(district = if_else(district == "Chanthaburi", "Pong Nam Ron & Soi Dao", district))

col_order <- c("district", "interface", "females_pro", "males_pro", "replicates_pro", 
               "females_lt",  "males_lt", "replicates_lt", "females_gat", "males_gat",
               "replicates_gat", "females_pkin" , "males_pkin", "replicates_pkin",
               "females_pkout", "males_pkout", "replicates_pkout", "Total")

full_summary_2024 <- full_summary_2024[, col_order]

#set levels
full_summary_2024$interface <- factor(full_summary_2024$interface, levels = c("House", "Temple", "Dumpsite", "Orchard", "Forest edge",  "Forest interior", "Bat cave"))
full_summary_2024$district <- factor(full_summary_2024$district, levels = c("Pong Nam Ron & Soi Dao", "Mae Fah Luang", "Wiang Kaen")) 
full_summary_2024 <- full_summary_2024[order(full_summary_2024$district, full_summary_2024$interface),]
  
#export overall counts and replicates summary
write.xlsx(full_summary_2024, file = "Datafiles\\Table_2_Mosquito_counts_at_interfaces.xlsx", sheetName = 'PANDASIA', append = FALSE)

rm(replicates_2024, f_sum_d, f_sum_w, f_temp, m_sum_d, m_sum_w, m_summary_2024, m_temp, mq_males)
rm(f_summary_2024)

```

```{r, f genus counts}
#remove whitespace
mq_data$genus <- gsub(" ", "", mq_data$genus, fixed = TRUE)

#select females only

genus_counts_f <- mq_data %>%
  filter(sex == "F") %>%
  group_by(district, interface, genus) %>%
  summarise(total = sum(count)) %>%
  filter(!genus %in% c("Unknown", "empty")) %>%
  pivot_wider(names_from = c(district,interface),
              values_from = c(total)) %>%
  replace(is.na(.), 0)

#set levels
genus_counts_f <- genus_counts_f[,c("genus", "Pong Nam Ron & Soi Dao_House", "Pong Nam Ron & Soi Dao_Temple", "Pong Nam Ron & Soi Dao_Dumpsite", "Pong Nam Ron & Soi Dao_Orchard", "Pong Nam Ron & Soi Dao_Forest edge", 
                                    "Pong Nam Ron & Soi Dao_Forest interior", "Pong Nam Ron & Soi Dao_Bat cave", "Mae Fah Luang_House", "Mae Fah Luang_Dumpsite", "Mae Fah Luang_Orchard",
                                    "Mae Fah Luang_Forest edge", "Mae Fah Luang_Forest interior", "Wiang Kaen_House", "Wiang Kaen_Orchard", "Wiang Kaen_Forest edge",
                                    "Wiang Kaen_Forest interior", "Wiang Kaen_Bat cave")]

#create totals column
genus_counts_f <- genus_counts_f %>% ungroup() %>%
  mutate(Total = rowSums(.[2:18])) %>%
  #arrange(desc(Total)) %>%
  bind_rows(summarise(., across(where(is.numeric), sum),
                         across(where(is.character), ~'Total')))

# export overall counts by genus
write.xlsx(genus_counts_f, file = "Datafiles\\Table_3.1_Female_mosquito_genera_at_interfaces.xlsx", sheetName = 'PANDASIA', append = FALSE)

```

```{r, f species counts}

species_counts_f <- mq_data %>% 
  filter(sex == "F") %>%
  group_by(district, interface, closest_id) %>%
  summarise(total = sum(count)) %>%
  filter(!closest_id %in% c("empty", "Unknown  Unknown", "Unknown Unknown Unknown")) %>%
  pivot_wider(names_from = c(district, interface),
              values_from = c(total)) %>%
  as.data.frame() %>%
  mutate_at(c(2:18), ~replace(., is.na(.), 0)) %>%
  mutate(Total = rowSums(.[2:18])) %>%
  bind_rows(summarise(., across(where(is.numeric), sum),
                         across(where(is.character), ~'Total')))
 
#set levels
species_counts_f <- species_counts_f[,c("closest_id", "Pong Nam Ron & Soi Dao_House", "Pong Nam Ron & Soi Dao_Temple",
                                        "Pong Nam Ron & Soi Dao_Dumpsite", "Pong Nam Ron & Soi Dao_Orchard", "Pong Nam Ron & Soi Dao_Forest edge", 
                                        "Pong Nam Ron & Soi Dao_Forest interior", "Pong Nam Ron & Soi Dao_Bat cave", "Mae Fah Luang_House", 
                                        "Mae Fah Luang_Dumpsite", "Mae Fah Luang_Orchard", "Mae Fah Luang_Forest edge", "Mae Fah Luang_Forest interior", 
                                        "Wiang Kaen_House", "Wiang Kaen_Orchard", "Wiang Kaen_Forest edge", "Wiang Kaen_Forest interior", "Wiang Kaen_Bat cave", "Total")]

#export overall counts by species
write.xlsx(species_counts_f, file = "Datafiles\\Female_mosquito_species_at_interfaces.xlsx", sheetName = 'PANDASIA', append = FALSE)

#combining both species and genus

full_id_counts_f <- genus_counts_f %>% 
  dplyr::rename(closest_id = genus) %>% filter(!closest_id == "Total")

#full_id_counts$closest_id <- paste0(full_id_counts$closest_id, " genus")
full_id_counts_f <- rbind(full_id_counts_f, species_counts_f)
full_id_counts_f <- arrange(full_id_counts_f, closest_id)
write.xlsx(full_id_counts_f, file = "Datafiles\\Mosquito_female_species_genus_at_interfaces.xlsx", sheetName = 'PANDASIA', append = FALSE)

```

```{r, f vector counts}
#select important arboviral vectors (females) to scan their distribution across study interfaces

vector_species_f <- species_counts_f %>%
  filter(closest_id %in% c("Aedes (Ste.) aegypti", "Aedes (Ste.) albopictus", "Culex (Cux.) fuscocephala",
                           "Culex (Cux.) gelidus", "Culex  (Cux.) pseudovishnui", "Culex (Cux.) tritaeniorhynchus", 
                           "Culex (Cux.) vishnui", "Culex (Cux.) quinquefasciatus")) %>%
   mutate(Total = rowSums(.[2:18])) %>%
   bind_rows(summarise(., across(where(is.numeric), sum),
                         across(where(is.character), ~'Total')))

write.xlsx(vector_species_f, file = "Datafiles\\Female_vectors_at_interfaces.xlsx", sheetName = 'PANDASIA', append = FALSE)
```
  
```{r bloodfed females}

feeding_counts <- mq_data %>%
  filter(sex == "F") %>%
  group_by(closest_id, interface, feeding_status) %>%
  summarise(total = sum(count)) %>% 
  filter(feeding_status %in% c("BF", "SG", "BF & SG")) %>%
  pivot_wider(names_from = c(interface),
              values_from = c(total)) %>%
  mutate_at(c(3:8), ~replace(., is.na(.), 0)) %>%
  ungroup() %>%
  mutate(Total = rowSums(.[3:8])) %>% arrange(desc(Total)) %>%
    bind_rows(summarise(., across(where(is.numeric), sum),
                         across(where(is.character), ~'Total'))) 

#set levels
feeding_counts <- feeding_counts[,c("closest_id", "feeding_status", "House", "Temple", "Dumpsite", "Orchard", "Forest edge", "Forest interior", "Total")]

kable(feeding_counts, align = 'l', col.names = c("Species", "Feeding status", "House", "Temple", "Dumpsite", "Orchard", "Forest edge", "Forest interior", "Total"), caption = "Blood-fed and semi-gravid female mosquitoes caught in study interfaces") %>%
kable_styling(bootstrap_options = "condensed", font_size = 11)

write.xlsx(feeding_counts, file = "Datafiles\\Bloodfed_species.xlsx", sheetName = 'PANDASIA', append = FALSE)
```

```{r, hs abundances}

#abundance of females of different behaviours at interfaces
#includes empty collections and collections with males and no females
#group data by catch

catch_females <- mq_females %>%
  group_by(start_date, start_time, district, interface, group, trap, season, site) %>%
  summarise(count_f = sum(females)) %>% ungroup()

#set categorical variables as factors
catch_females$district <-as.factor(catch_females$district)
catch_females$interface <-as.factor(catch_females$interface)
catch_females$trap <-as.factor(catch_females$trap)
catch_females$season <-as.factor(catch_females$season)
catch_females$start_date <-as.Date(catch_females$start_date)

hs_catch_f <- catch_females %>% filter(trap %in% c("pro", "lt"))
hs_catch_f$site <- as.factor(hs_catch_f$site)

hs_1 <- glmmTMB(count_f ~ interface + season + district + trap + (1|site),
            data = hs_catch_f, family = nbinom2)
summary(hs_1) #AIC 2154.1

hs_2 <- glmmTMB(count_f ~ interface + district + trap + (1|site),
           data = hs_catch_f, family = nbinom2)
summary(hs_2) #AIC 2196.6

hs_3 <- glmmTMB(count_f ~ season + district + trap + (1|site),
            data = hs_catch_f, family = nbinom2)
summary(hs_3) #aic 2168.0

hs_4 <- glmmTMB(count_f ~ interface + season + trap + (1|site),
            data = hs_catch_f, family = nbinom2)
summary(hs_4) #AIC 2160.6

hs_5 <- glmmTMB(count_f ~ interface + district + season + (1|site),
            data = hs_catch_f, family = nbinom2)
summary(hs_5) #AIC 2154.2

hs_6 <- glmmTMB(count_f ~ interface * trap + district + season + (1|site),
            data = hs_catch_f, family = nbinom2)
summary(hs_6) #AIC 2154.5

anova(hs_1, hs_2) #season sig. (Chisq = 44.5, Df = 1, p < 0.001)
anova(hs_1, hs_3) #interface sig.(Chisq = 26.0, Df = 6, p < 0.001)
anova(hs_1, hs_4) #district sig. (Chisq = 10.5, Df = 2, p < 0.01)
anova(hs_1, hs_5) #trap ns (Chisq = 2.1, Df = 1, p > 0.05)
anova(hs_1, hs_6) #interaction trap v interface (Chisq = 11.6, Df = 6, p > 0.05)

#dharma package, model diagnostics
simulationOutput <- simulateResiduals(fittedModel = hs_1)
plot(simulationOutput)
testZeroInflation(simulationOutput) #ratio greater than 1 means more zeros than expected but when p is not significant not critical to need a zeroinfl model
testUniformity(simulationOutput)
testDispersion(simulationOutput, alternative = "greater") #tests for only overdispersion
testDispersion(simulationOutput)

summary(glht(hs_1, linfct=mcp(season='Tukey')))
summary(glht(hs_1, linfct=mcp(interface='Tukey')))
summary(glht(hs_1, linfct=mcp(district='Tukey')))

#make predictions from best model
hsf_ab_pred <- ggpredict(hs_1, c("interface", "season", "district"))
hsf_ab_pred <- hsf_ab_pred %>%
  dplyr::rename(interface = x, mean = predicted, season = group, district = facet) %>%
  as.data.frame()

#remove predictions that dont make sense
hsf_ab_pred <- hsf_ab_pred %>%
  filter(!(interface == "Dumpsite" & district == "Wiang Kaen")) %>%
  filter(!(interface == "Temple" & district == "Wiang Kaen")) %>%
  filter(!(interface == "Temple" & district == "Mae Fah Luang")) %>%
  filter(!(interface == "Bat cave" & district == "Mae Fah Luang"))

hsf_ab_pred$interface <- factor(hsf_ab_pred$interface, levels = c("House", "Temple", "Dumpsite", "Orchard", "Forest edge",  "Forest interior", "Bat cave"))
hsf_ab_pred$district <- factor(hsf_ab_pred$district, levels = c("Pong Nam Ron & Soi Dao", "Mae Fah Luang", "Wiang Kaen"))

f_hs_plot <- ggplot(hsf_ab_pred, aes(fill=interface, y=mean, x=season)) + 
    geom_bar(position="dodge", color = "grey", stat="identity") + 
    geom_errorbar(hsf_ab_pred, mapping = aes(ymin=conf.low, ymax=conf.high), linewidth=0.1, position = "dodge") +
    scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(0, 1, 2, 5, 10, 20, 50, 100, 250, 500, 1000)) +
    scale_fill_brewer(palette="PRGn", name="Interface") +
    labs(fill = "Season", x = "Season", y = "Mean per collection (95% CI)", title = "") +
    facet_wrap(~district) +
  theme(strip.text.x = element_text(size=14, color="black",face="bold"),
        strip.text.y = element_text(size=14, color="black",face="bold"),
        strip.background = element_rect(fill = "white", colour = "black"),
        panel.background = element_rect(fill = "aliceblue", color = "black"),
        panel.grid.major = element_line(color = "gray", linetype = "dotted"),
        panel.grid.minor = element_line(color = "gray", linetype = "dotted"),
        axis.text = element_text(size=12, color="black"),
        axis.title = element_text(size=14, color="black"),
        legend.title = element_text(size = 14),
        legend.text=element_text(size=12),
        plot.title = element_text(size = 16),
        plot.caption = element_text(size=12, hjust = 0))

f_hs_plot

ggsave(f_hs_plot, file = "Plots\\Female_hs_ab.jpg", width = 10, height = 6)

```

```{r, pk abundances}

rs_catch_f <- catch_females %>% filter(trap %in% c("pkin", "pkout"))
rs_catch_f$site <- as.factor(rs_catch_f$site)

rs_1 <- glmmTMB(count_f ~ interface + season + district + trap + (1|site),
            data = rs_catch_f, family = nbinom2)
summary(rs_1) #AIC 710.5

rs_2 <- glmmTMB(count_f ~ interface + district + trap + (1|site),
           data = rs_catch_f, family = nbinom2)
summary(rs_2) #AIC 708.8

rs_3 <- glmmTMB(count_f ~ season + district + trap + (1|site),
            data = rs_catch_f, family = nbinom2)
summary(rs_3) #AIC 710.5

rs_4 <- glmmTMB(count_f ~ interface + season + trap + (1|site),
            data = rs_catch_f, family = nbinom2)
summary(rs_4) #AIC 707.7

rs_5 <- glmmTMB(count_f ~ interface + district + season + (1|site),
            data = rs_catch_f, family = nbinom2)
summary(rs_5) #AIC 710.1

rs_6 <- glmmTMB(count_f ~ interface * trap + district + season + (1|site),
            data = rs_catch_f, family = nbinom2)
summary(rs_6) #AIC 711.1

anova(rs_1, rs_2) #season ns
anova(rs_1, rs_3) #interface ns (borderline)
anova(rs_1, rs_4) #district ns
anova(rs_1, rs_5) #trap ns
anova(rs_1, rs_6) #interaction trap v interface ns

#dharma package, model diagnostics
simulationOutput <- simulateResiduals(fittedModel = rs_1)
plot(simulationOutput)
testZeroInflation(simulationOutput) #ratio greater than 1 means more zeros than expected but when p is not significant not critical to need a zeroinfl model
testUniformity(simulationOutput)
testDispersion(simulationOutput, alternative = "greater") #tests for only overdispersion
testDispersion(simulationOutput)

# summary(glht(rs_1, linfct=mcp(season='Tukey')))
# summary(glht(rs_1, linfct=mcp(interface='Tukey')))
# summary(glht(rs_1, linfct=mcp(district='Tukey')))
# summary(glht(rs_1, linfct=mcp(trap='Tukey')))

#make predictions from best model
rsf_ab_pred <- ggpredict(rs_1, c("interface", "season", "district", "trap"))
rsf_ab_pred <- rsf_ab_pred %>% dplyr::rename(interface = x, mean = predicted, season = group, district = facet, trap = panel)

#remove predictions that dont make sense
rsf_ab_pred <- rsf_ab_pred %>%
  filter(!(interface == "Dumpsite" & district == "Wiang Kaen")) %>%
  filter(!(interface == "Temple" & district == "Wiang Kaen")) %>%
  filter(!(interface == "Temple" & district == "Mae Fah Luang")) %>%
  filter(!(interface == "Bat cave" & district == "Mae Fah Luang")) %>%
  filter(!(trap == "pkin" & interface %in% c("Dumpsite", "Orchard", "Forest edge", "Forest interior"))) %>%
  filter(!(trap == "pkout" & interface == "Bat cave"))

rsf_ab_pred$interface <- factor(rsf_ab_pred$interface, levels = c("House", "Temple", "Dumpsite", "Orchard", "Forest edge",  "Forest interior", "Bat cave"))
rsf_ab_pred$district <- factor(rsf_ab_pred$district, levels = c("Pong Nam Ron & Soi Dao", "Mae Fah Luang", "Wiang Kaen"))

pkin <- rsf_ab_pred %>% filter(trap == "pkin")
pkout <- rsf_ab_pred %>% filter(trap == "pkout")

f_pkin_plot <- ggplot(pkin, aes(fill=interface, y=mean, x=season)) + 
    geom_bar(position="dodge", color = "grey", stat="identity") + 
    geom_errorbar(pkin, mapping = aes(ymin=conf.low, ymax=conf.high), linewidth=0.1, position = "dodge") +
    scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(0, 1, 2, 5, 10, 20, 50, 100, 250, 500, 1000)) +
    scale_fill_brewer(palette="PRGn", name="Interface") +
    labs(fill = "Season", x = "Season", y = "Mean per collection (95% CI)", title = "Female indoor resting abundances at study interfaces") +
    facet_wrap(~district) + my_theme
  
f_pkin_plot

f_pkout_plot <- ggplot(pkout, aes(fill=interface, y=mean, x=season)) + 
    geom_bar(position="dodge", color = "grey", stat="identity") + 
    geom_errorbar(pkout, mapping = aes(ymin=conf.low, ymax=conf.high), linewidth=0.1, position = "dodge") +
    scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(0, 1, 2, 5, 10, 20, 50, 100, 250, 500, 1000)) +
    scale_fill_brewer(palette="PRGn", name="Interface") +
    labs(fill = "Season", x = "Season", y = "Mean per collection (95% CI)", title = "Female outdoor resting abundances at study interfaces") +
    facet_wrap(~district) + my_theme
  
f_pkout_plot

ggsave(f_pkin_plot, file = "Plots\\Female_rs_in_ab.jpg", width = 10, height = 6)
ggsave(f_pkout_plot, file = "Plots\\Female_rs_out.jpg", width = 10, height = 6)

```

```{r, ovi abundances}

ovi_catch_f <- catch_females %>% filter(trap %in% c("gat"))
ovi_catch_f$site <- as.factor(ovi_catch_f$site)

ovi_1 <- glmmTMB(count_f ~ interface + season + district + (1|site),
            data = ovi_catch_f, family = nbinom2)
summary(ovi_1) #AIC 418.5

ovi_2 <- glmmTMB(count_f ~ interface + district + (1|site),
           data = ovi_catch_f, family = nbinom2)
summary(ovi_2) #AIC 423.0

ovi_3 <- glmmTMB(count_f ~ season + district + (1|site),
            data = ovi_catch_f, family = nbinom2)
summary(ovi_3) #AIC 411.1

ovi_4 <- glmmTMB(count_f ~ interface + season + (1|site),
            data = ovi_catch_f, family = nbinom2)
summary(ovi_4) #AIC 423.2

anova(ovi_1, ovi_2) #season sig. (Chisq = 6.5, Df = 1, p < 0.05)
anova(ovi_1, ovi_3) #interface ns (Chisq = 4.6, Df = 6, p > 0.05)
anova(ovi_1, ovi_4) #district sig. (Chisq = 8.7, Df = 2, p < 0.05)

#dharma package, model diagnostics
simulationOutput <- simulateResiduals(fittedModel = ovi_1)
plot(simulationOutput)
testZeroInflation(simulationOutput) #ratio greater than 1 means more zeros than expected but when p is not significant not critical to need a zeroinfl model
testUniformity(simulationOutput)
testDispersion(simulationOutput, alternative = "greater") #tests for only overdispersion
testDispersion(simulationOutput)

summary(glht(ovi_1, linfct=mcp(season='Tukey'))) # p < 0.01
summary(glht(ovi_1, linfct=mcp(interface='Tukey'))) # p > 0.05
summary(glht(ovi_1, linfct=mcp(district='Tukey'))) # sig.

#make predictions from best model
ovif_ab_pred <- ggpredict(ovi_1, c("interface", "season", "district"))
ovif_ab_pred <- ovif_ab_pred %>% dplyr::rename(interface = x, mean = predicted, season = group, district = facet)

#remove predictions that dont make sense
ovif_ab_pred <- ovif_ab_pred %>%
  filter(!(interface == "Dumpsite" & district == "Wiang Kaen")) %>%
  filter(!(interface == "Temple" & district == "Wiang Kaen")) %>%
  filter(!(interface == "Temple" & district == "Mae Fah Luang")) %>%
  filter(!(interface == "Bat cave" & district == "Mae Fah Luang"))

ovif_ab_pred$interface <- factor(ovif_ab_pred$interface, levels = c("House", "Temple", "Dumpsite", "Orchard", "Forest edge",  "Forest interior", "Bat cave"))
ovif_ab_pred$district <- factor(ovif_ab_pred$district, levels = c("Pong Nam Ron & Soi Dao", "Mae Fah Luang", "Wiang Kaen"))

f_ovi_plot <- ggplot(ovif_ab_pred, aes(fill=interface, y=mean, x=season)) + 
    geom_bar(position="dodge", color = "grey", stat="identity") + 
    geom_errorbar(ovif_ab_pred, mapping = aes(ymin=conf.low, ymax=conf.high), linewidth=0.1, position = "dodge") +
    scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(0, 1, 2, 5, 10, 20, 50, 100, 250, 500, 1000)) +
    scale_fill_brewer(palette="PRGn", name="Interface") +
    labs(fill = "Season", x = "Season", y = "Mean per collection (95% CI)", title = "Female ovipositing abundances at study interfaces") +
    facet_wrap(~district) + my_theme

f_ovi_plot

ggsave(f_ovi_plot, file = "Plots\\Female_ovi_ab.jpg", width = 10, height = 6)

```

```{r, restructure vector data} 

female_wide <- mq_females %>% 
  group_by(start_date, start_time, district, interface, group, trap, season, site, closest_id) %>%
  summarise(sum_f = sum(females)) %>%
  pivot_wider(names_from = c(closest_id), values_from = c(sum_f)) %>%
  mutate_at(c(9:93), ~replace(., is.na(.), 0)) %>%
  as.data.frame() %>%
  dplyr::select("start_date", "start_time", "district", "interface", "group", "trap", "season", "site", "empty", "Aedes (Ste.) albopictus", "Aedes (Ste.) aegypti", "Culex (Cux.) vishnui", "Culex (Cux.) gelidus", "Culex (Cux.) quinquefasciatus", "Culex (Cux.) tritaeniorhynchus", "Culex (Cux.) fuscocephala", "Culex  (Cux.) pseudovishnui") %>%
  dplyr::rename("albo_ab" ="Aedes (Ste.) albopictus", 
                "aegy_ab" = "Aedes (Ste.) aegypti",
                "vish_ab" = "Culex (Cux.) vishnui",
                "gel_ab" = "Culex (Cux.) gelidus",
                "quin_ab" = "Culex (Cux.) quinquefasciatus",
                "tri_ab" = "Culex (Cux.) tritaeniorhynchus",
                "fusco_ab" = "Culex (Cux.) fuscocephala",
                "psvish_ab" = "Culex  (Cux.) pseudovishnui") %>%
  mutate(albo_p = ifelse(albo_ab > 0, 1, 0),
         aegy_p = ifelse(aegy_ab > 0, 1, 0),
         vish_p = ifelse(vish_ab > 0, 1, 0),
         gel_p = ifelse(gel_ab > 0, 1, 0),
         quin_p = ifelse(quin_ab > 0, 1, 0),
         tri_p = ifelse(tri_ab > 0, 1, 0),
         fusco_p = ifelse(fusco_ab > 0, 1, 0),
         psvish_p = ifelse(psvish_ab > 0, 1, 0))

cols <- c("interface", "district", "season", "site", "trap")
vector_hs <- female_wide %>%
  filter(trap %in% c("lt", "pro")) %>%
  mutate_at(cols, factor)

```

```{r, vector probability} 

#ae.albo
albo_p_1 <- glmmTMB(albo_p ~ district + interface + season, data = vector_hs, family = "binomial")
summary(albo_p_1) #AIC 401.9
albo_p_2 <- glmmTMB(albo_p ~ interface + season, data = vector_hs, family = "binomial")
summary(albo_p_2) #AIC 400.0
albo_p_3 <- glmmTMB(albo_p ~ district + season, data = vector_hs, family = "binomial")
summary(albo_p_3) #AIC 402.8
albo_p_4 <- glmmTMB(albo_p ~ interface + district, data = vector_hs, family = "binomial")
summary(albo_p_4) #AIC 405.5

anova(albo_p_1, albo_p_2) #district, ns but keep in as want to plot values by district (Chisq = 2.7, Df = 2, p > 0.05)
anova(albo_p_1, albo_p_3) #interface sig *** (Chisq = 36.1, Df = 6, p < 0.001)
anova(albo_p_1, albo_p_4) #season sig * (Chisq = 5.7, Df = 1, p < 0.05)

#dharma check model
simulationOutput <- simulateResiduals(fittedModel = albo_p_1)
plot(simulationOutput)
testZeroInflation(simulationOutput)
testUniformity(simulationOutput)
testDispersion(simulationOutput, alternative = "greater") #tests for only overdispersion
testDispersion(simulationOutput)

#make predictions from best model
albo_pred <- ggpredict(albo_p_1, c("district", "interface", "season"), bias_correction = TRUE)
summary(glht(albo_p_1, linfct=mcp(interface='Tukey'))) # sig DU and FI more abundant than HO
summary(glht(albo_p_1, linfct=mcp(season='Tukey'))) # sig *
summary(glht(albo_p_1, linfct=mcp(district='Tukey'))) #ns

#cx.vish
#remove bat cave from df as 0 vishnui caught there
vector_hs_bc <- vector_hs %>% filter(!interface == "Bat cave")
vish_p_1 <- glmmTMB(vish_p ~ district + interface + season, data = vector_hs_bc, family = "binomial")
summary(vish_p_1) #AIC 282
vish_p_2 <- glmmTMB(vish_p ~ interface + season, data = vector_hs_bc, family = "binomial")
summary(vish_p_2) #AIC 342.5
vish_p_3 <- glmmTMB(vish_p ~ district + season, data = vector_hs_bc, family = "binomial")
summary(vish_p_3) #AIC 299.3
vish_p_4 <- glmmTMB(vish_p ~ interface + district, data = vector_hs_bc, family = "binomial")
summary(vish_p_4) #AIC 290.6

anova(vish_p_1, vish_p_2) #district sig *** (Chisq = 64.5, Df = 2, p < 0.001)
anova(vish_p_1, vish_p_3) #interface sig *** (Chisq = 27.3, Df = 5, p < 0.001)
anova(vish_p_1, vish_p_4) #season sig ** (Chisq = 10.6, Df = 1, p < 0.01)

#dharma check model
simulationOutput <- simulateResiduals(fittedModel = vish_p_1)
plot(simulationOutput)
testZeroInflation(simulationOutput)
testUniformity(simulationOutput)
testDispersion(simulationOutput, alternative = "greater") #tests for only overdispersion
testDispersion(simulationOutput)

#make predictions from best model
vish_pred <- ggpredict(vish_p_1, c("district", "interface", "season"))

summary(glht(vish_p_1, linfct=mcp(interface='Tukey'))) #sig ***
summary(glht(vish_p_1, linfct=mcp(season='Tukey'))) #sig **
summary(glht(vish_p_1, linfct=mcp(district='Tukey'))) #sig ***

#cx.gel
#remove WK and bat cave from df as 0 gelidus caught there
vector_hs_wk <- vector_hs_bc %>% filter(!district == "Wiang Kaen")
gel_p_1 <- glmmTMB(gel_p ~ district + interface + season, data = vector_hs_wk, family = "binomial")
summary(gel_p_1) #AIC 151.2
gel_pred <- ggpredict(gel_p_1, c("district", "interface", "season"))

gel_p_2 <- glmmTMB(gel_p ~ interface + season, data = vector_hs_wk, family = "binomial")
summary(gel_p_2) #AIC 149.6
gel_pred <- ggpredict(gel_p_2, c("interface", "season"))

gel_p_3 <- glmmTMB(gel_p ~ district + season, data = vector_hs_wk, family = "binomial")
summary(gel_p_3) #AIC 178.6
gel_pred <- ggpredict(gel_p_3, c("district", "season"))

gel_p_4 <- glmmTMB(gel_p ~ interface + district, data = vector_hs_wk, family = "binomial")
summary(gel_p_4) #AIC 238.4
gel_pred <- ggpredict(gel_p_4, c("district", "interface"))

anova(gel_p_1, gel_p_2) #district ns (Chisq = 0.4, Df = 1, p > 0.05)
anova(gel_p_1, gel_p_3) #interface sig *** (Chisq = 37.4, Df = 5, p < 0.001)
anova(gel_p_1, gel_p_4) #season sig *** (Chisq = 89.2, Df = 1, p < 0.001)

#dharma check model
simulationOutput <- simulateResiduals(fittedModel = gel_p_1)
plot(simulationOutput)
testZeroInflation(simulationOutput)
testUniformity(simulationOutput)
testDispersion(simulationOutput, alternative = "greater") #tests for only overdispersion
testDispersion(simulationOutput)

#make predictions from best model
gel_pred <- ggpredict(gel_p_1, c("district", "interface", "season"))

summary(glht(gel_p_1, linfct=mcp(interface='Tukey'))) #sig ***
summary(glht(gel_p_1, linfct=mcp(season='Tukey'))) #sig ***
summary(glht(gel_p_1, linfct=mcp(district='Tukey'))) #ns

#cx.quin
quin_p_1 <- glmmTMB(quin_p ~ district + interface + season, data = vector_hs, family = "binomial")
summary(quin_p_1) #AIC 358.5
quin_pred <- ggpredict(quin_p_1, c("district", "interface", "season"))

quin_p_2 <- glmmTMB(quin_p ~ interface + season, data = vector_hs, family = "binomial")
summary(quin_p_2) #AIC 363.0
quin_pred <- ggpredict(quin_p_2, c("interface", "season"))

quin_p_3 <- glmmTMB(quin_p ~ district + season, data = vector_hs, family = "binomial")
summary(quin_p_3) #AIC 355.2
quin_pred <- ggpredict(quin_p_3, c("district", "season"))

quin_p_4 <- glmmTMB(quin_p ~ interface + district, data = vector_hs, family = "binomial")
summary(quin_p_4) #AIC 356.6
quin_pred <- ggpredict(quin_p_4, c("district", "interface"))

anova(quin_p_1, quin_p_2) #district sig * (Chisq = 8.5, Df = 2, p < 0.01)
anova(quin_p_1, quin_p_3) #interface ns (Chisq = 8.7, Df = 6, p > 0.05)
anova(quin_p_1, quin_p_4) #season ns (Chisq = 0.1, Df = 1, p > 0.05)

#dharma check model
simulationOutput <- simulateResiduals(fittedModel = quin_p_1)
plot(simulationOutput)
testZeroInflation(simulationOutput)
testUniformity(simulationOutput)
testDispersion(simulationOutput, alternative = "greater") #tests for only overdispersion
testDispersion(simulationOutput)

#make predictions from best model
quin_pred <- ggpredict(quin_p_1, c("district", "interface", "season"))

summary(glht(quin_p_1, linfct=mcp(interface='Tukey'))) #ns
summary(glht(quin_p_1, linfct=mcp(season='Tukey'))) #ns
summary(glht(quin_p_1, linfct=mcp(district='Tukey'))) #sig *

#cx.tri
#remove WK and bat cave from df as 0 tritaeniorhynchus caught there
tri_p_1 <- glmmTMB(tri_p ~ district + interface + season, data = vector_hs_wk, family = "binomial")
summary(tri_p_1) #AIC 152.6
tri_pred <- ggpredict(tri_p_1, c("district", "interface", "season"))

tri_p_2 <- glmmTMB(tri_p ~ interface + season, data = vector_hs_wk, family = "binomial")
summary(tri_p_2) #AIC 163.6
tri_pred <- ggpredict(tri_p_2, c("interface", "season"))

tri_p_3 <- glmmTMB(tri_p ~ district + season, data = vector_hs_wk, family = "binomial")
summary(tri_p_3) #AIC 159.2
tri_pred <- ggpredict(tri_p_3, c("district", "season"))

tri_p_4 <- glmmTMB(tri_p ~ interface + district, data = vector_hs_wk, family = "binomial")
summary(tri_p_4) #AIC 155.8
tri_pred <- ggpredict(tri_p_4, c("district", "interface"))

anova(tri_p_1, tri_p_2) #district sig *** (Chisq = 13.1, Df = 1, p < 0.001)
anova(tri_p_1, tri_p_3) #interface sig ** (Chisq = 16.6, Df = 5, p < 0.01)
anova(tri_p_1, tri_p_4) #season sig * (Chisq = 5.3, Df = 1, p < 0.05)

#dharma check model
simulationOutput <- simulateResiduals(fittedModel = tri_p_1)
plot(simulationOutput)
testZeroInflation(simulationOutput)
testUniformity(simulationOutput)
testDispersion(simulationOutput, alternative = "greater") #tests for only overdispersion
testDispersion(simulationOutput)

#make predictions from best model
tri_pred <- ggpredict(tri_p_1, c("district", "interface", "season"))

summary(glht(tri_p_1, linfct=mcp(interface='Tukey'))) #ns
summary(glht(tri_p_1, linfct=mcp(season='Tukey'))) #sig *
summary(glht(tri_p_1, linfct=mcp(district='Tukey'))) #sig **

#ae.aegy #not enough data 18/310
#aegy_p_1 <- glmmTMB(aegy_p ~ district + interface + season + (1|site), data = vector_hs, family = "binomial")
#cx.fusco ## not enough data
#fusco_p_1 <- glmmTMB(fusco_p ~ district + interface + season, data = vector_hs, family = "binomial")
#cx.pseudo ## not enough data
#psvish_p_1 <- glmmTMB(psvish_p ~ district + interface + season, data = vector_hs, family = "binomial")

```

```{r, plot vector prob} 
#create vector name col
albo_pred$Species <- "Ae.albopictus"
vish_pred$Species <- "Cx.vishnui"
gel_pred$Species <- "Cx.gelidus"
quin_pred$Species <- "Cx.quinquefasciatus"
tri_pred$Species <- "Cx.tritaeniorhynchus"

#bring model predictions to one df
vector_prob <- rbind(albo_pred, vish_pred, gel_pred, quin_pred, tri_pred)
vector_prob <- vector_prob %>%
  dplyr::rename(District = x, Probability = predicted, Interface = group, Season = facet) %>%
  as.data.frame() 

pnr_vectors <- vector_prob %>%
  filter(District == "Pong Nam Ron & Soi Dao")
wk_vectors <- vector_prob %>%
  filter(District == "Wiang Kaen") %>%
  filter(!(Interface %in% c("Temple", "Dumpsite")))
mfl_vectors <- vector_prob %>%
  filter(District == "Mae Fah Luang")  %>%
  filter(!(Interface %in% c("Temple", "Bat cave")))

pnr_vectors$Interface <- factor(pnr_vectors$Interface, levels = c("House", "Temple", "Dumpsite", "Orchard", "Forest edge",  "Forest interior", "Bat cave"))
wk_vectors$Interface <- factor(wk_vectors$Interface, levels = c("House", "Orchard", "Forest edge",  "Forest interior", "Bat cave"))
mfl_vectors$Interface <- factor(mfl_vectors$Interface, levels = c("House", "Dumpsite", "Orchard", "Forest edge",  "Forest interior"))

pnr_plot <- ggplot(pnr_vectors, aes(x = Interface, y = Probability, group = Season, color = Season)) +
  geom_line() +
  geom_point() +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill = Season), alpha=0.3) +
  ylim(0,1) +
  labs(title = "Pong Nam Ron & Soi Dao", y = "Probability (95% CI)") +
  facet_wrap(~Species, nrow=5) + 
  theme_classic() +
  theme(axis.text.x = element_text(size=13, color="black"),
        axis.text.y = element_text(size=13, color="black"),
        axis.title = element_text(size=15, color="black"),
        legend.title = element_text(size = 15),
        legend.text=element_text(size=13),
        plot.title = element_text(size = 17, face = "bold"),
        strip.text.x = element_text(size=15, color="black", face ="italic")) +
   scale_x_discrete(labels = function(x) 
    stringr::str_wrap(x, width = 6))
pnr_plot

wk_plot <- ggplot(wk_vectors, aes(x = Interface, y = Probability, group = Season, color = Season)) +
  geom_line() +
  geom_point() +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill = Season), alpha=0.3) +
  ylim(0,1) +
  labs(title = "Wiang Kaen", y = "Probability (95% CI)") +
  facet_wrap(~Species, nrow=5) + 
  theme_classic() +
  theme(axis.text.x = element_text(size=13, color="black"),
        axis.text.y = element_text(size=13, color="black"),
        axis.title = element_text(size=15, color="black"),
        legend.title = element_text(size = 15),
        legend.text=element_text(size=13),
        plot.title = element_text(size = 17, face = "bold"),
        strip.text.x = element_text(size=15, color="black", face ="italic")) +
   scale_x_discrete(labels = function(x) 
    stringr::str_wrap(x, width = 6))
wk_plot

mfl_plot <- ggplot(mfl_vectors, aes(x = Interface, y = Probability, group = Season, color = Season)) +
  geom_line() +
  geom_point() +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill = Season), alpha=0.3) +
  ylim(0,1) +
  labs(title = "Mae Fah Luang", y = "Probability (95% CI)") +
  facet_wrap(~Species, nrow=5) + theme_classic() +
  theme(axis.text.x = element_text(size=13, color="black"),
        axis.text.y = element_text(size=13, color="black"),
        axis.title = element_text(size=15, color="black"),
        legend.title = element_text(size = 15),
        legend.text=element_text(size=13),
        plot.title = element_text(size = 17, face = "bold"),
        strip.text.x = element_text(size=15, color="black", face ="italic")) +
   scale_x_discrete(labels = function(x) 
    stringr::str_wrap(x, width = 6))
mfl_plot

#make a single plot
title_gg <- ggplot() +
  labs(title = "") +
  theme(panel.background = element_rect(fill = "white", color = "white"),
        plot.title = element_text(size = 20, hjust = 0.5))
prob_plot <- plot_grid(pnr_plot, mfl_plot, wk_plot, ncol = 3, rel_widths = c(1, 0.9, 0.8))
all_prob_plot <- plot_grid(title_gg, prob_plot, ncol = 1, rel_heights = c(0.05, 1))
all_prob_plot
ggsave(all_prob_plot, file = "Plots\\Arbo_vector_presence_all_districts.jpg", width = 20, height = 10)

ggsave(pnr_plot, file = "Plots\\Arbo_vector_presence_cth.jpg", width = 8, height = 9)
ggsave(mfl_plot, file = "Plots\\Arbo_vector_presence_mfl.jpg", width = 8, height = 9)
ggsave(wk_plot, file = "Plots\\Arbo_vector_presence_wk.jpg", width = 7, height = 6)

```

```{r, vector abundance}

#use vector_hs dataset to look at abundance
#start with ae.albo
  
albo_ab_1 <- glmmTMB(albo_ab ~ district + interface + season, vector_hs, family = nbinom2)
summary(albo_ab_1) #AIC = 1078 *best mod, tried trap as RE but got singularity
albo_ab_2 <- glmmTMB(albo_ab ~ district + interface, vector_hs, family = nbinom2)
summary(albo_ab_2) #AIC = 1086
albo_ab_3 <- glmmTMB(albo_ab ~ district + season, vector_hs, family = nbinom2)
summary(albo_ab_3) #AIC = 1122
albo_ab_4 <- glmmTMB(albo_ab ~ interface + season, vector_hs, family = nbinom2)
summary(albo_ab_4) #AIC = 1079

anova(albo_ab_1, albo_ab_2)  #season sig *** (Chisq = 9.7, Df = 1, p < 0.01)
anova(albo_ab_1, albo_ab_3)  #interface sig *** (Chisq = 56.0, Df = 6, p < 0.001)
anova(albo_ab_1, albo_ab_4)  #district n sig (Chisq = 5.0, Df = 2, p > 0.05)

simulationOutput <- simulateResiduals(fittedModel = albo_ab_1, plot = T)
plot(simulationOutput)
testZeroInflation(simulationOutput) #compares observed number of 0s with 0s expected from simulations, ns
testDispersion(simulationOutput)

#make predictions from albo_ab_1
albo_ab_pred <- ggpredict(albo_ab_1, c("district", "interface", "season"))
albo_ab_pred <- albo_ab_pred %>%
  dplyr::rename(District = x, Mean = predicted, Interface = group, Season = facet) %>%
  as.data.frame() %>%
  mutate(District = if_else(District == "Chanthaburi", "Pong Nam Ron & Soi Dao", District))

summary(glht(albo_ab_1,linfct=mcp(season='Tukey')))
summary(glht(albo_ab_1,linfct=mcp(interface='Tukey')))
summary(glht(albo_ab_1,linfct=mcp(district='Tukey')))

#next cx.vish
#using dataset without BC as no vishnui caught there
  
vish_ab_1 <- glmmTMB(vish_ab ~ district + interface + season, vector_hs_bc, family = nbinom2)
summary(vish_ab_1) #AIC 977
vish_ab_2 <- glmmTMB(vish_ab ~ district + interface, vector_hs_bc, family = nbinom2)
summary(vish_ab_2) #AIC 1021
vish_ab_3 <- glmmTMB(vish_ab ~ district + season, vector_hs_bc, family = nbinom2)
summary(vish_ab_3) #AIC 1039
vish_ab_4 <- glmmTMB(vish_ab ~ interface + season, vector_hs_bc, family = nbinom2)
summary(vish_ab_4) #AIC 1046

anova(vish_ab_1, vish_ab_2)  #season sig *** (Chisq = 45.1, Df = 1, p < 0.001)
anova(vish_ab_1, vish_ab_3)  #interface sig *** (Chisq = 71.0, Df = 5, p < 0.001)
anova(vish_ab_1, vish_ab_4)  #district sig *** (Chisq = 72.6, Df = 2, p < 0.001)

simulationOutput <- simulateResiduals(fittedModel = vish_ab_1, plot = T)
plot(simulationOutput)
testZeroInflation(simulationOutput) #compares observed number of 0s with 0s expected from simulations, ns
testDispersion(simulationOutput)

#make predictions from vish_ab_1
vish_ab_pred <- ggpredict(vish_ab_1, c("district", "interface", "season"))
vish_ab_pred <- vish_ab_pred %>%
  dplyr::rename(District = x, Mean = predicted, Interface = group, Season = facet) %>%
  as.data.frame() %>%
  mutate(District = if_else(District == "Chanthaburi", "Pong Nam Ron & Soi Dao", District))

summary(glht(vish_ab_1,linfct=mcp(season='Tukey'))) #sig ***
summary(glht(vish_ab_1,linfct=mcp(interface='Tukey'))) #sig ***
summary(glht(vish_ab_1,linfct=mcp(district='Tukey'))) #sig ***

#next cx.gel
#using dataset without BC and WK as no gelidus caught there
  
gel_ab_1 <- glmmTMB(gel_ab ~ district + interface + season, vector_hs_wk, family = nbinom2)
summary(gel_ab_1) #AIC 720
gel_ab_2 <- glmmTMB(gel_ab ~ district + interface, vector_hs_wk, family = nbinom2)
summary(gel_ab_2) #AIC 806
gel_ab_3 <- glmmTMB(gel_ab ~ district + season, vector_hs_wk, family = nbinom2)
summary(gel_ab_3) #AIC 757
gel_ab_4 <- glmmTMB(gel_ab ~ interface + season, vector_hs_wk, family = nbinom2)
summary(gel_ab_4) #AIC 720

anova(gel_ab_1, gel_ab_2)  #season sig *** (Chisq = 88.4, Df = 1, p < 0.001)
anova(gel_ab_1, gel_ab_3)  #interface sig *** (Chisq = 47.4, Df = 5, p < 0.001)
anova(gel_ab_1, gel_ab_4)  #district ns (Chisq = 1.9, Df = 1, p > 0.05)

simulationOutput <- simulateResiduals(fittedModel = gel_ab_1, plot = T)
plot(simulationOutput)
testZeroInflation(simulationOutput) #compares observed number of 0s with 0s expected from simulations, ns
testDispersion(simulationOutput)

#make predictions from gel_ab_1
gel_ab_pred <- ggpredict(gel_ab_1, c("district", "interface", "season"))
gel_ab_pred <- gel_ab_pred %>%
  dplyr::rename(District = x, Mean = predicted, Interface = group, Season = facet) %>%
  as.data.frame() %>%
  mutate(District = if_else(District == "Chanthaburi", "Pong Nam Ron & Soi Dao", District))

summary(glht(gel_ab_1,linfct=mcp(season='Tukey'))) #sig ***
summary(glht(gel_ab_1,linfct=mcp(interface='Tukey'))) #sig ***
summary(glht(gel_ab_1,linfct=mcp(district='Tukey'))) #ns

#cx.quin
  
quin_ab_1 <- glmmTMB(quin_ab ~ district + interface + season, vector_hs, family = nbinom2)
summary(quin_ab_1) #AIC = 675 
quin_ab_2 <- glmmTMB(quin_ab ~ district + interface, vector_hs, family = nbinom2)
summary(quin_ab_2) #AIC = 677
quin_ab_3 <- glmmTMB(quin_ab ~ district + season, vector_hs, family = nbinom2)
summary(quin_ab_3) #AIC = 678
quin_ab_4 <- glmmTMB(quin_ab ~ interface + season, vector_hs, family = nbinom2)
summary(quin_ab_4) #AIC = 688
quin_ab_5 <- glmmTMB(quin_ab ~ district + interface + season, vector_hs, family = nbinom2)
summary(quin_ab_5) #AIC = 675

anova(quin_ab_1, quin_ab_2)  #season sig * (Chisq = 4.1, Df = 1, p < 0.05)
anova(quin_ab_1, quin_ab_3)  #interface sig * (Chisq = 14.7, Df = 6, p < 0.05)
anova(quin_ab_1, quin_ab_4)  #district sig *** (Chisq = 17.5, Df = 2, p > 0.001)

simulationOutput <- simulateResiduals(fittedModel = quin_ab_1, plot = T)
plot(simulationOutput)
testZeroInflation(simulationOutput) #compares observed number of 0s with 0s expected from simulations, ns
testDispersion(simulationOutput)

#make predictions from quin_ab_1
quin_ab_pred <- ggpredict(quin_ab_1, c("district", "interface", "season"))
quin_ab_pred <- quin_ab_pred %>%
  dplyr::rename(District = x, Mean = predicted, Interface = group, Season = facet) %>%
  as.data.frame() %>%
  mutate(District = if_else(District == "Chanthaburi", "Pong Nam Ron & Soi Dao", District))

summary(glht(quin_ab_1,linfct=mcp(season='Tukey'))) #sig *
summary(glht(quin_ab_1,linfct=mcp(interface='Tukey'))) # sig *
summary(glht(quin_ab_1,linfct=mcp(district='Tukey'))) #sig ***

#next cx.tri
#using dataset without BC and WK as no tritaeniorhynchus caught there
  
tri_ab_1 <- glmmTMB(tri_ab ~ district + interface + season, vector_hs_wk, family = nbinom2)
summary(tri_ab_1) #AIC 296
tri_ab_2 <- glmmTMB(tri_ab ~ district + interface, vector_hs_wk, family = nbinom2)
summary(tri_ab_2) #AIC 299
tri_ab_3 <- glmmTMB(tri_ab ~ district + season, vector_hs_wk, family = nbinom2)
summary(tri_ab_3) #AIC 305
tri_ab_4 <- glmmTMB(tri_ab ~ interface + season, vector_hs_wk, family = nbinom2)
summary(tri_ab_4) #AIC 294

anova(tri_ab_1, tri_ab_2)  #season sig * (Chisq = 4.7, Df = 1, p < 0.05)
anova(tri_ab_1, tri_ab_3)  #interface sig ** (Chisq = 18.8, Df = 5, p < 0.01)
anova(tri_ab_1, tri_ab_4)  #district ns (Chisq = 0.3, Df = 1, p > 0.05)

simulationOutput <- simulateResiduals(fittedModel = tri_ab_1, plot = T)
plot(simulationOutput)
testZeroInflation(simulationOutput) #compares observed number of 0s with 0s expected from simulations, ns
testDispersion(simulationOutput)

#make predictions from tri_ab_1
tri_ab_pred <- ggpredict(tri_ab_1, c("district", "interface", "season"))
tri_ab_pred <- tri_ab_pred %>%
  dplyr::rename(District = x, Mean = predicted, Interface = group, Season = facet) %>%
  as.data.frame() %>%
  mutate(District = if_else(District == "Chanthaburi", "Pong Nam Ron & Soi Dao", District))

summary(glht(tri_ab_1,linfct=mcp(season='Tukey'))) #sig *
summary(glht(tri_ab_1,linfct=mcp(interface='Tukey'))) #sig **
summary(glht(tri_ab_1,linfct=mcp(district='Tukey'))) #ns

```

```{r, vector ab plots}

albo_ab_pred$Interface <- factor(albo_ab_pred$Interface, levels = c("House", "Temple", "Dumpsite", "Orchard", "Forest edge",  "Forest interior", "Bat cave"))
albo_ab_pred$District <- factor(albo_ab_pred$District, levels = c("Pong Nam Ron & Soi Dao", "Mae Fah Luang", "Wiang Kaen"))

vish_ab_pred$Interface <- factor(vish_ab_pred$Interface, levels = c("House", "Temple", "Dumpsite", "Orchard", "Forest edge",  "Forest interior"))
vish_ab_pred$District <- factor(vish_ab_pred$District, levels = c("Pong Nam Ron & Soi Dao", "Mae Fah Luang", "Wiang Kaen"))

gel_ab_pred$Interface <- factor(gel_ab_pred$Interface, levels = c("House", "Temple", "Dumpsite", "Orchard", "Forest edge",  "Forest interior"))
gel_ab_pred$District <- factor(gel_ab_pred$District, levels = c("Pong Nam Ron & Soi Dao", "Mae Fah Luang"))

quin_ab_pred$Interface <- factor(quin_ab_pred$Interface, levels = c("House", "Temple", "Dumpsite", "Orchard", "Forest edge",  "Forest interior", "Bat cave"))
quin_ab_pred$District <- factor(quin_ab_pred$District, levels = c("Pong Nam Ron & Soi Dao", "Mae Fah Luang", "Wiang Kaen"))

tri_ab_pred$Interface <- factor(tri_ab_pred$Interface, levels = c("House", "Temple", "Dumpsite", "Orchard", "Forest edge",  "Forest interior", "Bat cave"))
tri_ab_pred$District <- factor(tri_ab_pred$District, levels = c("Pong Nam Ron & Soi Dao", "Mae Fah Luang"))

interface_colours <- c("House" = "#762A83", "Temple" = "#AF8DC3", "Dumpsite" = "#E7D4E8", "Orchard" =  "#F7F7F7",  "Forest edge" = "#D9F0D3", "Forest interior" = "#7FBF7B", "Bat cave" = "#1B7837")

# district.labs <- c("Wiang-Kaen (CR)", "Chanthaburi", "Mae Fah Luang (CR)")
# names(district.labs) <- c("Wiang Kaen", "Chanthaburi", "Mae Fah Luang")

ae.albo_plot <- ggplot(albo_ab_pred, aes(x = Season, y = Mean, fill = Interface)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "grey") +
  geom_errorbar(albo_ab_pred, mapping = aes(ymin=conf.low, ymax=conf.high), linewidth=0.1, position = position_dodge(.9), width = 0.2) +
  labs(fill = "Interface", x = "", y = "Mean per collection (95% CI)", title = "Aedes albopictus", caption = "") +
  # scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(0, 1, 2, 5, 10, 20, 50, 100, 250, 500, 1000)) +
  scale_y_continuous(breaks=seq(0,25,5), limits=c(0, 25)) +
  scale_fill_manual(values = interface_colours) + my_theme +
  theme(plot.title = element_text(face = "bold.italic"),
        legend.title = element_text(size = 16),
        legend.text=element_text(size=14)) +
 facet_grid(~District) +
  scale_x_discrete(labels = function(x) 
    stringr::str_wrap(x, width = 6)) +
  guides(color = FALSE)

ae.albo_plot

ggsave(ae.albo_plot, file = "Plots\\Ae.albo_ab_plot.jpg", width = 10, height = 6)

# cx.vish 

cx.vish_plot <- ggplot(vish_ab_pred, aes(x = Season, y = Mean, fill = Interface)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "grey") +
  geom_errorbar(vish_ab_pred, mapping = aes(ymin=conf.low, ymax=conf.high), linewidth=0.1, position = position_dodge(.9), width = 0.2) +
  # ylim(0, 12) +
  labs(fill = "Interface", x = "Season", y = "Mean per collection (95% CI)", title = "Culex vishnui", caption = "") +
  scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(0, 1, 2, 5, 10, 20, 50, 100, 250, 500, 1000)) +
  scale_fill_manual(values = interface_colours) + my_theme +
  theme(plot.title = element_text(face = "bold.italic")) +
  facet_grid(~District) +
  scale_x_discrete(labels = function(x) 
    stringr::str_wrap(x, width = 6)) +
  guides(color = FALSE)

cx.vish_plot

ggsave(cx.vish_plot, file = "Plots\\Cx.vish_ab_plot.jpg", width = 38, height = 14)

#cx.gel

cx.gel_plot <- ggplot(gel_ab_pred, aes(x = Season, y = Mean, fill = Interface)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "grey") +
  geom_errorbar(gel_ab_pred, mapping = aes(ymin=conf.low, ymax=conf.high), linewidth=0.1, position = position_dodge(.9), width = 0.2) +
  labs(fill = "Interface", x = "", y = "Mean per collection (95% CI)", title = "Culex gelidus", caption = "") +
  scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(0, 2, 5, 10, 25, 50, 100, 250, 500, 1000)) +
  scale_fill_manual(values = interface_colours) + my_theme +
  theme(legend.position="none",
        plot.title = element_text(face = "bold.italic")) +
  facet_grid(~District) +
  scale_x_discrete(labels = function(x) 
    stringr::str_wrap(x, width = 6)) +
  guides(color = FALSE)

cx.gel_plot

ggsave(cx.gel_plot, file = "Plots\\Cx.gel_plot.jpg", width = 38, height = 14)

#cx.quin

cx.quin_plot <- ggplot(quin_ab_pred, aes(x = Season, y = Mean, fill = Interface)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "grey") +
  geom_errorbar(quin_ab_pred, mapping = aes(ymin=conf.low, ymax=conf.high), linewidth=0.1, position = position_dodge(.9), width = 0.2) +
  labs(fill = "Interface", x = "", y = "Mean per collection (95% CI)", title = "Culex quinquefasciatus") +
  #scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(0, 1, 2, 5, 10, 20, 50, 100, 250, 500, 1000)) +
  # scale_y_continuous(breaks=seq(0,25,5), limits=c(0, 25)) +
  scale_fill_manual(values = interface_colours) +
  my_theme +
  theme(plot.title = element_text(face = "bold.italic")) +
  facet_grid(~District) +
  scale_x_discrete(labels = function(x) 
    stringr::str_wrap(x, width = 6)) +
  guides(color = FALSE)

cx.quin_plot

ggsave(cx.quin_plot, file = "Plots\\Cx.quin_plot.jpg", width = 10, height = 6)

#cx.tri

cx.tri_plot <- ggplot(tri_ab_pred, aes(x = Season, y = Mean, fill = Interface)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "grey") +
  geom_errorbar(tri_ab_pred, mapping = aes(ymin=conf.low, ymax=conf.high), linewidth=0.1, position = position_dodge(.9), width = 0.2) +
  labs(fill = "Interface", x = "", y = "Mean per collection (95% CI)", title = "Culex tritaeniorhynchus") +
  # scale_y_continuous(trans=scales::pseudo_log_trans(base = 10), breaks = c(0, 1, 2, 5, 10, 20, 50, 100, 250, 500, 1000)) +
  scale_y_continuous(breaks=seq(0,12,2), limits=c(0, 12)) +
  scale_fill_manual(values = interface_colours) +
  my_theme +
  theme(legend.position="none",
        plot.title = element_text(face = "bold.italic")) +
  facet_grid(~District) +
  scale_x_discrete(labels = function(x) 
    stringr::str_wrap(x, width = 6)) +
  guides(color = FALSE)

cx.tri_plot

ggsave(cx.tri_plot, file = "Plots\\Cx.tri_plot.jpg", width = 10, height = 6)

#bring all 5 plots in one

all_vector_plot <- plot_grid(ae.albo_plot + theme(legend.position = "none"),
                             cx.gel_plot + theme(legend.position = "none"),
                             cx.quin_plot + theme(legend.position = "none"),
                             cx.tri_plot + theme(legend.position = "none"),
                             cx.vish_plot + theme(legend.position = "none"),
                             align = 'v',
                             labels = c("A", "B", "C", "D", "E"),
                             ncol = 1)

legend <- get_legend(ae.albo_plot + theme(legend.position = "left"))    

full_vector_plot <- plot_grid(all_vector_plot, legend, ncol = 2, rel_widths = c(0.8, 0.2))
full_vector_plot           

ggsave(full_vector_plot, file = "Plots\\vector_ab_plot.jpg", width = 12, height = 18)

```

```{r, pnr diversity}

#format df for iNEXT

f_species_pnr <- mq_females %>% 
  filter(district == "Pong Nam Ron & Soi Dao") %>%
  dplyr::select(closest_id, interface, females) %>%
  group_by(closest_id, interface) %>%
  summarise(sum = sum(females)) %>% ungroup() %>%
  pivot_wider(names_from = "interface", 
              values_from = "sum") %>%
  mutate_at(c(2:8), ~replace(., is.na(.), 0)) %>%
  as.data.frame()
  
rownames(f_species_pnr) <- f_species_pnr$closest_id
f_species_pnr <- f_species_pnr[,-1]

# remove unidentified sp, kept Arm (Lei) spp. as represents a different species
f_species_pnr <- f_species_pnr %>%
  filter(!(rownames(.) %in% c("Aedes  spp.", "Aedes (Ste.) spp.", "Anopheles  spp.", "Anopheles (Cel.) spp.", "Armigeres  spp.", "Armigeres (Arm.) spp.", "Coquillettidia  spp.", "Coquillettidia (Coq.) spp.", "Culex   spp.", "Culex  (Cux.) spp.", "Culex  (Cux.) spp.", "Culex  spp.", "Culex (Cux.) spp.",  "Mansonia  spp.", "Mansonia (Mnd.) spp.", "Unknown Unknown Unknown", "Culex (Cux.) Vishnui subgroup", "empty")))

# reorder columns
f_species_pnr <- f_species_pnr[, c("House", "Temple", "Dumpsite", "Orchard", "Forest edge",  "Forest interior", "Bat cave")]

#iNEXT creates an iNEXT object of a list of df, q = 0, species richness, q = 1, shannon, q = 2 simpson
#diversity estimates with rarefied and extrapolated samples, full estimates, extrapolates to double sample size, unless specified

# richness pnr only
# pnr_iobj_rich <- iNEXT(f_species_pnr, q=0, datatype="abundance")
# pnr_rich_plot <- ggiNEXT(pnr_iobj_rich, type=1, se=TRUE, color.var="Assemblage", grey=FALSE)
# pnr_rich_plot <- pnr_rich_plot + labs(x = "Mosquito count", y = "Species richness", title = "") + my_theme + scale_shape_manual(values = c("House" = 20, "Temple" = 20, "Dumpsite" = 20, "Orchard" = 20, "Forest edge" = 20,  "Forest interior" = 20, "Bat cave" = 20)) + theme(legend.position = "right")
# pnr_rich_plot

# shannon pnr only
pnr_iobj_shan <- iNEXT(f_species_pnr, q=1, datatype="abundance", endpoint=750) 
pnr_shan_plot <- ggiNEXT(pnr_iobj_shan, type=1, se=TRUE, color.var="Assemblage", grey=FALSE) # makes classic iNEXT plot
pnr_shan_plot

#want to control design on plot
pnr_shan_data <- fortify(pnr_iobj_shan, type=1)
pnr_shan.point <- pnr_shan_data[which(pnr_shan_data$Method=="Observed"),]
pnr_shan.line <- pnr_shan_data[which(pnr_shan_data$Method!="Observed"),]
pnr_shan.line$Method <- as.factor(pnr_shan.line$Method)

pnr_shan <- ggplot(pnr_shan_data, aes(x=x, y=y, colour=Assemblage)) + 
  geom_point(aes(shape=Assemblage), size=4, data=pnr_shan.point, show.legend = F) +
  geom_line(aes(linetype=Method), lwd=1, data=pnr_shan.line) +
  geom_ribbon(aes(ymin=y.lwr, ymax=y.upr,
                  fill=Assemblage, colour=NULL), alpha=0.2) +
  labs(x="Number of individuals", y="Shannon diversity") +
  theme(legend.title=element_blank(),
        legend.position = "bottom", 
        text=element_text(size=18),
        legend.box = "vertical",
        axis.text = element_text(size=12, color="black"),
        axis.title = element_text(size=14, color="black"),
        legend.text=element_text(size=12),
        strip.background = element_rect(fill = "white", colour = "black"),
        panel.background = element_rect(fill = "aliceblue", color = "black"),
        panel.grid.major = element_line(color = "gray", linetype = "dotted"),
        panel.grid.minor = element_line(color = "gray", linetype = "dotted")) +
  scale_shape_manual(values = c(16,16,16,16,16,16,16)) +
  scale_linetype_manual(values = c("twodash", "solid")) 

# simpson pnr only
pnr_iobj_simp <- iNEXT(f_species_pnr, q=2, datatype="abundance", endpoint=750) 
pnr_simp_plot <- ggiNEXT(pnr_iobj_simp, type=1, se=TRUE, color.var="Assemblage", grey=FALSE)
pnr_simp_plot 

#want to control design on plot
pnr_simp_data <- fortify(pnr_iobj_simp, type=1)
pnr_simp.point <- pnr_simp_data[which(pnr_simp_data$Method=="Observed"),]
pnr_simp.line <- pnr_simp_data[which(pnr_simp_data$Method!="Observed"),]
pnr_simp.line$Method <- as.factor(pnr_simp.line$Method)

pnr_simp <- ggplot(pnr_simp_data, aes(x=x, y=y, colour=Assemblage)) + 
  geom_point(aes(shape=Assemblage), size=4, data=pnr_simp.point, show.legend = F) +
  geom_line(aes(linetype=Method), lwd=1, data=pnr_simp.line) +
  geom_ribbon(aes(ymin=y.lwr, ymax=y.upr,
                  fill=Assemblage, colour=NULL), alpha=0.2) +
  labs(x="Number of individuals", y="Simpson diversity") +
  # ylim(0,7) +
  theme(legend.title=element_blank(),
        legend.position = "none", 
        text=element_text(size=18),
        legend.box = "vertical",
        axis.text = element_text(size=12, color="black"),
        axis.title = element_text(size=14, color="black"),
        legend.text=element_text(size=12),
        strip.background = element_rect(fill = "white", colour = "black"),
        panel.background = element_rect(fill = "aliceblue", color = "black"),
        panel.grid.major = element_line(color = "gray", linetype = "dotted"),
        panel.grid.minor = element_line(color = "gray", linetype = "dotted")) +
  scale_shape_manual(values = c(16,16,16,16,16,16,16)) +
  scale_linetype_manual(values = c("twodash", "solid")) 

pnr_plots <- grid_arrange_shared_legend(pnr_shan, pnr_simp, ncol=2, nrow=1, position='right')
pnr_title <- ggdraw() + draw_label("Pong Nam Ron & Soi Dao", fontface = 'bold')
pnr_div_plot <- plot_grid(pnr_title, pnr_plots, ncol=1, rel_heights = c(0.1,1))
pnr_div_plot

ggsave(pnr_div_plot, file = "Plots\\pnr_diversity.jpg", width = 10, height = 4)

```  

```{r, mfl diversity}

#format df for iNEXT

f_species_mfl <- mq_females %>% 
  filter(district == "Mae Fah Luang") %>%
  dplyr::select(closest_id, interface, females) %>%
  group_by(closest_id, interface) %>%
  summarise(sum = sum(females)) %>% ungroup() %>%
  pivot_wider(names_from = "interface", 
              values_from = "sum") %>%
  mutate_at(c(2:6), ~replace(., is.na(.), 0)) %>%
  as.data.frame()

rownames(f_species_mfl) <- f_species_mfl$closest_id
f_species_mfl <- f_species_mfl[,-1]

# remove unidentified sp (Arm. aureolineatus - this was male)
f_species_mfl <- f_species_mfl %>%
  filter(!(rownames(.) %in% c("Aedes  spp.", "Aedes (Ste.) spp.", "Armigeres  spp.", "Armigeres (Arm.) spp.", "Armigeres (Lei.) spp.", "Coquillettidia (Coq.) spp.",  "Culex  (Cux.) spp.", "Culex  spp.", "Culex (Cui.) spp.", "Culex (Cux.) spp.", "Culex (Eum.) spp.", "Mansonia  spp.",  "Tripteroides  spp.", "Tripteroides (Rah.) spp.", "Armigeres (Lei.) aureolineatus", "empty")))

# reorder columns
f_species_mfl <- f_species_mfl[, c("House", "Dumpsite", "Orchard", "Forest edge",  "Forest interior")]

# shannon mfl only
mfl_iobj_shan <- iNEXT(f_species_mfl, q=1, datatype="abundance", endpoint=750) 
# mfl_shan_plot <- ggiNEXT(mfl_iobj_shan, type=1, se=TRUE, color.var="Assemblage", grey=FALSE) # makes classic iNEXT plot
# mfl_shan_plot

#want to control design on plot
mfl_shan_data <- fortify(mfl_iobj_shan, type=1)
mfl_shan.point <- mfl_shan_data[which(mfl_shan_data$Method=="Observed"),]
mfl_shan.line <- mfl_shan_data[which(mfl_shan_data$Method!="Observed"),]
mfl_shan.line$Method <- as.factor(mfl_shan.line$Method)

mfl_shan <- ggplot(mfl_shan_data, aes(x=x, y=y, colour=Assemblage)) + 
  geom_point(aes(shape=Assemblage), size=4, data=mfl_shan.point, show.legend = F) +
  geom_line(aes(linetype=Method), lwd=1, data=mfl_shan.line) +
  geom_ribbon(aes(ymin=y.lwr, ymax=y.upr,
                  fill=Assemblage, colour=NULL), alpha=0.2) +
  labs(x="Number of individuals", y="Shannon diversity") +
  theme(legend.title=element_blank(),
        legend.position = "bottom", 
        text=element_text(size=18),
        legend.box = "vertical",
        axis.text = element_text(size=12, color="black"),
        axis.title = element_text(size=14, color="black"),
        legend.text=element_text(size=12),
        strip.background = element_rect(fill = "white", colour = "black"),
        panel.background = element_rect(fill = "aliceblue", color = "black"),
        panel.grid.major = element_line(color = "gray", linetype = "dotted"),
        panel.grid.minor = element_line(color = "gray", linetype = "dotted")) +
  scale_shape_manual(values = c(16,16,16,16,16,16,16)) +
  scale_linetype_manual(values = c("twodash", "solid")) 

# simpson mfl only
mfl_iobj_simp <- iNEXT(f_species_mfl, q=2, datatype="abundance", endpoint=750) 
mfl_simp_plot <- ggiNEXT(mfl_iobj_simp, type=1, se=TRUE, color.var="Assemblage", grey=FALSE)
mfl_simp_plot 

#want to control design on plot
mfl_simp_data <- fortify(mfl_iobj_simp, type=1)
mfl_simp.point <- mfl_simp_data[which(mfl_simp_data$Method=="Observed"),]
mfl_simp.line <- mfl_simp_data[which(mfl_simp_data$Method!="Observed"),]
mfl_simp.line$Method <- as.factor(mfl_simp.line$Method)

mfl_simp <- ggplot(mfl_simp_data, aes(x=x, y=y, colour=Assemblage)) + 
  geom_point(aes(shape=Assemblage), size=4, data=mfl_simp.point, show.legend = F) +
  geom_line(aes(linetype=Method), lwd=1, data=mfl_simp.line) +
  geom_ribbon(aes(ymin=y.lwr, ymax=y.upr,
                  fill=Assemblage, colour=NULL), alpha=0.2) +
  labs(x="Number of individuals", y="Simpson diversity") +
  # ylim(0,7) +
  theme(legend.title=element_blank(),
        legend.position = "none", 
        text=element_text(size=18),
        legend.box = "vertical",
        axis.text = element_text(size=12, color="black"),
        axis.title = element_text(size=14, color="black"),
        legend.text=element_text(size=12),
        strip.background = element_rect(fill = "white", colour = "black"),
        panel.background = element_rect(fill = "aliceblue", color = "black"),
        panel.grid.major = element_line(color = "gray", linetype = "dotted"),
        panel.grid.minor = element_line(color = "gray", linetype = "dotted")) +
  scale_shape_manual(values = c(16,16,16,16,16,16,16)) +
  scale_linetype_manual(values = c("twodash", "solid")) 

mfl_plots <- grid_arrange_shared_legend(mfl_shan, mfl_simp, ncol=2, nrow=1, position='right')
mfl_title <- ggdraw() + draw_label("Mae Fah Luang", fontface = 'bold')
mfl_div_plot <- plot_grid(mfl_title, mfl_plots, ncol=1, rel_heights = c(0.1,1))
mfl_div_plot

ggsave(mfl_div_plot, file = "Plots\\mfl_diversity.jpg", width = 10, height = 4)


```

```{r, wk diversity}

#format df for iNEXT

f_species_wk <- mq_females %>% 
  filter(district == "Wiang Kaen") %>%
  dplyr::select(closest_id, interface, females) %>%
  group_by(closest_id, interface) %>%
  summarise(sum = sum(females)) %>% ungroup() %>%
  pivot_wider(names_from = "interface", 
              values_from = "sum") %>%
  mutate_at(c(2:6), ~replace(., is.na(.), 0)) %>%
  as.data.frame()

rownames(f_species_wk) <- f_species_wk$closest_id
f_species_wk <- f_species_wk[,-1]

# remove unidentified sp
f_species_wk <- f_species_wk %>%
  filter(!(rownames(.) %in% c("Aedes  spp.", "Anopheles (Cel.) spp.", "Armigeres  spp.", "Armigeres (Arm.) spp.", "Culex (Cux.) spp.", "Tripteroides (Rah.) spp.", "Tripteroides (Trp.) spp.", "Unknown  Unknown", "empty")))

# reorder columns
f_species_wk <- f_species_wk[, c("House", "Orchard", "Forest edge",  "Forest interior", "Bat cave")]

# shannon wk only
wk_iobj_shan <- iNEXT(f_species_wk, q=1, datatype="abundance", endpoint=750) 
# wk_shan_plot <- ggiNEXT(wk_iobj_shan, type=1, se=TRUE, color.var="Assemblage", grey=FALSE) # makes classic iNEXT plot
# wk_shan_plot

#want to control design on plot
wk_shan_data <- fortify(wk_iobj_shan, type=1)
wk_shan.point <- wk_shan_data[which(wk_shan_data$Method=="Observed"),]
wk_shan.line <- wk_shan_data[which(wk_shan_data$Method!="Observed"),]
wk_shan.line$Method <- as.factor(wk_shan.line$Method)

wk_shan <- ggplot(wk_shan_data, aes(x=x, y=y, colour=Assemblage)) + 
  geom_point(aes(shape=Assemblage), size=4, data=wk_shan.point, show.legend = F) +
  geom_line(aes(linetype=Method), lwd=1, data=wk_shan.line) +
  geom_ribbon(aes(ymin=y.lwr, ymax=y.upr,
                  fill=Assemblage, colour=NULL), alpha=0.2) +
  labs(x="Number of individuals", y="Shannon diversity") +
  theme(legend.title=element_blank(),
        legend.position = "bottom", 
        text=element_text(size=18),
        legend.box = "vertical",
        axis.text = element_text(size=12, color="black"),
        axis.title = element_text(size=14, color="black"),
        legend.text=element_text(size=12),
        strip.background = element_rect(fill = "white", colour = "black"),
        panel.background = element_rect(fill = "aliceblue", color = "black"),
        panel.grid.major = element_line(color = "gray", linetype = "dotted"),
        panel.grid.minor = element_line(color = "gray", linetype = "dotted")) +
  scale_shape_manual(values = c(16,16,16,16,16,16,16)) +
  scale_linetype_manual(values = c("twodash", "solid")) 

# simpson wk only
wk_iobj_simp <- iNEXT(f_species_wk, q=2, datatype="abundance", endpoint=750) 
wk_simp_plot <- ggiNEXT(wk_iobj_simp, type=1, se=TRUE, color.var="Assemblage", grey=FALSE)
wk_simp_plot 

#want to control design on plot
wk_simp_data <- fortify(wk_iobj_simp, type=1)
wk_simp.point <- wk_simp_data[which(wk_simp_data$Method=="Observed"),]
wk_simp.line <- wk_simp_data[which(wk_simp_data$Method!="Observed"),]
wk_simp.line$Method <- as.factor(wk_simp.line$Method)

wk_simp <- ggplot(wk_simp_data, aes(x=x, y=y, colour=Assemblage)) + 
  geom_point(aes(shape=Assemblage), size=4, data=wk_simp.point, show.legend = F) +
  geom_line(aes(linetype=Method), lwd=1, data=wk_simp.line) +
  geom_ribbon(aes(ymin=y.lwr, ymax=y.upr,
                  fill=Assemblage, colour=NULL), alpha=0.2) +
  labs(x="Number of individuals", y="Simpson diversity") +
  ylim(0,7) +
  theme(legend.title=element_blank(),
        legend.position = "none", 
        text=element_text(size=18),
        legend.box = "vertical",
        axis.text = element_text(size=12, color="black"),
        axis.title = element_text(size=14, color="black"),
        legend.text=element_text(size=12),
        strip.background = element_rect(fill = "white", colour = "black"),
        panel.background = element_rect(fill = "aliceblue", color = "black"),
        panel.grid.major = element_line(color = "gray", linetype = "dotted"),
        panel.grid.minor = element_line(color = "gray", linetype = "dotted")) +
  scale_shape_manual(values = c(16,16,16,16,16,16,16)) +
  scale_linetype_manual(values = c("twodash", "solid")) 

wk_plots <- grid_arrange_shared_legend(wk_shan, wk_simp, ncol=2, nrow=1, position='right')
wk_title <- ggdraw() + draw_label("Wiang Kaen", fontface = 'bold')
wk_div_plot <- plot_grid(wk_title, wk_plots, ncol=1, rel_heights = c(0.1,1))

ggsave(wk_div_plot, file = "Plots\\wk_diversity.jpg", width = 10, height = 4)

all_title <- ggdraw() + draw_label("Mosquito species diversity indices in study interfaces", size=16)
div_all_districts <- plot_grid(all_title, pnr_div_plot, mfl_div_plot, wk_div_plot, nrow=4, rel_heights = c(0.1,1,1,1))
div_all_districts

ggsave(div_all_districts, file = "Plots\\diversity_all_districts.jpg", width = 14, height = 8)

```


```{r, shannon diversity}

#for shannon plot
#bring all data together
#pnr
pnr_shan_data <- fortify(pnr_iobj_shan, type=1)
pnr_shan_data$district <- "Pong Nam Ron & Soi Dao"
pnr_shan.point <- pnr_shan_data[which(pnr_shan_data$Method=="Observed"),]
pnr_shan.point$district <- "Pong Nam Ron & Soi Dao"
pnr_shan.line <- pnr_shan_data[which(pnr_shan_data$Method!="Observed"),]
pnr_shan.line$Method <- as.factor(pnr_shan.line$Method)
pnr_shan.line$district <- "Pong Nam Ron & Soi Dao"

#mfl
mfl_shan_data <- fortify(mfl_iobj_shan, type=1)
mfl_shan_data$district <- "Mae Fah Luang"
mfl_shan.point <- mfl_shan_data[which(mfl_shan_data$Method=="Observed"),]
mfl_shan.point$district <- "Mae Fah Luang"
mfl_shan.line <- mfl_shan_data[which(mfl_shan_data$Method!="Observed"),]
mfl_shan.line$Method <- as.factor(mfl_shan.line$Method)
mfl_shan.line$district <- "Mae Fah Luang"

#wk
wk_shan_data <- fortify(wk_iobj_shan, type=1)
wk_shan_data$district <- "Wiang Kaen" 
wk_shan.point <- wk_shan_data[which(wk_shan_data$Method=="Observed"),]
wk_shan.point$district <- "Wiang Kaen"  
wk_shan.line <- wk_shan_data[which(wk_shan_data$Method!="Observed"),]
wk_shan.line$Method <- as.factor(wk_shan.line$Method)
wk_shan.line$district <- "Wiang Kaen" 

#combine districts in one df
shan_data <- rbind(pnr_shan_data, mfl_shan_data, wk_shan_data)
shan.point <- rbind(pnr_shan.point, mfl_shan.point, wk_shan.point)
shan.line <- rbind(pnr_shan.line, mfl_shan.line, wk_shan.line)

#set levels
shan_data$Assemblage <- factor(shan_data$Assemblage, levels = c("House", "Temple", "Dumpsite", "Orchard", "Forest edge",  "Forest interior", "Bat cave"))
shan_data$district <- factor(shan_data$district, levels = c("Pong Nam Ron & Soi Dao", "Mae Fah Luang", "Wiang Kaen")) 
shan.point$Assemblage <- factor(shan.point$Assemblage, levels = c("House", "Temple", "Dumpsite", "Orchard", "Forest edge",  "Forest interior", "Bat cave"))
shan.point$district <- factor(shan.point$district, levels = c("Pong Nam Ron & Soi Dao", "Mae Fah Luang", "Wiang Kaen")) 
shan.line$Assemblage <- factor(shan.line$Assemblage, levels = c("House", "Temple", "Dumpsite", "Orchard", "Forest edge",  "Forest interior", "Bat cave"))
shan.line$district <- factor(shan.line$district, levels = c("Pong Nam Ron & Soi Dao", "Mae Fah Luang", "Wiang Kaen")) 


```

```{r, simpson diversity}

#for simpson plot
#bring all data together
pnr_simp_data <- fortify(pnr_iobj_simp, type=1)
pnr_simp_data$district <- "Pong Nam Ron & Soi Dao"
pnr_simp.point <- pnr_simp_data[which(pnr_simp_data$Method=="Observed"),]
pnr_simp.point$district <- "Pong Nam Ron & Soi Dao"
pnr_simp.line <- pnr_simp_data[which(pnr_simp_data$Method!="Observed"),]
pnr_simp.line$Method <- as.factor(pnr_simp.line$Method)
pnr_simp.line$district <- "Pong Nam Ron & Soi Dao" 

#mfl
mfl_simp_data <- fortify(mfl_iobj_simp, type=1)
mfl_simp_data$district <- "Mae Fah Luang"
mfl_simp.point <- mfl_simp_data[which(mfl_simp_data$Method=="Observed"),]
mfl_simp.point$district <- "Mae Fah Luang"
mfl_simp.line <- mfl_simp_data[which(mfl_simp_data$Method!="Observed"),]
mfl_simp.line$Method <- as.factor(mfl_simp.line$Method)
mfl_simp.line$district <- "Mae Fah Luang"

#wk
wk_simp_data <- fortify(wk_iobj_simp, type=1)
wk_simp_data$district <- "Wiang Kaen" 
wk_simp.point <- wk_simp_data[which(wk_simp_data$Method=="Observed"),]
wk_simp.point$district <- "Wiang Kaen" 
wk_simp.line <- wk_simp_data[which(wk_simp_data$Method!="Observed"),]
wk_simp.line$Method <- as.factor(wk_simp.line$Method)
wk_simp.line$district <- "Wiang Kaen" 


#combine districts in one df
simp_data <- rbind(pnr_simp_data, mfl_simp_data, wk_simp_data)
simp.point <- rbind(pnr_simp.point, mfl_simp.point, wk_simp.point)
simp.line <- rbind(pnr_simp.line, mfl_simp.line, wk_simp.line)

#set levels
simp_data$Assemblage <- factor(simp_data$Assemblage, levels = c("House", "Temple", "Dumpsite", "Orchard", "Forest edge",  "Forest interior", "Bat cave"))
simp_data$district <- factor(simp_data$district, levels = c("Pong Nam Ron & Soi Dao", "Mae Fah Luang", "Wiang Kaen")) 
simp.point$Assemblage <- factor(simp.point$Assemblage, levels = c("House", "Temple", "Dumpsite", "Orchard", "Forest edge",  "Forest interior", "Bat cave"))
simp.point$district <- factor(simp.point$district, levels = c("Pong Nam Ron & Soi Dao", "Mae Fah Luang", "Wiang Kaen")) 
simp.line$Assemblage <- factor(simp.line$Assemblage, levels = c("House", "Temple", "Dumpsite", "Orchard", "Forest edge",  "Forest interior", "Bat cave"))
simp.line$district <- factor(simp.line$district, levels = c("Pong Nam Ron & Soi Dao", "Mae Fah Luang", "Wiang Kaen")) 

```

```{r, diversity combined}

#try out combined plot with district~diversity facet_grid

shan_data$diversity <- "Shannon"
simp_data$diversity <- "Simpson"
shan.point$diversity <- "Shannon"
simp.point$diversity <- "Simpson"
shan.line$diversity <- "Shannon"
simp.line$diversity <- "Simpson"

div_data <- rbind(shan_data, simp_data)
div_data <- div_data %>% dplyr::rename("Interface" = "Assemblage")
point_data <- rbind(shan.point, simp.point)
point_data <- point_data %>% dplyr::rename("Interface" = "Assemblage")
line_data <- rbind(shan.line, simp.line)
line_data <- line_data %>% dplyr::rename("Interface" = "Assemblage")

div_plot <- ggplot(div_data, aes(x=x, y=y, colour=Interface)) + 
  geom_point(aes(shape=Interface), size=4, data=point_data, show.legend = F) +
  geom_line(aes(linetype=Method), lwd=1, data=line_data) +
  geom_ribbon(aes(ymin=y.lwr, ymax=y.upr,
                  fill=Interface, colour=NULL), alpha=0.2, show.legend = F) +
  labs(x="Number of mosquitoes", y="Diversity", legend = "test") +
  facet_grid(diversity~district) +
  theme(legend.title=element_blank(),
        legend.position = "right", 
        text=element_text(size=18),
        legend.box = "vertical",
        axis.text = element_text(size=12, color="black"),
        axis.title = element_text(size=14, color="black"),
        legend.text=element_text(size=12),
        strip.background = element_rect(fill = "white", colour = "black"),
        panel.background = element_rect(fill = "aliceblue", color = "black"),
        panel.grid.major = element_line(color = "gray", linetype = "dotted"),
        panel.grid.minor = element_line(color = "gray", linetype = "dotted")) +
  scale_shape_manual(values = c(16,16,16,16,16,16,16)) +
  scale_linetype_manual(values = c("twodash", "solid")) +
  scale_y_continuous(breaks = c(0,2,4,6,8,10,12), limits = c(0,13)) +
  my_theme

div_plot

ggsave(div_plot, file = "Plots\\diversity_indices.jpg", width = 14, height = 8)

```