Commit 234f7026 authored by de Freitas Costa, Eduardo's avatar de Freitas Costa, Eduardo
Browse files

First commmit

parent 084419c9
#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
This diff is collapsed.
---
title: "Summary report for dairy farms"
author: "Eduardo de Freitas Costa"
date: '`r format(Sys.Date(), "%Y-%m-%d")`'
output:
html_document:
df_print: paged
code_folding: hide
toc: yes
toc_depth: 3
---
```{r}
#header #################################################################################
#'veal_calves.R'
#Title: MADRA
#Project ID: 1600001916
#Client: client
#Author: <Eduardo> <Costa>, Wageningen Bioveterinary Research
#Description: Factors associated to ESBL in dayry farms. Results not used
#Start date: date
#Last Update: {6:date}
#R version: r.version
#Scriptversion: version
#Dependencies
#<-Downstream
#- none
#->Upstream
#- none
#Input:
#- Data/Raw/veal_calves.xlsx
#Output:
#- none
#Peer reviewer(s)
#Please ensure directories are relative. Please comment code sufficiently.
#Script start#############################################################################
```
```{r setup, include=FALSE}
# Library and source statements ###############################################
rm(list = ls())
#Packages to be used
packages<-c("readxl","here","tidyverse","ggplot2","lme4","logistf","car","knitr","glmmsr","plotly","gridExtra","grid","ggridges","ggthemes","glmmTMB",
"bbmle","multcomp","multcompView","lsmeans","lattice")
# Install packages not yet installed
installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
install.packages(packages[!installed_packages])
}
# Packages loading
invisible(lapply(packages, library, character.only = TRUE))
# Create folders ##############################################################
dir.create(here("Figures"))
dir.create(here("Outputs"))
# Define knitr options ########################################################
# global options
knitr::opts_chunk$set(message = FALSE) # code default outputs
knitr::opts_chunk$set(warning = FALSE) # code warnings
knitr::opts_chunk$set(echo = TRUE) # code chunks show
```
```{r}
#Importing data form veal calve studies
veal<-rbind(
read_excel(here("Data","Raw","veal_calves.xlsx"),sheet = "TRNS")
)
colnames(veal)[8]<-"output"
veal<-subset(veal,`Maldi-TOF` %in% c(0,1)) #Only E. coli
```
#Descriptives
```{r}
jpeg(file=here("Figures","FigS2.jpg"), units ="in",width = 10, height = 8,res = 1000)
ggplot(veal,aes(x =factor(DF),fill=factor(output))) +
theme(legend.title = element_blank())+
scale_fill_tableau(labels=c("Negative","Positive")) +
scale_color_manual(labels=c("Negative","Positive")) +
theme(legend.text = element_text(size = 10),
legend.margin = margin(6, 2, 6, 2),legend.position = "bottom")+
geom_bar(position="fill")+
theme(plot.title = element_text(hjust = 0.5,size=12))+
ggtitle(expression(paste("Percentage of ESC-R ",italic("E. coli"), "at dairy farms")))+
labs(y = expression(paste("Percentage of ESC-R ",italic("E. coli")," (%)"), x = "Dairy farms"))+
theme(axis.title.x = element_blank())+
scale_y_continuous(labels = scales::percent_format(accuracy = 1))+
scale_x_discrete(labels=c("DF1","DF2","DF3","DF4","DF5","DF6","DF7","DF8","DF9","DF10","DF11","DF12","DF13"),element_text(size=10))
dev.off()
ggplot(veal,aes(x =factor(treatment),fill=factor(output))) +
labs(fill="ESBL")+
geom_bar(position="fill")+
scale_fill_tableau() +
scale_color_tableau() +
geom_text(aes(label = ..count..), stat = "count", position = position_fill(.5))+
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 0.5))+
ggtitle("Ages fo the veal calves")
ggplot(veal,aes(x =factor(`#Cows`),fill=factor(output))) +
labs(fill="ESBL")+
geom_bar(position="fill")+
scale_fill_tableau() +
scale_color_tableau() +
geom_text(aes(label = ..count..), stat = "count", position = position_fill(.5))+
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 0.5))+
ggtitle("Number of lactating cows at the farm")
ggplot(veal,aes(x =factor(HousingCows),fill=factor(output))) +
labs(fill="ESBL")+
geom_bar(position="fill")+
scale_fill_tableau() +
scale_color_tableau() +
geom_text(aes(label = ..count..), stat = "count", position = position_fill(.5))+
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 0.5))+
ggtitle("Housing of the cows")
ggplot(veal,aes(x =factor(OtherAnimals),fill=factor(output))) +
labs(fill="ESBL")+
geom_bar(position="fill")+
scale_fill_tableau() +
scale_color_tableau() +
geom_text(aes(label = ..count..), stat = "count", position = position_fill(.5))+
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 0.5))+
ggtitle("Presence of other animals 0=no; 1=dog(s); 2= cat(s), 3= dog(s)/cat(s) and other")
ggplot(veal,aes(x =factor(OtherJob),fill=factor(output))) +
labs(fill="ESBL")+
geom_bar(position="fill")+
scale_fill_tableau() +
scale_color_tableau() +
geom_text(aes(label = ..count..), stat = "count", position = position_fill(.5))+
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 0.5))+
ggtitle("If the farmer has another job next to taking care of the dairy farm 0=no; 1= yes ")
ggplot(veal,aes(x =factor(HealthStatus),fill=factor(output))) +
labs(fill="ESBL")+
geom_bar(position="fill")+
scale_fill_tableau() +
scale_color_tableau() +
geom_text(aes(label = ..count..), stat = "count", position = position_fill(.5))+
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 0.5))+
ggtitle("What is the health status of your farm?")
ggplot(veal,aes(x =factor(HousingDryPeriod),fill=factor(output))) +
labs(fill="ESBL")+
geom_bar(position="fill")+
scale_fill_tableau() +
scale_color_tableau() +
geom_text(aes(label = ..count..), stat = "count", position = position_fill(.5))+
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 0.5))+
ggtitle("Housing of the dry cows")
ggplot(veal,aes(x =factor(SeperateCalves),fill=factor(output))) +
labs(fill="ESBL")+
geom_bar(position="fill")+
scale_fill_tableau() +