PART 3 - Subclinical mastitis analysis 




library(knitr)
library(readr)
library(ggplot2)
library(survminer)
## Loading required package: ggpubr
## Loading required package: magrittr
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.0.2
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(survival)
library(broom)
library(forestplot)
## Loading required package: grid
## Loading required package: checkmate
library(table1)
## 
## Attaching package: 'table1'
## The following objects are masked from 'package:base':
## 
##     units, units<-




knitr::opts_chunk$set(echo=TRUE, warning=FALSE, message=FALSE)


#Set working directory ----
#PC
#wd <- function(){setwd("C:/Users/rowe0122/Dropbox/R backup/2020/IMI cohort"")}
#Mac
wd <- function(){setwd("~/Dropbox/R backup/2020/IMI cohort")}


#Set WD ----
wd()

load(file="IMIcohortDHI.Rdata")
load(file="IMIcohortfull.Rdata")
load(file="IMIcohortpostcalving.Rdata")
ABX <- read_csv("ABXkeep.csv", col_types = cols(X1 = col_skip()))
## Warning: Missing column names filled in: 'X1' [1]
cow$SSLO <- 0
cow$SSLO[cow$StrepALL >=1 | cow$LactoALL >= 1 | cow$Entero >=1 | cow$Aerococ >=1] <- 1

postcalve$SSLO <- 0
postcalve$SSLO[postcalve$StrepALL >=1 | postcalve$LactoALL >= 1 | postcalve$Entero >=1 | postcalve$Aerococ >=1] <- 1

dhi$SSLO <- 0
dhi$SSLO[dhi$StrepALL >=1 | dhi$LactoALL >= 1 | dhi$Entero >=1 | dhi$Aerococ >=1] <- 1

postcalve <- merge(ABX,postcalve,by=c("HerdID","Cow"))
dhi <- merge(ABX,dhi,by=c("HerdID","Cow"))

Create functions to simplify the modelling syntax

Sort data

dhi <- dhi[order(dhi$TestNo),]
dhi <- dhi[order(dhi$Herdcow),]
dhi <- dhi[order(dhi$HerdID),]
dhi2 <- dhi
dhi <- dhi %>% filter(!is.na(MYDO))
library(lme4)
library(geepack)
library(emmeans)

Outcome 4 - Subclinical mastitis (1-120 DIM) (SCC > 200,000 cells/ml)

Model building plan

Model type: Generalized estimating equations (~GLM). Clustering indicated at the herd level. Working covariance structure = independence.

Step 1: Look for interactions with DIM by fitting a interaction model (Exposure*DIM) and plotting infection risk using estimated marginal means. If there is no biologically relevant interaction, then a main effects model will be reported.

Step 2: Create model with confounders as covariates (using DAG from earlier) +/- interaction term.

Step 3: Report final model



Sort data

dhi <- dhi[order(dhi$TestNo),]
dhi <- dhi[order(dhi$Herdcow),]
dhi <- dhi[order(dhi$HerdID),]
dhi <- dhi2

Create functions to simplify the modelling syntax

library(emmeans)

gee.plot <- function(expo,exponame){
  mm0 <-  geeglm(SCM ~ factor(expo)*dhi$TestDIMcat20, data=dhi, id = HerdID,family=binomial(link="identity"))
  atx <- unique(dhi$TestDIMcat20)
  emm <- emmeans(mm0, ~expo*TestDIMcat20) %>% tidy()
  emm$LCL <- emm$asymp.LCL
  emm$UCL <- emm$asymp.UCL
  emm$Exposure <- emm$expo %>% as.factor
  emm <- emm %>% subset(select=c(Exposure,TestDIMcat20,estimate,LCL,UCL))
  curve <- ggplot(emm) + coord_cartesian(ylim = (c(0,.5))) + aes(x=TestDIMcat20, y=estimate, group=Exposure, colour=Exposure) +
    labs(y = "Subclinical mastitis risk", x = "Days in milk") + ggtitle(exponame) + geom_point(size=3) + geom_ribbon(aes(ymin=emm$LCL, ymax=emm$UCL,colour=Exposure,fill=Exposure), linetype=0, alpha=0.1) + geom_line(aes(colour=Exposure,linetype=Exposure),size=1.1) + theme(panel.border = element_rect(colour = "black", fill=NA, size=1),axis.text=element_text(size=12,family="Times"),axis.title=element_text(size=14,face="bold",family="Times"),panel.background = element_rect(fill = "white", colour = "white",size = 0.5, linetype = "solid"),panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "grey"), panel.grid.minor = element_line(size = 0.25, linetype = 'solid',colour = "grey"),legend.position="bottom",legend.text = element_text(colour="black", size=12,face="bold",family="Times")) + scale_color_brewer(palette="Set1") 
  curve
}



