#Loading necessary packages
library(MASS)
library(ggplot2)
library(ggrepel)

#Loading necessart data
Mort_Inf_NATURE_ID50_LD50 <- read.csv("Mort_Inf_NATURE_ID50_LD50.csv")


#Creating dataframes for each species
New<-split(Mort_Inf_NATURE_ID50_LD50,Mort_Inf_NATURE_ID50_LD50$Species)
Y <- lapply(seq_along(New), function(x) as.data.frame(New[[x]])[, 1:10]) 
Y

#Creating New dataframes
AMME <- Y[[1]]
AMOP <- Y[[2]]
ANAE <- Y[[3]]
AQCE <- Y[[4]]
CHPE <- Y[[5]]
CRAL<- Y[[6]]
DEOC<- Y[[7]]
ENKL<- Y[[8]]
ENXA<- Y[[9]]
EUBI<- Y[[10]]
EULU<- Y[[11]]
EUWI<- Y[[12]]
HYCH<- Y[[13]]
LICH<- Y[[14]]
NOME<- Y[[15]]
NOPE<- Y[[16]]
OSSE<- Y[[17]]
PLME<- Y[[18]]
PLSH<- Y[[19]]
PSRU<- Y[[20]]
PSST<- Y[[21]]
SCHO<- Y[[22]]
SILA<- Y[[23]]
TAGR<- Y[[24]]
#AMME
AMMEprop<-AMME$Infected/(AMME$Total)
AMMEprop
AMMEtot<-AMME$Total
AMMEtot
AMMEmodel<-glm(AMMEprop~log(Treatment),family=binomial(link="probit"),data=AMME,weight=AMMEtot)
dose.p(AMMEmodel,p=0.5)
exp(17.43764)
exp(4.868299)


####ID-50 Estimations
#AMOP
AMOPprop<-AMOP$Infected/(AMOP$Total)
AMOPprop
AMOPtot<-AMOP$Total
AMOPtot
AMOPmodel<-glm(AMOPprop~log(Treatment),family=binomial(link="probit"),data=AMOP,weight=AMOPtot)
dose.p(AMOPmodel,p=0.5)
exp(15.9707)
exp(857.8431)

#ANAE
ANAEprop<-ANAE$Infected/(ANAE$Total)
ANAEprop
ANAEtot<-ANAE$Total
ANAEtot
ANAEmodel<-glm(ANAEprop~log(Treatment),family=binomial(link="probit"),data=ANAE,weight=ANAEtot)
dose.p(ANAEmodel,p=0.5)
exp(9.162012)
exp(0.7165214)
#AQCE
AQCEprop<-AQCE$Infected/(AQCE$Total)
AQCEprop
AQCEtot<-AQCE$Total
AQCEtot
AQCEmodel<-glm(AQCEprop~log(Treatment),family=binomial(link="probit"),data=AQCE,weight=AQCEtot)
dose.p(AQCEmodel,p=0.5)
exp(11.07163)
exp(0.6336817)

#CHPE
CHPEprop<-CHPE$Infected/(CHPE$Total)
CHPEprop
CHPEtot<-CHPE$Total
CHPEtot
CHPEmodel<-glm(CHPEprop~log(Treatment),family=binomial(link="probit"),data=CHPE,weight=CHPEtot)
dose.p(CHPEmodel,p=0.5)
exp(9.707659)
exp(1.245197)
#CRAL
CRALprop<-CRAL$Infected/(CRAL$Total)
CRALprop
CRALtot<-CRAL$Total
CRALtot
CRALmodel<-glm(CRALprop~log(Treatment),family=binomial(link="probit"),data=CRAL,weight=CRALtot)
dose.p(CRALmodel,p=0.5)
exp(-143.6688)

#DEOC
DEOCprop<-DEOC$Infected/(DEOC$Total)
DEOCprop
DEOCtot<-DEOC$Total
DEOCtot
DEOCmodel<-glm(DEOCprop~log(Treatment),family=binomial(link="probit"),data=DEOC,weight=DEOCtot)
dose.p(DEOCmodel,p=0.5)
exp(16.57926)
exp(1.405215)
#EUBI
EUBIprop<-EUBI$Infected/(EUBI$Total)
EUBIprop
EUBItot<-EUBI$Total
EUBItot
EUBImodel<-glm(EUBIprop~log(Treatment),family=binomial(link="probit"),data=EUBI,weight=EUBItot)
dose.p(EUBImodel,p=0.5)
exp(8.001905)
exp(810.0721)
#EULU
EULUprop<-EULU$Infected/(EULU$Total)
EULUprop
EULUtot<-EULU$Total
EULUtot
EULUmodel<-glm(EULUprop~log(Treatment),family=binomial(link="probit"),data=EULU,weight=EULUtot)
dose.p(EULUmodel,p=0.5)
exp(11.56568)
exp( 0.934861)
#LICH
LICHprop<-LICH$Infected/(LICH$Total)
LICHprop
LICHtot<-LICH$Total
LICHtot
LICHmodel<-glm(LICHprop~log(Treatment),family=binomial(link="probit"),data=LICH,weight=LICHtot)
dose.p(LICHmodel,p=0.5)
exp(15.31872)
exp(105.5147)

#EUWI
EUWIprop<-EUWI$Infected/(EUWI$Total)
EUWIprop
EUWItot<-EUWI$Total
EUWItot
EUWImodel<-glm(EUWIprop~log(Treatment),family=binomial(link="probit"),data=EUWI,weight=EUWItot)
dose.p(EUWImodel,p=0.5)
exp(0.3445269)

#HYCH
HYCHprop<-HYCH$Infected/(HYCH$Total)
HYCHprop
HYCHtot<-HYCH$Total
HYCHtot
HYCHmodel<-glm(HYCHprop~log(Treatment),family=binomial(link="probit"),data=HYCH,weight=HYCHtot)
dose.p(HYCHmodel,p=0.5)
exp(15.42495)
exp(0.1378277)
#NOME
NOMEprop<-NOME$Infected/(NOME$Total)
NOMEprop
NOMEtot<-NOME$Total
NOMEtot
NOMEmodel<-glm(NOMEprop~log(Treatment),family=binomial(link="probit"),data=NOME,weight=NOMEtot)
dose.p(NOMEmodel,p=0.5)
exp(8.218279)

#NOPE
NOPEprop<-NOPE$Infected/(NOPE$Total)
NOPEprop
NOPEtot<-NOPE$Total
NOPEtot
NOPEmodel<-glm(NOPEprop~log(Treatment),family=binomial(link="probit"),data=NOPE,weight=NOPEtot)
dose.p(NOPEmodel,p=0.5)


#OSSE
OSSEprop<-OSSE$Infected/(OSSE$Total)
OSSEprop
OSSEtot<-OSSE$Total
OSSEtot
OSSEmodel<-glm(OSSEprop~log(Treatment),family=binomial(link="probit"),data=OSSE,weight=OSSEtot)
dose.p(OSSEmodel,p=0.5)
exp(8.73286)

#PLME
PLMEprop<-PLME$Infected/(PLME$Total)
PLMEprop
PLMEtot<-PLME$Total
PLMEtot
PLMEmodel<-glm(PLMEprop~log(Treatment),family=binomial(link="probit"),data=PLME,weight=PLMEtot)
dose.p(PLMEmodel,p=0.5)
exp(17.7582)
exp(2.437811)
#PLSH
PLSHprop<-PLSH$Infected/(PLSH$Total)
PLSHprop
PLSHtot<-PLSH$Total
PLSHtot
PLSHmodel<-glm(PLSHprop~log(Treatment),family=binomial(link="probit"),data=PLSH,weight=PLSHtot)
dose.p(PLSHmodel,p=0.5)
exp(15.48947)
exp(77.52465)
#PSRU
PSRUprop<-PSRU$Infected/(PSRU$Total)
PSRUprop
PSRUtot<-PSRU$Total
PSRUtot
PSRUmodel<-glm(PSRUprop~log(Treatment),family=binomial(link="probit"),data=PSRU,weight=PSRUtot)
dose.p(PSRUmodel,p=0.5)
exp(9.652154)
exp(0.5912302)
#PSST
PSSTprop<-PSST$Infected/(PSST$Total)
PSSTprop
PSSTtot<-PSST$Total
PSSTtot
PSSTmodel<-glm(PSSTprop~log(Treatment),family=binomial(link="probit"),data=PSST,weight=PSSTtot)
dose.p(PSSTmodel,p=0.5)
exp(13.8674)
exp(0.6557128)
#SCHO
SCHOprop<-SCHO$Infected/(SCHO$Total)
SCHOprop
SCHOtot<-SCHO$Total
SCHOtot
SCHOmodel<-glm(SCHOprop~log(Treatment),family=binomial(link="probit"),data=SCHO,weight=SCHOtot)
dose.p(SCHOmodel,p=0.5)
exp(12.9297)
exp(1.039396)
#SILA
SILAprop<-SILA$Infected/(SILA$Total)
SILAprop
SILAtot<-SILA$Total
SILAtot
SILAmodel<-glm(SILAprop~log(Treatment),family=binomial(link="probit"),data=SILA,weight=SILAtot)
dose.p(SILAmodel,p=0.5)
exp(17.78375)
exp(3.564534)
#TAGR
TAGRprop<-TAGR$Infected/(TAGR$Total)
TAGRprop
TAGRtot<-TAGR$Total
TAGRtot
TAGRmodel<-glm(TAGRprop~log(Treatment),family=binomial(link="probit"),data=TAGR,weight=TAGRtot)
dose.p(TAGRmodel,p=0.5)
exp(8.686874)

