# Loaded package ---------------------------------------------------------
library(tidyverse)

#the line below removes everything from the R environment.
#This ensures there is no other data etc that could interfere with the analysis. 
rm(list=ls())

# Fungi virulence analysis ------------------------------------------------

## Data importation -------------------------------------------------------
mortality<-read.csv("Mortality.csv", 
                    header = T, 
                    stringsAsFactors = T)

## Data set exploration -----------------------------------------------------
mortality
str(mortality)
levels(mortality$Treatment)

# Let's get mosquitoes mortality regarding the fungi strains --------------
prop.table<- aggregate(mortality$Death, by=list(mortality$Replicate,mortality$Treatment), FUN='sum')
colnames(prop.table)<- c('Replicate','Treatment','Dead')
prop.table$prop.dead<- prop.table$Dead/25
prop.table$Alive<- 25-prop.table$Dead
prop.table$prop.alive<- (25-prop.table$Dead)/25
prop.table

# Plot the mortality data -------------------------------------------------

ggplot(data = prop.table, aes(x = Treatment, y = prop.dead, fill=Treatment)) +
  geom_boxplot() +
  labs(y = "Mortality",
       x = "Treatments",
       title = "Fig: Native Metarhizium strains virulence against insecticide resistant malaria vector Anopheles coluzzii VKPR") +
  theme_classic() +
  theme(strip.text = element_text(size = rel(1)),
        panel.grid = element_blank())

# Test efficacy of fungi strains (Treatment) on mosquitoes mortality ------------------
m0<-glm(cbind(Dead,Alive) ~ Treatment, data = prop.table,
        family = binomial(link = "logit"))
m1<-glm(cbind(Dead,Alive) ~ 1, data = prop.table,
        family = binomial(link = "logit"))

anova(m0,m1, test = "LRT")
#the model m0 with treatment is better so we keep that one

summary(m0)
#All treatments significantly increase mortality of mosquitoes compared to control (about 2 log odds)
exp(m0$coefficients)
#Mosquitoes have 7.9, 12.2 and 16.6 times the odds of dying with S10, S26 and S31 respectively, compared to control) .