geefunc <- function(expo,exposure,exponame,expo2,Other){
  OtherCov <- "+ factor(ABX) + factor(Summer) + factor(Dryfreestall) + factor(calvedalone) + factor(Bedtype) + factor(Parity) + factor(Breed) + factor(East) + factor(TestDIMcat20)"
  OtherIMI <- c(expo,Other)
  Formula <- formula(paste("SCM ~  ",paste(OtherIMI, collapse=" + "),OtherCov))
  mm1 <-  geeglm(Formula, id=HerdID,family=poisson(link="log"), data=dhi)
  mm0 <- mm1 %>% tidy(conf.int=T,exp=T) %>% select(estimate,conf.low,conf.high) %>% slice(2)
  mm0$term <- exponame
  mm0$Exposure <- mm0$term
  mm0$term <- NULL
  
  a <- table(expo2,dhi$SCM) %>% as.data.frame
  names(a) <- c("Exp","CM","Freq")
  b <- a %>% filter(CM==0)
  c <- a %>% filter(CM==1)
  names(b) <- c("Exp","CM","Freq.CM0")
  names(c) <- c("Exp","CM","Freq.CM1")
  b$CM <- NULL
  c$CM <- NULL
  a <- merge(b,c)
  a$total <- a$Freq.CM0 + a$Freq.CM1
  a$Freq.CM0 <- NULL
  a$Percent <- a$Freq.CM1/a$total
  a$Risk <- paste0(sprintf("%.1f",round(100*a$Percent, 1)), sep="","% (",a$Freq.CM1,"/",a$total,")")
  a <- a %>% select(Exp,Risk)
  b <- a %>% filter(Exp==0)
  c <- a %>% filter(Exp==1)
  b$RiskExp0 <- b$Risk
  b <- b %>% select(RiskExp0)
  c$RiskExp1 <- c$Risk
  c <- c %>% select(RiskExp1)
  a <- merge(b,c)
  a$Exposure <- exponame
  mm0 <- merge(a,mm0,by=c("Exposure"))
  mm0$RR <- paste0(sprintf("%.1f",round(mm0$estimate,1))," (",sprintf("%.1f",round(mm0$conf.low,1)),
                   ", ",sprintf("%.1f",round(mm0$conf.high,1)),")")
  mm0 %>% select(Exposure,RiskExp0,RiskExp1,RR,estimate,conf.low,conf.high)

}

Exposure = Staphylococcus aureus

Run model and check residuals

gee.plot(dhi$Staaur, "Staphylococcus aureus")

#Staaur <- geefunc(dhi$Staaur,"Staphylococcus aureus")
Staaur <- geefunc(expo="factor(Staaur)",
          exposure="Staaur",
          exponame="Staphylococcus aureus",
          expo2=dhi$Staaur,
          Other=c("Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","Staspp",
                          "LactoALL","StrepALL",
                          "Bac","Coryn"))

Staaur

Exposure = Non aureus Staphylococcus sp. (whole group)

Run model and check residuals

gee.plot(dhi$NAS,"Non aureus Staphylococcus sp.")

#NAS <- geefunc(dhi$NAS,"Non aureus Staphylococcus sp.")

NAS<- geefunc("factor(NAS)",
        exposure="NAS",
        exponame="Non aureus Staphylococcus sp.",
        expo2=dhi$NAS,
        Other=c("Aerococ","Entero","LactoALL","StrepALL",
                "Bac","Coryn"))
NAS



Exposure = Staphylococcus chromogenes

Run model and check residuals

gee.plot(dhi$Stachr,"Staphylococcus chromogenes")

#Stachr <- geefunc(dhi$Stachr,"Staphylococcus chromogenes")
Stachr <- geefunc(expo="factor(Stachr)",
        exposure="Stachr",
        exponame="Staphylococcus chromogenes",
                expo2=dhi$Stachr ,
        Other=c("Staepi","Stahaem","Stasci","Stasim","Staxyl","Staspp",
                "Aerococ","Entero","LactoALL","StrepALL",
                "Bac","Coryn"))
Stachr



Exposure = Staphylococcus epidermidis

Run model and check residuals

gee.plot(dhi$Staepi,"Staphylococcus epidermidis")

