ANALYSIS FOR PRODUCTION SYSTEM AMR DRIVERS IN UGANDA
Adrian Muwonge
2024-08-08
R Markdown
CONTEXT AND DATA DICTIONANRY
This code is written to analyse the Production associated AMR data collected in Uganda as part of BBSRC funded fellowship. There are three datasets, Metadata, Phenotype and Qpcr datases METADATA.csv has 476 rows and 14 colums
This data set is named FLF_PHENOTYPE DATA It contains 8 columns and 2830 rows This data is output from analysis of bacteria culture from farmers and thier pigs in the one year longitudinal study. We collected these sample 6 times growing three bacteria E.coli, Klebsiella and Salmonella
Pig and Human ID
SAMPLE ID E.G KLAMAK011P The first three letters District, Second three letter Subcounty The two number digit is for the household next number is for the visit number the last number is for the host (Pig or Human)
Exposure references (can be viewed as controls) KLANANCTR021 CRT means References 02 = Low exposure references (LEFs) Little to no conctact with pigs because of region or cultural taboos 01= High exposure references (HEFs) high frequency of contact such as people who eviscerate pig guts in the abattoirs The last number is for the visit number
For analysis use column 3…. SAMPLE_ID2
We often generate variables out of the ID as data intergrity validation process.. usually in a data set as ID3
PRODUCTION SYSTEM…. captured as District (Urban =Kampala Semi intensive pigger), Rural=Mubende(free range)) Subcounty….. is an administrative division in the district SAMPLING…. is the visit time (1-6) Species…… is the Host (Human or Pig) Isolate…… Is a tracking number for the laboratory Bacteria…….Is the Bacteria isolated
Antibiotics, these will have three states(R,I,S) resistant, intermediet and Susceptible Note that, we did not have results for Imipenem for the first sampling point
Reading in packages
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(ggplot2)
library(reshape)##
## Attaching package: 'reshape'
##
## The following object is masked from 'package:lubridate':
##
## stamp
##
## The following object is masked from 'package:dplyr':
##
## rename
##
## The following objects are masked from 'package:tidyr':
##
## expand, smithslibrary(reshape2)##
## Attaching package: 'reshape2'
##
## The following objects are masked from 'package:reshape':
##
## colsplit, melt, recast
##
## The following object is masked from 'package:tidyr':
##
## smithslibrary(tidyr)
library(tidyverse)
library(ggpubr)
library(lme4)## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following object is masked from 'package:reshape':
##
## expand
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacklibrary(jtools)
library(sjPlot)## #refugeeswelcomelibrary(stringr)
library(car)## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## somelibrary(broom)
library(ggeffects)##
## Attaching package: 'ggeffects'
##
## The following object is masked from 'package:jtools':
##
## johnson_neymanlibrary(jcolors)
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
library(grid)
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots = length(plots)
# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel
# ncol: Number of columns of plots
# nrow: Number of rows needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
ncol = cols, nrow = ceiling(numPlots/cols))
}
if (numPlots==1) {
print(plots[[1]])
} else {
# Set up the page
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col))
}
}
}Reading in Datasets
METADATA<- read.csv("../META_DATA.csv",header = T,sep = ",")
PHENODATA1<- read.csv("../For_submission/Datasets/PHENOTYPE_DB.csv",header = T,sep = ",")
QPCRDATA<- read.csv("../For_submission/Datasets/FLF_QPCR_DB.csv",header = T,sep = ",")
HH_data<- read.csv("../For_submission/Datasets/HOUSE_HOLD_PAIRS.csv",header = T,sep = ",")## The household pairs in this
HH_data_MDR<- read.csv("../For_submission/Datasets/HOUSE_HOLD_PAIRS_MDR_SHARED.csv",header = T,sep = ",")
Pattern_data<- read.csv("../For_submission/Datasets/E_COLI_MDR_SHARED.csv",header = T,sep = ",")
staggered<-read.csv("../For_submission/Datasets/Farmer_pig_carriage.csv",header = T,sep = ',')# names(PHENODATA1)[5]<-"Visits"
# names(PHENODATA1)[4]<-"setting"
# PHENODATA1<-select(PHENODATA1, -Amoxy_clav)
PHENODATA1<-select(PHENODATA1, -X)
Controlsq<-QPCRDATA[str_detect(QPCRDATA$SAMPLE_ID, "CTR0"),]
# #PHENODATA1<-PHENODATA1[!str_detect(PHENODATA1$SAMPLE_ID2,"CTR"),]
# PHENODATA1$setting[PHENODATA1$setting=="KAMPALA"]<-"Peri_urban"
# PHENODATA1$setting[PHENODATA1$setting=="MUBENDE"]<-"Rural"
#
# PHENODATA1$Ciproflaxacine[PHENODATA1$Ciproflaxacine=="I"]<-"S"
# PHENODATA1$Ciproflaxacine[PHENODATA1$Ciproflaxacine=="1"]<-"S"
# PHENODATA1$Ciproflaxacine[PHENODATA1$Ciproflaxacine=="SI"]<-"S"
#
# PHENODATA1$Ciproflaxacine[PHENODATA1$Ciproflaxacine=="S"]<-0
# PHENODATA1$Ciproflaxacine[PHENODATA1$Ciproflaxacine=="R"]<-1
#
#
# PHENODATA1$Sulfa_Trimethoprim[PHENODATA1$Sulfa_Trimethoprim=="I"]<-"S"
# PHENODATA1$Sulfa_Trimethoprim[PHENODATA1$Sulfa_Trimethoprim=="1"]<-"S"
# PHENODATA1$Sulfa_Trimethoprim[PHENODATA1$Sulfa_Trimethoprim=="SI"]<-"S"
#
# PHENODATA1$Sulfa_Trimethoprim[PHENODATA1$Sulfa_Trimethoprim=="S"]<-0
# PHENODATA1$Sulfa_Trimethoprim[PHENODATA1$Sulfa_Trimethoprim=="R"]<-1
#
#
# PHENODATA1$Strepptomycine[PHENODATA1$Strepptomycine=="I"]<-"S"
# PHENODATA1$Strepptomycine[PHENODATA1$Strepptomycine=="1"]<-"S"
# PHENODATA1$Strepptomycine[PHENODATA1$Strepptomycine=="SI"]<-"S"
#
# PHENODATA1$Strepptomycine[PHENODATA1$Strepptomycine=="S"]<-0
# PHENODATA1$Strepptomycine[PHENODATA1$Strepptomycine=="R"]<-1
#
# PHENODATA1$Nalidixic_acid[PHENODATA1$Nalidixic_acid=="I"]<-"S"
# PHENODATA1$Nalidixic_acid[PHENODATA1$Nalidixic_acid=="1"]<-"S"
# PHENODATA1$Nalidixic_acid[PHENODATA1$Nalidixic_acid=="SI"]<-"S"
#
# PHENODATA1$Nalidixic_acid[PHENODATA1$Nalidixic_acid=="S"]<-0
# PHENODATA1$Nalidixic_acid[PHENODATA1$Nalidixic_acid=="R"]<-1
#
#
# PHENODATA1$Tetracycline[PHENODATA1$Tetracycline=="I"]<-"S"
# PHENODATA1$Tetracycline[PHENODATA1$Tetracycline=="1"]<-"S"
# PHENODATA1$Tetracycline[PHENODATA1$Tetracycline=="SI"]<-"S"
#
# PHENODATA1$Tetracycline[PHENODATA1$Tetracycline=="S"]<-0
# PHENODATA1$Tetracycline[PHENODATA1$Tetracycline=="R"]<-1
#
#
# PHENODATA1$Gentamycine[PHENODATA1$Gentamycine=="I"]<-"S"
# PHENODATA1$Gentamycine[PHENODATA1$Gentamycine=="1"]<-"S"
# PHENODATA1$Gentamycine[PHENODATA1$Gentamycine=="SI"]<-"S"
#
# PHENODATA1$Gentamycine[PHENODATA1$Gentamycine=="S"]<-0
# PHENODATA1$Gentamycine[PHENODATA1$Gentamycine=="R"]<-1
#
#
# PHENODATA1$Chloroamphenicol[PHENODATA1$Chloroamphenicol=="I"]<-"S"
# PHENODATA1$Chloroamphenicol[PHENODATA1$Chloroamphenicol=="1"]<-"S"
# PHENODATA1$Chloroamphenicol[PHENODATA1$Chloroamphenicol=="SI"]<-"S"
#
# PHENODATA1$Chloroamphenicol[PHENODATA1$Chloroamphenicol=="S"]<-0
# PHENODATA1$Chloroamphenicol[PHENODATA1$Chloroamphenicol=="R"]<-1
#
# PHENODATA2<-melt(PHENODATA1, id.vars = c("Sample.ID","Bacteria","Visits","Species","setting"), measure.vars = 9:16)
#
# names(PHENODATA2)[5]<-"setting"
# names(PHENODATA2)[6]<-"Antibiotics"
# names(PHENODATA2)[7]<-"Resistance"
# names(PHENODATA2)[3]<-"Visits"
# PHENODATA2$setting[PHENODATA2$setting=="KAMPALA"]<-"Peri_urban"
# PHENODATA2$setting[PHENODATA2$setting!="KAMPALA"]<-"Rural"
# PHENODATA2$Resistance[PHENODATA2$Resistance=="I"]<-"S"
# PHENODATA2$Resistance[PHENODATA2$Resistance=="1"]<-"S"
# PHENODATA2$Resistance[PHENODATA2$Resistance=="SI"]<-"S"
# PHENODATA2<-PHENODATA2[PHENODATA2$Resistance!="",]
#
#
# PHENODATA1<- PHENODATA1[,c(2:15),]
# names(PHENODATA1)[2]<-"sampleid"
#
PHENODATA1B<-melt(PHENODATA1, id.vars = c("sampleid","Bacteria","Visits","Species","setting"), measure.vars = 8:14)
names(PHENODATA1B)[5]<-"setting"
names(PHENODATA1B)[6]<-"Antibiotics"
names(PHENODATA1B)[7]<-"Resistance"
names(PHENODATA1B)[3]<-"Visits"
PHENODATA1B<-PHENODATA1B[PHENODATA1B$Antibiotics!="Amoxy_clav",]
PHENODATA1B$setting[PHENODATA1B$setting=="KAMPALA"]<-"Peri_urban"
PHENODATA1B$setting[PHENODATA1B$setting=="MUBENDE"]<-"Rural"
PHENODATA1B<-PHENODATA1B[PHENODATA1B$Resistance!="",]SUMMARY TABLE
DAT1<-METADATA[METADATA$experiment!="HUMAN",] ### keep only References
DAT1<-DAT1[DAT1$experiment!="PIGS",] ## Keep only References
DAT1<-droplevels.data.frame(DAT1) ## Keep only References
PIGZ<-METADATA[METADATA$experiment!="HUMAN",]
PIGZ<-PIGZ[PIGZ$experiment!="NEGATIVE CONTROL",]
PIGZ<-PIGZ[PIGZ$experiment!="POSITIVE CONTROL",]
CTRS<-METADATA[METADATA$experiment!="HUMAN",]
CTRS<-CTRS[CTRS$experiment!="PIGS",]
META_PGHU<-METADATA[!str_detect(METADATA$sampleid,"CTR"),]
CTRS<- CTRS[,c(1:9,13:17),]
CTRS$hh<-substr(CTRS$sampleid,1,11)
CTRS$experiment[CTRS$experiment=="NEGATIVE CONTROL"]<-"Negative_control" ## Low exposure reference
CTRS$experiment[CTRS$experiment=="POSITIVE CONTROL"]<-"Positive_control" ## Highe exposure reference
CTRS_pheno<-PHENODATA1[PHENODATA1$Species!="Human",]
CTRS_pheno<-CTRS_pheno[CTRS_pheno$Species!="Pig",]
names(CTRS_pheno)[3]<-"Setting"
CTR_META<-CTRS_pheno%>% inner_join(CTRS,by="sampleid")# joining
CTR_META<-CTR_META[CTR_META$Bacteria!="Salmonella",]
table(CTR_META$Bacteria,CTR_META$maritat_status)##
## divorced married single widow
## E.coli 0 31 8 1
## Klebsiella 2 19 3 2names(PHENODATA1)[3]<-"Setting"
Phe_META<-PHENODATA1%>% inner_join(META_PGHU,by="sampleid")
Phe_META<-Phe_META[Phe_META$Bacteria!="Salmonella",]
PHENODATA1_pig<- PHENODATA1[PHENODATA1$Species=="Pig",]
Phe_META_pig<-PHENODATA1_pig%>% inner_join(PIGZ,by="sampleid")
Phe_META_pig<-Phe_META_pig[Phe_META_pig$Bacteria!="Salmonella",]
## FARMER SUMMARY STATS
table(Phe_META[Phe_META$Species=="Human",]$Bacteria,Phe_META[Phe_META$Species=="Human",]$maritat_status)##
## divorced married single widow
## E.coli 5 130 41 7
## Klebsiella 5 72 25 8 table(Phe_META[Phe_META$Species=="Human",]$Bacteria,Phe_META[Phe_META$Species=="Human",]$education)##
## PRIMARY SECONDARY TERTIARY
## E.coli 69 89 25
## Klebsiella 43 54 13 table(Phe_META[Phe_META$Species=="Human",]$Bacteria,Phe_META[Phe_META$Species=="Human",]$smoker_status)##
## non-smoker smoker(BID) smoker(SID) smoker(TID or more)
## E.coli 175 2 5 1
## Klebsiella 106 1 3 0 table(Phe_META[Phe_META$Species=="Human",]$Bacteria,Phe_META[Phe_META$Species=="Human",]$alcohol_status)##
## NO YES
## E.coli 135 48
## Klebsiella 84 26 table(Phe_META[Phe_META$Species=="Human",]$Bacteria,Phe_META[Phe_META$Species=="Human",]$medication_human)##
## NO YES
## E.coli 146 37
## Klebsiella 86 24 table(Phe_META$Bacteria,Phe_META$medication_pigs)##
## NO YES
## E.coli 132 51
## Klebsiella 79 31 table(Phe_META$Bacteria,Phe_META$breed_verif)##
## EXOTIC_BREED LOCAL MIXED_BREED
## E.coli 33 33 117
## Klebsiella 22 25 63 table(Phe_META$Bacteria,Phe_META$pig_housing)##
## INDOOR OUTDOOR
## E.coli 145 38
## Klebsiella 88 22 table(Phe_META$Bacteria,Phe_META$pig.breeds)##
## EXOTIC_BREED LOCAL MIXED_BREED
## E.coli 33 45 105
## Klebsiella 22 24 64 table(Phe_META$Bacteria,Phe_META$cleaning_frequency)##
## DAILY NEVER TWICE_AWEEK WEEKLY
## E.coli 109 13 24 37
## Klebsiella 65 11 14 20 ## PIG SUMMARY STATS
PIGZ$MERGE_ID<-substr(PIGZ$sampleid,1,9)
METADATA$MERGE_ID<-substr(METADATA$sampleid,1,9)
PIGZ$medication_pigs<- METADATA[match(PIGZ$MERGE_ID,METADATA$MERGE_ID),]$medication_pigs
PIGZ$breed_verif<- METADATA[match(PIGZ$MERGE_ID,METADATA$MERGE_ID),]$breed_verif
PIGZ$vaccination<- METADATA[match(PIGZ$MERGE_ID,METADATA$MERGE_ID),]$vaccination
PIGZ$cleaning.frequency<- METADATA[match(PIGZ$MERGE_ID,METADATA$MERGE_ID),]$cleaning_frequency
PIGZ$pig.feed<- METADATA[match(PIGZ$MERGE_ID,METADATA$MERGE_ID),]$pig_feed
PIGZ$pig.housing<- METADATA[match(PIGZ$MERGE_ID,METADATA$MERGE_ID),]$pig_housing
table(PIGZ$medication_pigs,PIGZ$sampling.point)##
## 1 2 3 4 5 6
## NO 23 29 45 43 34 28
## YES 39 25 13 12 11 10 table(PIGZ$breed_verif,PIGZ$sampling.point)##
## 1 2 3 4 5 6
## EXOTIC_BREED 13 15 9 10 5 4
## LOCAL 8 4 11 11 9 9
## MIXED_BREED 41 35 38 34 31 25 table(PIGZ$pig.housing,PIGZ$sampling.point)##
## 1 2 3 4 5 6
## INDOOR 49 46 43 43 34 31
## OUTDOOR 13 8 15 12 11 7 table(PIGZ$cleaning.frequency,PIGZ$sampling.point)##
## 1 2 3 4 5 6
## DAILY 31 36 28 36 25 16
## NEVER 11 6 9 1 5 0
## TWICE_AWEEK 10 8 11 5 7 14
## WEEKLY 10 4 10 13 8 8 table(PIGZ$pig.feed,PIGZ$sampling.point)##
## 1 2 3 4 5 6
## COMMERCIAL_FEED 2 6 5 5 1 0
## MIXED 53 42 48 49 44 37
## SWILL 7 6 5 1 0 1## RERERENCE SUMMARY STATS
table(CTR_META$Setting,CTR_META$maritat_status)##
## divorced married single widow
## Peri_urban 0 25 11 0
## Rural 2 25 0 3 table(CTR_META$Bacteria,CTR_META$education)##
## NONE PRIMARY SECONDARY TERTIARY
## E.coli 2 14 13 11
## Klebsiella 1 9 12 4 table(CTR_META$Setting,CTR_META$smoker_status)##
## non-smoker smoker(BID)
## Peri_urban 36 0
## Rural 25 5 table(CTR_META$Setting,CTR_META$alcohol_status)##
## NO YES
## Peri_urban 9 27
## Rural 22 8 table(CTR_META$Setting,CTR_META$medication_human)##
## NO YES
## Peri_urban 31 5
## Rural 26 4 table(CTR_META$Bacteria,CTR_META$Setting)##
## Peri_urban Rural
## E.coli 32 15
## Klebsiella 14 15 table(CTR_META$Setting,CTR_META$education)##
## NONE PRIMARY SECONDARY TERTIARY
## Peri_urban 3 10 15 8
## Rural 0 13 10 7summary stats for Phenotypic and genetic AMR characteristic
names(METADATA)<- tolower(names(METADATA))
names(PHENODATA1B)<-tolower(names(PHENODATA1B))
names(QPCRDATA)<-tolower(names(QPCRDATA))
##################################################################################
#creating a new object phen.
names(PHENODATA1B)[5]<-"Visit"
names(PHENODATA1B)[3]<-"Visit"
names(PHENODATA1B)[2]<-"bacteria"
names(PHENODATA1B)[4]<-"species_id"
names(PHENODATA1B)[5]<-"location_id"
PHENODATA1B<- PHENODATA1B[PHENODATA1B$bacteria!="Salmonella",]
phen<-dcast(PHENODATA1B,antibiotics+Visit+bacteria+species_id+location_id~resistance)## Using resistance as value column: use value.var to override.## Aggregation function missing: defaulting to length names(phen)[6]<-"S"
names(phen)[7]<-"R"
phen$N <- phen$R + phen$S # add column N, Sum of R & S.
phen$proportion<-phen$R/phen$N # Add column proportion .
phenx<-phen %>%
group_by(antibiotics,bacteria,Visit,species_id,location_id) %>%
summarise(med = mean(proportion))## `summarise()` has grouped output by 'antibiotics', 'bacteria', 'Visit',
## 'species_id'. You can override using the `.groups` argument. # Generating figure 3A (as bar graph)
#'## Plot of Proportion of resistance
phenx$antibiotics<- as.character(phenx$antibiotics)
phenx$antibiotics[phenx$antibiotics=="gentamycine"]<- "Gentamicin"
phenx$antibiotics[phenx$antibiotics=="Strepptomycine"]<- "Streptomycin"
phenx$antibiotics[phenx$antibiotics=="Gentamycine"]<- "Gentamicin"
phenx$antibiotics[phenx$antibiotics=="Nalidixic_acid"]<- "Nalidixic acid"
phenx$antibiotics[phenx$antibiotics=="Ciproflaxacine"]<- "Ciprofloxacin"
phenx$antibiotics[phenx$antibiotics=="Tetracycline"]<- "Tetracycline"
phenx$antibiotics[phenx$antibiotics=="Sulfa_Trimethoprim"]<- "Trimethoprim"
names(phenx)[4]<-"Host"
phenx$Host<- as.character(phenx$Host)
phenx$Host[phenx$Host=="Human"]<- "Farmer"
phenx$Host[phenx$Host=="Pig"]<- "Pig"
phenx<-phenx[!phenx$Host=="",]
phenx<-phenx[!phenx$Host=="Negative_control",]
phenx<-phenx[!phenx$Host=="Positive_control",]
names(phenx)[5]<-"Production"
phenx$Production<- as.character(phenx$Production)
phenx$Production[phenx$Production=="Rural"]<- "Free range"
phenx$Production[phenx$Production=="Peri_urban"]<- "Semi_intensive"
## Temporal relation ship
ggscatter(phenx, x = "Visit", y = "med",
color = "Production", # Points color, shape and size
add = "reg.line", # Add regressin line
palette=c("#000080","#FC4E07","#800000"),
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE,point = F, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n") ) +
facet_wrap(~antibiotics,scales = "free_y") + xlab("Visits") + ylab("Proportion of resistant bacteria")## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_cor()`).
ggscatter(phenx, x = "Visit", y = "med",
color = "Host", # Points color, shape and size
add = "reg.line", # Add regressin line
palette=c("#69b3a2", "#404080"),
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE,point = F, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n") ) +
facet_wrap(~antibiotics,scales = "free_y") + xlab("Visits") + ylab("Proportion of resistant bacteria")## Warning: Removed 1 row containing non-finite outside the scale range (`stat_smooth()`).
## Removed 1 row containing non-finite outside the scale range (`stat_cor()`).
## boxplot of proportion of resistance to each antibiotic in the two production setting
ggplot(phenx,aes(x=antibiotics, y=med,fill=Production))+
geom_boxplot()+ facet_wrap(~Host)+theme_bw()+
theme(axis.text.x = element_text(angle = 90))+xlab("Antibiotics") +
scale_fill_manual(values=c("#000080","#FC4E07","#800000")) + ylab("proportion of resistant sentinel")## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_boxplot()`).
## SPREAD DATA TO PLOT MULTI REISTANCE
PHENODATA1$household_id<-substring(PHENODATA1$sampleid, 1,8)
names(PHENODATA1)[3]<-"location_id"
names(PHENODATA1)[4]<-"Visit"
names(PHENODATA1)[5]<-"species_id"
names(PHENODATA1)[7]<-"bacteria"
PHENODATA1$Tetracycline[is.na(PHENODATA1$Tetracycline)]<-0
PHENODATA1$Tetracycline<- as.numeric(PHENODATA1$Tetracycline)
PHENODATA1$Sulfa_Trimethoprim<- as.numeric(PHENODATA1$Sulfa_Trimethoprim)
PHENODATA1$Strepptomycine<- as.numeric(PHENODATA1$Strepptomycine)
PHENODATA1$Nalidixic_acid<- as.numeric(PHENODATA1$Nalidixic_acid)
PHENODATA1$Gentamycine<- as.numeric(PHENODATA1$Gentamycine)
PHENODATA1$Ciproflaxacine<- as.numeric(PHENODATA1$Ciproflaxacine)
PHENODATA1$MDR_count<- PHENODATA1$Ciproflaxacine + PHENODATA1$Gentamycine +
PHENODATA1$Nalidixic_acid+PHENODATA1$Strepptomycine+PHENODATA1$Sulfa_Trimethoprim+
PHENODATA1$Tetracycline
PHENODATA1<-PHENODATA1[!is.na(PHENODATA1$MDR_count),]
phMDR<-dcast(PHENODATA1,bacteria+species_id+location_id~MDR_count)## Using MDR_count as value column: use value.var to override.
## Aggregation function missing: defaulting to lengthphMDRx<-phMDR
phMDRx$N <- phMDRx$"0" + phMDRx$"1" +phMDRx$"2"+phMDRx$"3"+phMDRx$"4" + phMDRx$"5" + phMDRx$"6"
phMDRx$"00" <- phMDRx$"0"/phMDRx$N
phMDRx$"01" <- phMDRx$"1"/phMDRx$N
phMDRx$"02" <- phMDRx$"2"/phMDRx$N
phMDRx$"03" <- phMDRx$"3"/phMDRx$N
phMDRx$"04" <- phMDRx$"4"/phMDRx$N
phMDRx$"05" <- phMDRx$"5"/phMDRx$N
phMDRx$"06" <- phMDRx$"6"/phMDRx$N
phMDRx<-phMDRx[phMDRx$species_id!="Positive_control",]
phMDRx<-phMDRx[phMDRx$species_id!="Negative_control",]
phMDRx<-phMDRx[phMDRx$bacteria!="Salmonella",]
phMDRx<-phMDRx[phMDRx$species_id!="",]
phMDRx1<- phMDRx[,c(1:3,12:18),]
phMDRx1_melt<-melt(phMDRx1, id.vars = c("bacteria","species_id","location_id"), measure.vars = 4:9)
names(phMDRx1_melt)[4]<-"MDR"
names(phMDRx1_melt)[5]<-"Proportion"
phMDRx1_melt$location_id[phMDRx1_melt$location_id=="Peri_urban"]<-"Semi-intensive"
phMDRx1_melt$location_id[phMDRx1_melt$location_id=="Rural"]<-"Free-range"
names(phMDRx1_melt)[3]<-"Production"
# FINAL PLOTE PUBLISHED
phMDRx1_melt$species_id[phMDRx1_melt$species_id=="Human"]<-"Farmer"
ggplot(phMDRx1_melt,aes(MDR,Proportion,fill=Production)) + geom_bar(stat = "identity",position = "dodge") +
facet_wrap(~species_id) + theme_bw() + ylab("Proportion of isolate") + xlab("Multi-resistance") +
scale_fill_manual(values=c("#FC4E07","#000080","#800000")) + theme(legend.position="bottom") +
guides(fill=guide_legend(title="Production"))
##MERGE GENO AND PHENO
PHENO_SPREAD_1<-PHENODATA1
PHENO_SPREAD_1$species_id<- as.character(PHENO_SPREAD_1$species_id)
PHENO_SPREAD_1$species_id[PHENO_SPREAD_1$species_id=="Human"]<-"H"
PHENO_SPREAD_1$species_id[PHENO_SPREAD_1$species_id=="Pig"]<-"P"
PHENO_SPREAD_1<-PHENO_SPREAD_1[!PHENO_SPREAD_1$species_id=="",]
PHENO_SPREAD_1<-PHENO_SPREAD_1[!PHENO_SPREAD_1$species_id=="Positive_control",]
names(PHENO_SPREAD_1)[7]<-"Bacteria"
names(PHENO_SPREAD_1)[3]<-"setting"
names(PHENO_SPREAD_1)[4]<-"time"
names(PHENO_SPREAD_1)[2]<-"sample_id"
names(PHENO_SPREAD_1)[5]<-"host"
PHENO_SPREAD_1$SAMP_ID<- paste(PHENO_SPREAD_1$household_id,PHENO_SPREAD_1$time,PHENO_SPREAD_1$host,sep = "")
PHENO_SPREAD_EC<-PHENO_SPREAD_1[PHENO_SPREAD_1$Bacteria=="E.coli",]
PHENO_SPREAD_KB<-PHENO_SPREAD_1[PHENO_SPREAD_1$Bacteria=="Klebsiella",]
genotype_EC<-QPCRDATA # E.coli
genotype_KB<-QPCRDATA # Klebsiella
names(genotype_KB)[2]<-"sample_id"
names(genotype_KB)[14]<-"host"
names(genotype_KB)[8]<-"setting"
names(genotype_KB)[15]<-"time"
names(genotype_EC)[2]<-"sample_id"
names(genotype_EC)[14]<-"host"
names(genotype_EC)[8]<-"setting"
names(genotype_EC)[15]<-"time"
genotype_EC$Bacteria<-"E.coli"
genotype_KB$Bacteria<-"Klebsiella"
genotype_EC$MDR_count<-PHENO_SPREAD_EC[match(genotype_EC$sample_id,PHENO_SPREAD_EC$SAMP_ID),]$MDR_count
genotype_KB$MDR_count<-PHENO_SPREAD_KB[match(genotype_KB$sample_id,PHENO_SPREAD_KB$SAMP_ID),]$MDR_count
genotype_EC<-genotype_EC[!is.na(genotype_EC$MDR_count),]
genotype_KB<-genotype_KB[!is.na(genotype_KB$MDR_count),]
# transform, dcasr abnd plot genes
phMDR<-dcast(genotype_EC,Bacteria+host+time~MDR_count)## Using MDR_count as value column: use value.var to override.## Aggregation function missing: defaulting to lengthGeno_db_melt<-melt(genotype_EC, id.vars = c("sample_id","Bacteria","setting","host","time","MDR_count"), measure.vars = 9:12)
Geno_db_meltkb<-melt(genotype_KB, id.vars = c("sample_id","Bacteria","setting","host","time","MDR_count"), measure.vars = 9:12)
names(Geno_db_melt)[7]<-"Gene"
names(Geno_db_melt)[8]<-"Abundance"
names(Geno_db_meltkb)[7]<-"Gene"
names(Geno_db_meltkb)[8]<-"Abundance"
# Numeric value of the visit
Geno_db_melt$time<- as.character(Geno_db_melt$time)
Geno_db_melt$time[Geno_db_melt$time=="ONE"]<-1
Geno_db_melt$time[Geno_db_melt$time=="TWO"]<-2
Geno_db_melt$time[Geno_db_melt$time=="THREE"]<-3
Geno_db_melt$time[Geno_db_melt$time=="FOUR"]<-4
Geno_db_melt$time[Geno_db_melt$time=="FIVE"]<-5
Geno_db_melt$time[Geno_db_melt$time=="SIX"]<-6
Geno_db_meltkb$time<- as.character(Geno_db_meltkb$time)
Geno_db_meltkb$time[Geno_db_meltkb$time=="ONE"]<-1
Geno_db_meltkb$time[Geno_db_meltkb$time=="TWO"]<-2
Geno_db_meltkb$time[Geno_db_meltkb$time=="THREE"]<-3
Geno_db_meltkb$time[Geno_db_meltkb$time=="FOUR"]<-4
Geno_db_meltkb$time[Geno_db_meltkb$time=="FIVE"]<-5
Geno_db_meltkb$time[Geno_db_meltkb$time=="SIX"]<-6
Geno_db_melt$Genes<-sapply(strsplit(as.character(Geno_db_melt$Gene), "\\_N"), function(oo) oo[1])
Geno_db_meltkb$Genes<-sapply(strsplit(as.character(Geno_db_meltkb$Gene), "\\_N"), function(oo) oo[1])
Geno_db_melt<- Geno_db_melt[!is.na(Geno_db_melt$host),]
Geno_db_meltkb<- Geno_db_meltkb[!is.na(Geno_db_meltkb$host),]
Geno_db<- rbind(Geno_db_melt,Geno_db_meltkb)
Geno_db$setting[Geno_db$setting=="Rural"]<-"Free-range"
Geno_db$setting[Geno_db$setting=="Urban"]<-"Semi-intensive"
names(Geno_db)[3]<-"Production"
Geno_db$lgAbund<-log10(Geno_db$Abundance)
### NEW PLOT MDR & GENE
# Renaming
Geno_db$Genes<- as.character(Geno_db$Genes)
Geno_db$Genes[Geno_db$Genes=="dfra_1_norm"]<-"dfrA1"
Geno_db$Genes[Geno_db$Genes=="erm_b_norm"]<-"ermB"
Geno_db$Genes[Geno_db$Genes=="teb_b_norm"]<-"tetB"
Geno_db$Genes[Geno_db$Genes=="tet_q_norm"]<-"tetQ"
names(Geno_db)[3]<-"Production system"
Geno_db$`Production system`[Geno_db$`Production system`=="Rural"]<-"Free range"
Geno_db$`Production system`[Geno_db$`Production system`=="Urban"]<-"Semi-intensive"
ggscatter(Geno_db, x = "MDR_count", y = "lgAbund",
color = "Production system", # Points color, shape and size
add = "reg.line", # Add regressin line
palette=c("#000080","#FC4E07","#800000"),
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE,point = F, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n")) +
facet_wrap(~Genes,scales = "free_y") + xlab("Multi-antibiotic count") + ylab("log10(Normalised copies)")## Warning: Removed 100 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 100 rows containing non-finite outside the scale range
## (`stat_cor()`).
ggscatter(Geno_db, x = "MDR_count", y = "lgAbund",
color = "host", # Points color, shape and size
add = "reg.line", # Add regressin line
palette=c("#69b3a2", "#404080"),
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE,point = F, # Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 3, label.sep = "\n") ) +
facet_wrap(~Genes,scales = "free_y") + xlab("Multi-antibiotic count") + ylab("log10(Normalised copies)")## Warning: Removed 100 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 100 rows containing non-finite outside the scale range (`stat_cor()`).
## TEMPORAL SIGNAL OF VARIATION FOR AMR GENES
Geno_db$time<-as.numeric(Geno_db$time)
ggscatter(Geno_db, x = "time", y = "lgAbund",
color = "Production system", # Points color, shape and size
add = "reg.line", # Add regressin line
palette=c("#000080","#FC4E07","#800000"),
conf.int = TRUE, # Add confidence interval
cor.coef = TRUE, point = F,
# Add correlation coefficient. see ?stat_cor
cor.coeff.args = list(method = "pearson", label.x = 0.1, label.sep = "\n")
) + facet_wrap(~Genes,scales = "free_y") + ylab("log10(Normalised copies)") + xlab("Farm visits")## Warning: Removed 100 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 100 rows containing non-finite outside the scale range (`stat_cor()`).
ggplot(Geno_db,aes (x=time, y=log10(Abundance), group=host,color=host)) +
geom_smooth(method = "lm") + theme_bw()+
theme(axis.text.x = element_text(angle = 0))+ylab("log10(Normalised gene copies)") +
scale_color_manual(values=c("#69b3a2", "#404080")) +
facet_wrap(~Genes,scales = "free_y") + xlab("Visits")## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 100 rows containing non-finite outside the scale range
## (`stat_smooth()`).
my_comparisons <- list( c("0", "1" ), c("2" , "3"), c("0", "3") ) ## comparisons between MDR counts
names(Geno_db_melt)[7]<-"AMRgenes"
names(Geno_db_melt)[8]<-"Normalized_Genecopycount"
## What is the relatioship between gene carriage and multi-resistance?
Geno_db_melt$Genes<- as.character(Geno_db_melt$Genes)
Geno_db_melt$Genes[Geno_db_melt$Genes=="dfra_1_norm"]<-"dfrA1"
Geno_db_melt$Genes[Geno_db_melt$Genes=="erm_b_norm"]<-"ermB"
Geno_db_melt$Genes[Geno_db_melt$Genes=="teb_b_norm"]<-"tetB"
Geno_db_melt$Genes[Geno_db_melt$Genes=="tet_q_norm"]<-"tetQ"
ggplot(Geno_db_melt,aes(factor(MDR_count),log10(Normalized_Genecopycount),fill=host)) + geom_boxplot() + theme_bw() +
ylab("log10(Gene copy count)") + xlab("Multi-resistance") +
stat_compare_means(comparisons = my_comparisons) +
scale_fill_manual(values=c("#69b3a2", "#404080")) + facet_wrap(~Genes)## Warning: Removed 64 rows containing non-finite outside the scale range
## (`stat_boxplot()`).## Warning: Removed 64 rows containing non-finite outside the scale range
## (`stat_signif()`).
Geno_db_melt$loggene<-log10(Geno_db_melt$Normalized_Genecopycount)
ggscatter(Geno_db_melt, x = "MDR_count", y = "loggene", add = "reg.line", conf.int = T,point = F,
color="host",size=4, alpha=0.7) +
ggpubr::stat_cor(label.x = 2, label.y = log10(4)) + ## Change the label.x value til u get a plot
scale_color_manual(values=c("#69b3a2", "#404080")) +
xlab("MDR count") + ylab("log10(Normalised gene copy)") + facet_wrap(~Genes,scales = "free_y")## Warning: Removed 64 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 64 rows containing non-finite outside the scale range
## (`stat_cor()`).
Geno_db_melt$loggene<-log10(Geno_db_melt$Normalized_Genecopycount)
ggpubr::ggscatter(Geno_db_melt, x = "MDR_count", y = "loggene", add = "reg.line", conf.int = T,point = F,
color="setting",size=4, alpha=0.7) +
ggpubr::stat_cor(label.x = 1, label.y = log10(9)) + ## Change the label.x value til u get a plot
scale_color_manual(values=c("#000080","#FC4E07","#800000")) +
xlab("MDR count") + ylab("log10(Normalised gene copies") + facet_wrap(~Genes,scales = "free_y") ## Warning: Removed 64 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 64 rows containing non-finite outside the scale range (`stat_cor()`).
examine associations with Phenotypic and genetics AMR
## MODELLING MONORESISTANCE
# We Combine phenotypic data and meta data, then fit models
#Rename phenotype object, df_phenotype
PHENO_SPREAD_1$Sample_ID<- paste(PHENO_SPREAD_1$household_id,PHENO_SPREAD_1$Visit,PHENO_SPREAD_1$species_id,sep = "")
phenotype <- PHENO_SPREAD_1
names(phenotype)[5]<-"species_id"
names(phenotype)[4]<-"Visit"
met_df<-METADATA
#characters of id2
names(met_df)[7]<-"visit"
met_df$hh_time<-paste(met_df$hh,met_df$visit, sep="" ) # add hh_time column
met_df$species_id[met_df$host=="HUMAN"]<-"H"
met_df$species_id[met_df$host=="PIGS"]<-"P"
#
#made up of hh(household id) and sampling time.
# combining phenotypic data with meta data , first we split phenotype data
# to create a data set for each species.
pig_phenotype<- phenotype%>%filter(species_id=="P") # to create pig
#phenotype object
human_phenotype<- phenotype%>%filter(species_id=="H")# to create human
# phenotype data object .
# human_phenotype$species2_id <- sub("human", "H", human_phenotype$species_id) # to
# add species2_id composed of H for human.
met_df$id2<-paste(met_df$hh_time,met_df$species_id,sep = "")
human_phenotype$id2<-paste(human_phenotype$household_id,
human_phenotype$Visit,
human_phenotype$species_id, sep="" )# to
human_phenotype$id2<-paste(human_phenotype$household_id,
human_phenotype$Visit,
human_phenotype$species_id, sep="" )# to
#add column id2 , so that it's the same as in met-df to be used for merging.
combined_df<-human_phenotype%>% inner_join(met_df,by="id2") # joining the two## Warning in inner_join(., met_df, by = "id2"): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 162 of `x` matches multiple rows in `y`.
## ℹ Row 197 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.#data frames by column id2.
combined_df$hh_time<-paste(combined_df$household_id,
combined_df$Visit, sep="" ) # add hh_time
# to be used in merging .
pig_phenotype$hh_time <-paste(pig_phenotype$household_id,
pig_phenotype$Visit, sep="" )# to
#add hh_time column to be used in merging .
combined_df_pig<-pig_phenotype%>% inner_join(met_df,by="hh_time")# joining## Warning in inner_join(., met_df, by = "hh_time"): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 197 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.#pig_phenotype with meta data via column hh_time.
#combined_df_pig<-combined_df_pig[,c(1:22),] # select column 1:22 to be included
# in the data frame .
combined_df_pig$hh_time<- paste(combined_df_pig$household_id,
combined_df_pig$Visit,sep = "") #adding
#Remember the metadata from the pig was asked from the farmer....
# no talking pigs
combined_df_pig<-combined_df_pig[!is.na(combined_df_pig$alcohol_status),]
#names(combined_df_pig)[15]<-"id2"
combined_df_pig<-combined_df_pig[names(combined_df)] ## sort dataframes
# the merging .
final_df<-rbind(combined_df,combined_df_pig)# construct the final data object
#with phenotypic and meta data MELT BY ANTIBIOTICS AND REISTANCE STATUS
final_df_melt<-melt(final_df, id.vars = c(1:7,15,16:40), measure.vars = 8:14)
names(final_df_melt)[34]<-"Antibtiotic"
names(final_df_melt)[35]<-"resistance"
final_df_melt$resistance<- as.character(final_df_melt$resistance)
final_df_melt$resistance[final_df_melt$resistance=="0"]<-"S"
final_df_melt$resistance[final_df_melt$resistance=="1"]<-"R"
final_df_melt$resistance<- factor(final_df_melt$resistance, levels=c("S","R"))
# defining variables to be considered as factors .
final_df_melt$Visit<-as.factor(final_df_melt$Visit)
names(final_df_melt)[3]<-"location_id"
names(final_df_melt)[5]<-"species_id"
final_df_melt$species_id[final_df_melt$species_id=="H"]<-"Farmer"
final_df_melt$species_id[final_df_melt$species_id=="P"]<-"Pig"
final_df_melt$location_id<-as.factor(final_df_melt$location_id)
final_df_melt$Antibtiotic<-as.factor(final_df_melt$Antibtiotic)
final_df_melt$species_id<-as.factor(final_df_melt$species_id)
final_df_melt$resistance<-as.factor(final_df_melt$resistance)
final_df_melt<- final_df_melt[!final_df_melt$Bacteria=="Salmonella",]
## CREATING TABLE 2 This model shows exploration of host and production system
table(final_df_melt$Bacteria,final_df_melt$maritat_status)##
## divorced married single widow
## E.coli 42 1890 595 112
## Klebsiella 56 938 315 98table(final_df_melt$resistance,final_df_melt$location_id)##
## Peri_urban Rural
## S 1718 1348
## R 704 318md_phenotype<-glmer(data = final_df_melt,resistance~location_id+species_id+
(1|household_id)+(1|Visit)+(1|Antibtiotic),
family = binomial,control = glmerControl(optimizer ="bobyqa"))
summary(md_phenotype)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: resistance ~ location_id + species_id + (1 | household_id) +
## (1 | Visit) + (1 | Antibtiotic)
## Data: final_df_melt
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 3963.0 4000.9 -1975.5 3951.0 4082
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7081 -0.5390 -0.3450 0.0296 6.1845
##
## Random effects:
## Groups Name Variance Std.Dev.
## household_id (Intercept) 0.1096 0.3311
## Antibtiotic (Intercept) 0.8761 0.9360
## Visit (Intercept) 0.1483 0.3851
## Number of obs: 4088, groups: household_id, 70; Antibtiotic, 7; Visit, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.00636 0.39664 -2.537 0.01117 *
## location_idRural -0.69694 0.11768 -5.922 3.18e-09 ***
## species_idPig -0.21207 0.08151 -2.602 0.00927 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lctn_R
## locatn_dRrl -0.119
## species_dPg -0.095 0.018plot_model(md_phenotype, type = "est") + theme_bw()
plot_model(md_phenotype, type = "eff", terms = c("location_id","species_id"),colors = c("#69b3a2", "#404080")) +
theme_bw() + theme(axis.text.x = element_text(angle = 90, hjust = 1,size = 10))
##################################################################################
# If there's a difference in levels of resistance based on the location, does it differ based on bacteria?
#If there's a difference in levels of resistance based on the host, does it differ based on bacteria?
#if there's difference in levels of resistance based on host and location , does that differ by bacteria?
md_phenotype1<-glmer(data = final_df_melt,resistance~location_id*Bacteria+species_id*Bacteria+location_id*species_id*Bacteria+
(1|household_id)+(1|Visit)+(1|Antibtiotic),
family = binomial,control = glmerControl(optimizer ="bobyqa"))
summ(md_phenotype1)## MODEL INFO:
## Observations: 4088
## Dependent Variable: resistance
## Type: Mixed effects generalized linear regression
## Error Distribution: binomial
## Link function: logit
##
## MODEL FIT:
## AIC = 3890.75, BIC = 3960.22
## Pseudo-R² (fixed effects) = 0.06
## Pseudo-R² (total) = 0.31
##
## FIXED EFFECTS:
## -------------------------------------------------------------------------------
## Est. S.E. z val.
## ------------------------------------------------------- ------- ------ --------
## (Intercept) -0.87 0.40 -2.15
## location_idRural -0.63 0.17 -3.82
## BacteriaKlebsiella -0.74 0.16 -4.77
## species_idPig -0.05 0.12 -0.43
## location_idRural:BacteriaKlebsiella 0.59 0.24 2.42
## BacteriaKlebsiella:species_idPig -0.22 0.24 -0.91
## location_idRural:species_idPig -0.21 0.21 -1.01
## location_idRural:BacteriaKlebsiella:species_idPig -0.67 0.39 -1.73
## -------------------------------------------------------------------------------
##
## --------------------------------------------------------------
## p
## ------------------------------------------------------- ------
## (Intercept) 0.03
## location_idRural 0.00
## BacteriaKlebsiella 0.00
## species_idPig 0.67
## location_idRural:BacteriaKlebsiella 0.02
## BacteriaKlebsiella:species_idPig 0.36
## location_idRural:species_idPig 0.31
## location_idRural:BacteriaKlebsiella:species_idPig 0.08
## --------------------------------------------------------------
##
## RANDOM EFFECTS:
## ----------------------------------------
## Group Parameter Std. Dev.
## -------------- ------------- -----------
## household_id (Intercept) 0.34
## Antibtiotic (Intercept) 0.96
## Visit (Intercept) 0.36
## ----------------------------------------
##
## Grouping variables:
## --------------------------------
## Group # groups ICC
## -------------- ---------- ------
## household_id 70 0.03
## Antibtiotic 7 0.21
## Visit 6 0.03
## --------------------------------plot_model(md_phenotype1, type = "est") + theme_bw()
plot_model(md_phenotype1, type = "eff", terms = c("location_id","species_id"),colors = c("#69b3a2", "#404080")) +
theme_bw() + theme(axis.text.x = element_text(angle = 90, hjust = 1,size = 10))
drop1(md_phenotype1,test="Chisq") # dropped the three-way interaction ## Single term deletions
##
## Model:
## resistance ~ location_id * Bacteria + species_id * Bacteria +
## location_id * species_id * Bacteria + (1 | household_id) +
## (1 | Visit) + (1 | Antibtiotic)
## npar AIC LRT Pr(Chi)
## <none> 3890.7
## location_id:Bacteria:species_id 1 3891.7 2.9792 0.08434 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##################################################################################################
md_phenotype2<-glmer(data = final_df_melt,resistance~location_id*Bacteria+species_id*Bacteria+location_id*species_id+
(1|household_id)+(1|Visit)+(1|Antibtiotic),
family = binomial,control = glmerControl(optimizer ="bobyqa"))
summary(md_phenotype2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: resistance ~ location_id * Bacteria + species_id * Bacteria +
## location_id * species_id + (1 | household_id) + (1 | Visit) +
## (1 | Antibtiotic)
## Data: final_df_melt
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 3891.7 3954.9 -1935.9 3871.7 4078
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7903 -0.5350 -0.3350 0.0487 8.2632
##
## Random effects:
## Groups Name Variance Std.Dev.
## household_id (Intercept) 0.1117 0.3342
## Antibtiotic (Intercept) 0.9119 0.9549
## Visit (Intercept) 0.1296 0.3600
## Number of obs: 4088, groups: household_id, 70; Antibtiotic, 7; Visit, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.89793 0.40278 -2.229 0.025793 *
## location_idRural -0.54169 0.15554 -3.483 0.000496 ***
## BacteriaKlebsiella -0.63616 0.14178 -4.487 7.23e-06 ***
## species_idPig 0.01094 0.11115 0.098 0.921565
## location_idRural:BacteriaKlebsiella 0.32839 0.19123 1.717 0.085944 .
## BacteriaKlebsiella:species_idPig -0.48577 0.18823 -2.581 0.009862 **
## location_idRural:species_idPig -0.40760 0.17619 -2.313 0.020696 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lctn_R BctrKl spcs_P l_R:BK BcK:_P
## locatn_dRrl -0.150
## BactrKlbsll -0.106 0.187
## species_dPg -0.144 0.304 0.309
## lctn_dRr:BK 0.061 -0.422 -0.546 -0.081
## BctrKlbs:_P 0.056 0.040 -0.513 -0.384 -0.015
## lctn_dRr:_P 0.071 -0.519 -0.003 -0.488 0.089 -0.124drop1(md_phenotype2,test = "Chisq")# dropped bacteria*species interaction .## Single term deletions
##
## Model:
## resistance ~ location_id * Bacteria + species_id * Bacteria +
## location_id * species_id + (1 | household_id) + (1 | Visit) +
## (1 | Antibtiotic)
## npar AIC LRT Pr(Chi)
## <none> 3891.7
## location_id:Bacteria 1 3892.6 2.9193 0.087527 .
## Bacteria:species_id 1 3896.4 6.6982 0.009651 **
## location_id:species_id 1 3895.1 5.3311 0.020948 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1####################################################################################################
md_phenotype3<-glmer(data = final_df_melt,resistance~location_id*Bacteria+species_id+location_id*species_id+
(1|household_id)+(1|Visit)+(1|Antibtiotic),
family = binomial,control = glmerControl(optimizer ="bobyqa"))
summary(md_phenotype3)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: resistance ~ location_id * Bacteria + species_id + location_id *
## species_id + (1 | household_id) + (1 | Visit) + (1 | Antibtiotic)
## Data: final_df_melt
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 3896.4 3953.3 -1939.2 3878.4 4079
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8246 -0.5383 -0.3340 0.0372 7.5155
##
## Random effects:
## Groups Name Variance Std.Dev.
## household_id (Intercept) 0.1065 0.3263
## Antibtiotic (Intercept) 0.9091 0.9535
## Visit (Intercept) 0.1297 0.3601
## Number of obs: 4088, groups: household_id, 70; Antibtiotic, 7; Visit, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8395 0.4014 -2.091 0.036484 *
## location_idRural -0.5267 0.1532 -3.438 0.000585 ***
## BacteriaKlebsiella -0.8314 0.1208 -6.885 5.77e-12 ***
## species_idPig -0.1015 0.1020 -0.995 0.319782
## location_idRural:BacteriaKlebsiella 0.3236 0.1900 1.704 0.088474 .
## location_idRural:species_idPig -0.4655 0.1740 -2.675 0.007474 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lctn_R BctrKl spcs_P l_R:BK
## locatn_dRrl -0.151
## BactrKlbsll -0.085 0.234
## species_dPg -0.130 0.345 0.105
## lctn_dRr:BK 0.059 -0.407 -0.644 -0.071
## lctn_dRr:_P 0.078 -0.509 -0.056 -0.584 0.042drop1(md_phenotype3,test = "Chisq") # dropped location* bacteria interaction ## Single term deletions
##
## Model:
## resistance ~ location_id * Bacteria + species_id + location_id *
## species_id + (1 | household_id) + (1 | Visit) + (1 | Antibtiotic)
## npar AIC LRT Pr(Chi)
## <none> 3896.4
## location_id:Bacteria 1 3897.3 2.8731 0.090070 .
## location_id:species_id 1 3901.6 7.1462 0.007512 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#############################################################################################
md_phenotype4<-glmer(data = final_df_melt,resistance~relevel(location_id, ref = "Peri_urban")+Bacteria+species_id+relevel(location_id, ref = "Peri_urban")*species_id+
(1|household_id)+(1|Visit)+(1|Antibtiotic),
family = binomial,control = glmerControl(optimizer ="bobyqa"))
summary(md_phenotype4)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: resistance ~ relevel(location_id, ref = "Peri_urban") + Bacteria +
## species_id + relevel(location_id, ref = "Peri_urban") * species_id +
## (1 | household_id) + (1 | Visit) + (1 | Antibtiotic)
## Data: final_df_melt
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 3897.3 3947.8 -1940.6 3881.3 4080
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8046 -0.5344 -0.3368 0.0444 7.9593
##
## Random effects:
## Groups Name Variance Std.Dev.
## household_id (Intercept) 0.1049 0.3239
## Antibtiotic (Intercept) 0.9075 0.9526
## Visit (Intercept) 0.1381 0.3717
## Number of obs: 4088, groups: household_id, 70; Antibtiotic, 7; Visit, 6
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) -0.88020 0.40204
## relevel(location_id, ref = "Peri_urban")Rural -0.42146 0.13955
## BacteriaKlebsiella -0.70129 0.09220
## species_idPig -0.08928 0.10152
## relevel(location_id, ref = "Peri_urban")Rural:species_idPig -0.47844 0.17400
## z value Pr(>|z|)
## (Intercept) -2.189 0.02857
## relevel(location_id, ref = "Peri_urban")Rural -3.020 0.00253
## BacteriaKlebsiella -7.606 2.82e-14
## species_idPig -0.879 0.37915
## relevel(location_id, ref = "Peri_urban")Rural:species_idPig -2.750 0.00596
##
## (Intercept) *
## relevel(location_id, ref = "Peri_urban")Rural **
## BacteriaKlebsiella ***
## species_idPig
## relevel(location_id, ref = "Peri_urban")Rural:species_idPig **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) rl(_,r="P_")R BctrKl spcs_P
## rl(_,r="P_")R -0.139
## BactrKlbsll -0.063 -0.024
## species_dPg -0.126 0.348 0.078
## r(_,r="P_")R: 0.075 -0.542 -0.035 -0.581drop1(md_phenotype4,test = "Chisq") # final model the final model(Table 2) examining if the trends and factors seen in figure 3 are statistically significant## Single term deletions
##
## Model:
## resistance ~ relevel(location_id, ref = "Peri_urban") + Bacteria +
## species_id + relevel(location_id, ref = "Peri_urban") * species_id +
## (1 | household_id) + (1 | Visit) + (1 | Antibtiotic)
## npar AIC LRT
## <none> 3897.3
## Bacteria 1 3955.3 59.999
## relevel(location_id, ref = "Peri_urban"):species_id 1 3902.8 7.552
## Pr(Chi)
## <none>
## Bacteria 9.489e-15 ***
## relevel(location_id, ref = "Peri_urban"):species_id 0.005994 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(md_phenotype4,md_phenotype3)## Data: final_df_melt
## Models:
## md_phenotype4: resistance ~ relevel(location_id, ref = "Peri_urban") + Bacteria + species_id + relevel(location_id, ref = "Peri_urban") * species_id + (1 | household_id) + (1 | Visit) + (1 | Antibtiotic)
## md_phenotype3: resistance ~ location_id * Bacteria + species_id + location_id * species_id + (1 | household_id) + (1 | Visit) + (1 | Antibtiotic)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## md_phenotype4 8 3897.3 3947.8 -1940.7 3881.3
## md_phenotype3 9 3896.4 3953.3 -1939.2 3878.4 2.8731 1 0.09007 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tab_model(md_phenotype4)| resistance | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 0.41 | 0.19 – 0.91 | 0.029 |
|
relevel(location id, ref = “Peri urban”)Rural |
0.66 | 0.50 – 0.86 | 0.003 |
| Bacteria [Klebsiella] | 0.50 | 0.41 – 0.59 | <0.001 |
| species id [Pig] | 0.91 | 0.75 – 1.12 | 0.379 |
|
relevel(location id, ref = “Peri urban”)Rural × species id [Pig] |
0.62 | 0.44 – 0.87 | 0.006 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 household_id | 0.10 | ||
| τ00 Antibtiotic | 0.91 | ||
| τ00 Visit | 0.14 | ||
| ICC | 0.26 | ||
| N household_id | 70 | ||
| N Visit | 6 | ||
| N Antibtiotic | 7 | ||
| Observations | 4088 | ||
| Marginal R2 / Conditional R2 | 0.057 / 0.301 | ||
#############################################################################################
## GENE INFORM INFORMATION
###***********************************
genotype_EC$MDR_count<-PHENO_SPREAD_EC[match(genotype_EC$sample_id,PHENO_SPREAD_EC$SAMP_ID),]$MDR_count
genotype_KB$MDR_count<-PHENO_SPREAD_KB[match(genotype_KB$sample_id,PHENO_SPREAD_KB$SAMP_ID),]$MDR_count
genotype_EC<-genotype_EC[!is.na(genotype_EC$MDR_count),]
genotype_KB<-genotype_KB[!is.na(genotype_KB$MDR_count),]
# genotype_EC<-genotype_EC[,c(1:16,17),]
# genotype_KB<-genotype_KB[,c(2:16,17),]
Geno_db<- rbind(genotype_EC,genotype_KB)
genotype_melt<-melt(Geno_db, id.vars = c(2:8,13:17), measure.vars = 9:12)
names(genotype_melt)[13]<-"gene"
names(genotype_melt)[14]<-"norm_copy_num"
genotype_melt$gene<- as.character(genotype_melt$gene)
### distribution of gene copy number normalised by 16s
#****************************************************
ggplot(genotype_melt,aes(log10(norm_copy_num),fill=setting)) + geom_density(alpha=0.7) +
facet_wrap(~gene) + theme_bw() + ylab("Frequency") + xlab("log10(Gene copy number)") +
scale_fill_manual(values=c("#000080","#FC4E07","#800000")) + theme(legend.position="bottom") +
guides(fill=guide_legend(title="Setting"))## Warning: Removed 100 rows containing non-finite outside the scale range
## (`stat_density()`).
genotype_melt$gene[genotype_melt$gene=="dfra_1_norm"]<-"dfrA1"
genotype_melt$gene[genotype_melt$gene=="erm_b_norm"]<-"ermB"
genotype_melt$gene[genotype_melt$gene=="teb_b_norm"]<-"tetB"
genotype_melt$gene[genotype_melt$gene=="tet_q_norm"]<-"tetQ"
genotype_melt$host[genotype_melt$host=="FARMER"]<-"Farmer"
genotype_melt$host[genotype_melt$host=="PIG"]<-"Pig"
ggplot(genotype_melt,aes(log10(norm_copy_num),fill=host)) + geom_density(alpha=0.7) +
facet_wrap(~gene) + theme_bw() + ylab("Frequency") + xlab("log10(Normalised copies)") +
scale_fill_manual(values=c("#69b3a2", "#404080")) + theme(legend.position="bottom") +
guides(fill=guide_legend(title="Host"))## Warning: Removed 100 rows containing non-finite outside the scale range
## (`stat_density()`).
# from anova md_phenotype1 is afar better model than md_phenotype
######################################################################################
# Modeling genotypic resistance
df_genotype<-genotype_melt
df_genotype$genl0g<- log10(df_genotype$norm_copy_num)
df_genotype<-df_genotype[df_genotype$genl0g!="Inf",]
df_genotype<-df_genotype[df_genotype$genl0g!="-Inf",]
df_genotype$setting<- as.factor(df_genotype$setting)
df_genotype$host<-as.factor(df_genotype$host)
df_genotype$time<-as.factor(df_genotype$time)
df_genotype$hh_id<- substring(df_genotype$sample_id, 1,8)
md_genotype<-lmer(data=df_genotype,genl0g~host+setting+host*setting+
(1|hh_id)+(1|time)+(1|gene))
summ(md_genotype)## MODEL INFO:
## Observations: 2120
## Dependent Variable: genl0g
## Type: Mixed effects linear regression
##
## MODEL FIT:
## AIC = 6916.69, BIC = 6961.96
## Pseudo-R² (fixed effects) = 0.02
## Pseudo-R² (total) = 0.78
##
## FIXED EFFECTS:
## -------------------------------------------------------------------
## Est. S.E. t val. d.f. p
## -------------------------- ------- ------ -------- --------- ------
## (Intercept) -2.56 1.08 -2.37 3.16 0.09
## hostPig -0.52 0.09 -5.78 1880.45 0.00
## settingUrban 0.61 0.11 5.54 114.94 0.00
## hostPig:settingUrban 0.07 0.11 0.65 2041.04 0.51
## -------------------------------------------------------------------
##
## p values calculated using Satterthwaite d.f.
##
## RANDOM EFFECTS:
## ------------------------------------
## Group Parameter Std. Dev.
## ---------- ------------- -----------
## hh_id (Intercept) 0.32
## time (Intercept) 0.38
## gene (Intercept) 2.14
## Residual 1.20
## ------------------------------------
##
## Grouping variables:
## -------------------------
## Group # groups ICC
## ------- ---------- ------
## hh_id 70 0.02
## time 6 0.02
## gene 4 0.73
## -------------------------########################################################################################
df_genotype$gene<- as.character(df_genotype$gene)
df_genotype$gene[df_genotype$gene=="dfra_1_norm"]<-"dfra1"
df_genotype$gene[df_genotype$gene=="erm_b_norm"]<-"ermB"
df_genotype$gene[df_genotype$gene=="teb_b_norm"]<-"tetB"
df_genotype$gene[df_genotype$gene=="tet_q_norm"]<-"tetQ"
df_genotype<- df_genotype[!is.na(df_genotype$gene),]
md_genotype1<-lmer(data=df_genotype,genl0g~host*gene+setting*gene+host*setting*gene+
(1|hh_id)+(1|time))
summary(md_genotype1)## Linear mixed model fit by REML ['lmerMod']
## Formula: genl0g ~ host * gene + setting * gene + host * setting * gene +
## (1 | hh_id) + (1 | time)
## Data: df_genotype
##
## REML criterion at convergence: 6814.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5139 -0.6398 -0.0593 0.6147 3.6572
##
## Random effects:
## Groups Name Variance Std.Dev.
## hh_id (Intercept) 0.1027 0.3205
## time (Intercept) 0.1470 0.3834
## Residual 1.3751 1.1727
## Number of obs: 2120, groups: hh_id, 70; time, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -3.74998 0.20390 -18.391
## hostPig -0.02712 0.16706 -0.162
## geneermB 0.21267 0.16554 1.285
## genetetB -0.05509 0.16356 -0.337
## genetetQ 4.54976 0.16279 27.949
## settingUrban 0.76337 0.16994 4.492
## hostPig:geneermB -0.42298 0.23054 -1.835
## hostPig:genetetB -1.03482 0.23418 -4.419
## hostPig:genetetQ -0.49107 0.22832 -2.151
## geneermB:settingUrban 0.25396 0.21065 1.206
## genetetB:settingUrban -0.24540 0.20958 -1.171
## genetetQ:settingUrban -0.58724 0.20837 -2.818
## hostPig:settingUrban -0.30851 0.21345 -1.445
## hostPig:geneermB:settingUrban 0.62578 0.29541 2.118
## hostPig:genetetB:settingUrban 0.29503 0.30027 0.983
## hostPig:genetetQ:settingUrban 0.54412 0.29334 1.855##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itdrop1(md_genotype1,test="Chisq") # drop interaction host_id1:gene:setting_id1## Single term deletions
##
## Model:
## genl0g ~ host * gene + setting * gene + host * setting * gene +
## (1 | hh_id) + (1 | time)
## npar AIC LRT Pr(Chi)
## <none> 6813.9
## host:gene:setting 3 6813.4 5.4921 0.1391df_genotype$setting<- factor(df_genotype$setting)
df_genotype$gene<- factor(df_genotype$gene)
df_genotype$host<- factor(df_genotype$host)
md_genotype2<-lmer(data=df_genotype,genl0g~ relevel(gene, ref = "dfrA1") + relevel(setting, ref = "Urban") + relevel(host, ref = "Farmer") +
(1|hh_id)+(1|time))
md_genotype3<-lmer(data=df_genotype,genl0g~ relevel(gene, ref = "dfrA1") + relevel(setting, ref = "Urban")*relevel(host, ref = "Farmer") +
(1|hh_id)+(1|time))
md_genotype4<-lmer(data=df_genotype,genl0g~ gene + setting*host +
(1|hh_id)+(1|time))
tab_model(md_genotype3)| genl 0 g | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | -2.92 | -3.26 – -2.57 | <0.001 |
|
relevel(gene, ref = “dfrA1”)ermB |
0.35 | 0.20 – 0.49 | <0.001 |
|
relevel(gene, ref = “dfrA1”)tetB |
-0.62 | -0.76 – -0.47 | <0.001 |
|
relevel(gene, ref = “dfrA1”)tetQ |
4.11 | 3.97 – 4.25 | <0.001 |
|
relevel(setting, ref = “Urban”)Rural |
-0.61 | -0.82 – -0.39 | <0.001 |
|
relevel(host, ref = “Farmer”)Pig |
-0.44 | -0.58 – -0.31 | <0.001 |
|
relevel(setting, ref = “Urban”)Rural × relevel(host, ref = “Farmer”)Pig |
-0.07 | -0.29 – 0.15 | 0.513 |
| Random Effects | |||
| σ2 | 1.43 | ||
| τ00 hh_id | 0.10 | ||
| τ00 time | 0.15 | ||
| ICC | 0.15 | ||
| N hh_id | 70 | ||
| N time | 6 | ||
| Observations | 2120 | ||
| Marginal R2 / Conditional R2 | 0.685 / 0.731 | ||
tab_model(md_genotype4)| genl 0 g | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | -3.52 | -3.88 – -3.17 | <0.001 |
| gene [ermB] | 0.35 | 0.20 – 0.49 | <0.001 |
| gene [tetB] | -0.62 | -0.76 – -0.47 | <0.001 |
| gene [tetQ] | 4.11 | 3.97 – 4.25 | <0.001 |
| setting [Urban] | 0.61 | 0.39 – 0.82 | <0.001 |
| host [Pig] | -0.52 | -0.69 – -0.34 | <0.001 |
|
setting [Urban] × host [Pig] |
0.07 | -0.15 – 0.29 | 0.513 |
| Random Effects | |||
| σ2 | 1.43 | ||
| τ00 hh_id | 0.10 | ||
| τ00 time | 0.15 | ||
| ICC | 0.15 | ||
| N hh_id | 70 | ||
| N time | 6 | ||
| Observations | 2120 | ||
| Marginal R2 / Conditional R2 | 0.685 / 0.731 | ||
plot_model(md_genotype4, type = "est") + theme_bw()
plot_summs(md_genotype1,show.values = TRUE, scale = TRUE, plot.distributions = TRUE)## Loading required namespace: broom.mixed
plot_model(md_genotype1,show.values = TRUE, type = "est") + theme_bw()
#Diagnostics
###################################################################################
#Creating TET_B_Norm only object
df_TETB<-df_genotype%>%filter(gene=="tetB")
md_TETB<-lmer(data=df_TETB,genl0g~host+setting+host*setting+
(1|hh_id)+(1|time))
summary(md_TETB)## Linear mixed model fit by REML ['lmerMod']
## Formula: genl0g ~ host + setting + host * setting + (1 | hh_id) + (1 |
## time)
## Data: df_TETB
##
## REML criterion at convergence: 1726.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.83857 -0.66968 0.02178 0.70233 2.87527
##
## Random effects:
## Groups Name Variance Std.Dev.
## hh_id (Intercept) 0.09388 0.3064
## time (Intercept) 0.23465 0.4844
## Residual 1.62372 1.2743
## Number of obs: 507, groups: hh_id, 70; time, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -3.85047 0.24147 -15.946
## hostPig -1.01541 0.18793 -5.403
## settingUrban 0.53220 0.17920 2.970
## hostPig:settingUrban -0.03738 0.23837 -0.157
##
## Correlation of Fixed Effects:
## (Intr) hostPg sttngU
## hostPig -0.351
## settingUrbn -0.443 0.473
## hstPg:sttnU 0.282 -0.786 -0.618summ(md_TETB)## MODEL INFO:
## Observations: 507
## Dependent Variable: genl0g
## Type: Mixed effects linear regression
##
## MODEL FIT:
## AIC = 1740.18, BIC = 1769.78
## Pseudo-R² (fixed effects) = 0.15
## Pseudo-R² (total) = 0.29
##
## FIXED EFFECTS:
## ------------------------------------------------------------------
## Est. S.E. t val. d.f. p
## -------------------------- ------- ------ -------- -------- ------
## (Intercept) -3.85 0.24 -15.95 9.46 0.00
## hostPig -1.02 0.19 -5.40 497.37 0.00
## settingUrban 0.53 0.18 2.97 158.76 0.00
## hostPig:settingUrban -0.04 0.24 -0.16 496.86 0.88
## ------------------------------------------------------------------
##
## p values calculated using Satterthwaite d.f.
##
## RANDOM EFFECTS:
## ------------------------------------
## Group Parameter Std. Dev.
## ---------- ------------- -----------
## hh_id (Intercept) 0.31
## time (Intercept) 0.48
## Residual 1.27
## ------------------------------------
##
## Grouping variables:
## -------------------------
## Group # groups ICC
## ------- ---------- ------
## hh_id 70 0.05
## time 6 0.12
## -------------------------plot_model(md_TETB, type = "eff", terms = c("setting","host"),colors = c("#69b3a2", "#404080")) +
theme_bw() + theme(axis.text.x = element_text(angle = 90, hjust = 1,size = 10))
#############################################################################
#Creating ERM_B_Norm only object
df_ERM_B<-df_genotype%>%filter(gene=="ermB")
md_ERM_B<-lmer(data = df_ERM_B,genl0g~host+setting+host*setting+
(1|hh_id)+(1|time))
summary(md_ERM_B)## Linear mixed model fit by REML ['lmerMod']
## Formula: genl0g ~ host + setting + host * setting + (1 | hh_id) + (1 |
## time)
## Data: df_ERM_B
##
## REML criterion at convergence: 1767.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5673 -0.6557 -0.0435 0.6105 3.5163
##
## Random effects:
## Groups Name Variance Std.Dev.
## hh_id (Intercept) 0.2100 0.4583
## time (Intercept) 0.1387 0.3725
## Residual 1.3576 1.1652
## Number of obs: 540, groups: hh_id, 70; time, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -3.4631 0.2106 -16.447
## hostPig -0.5123 0.1723 -2.974
## settingUrban 0.8952 0.1891 4.735
## hostPig:settingUrban 0.4413 0.2161 2.042
##
## Correlation of Fixed Effects:
## (Intr) hostPg sttngU
## hostPig -0.419
## settingUrbn -0.532 0.466
## hstPg:sttnU 0.338 -0.797 -0.582summ(md_ERM_B)## MODEL INFO:
## Observations: 540
## Dependent Variable: genl0g
## Type: Mixed effects linear regression
##
## MODEL FIT:
## AIC = 1781.18, BIC = 1811.23
## Pseudo-R² (fixed effects) = 0.16
## Pseudo-R² (total) = 0.33
##
## FIXED EFFECTS:
## ------------------------------------------------------------------
## Est. S.E. t val. d.f. p
## -------------------------- ------- ------ -------- -------- ------
## (Intercept) -3.46 0.21 -16.45 14.41 0.00
## hostPig -0.51 0.17 -2.97 530.95 0.00
## settingUrban 0.90 0.19 4.74 127.96 0.00
## hostPig:settingUrban 0.44 0.22 2.04 525.49 0.04
## ------------------------------------------------------------------
##
## p values calculated using Satterthwaite d.f.
##
## RANDOM EFFECTS:
## ------------------------------------
## Group Parameter Std. Dev.
## ---------- ------------- -----------
## hh_id (Intercept) 0.46
## time (Intercept) 0.37
## Residual 1.17
## ------------------------------------
##
## Grouping variables:
## -------------------------
## Group # groups ICC
## ------- ---------- ------
## hh_id 70 0.12
## time 6 0.08
## -------------------------plot_model(md_ERM_B, type = "eff", terms = c("setting","host"),colors = c("#69b3a2", "#404080")) +
theme_bw() + theme(axis.text.x = element_text(angle = 90, hjust = 1,size = 10))
############################################################################
#Creating DFRA_1_Norm only object
df_DFRA<-df_genotype%>%filter(gene=="dfrA1")
md_DFRA<-lmer(data = df_DFRA,genl0g~host+setting+host*setting+
(1|hh_id)+(1|time))
summ(md_DFRA)## MODEL INFO:
## Observations: 521
## Dependent Variable: genl0g
## Type: Mixed effects linear regression
##
## MODEL FIT:
## AIC = 1800.54, BIC = 1830.33
## Pseudo-R² (fixed effects) = 0.06
## Pseudo-R² (total) = 0.24
##
## FIXED EFFECTS:
## ------------------------------------------------------------------
## Est. S.E. t val. d.f. p
## -------------------------- ------- ------ -------- -------- ------
## (Intercept) -3.75 0.24 -15.58 10.75 0.00
## hostPig -0.08 0.19 -0.42 511.94 0.68
## settingUrban 0.78 0.20 4.01 153.79 0.00
## hostPig:settingUrban -0.27 0.24 -1.12 508.05 0.26
## ------------------------------------------------------------------
##
## p values calculated using Satterthwaite d.f.
##
## RANDOM EFFECTS:
## ------------------------------------
## Group Parameter Std. Dev.
## ---------- ------------- -----------
## hh_id (Intercept) 0.41
## time (Intercept) 0.46
## Residual 1.27
## ------------------------------------
##
## Grouping variables:
## -------------------------
## Group # groups ICC
## ------- ---------- ------
## hh_id 70 0.09
## time 6 0.11
## -------------------------plot_model(md_DFRA, type = "eff", terms = c("setting","host"),colors = c("#69b3a2", "#404080")) +
theme_bw() + theme(axis.text.x = element_text(angle = 90, hjust = 1,size = 10))
######################################################################################################################
#Creating TETQ_Norm only object
df_TETQ<-df_genotype%>%filter(gene=="tetQ")
md_TETQ<-lmer(data = df_TETQ,genl0g~host+setting+host*setting+
(1|hh_id)+(1|time))
summ(md_TETQ)## MODEL INFO:
## Observations: 552
## Dependent Variable: genl0g
## Type: Mixed effects linear regression
##
## MODEL FIT:
## AIC = 1343.48, BIC = 1373.68
## Pseudo-R² (fixed effects) = 0.06
## Pseudo-R² (total) = 0.46
##
## FIXED EFFECTS:
## ------------------------------------------------------------------
## Est. S.E. t val. d.f. p
## -------------------------- ------- ------ -------- -------- ------
## (Intercept) 0.80 0.27 3.01 5.94 0.02
## hostPig -0.51 0.11 -4.75 542.84 0.00
## settingUrban 0.20 0.11 1.87 154.79 0.06
## hostPig:settingUrban 0.21 0.14 1.53 537.71 0.13
## ------------------------------------------------------------------
##
## p values calculated using Satterthwaite d.f.
##
## RANDOM EFFECTS:
## ------------------------------------
## Group Parameter Std. Dev.
## ---------- ------------- -----------
## hh_id (Intercept) 0.21
## time (Intercept) 0.62
## Residual 0.76
## ------------------------------------
##
## Grouping variables:
## -------------------------
## Group # groups ICC
## ------- ---------- ------
## hh_id 70 0.04
## time 6 0.38
## -------------------------plot_model(md_TETQ, type = "eff", terms = c("setting","host"),colors = c("#69b3a2", "#404080")) +
theme_bw() + theme(axis.text.x = element_text(angle = 90, hjust = 1,size = 10))
######################################################################################################################
### compare gene models using jtools
plot_summs(md_TETQ,md_TETB,md_ERM_B,md_DFRA,plot.distributions = T, model.names = c("TETQ","TETB","ERMB","DFRA1"))
##**********************************************************************
tetq<-spread(df_TETQ,value =genl0g,key=gene)
ermb<-spread(df_ERM_B,value =genl0g,key=gene)
tetb<-spread(df_TETB,value =genl0g,key=gene)
Dfra1<-spread(df_DFRA,value =genl0g,key=gene)
# Adding the spread out gene column back to the phenotypic data set .
combined_df$TET_Q_Norm<-tetq[match(combined_df$id2,tetq$sample_id),]$tetQ
combined_df$TET_B_Norm<-tetb[match(combined_df$id2,tetb$sample_id),]$tetB
combined_df$ERM_B_Norm<-ermb[match(combined_df$id2,ermb$sample_id),]$ermB
combined_df$DFRA_1_Norm<-Dfra1[match(combined_df$id2,Dfra1$sample_id),]$dfrA1
combined_df$id3<- paste(combined_df$hh_time,combined_df$species_id,sep = "")
## pig combine
combined_df_pig$TET_Q_Norm<-tetq[match(combined_df_pig$id2,tetq$sample_id),]$tetQ
combined_df_pig$TET_B_Norm<-tetb[match(combined_df_pig$id2,tetb$sample_id),]$tetB
combined_df_pig$ERM_B_Norm<-ermb[match(combined_df_pig$id2,ermb$sample_id),]$ermB
combined_df_pig$DFRA_1_Norm<-Dfra1[match(combined_df_pig$id2,Dfra1$sample_id),]$dfrA1
combined_df_pig$id3<- paste(combined_df_pig$hh_time,combined_df_pig$species_id,sep = "")
final_df<-rbind(combined_df,combined_df_pig)# construct the final data object | Datas
final_df_melt2<-melt(final_df, id.vars = c(1:5,7,13,16:19,20:47), measure.vars = 8:14)
names(final_df_melt2)[40]<-"Antibtiotic"
names(final_df_melt2)[41]<-"resistance"
final_df_melt2$resistance<- as.character(final_df_melt2$resistance)
final_df_melt2$resistance[final_df_melt2$resistance=="0"]<-"S"
final_df_melt2$resistance[final_df_melt2$resistance=="1"]<-"R"
final_df_melt2$resistance<- factor(final_df_melt2$resistance, levels=c("S","R"))
names(final_df_melt2)[36]<-"tetQ"
names(final_df_melt2)[37]<-"tetB"
names(final_df_melt2)[38]<-"ermB"
names(final_df_melt2)[39]<-"dfra1"
names(final_df_melt2)[30]<-"breed_verif"
names(final_df_melt2)[3]<-"location_id"
final_df_melt2<- final_df_melt2[!is.na(final_df_melt2$tetQ),]
final_df_melt2<- final_df_melt2[!is.na(final_df_melt2$tetB),]
final_df_melt2<- final_df_melt2[!is.na(final_df_melt2$ermB),]
final_df_melt2<- final_df_melt2[!is.na(final_df_melt2$dfra1),]
final_df_melt2$VIZIT[final_df_melt2$Visit==1]<-"ONE" ## Categorical version of the temporal variable
final_df_melt2$VIZIT[final_df_melt2$Visit==2]<-"TWO"
final_df_melt2$VIZIT[final_df_melt2$Visit==3]<-"THREE"
final_df_melt2$VIZIT[final_df_melt2$Visit==4]<-"FOUR"
final_df_melt2$VIZIT[final_df_melt2$Visit==5]<-"FIVE"
final_df_melt2$VIZIT[final_df_melt2$Visit==6]<-"SIX"
final_df_melt2$species_id.x[final_df_melt2$species_id.x=="H"]<-"Farmer"
final_df_melt2$species_id.x[final_df_melt2$species_id.x=="P"]<-"Pig"
names(final_df_melt2)[5]<-"species_id"
## FINAL MODEL
final_df_melt2$location_id<- as.character(final_df_melt2$l)
final_df_melt2$location_id[final_df_melt2$location_id=="Rural"]<-"Free_range"
final_df_melt2$location_id[final_df_melt2$location_id=="Peri_urban"]<-"Semi_intensive"
final_df_melt2$location_id<- factor(final_df_melt2$location_id,levels =c("Semi_intensive","Free_range") )
final_df_melt2$breed_verif<-factor(final_df_melt2$breed_verif,levels = c("Exotic","Local","Mixed"))
names(final_df_melt2)[26]<-"household_id"
final_df_melt2$Antibtiotic<- as.character(final_df_melt2$Antibtiotic)
final_df_melt2$Antibtiotic[final_df_melt2$Antibtiotic=="Sulfa_Trimethoprim"]<-"TMPS"
final_df_melt2$Antibtiotic[final_df_melt2$Antibtiotic=="ciproflaxacine"]<-"Ciproflaxacine"
final_df_melt2$Antibtiotic[final_df_melt2$Antibtiotic=="gentamycine"]<-"Gentamycine"
final_df_melt2$Antibtiotic[final_df_melt2$Antibtiotic=="nalidixic_acid"]<-"Nalidixic_acid"
final_df_melt2$Antibtiotic[final_df_melt2$Antibtiotic=="strepptomycine"]<-"Strepptomycine"
final_df_melt2$Antibtiotic[final_df_melt2$Antibtiotic=="tetracycline"]<-"Tetracycline"
final_df_melt2$Antibtiotic<-factor(final_df_melt2$Antibtiotic,levels = c("Ciproflaxacine","Gentamycine","Nalidixic_acid","Chloroamphenicol","Strepptomycine","Tetracycline","TMPS" ))
final_df_melt2$breed_verif<-factor(final_df_melt2$breed_verif,levels = c("Exotic","Mixed","Local"))
final_df_melt2<-final_df_melt2[!final_df_melt2$Bacteria=="Salmonella",]
names(final_df_melt2)[30]<-"Breeds"
final_df_melt2$Breeds<- as.character(final_df_melt2$Breeds)
final_df_melt2$Breeds[final_df_melt2$Breeds=="LOCAL"]<-"Local"
final_df_melt2$Breeds[final_df_melt2$Breeds=="MIXED_BREED"]<-"Mixed"
final_df_melt2$Breeds[final_df_melt2$Breeds=="EXOTIC_BREED"]<-"Exotic"
names(final_df_melt2)[3]<-"Production"
final_df_melt2$Antibtiotic<- as.character(final_df_melt2$Antibtiotic)
final_df_melt2$Antibtiotic[final_df_melt2$Antibtiotic=="Sulfa_Trimethoprim"]<-"Trimethoprim/sulfamethoxazole"
names(final_df_melt2)[5]<-"Host"
phen_gen_mdl13<-glmer(data = final_df_melt2,resistance~Production+Bacteria+
Breeds+Antibtiotic+Host+Visit+tetQ+
Production*Host +
(1|Sample_ID)+(1|VIZIT),
family = binomial,control = glmerControl(optimizer ="bobyqa"))
summary(phen_gen_mdl13)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: resistance ~ Production + Bacteria + Breeds + Antibtiotic + Host +
## Visit + tetQ + Production * Host + (1 | Sample_ID) + (1 | VIZIT)
## Data: final_df_melt2
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 2753.3 2855.1 -1359.6 2719.3 2930
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8151 -0.5243 -0.3175 0.5167 8.0297
##
## Random effects:
## Groups Name Variance Std.Dev.
## Sample_ID (Intercept) 0.09673 0.3110
## VIZIT (Intercept) 0.03336 0.1827
## Number of obs: 2947, groups: Sample_ID, 65; VIZIT, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.07278 0.34349 -6.034 1.60e-09 ***
## ProductionFree_range -0.35402 0.15698 -2.255 0.02412 *
## BacteriaKlebsiella -0.64491 0.11074 -5.824 5.75e-09 ***
## BreedsLocal 0.31696 0.17644 1.796 0.07244 .
## BreedsMixed 0.28404 0.15342 1.851 0.06411 .
## AntibtioticCiproflaxacine -0.57243 0.23693 -2.416 0.01569 *
## AntibtioticGentamycine -0.35372 0.22562 -1.568 0.11694
## AntibtioticNalidixic_acid 0.30410 0.20169 1.508 0.13162
## AntibtioticStrepptomycine 1.34150 0.18509 7.248 4.23e-13 ***
## AntibtioticTetracycline 1.91948 0.18324 10.475 < 2e-16 ***
## AntibtioticTMPS 2.04775 0.18340 11.166 < 2e-16 ***
## HostPig -0.09706 0.11779 -0.824 0.40994
## Visit 0.17082 0.05584 3.059 0.00222 **
## tetQ 0.11956 0.04249 2.814 0.00489 **
## ProductionFree_range:HostPig -0.71870 0.23105 -3.111 0.00187 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itsumm(phen_gen_mdl13)## MODEL INFO:
## Observations: 2947
## Dependent Variable: resistance
## Type: Mixed effects generalized linear regression
## Error Distribution: binomial
## Link function: logit
##
## MODEL FIT:
## AIC = 2753.29, BIC = 2855.09
## Pseudo-R² (fixed effects) = 0.29
## Pseudo-R² (total) = 0.32
##
## FIXED EFFECTS:
## -----------------------------------------------------------------
## Est. S.E. z val. p
## ---------------------------------- ------- ------ -------- ------
## (Intercept) -2.07 0.34 -6.03 0.00
## ProductionFree_range -0.35 0.16 -2.26 0.02
## BacteriaKlebsiella -0.64 0.11 -5.82 0.00
## BreedsLocal 0.32 0.18 1.80 0.07
## BreedsMixed 0.28 0.15 1.85 0.06
## AntibtioticCiproflaxacine -0.57 0.24 -2.42 0.02
## AntibtioticGentamycine -0.35 0.23 -1.57 0.12
## AntibtioticNalidixic_acid 0.30 0.20 1.51 0.13
## AntibtioticStrepptomycine 1.34 0.19 7.25 0.00
## AntibtioticTetracycline 1.92 0.18 10.48 0.00
## AntibtioticTMPS 2.05 0.18 11.17 0.00
## HostPig -0.10 0.12 -0.82 0.41
## Visit 0.17 0.06 3.06 0.00
## tetQ 0.12 0.04 2.81 0.00
## ProductionFree_range:HostPig -0.72 0.23 -3.11 0.00
## -----------------------------------------------------------------
##
## RANDOM EFFECTS:
## -------------------------------------
## Group Parameter Std. Dev.
## ----------- ------------- -----------
## Sample_ID (Intercept) 0.31
## VIZIT (Intercept) 0.18
## -------------------------------------
##
## Grouping variables:
## -----------------------------
## Group # groups ICC
## ----------- ---------- ------
## Sample_ID 65 0.03
## VIZIT 6 0.01
## -----------------------------tab_model(phen_gen_mdl13)| resistance | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 0.13 | 0.06 – 0.25 | <0.001 |
| Production [Free_range] | 0.70 | 0.52 – 0.95 | 0.024 |
| Bacteria [Klebsiella] | 0.52 | 0.42 – 0.65 | <0.001 |
| Breeds [Local] | 1.37 | 0.97 – 1.94 | 0.072 |
| Breeds [Mixed] | 1.33 | 0.98 – 1.79 | 0.064 |
|
Antibtiotic [Ciproflaxacine] |
0.56 | 0.35 – 0.90 | 0.016 |
| Antibtiotic [Gentamycine] | 0.70 | 0.45 – 1.09 | 0.117 |
|
Antibtiotic [Nalidixic_acid] |
1.36 | 0.91 – 2.01 | 0.132 |
|
Antibtiotic [Strepptomycine] |
3.82 | 2.66 – 5.50 | <0.001 |
|
Antibtiotic [Tetracycline] |
6.82 | 4.76 – 9.76 | <0.001 |
| Antibtiotic [TMPS] | 7.75 | 5.41 – 11.10 | <0.001 |
| Host [Pig] | 0.91 | 0.72 – 1.14 | 0.410 |
| Visit | 1.19 | 1.06 – 1.32 | 0.002 |
| tetQ | 1.13 | 1.04 – 1.22 | 0.005 |
|
Production [Free_range] × Host [Pig] |
0.49 | 0.31 – 0.77 | 0.002 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 Sample_ID | 0.10 | ||
| τ00 VIZIT | 0.03 | ||
| ICC | 0.04 | ||
| N Sample_ID | 65 | ||
| N VIZIT | 6 | ||
| Observations | 2947 | ||
| Marginal R2 / Conditional R2 | 0.295 / 0.322 | ||
plot_model(phen_gen_mdl13,show.values = TRUE, type = "est") + theme_bw()
plot_model(phen_gen_mdl13, type = "eff", terms = c("Antibtiotic","Production"),colors =c("#FC4E07","#000080","#800000")) +
theme_bw() + theme(axis.text.x = element_text(angle = 90, hjust = 1,size = 10)) + ylab("Resistance") + xlab("Antibiotics")
## Figures for paper
plot_model(phen_gen_mdl13, type = "eff",title = "", terms = c("Host","Production"),colors = c("#FC4E07","#000080","#800000")) +
theme_bw() + ylab("Resistance") + xlab("Hosts") + ylab("Adjusted ABR prevalence")
plot_model(phen_gen_mdl13, type = "eff", title = "", terms = c("Bacteria","Production"),colors = c("#FC4E07","#000080","#800000")) +
theme_bw() + ylab("Adjusted ABR prevalence")
plot_model(phen_gen_mdl13, type = "eff", title = "", terms = c("Breeds","Production"),colors = c("#FC4E07","#000080","#800000")) +
theme_bw() +xlab("Hosts") + ylab("Adjusted ABR prevalence") + xlab("Pig breeds")
plot_model(phen_gen_mdl13, type = "eff", terms = c("Visit","Production"),colors = c("#FC4E07","#000080","#800000")) +
theme_bw() + xlab("Visit every two months")
plot_model(phen_gen_mdl13, type = "eff", terms = c("tetQ [all]","Production"),colors = c("#FC4E07","#000080","#800000")) +
theme_bw() + xlab("tetQ Normalise copy number")
plot_model(phen_gen_mdl13, type = "eff", terms = c("tetQ [all]","Host"),colors = c("#69b3a2", "#404080")) +
theme_bw() 
plot_model(phen_gen_mdl13, type = "eff", terms = c("tetQ [all]","Host","Production"),colors = c("#69b3a2", "#404080")) +
theme_bw() 
plot_model(phen_gen_mdl13, type = "eff", title = "", terms = c("Visit [all]","Production"),colors = c("#FC4E07","#000080","#800000")) +
theme_bw() + xlab("Visits") + ylab("Adjusted ABR prevalence")
plot_model(phen_gen_mdl13, type = "eff", terms = c("tetQ [all]","Host"),colors = c("#69b3a2", "#404080")) +
theme_bw() + xlab("Log10(Copy number of tetQ)") + ylab("Antibiotic resistance")
plot_model(phen_gen_mdl13, type = "eff", title = "", terms = c("tetQ [all]","Production"),colors = c("#FC4E07","#000080","#800000")) +
theme_bw() + xlab("Log10(Normalized tetQ copies)") + ylab("Antibiotic resistance")
plot_model(phen_gen_mdl13, type = "eff", terms = c("Visit [all]","Host"),colors = c("#FC4E07","#000080","#800000")) +
theme_bw() + xlab("Time in months") + ylab("Antibiotic resistance")
## DOES THE CARRIAGE OF A FARMER AT time T AFFECT THAT OF thier pig at
time T+1?
##############################################################################
staggered$GENE_ID1[staggered$GENE_ID1==Inf]<-0
staggered$GENE_ID2[staggered$GENE_ID2==Inf]<-0
cor.test(staggered$GENE_ID1,staggered$GENE_ID2)##
## Pearson's product-moment correlation
##
## data: staggered$GENE_ID1 and staggered$GENE_ID2
## t = -0.10281, df = 748, p-value = 0.9181
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.07532733 0.06784797
## sample estimates:
## cor
## -0.003758944# Ensure all names of columns in the data frame are in lower case
phenotype_data<- PHENODATA1
phenotype_data$host[phenotype_data$species_id=="Human"]<- 'H' # H for human
phenotype_data$host[phenotype_data$species_id=="Pig"]<-'P' # P for pig host .
phenotype_data$id1<-paste(phenotype_data$household_id,phenotype_data$Visit,phenotype_data$host,sep="" ) # add id1 column
# to be used for merging with staggered data .
names(staggered)<-tolower(names(staggered))
Combined_df<-phenotype_data%>% inner_join(staggered,by="id1") # joining the two## Warning in inner_join(., staggered, by = "id1"): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 2 of `x` matches multiple rows in `y`.
## ℹ Row 55 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.Combined_df$Phen_id2<- phenotype_data[match(Combined_df$id2,phenotype_data$id1),]$resistance
Combined_df$Phen_id1<- phenotype_data[match(Combined_df$id1,phenotype_data$id1),]$resistance
combined_df3<-Combined_df[!is.na(Combined_df$Phen_id2),]
#data frames by column id1.
Combined_df$household_id<- substr(Combined_df$sampleid, 1,8)
Combined_df<-mutate(Combined_df,loggene_id1=log(gene_id1))# add column
#log.gene_id1
Combined_df<-filter(Combined_df,loggene_id1!="-Inf") # remove -inf
Combined_df$time_id1<-as.factor(Combined_df$time_id1)
Combined_df$host_id1<-as.factor(Combined_df$host_id1)
Combined_df$loggene_id2<- log10(Combined_df$gene_id2)
Combined_df$loggene_id1<- log10(Combined_df$gene_id1)
Combined_df<-filter(Combined_df,loggene_id1!="-Inf") # remove -inf
Combined_df<-filter(Combined_df,loggene_id2!="-Inf") # remove -inf
combined_md<-glm(data = Combined_df, loggene_id1~ loggene_id2 + location_id+ gene)
## Farmer at T and pig at T2 (Final)
Combined_df$location_id[Combined_df$location_id=="Peri_urban"]<-"Semi-intensive"
Combined_df$location_id[Combined_df$location_id=="Rural"]<-"Free-Range"
Combined_df$gene[Combined_df$gene=="DFRA_1_Norm"]<-"dfra1"
Combined_df$gene[Combined_df$gene=="ERM_B_Norm"]<-"ermB"
Combined_df$gene[Combined_df$gene=="TET_B_Norm"]<-"tetB"
Combined_df$gene[Combined_df$gene=="TET_Q_Norm"]<-"tetQ"
names(Combined_df)[3]<-"Production"
combined_md2<-lmer(loggene_id1 ~ loggene_id2 + Production + gene + Production*gene + (1| household_id), Combined_df)
#levels(combined_df$time_id1)
summary(combined_md)##
## Call:
## glm(formula = loggene_id1 ~ loggene_id2 + location_id + gene,
## data = Combined_df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.6277 -0.8478 -0.0385 0.7751 4.6040
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.80279 0.20946 -13.381 < 2e-16 ***
## loggene_id2 0.03383 0.04938 0.685 0.493643
## location_idRural -0.45621 0.12335 -3.698 0.000239 ***
## geneERM_B_Norm 0.17744 0.16316 1.088 0.277294
## geneTET_B_Norm -0.32889 0.17584 -1.870 0.061957 .
## geneTET_Q_Norm 3.83836 0.26582 14.440 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 1.61659)
##
## Null deviance: 2772.80 on 551 degrees of freedom
## Residual deviance: 882.66 on 546 degrees of freedom
## AIC: 1839.6
##
## Number of Fisher Scoring iterations: 2summ(combined_md2)## MODEL INFO:
## Observations: 552
## Dependent Variable: loggene_id1
## Type: Mixed effects linear regression
##
## MODEL FIT:
## AIC = 1812.86, BIC = 1860.31
## Pseudo-R² (fixed effects) = 0.68
## Pseudo-R² (total) = 0.74
##
## FIXED EFFECTS:
## -------------------------------------------------------------------------------
## Est. S.E. t val. d.f. p
## --------------------------------------- ------- ------ -------- -------- ------
## (Intercept) -3.45 0.37 -9.25 463.02 0.00
## loggene_id2 -0.02 0.05 -0.36 542.63 0.72
## ProductionSemi-intensive 0.46 0.35 1.29 317.70 0.20
## geneermB -0.16 0.33 -0.49 531.63 0.62
## genetetB -0.33 0.34 -1.00 524.27 0.32
## genetetQ 4.34 0.40 10.93 538.01 0.00
## ProductionSemi-intensive:geneermB 0.55 0.37 1.47 524.74 0.14
## ProductionSemi-intensive:genetetB -0.07 0.38 -0.18 518.35 0.85
## ProductionSemi-intensive:genetetQ -0.43 0.37 -1.18 523.32 0.24
## -------------------------------------------------------------------------------
##
## p values calculated using Satterthwaite d.f.
##
## RANDOM EFFECTS:
## ----------------------------------------
## Group Parameter Std. Dev.
## -------------- ------------- -----------
## household_id (Intercept) 0.59
## Residual 1.14
## ----------------------------------------
##
## Grouping variables:
## --------------------------------
## Group # groups ICC
## -------------- ---------- ------
## household_id 58 0.21
## --------------------------------tab_model(combined_md2)| loggene id 1 | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | -3.45 | -4.18 – -2.72 | <0.001 |
| loggene id2 | -0.02 | -0.12 – 0.08 | 0.721 |
|
Production [Semi-intensive] |
0.46 | -0.24 – 1.15 | 0.197 |
| gene [ermB] | -0.16 | -0.81 – 0.49 | 0.623 |
| gene [tetB] | -0.33 | -0.99 – 0.32 | 0.318 |
| gene [tetQ] | 4.34 | 3.56 – 5.12 | <0.001 |
|
Production [Semi-intensive] × gene [ermB] |
0.55 | -0.18 – 1.28 | 0.141 |
|
Production [Semi-intensive] × gene [tetB] |
-0.07 | -0.81 – 0.67 | 0.855 |
|
Production [Semi-intensive] × gene [tetQ] |
-0.43 | -1.15 – 0.29 | 0.238 |
| Random Effects | |||
| σ2 | 1.31 | ||
| τ00 household_id | 0.35 | ||
| ICC | 0.21 | ||
| N household_id | 58 | ||
| Observations | 552 | ||
| Marginal R2 / Conditional R2 | 0.675 / 0.744 | ||
plot_model(combined_md2, type = "eff", terms = c("gene [all]","Production"),colors = c("#FC4E07","#000080","#800000")) + theme_bw() + xlab("Time in months") + ylab("Log10(gene copy number)")
plot_model(combined_md2, type = "eff", terms = c("loggene_id2 [all]","Production","gene"),colors = c("#FC4E07","#000080","#800000")) + theme_bw() + xlab("Log10(gene copy number) for Pig at time T+1") + ylab("Log10(gene copy number) for farmer at time T")
## GENERATE TRANSMISSION DATA
Pattern_data$pattern<-paste(Pattern_data$ciproflaxacine,Pattern_data$gentamycine,Pattern_data$nalidixic_acid,
Pattern_data$strepptomycine,Pattern_data$sulfa_trimethoprim,Pattern_data$tetracycline,sep = "")
Pattern_data$patternSUM<-as.numeric(Pattern_data$ciproflaxacine+Pattern_data$gentamycine+Pattern_data$nalidixic_acid +
Pattern_data$strepptomycine+Pattern_data$sulfa_trimethoprim+Pattern_data$tetracycline)## This creates a column with the MDR patterns
HH_data<-unique(HH_data)
HH_data$Pattern_farmer<- Pattern_data[match(HH_data$Farmer_ID,Pattern_data$ID3),]$pattern
HH_data$Pattern_pig<- Pattern_data[match(HH_data$Pig_ID,Pattern_data$ID3),]$pattern
Pattern_data$id3<- paste(Pattern_data$household_id,Pattern_data$SAMPLING.POINT,Pattern_data$Hostcode,sep = "")
Pattern_data<-Pattern_data[Pattern_data$bacteria=="E.coli",] ## We are only using E.coli here
Pattern_data_MDR<-Pattern_data[Pattern_data$patternSUM>2,]
## PLOT FOR FREQUENCY OF UNIQUE PATTERSN
###***************************************
Summarised_pattern<-Pattern_data_MDR %>%
dplyr::count(pattern,location_id,bacteria,species_id) %>%
mutate(Proportion = (n/110)*100)
## PLOT FOR FREQUENCY BY LOCATIONS
ggplot(Summarised_pattern,aes(x= fct_reorder(pattern,Proportion,.desc = F),Proportion,fill=location_id)) + geom_bar(stat="identity") +
guides(fill=guide_legend(title="Setting")) + coord_flip() + theme_bw()+ scale_fill_manual(values=c("#000080","#FC4E07","#800000")) + xlab("MDR patterns")
# As you can see the prevalence of these MDR is higher in urban
# To create the conservative dataset, we remove the most abundant pattern, so that the prevalens is less or equal to 10%(see line 1277-1279)
## PLOT FOR FREQUENCY BY HOST
ggplot(Summarised_pattern,aes(x= fct_reorder(pattern,Proportion,.desc = F),Proportion,fill=species_id)) + geom_bar(stat="identity") +
guides(fill=guide_legend(title="Host")) + coord_flip() + theme_bw() +
scale_fill_manual(values=c("#69b3a2", "#404080")) + xlab("MDR patterns")
# But pigs and humans seem to carry then in comparable propotions, in general farmers carry more
#**************************************************************************************************
HH_data$Pattern_farmer<- Pattern_data_MDR[match(HH_data$Farmer_ID,Pattern_data_MDR$ID3),]$pattern
HH_data$Pattern_pig<- Pattern_data_MDR[match(HH_data$Pig_ID,Pattern_data_MDR$ID3),]$pattern
HH_data$Transmission[HH_data$Pattern_farmer==HH_data$Pattern_pig]<-"Yes"
HH_data$Transmission[is.na(HH_data$Transmission)]<-"NO"
NAs<-HH_data[is.na(HH_data$Pattern_farmer) & is.na(HH_data$Pattern_pig),]
HH_data<- HH_data[!HH_data$Farmer_ID %in%NAs$Farmer_ID,]
### Create a matrix with all
MAT_HUM_PIG_same_HHx<-HH_data ## HUMAN AND PIGS IN THE SAME HOUSEHOLD
HH_data_P<-HH_data[str_detect(HH_data$Pig_ID, "P"),]
HH_data_PP<-HH_data[str_detect(HH_data$Farmer_ID, "P"),]
HH_data_PP<-HH_data_PP[str_detect(HH_data_PP$Farmer_ID, "H"),]
unique(HH_data_P$Pig_ID) ## double checking ## [1] "KLAKIR014P" "KLAKIR021P" "KLAKIR023P" "KLAKIR026P" "KLAKIR034P"
## [6] "KLAKIR036P" "KLAKIR051P" "KLAKIR055P" "KLAKIR073P" "KLAKIR074P"
## [11] "KLAKYN013P" "KLAKYN014P" "KLAKYN016P" "KLAKYN023P" "KLAKYN024P"
## [16] "KLAKYN034P" "KLAKYN044P" "KLAKYN045P" "KLAKYN055P" "KLAKYN056P"
## [21] "KLAKYN064P" "KLAKYN071P" "KLAKYN074P" "KLAMAK011P" "KLAMAK013P"
## [26] "KLAMAK016P" "KLAMAK023P" "KLAMAK025P" "KLAMAK035P" "KLAMAK036P"
## [31] "KLAMAK044P" "KLAMAK053P" "KLAMAK062P" "KLAMAK074P" "KLANAK023P"
## [36] "KLANAK024P" "KLANAK025P" "KLANAK034P" "KLANAK044P" "KLANAK054P"
## [41] "KLANAK075P" "KLANAK076P" "KLANAN012P" "KLANAN013P" "KLANAN014P"
## [46] "KLANAN015P" "KLANAN025P" "KLANAN026P" "KLANAN033P" "KLANAN044P"
## [51] "KLANAN045P" "KLANAN066P" "KLANAN072P" "MBDBAG011P" "MBDBAG013P"
## [56] "MBDBAG026P" "MBDBAG043P" "MBDBAG045P" "MBDBAG052P" "MBDKIG015P"
## [61] "MBDKIG043P" "MBDKIG045P" "MBDKIG076P" "MBDKIT012P" "MBDKIT014P"
## [66] "MBDKIT053P" "MBDKIT062P" "MBDKIT074P" "MBDKIY032P" "MBDKIY042P"
## [71] "MBDKIY072P" "MBDMAD062P" "MBDMAD071P" "KLAMAK063P" "KLANAN054P"unique(HH_data_P$Farmer_ID)## double checking ## [1] "KLAKIR014H" "KLAKIR021H" "KLAKIR023H" "KLAKIR026H" "KLAKIR034H"
## [6] "KLAKIR036H" "KLAKIR051H" "KLAKIR055H" "KLAKIR073H" "KLAKIR074H"
## [11] "KLAKYN013H" "KLAKYN014H" "KLAKYN016H" "KLAKYN023H" "KLAKYN024H"
## [16] "KLAKYN034H" "KLAKYN044H" "KLAKYN045H" "KLAKYN055H" "KLAKYN056H"
## [21] "KLAKYN064H" "KLAKYN071H" "KLAKYN074H" "KLAMAK011H" "KLAMAK013H"
## [26] "KLAMAK016H" "KLAMAK023H" "KLAMAK025H" "KLAMAK035H" "KLAMAK036H"
## [31] "KLAMAK044H" "KLAMAK053H" "KLAMAK062H" "KLAMAK074H" "KLANAK023H"
## [36] "KLANAK024H" "KLANAK025H" "KLANAK034H" "KLANAK044H" "KLANAK054H"
## [41] "KLANAK075H" "KLANAK076H" "KLANAN012H" "KLANAN013H" "KLANAN014H"
## [46] "KLANAN015H" "KLANAN025H" "KLANAN026H" "KLANAN033H" "KLANAN044H"
## [51] "KLANAN045H" "KLANAN066H" "KLANAN072H" "MBDBAG011H" "MBDBAG013H"
## [56] "MBDBAG026H" "MBDBAG043H" "MBDBAG045H" "MBDBAG052H" "MBDKIG015H"
## [61] "MBDKIG043H" "MBDKIG045H" "MBDKIG076H" "MBDKIT012H" "MBDKIT014H"
## [66] "MBDKIT053H" "MBDKIT062H" "MBDKIT074H" "MBDKIY032H" "MBDKIY042H"
## [71] "MBDKIY072H" "MBDMAD062H" "MBDMAD071H" "KLAMAK063H" "KLANAN054H"HH_data_F<-HH_data[str_detect(HH_data$Farmer_ID, "H"),]
HH_data_FF<-HH_data[str_detect(HH_data$Pig_ID, "H"),]
MAT_HUM_PIG<- matrix(data = NA, nrow = length(HH_data_MDR$Host), ncol = length(HH_data_MDR$Reference),
dimnames = list(HH_data_MDR$Host,HH_data_MDR$Reference))
MAT_HUM_PIG<-melt(as.matrix(MAT_HUM_PIG))
MAT_HUM_PIG$Var1<- as.character(MAT_HUM_PIG$Var1)
MAT_HUM_PIG$Var2<- as.character(MAT_HUM_PIG$Var2)
MAT_HUM_PIG<-MAT_HUM_PIG[!MAT_HUM_PIG$Var1==MAT_HUM_PIG$Var2,]
MAT_HUM_PIG<- MAT_HUM_PIG[,c(1:2),]
names(MAT_HUM_PIG)[1]<-"HOST1"
names(MAT_HUM_PIG)[2]<-"HOST2"
MAT_HUM_PIG$H1_pattern<-HH_data_MDR[match(MAT_HUM_PIG$HOST1,HH_data_MDR$Host),]$Host_Pattern
MAT_HUM_PIG$H2_pattern<-HH_data_MDR[match(MAT_HUM_PIG$HOST2,HH_data_MDR$Reference),]$Ref_pattern
MAT_HUM_PIG$SET1[str_detect(MAT_HUM_PIG$HOST1, "KLA")]<-"Urban"
MAT_HUM_PIG$SET2[str_detect(MAT_HUM_PIG$HOST2, "KLA")]<-"Urban"
MAT_HUM_PIG$SET1[str_detect(MAT_HUM_PIG$HOST1, "MBD")]<-"Rural"
MAT_HUM_PIG$SET2[str_detect(MAT_HUM_PIG$HOST2, "MBD")]<-"Rural"
MAT_HUM_PIG$same_dis <- MAT_HUM_PIG$SET1 == MAT_HUM_PIG$SET2
MAT_HUM_PIG$HOST1_SBC<-substr(MAT_HUM_PIG$HOST1,4,6)
MAT_HUM_PIG$HOST2_SBC<-substr(MAT_HUM_PIG$HOST2,4,6)
MAT_HUM_PIG$same_SBC <- MAT_HUM_PIG$HOST1_SBC == MAT_HUM_PIG$HOST2_SBC
MAT_HUM_PIG$HOST1_HH<-substr(MAT_HUM_PIG$HOST1,1,8)
MAT_HUM_PIG$HOST2_HH<-substr(MAT_HUM_PIG$HOST2,1,8)
MAT_HUM_PIG$same_HH <- MAT_HUM_PIG$HOST1_HH == MAT_HUM_PIG$HOST2_HH
MAT_HUM_PIG_same_HH<- MAT_HUM_PIG[MAT_HUM_PIG$same_HH==T,]
NAz<-MAT_HUM_PIG_same_HH[is.na(MAT_HUM_PIG_same_HH$H2_pattern) & is.na(MAT_HUM_PIG_same_HH$H1_pattern),]
MAT_HUM_PIG_same_HH<- MAT_HUM_PIG_same_HH[!MAT_HUM_PIG_same_HH$HOST1 %in%NAz$HOST1,]
MAT_HUM_PIG_same_HH$same <- MAT_HUM_PIG_same_HH$H1_pattern == MAT_HUM_PIG_same_HH$H2_pattern
MAT_HUM_PIG_same_HH<-MAT_HUM_PIG_same_HH[!is.na(MAT_HUM_PIG_same_HH$same),]
MAT_HUM_PIGNONHH<- MAT_HUM_PIG[MAT_HUM_PIG$same_HH!=T,]
MAT_HUM_PIG_SBC<-MAT_HUM_PIGNONHH[MAT_HUM_PIGNONHH$same_SBC==T,]
MAT_HUM_PIG_SBC$same <- MAT_HUM_PIG_SBC$H1_pattern == MAT_HUM_PIG_SBC$H2_pattern
MAT_HUM_PIG_SBC<-MAT_HUM_PIG_SBC[!is.na(MAT_HUM_PIG_SBC$same),]
MAT_HUM_PIG_SBC_BTN<-MAT_HUM_PIGNONHH[MAT_HUM_PIGNONHH$same_dis==F,]
MAT_HUM_PIG_SBC_BTN$same <- MAT_HUM_PIG_SBC_BTN$H1_pattern == MAT_HUM_PIG_SBC_BTN$H2_pattern
MAT_HUM_PIG_SBC_BTN<-MAT_HUM_PIG_SBC_BTN[!is.na(MAT_HUM_PIG_SBC_BTN$same),]
MAT_HUM_PIG_DIS<-MAT_HUM_PIG[MAT_HUM_PIG$same_dis!=T,]
MAT_HUM_PIG_DIS$same <- MAT_HUM_PIG_DIS$H1_pattern == MAT_HUM_PIG_DIS$H2_pattern
MAT_HUM_PIG_DIS<-MAT_HUM_PIG_DIS[!is.na(MAT_HUM_PIG_DIS$same),]
MAT_HUM_PIG_same_HH$LOCATION<-"SAME_HOUSEHOLD"
## P&F WITHIN SUBCOUNTY
MAT_HUM_PIG_SBC$LOCATION<-"SAME_SUBCOUNTY"
## P&F BETWEEN SUBCOUNTY BUT SAME DISTRICT
MAT_HUM_PIG_SBC_BTN$LOCATION<-"DIFF_SUBCOUNTY"
## P&FBETWEEN DISTRICS
MAT_HUM_PIG_DIS$LOCATION<-"DIFF_DISTRICT"
DB_MODEL_COMB<- rbind(MAT_HUM_PIG_same_HH,MAT_HUM_PIG_SBC,MAT_HUM_PIG_DIS) ## COMPARISONS COMBINED DB
names(DB_MODEL_COMB)[14]<-"Transmission"
DB_MODEL_COMB$Transmission<- as.character(DB_MODEL_COMB$Transmission)
DB_MODEL_COMB$Transmission[DB_MODEL_COMB$Transmission=="TRUE"]<-"Yes"
DB_MODEL_COMB$Transmission[DB_MODEL_COMB$Transmission=="FALSE"]<-"No"
DB_MODEL_COMB$Transmission<- as.factor(DB_MODEL_COMB$Transmission)
DB_MODEL_COMB$TIME_HOST1<-substr(DB_MODEL_COMB$HOST1,9,9)
DB_MODEL_COMB$TIME_HOST1<- as.numeric(DB_MODEL_COMB$TIME_HOST1)
DB_MODEL_COMB$TIME_HOST2<-substr(DB_MODEL_COMB$HOST2,9,9)
DB_MODEL_COMB$TIME_HOST2<- as.numeric(DB_MODEL_COMB$TIME_HOST2)## Warning: NAs introduced by coercionDB_MODEL_COMB$LAG<-DB_MODEL_COMB$TIME_HOST2-DB_MODEL_COMB$TIME_HOST1
DB_MODEL_COMB$LAG2<- abs(DB_MODEL_COMB$LAG)
DB_MODEL_COMB$LAG3<- factor(DB_MODEL_COMB$LAG2)
CONTS<-DB_MODEL_COMB[str_detect(DB_MODEL_COMB$HOST2, "CT"),]
DB_MODEL_COMB<-DB_MODEL_COMB[!DB_MODEL_COMB$HOST2%in%CONTS$HOST2,]
DB_MODEL_COMB_cons<-DB_MODEL_COMB[DB_MODEL_COMB$H1_pattern!="111",] ## conservative estimate with patterns that occur at less than 10%,
#so sharing these would be even more less likely
DB_MODEL_COMB_cons<-DB_MODEL_COMB_cons[DB_MODEL_COMB_cons$H2_pattern!="111",]
#*************************************************************************
## VISUALISING TRANSMISSION INFERENCE OUTPUT
#*************************************************************************
names(DB_MODEL_COMB)[15]<-"Comparison"
names(DB_MODEL_COMB)[20]<-"Time_Lag"
names(DB_MODEL_COMB_cons)[15]<-"Comparison"
names(DB_MODEL_COMB_cons)[20]<-"Time_Lag"
ggplot(DB_MODEL_COMB, aes(LAG2, color=Comparison)) + geom_density()
ggplot(DB_MODEL_COMB, aes(LAG2, color=Comparison)) + geom_density() +
facet_wrap(~Transmission) + theme_bw()
DB_MODEL_COMB %>%
group_by(Comparison,Time_Lag,Transmission) %>%
dplyr::count(Comparison,Time_Lag,Transmission) %>%
ggplot(aes(Time_Lag,log10(n), fill=Transmission)) + geom_bar(stat = "identity",position = "dodge") +
theme_bw() + facet_wrap(~Comparison)
DB_MODEL_COMB_cons %>%
group_by(Comparison,Time_Lag,Transmission) %>%
dplyr::count(Comparison,Time_Lag,Transmission) %>%
ggplot(aes(Time_Lag,log10(n), fill=Transmission)) + geom_bar(stat = "identity",position = "dodge") +
theme_bw() + facet_wrap(~Comparison)
DB_MODEL_COMB$Comparison<- as.character(DB_MODEL_COMB$Comparison)
DB_MODEL_COMB$Comparison[DB_MODEL_COMB$Comparison=="DIFF_DISTRICT"]<- "Between production system"
DB_MODEL_COMB$Comparison[DB_MODEL_COMB$Comparison=="DIFF_SUBCOUNTY"]<- "Farmers and pigs between Subcounties in same district"
DB_MODEL_COMB$Comparison[DB_MODEL_COMB$Comparison=="SAME_HOUSEHOLD"]<- "Within Farm"
DB_MODEL_COMB$Comparison[DB_MODEL_COMB$Comparison=="SAME_SUBCOUNTY"]<- "Within Subcounties"
DB_MODEL_COMB_cons$Comparison<- as.character(DB_MODEL_COMB_cons$Comparison)
DB_MODEL_COMB_cons$Comparison[DB_MODEL_COMB_cons$Comparison=="DIFF_DISTRICT"]<- "Between production system"
DB_MODEL_COMB_cons$Comparison[DB_MODEL_COMB_cons$Comparison=="DIFF_SUBCOUNTY"]<- "Farmers and pigs between Subcounties in same district"
DB_MODEL_COMB_cons$Comparison[DB_MODEL_COMB_cons$Comparison=="SAME_HOUSEHOLD"]<- "Within Subcounties"
DB_MODEL_COMB_cons$Comparison[DB_MODEL_COMB_cons$Comparison=="SAME_SUBCOUNTY"]<- "Within Farm"
names(DB_MODEL_COMB)[15]<-"Comparison"
names(DB_MODEL_COMB)[14]<-"Transmission"
names(DB_MODEL_COMB_cons)[15]<-"Comparison" ## conservative estimate based on pattern with less than 10% frequency
names(DB_MODEL_COMB_cons)[14]<-"Transmission"
prop.test(x = (68+330+66),
n = 1850,
p = (68+330+66)/1850)##
## 1-sample proportions test without continuity correction
##
## data: (68 + 330 + 66) out of 1850, null probability (68 + 330 + 66)/1850
## X-squared = 0, df = 1, p-value = 1
## alternative hypothesis: true p is not equal to 0.2508108
## 95 percent confidence interval:
## 0.2315880 0.2710664
## sample estimates:
## p
## 0.2508108plot_xtab(
x = DB_MODEL_COMB$Comparison,
grp = DB_MODEL_COMB$Transmission,
margin = "row",
bar.pos = "stack",
show.summary = TRUE,
coord.flip = TRUE) ## conservative transmission estimate
plot_xtab(
x = DB_MODEL_COMB_cons$Comparison,
grp = DB_MODEL_COMB_cons$Transmission,
margin = "row",
bar.pos = "stack",
show.summary = TRUE,
coord.flip = TRUE) ## Conservative transmission estimate## Warning in stats::chisq.test(ftab): Chi-squared approximation may be incorrect
tab_xtab(
var.row = DB_MODEL_COMB$Comparison,
var.col = DB_MODEL_COMB$Transmission,
show.row.prc = T)## Liberal transmission estimate| Comparison | Transmission | Total | |
|---|---|---|---|
| No | Yes | ||
|
Between production system |
530 88.9 % |
66 11.1 % |
596 100 % |
| Within Farm |
110 61.8 % |
68 38.2 % |
178 100 % |
| Within Subcounties |
746 69.3 % |
330 30.7 % |
1076 100 % |
| Total |
1386 74.9 % |
464 25.1 % |
1850 100 % |
χ2=96.422 · df=2 · Cramer’s V=0.228 · p=0.000 |
tab_xtab(
var.row = DB_MODEL_COMB_cons$Comparison,
var.col = DB_MODEL_COMB_cons$Transmission,
show.row.prc = T)## Conservative transmission estimate| Comparison | Transmission | Total | |
|---|---|---|---|
| No | Yes | ||
|
Between production system |
210 89 % |
26 11 % |
236 100 % |
| Within Farm |
200 94.3 % |
12 5.7 % |
212 100 % |
| Within Subcounties |
32 76.2 % |
10 23.8 % |
42 100 % |
| Total |
442 90.2 % |
48 9.8 % |
490 100 % |
χ2=13.836 · df=2 · Cramer’s V=0.168 · Fisher’s p=0.001 |
DB_MODEL_COMy<- DB_MODEL_COMB %>%
group_by(Comparison,Time_Lag,Transmission) %>%
dplyr::count(Time_Lag,Transmission) %>%
spread(key = Transmission, value = n)
DB_MODEL_COMz<- DB_MODEL_COMB_cons %>%
group_by(Comparison,Time_Lag,Transmission) %>%
dplyr::count(Time_Lag,Transmission) %>%
spread(key = Transmission, value = n)
DB_MODEL_COMy$Yes[is.na(DB_MODEL_COMy$Yes)]<-0
DB_MODEL_COMy$No[is.na(DB_MODEL_COMy$No)]<-0
DB_MODEL_COMz$Yes[is.na(DB_MODEL_COMz$Yes)]<-0
DB_MODEL_COMz$No[is.na(DB_MODEL_COMz$No)]<-0
DB_MODEL_COMy$Total<- DB_MODEL_COMy$No + DB_MODEL_COMy$Yes
DB_MODEL_COMy$Proportion<- DB_MODEL_COMy$Yes/DB_MODEL_COMy$Total
DB_MODEL_COMz$Total<- DB_MODEL_COMz$No + DB_MODEL_COMz$Yes
DB_MODEL_COMz$Proportion<- DB_MODEL_COMz$Yes/DB_MODEL_COMz$Total
ggplot(DB_MODEL_COMy,aes(Time_Lag,Proportion,fill=Comparison)) +
geom_bar(stat = "identity",position = "dodge")
ggplot(DB_MODEL_COMz,aes(Time_Lag,Proportion,fill=Comparison)) +
geom_bar(stat = "identity",position = "dodge")
ggplot(DB_MODEL_COMy,aes(Time_Lag,Proportion,fill=Comparison)) +
geom_bar(stat = "identity",position = "dodge") + theme_bw()
ggplot(DB_MODEL_COMz,aes(Time_Lag,Proportion,fill=Comparison)) +
geom_bar(stat = "identity",position = "dodge") + theme_bw()
Prop_AVE<-DB_MODEL_COMy %>%
group_by(Time_Lag,Comparison) %>%
mutate(Aveper= mean(Proportion))
Prop_AVE_con<-DB_MODEL_COMz %>%
group_by(Time_Lag,Comparison) %>%
mutate(Aveper= mean(Proportion))
ggplot(Prop_AVE,aes(Time_Lag,Aveper,color=Comparison)) + geom_point(size=3) +
theme_bw() 
ggplot(Prop_AVE_con,aes(Time_Lag,Aveper,color=Comparison)) + geom_point(size=3) +
theme_bw() 
ggplot(Prop_AVE, aes(x=Time_Lag, y=Aveper, group=Comparison, color=Comparison)) +
geom_line(size=1)+
geom_point(size=3) +
scale_fill_manual(values=c("#000080","#FC4E07","#800000")) + theme_bw() +
ylab("Probability of transmission event") +
xlab("Time lag")## Liberal transmission estimate## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ggplot(Prop_AVE_con, aes(x=Time_Lag, y=Aveper, group=Comparison, color=Comparison)) +
geom_line(size=1)+
geom_point(size=3) +
scale_fill_manual(values=c("#000080","#FC4E07","#800000")) + theme_bw() +
ylab("Probability of transmission event") +
xlab("Time lag")## Cons transmission estimate
DB_MODEL_COMy$Time_Lag <-as.numeric(DB_MODEL_COMy$Time_Lag)
DB_MODEL_COMy$merger<- paste(DB_MODEL_COMy$Comparison,DB_MODEL_COMy$Time_Lag,sep = "_")
ggplot(DB_MODEL_COMy,aes(Time_Lag,Proportion,color=Comparison)) +
geom_smooth(size=1.5) + theme_bw() + xlab("Time lag") + ylab("Percentage MDR share by Farmer and Pig")## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : span too small. fewer data values than degrees of freedom.## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 0.98## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 2.02## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 4.0804## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small. fewer
## data values than degrees of freedom.## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 0.98## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.02## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0804## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : span too small. fewer data values than degrees of freedom.## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 0.985## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 2.015## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 4.0602## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small. fewer
## data values than degrees of freedom.## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 0.985## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.015## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0602## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : span too small. fewer data values than degrees of freedom.## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 0.985## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 2.015## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 4.0602## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small. fewer
## data values than degrees of freedom.## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 0.985## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 2.015## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0602## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
##************************************
DB_MODEL_COM3<- DB_MODEL_COMB[DB_MODEL_COMB$Comparison!="Farmers and pigs between Subcounties in same district",]
names(DB_MODEL_COMB)[15]<-"Comparison"
names(DB_MODEL_COMB)[20]<-"Time_Lag"
DB_MODEL_COMB_cons$Comparison[DB_MODEL_COMB_cons$Comparison=="DIFF_DISTRICT"]<- "Between production system"
DB_MODEL_COMB_cons$Comparison[DB_MODEL_COMB_cons$Comparison=="DIFF_SUBCOUNTY"]<- "Farmers and pigs between Subcounties in same district"
DB_MODEL_COMB_cons$Comparison[DB_MODEL_COMB_cons$Comparison=="SAME_HOUSEHOLD"]<- "Within Farm"
DB_MODEL_COMB_cons$Comparison[DB_MODEL_COMB_cons$Comparison=="SAME_SUBCOUNTY"]<- "Within Subcounties"
DB_MODEL_COMB$Comparison<- factor(DB_MODEL_COMB$Comparison,levels = c("Between production system",
"Within Subcounties",
"Within Farm"))
## PLOT COMPARISONS
DB_MODEL_COMB %>%
plot_gpt(x = Transmission, y = Time_Lag, grp = Comparison) 
DB_MODEL_COMB_cons %>%
plot_gpt(x = Transmission, y = Time_Lag, grp = Comparison) 
#************************************************************************
### TRANSMISSION INFERENCES MODEL WITHOUT THE GENES
#************************************************************************
Model <- glm(Transmission ~ Comparison + Time_Lag , data = DB_MODEL_COMB, family = "binomial")
tab_model(Model)| Transmission | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 0.32 | 0.23 – 0.45 | <0.001 |
|
Comparison [Within Subcounties] |
3.43 | 2.53 – 4.72 | <0.001 |
| Comparison [Within Farm] | 3.08 | 1.98 – 4.79 | <0.001 |
| Time Lag [1] | 0.24 | 0.18 – 0.32 | <0.001 |
| Time Lag [2] | 0.40 | 0.29 – 0.53 | <0.001 |
| Time Lag [3] | 0.13 | 0.06 – 0.24 | <0.001 |
| Time Lag [4] | 1.20 | 0.52 – 2.62 | 0.652 |
| Observations | 1850 | ||
| R2 Tjur | 0.128 | ||
DB_MODEL_COMB_cons$Time_Lag[DB_MODEL_COMB_cons$Time_Lag==3]<-4
Model_conz <- glm(Transmission ~ Comparison + Time_Lag , data = DB_MODEL_COMB_cons, family = "binomial")
tab_model(Model)| Transmission | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 0.32 | 0.23 – 0.45 | <0.001 |
|
Comparison [Within Subcounties] |
3.43 | 2.53 – 4.72 | <0.001 |
| Comparison [Within Farm] | 3.08 | 1.98 – 4.79 | <0.001 |
| Time Lag [1] | 0.24 | 0.18 – 0.32 | <0.001 |
| Time Lag [2] | 0.40 | 0.29 – 0.53 | <0.001 |
| Time Lag [3] | 0.13 | 0.06 – 0.24 | <0.001 |
| Time Lag [4] | 1.20 | 0.52 – 2.62 | 0.652 |
| Observations | 1850 | ||
| R2 Tjur | 0.128 | ||
tab_model(Model_conz)| Transmission | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 0.54 | 0.25 – 1.11 | 0.097 |
| Comparison [Within Farm] | 0.45 | 0.21 – 0.92 | 0.035 |
|
Comparison [Within Subcounties] |
0.81 | 0.29 – 2.15 | 0.684 |
| Time Lag [1] | 0.17 | 0.07 – 0.41 | <0.001 |
| Time Lag [2] | 0.18 | 0.07 – 0.46 | <0.001 |
| Time Lag [4] | 0.19 | 0.05 – 0.59 | 0.007 |
| Observations | 490 | ||
| R2 Tjur | 0.077 | ||
plot_model(Model, show.values = TRUE, width = 0.1) + theme_bw()
plot_model(Model_conz, show.values = TRUE, width = 0.1) + theme_bw() ## conservative is less precise
plot_model(Model,type="pred",
terms=c("Time_Lag","Comparison"),show.legend = F) + theme_bw()
plot_model(Model_conz,type="pred",
terms=c("Time_Lag","Comparison"),show.legend = F) + theme_bw()
plot_model(Model,type="pred",
terms=c("Time_Lag","Comparison"),show.legend = T) + theme_bw()
plot_model(Model_conz,type="pred",
terms=c("Time_Lag","Comparison"),show.legend = T) + theme_bw()
#************************************************************************
## TRANSMISSION INFERENCES MODELS WITH THE GENES
#***********************************************************************
GENE_COUNT<-QPCRDATA
names(GENE_COUNT)[9]<-"teb_b_norm"
names(GENE_COUNT)[10]<-"tet_q_norm"
names(GENE_COUNT)[11]<-"erm_b_norm"
names(GENE_COUNT)[12]<-"dfra_1_norm"
DB_MODEL_COMB$tetQ_H1<-GENE_COUNT[match(DB_MODEL_COMB$HOST1,GENE_COUNT$sample_id),]$tet_q_norm
DB_MODEL_COMB$tetB_H1<-GENE_COUNT[match(DB_MODEL_COMB$HOST1,GENE_COUNT$sample_id),]$teb_b_norm
DB_MODEL_COMB$ermB_H1<-GENE_COUNT[match(DB_MODEL_COMB$HOST1,GENE_COUNT$sample_id),]$erm_b_norm
DB_MODEL_COMB$dfra_H1<-GENE_COUNT[match(DB_MODEL_COMB$HOST1,GENE_COUNT$sample_id),]$dfra_1_norm
DB_MODEL_COMB$tetQ_H2<-GENE_COUNT[match(DB_MODEL_COMB$HOST2,GENE_COUNT$sample_id),]$tet_q_norm
DB_MODEL_COMB$tetB_H2<-GENE_COUNT[match(DB_MODEL_COMB$HOST2,GENE_COUNT$sample_id),]$teb_b_norm
DB_MODEL_COMB$ermB_H2<-GENE_COUNT[match(DB_MODEL_COMB$HOST2,GENE_COUNT$sample_id),]$erm_b_norm
DB_MODEL_COMB$dfra_H2<-GENE_COUNT[match(DB_MODEL_COMB$HOST2,GENE_COUNT$sample_id),]$dfra_1_norm
DB_MODEL_COM4<-DB_MODEL_COMB[DB_MODEL_COMB$tetQ_H2!=Inf,]
DB_MODEL_COM4$tetQ_ratio<- DB_MODEL_COM4$tetQ_H1/ DB_MODEL_COM4$tetQ_H2
DB_MODEL_COM4$tetB_ratio<- DB_MODEL_COM4$tetB_H1/ DB_MODEL_COM4$tetB_H2
DB_MODEL_COM4$ermB_ratio<- DB_MODEL_COM4$ermB_H1/ DB_MODEL_COM4$ermB_H2
DB_MODEL_COM4$dfra_ratio<- DB_MODEL_COM4$dfra_H1/ DB_MODEL_COM4$dfra_H2
DB_MODEL_COM4<- DB_MODEL_COM4[ DB_MODEL_COM4$tetB_ratio!=Inf,]
DB_MODEL_COM4<- DB_MODEL_COM4[ DB_MODEL_COM4$ermB_ratio!=Inf,]
DB_MODEL_COM4<- DB_MODEL_COM4[ DB_MODEL_COM4$dfra_ratio!=Inf,]
DB_MODEL_COM4$lgtetq<- log10( DB_MODEL_COM4$tetQ_ratio)
DB_MODEL_COM4$lgtetb<- log10( DB_MODEL_COM4$tetB_ratio)
DB_MODEL_COM4$lgwemB<- log10( DB_MODEL_COM4$ermB_ratio)
DB_MODEL_COM4$lgdfra<- log10( DB_MODEL_COM4$dfra_ratio)
DB_MODEL_COM4 %>%
group_by(Comparison) %>%
plot_frq(lgtetq, type = "density", show.mean = TRUE, normal.curve = TRUE) %>%
plot_grid()## Warning in plot_grid(.): Not enough tags labels in list. Using letters instead.## Warning: `stat(density)` was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## ℹ The deprecated feature was likely used in the sjPlot package.
## Please report the issue at <https://github.com/strengejacke/sjPlot/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
DB_MODEL_COM4_x<- DB_MODEL_COM4[ DB_MODEL_COM4$lgdfra!=(-Inf),]
DB_MODEL_COM4_x<- DB_MODEL_COM4_x[ DB_MODEL_COM4_x$lgwemB!=(-Inf),]
DB_MODEL_COM4_x<- DB_MODEL_COM4_x[ DB_MODEL_COM4_x$lgtetb!=(-Inf),]
DB_MODEL_COM4_x<- DB_MODEL_COM4_x[ DB_MODEL_COM4_x$lgtetq!=(-Inf),]
## The transmission inference model including the effects of AMRgenes counts
Model2 <- glm(Transmission ~ Comparison + Time_Lag + lgdfra + lgwemB + lgtetb + lgtetq, data = DB_MODEL_COM4_x, family = "binomial")
tab_model(Model2)| Transmission | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 0.31 | 0.21 – 0.47 | <0.001 |
|
Comparison [Within Subcounties] |
3.13 | 2.20 – 4.53 | <0.001 |
| Comparison [Within Farm] | 3.16 | 1.83 – 5.43 | <0.001 |
| Time Lag [1] | 0.23 | 0.16 – 0.33 | <0.001 |
| Time Lag [2] | 0.54 | 0.37 – 0.79 | 0.001 |
| Time Lag [3] | 0.12 | 0.05 – 0.25 | <0.001 |
| Time Lag [4] | 0.98 | 0.37 – 2.41 | 0.972 |
| lgdfra | 1.24 | 1.12 – 1.38 | <0.001 |
| lgwemB | 0.79 | 0.72 – 0.87 | <0.001 |
| lgtetb | 0.99 | 0.89 – 1.10 | 0.869 |
| lgtetq | 1.25 | 1.09 – 1.43 | 0.002 |
| Observations | 1276 | ||
| R2 Tjur | 0.156 | ||
Model2 %>% ggpredict() ## Data were 'prettified'. Consider using `terms="lgdfra [all]"` to get
## smooth plots.
## Data were 'prettified'. Consider using `terms="lgwemB [all]"` to get
## smooth plots.
## Data were 'prettified'. Consider using `terms="lgtetb [all]"` to get
## smooth plots.
## Data were 'prettified'. Consider using `terms="lgtetq [all]"` to get
## smooth plots.## $Comparison
## # Predicted probabilities of Transmission
##
## Comparison | Predicted | 95% CI
## ----------------------------------------------------
## Between production system | 0.23 | [0.17, 0.31]
## Within Subcounties | 0.48 | [0.41, 0.55]
## Within Farm | 0.48 | [0.38, 0.59]
##
## Adjusted for:
## * Time_Lag = 0
## * lgdfra = -0.03
## * lgwemB = 0.17
## * lgtetb = 0.37
## * lgtetq = -0.05
##
## $Time_Lag
## # Predicted probabilities of Transmission
##
## Time_Lag | Predicted | 95% CI
## -----------------------------------
## 0 | 0.23 | [0.17, 0.31]
## 1 | 0.06 | [0.04, 0.09]
## 2 | 0.14 | [0.10, 0.19]
## 3 | 0.03 | [0.02, 0.07]
## 4 | 0.23 | [0.11, 0.40]
##
## Adjusted for:
## * Comparison = Between production system
## * lgdfra = -0.03
## * lgwemB = 0.17
## * lgtetb = 0.37
## * lgtetq = -0.05
##
## $lgdfra
## # Predicted probabilities of Transmission
##
## lgdfra | Predicted | 95% CI
## ---------------------------------
## -6 | 0.07 | [0.04, 0.14]
## -4 | 0.11 | [0.07, 0.18]
## -2 | 0.16 | [0.11, 0.23]
## 0 | 0.23 | [0.17, 0.31]
## 2 | 0.32 | [0.23, 0.42]
## 4 | 0.42 | [0.29, 0.56]
## 6 | 0.52 | [0.35, 0.70]
##
## Adjusted for:
## * Comparison = Between production system
## * Time_Lag = 0
## * lgwemB = 0.17
## * lgtetb = 0.37
## * lgtetq = -0.05
##
## $lgwemB
## # Predicted probabilities of Transmission
##
## lgwemB | Predicted | 95% CI
## ---------------------------------
## -10 | 0.76 | [0.54, 0.89]
## -5 | 0.49 | [0.35, 0.64]
## 0 | 0.24 | [0.17, 0.32]
## 5 | 0.09 | [0.05, 0.15]
## 10 | 0.03 | [0.01, 0.08]
##
## Adjusted for:
## * Comparison = Between production system
## * Time_Lag = 0
## * lgdfra = -0.03
## * lgtetb = 0.37
## * lgtetq = -0.05
##
## $lgtetb
## # Predicted probabilities of Transmission
##
## lgtetb | Predicted | 95% CI
## ---------------------------------
## -6 | 0.24 | [0.13, 0.41]
## -4 | 0.24 | [0.14, 0.36]
## -2 | 0.23 | [0.16, 0.33]
## 0 | 0.23 | [0.17, 0.31]
## 2 | 0.23 | [0.16, 0.31]
## 4 | 0.22 | [0.14, 0.33]
## 6 | 0.22 | [0.12, 0.37]
##
## Adjusted for:
## * Comparison = Between production system
## * Time_Lag = 0
## * lgdfra = -0.03
## * lgwemB = 0.17
## * lgtetq = -0.05
##
## $lgtetq
## # Predicted probabilities of Transmission
##
## lgtetq | Predicted | 95% CI
## ---------------------------------
## -6 | 0.07 | [0.03, 0.17]
## -4 | 0.11 | [0.06, 0.20]
## -2 | 0.16 | [0.11, 0.24]
## 0 | 0.23 | [0.17, 0.31]
## 2 | 0.32 | [0.22, 0.43]
## 4 | 0.42 | [0.27, 0.59]
##
## Adjusted for:
## * Comparison = Between production system
## * Time_Lag = 0
## * lgdfra = -0.03
## * lgwemB = 0.17
## * lgtetb = 0.37
##
## attr(,"class")
## [1] "ggalleffects" "list"
## attr(,"model.name")
## [1] "."drop1(Model2,test="Chisq") ## Single term deletions
##
## Model:
## Transmission ~ Comparison + Time_Lag + lgdfra + lgwemB + lgtetb +
## lgtetq
## Df Deviance AIC LRT Pr(>Chi)
## <none> 1192.3 1214.3
## Comparison 2 1237.5 1255.5 45.139 1.578e-10 ***
## Time_Lag 4 1279.0 1293.0 86.725 < 2.2e-16 ***
## lgdfra 1 1210.5 1230.5 18.191 1.999e-05 ***
## lgwemB 1 1217.9 1237.9 25.553 4.304e-07 ***
## lgtetb 1 1192.3 1212.3 0.027 0.868639
## lgtetq 1 1202.4 1222.4 10.061 0.001514 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Visualising of the model output
plot_model(Model2, show.values = TRUE, width = 0.1) + theme_bw()
plot_model(Model2,type="pred",
terms=c("Time_Lag","Comparison"),show.legend = F) + theme_bw()
plot_model(Model2,type="pred",
terms=c("lgtetq [all]","Comparison"),show.legend = F) + theme_bw()
plot_model(Model2,type="pred",
terms=c("lgwemB [all]","Comparison"),show.legend = F) + theme_bw()
plot_model(Model2,type="pred",
terms=c("lgdfra [all]","Comparison"),show.legend = F) + theme_bw()
plot_model(Model2,type="pred",
terms=c("lgtetq [all]","Comparison"),show.legend = F) + theme_bw()
plot_model(Model2,type="pred",
terms=c("lgtetq [all]","Comparison","Time_Lag"),show.legend = F) + theme_bw()
plot_model(Model2,type="pred",
terms=c("lgdfra [all]","Comparison","Time_Lag"),show.legend = F) + theme_bw()
DB_MODEL_COM5<- DB_MODEL_COM4[ DB_MODEL_COM4$lgtetb!=Inf,]
DB_MODEL_COM5<-DB_MODEL_COM5[DB_MODEL_COM5$lgtetb!=(-Inf),]
DB_MODEL_COM5<-DB_MODEL_COM5[DB_MODEL_COM5$lgwemB!=Inf,]
DB_MODEL_COM5<-DB_MODEL_COM5[DB_MODEL_COM5$lgwemB!=(-Inf),]
DB_MODEL_COM5<-DB_MODEL_COM5[DB_MODEL_COM5$lgdfra!=(-Inf),]
DB_MODEL_COM5<-DB_MODEL_COM5[DB_MODEL_COM5$lgdfra!=Inf,]
DB_MODEL_COM5<-DB_MODEL_COM5[DB_MODEL_COM5$lgtetq!=(-Inf),]
DB_MODEL_COM5<-DB_MODEL_COM5[DB_MODEL_COM5$lgtetq!=Inf,]
DB_MODEL_COM5<-DB_MODEL_COM5[!is.na(DB_MODEL_COM5$lgdfra),]
names(DB_MODEL_COM5)[33]<-"tetQ"
names(DB_MODEL_COM5)[34]<-"tetB"
names(DB_MODEL_COM5)[35]<-"ermB"
names(DB_MODEL_COM5)[36]<-"defra1"
Model3 <- glm(Transmission ~ Comparison + Time_Lag + tetQ + tetB + ermB + defra1, data = DB_MODEL_COM5, family = "binomial")
tab_model(Model3)| Transmission | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 0.31 | 0.21 – 0.47 | <0.001 |
|
Comparison [Within Subcounties] |
3.13 | 2.20 – 4.53 | <0.001 |
| Comparison [Within Farm] | 3.16 | 1.83 – 5.43 | <0.001 |
| Time Lag [1] | 0.23 | 0.16 – 0.33 | <0.001 |
| Time Lag [2] | 0.54 | 0.37 – 0.79 | 0.001 |
| Time Lag [3] | 0.12 | 0.05 – 0.25 | <0.001 |
| Time Lag [4] | 0.98 | 0.37 – 2.41 | 0.972 |
| tetQ | 1.25 | 1.09 – 1.43 | 0.002 |
| tetB | 0.99 | 0.89 – 1.10 | 0.869 |
| ermB | 0.79 | 0.72 – 0.87 | <0.001 |
| defra1 | 1.24 | 1.12 – 1.38 | <0.001 |
| Observations | 1276 | ||
| R2 Tjur | 0.156 | ||
ERMBPLOT<-plot_model(Model3,type="pred", title = "",
terms=c("ermB [all]","Comparison"),show.legend = F) + theme_bw()
TETBPLOT<-plot_model(Model3,type="pred", title = "",
terms=c("tetB [all]","Comparison"),show.legend = F) + theme_bw()
DFRAPLOT<-plot_model(Model3,type="pred", title = "",
terms=c("defra1 [all]","Comparison"),show.legend = F) + theme_bw()
TETQPLOT<-plot_model(Model3,type="pred", title = "",
terms=c("tetQ [all]","Comparison"),show.legend = F) + theme_bw()
TIMEPLOT<-plot_model(Model3,type="pred", title = "",
terms=c("Time_Lag","Comparison"),show.legend = F) + theme_bw()
multiplot(ERMBPLOT,DFRAPLOT,TETQPLOT,TETBPLOT,cols = 4 )
plot_model(Model3, show.values = TRUE, width = 0.1) + theme_bw()
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.