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