Commit 4f4414cd authored by de Freitas Costa, Eduardo's avatar de Freitas Costa, Eduardo
Browse files

Sensitivity analysis

parent 785450e3
......@@ -29,13 +29,9 @@ set.seed(13)
wd.dados <- paste(wd,"/output","/output_scenario",j,sep="")
dir.create(wd.dados,showWarnings = F)
nome1<-c("Baseline","No Farm exposure","High beta_ch", "Low beta_ch", "No food", "No cross-cont. kitchen","Full cross-cont kitchen", "High cont at SH",
"no name","no name" ,"no name","no name","no name","no name", "no name", "no name")
nome<-c("No Farm exposure","High beta_ch", "Low beta_ch", "No food", "No cross-cont. kitchen","Full cross-cont kitchen", "High cont at SH",
"no name","no name" ,"no name","no name","no name","no name", "no name", "no name")
nome1<-c("Baseline","No Farm exposure","High beta_ch", "Low beta_ch", "No food", "No cross-cont. kitchen","Full cross-cont kitchen", "High cont at SH","High ATM use",
"Low col_prob1" ,"High col_prob1","Low fade_h","High fade_h","Low beta_c", "high beta_c", "Low rate1", "High rate1",
"Low rate2", "High rate2", "Low heavy_prop", "High heavy_prop","Low heavy_effect", "High heavy_effect")
par(mar=c(5,5,4,3))
......@@ -61,8 +57,8 @@ set.seed(13)
amu[[1]]<-rbinom(flocks,1,heavy_prop[j])
r1[1]<- (H[1]*(1-H[1]))*col_proc
r2[1]<- (Fa[1])*(1-H[1])*expo_fa[j]*col_proc
r1[1]<- (H[1]*(1-H[1]))*(1-exp(-1/col_proc[j]))
r2[1]<- (Fa[1])*(1-H[1])*expo_fa[j]*(1-exp(-1/col_proc[j]))
r3[1]<-(1 - (1 + ((dose_chi[1] + dose_veg[1])/alfa))^-gama)*non_veg*prob_consu[j]*prev_chi[1]
r4[1]<-(1 - (1 + ((dose_veg[1])/alfa))^-gama)*non_veg*prob_consu[j]*prev_veg[1]
r5[1]<- (1-exp(-1/fade_h[j])) # decolonization probability in one week
......@@ -70,7 +66,7 @@ set.seed(13)
p1[1]<-1-((1-r1[1])*(1-r2[1])*(1-r3[1])*(1-r4[1]))
p2[1]<-r5[1]
r6[1]<-(1-Fa[1])*(H[1])*col_proc2
r6[1]<-(1-Fa[1])*(H[1])**(1-exp(-1/col_proc[j]))
r7[1]<-(1-Fa[1])*(C[1])*rate1[j]
r8[1]<-(1-exp(-1/fade_fa[j]))
......@@ -214,8 +210,8 @@ set.seed(13)
#Frac of open community infec by open community + Frac of infec by farmers +Frac. of open community infec by food.
r1[t]<- H[t-1]*((1-H[t-1]))*col_proc
r2[t]<- (Fa[t-1])*(1-H[t-1])*expo_fa[j]*col_proc
r1[t]<- H[t-1]*((1-H[t-1]))*(1-exp(-1/col_proc[j]))
r2[t]<- (Fa[t-1])*(1-H[t-1])*expo_fa[j]*(1-exp(-1/col_proc[j]))
r3[t]<-(1 - (1 + ((dose_chi[t-1])/alfa))^-gama)*non_veg*prob_consu[j]*prev_chi[t-1]
r4[t]<-(1 - (1 + ((dose_veg[t-1])/alfa))^-gama)*non_veg*prob_consu[j]*prev_veg[t-1]
r5[t]<- (1-exp(-1/fade_h[j])) # decolonization probability in one week
......@@ -231,7 +227,7 @@ set.seed(13)
###############################
#Frac of farmers infec by open comunity+ Frac of farmers infec by chickens + Frac of farmers infec by food
r6[t]<-(H[t-1])*(1-Fa[t-1])*col_proc2
r6[t]<-(H[t-1])*(1-Fa[t-1])*(1-exp(-1/col_proc[j]))
r7[t]<-(1-Fa[t-1])*(C[t-1])*rate1[j]
r8[t]<-(1-exp(-1/fade_fa[j]))
......@@ -251,7 +247,7 @@ set.seed(13)
r9[[t]]<-ifelse(time_flock[[t]]==1,(1-exp(-beta_ch[j]*prev3)),(1-exp(-beta_ch[j]*flock[[t-1]])) )
r10[[t]]<-ifelse(time_flock[[t]]!=1,(Fa[t-1])*(1-flock[[t-1]])*rate2[j], (Fa[t-1])*(1-prev3)*rate2[j])
r11[[t]]<-ifelse(amu[[t]]!=1,(1-exp(-(1/(fade_ch[j])))), (1-exp(-(1/(fade_ch[j]*heavy)))))
r11[[t]]<-ifelse(amu[[t]]!=1,(1-exp(-(1/(fade_ch[j])))), (1-exp(-(1/(fade_ch[j]*heavy[j])))))
p5[[t]]<-1-((1-r9[[t]])*(1-r10[[t]]))
p6[[t]]<-ifelse(time_flock[[t]]!=1,r11[[t]], r11[[t]])
......@@ -415,7 +411,7 @@ set.seed(13)
par(xpd=NA)
legend('center', legend=c("Ope. Commu.", "Chicken","Farmer","Slaug_prev"),
legend(-1.2,-1, legend=c("Ope. Commu.", "Chicken","Farmer","Slaug_prev"),
horiz=TRUE,col=c("red","green", "blue","yellow"), lty=1:1, cex=1.2)
......@@ -428,28 +424,51 @@ cena<-readRDS(paste(wd,"/output/cena.rds",sep=""))
tornado<-function(){
par(mfrow=c(1,1))
par(mfrow=c(1,1))
base<-cena[1]
nome<-c("No Farm exposure","High beta_ch", "Low beta_ch", "No food", "No cross-cont. kitchen","Full cross-cont kitchen", "High cont at SH",
"no name","no name" ,"no name","no name","no name","no name", "no name", "no name")
nome2<-nome1[2:9]
nome3<-nome1[10:17]
nome4<-nome1[18:simtable]
dado<-cbind.data.frame(cena[2:simtable],nome)
dado1<-cbind.data.frame(cena[2:9],nome2)
dado2<-cbind.data.frame(cena[10:17],nome3)
dado3<-cbind.data.frame(cena[18:simtable],nome4)
sorted1<-dado1[order(abs(log(cena[2:9]/base))),]
sorted<-dado[order(abs(log(cena[2:simtable]/base))),]
sorted2<-dado1[order(abs(((cena[2:9]-base)))),]
sorted3<-dado2[order(abs(log(cena[10:17]/base))),]
sorted4<-dado2[order(abs(((cena[10:17]-base)))),]
sorted5<-dado3[order(abs(log(cena[18:simtable]/base))),]
sorted6<-dado3[order(abs(((cena[18:simtable]-base)))),]
sorted1<-dado[order(abs(((cena[2:simtable]-base)))),]
saveRDS(sorted,paste(wd,"/output/sorted.rds",sep=""))
saveRDS(sorted1,paste(wd,"/output/sorted1.rds",sep=""))
saveRDS(sorted2,paste(wd,"/output/sorted2.rds",sep=""))
saveRDS(sorted3,paste(wd,"/output/sorted3.rds",sep=""))
saveRDS(sorted4,paste(wd,"/output/sorted4.rds",sep=""))
saveRDS(sorted5,paste(wd,"/output/sorted5.rds",sep=""))
saveRDS(sorted6,paste(wd,"/output/sorted6.rds",sep=""))
png(file=here("Figures", "tornado1.png"), width = 465, height = 225, units='mm', res = 300)
par(mar=c(5,13,1,2))
barplot(log(sorted$cena/base), horiz = T,
names.arg = sorted$nome,col = "darkgray",las=1,xlab="Log difference",
ann = FALSE,xlim=c(min(log(sorted$cena/base)),4+max(log(sorted$cena/base))))
barplot(log(sorted1$cena/base), horiz = T,
names.arg = sorted1$nome2,col = "darkgray",las=1,xlab="Log difference",
ann = FALSE,xlim=c(min(log(sorted1$cena/base)),max(log(sorted1$cena/base))))
mtext(side = 1, text = "log difference", line = 3)
......@@ -457,13 +476,66 @@ mtext(side = 1, text = "log difference", line = 3)
dev.off()
png(file=here("Figures", "tornado2.png"), width = 465, height = 225, units='mm', res = 300)
par(mar=c(5,13,1,2))
barplot((sorted1$cena-base), horiz = T,
names.arg = sorted1$nome,col = "darkgray",las=1,xaxt="n",
ann = FALSE,xlim=c(min((sorted1$cena-base)),max((sorted1$cena-base))))
axis(side=1, at=seq(-0.002,0.032,0.001), labels = seq(-0.002,0.032,0.001))
barplot(log(sorted3$cena/base), horiz = T,
names.arg = sorted3$nome3,col = "darkgray",las=1,xlab="Log difference",
ann = FALSE,xlim=c(min(log(sorted3$cena/base)),max(log(sorted3$cena/base))))
mtext(side = 1, text = "log difference", line = 3)
#mtext(side = 2, text = "Scenarios", line = 6.5)
dev.off()
png(file=here("Figures", "tornado3.png"), width = 465, height = 225, units='mm', res = 300)
par(mar=c(5,13,1,2))
barplot(log(sorted5$cena/base), horiz = T,
names.arg = sorted5$nome4,col = "darkgray",las=1,xlab="Log difference",
ann = FALSE,xlim=c(min(log(sorted5$cena/base)),max(log(sorted5$cena/base))))
mtext(side = 1, text = "log difference", line = 3)
#mtext(side = 2, text = "Scenarios", line = 6.5)
dev.off()
png(file=here("Figures", "tornado4.png"), width = 465, height = 225, units='mm', res = 300)
par(mar=c(5,13,1,2))
barplot((sorted2$cena-base), horiz = T,
names.arg = sorted2$nome2,col = "darkgray",las=1,xaxt="n",
ann = FALSE,xlim=c(min((sorted2$cena-base)),max((sorted2$cena-base))))
#axis(side=1, at=seq(-0.002,0.032,0.001), labels = seq(-0.002,0.032,0.001))
mtext(side = 1, text = "Percentage difference", line = 3)
dev.off()
png(file=here("Figures", "tornado5.png"), width = 465, height = 225, units='mm', res = 300)
par(mar=c(5,13,1,2))
barplot((sorted4$cena-base), horiz = T,
names.arg = sorted4$nome3,col = "darkgray",las=1,xaxt="n",
ann = FALSE,xlim=c(min((sorted4$cena-base)),max((sorted4$cena-base))))
#axis(side=1, at=seq(-0.002,0.032,0.001), labels = seq(-0.002,0.032,0.001))
mtext(side = 1, text = "Percentage difference", line = 3)
dev.off()
png(file=here("Figures", "tornado6.png"), width = 465, height = 225, units='mm', res = 300)
par(mar=c(5,13,1,2))
barplot((sorted6$cena-base), horiz = T,
names.arg = sorted6$nome4,col = "darkgray",las=1,xaxt="n",
ann = FALSE,xlim=c(min((sorted6$cena-base)),max((sorted6$cena-base))))
#axis(side=1, at=seq(-0.002,0.032,0.001), labels = seq(-0.002,0.032,0.001))
mtext(side = 1, text = "Percentage difference", line = 3)
dev.off()
......@@ -473,7 +545,7 @@ dev.off()
tornado()
}#end function
human_prev<-list()
......@@ -500,18 +572,22 @@ for (i in 1:simtable){
png(file=here("Figures", "human_prev.png"), width = 465, height = 225, units='mm', res = 300)
par(mar=c(4,4,1,1))
plot(1:200,human_prev[[1]]$H,type="l",col=1,lwd=2,ylab="ESBL prevalence in the open community",xlab="Time in weeks")
lines(1:200,human_prev[[2]]$H,type="l",col=15)
lines(1:200,human_prev[[3]]$H,type="l",col=9)
lines(1:200,human_prev[[4]]$H,type="l",col=14)
lines(1:200,human_prev[[5]]$H,type="l",col=12)
lines(1:200,human_prev[[6]]$H,type="l",col=4)
lines(1:200,human_prev[[7]]$H,type="l",col=5)
lines(1:200,human_prev[[8]]$H,type="l",col=2)
lines(1:200,human_prev[[2]]$H,type="l",col=2)
lines(1:200,human_prev[[3]]$H,type="l",col=3)
lines(1:200,human_prev[[4]]$H,type="l",col=4)
lines(1:200,human_prev[[5]]$H,type="l",col=5)
lines(1:200,human_prev[[6]]$H,type="l",col=6)
lines(1:200,human_prev[[7]]$H,type="l",col=7)
lines(1:200,human_prev[[8]]$H,type="l",col=8)
lines(1:200,human_prev[[9]]$H,type="l",col=9)
nome2<-c("No Farm exposure","High beta_ch", "Low beta_ch", "No food", "No cross-cont. kitchen","Full cross-cont kitchen", "High cont at SH")
nome2<-nome1[1:9]
legend(130, 0.05, legend=nome2,
col=c(1,15,9,14,12,4,5,2), lty=1,lwd=c(2,1,1,1,1,1,1,1), cex=0.8, text.font=4)
col=c(1,2,3,4,5,6,7,8,9), lty=1,lwd=c(2,1,1,1,1,1,1,1,1), cex=0.8, text.font=4)
dev.off()
dev.off()
\ No newline at end of file
}#end function
\ No newline at end of file
......@@ -40,5 +40,5 @@ dir.create("Output",showWarnings = F)
dir.create("Figures",showWarnings = F)
source(file.path("Scripts","dynamic.R"))
model(10,200,500,16)
model(10,200,500,23)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment