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

new uncertainty scenario

parent c0034e94
This diff is collapsed.
#Final colonization week 200
colo <- read.csv(here("Output","multi","scenario1","data1.txt"), sep="")
colo$week<-seq(1,200,1)
base<-colo[200,1]
colo1.2<- gather(colo,key="pop",value="prop",-week)
colo2<-subset(colo1.2,pop %in% c("H","media","media1","Fa","Fa1"))
colo2$pop<-as.factor(colo2$pop)
ggplot(colo2, aes(x=week, 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.55,0.05))+
scale_linetype_manual(values=c(1, 2, 3, 4,5),breaks =c("Fa", "Fa1","H", "media","media1") ,
labels=c("Broiler Farmers","Veal Farmers", "Open community", "Flocks","Herds"),element_text(size=10))
#Source attribution
# Final risk/week
risk <- read.csv(here("Output","multi_risk","scenario1","data1.txt"), sep="")
##Dataset with the sources according Thomas's eq
#For OC
h1<-(-log(1-risk$r1)) #OC
h2<-(-log(1-risk$r2)) #chicken_Farmer
h2.1<-(-log(1-risk$r2.1)) #veal_Farmer
h3<-(-log(1-risk$r3)) #Chicken_meat
h3.1<-(-log(1-risk$r3.1)) #Beef_meat
h4<-(-log(1-risk$r4)) #Chicken_veg
h4.1<-(-log(1-risk$r4.1)) #Beef_veg
pcol_oc<-tail((h1/(h1+h2+h2.1+h3+h3.1+h4+h4.1)),1)
pcol_cfa<-tail(h2/(h1+h2+h2.1+h3+h3.1+h4+h4.1),1)
pcol_vfa<-tail(h2.1/(h1+h2+h2.1+h3+h3.1+h4+h4.1),1)
pcol_cmeat<-tail(h3/(h1+h2+h2.1+h3+h3.1+h4+h4.1),1)
pcol_vmeat<-tail(h3.1/(h1+h2+h2.1+h3+h3.1+h4+h4.1),1)
pcol_chve<-tail(h4/(h1+h2+h2.1+h3+h3.1+h4+h4.1),1)
pcol_beefve<-tail(h4.1/(h1+h2+h2.1+h3+h3.1+h4+h4.1),1)
pcol_oc;pcol_cfa;pcol_vfa;pcol_cmeat;pcol_vmeat;pcol_chve;pcol_beefve
#For chicken farmers
h5<-(-log(1-risk$r6)) #OC
h6<-(-log(1-risk$r6.1)) #Farmer
h7<-(-log(1-risk$r7)) #Flock
pcol_oc1<-tail(h5/(h5+h6+h7+h3+h3.1+h4+h4.1),1)
pcol_fa<-tail(h6/(h5+h6+h7+h3+h3.1+h4+h4.1),1)
pcol_fl<-tail(h7/(h5+h6+h7+h3+h3.1+h4+h4.1),1)
pcol_cmeat1<-tail(h3/(h5+h6+h7+h3+h3.1+h4+h4.1),1)
pcol_vmeat1<-tail(h3.1/(h5+h6+h7+h3+h3.1+h4+h4.1),1)
pcol_chve1<-tail(h4/(h5+h6+h7+h3+h3.1+h4+h4.1),1)
pcol_beefve1<-tail(h4.1/(h5+h6+h7+h3+h3.1+h4+h4.1),1)
pcol_oc1;pcol_fa;pcol_fl;pcol_cmeat1;pcol_vmeat1;pcol_chve1;pcol_beefve1
#For veal farmers
h8<-(-log(1-risk$r6.2)) #OC
h9<-(-log(1-risk$r6.3)) #Farmer
h10<-(-log(1-risk$r7.1)) #Veal
pcol_oc2<-tail(h8/(h8+h9+h10+h3+h3.1+h4+h4.1),1)
pcol_fa2<-tail(h9/(h8+h9+h10+h3+h3.1+h4+h4.1),1)
pcol_vc2<-tail(h10/(h8+h9+h10+h3+h3.1+h4+h4.1),1)
pcol_cmeat2<-tail(h3/(h8+h9+h10+h3+h3.1+h4+h4.1),1)
pcol_vmeat2<-tail(h3.1/(h8+h9+h10+h3+h3.1+h4+h4.1),1)
pcol_chve2<-tail(h4/(h8+h9+h10+h3+h3.1+h4+h4.1),1)
pcol_beefve2<-tail(h4.1/(h8+h9+h10+h3+h3.1+h4+h4.1),1)
pcol_oc2;pcol_fa2;pcol_vc2;pcol_cmeat2;pcol_vmeat2;pcol_chve2;pcol_beefve2
#For flocks
h11<-(-log(1-risk$ere9)) #flock
h12<-(-log(1-risk$ere10)) #farmer
pcol_fl<-tail(h11/(h11+h12),1)
pcol_fa<-tail(h12/(h11+h12),1)
pcol_fl;pcol_fa
#For herds
h10.1<-(-log(1-risk$ere9.1)) #herd
h11.1<-(-log(1-risk$ere10.1)) #farmer
pcol_herd<-tail(h10.1/(h10.1+h11.1),1)
pcol_vfa<-tail(h11.1/(h10.1+h11.1),1)
pcol_herd;pcol_vfa
long<-read_excel(here("Data","Raw","longitudinal_veal.xlsx"),sheet="Rsheet")
long$week1<-as.numeric(sapply(strsplit(long$Week, "wk"), "[", 2))
# For ESBL
ggplot(subset(long,ESBL!="none"),aes(x =factor(week1),fill=ESBL)) +
geom_bar(position="fill")+
scale_fill_tableau() +
scale_color_tableau() +
geom_text(aes(label = ..count..), stat = "count" ,position = position_fill(.5),size=3.8)+
theme(axis.text.x = element_text(angle = 0, hjust = 1, vjust = 0.5))+
xlab("Weeks") + ylab("Relative frequency")+
theme(
axis.title.x = element_text( size=13),
axis.title.y = element_text( size=13)
)+
labs(fill =substitute(paste("ESC-R",italic(' E.coli'))))
#For MLST
ggplot(subset(long,!is.na(MLST)),aes(x =factor(week1),fill=MLST)) +
geom_bar(position="fill")+
scale_fill_tableau() +
scale_color_tableau() +
geom_text(aes(label = ..count..), stat = "count" ,position = position_fill(.5),size=3.8)+
theme(axis.text.x = element_text(angle = 0, hjust = 1, vjust = 0.5))+
xlab("Weeks") + ylab("Relative frequency")+
theme(
axis.title.x = element_text( size=13),
axis.title.y = element_text( size=13)
)+
labs(fill =substitute(paste("ESC-R",italic(' E.coli'))))
# substitution
su<-long[,c(2,8,9)]%>%
group_by(Eartag,ESBL,week1)
ggplot(subset(su,ESBL!="none" & Eartag!="-"),mapping=aes(x=as.numeric(week1),y=Eartag,colour=ESBL))+
geom_point(size=2)+
geom_line(size=1.5)+
xlab("Weeks") + ylab("Eartag")+labs(colour = "ESBL")+
#scale_y_continuous(minor_breaks = seq(0, 4, 1))+
scale_x_continuous(minor_breaks = c(0,2,6,10,18,24))+
theme(text = element_text(size = 20))
# ESBL
su1<-spread(subset(su,ESBL!="none" & Eartag!="-"),key=week1,value=ESBL)
subs<-NULL
for (i in 1:dim(su1)[1]){
subs[i]<-length(table(unique(unlist(su1[i,2:7]))))
}
mean(subs>1)
mean(subset(subs, subs>1))
# SNP
su1.2<-subset(long[,c(2,7,9)],Eartag!="-" & !is.na(SNP_clones))
su2<-spread(su1.2,key=week1,value=SNP_clones)
subs2<-NULL
for (i in 1:dim(su2)[1]){
subs2[i]<-length(table(unique(unlist(su2[i,2:7]))))
}
mean(subs2>1)
mean(subset(subs2, subs2>1))
# MLST
su1.3<-subset(long[,c(2,6,9)],Eartag!="-" & MLST!="-")
su3<-spread(su1.3,key=week1,value=MLST)
subs3<-NULL
for (i in 1:dim(su3)[1]){
subs3[i]<-length(table(unique(unlist(su3[i,2:7]))))
}
mean(subs3>1)
# SNP ->ESBL
su1.4<-subset(long[,c(2,7,8,9)],!is.na(SNP_clones)&ESBL!="none"&Eartag!="-")
su4<-spread(su1.4,key=week1,value=ESBL)
subs4<-NULL
for (i in 1:dim(su4)[1]){
subs4[i]<-length(table(unique(unlist(su4[i,3:8]))))
}
mean(subs4>1)
# MLST ->ESBL
su1.5<-subset(long[,c(2,6,8,9)],MLST!="-" & ESBL!="none"& Eartag!="-")
su5<-spread(su1.5,key=week1,value=ESBL)
subs5<-NULL
for (i in 1:dim(su5)[1]){
subs5[i]<-length(table(unique(unlist(su5[i,3:8]))))
}
mean(subs5>1)
# SNP vs ESBL
table(long$SNP_clones,long$ESBL)
# MLST vs ESBL
table(long$MLST,long$ESBL)
# MLST vs SNP
table(long$MLST,long$SNP_clones)
library(ggplot2)
library(gridExtra)
library(tidyverse)
pop1<-list()
pop2<-list()
p1<-list()
p2<-list()
r1<-NULL
r2<-NULL
r3<-NULL
r4<-NULL
wp1<-NULL
wp2<-NULL
wp3<-NULL
wp4<-NULL
pop1[[1]]<-c(0,0.25,0.75)
pop2[[1]]<-c(0,0.5,0.5)
r1[1]<-1-exp(-(1-pop1[[1]][1]))
r2[1]<-1-exp(-(1-pop2[[1]][1]))
r3[1]<-1-exp(-(1-pop2[[1]][1]))
r4[1]<-1-exp(-(1-pop1[[1]][1]))
if(r1[1]+r2[1]==0){wp1[1]=wp2[1]=0}else{
wp1[1]<-r1[1]/(r1[1]+r2[1])*(1-(1-r1[1])*(1-r2[1]))
wp2[1]<-r2[1]/(r1[1]+r2[1])*(1-(1-r1[1])*(1-r2[1]))
}
if(r3[1]+r4[1]==0){wp3[1]=wp4[1]=0}else{
wp3[1]<-r3[1]/(r3[1]+r4[1])*(1-(1-r3[1])*(1-r4[1]))
wp4[1]<-r4[1]/(r3[1]+r4[1])*(1-(1-r3[1])*(1-r4[1]))
}
p1[[1]]=rbind(c((1-sum(wp1[1]*pop1[[1]][2:3]+wp2[1]*pop2[[1]][2:3])), wp1[1]*pop1[[1]][2:3]+wp2[1]*pop2[[1]][2:3] ),
c(0.75,0.25,0),
c(0.6,0,0.4)
)
p2[[1]]=rbind(c((1-sum(wp3[1]*pop2[[1]][2:3]+wp4[1]*pop1[[1]][2:3])), wp3[1]*pop2[[1]][2:3]+wp4[1]*pop1[[1]][2:3] ),
c(0.2,0.8,0),
c(0.3,0,0.7)
)
for(i in 2:100){
r1[i]<-1-exp(-(1-pop1[[i-1]][1]))
r2[i]<-1-exp(-(1-pop2[[i-1]][1]))
if(r1[i]+r2[i]==0){wp1[i]=wp2[i]=0}else{
wp1[i]<-r1[i]/(r1[i]+r2[i])*(1-(1-r1[1])*(1-r2[1]))
wp2[i]<-r2[i]/(r1[i]+r2[i])*(1-(1-r1[1])*(1-r2[1]))
}
p1[[i]]=rbind(c((1-sum(wp1[i-1]*pop1[[i-1]][2:3]+wp2[1]*pop2[[i-1]][2:3])), wp1[i-1]*pop1[[i-1]][2:3]+wp2[i-1]*pop2[[i-1]][2:3] ),
c(0.2,0.8,0),
c(0.3,0,0.7)
)
pop1[[i]]<-pop1[[i-1]]%*%p1[[1]]
r3[i]<-1-exp(-(1-pop2[[i-1]][1]))
r4[i]<-1-exp(-(1-pop1[[i-1]][1]))
if(r3[i]+r4[i]==0){wp3[i]=wp4[i]=0}else{
wp3[i]<-r3[i]/(r3[1]+r4[i])*(1-(1-r3[1])*(1-r4[1]))
wp4[i]<-r4[i]/(r3[1]+r4[i])*(1-(1-r3[1])*(1-r4[1]))
}
p2[[i]]=rbind(c((1-sum(wp3[i-1]*pop2[[i-1]][2:3]+wp4[i-1]*pop1[[i-1]][2:3])), wp3[i-1]*pop2[[i-1]][2:3]+wp4[i-1]*pop1[[i-1]][2:3] ),
c(0.75,0.25,0),
c(0.6,0,0.4)
)
pop2[[i]]<-pop2[[i-1]]%*%p2[[1]]
}
#Organizing data frames
population1<-as.data.frame(matrix(unlist(pop1),byrow = T,ncol=3,dimnames = list(NULL, c("Neg", "G1","G2"))))
population1$time<-1:100
population2<-as.data.frame(matrix(unlist(pop2),byrow = T,ncol=3,dimnames = list(NULL, c("Neg", "G1","G2"))))
population2$time<-1:100
#Tidy data to plot the results
plot_pop1<-population1%>%
gather(key="Status",value="Prev",-time)
plot_pop2<-population2%>%
gather(key="Status",value="Prev",-time)
#Plot the results
plot1<-ggplot() + geom_area(aes(y = Prev, x = time, fill = Status), data = plot_pop1,
stat="identity")+
theme_minimal()+
labs(title = "Population 1")
plot2<-ggplot() + geom_area(aes(y = Prev, x = time, fill = Status), data = plot_pop2,
stat="identity")+
theme_minimal()+
labs(title = "Population 2")
grid.arrange(plot1, plot2, ncol=1,nrow = 2)
\ No newline at end of file
......@@ -208,7 +208,7 @@ days_ab2<-read_excel(here("Data","Raw","Group_abuse2.xlsx"),sheet="#Times")
# Descriptive of the animals and atb
```{r}
veal2<-subset(veal1,time==1 & `Maldi-TOF` %in% c(0,1))
veal2<-subset(veal1,time==1)
##Population description
......
......@@ -36,8 +36,7 @@
#Script start#############################################################################
library(ggplot2)
library(tidyverse)
# Final risk/week
risk <- read.csv(here("Output","chicken_risk","scenario1","data1.txt"), sep="")
......@@ -49,6 +48,8 @@ risk <- read.csv(here("Output","chicken_risk","scenario1","data1.txt"), sep="")
1-((1-risk$ere9[dim(risk)[1]])*(1-risk$ere10[dim(risk)[1]]))
# Table with source attribution
rbind(
cbind(tail(risk$r1,n=1),tail(risk$r2,n=1),"NA",tail(risk$r3,n=1),tail(risk$r4,n=1)),
......@@ -58,10 +59,131 @@ cbind("NA",tail(risk$ere10,n=1),tail(risk$ere9,n=1),"NA","NA")
)
##Dataset with the sources according Thomas's eq
#For OC
h1<-(-log(1-risk$r1)) #OC
h2<-(-log(1-risk$r2)) #Farmer
h3<-(-log(1-risk$r3)) #Meat
h4<-(-log(1-risk$r4)) #Veg
pcol_oc<-tail((h1/(h1+h2+h3+h4)),1)
pcol_fa<-tail(h2/(h1+h2+h3+h4),1)
pcol_me<-tail(h3/(h1+h2+h3+h4),1)
pcol_ve<-tail(h4/(h1+h2+h3+h4),1)
pcol_oc;pcol_fa;pcol_me;pcol_ve
#For farmers
h5<-(-log(1-risk$r6)) #OC
h6<-(-log(1-risk$r6.1)) #Farmer
h7<-(-log(1-risk$r7)) #Flock
h8<-(-log(1-risk$r3)) #Meat
h9<-(-log(1-risk$r4)) #Veg
pcol_oc<-tail(h5/(h5+h6+h7+h8+h9),1)
pcol_fa<-tail(h6/(h5+h6+h7+h8+h9),1)
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)
#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)
# Graph for source attribution
attr<-cbind.data.frame(Time=1:200,OC=risk$r1,Farmer=risk$r2,Meat=risk$r3,Veg=risk$r4)
attr.1<-cbind.data.frame(Time=1:200,OC=risk$r6,Farmer=risk$r6.1,Flock=risk$r7,Meat=risk$r3,Veg=risk$r4)
attr.2<-cbind.data.frame(Time=1:200,Farmer=risk$ere10,Flock=risk$ere9)
attr1<-attr%>%
gather(key="Source",value="Attribution",-Time)
attr2<-attr.1%>%
gather(key="Source",value="Attribution",-Time)
attr3<-attr.2%>%
gather(key="Source",value="Attribution",-Time)
at1<-ggplot(attr1, aes(x=Time, y=Attribution, linetype=Source)) +
geom_line(size=1.01)+
ylab("Attribution")+
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,1.850138e-03,0.0006))+
scale_linetype_manual(values=c(1, 3, 5, 10),element_text(size=10))+
ggtitle("Source attribution for Open community")
at2<-ggplot(attr2, aes(x=Time, y=Attribution, linetype=Source)) +
geom_line(size=1.01)+
ylab("Attribution")+
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.02,0.005))+
scale_linetype_manual(values=c(1, 3, 5, 8,10),element_text(size=10))+
ggtitle("Source attribution for Farmers")
at3<-ggplot(attr3, aes(x=Time, y=Attribution, linetype=Source)) +
geom_line(size=1.01)+
ylab("Attribution")+
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.15,0.025))+
scale_linetype_manual(values=c(1, 3),element_text(size=10))+
ggtitle("Source attribution for Flocks")
grid.arrange(at1, at2,at3, ncol=2,nrow = 2)
#Final colonization week 200
colo <- read.csv(here("Output","chicken","scenario1","data1.txt"), sep="")
base<-colo[200,1]
colo1.2<- gather(colo,key="pop",value="prop")
colo1.2$time<-rep(seq(1,200,1),4)
......@@ -89,51 +211,312 @@ colo2$pop<-as.factor(colo2$pop)
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