#Staepi <- geefunc(dhi$Staepi,"Staphylococcus epidermidis")
Staepi <- geefunc("factor(Staepi)",
        exposure="Staepi",
        exponame="Staphylococcus epidermidis",
                        expo2=dhi$Staepi ,
        Other=c("Stachr","Stahaem","Stasci","Stasim","Staxyl","Staspp",
                "Aerococ","Entero","LactoALL","StrepALL",
                "Bac","Coryn"))
Staepi



Exposure = Staphylococcus haemolyticus

Run model and check residuals

gee.plot(dhi$Stahaem,"Staphylococcus haemolyticus")

#Stahaem <- geefunc(dhi$Stahaem,"Staphylococcus haemolyticus")
Stahaem <- geefunc("factor(Stahaem)",
        exposure="Stahaem",
        exponame="Staphylococcus haemolyticus",
                        expo2=dhi$Stahaem ,
        Other=c("Stachr","Staepi","Stasci","Stasim","Staxyl","Staspp",
                "Aerococ","Entero","LactoALL","StrepALL",
                "Bac","Coryn"))
Stahaem



Exposure = Staphylococcus sciuri

Run model and check residuals

#gee.plot(dhi$Stasci,"Staphylococcus sciuri")
#Stasci <- geefunc(dhi$Stasci,"Staphylococcus sciuri")
Stasci <- geefunc("factor(Stasci)",
        exposure="Stasci",
        exponame="Staphylococcus sciuri",
                        expo2=dhi$Stasci ,
        Other=c("Stachr","Staepi","Stahaem","Stasim","Staxyl","Staspp",
                "Aerococ","Entero","LactoALL","StrepALL",
                "Bac","Coryn"))
Stasci



Exposure = Staphylococcus simulans

Run model and check residuals

gee.plot(dhi$Stasim,"Staphylococcus simulans")

Stasim <- geefunc("factor(Stasim)",
        exposure="Stasim",
        exponame="Staphylococcus simulans",
                        expo2=dhi$Stasim ,
        Other=c("Stachr","Staepi","Stahaem","Stasci","Staxyl","Staspp",
                "Aerococ","Entero","LactoALL","StrepALL",
                "Bac","Coryn"))
#Stasim <- geefunc(dhi$Stasim,"Staphylococcus simulans")
Stasim



Exposure = Staphylococcus xylosus

Run model and check residuals

gee.plot(dhi$Staxyl,"Staphylococcus xylosus")

#Staxyl <- geefunc(dhi$Staxyl,"Staphylococcus xylosus")
Staxyl <- geefunc("factor(Staxyl)",
        exposure="Staxyl",
        exponame="Staphylococcus xylosus",
                        expo2=dhi$Staxyl ,
        Other=c("Stachr","Staepi","Stahaem","Stasci","Stasim","Staspp",
                "Aerococ","Entero","LactoALL","StrepALL",
                "Bac","Coryn"))
Staxyl



Exposure = Staphylococcus species (i.e. species unable to be determined using MALDI-TOF)

Run model and check residuals

gee.plot(dhi$Staspp,"Staphylococcus species")

#Staspp <- geefunc(dhi$Staspp,"Staphylococcus species")
Staspp <- geefunc("factor(Staspp)",
        exposure="Staspp",
        exponame="Staphylococcus sp.",
                        expo2=dhi$Staspp ,
        Other=c("Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl",
                "Aerococ","Entero","LactoALL","StrepALL",
                "Bac","Coryn"))
Staspp



Exposure = SSLO

Run model and check residuals

gee.plot(dhi$SSLO,"SSLO")

#SSLO <- geefunc(dhi$SSLO,"SSLO")
SSLO <- geefunc("factor(SSLO)",
        exposure="SSLO",
        exponame="SSLO",
                        expo2=dhi$SSLO ,
        Other=c("Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","Staspp",
                "Bac","Coryn"))
SSLO



Exposure = Aerococcus species

Run model and check residuals

gee.plot(dhi$Aerococ,"Aerococcus sp.")

#Aerococ <- geefunc(dhi$Aerococ,"Aerococcus sp.")
Aerococ <- geefunc("factor(Aerococ)",
        exposure="Aerococ",
        exponame="Aerococcus spp.",
                        expo2=dhi$Aerococ ,
        Other=c("Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","Staspp",
                "Entero","LactoALL","StrepALL",
                "Bac","Coryn"))
Aerococ



Exposure = Enterococcus species

Run model and check residuals

gee.plot(dhi$Entero,"Enterococcus sp.")

