Commit 7c935dd1 authored by de Freitas Costa, Eduardo's avatar de Freitas Costa, Eduardo
Browse files

Changes in the Heavy AMU

parent cb330326
......@@ -29,7 +29,7 @@ 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","High ATM use",
nome1<-c("Baseline","No Farm exposure", "Low spread flock", "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")
......@@ -66,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])**(1-exp(-1/col_proc[j]))
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]))
......@@ -74,7 +74,7 @@ set.seed(13)
p4[1]<- r8[1]
r9[[1]]<- (1-exp(-beta_ch[j]*time_flock[[1]]))*prev3
r9[[1]]<- (1-exp(-(beta_ch[j]*prev3*time_flock[[1]])))
r10[[1]]<- (Fa[1])*(1-prev3)*rate2[j]
r11[[1]]<- (1-exp(-(1/fade_ch[j])))
p5[[1]]<- (r9[[1]]+r10[[1]])-(r9[[1]]*r10[[1]])
......@@ -245,8 +245,18 @@ set.seed(13)
amu[[t]]<-amu[[t-1]]
time_flock[[t]]<-ifelse (time_flock[[t-1]]<6,time_flock[[t-1]]+1,0)
r9[[t]]<-ifelse(time_flock[[t]]==1,(1-exp(-beta_ch[j]*prev3)),(1-exp(-beta_ch[j]*flock[[t-1]])) )
ifelse(amu[[t]]!=1,
r9[[t]]<-ifelse(time_flock[[t]]==1,(1-exp(-beta_ch[j]*prev3)),(1-exp(-beta_ch[j]*flock[[t-1]]))),
r9[[t]]<-ifelse(time_flock[[t]]==1,(1-exp(-beta_ch[j]*heavy[j]*prev3)),(1-exp(-beta_ch[j]*heavy[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[j])))))
p5[[t]]<-1-((1-r9[[t]])*(1-r10[[t]]))
......@@ -427,30 +437,30 @@ tornado<-function(){
par(mfrow=c(1,1))
base<-cena[1]
nome2<-nome1[2:9]
nome3<-nome1[10:17]
nome4<-nome1[18:simtable]
nome2<-nome1[2:8]
nome3<-nome1[9:16]
nome4<-nome1[17:simtable]
dado1<-cbind.data.frame(cena[2:9],nome2)
dado1<-cbind.data.frame(cena[2:8],nome2)
dado2<-cbind.data.frame(cena[10:17],nome3)
dado2<-cbind.data.frame(cena[9:16],nome3)
dado3<-cbind.data.frame(cena[18:simtable],nome4)
dado3<-cbind.data.frame(cena[17:simtable],nome4)
sorted1<-dado1[order(abs(log(cena[2:9]/base))),]
sorted1<-dado1[order(abs(log(cena[2:8]/base))),]
sorted2<-dado1[order(abs(((cena[2:9]-base)))),]
sorted2<-dado1[order(abs(((cena[2:8]-base)))),]
sorted3<-dado2[order(abs(log(cena[10:17]/base))),]
sorted3<-dado2[order(abs(log(cena[9:16]/base))),]
sorted4<-dado2[order(abs(((cena[10:17]-base)))),]
sorted4<-dado2[order(abs(((cena[9:16]-base)))),]
sorted5<-dado3[order(abs(log(cena[18:simtable]/base))),]
sorted5<-dado3[order(abs(log(cena[17:simtable]/base))),]
sorted6<-dado3[order(abs(((cena[18:simtable]-base)))),]
sorted6<-dado3[order(abs(((cena[17:simtable]-base)))),]
saveRDS(sorted1,paste(wd,"/output/sorted1.rds",sep=""))
......@@ -568,24 +578,26 @@ for (i in 1:simtable){
}
nome4<-nome1[c(1,2,3,4,5,6,8,9,10,11,12)]
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")
plot(1:200,human_prev[[1]]$H,type="l",col=1,lwd=2,ylab="ESBL prevalence in the open community",xlab="Time in weeks",ylim = c(0,0.15))
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)
lines(1:200,human_prev[[10]]$H,type="l",col=10)
lines(1:200,human_prev[[11]]$H,type="l",col=11)
lines(1:200,human_prev[[12]]$H,type="l",col=10)
nome2<-nome1[1:9]
legend(130, 0.05, legend=nome2,
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)
legend("topleft", legend=nome4,
col=c(1,2,3,4,5,6,8,9,10,11,12), lty=1,lwd=c(2,1,1,1,1,1,1,1,1,1,1,1), cex=0.8, text.font=4)
dev.off()
......
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