#ENKL
ENKLprop<-ENKL$Infected/(ENKL$Total)
ENKLprop
ENKLtot<-ENKL$Total
ENKLtot
ENKLmodel<-glm(ENKLprop~log(Treatment),family=binomial(link="probit"),data=ENKL,weight=ENKLtot)
dose.p(ENKLmodel,p=0.5)
exp(9.675397)

#ENXA
ENXAprop<-ENXA$Infected/(ENXA$Total)
ENXAprop
ENXAtot<-ENXA$Total
ENXAtot
ENXAmodel<-glm(ENXAprop~log(Treatment),family=binomial(link="probit"),data=ENXA,weight=ENXAtot)
dose.p(ENXAmodel,p=0.5)
exp( 12.43817)
exp( 1.255821)


#LD-50 Estimations
#ANAE
ANAEpropDead<-ANAE$Dead/(ANAE$Total)
ANAEpropDead
ANAEmodelDead<-glm(ANAEpropDead~log(Treatment),family=binomial(link="probit"),data=ANAE,weight=ANAEtot)
dose.p(ANAEmodelDead,p=0.5)
exp(14.15861)
exp(2.611782)

#AQCE
AQCEpropDead<-AQCE$Dead/(AQCE$Total)
AQCEmodelDead<-glm(AQCEpropDead~log(Treatment),family=binomial(link="probit"),data=AQCE,weight=AQCEtot)
dose.p(AQCEmodelDead,p=0.5)
exp(11.07163)
exp(0.6336817)
#
CHPEpropDead<-CHPE$Dead/(CHPE$Total)
CHPEmodelDead<-glm(CHPEpropDead~log(Treatment),family=binomial(link="probit"),data=CHPE,weight=CHPEtot)
dose.p(CHPEmodelDead,p=0.5)
exp(11.53898)
exp(0.6258333)
#EUBI
EUBIpropDead<-EUBI$Dead/(EUBI$Total)
EUBImodelDead<-glm(EUBIpropDead~log(Treatment),family=binomial(link="probit"),data=EUBI,weight=EUBItot)
dose.p(EUBImodelDead,p=0.5)
exp(14.86623)
exp(1.099715)

#EUWI
EUWIpropDead<-EUWI$Dead/(EUWI$Total)
EUWImodelDead<-glm(EUWIpropDead~log(Treatment),family=binomial(link="probit"),data=EUWI,weight=EUWItot)
dose.p(EUWImodelDead,p=0.5)
exp(14.91199)
exp(0.3357429)
#NOME
NOMEpropDead<-NOME$Dead/(NOME$Total)
NOMEmodelDead<-glm(NOMEpropDead~log(Treatment),family=binomial(link="probit"),data=NOME,weight=NOMEtot)
dose.p(NOMEmodelDead,p=0.5)
exp(9.921735)
exp( 0.4014456)
#NOPE
NOPEpropDead<-NOPE$Dead/(NOPE$Total)
NOPEpropDead
NOPEmodelDead<-glm(NOPEpropDead~log(Treatment),family=binomial(link="probit"),data=NOPE,weight=NOPEtot)
dose.p(NOPEmodelDead,p=0.5)
exp(11.18118)
exp(0.6144835)
#OSSE
OSSEpropDead<-OSSE$Dead/(OSSE$Total)