#Entero <- geefunc(dhi$Entero,"Enterococcus sp.")
Entero <- geefunc("factor(Entero)",
        exposure="Entero",
        exponame="Enterococcus spp.",
                        expo2=dhi$Entero ,
        Other=c("Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","Staspp",
                "Aerococ","LactoALL","StrepALL",
                "Bac","Coryn"))
Entero



Exposure = Lactococcus (all species combined)

Run model and check residuals

gee.plot(dhi$LactoALL,"Lactococcus spp.")

#LactoALL <- geefunc(dhi$LactoALL,"Lactococcus spp.")
LactoALL <- geefunc("factor(LactoALL)",
        exposure="LactoALL",
        exponame="Lactococcus spp.",
                        expo2=dhi$LactoALL ,
        Other=c("Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","Staspp",
                "Aerococ","Entero","StrepALL",
                "Bac","Coryn"))
LactoALL



Exposure = Lactococcus garviae

Run model and check residuals

gee.plot(dhi$Lactogar,"Lactococcus garviae")

#Lactogar <- geefunc(dhi$Lactogar,"Lactococcus garviae")
Lactogar <- geefunc("factor(Lactogar)",
        exposure="Lactogar",
        exponame="Lactococcus garviae",
                        expo2=dhi$Lactogar ,
        Other=c("Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","Staspp",
                "Aerococ","Entero","Lactolac","Lactosp","StrepALL",
                "Bac","Coryn"))
Lactogar



Exposure = Lactococcus lactis

Run model and check residuals

gee.plot(dhi$Lactolac,"Lactococcus lactis")

#Lactolac <- geefunc(dhi$Lactolac,"Lactococcus lactis")
Lactolac <- geefunc("factor(Lactolac)",
        exposure="Lactolac",
        exponame="Lactococcus lactis",
                        expo2=dhi$Lactolac ,
        Other=c("Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","Staspp",
                "Aerococ","Entero","Lactogar","StrepALL",
                "Bac","Coryn"))
Lactolac



Exposure = Streptococcus

Run model and check residuals

gee.plot(dhi$StrepALL,"Streptococcus spp.")

StrepALL <- geefunc("factor(StrepALL)",
        exposure="StrepALL",
        exponame="Streptococcus spp.",
                        expo2=dhi$StrepALL ,
        Other=c("Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","Staspp",
                "Aerococ","Entero","LactoALL",
                "Bac","Coryn"))
#StrepALL <- geefunc(dhi$StrepALL,"Streptococcus spp.")
StrepALL



Exposure = Streptococcus dysgalactiae

Run model and check residuals

gee.plot(dhi$Strepdys,"Streptococcus dysgalactiae")

#Strepdys <- geefunc(dhi$Strepdys,"Streptococcus dysgalactiae")
Strepdys <- geefunc("factor(Strepdys)",
        exposure="Strepdys",
        exponame="Streptococcus dysgalactiae",
                        expo2=dhi$Strepdys ,
        Other=c("Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","Staspp",
                "Aerococ","Entero","LactoALL","Strepub","Strepspp",
                "Bac","Coryn"))
Strepdys



Exposure = Streptococcus uberis

Run model and check residuals

gee.plot(dhi$Strepub,"Streptococcus uberis")

#Strepub <- geefunc(dhi$Strepub,"Streptococcus uberis")
Strepub <- geefunc("factor(Strepub)",
        exposure="Strepub",
        exponame="Streptococcus uberis",
                        expo2=dhi$Strepub ,
        Other=c("Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","Staspp",
                "Aerococ","Entero","LactoALL","Strepdys","Strepspp",
                "Bac","Coryn"))
Strepub



Exposure = Strep species (unable to be determined using MALDI-TOF)

Run model and check residuals

gee.plot(dhi$Strepspp,"Streptococcus sp.")

#Strepspp <- geefunc(dhi$Strepspp,"Streptococcus sp.")
Strepspp <- geefunc("factor(Strepspp)",
        exposure="Strepspp",
        exponame="Streptococcus sp.",
                        expo2=dhi$Strepspp ,
        Other=c("Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","Staspp",
                "Aerococ","Entero","LactoALL","Strepub","Strepdys",
                "Bac","Coryn"))
Strepspp



Exposure = Bacillus

Run model and check residuals

gee.plot(dhi$Bac,"Bacillus sp.")

#Bac <- geefunc(dhi$Bac,"Bacillus sp.")
Bac <- geefunc("factor(Bac)",
        exposure="Bac",
        exponame="Bacillus spp.",
                        expo2=dhi$Bac ,
        Other=c("Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","Staspp",
                "Aerococ","Entero","LactoALL","Strepub","Strepdys",
                "Coryn"))
Bac



