#Loading required packages
library(xlsx)
library(picante)
library(ggtree)

#Loading needed data files
Nature_Paper_Infection_Mort_Perc <- read.csv("Nature_Paper_Infection_Mortality_Perc.csv", row.names=1)
CaudataTree<-read.tree("SalamanderSpecies.nwk")

#Creating new species list that only includes salamanders
Nature_Paper_Infection_Mort_Perc<-subset(Nature_Paper_Infection_Mort_Perc,Order=Caudata)

#Reformatting data 
InfectionTrait<-Nature_Paper_Infection_Mort_Perc$Infection
MortTrait<-Nature_Paper_Infection_Mort_Perc$Mortality
names(InfectionTrait)<-Nature_Paper_Infection_Mort_Perc$Species.1
names(MortTrait)<-Nature_Paper_Infection_Mort_Perc$Species.1

#Using the phyEstimate function to predict infection and mortality for species not challenged with Bsal that are present in US
EstimatedInfectionData<-phyEstimate(CaudataTree,InfectionTrait,method="pic")
EstimatedInfectionData
EstimatedMortData<-phyEstimate(CaudataTree,MortTrait,method="pic")
EstimatedMortData

#Making sure the predicted values are numeric
Infection<-as.numeric(EstimatedInfectionData$estimate)
Mortality<-as.numeric(EstimatedMortData$estimate)

#Combining the infeciton and mortality values for each species into a dataframe
Infection_Mort_CatData<-as.data.frame(cbind(rownames(EstimatedInfectionData),Infection,Mortality))
str(Infection_Mort_CatData)
Infection_Mort_CatData$Species.1<-Infection_Mort_CatData$V1

#Reformatting original dataframe and new data so they can be combined
Infection_Mort_CatData<-Infection_Mort_CatData[,c(1:3)]
Infection_Mort_CatData$Infection<-as.numeric(Infection_Mort_CatData$Infection)
Infection_Mort_CatData$Mortality<-as.numeric(Infection_Mort_CatData$Mortality)
Nature_Paper_Infection_Mort_Perc_Reduced<-Nature_Paper_Infection_Mort_Perc[,c(1,4:5)]
#Renaming column from new dataframe
names(Infection_Mort_CatData)[names(Infection_Mort_CatData) == "V1"] <- "Species.1"


#Making sure the original infection and mortality values are numeric 
Nature_Paper_Infection_Mort_Perc_Reduced$Infection<-as.numeric(Nature_Paper_Infection_Mort_Perc_Reduced$Infection)
Nature_Paper_Infection_Mort_Perc_Reduced$Mortality<-as.numeric(Nature_Paper_Infection_Mort_Perc_Reduced$Mortality)


#Combining original and new dataframe into one dataframe for export
CombinedData<-dplyr::bind_rows(Infection_Mort_CatData,Nature_Paper_Infection_Mort_Perc_Reduced)
CombinedData
rownames(CombinedData) <- c()

#Creating a new xls file 
write.xlsx(CombinedData, "CombinedData.xls")


#Creating scaled data and a new dataframe for GIS Analysis. Data imported are those calculatedfor each US county (FIPS) where salamanders are present.
BsalRiskData <- read.csv("BsalRiskData.csv")
str(BsalRiskData)

#Scaling Variables from 1-4
BsalRiskData<-within(BsalRiskData,MortScale<-((BsalRiskData$AverageMort-min(BsalRiskData$AverageMort))*3/(max(BsalRiskData$AverageMort)-min(BsalRiskData$AverageMort))+1))

BsalRiskData<-within(BsalRiskData,InfScale<-((BsalRiskData$InfectionPerc-min(BsalRiskData$InfectionPerc))*3/(max(BsalRiskData$InfectionPerc)-min(BsalRiskData$InfectionPerc))+1))

BsalRiskData<-within(BsalRiskData,RichnessScale<-((BsalRiskData$Richness-min(BsalRiskData$Richness))*3/(max(BsalRiskData$Richness)-min(BsalRiskData$Richness))+1))

BsalRiskData<-within(BsalRiskData,ResistantSpeciesCountScale<-((BsalRiskData$ResistantSpeciesCount-min(BsalRiskData$ResistantSpeciesCount))*3/(max(BsalRiskData$ResistantSpeciesCount)-min(BsalRiskData$ResistantSpeciesCount))+1))

BsalRiskData<-within(BsalRiskData,MeanTempScale<-ifelse(BsalRiskData$AverageAnnualTemp<5.9,
                                                        ((BsalRiskData$AverageAnnualTemp-min(BsalRiskData$AverageAnnualTemp))*3/(5.9-min(BsalRiskData$AverageAnnualTemp))+1),
                                                        ifelse(BsalRiskData$AverageAnnualTemp>5.9&BsalRiskData$AverageAnnualTemp<=14,4,
                                                               ifelse(BsalRiskData$AverageAnnualTemp>14&BsalRiskData$AverageAnnualTemp<=22,
                                                                      ((((14-BsalRiskData$AverageAnnualTemp)/(22-14))*(4-1))+4),
                                                                      ifelse(BsalRiskData$AverageAnnualTemp>22,1,NA)))))

BsalRiskData<-within(BsalRiskData,MeanTempWarmestMonthScale<-ifelse(BsalRiskData$MeanTempWarmestMonth<6,
                                                                    ((BsalRiskData$MeanTempWarmestMonth-min(BsalRiskData$MeanTempWarmestMonth))*3/(6-min(BsalRiskData$MeanTempWarmestMonth))+1),
                                                                    ifelse(BsalRiskData$MeanTempWarmestMonth<=15,4,
                                                                           ifelse(BsalRiskData$MeanTempWarmestMonth>15&BsalRiskData$MeanTempWarmestMonth<=22,
                                                                                  ((((15-BsalRiskData$MeanTempWarmestMonth)/(22-15))*(4-1))+4),
                                                                                  ifelse(BsalRiskData$MeanTempWarmestMonth>22,1,NA)))))

#Creating Combined Scales
BsalRiskData$Temperature_Suitability_score<-(BsalRiskData$MeanTempScale+BsalRiskData$MeanTempWarmestMonthScale)/2
BsalRiskData$Invasion_Risk_score<-(BsalRiskData$Temperature_Suitability_score+BsalRiskData$InfScale)/2
BsalRiskData$Extinction_Risk_score<-(BsalRiskData$Temperature_Suitability_score+BsalRiskData$InfScale+BsalRiskData$MortScale)/3
BsalRiskData$Biodiversity_Risk_score<-(BsalRiskData$Temperature_Suitability_score+BsalRiskData$InfScale+BsalRiskData$MortScale+BsalRiskData$RichnessScale)/4


write.csv(BsalRiskData,"SalamanderRiskMod_GIS_Dataframe.csv", row.names = FALSE)