OSSEmodelDead<-glm(OSSEpropDead~log(Treatment),family=binomial(link="probit"),data=OSSE,weight=OSSEtot)
dose.p(OSSEmodelDead,p=0.5)
exp(12.90576)

#PSRU
PSRUpropDead<-PSRU$Dead/(PSRU$Total)
PSRUmodelDead<-glm(PSRUpropDead~log(Treatment),family=binomial(link="probit"),data=PSRU,weight=PSRUtot)
dose.p(PSRUmodelDead,p=0.5)
exp(14.26702)

#SCHO
SCHOpropDead<-SCHO$Dead/(SCHO$Total)
SCHOmodelDead<-glm(SCHOpropDead~log(Treatment),family=binomial(link="probit"),data=SCHO,weight=SCHOtot)
dose.p(SCHOmodelDead,p=0.5)
exp(15.7704)

#TAGR
TAGRpropDead<-TAGR$Dead/(TAGR$Total)
TAGRmodelDead<-glm(TAGRpropDead~log(Treatment),family=binomial(link="probit"),data=TAGR,weight=TAGRtot)
dose.p(TAGRmodelDead,p=0.5)
exp(15.39697)
exp(0.6518898)
#ENKL
ENKLpropDead<-ENKL$Dead/(ENKL$Total)
ENKLmodelDead<-glm(ENKLpropDead~log(Treatment),family=binomial(link="probit"),data=ENKL,weight=ENKLtot)
dose.p(ENKLmodelDead,p=0.5)
exp(13.789)
exp(0.7984993)
#ENXA
ENXApropDead<-ENXA$Dead/(ENXA$Total)
ENXAmodelDead<-glm(ENXApropDead~log(Treatment),family=binomial(link="probit"),data=ENXA,weight=ENXAtot)
dose.p(ENXAmodelDead,p=0.5)
exp(15.58563)
#NOVI
NOVIdata = data.frame(Treatment = c(5000000,500000,50000,5000),
                      alive = c(0,0,0,.1),
                      dead = c(1,1,1,.9),NOVItot=c(10,10,10,10))
NOVImodelDead<-glm(dead~log(Treatment),family=binomial(link="probit"),data=NOVIdata,weight=NOVItot)
dose.p(NOVImodelDead)
exp(7.97144)




#Creating ID50/LD50 Ordination Figure 3
ID50_LD50_Ordination_Data.csv <- read.csv("ID50_LD50_Ordination_Data.csv", row.names=1)
ID50_LD50_Ordination_Data.csv 

#Figure3
OrdinationPlot<-ggplot(ID50_LD50_Ordination_Data.csv, aes(x=log(ID50), y=log(LD50),color=ConservationRisk,fill=ConservationRisk,size=Amplification_Potential))+
  geom_label_repel(label=rownames(ID50_LD50_Ordination_Data.csv), 
                   nudge_x = 0.25, nudge_y = -0.25,color="black",fill="white",size=3.5)+ 
  geom_hline(yintercept=16.25,color="black",size=1,lty=2)+
  scale_x_continuous(expand=c(0.01,0.01)) + 
  scale_y_continuous(expand=c(0.01,0.01))+
  theme(axis.line.x=element_line(),
        axis.line.y=element_line())+geom_point(pch=21)+scale_size_continuous(range=c(2,8))+
  coord_cartesian(xlim=c(7.5,20),ylim=c(7.5,17.5))+
  scale_color_manual(values=c("black","black","black","black"))+
  scale_fill_manual(values=c("darkorange2","skyblue2","yellow2","firebrick4"))+
  theme(legend.position="none",axis.title.x = element_blank(),axis.title.y=element_blank(),
        axis.text.x=element_text(face="bold",color="black"),
        axis.text.y=element_text(face="bold",color="black"))
OrdinationPlot