Exposure = Corynebacterium

Run model and check residuals

gee.plot(dhi$Coryn,"Corynebacterium sp.")

Coryn <- geefunc("factor(Coryn)",
        exposure="Coryn",
        exponame="Corynebacterium spp.",
                        expo2=dhi$Coryn,
        Other=c("Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","Staspp",
                "Aerococ","Entero","LactoALL","Strepub","Strepdys",
                "Bac"))
#Coryn <- geefunc(dhi$Coryn,"Corynebacterium sp.")
Coryn



Exposure = All IMI

Run model and check residuals

gee.plot(dhi$IMI,"All pathogens")

#IMIALL <- geefunc(dhi$IMI,"All pathogens")
IMIALL <- geefunc("factor(IMI)",
        exposure="IMI",
        exponame="All pathogens",
                        expo2=dhi$IMI ,
        Other=c())
IMIALL



Exposure = Major pathogens

Run model and check residuals

gee.plot(dhi$IMIMAJ,"Major pathogens")

#IMIMAJ <- geefunc(dhi$IMIMAJ,"Major pathogens")
IMIMAJ <- geefunc("factor(IMIMAJ)",
        exposure="IMIMAJ",
        exponame="Major pathogens",
                        expo2=dhi$IMIMAJ ,
        Other=c("IMIMIN"))
IMIMAJ

Exposure = Minor pathogens (NAS, Bacillus and Corynebacterium)

Run model and check residuals

gee.plot(dhi$IMIMIN,"Minor pathogens")

IMIMIN <- geefunc("factor(IMIMIN)",
        exposure="IMIMIN",
        exponame="Minor pathogens",
                        expo2=dhi$IMIMIN ,
        Other=c("IMIMAJ"))
#IMIMIN <- geefunc(dhi$IMIMIN,"Minor pathogens")
IMIMIN

Export subclinical mastitis curves

dhi$TestDIMcat20 <- recode(dhi$TestDIMcat20,`0-20` = "1-20")
gee.plot <- function(expo,exponame){
  mm0 <-  geeglm(SCM ~ factor(expo) + factor(TestDIMcat20) + factor(Dryfreestall) +
                 factor(East) + factor(Parity) + factor(Breed) + factor(Bedtype) + 
                 factor(calvedalone) +  factor(Summer), id = HerdID, family=binomial(link="identity"), data=dhi)
  atx <- unique(dhi$TestDIMcat20)
  emm <- emmeans(mm0, ~expo*TestDIMcat20) %>% tidy()
  emm$LCL <- emm$asymp.LCL
  emm$UCL <- emm$asymp.UCL
  emm$Exposure <- emm$expo %>% as.factor
  emm <- emm %>% subset(select=c(Exposure,TestDIMcat20,estimate,LCL,UCL))
  curve <- ggplot(emm) + coord_cartesian(ylim = (c(0,.35))) + aes(x=TestDIMcat20, y=estimate, group=Exposure, colour=Exposure) +
    labs(y = "Subclinical mastitis risk", x = "Days in milk") + ggtitle(exponame) + geom_point(size=3) + geom_ribbon(aes(ymin=emm$LCL, ymax=emm$UCL,colour=Exposure,fill=Exposure), linetype=0, alpha=0.1) + geom_line(aes(colour=Exposure,linetype=Exposure),size=1.1) + theme(panel.border = element_rect(colour = "black", fill=NA, size=1),axis.text=element_text(size=12,family="Times"),axis.title=element_text(size=14,face="bold",family="Times"),panel.background = element_rect(fill = "white", colour = "white",size = 0.5, linetype = "solid"),panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "grey"), panel.grid.minor = element_line(size = 0.25, linetype = 'solid',colour = "grey"),legend.position = c(0.8, 0.2),legend.text = element_text(colour="black", size=12,face="bold",family="Times")) + scale_color_brewer(palette="Set1")
curve
}

mm0 <-  geeglm(SCM ~ factor(IMIMAJ) + factor(IMIMIN) + factor(TestDIMcat20) + factor(Dryfreestall) +
                 factor(East) + factor(Parity) + factor(Breed) + factor(Bedtype) + 
                 factor(calvedalone) +  factor(Summer), id = HerdID, family=binomial(link="identity"), data=dhi)
atx <- unique(dhi$TestDIMcat20)
emm <- emmeans(mm0, ~IMIMIN*TestDIMcat20) %>% tidy()
emm$LCL <- emm$asymp.LCL
emm$UCL <- emm$asymp.UCL
emm$Exposure <- emm$IMIMIN %>% as.factor
emm <- emm %>% subset(select=c(Exposure,TestDIMcat20,estimate,LCL,UCL))

