Commit 084419c9 authored by de Freitas Costa, Eduardo's avatar de Freitas Costa, Eduardo
Browse files

public commit

parent 7a08f60b
......@@ -114,20 +114,57 @@ dose_veg<-numeric()
dose_veg1<-numeric()
OC<-list()
FA<-list()
Ce<-list()
ME<-list()
media_H2<-numeric()
media_C2<-numeric()
media_F2<-numeric()
R1<-list()
R2<-list()
R3<-list()
R4<-list()
R6<-list()
R6.1<-list()
R7<-list()
ERE9<-list()
ERE10<-list()
A<-list()
Prev_slaug<-list()
Raw_af<-list()
Raw<-list()
Veg_cross<-list()
Raw_cross<-list()
Ready<-list()
Ready2<-list()
Dose_chi<-list()
Dose_veg<-list()
## Infection dynamics in Open commmunity
prev1<-0.05 #prev in human time zero
col_proc<- c(1/27,1/27,1/27,1/27,1/27,1/27,1/27,1/27,
1/57,1/12,1/27,1/27,1/27,1/27,1/27,1/27,
1/27,1/27,1/27,1/27,1/27,1/27,1/27,1/27,1/27) #rate of within-household colonization (week) Haverkate et al. (2017)
1/27,1/27,1/27,1/27,1/27,1/27,1/27,1/27,1/27,1/27,1/27,1/27,1/27) #rate of within-household colonization (week) Haverkate et al. (2017)
fade_h<-c(1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,
1/16,1/31,1/8,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16) # Rate of colonization clearance in healthy people Haverkate et al. (2017)
1/16,1/31,1/8,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16) # Rate of colonization clearance in healthy people Haverkate et al. (2017)
#beta(154,10625)=(0.012; 0.014; 0.016)/4
expo_fa<-c(0.003645,0.003092674,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645
,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.0042,0.003645,0.003645,0.003645,0.003645)*(500/1800) #Proportion of humans in open community which have contact with BF farmers Mughini-Gras et al. (2019) Supp appendix S2.
,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.0042,0.003645,0.003645,0.003645,0.003645)*(1800/1800) #Proportion of humans in open community which have contact with BF farmers Mughini-Gras et al. (2019) Supp appendix S2.
#beta(154,10625)=(0.012; 0.014; 0.016)/4
expo_fa1<-c(0.003645,0.003092674,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645,0.003645
......@@ -144,7 +181,7 @@ rate1<-c(0.19,0.19,0.19,0.19,0.19,0.19,0.19,0.19,0.19,0.19,0.19,0.19,0.19,0.19
col_proc2<- 1-exp(-1/27) #Colonization given the contact with open community Haverkate et al. (2017)
fade_fa<-c(1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,
1/16,1/31,1/8,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16)#average time (week) to colonization cleareance in healthy people Haverkate et al. (2017)
1/16,1/31,1/8,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16,1/16)#average time (week) to colonization cleareance in healthy people Haverkate et al. (2017)
......@@ -190,8 +227,8 @@ rate2.1<-c(0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.01,
#rate4<-c(0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.01,0.2,0.1,0.1,0.1,0.1) #Prob of colonization of a calve given the contact with OC: Guess
heavy_prop<-c(0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.1,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,
0.01,0.001,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01) #proportion of heavy ATM flocks No data, scenario created based on Luiken et al. (2019)
heavy_prop<-c(0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.015,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,
0.01,0.0075,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01) #proportion of heavy ATM flocks No data, scenario created based on Luiken et al. (2019)
heavy_b<-c(1.46,1.46,1.46,1.46,1.46,1.46,1.46,1.46,1.46,1.46,1.46,1.46,1.46,1.46,1.46,1.46,1.46,1.46,1.46,1.46,1.46,1.05,2,1.46,1.46,1.46,1.46,1.46,1.46) #effect of high ATU on the beta_ch Luiken et al. (2019) Tb2 and Fig2
......@@ -223,7 +260,7 @@ sig_conc1<-0.67 # SD concentration o ESBL in feces log10(CFU/gram) Ceccarel
# Scalding Defeathering
# Defeathering
aext2<-0.65 # tranfer from carcass to environ during scalding
cext2<-0 # inactivation/removal from carcass surface during scalding https://browser.combase.cc/ComBase_Predictor.aspx?model=2 & Pacholewicz et al. 2016
benv2<-c(0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02) # transfer from envir to carcass surface during scalding
......@@ -236,7 +273,7 @@ conc2<-5.68 # concentration o ESBL in feces log10(CFU/gram) Ceccarelli
sig_conc2<-0.67 # SD concentration o ESBL in feces log10(CFU/gram) Ceccarelli et al. (2017)
# Scalding evisceration
# Evisceration
aext3<-0.053 # tranfer from carcass to environ during scalding
cext3<-0.5 # inactivation/removal from carcass surface during scalding https://browser.combase.cc/ComBase_Predictor.aspx?model=2 & Pacholewicz et al. 2016
benv3<-c(0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02) # transfer from envir to carcass surface during scalding
......@@ -249,7 +286,7 @@ conc3<-5.68 # concentration o ESBL in feces log10(CFU/gram) Ceccarelli
sig_conc3<-0.67 # SD concentration o ESBL in feces log10(CFU/gram) Ceccarelli et al. (2017)
# Scalding chilling
# Chilling
aext5<-0.01 # transfer from carcass to environ during scalding
cext5<-0.74 # inactivation/removal from carcass surface during scalding https://browser.combase.cc/ComBase_Predictor.aspx?model=2 & Pacholewicz et al. 2016
benv5<-c(0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02) # transfer from envir to carcass surface during scalding
......
......@@ -36,16 +36,56 @@
#Script start#############################################################################
####################
# Baseline outputs #
####################
#Final prev/week
# Final risk/week
risk <- read.csv(here("Output","chicken_risk","scenario1","data1.txt"), sep="")
prev <- read.csv(here("Output","chicken","scenario1.txt"),sep=" ")
1-((1-risk$r1[dim(risk)[1]])*(1-risk$r2[dim(risk)[1]])*(1-risk$r3[dim(risk)[1]])*(1-risk$r4[dim(risk)[1]]))
colo1<-cbind.data.frame(prev,time=seq(1,200,1))
1-((1-risk$r6[dim(risk)[1]])*(1-risk$r6.1[dim(risk)[1]])*(1-risk$r7[dim(risk)[1]])*(1-risk$r3[dim(risk)[1]])*(1-risk$r4[dim(risk)[1]]))
colo2<-colo1[,c(1,2,4,5)]%>%gather(key="pop",value="prop",-time)
1-((1-risk$ere9[dim(risk)[1]])*(1-risk$ere10[dim(risk)[1]]))
# Plot with prevalence per week
ggplot(colo2, aes(x=time, y=prop, linetype=pop)) +
geom_line(size=1.01)+
ylab("Prevalence")+
xlab("Time (weeks)")+
theme_minimal()+
theme(legend.position="bottom",legend.title=element_blank(),
axis.text.x = element_text(colour = 'black',
size = 20),
axis.text.y = element_text(colour = 'black',
size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))+
theme(legend.text=element_text(size=15), legend.spacing.x = unit(0.5, 'cm'),
legend.key.height= unit(1, 'cm'),
legend.key.width= unit(1.3, 'cm'))+
scale_linetype_manual(values=c(1, 2, 3),breaks =c("Fa", "H", "media") ,
labels=c("Farmers","Open community", "Flocks"),element_text(size=10))+
scale_y_continuous(labels = function(x) paste0(x*100, "%"),breaks = seq(0,0.35,0.05))
ggsave(here("Figures","chicken",'Fig_2.png'), dpi = 300, height = 5, width = 7, unit = 'in',bg="white")
######################
# Source attribution #
######################
risk <- read.csv(here("Output","chicken_risk","scenario1.txt"),sep=" ")
1-((1-tail(risk$r1,1))*(1-tail(risk$r2,1))*(1-tail(risk$r3,1))*(1-tail(risk$r4,1)))
1-((1-tail(risk$r6,1))*(1-tail(risk$r6.1,1))*(1-tail(risk$r7,1))*(1-tail(risk$r3,1))*(1-tail(risk$r4,1)))
1-((1-tail(risk$ere9,1))*(1-tail(risk$ere10,1)))
# Table with source attribution
......@@ -88,15 +128,17 @@ pcol_fl<-tail(h7/(h5+h6+h7+h8+h9),1)
pcol_me<-tail(h8/(h5+h6+h7+h8+h9),1)
pcol_ve<-tail(h9/(h5+h6+h7+h8+h9),1)
pcol_oc;pcol_fa;pcol_fl;pcol_me;pcol_ve
#For flocks
h10<-(-log(1-risk$ere9)) #flock
h11<-(-log(1-risk$ere10)) #farmer
pcol_fl<-tail(h10/(h10+h11),1)
pcol_fa<-tail(h11/(h10+h11),1)
pcol_fl;pcol_fa
# Graph for source attribution
attr<-cbind.data.frame(Time=1:200,OC=risk$r1,Farmer=risk$r2,Meat=risk$r3,Veg=risk$r4)
......@@ -177,111 +219,170 @@ at3<-ggplot(attr3, aes(x=Time, y=Attribution, linetype=Source)) +
grid.arrange(at1, at2,at3, ncol=2,nrow = 2)
#################
# Contamination #
#################
conta <-read.csv(here("Output","chicken_cont","scenario1.txt"), encoding="UTF8", sep=" ")
#Final colonization week 200
colo <- read.csv(here("Output","chicken","scenario1","data1.txt"), sep="")
cont<-cbind.data.frame(conta,x=1:200)
base<-colo[200,1]
colo1.2<- gather(colo,key="pop",value="prop")
colo1.2$time<-rep(seq(1,200,1),4)
colo2<-subset(colo1.2,pop %in% c("H","media","Fa"))
colo2$pop<-as.factor(colo2$pop)
ggplot(colo2, aes(x=time, y=prop, linetype=pop)) +
geom_line(size=1.01)+
ylab("Prevalence")+
xlab("Time (weeks)")+
theme_minimal()+
theme(legend.position="bottom",legend.title=element_blank(),
axis.text.x = element_text(colour = 'black',
size = 20),
axis.text.y = element_text(colour = 'black',
size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))+
theme(legend.text=element_text(size=15), legend.spacing.x = unit(0.5, 'cm'),
legend.key.height= unit(1, 'cm'),
legend.key.width= unit(1.3, 'cm'))+
scale_y_continuous(labels = function(x) paste0(x*100, "%"),breaks = seq(0,0.35,0.05))+
scale_linetype_manual(values=c(1, 2, 3),breaks =c("Fa", "H", "media") ,
labels=c("Farmers", "Open community", "Flocks"),element_text(size=10))
# Contaminations
conta <-read.csv(here("Output","chicken_cont","scenario1","data1.txt"), encoding="UTF8", sep="")
a <- ggplot(data = conta, aes(x = 1:200, y = a)) +
a1 <- ggplot(data = cont, aes(x = x, y = a)) +
theme_minimal()+
geom_line() +
labs(y = "Concentration (log10 cfu/g)", x = "Time (weeks)")+
ggtitle("Carcasses surface")
ggtitle("Carcasses surface")+
scale_y_continuous(breaks = seq(0,1.2,0.2))
b <- ggplot(data = conta, aes(x = 1:200, y = raw)) +
b <- ggplot(data = cont, aes(x = x, y = raw)) +
theme_minimal()+
geom_line() +
labs(y = "Dose (cfu/portion)", x = "Time (weeks)")+
ggtitle("Meat portions before cross-contamination")
c <- ggplot(data = conta, aes(x = 1:200, y = raw_cross)) +
c <- ggplot(data = cont, aes(x = 1:200, y = raw_cross)) +
theme_minimal()+
geom_line() +
labs(y = "Dose (cfu/portion)", x = "Time (weeks)")+
ggtitle("Meat portions after cross-contamination")
d <- ggplot(data = conta, aes(x = 1:200, y = dose_veg)) +
d <- ggplot(data = cont, aes(x = 1:200, y = veg_cross)) +
theme_minimal()+
geom_line() +
labs(y = "Dose (cfu/portion)", x = "Time (weeks)")+
ggtitle("Vegetables portions after cross-contamination")
ggtitle("Vegetables portions after cross-contamination")+
scale_y_continuous(breaks = seq(0,0.35,0.05))
#tiff(file=paste(wd,"/Figure_simulation","/bias.tiff",sep=""), height = 4, width = 6, units = 'in', res=300)
gridExtra::grid.arrange(a,b,c,d,nrow=2)
#dev.off()
png(file=here("Figures","chicken","Fig_S1.png"), height = 4, width = 9, units = 'in', res=300)
gridExtra::grid.arrange(a1,b,c,d,nrow=2)
#ggsave(here("Figures","chicken","Fig_S1.png"), dpi = 300, height = 7, width = 7, unit = 'in')
dev.off()
############
# Tornados #
############
# Data gathering for tornados
cena <- readRDS(here("Output","chicken","cena.rds"))
simtable<-29
nome1<-c("Baseline","Low farm exposure", "Low spread flock", "Low chicken meat consumption", "Low cross-cont. kitchen","High cross-cont. kitchen", "High cont. SH*","High prop. heavy ATM** user",
nome1<-c("Baseline","Low farm exposure", "Low spread flock", "Low chicken meat consumption", "Low cross-cont. kitchen","High cross-cont. kitchen", "Low cont. SH*","High prop. heavy ATM** user",
"Low colh" ,"High colh","Low fadeh","High fadeh","Low colc", "High colc", "Low rate1", "High rate1",
"Low rate2", "High rate2", "Low prop. heavy ATM** user", "Low fadec","High facec","Low effect ATM**", "High effect ATM**","Low meat portion size","High farm exposure","Low delta", "High delta", "Low gamma","High gamma")
nome2<-nome1[c(3:8,19,24)]
nome3<-nome1[c(2,9:18,20:23,simtable)]
dado1<-cbind.data.frame(prev=cena[c(3:8,19,24)],name=nome2)
dado1<-cbind.data.frame(cena[c(3:8,19,24)],nome2)
dado2<-cbind.data.frame(cena[c(2,9:18,20:23,simtable)],nome3)
nome3<-nome1[c(2,9:18,20:23,25:simtable)]
dado2<-cbind.data.frame(prev=cena[c(2,9:18,20:23,25:simtable)],name=nome3)
base<-cena[1]
# Scenarios for humans
#percentage difference
teste<-cbind.data.frame(perdiff=(dado1$`cena[c(3:8, 19, 24)]`-base)/base,name=dado1$nome2)
ggplot(teste,aes(y=perdiff,x=reorder(name, abs(perdiff))))+
geom_bar(stat = "identity")+
ylab("Percentage difference")+
# Uncertainty for humans
#Percentage difference
teste2<-cbind.data.frame(perdiff=(dado2$prev-base)/base,name=dado2$name)
ggplot(teste2,aes(y=perdiff,x=reorder(name, abs(perdiff))))+
geom_bar(stat = "identity")+
ylab("Log difference")+
theme_minimal()+
theme(text = element_text(size=15))+
xlab(" ")+
scale_y_continuous(labels = scales::percent)+
xlab(" ")+
coord_flip()
#crude prevalence
teste3<-cbind.data.frame(prev=dado2$prev,name=dado2$name)
teste3<-rbind(teste3,teste3[14:15,])
teste3$name[14]<-"Low_heavy_effect_c"
teste3$name[15]<-"High_heavy_effect_c"
teste3$name[21]<-"Low_heavy_effect_f"
teste3$name[22]<-"High_heavy_effect_f"
teste3$sce<-c("Lower","Lower","Upper","Lower","Upper","Lower","Upper","Lower","Upper","Lower",
"Upper","Lower","Upper","Lower","Upper","Upper","Lower","Upper","Lower","Upper","Lower","Upper")
teste3<-spread(teste3,key = "sce","prev")
teste3<-cbind(teste3[12:22,c(1,2)],teste3[1:11,c(1,3)])
teste3$abs<-abs(teste3$Upper-teste3$Lower)
teste3<-teste3[with(teste3, order(-abs)), ]
teste3<-teste3[,-c(3)]
parameter<-c("beta^h","theta^h","lambda^beta^c","lambda^theta^c","beta^c","p^fc","theta^c","delta","gamma","p^cf","p^hf")
teste3$Parameter<-parameter
base.value<-base
order.parameters <- teste3 %>% arrange(abs) %>%
mutate(Parameter=factor(x=Parameter, levels=Parameter)) %>%
select(Parameter) %>% unlist() %>% levels()
width<-0.7
# get data frame in shape for ggplot and geom_rect
df.3 <- teste3 %>%
# gather columns Lower_Bound and Upper_Bound into a single column using gather
gather(key='type', value='output.value', Lower:Upper) %>%
# just reordering columns
select(Parameter, type, output.value, abs) %>%
# create the columns for geom_rect
mutate(Parameter=factor(Parameter, levels=order.parameters),
ymin=pmin(output.value, base.value),
ymax=pmax(output.value, base.value),
xmin=as.numeric(Parameter)-width/1.5,
xmax=as.numeric(Parameter)+width/1.5)
teste1<-cbind.data.frame(prev=(dado1$`cena[c(3:8, 19, 24)]`),name=dado1$nome2)
# create plot
# (use scale_x_continuous to change labels in y axis to name of parameters)
.expressions <- paste(order.parameters)
legend_expressions <-parse(text = .expressions)
ggplot() +
geom_rect(data = subset(df.3,abs>0.0024572273),
aes(ymax=ymax, ymin=ymin, xmax=xmax, xmin=xmin), fill=c(rep("gray",5),rep("black",5))) +
theme_minimal()+
theme(axis.title.y=element_blank(), legend.position = 'right',
legend.title = element_blank()) +
geom_hline(yintercept = base.value) +
scale_x_continuous(breaks = c(1:length(order.parameters)),
labels = legend_expressions) +
scale_y_continuous(breaks=seq(0,0.3,0.05),labels = scales::percent_format(accuracy=1))+
theme(text = element_text(size=15))+
ylab("Prevalence at week 200")+
coord_flip()
ggsave(here("Figures","chicken",'Fig_3.png'), dpi = 300, height = 5, width = 7, unit = 'in',bg="white")
## Scenarios for humans
#crude prevalence all scenarios
teste1<-cbind.data.frame(prev=dado1$prev,name=dado1$name)
teste1$Lower<-ifelse(teste1$prev<base,teste1$prev,base)
teste1$Upper<-ifelse(teste1$prev>base,teste1$prev,base)
teste1$UL_Difference<-abs(teste1$Upper-teste1$Lower)
......@@ -324,9 +425,14 @@ gridExtra::grid.arrange(a,b,c,d,nrow=2)
#Only reductions
teste1<-cbind.data.frame(prev=(dado1$`cena[c(3:8, 19, 24)]`),name=dado1$nome2)
teste2<-subset(teste1,name!="High cross-cont. kitchen" & name!="High prop. heavy AMU** user")
teste2<-subset(teste1,name!="High cross-cont. kitchen" & name!="High prop. heavy ATM** user")
teste2$name<-c("Reducing spread flock",
"Reducing meat consumption",
"Reducing cross-contamination kitchen",
"Reducing contamination slaughterhouse",
"Reducing proportion heavy AMU*",
"Reducing meat portion size")
teste2$Lower<-ifelse(teste2$prev<base,teste2$prev,base)
teste2$Upper<-ifelse(teste2$prev>base,teste2$prev,base)
......@@ -368,155 +474,7 @@ gridExtra::grid.arrange(a,b,c,d,nrow=2)
scale_y_continuous(breaks=seq(0.0032,base,0.001),labels = scales::percent)+
coord_flip()
ggsave(here("Figures","chicken",'Fig_4.png'), dpi = 300, height = 5, width = 7, unit = 'in',bg="white")
# Uncertianty for humans
#Percentage difference
teste2<-cbind.data.frame(perdiff=(dado2$`cena[c(2, 9:18, 20:23, simtable)]`-base)/base,name=dado2$nome3)
ggplot(teste2,aes(y=perdiff,x=reorder(name, abs(perdiff))))+
geom_bar(stat = "identity")+
ylab("Log difference")+
theme_minimal()+
theme(text = element_text(size=15))+
scale_y_continuous(labels = scales::percent)+
xlab(" ")+
coord_flip()
Parameter<-c("colh","fadeh","colc","Heavy_effect_c","Heavy_effect_f","rate2")
Lower_Bound<-c(-0.55079302,13.31194989 ,-0.54246917,-0.68866605,-0.68866605,-0.16532534)
Upper_Bound<-c(31.79585354,-0.71343386,1.24034604,0.91456462,0.91456462,0.15817645)
UL_Difference<-abs(Upper_Bound-Lower_Bound)
df<-cbind.data.frame(Parameter,Lower_Bound,Upper_Bound,UL_Difference)
base.value <- 0
order.parameters <- df %>% arrange(UL_Difference) %>%
mutate(Parameter=factor(x=Parameter, levels=Parameter)) %>%
select(Parameter) %>% unlist() %>% levels()
# width of columns in plot (value between 0 and 1)
width <- 0.7
# get data frame in shape for ggplot and geom_rect
df.2 <- df %>%
# gather columns Lower_Bound and Upper_Bound into a single column using gather
gather(key='type', value='output.value', Lower_Bound:Upper_Bound) %>%
# just reordering columns
select(Parameter, type, output.value, UL_Difference) %>%
# create the columns for geom_rect
mutate(Parameter=factor(Parameter, levels=order.parameters),
ymin=pmin(output.value, base.value),
ymax=pmax(output.value, base.value),
xmin=as.numeric(Parameter)-width/1.5,
xmax=as.numeric(Parameter)+width/1.5)
# create plot
# (use scale_x_continuous to change labels in y axis to name of parameters)
ggplot() +
geom_rect(data = df.2,
aes(ymax=ymax, ymin=ymin, xmax=xmax, xmin=xmin), fill=c(rep("gray",6),rep("black",6))) +
theme_minimal()+
theme(axis.title.y=element_blank(), legend.position = 'NULL',
legend.title = element_blank()) +
geom_hline(yintercept = base.value) +
scale_x_continuous(breaks = c(1:length(order.parameters)),
labels = order.parameters) +
scale_y_continuous(labels = scales::percent)+
theme(text = element_text(size=15))+
ylab("Percentage difference")+
coord_flip()
#crude prevalence
teste3<-cbind.data.frame(prev=(dado2$`cena[c(2, 9:18, 20:23, simtable)]`),name=dado2$nome3)
teste3<-rbind(teste3,teste3[14:15,])
teste3$name[14]<-"Low_heavy_effect_c"
teste3$name[15]<-"High_heavy_effect_c"
teste3$name[17]<-"Low_heavy_effect_f"
teste3$name[18]<-"High_heavy_effect_f"
teste3$sce<-c("Lower","Lower","Upper","Lower","Upper","Lower","Upper","Lower","Upper","Lower","Upper","Lower","Upper","Lower","Upper","Upper","Lower","Upper")
teste3<-spread(teste3,key = "sce","prev")
teste3<-cbind(teste3[10:18,c(1,2)],teste3[1:9,c(1,3)])
teste3$abs<-abs(teste3$Upper-teste3$Lower)
teste3<-teste3[with(teste3, order(-abs)), ]
teste3<-teste3[,-c(3)]
parameter<-c("beta^h","theta^h","lambda^beta^c","lambda^theta^c","beta^c","biosec^fc","biosec^cf","p^hf","theta^c")
teste3$Parameter<-parameter
base.value<-base
order.parameters <- teste3 %>% arrange(abs) %>%
mutate(Parameter=factor(x=Parameter, levels=Parameter)) %>%
select(Parameter) %>% unlist() %>% levels()
width<-0.7
# get data frame in shape for ggplot and geom_rect
df.3 <- teste3 %>%
# gather columns Lower_Bound and Upper_Bound into a single column using gather
gather(key='type', value='output.value', Lower:Upper) %>%
# just reordering columns
select(Parameter, type, output.value, abs) %>%
# create the columns for geom_rect
mutate(Parameter=factor(Parameter, levels=order.parameters),
ymin=pmin(output.value, base.value),
ymax=pmax(output.value, base.value),
xmin=as.numeric(Parameter)-width/1.5,
xmax=as.numeric(Parameter)+width/1.5)
# create plot
# (use scale_x_continuous to change labels in y axis to name of parameters)
.expressions <- paste(order.parameters)