library("fast")
library("sensitivity")
library("ggplot2")
library("RColorBrewer")

#Generate parameter values
para_model=fast_parameters(minimum = c(0,0,0,0,0,0,0,0), maximum = c(1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0), factor=9,names = c("betaha", "betaah", "betahh", "betaaa","lambdah", "lambdaa", "muh","mua"))

#Output to csv-file
write.csv(para_model,"m:/animal-human_model/paras_sens_analysis_manuscript.csv")

#Read in the output from equilibrium analysis
outputequilibrium=unlist(read.csv("m:/animal-human_model/outcomeequilibrium.csv", header = F), use.names = F)

#Sensitivity analysis
sens<-sensitivity(x=outputequilibrium, numberf=8, make.plot=T, names = c("betaha", "betaah","betahh", "betaaa", "lambdah", "lambdaa", "muh","mua"))

#Plot partial variances (Figure 2)
df.equilibrium <- data.frame(parameter=rbind("betaha", "betaah","betahh", "betaaa", "lambdah", "lambdaa", "muh","mua"), value=sens)
windows(12,8)
p <- ggplot(df.equilibrium, aes(parameter, value))
p + geom_bar(stat="identity", fill=brewer.pal(3,"Set1")[2]) + 
   scale_x_discrete("Parameter", waiver(), c(expression(paste(beta[AA])), expression(paste(beta[AH])),expression(paste(beta[HA])), expression(paste(beta[HH])), expression(paste(Lambda[A])), expression(paste(Lambda[H])), expression(paste(mu[A])),expression(paste(mu[H])))) +
  ylab("Partial Variance") + theme_bw()+theme(axis.text=element_text(size = 16, colour = "black"), axis.title=element_text(size=20))


#Read in the output from impact analysis
outputimpact=unlist(read.csv("m:/animal-human_model/outcomeimpact.csv", header = F), use.names = F)

#Sensitivity analysis
sens<-sensitivity(x=outputimpact, numberf=8, make.plot=T, names = c("betaha", "betaah","betahh", "betaaa", "lambdah", "lambdaa", "muh","mua"))

#Plot partial variances (Figure 4)
df.impact <- data.frame(parameter=rbind("betaha", "betaah","betahh", "betaaa", "lambdah", "lambdaa", "muh","mua"), value=sens)
windows(12,8)
p <- ggplot(df.impact, aes(parameter, value))
p + geom_bar(stat="identity", fill=brewer.pal(3,"Set1")[2]) + 
  scale_x_discrete("Parameter", waiver(), c(expression(paste(beta[AA])), expression(paste(beta[AH])),expression(paste(beta[HA])), expression(paste(beta[HH])), expression(paste(Lambda[A])), expression(paste(Lambda[H])), expression(paste(mu[A])),expression(paste(mu[H])))) +
  ylab("Partial Variance") + theme_bw()+theme(axis.text=element_text(size = 16, colour = "black"), axis.title=element_text(size=20))