MIN <- ggplot(emm) + coord_cartesian(ylim = (c(0,.35))) + aes(x=TestDIMcat20, y=estimate, group=Exposure, colour=Exposure) +
    labs(y = "Subclinical mastitis risk", x = "Days in milk") + ggtitle("Minor pathogens") + geom_point(size=3) + geom_ribbon(aes(ymin=emm$LCL, ymax=emm$UCL,colour=Exposure,fill=Exposure), linetype=0, alpha=0.1) + geom_line(aes(colour=Exposure,linetype=Exposure),size=1.1) + theme(panel.border = element_rect(colour = "black", fill=NA, size=1),axis.text=element_text(size=12,family="Times"),axis.title=element_text(size=14,face="bold",family="Times"),panel.background = element_rect(fill = "white", colour = "white",size = 0.5, linetype = "solid"),panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "grey"), panel.grid.minor = element_line(size = 0.25, linetype = 'solid',colour = "grey"),legend.position = c(0.8, 0.2),legend.text = element_text(colour="black", size=12,face="bold",family="Times")) + scale_color_brewer(palette="Set1")


emm <- emmeans(mm0, ~IMIMAJ*TestDIMcat20) %>% tidy()
emm$LCL <- emm$asymp.LCL
emm$UCL <- emm$asymp.UCL
emm$Exposure <- emm$IMIMAJ %>% as.factor
emm <- emm %>% subset(select=c(Exposure,TestDIMcat20,estimate,LCL,UCL))
MAJ <- ggplot(emm) + coord_cartesian(ylim = (c(0,.35))) + aes(x=TestDIMcat20, y=estimate, group=Exposure, colour=Exposure) +
    labs(y = "Subclinical mastitis risk", x = "Days in milk") + ggtitle("Major pathogens") + geom_point(size=3) + geom_ribbon(aes(ymin=emm$LCL, ymax=emm$UCL,colour=Exposure,fill=Exposure), linetype=0, alpha=0.1) + geom_line(aes(colour=Exposure,linetype=Exposure),size=1.1) + theme(panel.border = element_rect(colour = "black", fill=NA, size=1),axis.text=element_text(size=12,family="Times"),axis.title=element_text(size=14,face="bold",family="Times"),panel.background = element_rect(fill = "white", colour = "white",size = 0.5, linetype = "solid"),panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "grey"), panel.grid.minor = element_line(size = 0.25, linetype = 'solid',colour = "grey"),legend.position = c(0.8, 0.2),legend.text = element_text(colour="black", size=12,face="bold",family="Times")) + scale_color_brewer(palette="Set1")



MIN <- gee.plot(dhi$IMIMAJ,"Major IMI exposure")
MAJ <- gee.plot(dhi$IMIMIN,"Minor IMI exposure")

tiff("DIM1.tiff", units="in", width=4, height=9, res=300)
cowplot::plot_grid(MIN,MAJ,nrow=2)
dev.off()
## quartz_off_screen 
##                 2

Outcome 4: Final models

Blank <- Staaur
Blank
Blank$Exposure <- NA
Blank$RiskExp0 <- NA
Blank$RiskExp1 <- NA
Blank$RR <- NA
Blank$estimate <- NA
Blank$conf.low <- NA
Blank$conf.high <- NA


SSLO$Exposure <- "Strep and Strep-like organisms"
OGP <- Blank
OGP$Exposure <- "Other Gram-positives"

Table <- rbind(NAS,Stachr,Staepi,Stahaem,Stasci,Stasim,Staxyl,Staspp,Blank,
               SSLO,Aerococ,Entero,LactoALL,Lactogar,Lactolac,
               StrepALL,Strepdys,Strepub,Strepspp,Blank,OGP,
               Staaur,Bac,Coryn,Blank,IMIALL,IMIMAJ,IMIMIN)

FP <- Table
FP$estimate <- FP$estimate
FP$conf.low <- FP$conf.low
FP$conf.high <- FP$conf.high

## Labels defining subgroups that are indented
subgps <- c(2:8,11:19,22:24)
FP$Exposure[subgps] <- paste("     ",FP$Exposure[subgps]) 
subgps <- c(14:15,17:19)
FP$Exposure[subgps] <- paste("  ",FP$Exposure[subgps]) 

## The rest of the columns in the table. 
tabletext <- cbind(c("Exposure","\n",FP$Exposure), 
                   c("exposed","\n",FP$RiskExp1),
                   c("unexposed","\n",FP$RiskExp0),
                   c("Risk Ratio (95% CI)","\n",FP$RR))

library(forestplot)
tiff("FP4.tiff", units="in", width=8, height=6, res=300)
forestplot(labeltext=tabletext, graph.pos=4, 
                  graphwidth = unit(35,'mm'), 
           clip=c(0.2,5),
           xticks = c(.2,1,5),
           xlog= T,
           mean=c(NA,NA,FP$estimate), 
           lower=c(NA,NA,FP$conf.low), upper=c(NA,NA,FP$conf.high),
           title="",
                  xlab="Risk ratio for subclinical mastitis",
                  #hrzl_lines=list("8" = gpar(lwd=100, lineend="butt", columns=c(1:5), col="#99999922")),
                  #               ,"14" = gpar(lwd=110, lineend="butt", columns=c(1:5), col="#99999922")),
           txt_gp=fpTxtGp(label=gpar(fontfamily="Times", cex=.8),
                          ticks=gpar(fontfamily="Times", cex=.8),
                          xlab=gpar(fontfamily="Times", cex=.8),
                          title=gpar(fontfamily="Times")),
           col=fpColors(box="black", lines="black", zero = "gray50"),
           zero=1, cex=0.9, lineheight = "auto", boxsize=0.25, colgap=unit(6,"mm"),
           lwd.ci=2, ci.vertices=TRUE, ci.vertices.height = 0.2)

dev.off()
## quartz_off_screen 
##                 2

Exposure table

#Tables
countexp <- function(Expo,Exponame){
a <- NA %>% as.data.frame
a$Exposure <- Exponame
a$Count <- table(cow$HerdID,Expo) %>% as.data.frame %>% filter(Expo==1) %>% select(Freq) %>% filter(Freq>0) %>% nrow
a$Herd <- a$Count
b <- a %>% select(Exposure,Herd)
a$Exposure <- Exponame
a$Count <- table(cow$HerdID,Expo) %>% as.data.frame %>% filter(Expo==1) %>% select(Freq) %>% sum
a$Cow <- a$Count
a <- a %>% select(Exposure,Cow)
merge(a,b,by=c("Exposure")) }

countexp(cow$Lactolac,"Lactococcus lactis")
Staaur <- countexp(cow$Staaur,"Staphylococcus aureus")
NAS <- countexp(cow$NAS,"Non aureus Staphylococcus")
Stachr <- countexp(cow$Stachr,"Staphylococcus chromogenes")
Staepi <- countexp(cow$Staepi,"Staphylococcus epidermidis")
Stahaem <- countexp(cow$Stahaem,"Staphylococcus haemolyticus")
Stasci <- countexp(cow$Stasci,"Staphylococcus sciuri")
Stasim <- countexp(cow$Stasim,"Staphylococcus simulans")
Staxyl <- countexp(cow$Staxyl,"Staphylococcus xylosus")
Staspp <- countexp(cow$Staspp,"Staphylococcus sp.")
SSLO <- countexp(cow$SSLO,"Strep and Strep-like organisms.")
Aerococ <- countexp(cow$Aerococ,"Aerococcus spp.")
Entero <- countexp(cow$Entero,"Enterococcus spp.")
LactoALL <- countexp(cow$LactoALL,"Lactococcus spp.")
Lactolac <- countexp(cow$Lactolac,"Lactococcus lactis")
Lactogar <- countexp(cow$Lactogar,"Lactococcus garviae")
StrepALL <- countexp(cow$StrepALL,"Streptococcus spp.")
Strepdys <- countexp(cow$Strepdys,"Streptococcus dysgalactiae")
Strepub <- countexp(cow$Strepub,"Streptococcus uberis")
Strepspp <- countexp(cow$Strepspp,"Streptococcus sp.")
Bac <- countexp(cow$Bac,"Bacillus spp.")
Coryn <- countexp(cow$Coryn,"Corynebacterium spp.")
IMI <- countexp(cow$IMI,"All pathogens")
IMIMAJ <- countexp(cow$IMIMAJ,"Major pathogens")
IMIMIN <- countexp(cow$IMIMIN,"Minor pathogens")

coli <- countexp(cow$Coliform,"Coliform")
Pseudo <- countexp(cow$Pseudo,"Pseudo")
miscgn <- countexp(cow$MiscGN,"Misc GN")
miscgp <- countexp(cow$MiscGP,"Misc GP")

table <- rbind(NAS,Stachr,Staepi,Stahaem,Stasci,Stasim,Staxyl,Staspp,
     SSLO,Aerococ,Entero,LactoALL,Lactogar,Lactolac,
     StrepALL,Strepdys,Strepub,Strepspp,Staaur,
     Bac,Coryn,IMI,IMIMAJ,IMIMIN,coli,Pseudo,miscgn,miscgp)

cowsatrisk <- table(cow$IMI) %>% as.data.frame %>% select(Freq) %>% sum
herdsatrisk <- table(cow$HerdID,cow$IMI) %>% as.data.frame %>% filter(Var2==0) %>% select(Freq) %>% nrow

table$Herdtotal <- herdsatrisk
table$Cowstotal <- cowsatrisk
table$CowPer <- table$Cow/table$Cowstotal
table$HerdPer <- table$Herd/table$Herdtotal
table$CowExp <- paste0(table$Cow," (",sprintf("%.1f",round(100*table$CowPer, 1)),")")
table$HerdExp <- paste0(table$Herd," (",sprintf("%.1f",round(100*table$HerdPer,1)),")")
table <- table %>% select(Exposure,CowExp,HerdExp)
table
write.csv(table,'exposure.csv')

Herd table

load(file="Herds.Rdata")
a <- cow %>% select(HerdID,calvedalone)
b <- merge(Herd,a,by="HerdID") %>% filter(Season==1)
c <- b %>% distinct(HerdID, .keep_all= TRUE) %>% select(State,Breed,ProdL,Herdsize,Lacthouse,DryHouse,Bedtype,MilkSchedual,Predip,Postdip,Forestripping,TS,Housedrysum,Housedrywin,DryPaddock,DCT,ITS,ETS,Vax,Calvingarea,calvedalone)

c$ProdL <- c$ProdL*0.453592
table1(~factor(State) + factor(Breed) + ProdL + Herdsize + factor(Lacthouse) + factor(DryHouse) +  factor(Bedtype) + factor(MilkSchedual)  + factor(Predip)  + factor(Postdip)  + factor(Forestripping) + factor(TS) + factor(DryPaddock) + factor(DCT) + factor(ITS) + factor(ETS) +factor(Vax) + factor(Calvingarea) + factor(calvedalone), data=c)
Overall
(N=74)
factor(State)
California 16 (21.6%)
Idaho 5 (6.8%)
indiana 3 (4.1%)
Michigan 5 (6.8%)
minnesota 10 (13.5%)
New York 9 (12.2%)
oregon 1 (1.4%)
texas 2 (2.7%)
Washington 4 (5.4%)
Wisconsin 19 (25.7%)
factor(Breed)
h 64 (86.5%)
j 3 (4.1%)
MIX 7 (9.5%)
ProdL
Mean (SD) 37.6 (4.81)
Median [Min, Max] 38.6 [22.7, 47.2]
Herdsize
Mean (SD) 2330 (1890)
Median [Min, Max] 1820 [235, 9650]
factor(Lacthouse)
Dry lot 2 (2.7%)
Free stall 71 (95.9%)
Mixed 1 (1.4%)
factor(DryHouse)
Free stall 50 (67.6%)
Loose - bedded pack 2 (2.7%)
Loose - dry lot 18 (24.3%)
Mixed 4 (5.4%)
factor(Bedtype)
0 22 (29.7%)
1 19 (25.7%)
2 15 (20.3%)
3 18 (24.3%)
factor(MilkSchedual)
2x 21 (28.4%)
3x 53 (71.6%)
factor(Predip)
n 4 (5.4%)
y 70 (94.6%)
factor(Postdip)
y 74 (100%)
factor(Forestripping)
n 13 (17.6%)
y 61 (82.4%)
factor(TS)
n 63 (85.1%)
y 11 (14.9%)
factor(DryPaddock)
No 70 (94.6%)
Yes - summer 3 (4.1%)
Yes - summer and winter 1 (1.4%)
factor(DCT)
Never 1 (1.4%)
Some cows --> 7 (9.5%)
Yes, for all cows---> 66 (89.2%)
factor(ITS)
0 1 (1.4%)
n 21 (28.4%)
y 52 (70.3%)
factor(ETS)
0 1 (1.4%)
n 67 (90.5%)
y 6 (8.1%)
factor(Vax)
NO 8 (10.8%)
YES 66 (89.2%)
factor(Calvingarea)
In a dedicated calving area (separated from all other cows) 27 (36.5%)
In a dedicated calving area (with other cows that are close to calving) 44 (59.5%)
With the other dry cows in the dry cow area 3 (4.1%)
factor(calvedalone)
0 47 (63.5%)
1 27 (36.5%)