Post-calving udder health in cows approaching dry-off with intramammary infections caused by non-aureus Staphylococcus, Aerococcus, Enterococcus, Lactococcus and Streptococcus species
Sam Rowe (samrowe101@gmail.com)
April 2020
- Evaluate Data
- DAG
- Check for correlation between covariates
- Outcome 1: Clinical mastitis
- 120 d Clinical mastitis risk
- Kaplan Meier curve & log rank test for clinical mastitis
- Cox proportional hazards regression for clinical mastitis (1-120 DIM)
- Model building plan
- Check proprtional hazards assumption for covariates
- Create functions to simplify the modelling syntax
- Exposure = Staphylococcus aureus
- Exposure = Non aureus Staphylococcus sp. (whole group)
- Exposure = Staphylococcus chromogenes
- Exposure = Staphylococcus epidermidis
- Exposure = Staphylococcus haemolyticus
- Exposure = Staphylococcus sciuri
- Exposure = Staphylococcus simulans
- Exposure = Staphylococcus xylosus
- Exposure = Staphylococcus species (i.e. species unable to be determined using MALDI-TOF)
- Exposure = Aerococcus species
- Exposure = Enterococcus species
- Exposure = Lactococcus (all species combined)
- Exposure = Lactococcus garviae
- Exposure = Lactococcus lactis
- Exposure = Streptococcus
- Exposure = Streptococcus dysgalactiae
- Exposure = Streptococcus uberis
- Exposure = Strep species (unable to be determined using MALDI-TOF)
- Exposure = Bacillus
- Exposure = Corynebacterium
- Exposure = All IMI
- Exposure = Major pathogens
- Exposure = Minor pathogens (NAS, Bacillus and Corynebacterium)
- Outcome 1: Final models
- Outcome 2: Culling risk 0-120 DIM
- Kaplan Meier curve & log rank test for culling risk
- Cox proportional hazards regression for culling or death (1-120 DIM)
- Model building plan
- Check proprtional hazards assumption for covariates
- Create functions to simplify the modelling syntax
- Exposure = Staphylococcus aureus
- Exposure = Non aureus Staphylococcus sp. (whole group)
- Exposure = Staphylococcus chromogenes
- Exposure = Staphylococcus epidermidis
- Exposure = Staphylococcus haemolyticus
- Exposure = Staphylococcus sciuri
- Exposure = Staphylococcus simulans
- Exposure = Staphylococcus xylosus
- Exposure = Staphylococcus species (i.e. species unable to be determined using MALDI-TOF)
- Exposure = Aerococcus species
- Exposure = Enterococcus species
- Exposure = Lactococcus (all species combined)
- Exposure = Lactococcus garviae
- Exposure = Lactococcus lactis
- Exposure = Streptococcus
- Exposure = Streptococcus dysgalactiae
- Exposure = Streptococcus uberis
- Exposure = Strep species (unable to be determined using MALDI-TOF)
- Exposure = Bacillus
- Exposure = Corynebacterium
- Exposure = All IMI
- Exposure = Major pathogens
- Exposure = Minor pathogens (NAS, Bacillus and Corynebacterium)
- Outcome 2 Final models
PART 1: Initial analyses and clinical mastitis / removal from the herd
library(knitr)
library(readr)
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"))Evaluate Data
Show first 10 rows of data
library(dplyr)
Desc1 <- postcalve %>% subset(select=c(Herdcow,HerdID,Parity,DPlength,CM2,CM2TAR,Rem2,Rem2TAR))
head(Desc1, n=10)
Inspect Data
#summarytools::dfSummary(Baseline, style='grid')
print(summarytools::dfSummary(postcalve,valid.col=FALSE, graph.col=F, style="grid"))## Data Frame Summary
## postcalve
## Dimensions: 2581 x 95
## Duplicates: 0
##
## +----+---------------+---------------------------------+----------------------+----------+
## | No | Variable | Stats / Values | Freqs (% of Valid) | Missing |
## +====+===============+=================================+======================+==========+
## | 1 | HerdID | Mean (sd) : 40.8 (23.3) | 74 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 1 < 41 < 80 | | |
## | | | IQR (CV) : 41 (0.6) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 2 | Cow | Mean (sd) : 12281.4 (13316.5) | 2440 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 1 < 7967 < 82342 | | |
## | | | IQR (CV) : 9993 (1.1) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 3 | ABX | Min : 0 | 0 : 185 ( 7.2%) | 0 |
## | | [numeric] | Mean : 0.9 | 1 : 2396 (92.8%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 4 | Rem2 | Min : 0 | 0 : 2271 (88.0%) | 0 |
## | | [numeric] | Mean : 0.1 | 1 : 310 (12.0%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 5 | Rem2comment | 1. ··············· | 77 (11.4%) | 1908 |
## | | [character] | 2. CAR3··········· | 54 ( 8.0%) | (73.92%) |
## | | | 3. CAR7··········· | 31 ( 4.6%) | |
## | | | 4. CAR1··········· | 30 ( 4.5%) | |
## | | | 5. CAR5··········· | 27 ( 4.0%) | |
## | | | 6. MAST··········· | 25 ( 3.7%) | |
## | | | 7. CAR2··········· | 16 ( 2.4%) | |
## | | | 8. DCAR3·········· | 15 ( 2.2%) | |
## | | | 9. PROD··········· | 15 ( 2.2%) | |
## | | | 10. LAME··········· | 13 ( 1.9%) | |
## | | | [ 196 others ] | 370 (55.0%) | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 6 | RemDate | min : 2017-09-05 | 310 distinct values | 1908 |
## | | [Date] | med : 2018-06-12 | | (73.92%) |
## | | | max : 2019-05-15 | | |
## | | | range : 1y 8m 10d | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 7 | RemLate | 1 distinct value | 0 : 2581 (100.0%) | 0 |
## | | [numeric] | | | (0%) |
## +----+---------------+---------------------------------+----------------------+----------+
## | 8 | CM2 | Min : 0 | 0 : 2080 (86.1%) | 165 |
## | | [numeric] | Mean : 0.1 | 1 : 336 (13.9%) | (6.39%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 9 | CMDate | min : 2017-10-12 | 200 distinct values | 2242 |
## | | [Date] | med : 2018-03-14 | | (86.87%) |
## | | | max : 2018-10-28 | | |
## | | | range : 1y 0m 16d | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 10 | CMLate | Min : 0 | 0 : 2395 (99.1%) | 165 |
## | | [numeric] | Mean : 0 | 1 : 21 ( 0.9%) | (6.39%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 11 | Calv2Date | min : 2017-08-25 | 248 distinct values | 0 |
## | | [Date] | med : 2017-12-13 | | (0%) |
## | | | max : 2018-07-15 | | |
## | | | range : 10m 20d | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 12 | Herdcow | Mean (sd) : 4432.2 (3378.3) | 2580 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 38 < 3838.2 < 18014.3 | | |
## | | | IQR (CV) : 4136 (0.8) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 13 | PeakSCC42D | Mean (sd) : 4.4 (1.5) | 486 distinct values | 439 |
## | | [numeric] | min < med < max: | | (17.01%) |
## | | | 0 < 4.4 < 9.2 | | |
## | | | IQR (CV) : 1.6 (0.3) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 14 | PeakSCC3D | Mean (sd) : 4.8 (1.4) | 644 distinct values | 246 |
## | | [numeric] | min < med < max: | | (9.53%) |
## | | | 0 < 4.8 < 9.2 | | |
## | | | IQR (CV) : 1.6 (0.3) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 15 | PeakSCCD | Mean (sd) : 5.4 (1.5) | 858 distinct values | 246 |
## | | [numeric] | min < med < max: | | (9.53%) |
## | | | 0 < 5.3 < 9.2 | | |
## | | | IQR (CV) : 1.9 (0.3) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 16 | SCCDDO | Mean (sd) : 4.4 (1.5) | 516 distinct values | 246 |
## | | [numeric] | min < med < max: | | (9.53%) |
## | | | 0 < 4.4 < 9.2 | | |
## | | | IQR (CV) : 1.7 (0.3) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 17 | MYDDO | Mean (sd) : 60.4 (22.5) | 123 distinct values | 246 |
## | | [numeric] | min < med < max: | | (9.53%) |
## | | | 0 < 61 < 179 | | |
## | | | IQR (CV) : 29 (0.4) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 18 | MYMEDDO | 1. 30520 | 9 ( 0.4%) | 246 |
## | | [character] | 2. 29270 | 7 ( 0.3%) | (9.53%) |
## | | | 3. 25850 | 6 ( 0.3%) | |
## | | | 4. 29760 | 6 ( 0.3%) | |
## | | | 5. 30970 | 6 ( 0.3%) | |
## | | | 6. 25570 | 5 ( 0.2%) | |
## | | | 7. 27690 | 5 ( 0.2%) | |
## | | | 8. 27990 | 5 ( 0.2%) | |
## | | | 9. 28040 | 5 ( 0.2%) | |
## | | | 10. 28790 | 5 ( 0.2%) | |
## | | | [ 1402 others ] | 2276 (97.5%) | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 19 | CM1Dcount | Mean (sd) : 0.3 (0.6) | 0 : 1966 (82.4%) | 194 |
## | | [numeric] | min < med < max: | 1 : 301 (12.6%) | (7.52%) |
## | | | 0 < 0 < 7 | 2 : 84 ( 3.5%) | |
## | | | IQR (CV) : 0 (2.6) | 3 : 23 ( 1.0%) | |
## | | | | 4 : 8 ( 0.3%) | |
## | | | | 5 : 3 ( 0.1%) | |
## | | | | 6 : 1 ( 0.0%) | |
## | | | | 7 : 1 ( 0.0%) | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 20 | PeakSCC42 | Mean (sd) : 4.3 (1.4) | 461 distinct values | 455 |
## | | [numeric] | min < med < max: | | (17.63%) |
## | | | 0 < 4.3 < 8.9 | | |
## | | | IQR (CV) : 1.6 (0.3) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 21 | PeakSCC3 | Mean (sd) : 4.7 (1.4) | 621 distinct values | 238 |
## | | [numeric] | min < med < max: | | (9.22%) |
## | | | 0 < 4.7 < 9.2 | | |
## | | | IQR (CV) : 1.6 (0.3) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 22 | PeakSCC | Mean (sd) : 5.4 (1.5) | 842 distinct values | 238 |
## | | [numeric] | min < med < max: | | (9.22%) |
## | | | 0 < 5.2 < 9.2 | | |
## | | | IQR (CV) : 1.9 (0.3) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 23 | SCCDO | Mean (sd) : 4.3 (1.4) | 491 distinct values | 238 |
## | | [numeric] | min < med < max: | | (9.22%) |
## | | | 0 < 4.3 < 9 | | |
## | | | IQR (CV) : 1.6 (0.3) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 24 | MYDO | Mean (sd) : 64.8 (22.9) | 123 distinct values | 238 |
## | | [numeric] | min < med < max: | | (9.22%) |
## | | | 0 < 65 < 179 | | |
## | | | IQR (CV) : 29 (0.4) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 25 | MYMEDO | Mean (sd) : 13515.7 (2616.1) | 1422 distinct values | 238 |
## | | [numeric] | min < med < max: | | (9.22%) |
## | | | 4018.8 < 13580.5 < 23863.5 | | |
## | | | IQR (CV) : 3517.6 (0.2) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 26 | CM1count | Mean (sd) : 0.2 (0.6) | 0 : 2001 (82.8%) | 165 |
## | | [numeric] | min < med < max: | 1 : 299 (12.4%) | (6.39%) |
## | | | 0 < 0 < 7 | 2 : 83 ( 3.4%) | |
## | | | IQR (CV) : 0 (2.6) | 3 : 21 ( 0.9%) | |
## | | | | 4 : 8 ( 0.3%) | |
## | | | | 5 : 3 ( 0.1%) | |
## | | | | 7 : 1 ( 0.0%) | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 27 | DryDate | min : 2017-08-01 | 153 distinct values | 30 |
## | | [Date] | med : 2017-10-23 | | (1.16%) |
## | | | max : 2018-05-02 | | |
## | | | range : 9m 1d | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 28 | Calv1Date | min : 2015-08-17 | 446 distinct values | 38 |
## | | [Date] | med : 2016-12-10 | | (1.47%) |
## | | | max : 2017-08-09 | | |
## | | | range : 1y 11m 23d | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 29 | Parity | 1. 1 | 1119 (43.4%) | 0 |
## | | [character] | 2. 2 | 826 (32.0%) | (0%) |
## | | | 3. 3 | 390 (15.1%) | |
## | | | 4. 4 | 246 ( 9.5%) | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 30 | Breed | Min : 0 | 0 : 2217 (85.9%) | 0 |
## | | [numeric] | Mean : 0.1 | 1 : 364 (14.1%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 31 | Summer | Min : 0 | 0 : 1267 (49.1%) | 0 |
## | | [numeric] | Mean : 0.5 | 1 : 1314 (50.9%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 32 | Path1 | 1. Staphyloccocus chromogene | 184 (33.5%) | 2032 |
## | | [character] | 2. Bacillus sp. | 49 ( 8.9%) | (78.73%) |
## | | | 3. Staphylococcus sp. | 46 ( 8.4%) | |
## | | | 4. Aerococcus viridans | 40 ( 7.3%) | |
## | | | 5. Lactococcus garvieae | 33 ( 6.0%) | |
## | | | 6. Corynebacterium sp. | 22 ( 4.0%) | |
## | | | 7. Staphylococcus simulans | 19 ( 3.5%) | |
## | | | 8. Staphylococcus xylosus/sa | 14 ( 2.6%) | |
## | | | 9. Gram Positive Cocci | 13 ( 2.4%) | |
## | | | 10. Staphylococcus epidermidi | 12 ( 2.2%) | |
## | | | [ 22 others ] | 117 (21.3%) | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 33 | Path2 | 1. Staphyloccocus chromogene | 8 (20.0%) | 2541 |
## | | [character] | 2. Corynebacterium sp. | 6 (15.0%) | (98.45%) |
## | | | 3. Lactococcus garvieae | 4 (10.0%) | |
## | | | 4. Pseudomonas sp. | 3 ( 7.5%) | |
## | | | 5. Enterococcus sp. | 2 ( 5.0%) | |
## | | | 6. Gram Positive Cocci | 2 ( 5.0%) | |
## | | | 7. Staphylococcus xylosus/sa | 2 ( 5.0%) | |
## | | | 8. Aerococcus sp. | 1 ( 2.5%) | |
## | | | 9. Aerococcus viridans | 1 ( 2.5%) | |
## | | | 10. Bacillus sp. | 1 ( 2.5%) | |
## | | | [ 10 others ] | 10 (25.0%) | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 34 | Bedtype | Mean (sd) : 1.4 (1.2) | 0 : 792 (30.7%) | 0 |
## | | [numeric] | min < med < max: | 1 : 581 (22.5%) | (0%) |
## | | | 0 < 1 < 3 | 2 : 520 (20.2%) | |
## | | | IQR (CV) : 3 (0.8) | 3 : 688 (26.7%) | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 35 | East | Min : 0 | 0 : 1031 (40.0%) | 0 |
## | | [numeric] | Mean : 0.6 | 1 : 1550 (60.1%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 36 | Herdsize | Mean (sd) : 2366 (1909.8) | 64 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 235 < 1840 < 9650 | | |
## | | | IQR (CV) : 1550 (0.8) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 37 | Dryfreestall | Min : 0 | 0 : 928 (36.0%) | 0 |
## | | [numeric] | Mean : 0.6 | 1 : 1653 (64.0%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 38 | ITS | Min : 0 | 0 : 742 (29.2%) | 40 |
## | | [numeric] | Mean : 0.7 | 1 : 1799 (70.8%) | (1.55%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 39 | Vax | Min : 0 | 0 : 268 (10.4%) | 0 |
## | | [numeric] | Mean : 0.9 | 1 : 2313 (89.6%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 40 | calvedalone | Min : 0 | 0 : 1669 (64.7%) | 0 |
## | | [numeric] | Mean : 0.4 | 1 : 912 (35.3%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 41 | Aerococ | Min : 0 | 0 : 2416 (93.6%) | 0 |
## | | [numeric] | Mean : 0.1 | 1 : 165 ( 6.4%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 42 | Bac | Min : 0 | 0 : 2429 (94.1%) | 0 |
## | | [numeric] | Mean : 0.1 | 1 : 152 ( 5.9%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 43 | Coryn | Min : 0 | 0 : 2480 (96.1%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 101 ( 3.9%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 44 | Entero | Min : 0 | 0 : 2529 (98.0%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 52 ( 2.0%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 45 | LactoALL | Min : 0 | 0 : 2441 (94.6%) | 0 |
## | | [numeric] | Mean : 0.1 | 1 : 140 ( 5.4%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 46 | Lactogar | Min : 0 | 0 : 2487 (96.4%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 94 ( 3.6%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 47 | Lactolac | Min : 0 | 0 : 2538 (98.3%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 43 ( 1.7%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 48 | Lactosp | Min : 0 | 0 : 2571 (99.6%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 10 ( 0.4%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 49 | Pseudo | Min : 0 | 0 : 2566 (99.4%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 15 ( 0.6%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 50 | Staaur | Min : 0 | 0 : 2550 (98.8%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 31 ( 1.2%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 51 | NAS | Min : 0 | 0 : 1818 (70.4%) | 0 |
## | | [numeric] | Mean : 0.3 | 1 : 763 (29.6%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 52 | Stachr | Min : 0 | 0 : 2087 (80.9%) | 0 |
## | | [numeric] | Mean : 0.2 | 1 : 494 (19.1%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 53 | Staepi | Min : 0 | 0 : 2537 (98.3%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 44 ( 1.7%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 54 | Stahaem | Min : 0 | 0 : 2545 (98.6%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 36 ( 1.4%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 55 | Stasci | Min : 0 | 0 : 2548 (98.7%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 33 ( 1.3%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 56 | Stasim | Min : 0 | 0 : 2536 (98.3%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 45 ( 1.7%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 57 | Staxyl | Min : 0 | 0 : 2537 (98.3%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 44 ( 1.7%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 58 | Staspp | Min : 0 | 0 : 2404 (93.1%) | 0 |
## | | [numeric] | Mean : 0.1 | 1 : 177 ( 6.9%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 59 | StrepALL | Min : 0 | 0 : 2486 (96.3%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 95 ( 3.7%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 60 | Strepdys | Min : 0 | 0 : 2554 (99.0%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 27 ( 1.1%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 61 | Strepub | Min : 0 | 0 : 2553 (98.9%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 28 ( 1.1%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 62 | Strepspp | Min : 0 | 0 : 2538 (98.3%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 43 ( 1.7%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 63 | Yeast | Min : 0 | 0 : 2566 (99.4%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 15 ( 0.6%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 64 | Coliform | Min : 0 | 0 : 2567 (99.5%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 14 ( 0.5%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 65 | MiscGP | Min : 0 | 0 : 2509 (97.2%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 72 ( 2.8%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 66 | MiscGN | Min : 0 | 0 : 2541 (98.5%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 40 ( 1.6%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 67 | Enrolldate | min : 2017-08-01 | 60 distinct values | 0 |
## | | [Date] | med : 2017-09-13 | | (0%) |
## | | | max : 2018-03-27 | | |
## | | | range : 7m 26d | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 68 | IMI | Min : 0 | 0 : 1301 (50.4%) | 0 |
## | | [numeric] | Mean : 0.5 | 1 : 1280 (49.6%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 69 | IMIMAJ | Min : 0 | 0 : 2114 (81.9%) | 0 |
## | | [numeric] | Mean : 0.2 | 1 : 467 (18.1%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 70 | IMIMAJ2 | Min : 0 | 0 : 2019 (78.2%) | 0 |
## | | [numeric] | Mean : 0.2 | 1 : 562 (21.8%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 71 | IMIMIN | Min : 0 | 0 : 1647 (63.8%) | 0 |
## | | [numeric] | Mean : 0.4 | 1 : 934 (36.2%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 72 | DIMDO | Mean (sd) : 333.9 (60.2) | 292 distinct values | 56 |
## | | [numeric] | min < med < max: | | (2.17%) |
## | | | 212 < 314 < 823 | | |
## | | | IQR (CV) : 60 (0.2) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 73 | DIMEN | Mean (sd) : 321 (65.5) | 302 distinct values | 38 |
## | | [numeric] | min < med < max: | | (1.47%) |
## | | | 212 < 302 < 887 | | |
## | | | IQR (CV) : 63 (0.2) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 74 | RecCM3HT | Min : 0 | 0 : 2320 (96.0%) | 165 |
## | | [numeric] | Mean : 0 | 1 : 96 ( 4.0%) | (6.39%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 75 | RecCM14 | Min : 0 | 0 : 2409 (99.7%) | 165 |
## | | [numeric] | Mean : 0 | 1 : 7 ( 0.3%) | (6.39%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 76 | RecCM3HTD | Min : 0 | 0 : 2318 (95.9%) | 165 |
## | | [numeric] | Mean : 0 | 1 : 98 ( 4.1%) | (6.39%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 77 | RecCM14D | Min : 0 | 0 : 2405 (99.5%) | 165 |
## | | [numeric] | Mean : 0 | 1 : 11 ( 0.5%) | (6.39%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 78 | En2Calv | min : 3 | 110 distinct values | 0 |
## | | [difftime] | med : 70 | | (0%) |
## | | | max : 129 | | |
## | | | units : days | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 79 | DPlength | min : 3 | 109 distinct values | 30 |
## | | [difftime] | med : 56 | | (1.16%) |
## | | | max : 129 | | |
## | | | units : days | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 80 | CM1 | Min : 0 | 0 : 2001 (82.8%) | 165 |
## | | [numeric] | Mean : 0.2 | 1 : 415 (17.2%) | (6.39%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 81 | CM1D | Min : 0 | 0 : 1966 (82.4%) | 194 |
## | | [numeric] | Mean : 0.2 | 1 : 421 (17.6%) | (7.52%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 82 | RCens | min : 2017-12-23 | 248 distinct values | 0 |
## | | [Date] | med : 2018-04-12 | | (0%) |
## | | | max : 2018-11-12 | | |
## | | | range : 10m 20d | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 83 | RCens60 | min : 2017-10-24 | 248 distinct values | 0 |
## | | [Date] | med : 2018-02-11 | | (0%) |
## | | | max : 2018-09-13 | | |
## | | | range : 10m 20d | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 84 | CMCens | min : 2017-09-05 | 359 distinct values | 0 |
## | | [Date] | med : 2018-04-07 | | (0%) |
## | | | max : 2018-11-08 | | |
## | | | range : 1y 2m 3d | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 85 | CM2Cens60 | min : 2017-09-05 | 284 distinct values | 0 |
## | | [Date] | med : 2018-02-10 | | (0%) |
## | | | max : 2018-09-13 | | |
## | | | range : 1y 0m 8d | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 86 | CM2TAR | Mean (sd) : 103.2 (34.5) | 120 distinct values | 165 |
## | | [numeric] | min < med < max: | | (6.39%) |
## | | | 0 < 120 < 120 | | |
## | | | IQR (CV) : 0 (0.3) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 87 | CM2TAR60 | Mean (sd) : 55 (13.9) | 60 distinct values | 165 |
## | | [numeric] | min < med < max: | | (6.39%) |
## | | | 0 < 60 < 60 | | |
## | | | IQR (CV) : 0 (0.3) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 88 | CM260 | Min : 0 | 0 : 2225 (92.1%) | 165 |
## | | [numeric] | Mean : 0.1 | 1 : 191 ( 7.9%) | (6.39%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 89 | Rem2Cens | min : 2017-09-05 | 321 distinct values | 0 |
## | | [Date] | med : 2018-04-09 | | (0%) |
## | | | max : 2018-11-12 | | |
## | | | range : 1y 2m 7d | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 90 | RemTAR | Mean (sd) : 182.2 (32.1) | 191 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 18 < 187 < 249 | | |
## | | | IQR (CV) : 26 (0.2) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 91 | Rem2TAR | Mean (sd) : 111.1 (27.2) | 110 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 0 < 120 < 120 | | |
## | | | IQR (CV) : 0 (0.2) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 92 | RemALLRCens | min : 2017-09-05 | 254 distinct values | 0 |
## | | [Date] | med : 2018-05-16 | | (0%) |
## | | | max : 2018-12-02 | | |
## | | | range : 1y 2m 27d | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 93 | RemALLTAR | Mean (sd) : 232.8 (45.8) | 176 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 18 < 250 < 250 | | |
## | | | IQR (CV) : 0 (0.2) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 94 | DOSCM | 1. 0 | 1850 (79.0%) | 238 |
## | | [factor] | 2. 1 | 493 (21.0%) | (9.22%) |
## +----+---------------+---------------------------------+----------------------+----------+
## | 95 | SSLO | Min : 0 | 0 : 2171 (84.1%) | 0 |
## | | [numeric] | Mean : 0.2 | 1 : 410 (15.9%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
#Describe subjects in study
colnames(cow)## [1] "Cow" "RemALL" "RemALLDate" "Rem2" "Rem2comment"
## [6] "RemDate" "RemDP" "RemLate" "CMALL" "CMALLDate"
## [11] "CM2" "CMDate" "CMEnCa" "CMDry" "CMLate"
## [16] "Calv2Date" "Herdcow" "PeakSCC42D" "PeakSCC3D" "PeakSCCD"
## [21] "SCCDDO" "MYDDO" "MYMEDDO" "CM1Dcount" "PeakSCC42"
## [26] "PeakSCC3" "PeakSCC" "SCCDO" "MYDO" "MYMEDO"
## [31] "CM1count" "DryDate" "Calv1Date" "Parity" "HerdID"
## [36] "Breed" "Summer" "Path1" "Path2" "Bedtype"
## [41] "East" "Herdsize" "Dryfreestall" "ITS" "Vax"
## [46] "calvedalone" "Aerococ" "Bac" "Coryn" "Entero"
## [51] "LactoALL" "Lactogar" "Lactolac" "Lactosp" "Pseudo"
## [56] "Staaur" "NAS" "Stachr" "Staepi" "Stahaem"
## [61] "Stasci" "Stasim" "Staxyl" "Staspp" "StrepALL"
## [66] "Strepdys" "Strepub" "Strepspp" "Yeast" "Coliform"
## [71] "MiscGP" "MiscGN" "Enrolldate" "IMI" "IMIMAJ"
## [76] "IMIMAJ2" "IMIMIN" "DIMDO" "DIMEN" "RecCM3HT"
## [81] "RecCM14" "RecCM3HTD" "RecCM14D" "En2Calv" "Calvedlate"
## [86] "nocalf" "DPlength" "CM1" "CM1D" "RCens"
## [91] "RCens60" "CMCens" "CM2Cens60" "CM2TAR" "CM2TAR60"
## [96] "CM260" "CMALLRCens" "CMALLTAR" "Rem2Cens" "RemTAR"
## [101] "Rem2TAR" "RemALLRCens" "RemALLTAR" "DOSCM" "SSLO"library(table1)
cow$MYDO <- cow$MYDO * 0.453592
cow$DPlength <- cow$DPlength %>% as.numeric
table1(~ factor(Parity) + MYDO + SCCDO + factor(DOSCM) + PeakSCC + factor(CM1) + factor(Breed), data=cow)| Overall (N=2763) |
|
|---|---|
| factor(Parity) | |
| 1 | 1184 (42.9%) |
| 2 | 881 (31.9%) |
| 3 | 422 (15.3%) |
| 4 | 276 (10.0%) |
| MYDO | |
| Mean (SD) | 29.4 (10.4) |
| Median [Min, Max] | 29.5 [0, 81.2] |
| Missing | 277 (10.0%) |
| SCCDO | |
| Mean (SD) | 4.32 (1.42) |
| Median [Min, Max] | 4.32 [0, 8.99] |
| Missing | 277 (10.0%) |
| factor(DOSCM) | |
| 0 | 1955 (70.8%) |
| 1 | 531 (19.2%) |
| Missing | 277 (10.0%) |
| PeakSCC | |
| Mean (SD) | 5.35 (1.51) |
| Median [Min, Max] | 5.23 [0, 9.21] |
| Missing | 277 (10.0%) |
| factor(CM1) | |
| 0 | 2112 (76.4%) |
| 1 | 438 (15.9%) |
| Missing | 213 (7.7%) |
| factor(Breed) | |
| 0 | 2376 (86.0%) |
| 1 | 387 (14.0%) |
Cows that calved
postcalve$DPlength <- postcalve$DPlength %>% as.numeric
postcalve$En2Calv <- postcalve$En2Calv %>% as.numeric
table1(~ factor(Parity) + MYDO + SCCDO + DOSCM + PeakSCC + factor(CM1) + factor(Breed) + DPlength + En2Calv + factor(CM2) + factor(Rem2), data=postcalve)| Overall (N=2581) |
|
|---|---|
| factor(Parity) | |
| 1 | 1119 (43.4%) |
| 2 | 826 (32.0%) |
| 3 | 390 (15.1%) |
| 4 | 246 (9.5%) |
| MYDO | |
| Mean (SD) | 64.8 (22.9) |
| Median [Min, Max] | 65.0 [0, 179] |
| Missing | 238 (9.2%) |
| SCCDO | |
| Mean (SD) | 4.31 (1.42) |
| Median [Min, Max] | 4.32 [0, 8.99] |
| Missing | 238 (9.2%) |
| DOSCM | |
| 0 | 1850 (71.7%) |
| 1 | 493 (19.1%) |
| Missing | 238 (9.2%) |
| PeakSCC | |
| Mean (SD) | 5.35 (1.51) |
| Median [Min, Max] | 5.23 [0, 9.21] |
| Missing | 238 (9.2%) |
| factor(CM1) | |
| 0 | 2001 (77.5%) |
| 1 | 415 (16.1%) |
| Missing | 165 (6.4%) |
| factor(Breed) | |
| 0 | 2217 (85.9%) |
| 1 | 364 (14.1%) |
| DPlength | |
| Mean (SD) | 57.5 (12.9) |
| Median [Min, Max] | 56.0 [3.00, 129] |
| Missing | 30 (1.2%) |
| En2Calv | |
| Mean (SD) | 71.2 (16.8) |
| Median [Min, Max] | 70.0 [3.00, 129] |
| factor(CM2) | |
| 0 | 2080 (80.6%) |
| 1 | 336 (13.0%) |
| Missing | 165 (6.4%) |
| factor(Rem2) | |
| 0 | 2271 (88.0%) |
| 1 | 310 (12.0%) |
DAG
Below is a DAG that demonstrates my understanding of the causal relationships between observed (and 1 unobserved) variables. I will use this to determine what variables could function as confounders and thus be included as covariates in final models. I will use the same confounder list for all models.
Image.
According to this DAG, I will controll for the following variables in my analyses
Parity [“Parity”]
Herd breed [“Breed”]
Housing –> Calving pen managment [“calvedalone” (0/1)]
Housing –> Dry cow housing [“DryFreestall” (0/1)]
Region of the herd [“East” (0/1)]
Bedding type (Manure solids, New Sand, Other organic, Recycled Sand) [“Bedtype” (0/1/2/3)]
Season of enrollment [“Season”]
I decided to not control for udder health indicators during the lactation prior to enrollment as it is possible that they are on the causal pathway:
Clinical mastitis in previous lactation [“CM1”] Yield at most recent test before dry off [“DOMY”] Somatic cell count at last herd test during previous lactation [“DOSCC” or “PeakSCC”]
Check for correlation between covariates
Pearson correlation matrix among potential predictors
library(psych)
cor <- postcalve %>% select(Parity, East, Breed, Summer, Dryfreestall,calvedalone)
cor$Parity <- cor$Parity %>% as.numeric
corplot <- cor(cor, use = "complete.obs")
corPlot(corplot,numbers=TRUE)
Kendal’s tau
corplot <- cor(cor, use = "complete.obs", method="kendal")
corPlot(corplot,numbers=TRUE)
Spearman
corplot <- cor(cor, use = "complete.obs", method="spearman")
corPlot(corplot,numbers=TRUE)
No variables correlated > 0.7. Will include all covariates in all models.
Outcome 1: Clinical mastitis
120 d Clinical mastitis risk
Kaplan Meier curve & log rank test for clinical mastitis
Kaplan Meier curve
library(ggplot2)
library(survminer)
library(survival)
KM <- survfit(Surv(CM2TAR, CM2) ~ 1, data = postcalve)
summary(KM, times = c(1,30,60,90,120*(1:10)))## Call: survfit(formula = Surv(CM2TAR, CM2) ~ 1, data = postcalve)
##
## 165 observations deleted due to missingness
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 2411 5 0.998 0.000925 0.996 1.000
## 30 2217 93 0.958 0.004110 0.950 0.967
## 60 2063 95 0.917 0.005747 0.905 0.928
## 90 1939 91 0.876 0.006904 0.862 0.889
## 120 1845 52 0.852 0.007462 0.838 0.867ggsurvplot(KM, data = postcalve, title = "", pval = T, conf.int = F,risk.table = T, risk.table.y.text.col = TRUE , surv.plot.height = 5, tables.theme = theme_cleantable(), ggtheme = theme_bw(base_family = "Times"), palette = "Set1")
KMplot <- function(Survfuction){
ggsurvplot(KM, data = postcalve, title = "", pval = T, censor=F, conf.int = F,risk.table = T, surv.plot.height = 5, risk.table.pos="in", tables.theme = theme_cleantable(), ggtheme = theme_bw(base_family = "Times"), xlab="Days in Milk",palette = "Set1",risk.table.fontsize=3)}
KM <- survfit(Surv(CM2TAR, CM2) ~ Staaur, data = postcalve)
Staaur <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ NAS, data = postcalve)
NAS <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ Stachr, data = postcalve)
Stachr <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ Staepi, data = postcalve)
Staepi <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ Stahaem, data = postcalve)
Stahaem <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ Stasci, data = postcalve)
Stasci <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ Stasim, data = postcalve)
Stasim <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ Staxyl, data = postcalve)
Staxyl <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ Staspp, data = postcalve)
Staspp <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ SSLO, data = postcalve)
SSLO <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ Aerococ, data = postcalve)
Aerococ <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ Entero, data = postcalve)
Entero <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ LactoALL, data = postcalve)
LactoALL <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ Lactolac, data = postcalve)
Lactolac <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ Lactogar, data = postcalve)
Lactogar <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ StrepALL, data = postcalve)
StrepALL <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ Strepdys, data = postcalve)
Strepdys <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ Strepub, data = postcalve)
Strepub <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ Strepspp, data = postcalve)
Strepspp <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ Bac, data = postcalve)
Bac <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ Coryn, data = postcalve)
Coryn <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ IMI, data = postcalve)
IMI <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ IMIMAJ, data = postcalve)
IMIMAJ <- KMplot(KM)
KM <- survfit(Surv(CM2TAR, CM2) ~ IMIMIN, data = postcalve)
IMIMIN <- KMplot(KM)
splots <- list(Staaur,NAS,Stachr,Staepi,Stahaem,Stasci,Stasim,Staxyl,Staspp,SSLO,
Aerococ,Entero,LactoALL,Lactogar,Lactolac,
StrepALL,Strepdys,Strepub,Strepspp,
Bac,Coryn,IMI,IMIMAJ,IMIMIN)
arrange_ggsurvplots(splots,ncol=2)
Export KM for IMIMAJ and IMIMIN
library(dplyr)
mm0 <- coxph(Surv(CM2TAR, CM2) ~ Staaur + Stachr + Staepi + Stahaem + Stasci + Stasim + Staxyl + Staspp +
Aerococ + Entero + Lactogar + Lactolac + Strepdys + Strepub + Strepspp +
Bac + Coryn + factor(Summer) + factor(Dryfreestall) + factor(calvedalone) + factor(Bedtype) + strata(Parity) + strata(Breed) + strata(East) + cluster(HerdID), data=postcalve)
mm0 <- coxph(Surv(CM2TAR, CM2) ~ NAS + SSLO + Staaur + Bac + Coryn + factor(Summer) + factor(Dryfreestall) + factor(calvedalone) + factor(Bedtype) + strata(Parity) + strata(Breed) + strata(East) + cluster(HerdID), data=postcalve)
mm0 %>% summary## Call:
## coxph(formula = Surv(CM2TAR, CM2) ~ NAS + SSLO + Staaur + Bac +
## Coryn + factor(Summer) + factor(Dryfreestall) + factor(calvedalone) +
## factor(Bedtype) + strata(Parity) + strata(Breed) + strata(East),
## data = postcalve, cluster = HerdID)
##
## n= 2416, number of events= 336
## (165 observations deleted due to missingness)
##
## coef exp(coef) se(coef) robust se z Pr(>|z|)
## NAS 0.163021 1.177061 0.120292 0.125209 1.302 0.192922
## SSLO 0.393005 1.481426 0.139523 0.139353 2.820 0.004799
## Staaur 1.319522 3.741632 0.313981 0.395239 3.339 0.000842
## Bac -0.318128 0.727510 0.259523 0.257243 -1.237 0.216204
## Coryn 0.241710 1.273425 0.255042 0.285895 0.845 0.397859
## factor(Summer)1 -0.038827 0.961917 0.112020 0.097945 -0.396 0.691796
## factor(Dryfreestall)1 -0.388144 0.678314 0.159974 0.234406 -1.656 0.097750
## factor(calvedalone)1 0.102107 1.107502 0.144818 0.196890 0.519 0.604040
## factor(Bedtype)1 0.011499 1.011566 0.183832 0.252355 0.046 0.963655
## factor(Bedtype)2 -0.092145 0.911973 0.197118 0.258548 -0.356 0.721546
## factor(Bedtype)3 0.001177 1.001178 0.201921 0.259202 0.005 0.996376
##
## NAS
## SSLO **
## Staaur ***
## Bac
## Coryn
## factor(Summer)1
## factor(Dryfreestall)1 .
## factor(calvedalone)1
## factor(Bedtype)1
## factor(Bedtype)2
## factor(Bedtype)3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## NAS 1.1771 0.8496 0.9209 1.504
## SSLO 1.4814 0.6750 1.1274 1.947
## Staaur 3.7416 0.2673 1.7244 8.119
## Bac 0.7275 1.3746 0.4394 1.204
## Coryn 1.2734 0.7853 0.7271 2.230
## factor(Summer)1 0.9619 1.0396 0.7939 1.165
## factor(Dryfreestall)1 0.6783 1.4742 0.4285 1.074
## factor(calvedalone)1 1.1075 0.9029 0.7529 1.629
## factor(Bedtype)1 1.0116 0.9886 0.6169 1.659
## factor(Bedtype)2 0.9120 1.0965 0.5494 1.514
## factor(Bedtype)3 1.0012 0.9988 0.6024 1.664
##
## Concordance= 0.588 (se = 0.024 )
## Likelihood ratio test= 31.1 on 11 df, p=0.001
## Wald test = 37.34 on 11 df, p=1e-04
## Score (logrank) test = 39.76 on 11 df, p=4e-05, Robust = 16.14 p=0.1
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).KM <- survfit(Surv(CM2TAR, CM2) ~ IMIMAJ, data = postcalve)
IMIMAJ <- ggsurvplot(KM, data = postcalve, title = "Major IMI exposure", pval = T, censor=F, conf.int = F,risk.table = T, surv.plot.height = 5, risk.table.pos="in", tables.theme = theme_cleantable(), ggtheme = theme_bw(base_family = "Times"), xlab="Days in Milk",palette = "Set1",risk.table.fontsize=3,legend.title = "",legend.labs = c("Unexposed", "Exposed"),legend = c(0.82, 0.25))
KM <- survfit(Surv(CM2TAR, CM2) ~ IMIMIN, data = postcalve)
IMIMIN <- ggsurvplot(KM, data = postcalve, title = "Minor IMI exposure", pval = T, censor=F, conf.int = F,risk.table = T, surv.plot.height = 5, risk.table.pos="in", tables.theme = theme_cleantable(), ggtheme = theme_bw(base_family = "Times"), xlab="Days in Milk",palette = "Set1",risk.table.fontsize=3,legend.title = "",legend.labs = c("Unexposed", "Exposed"),legend = c(0.82, 0.25))
splots <- list(IMIMAJ,IMIMIN)
#tiff("KM1.tiff", units="in", width=8, height=5, res=300)
arrange_ggsurvplots(splots,ncol=2,nrow=1)
#dev.off()
tiff("KM1.tiff", units="in", width=4, height=9, res=300)
arrange_ggsurvplots(splots,ncol=1,nrow=2)
dev.off()## quartz_off_screen
## 2
Cox proportional hazards regression for clinical mastitis (1-120 DIM)
Model building plan
Model type: Cox proportional hazards analysis with farm included as a cluster variable (robust sandwich standard error estimator) to account for lack of indepedence.
Step 1: Create model with confounders as covariates
Step 2: Investigate if covariates meet proportional hazards assumption for model without IMI as a predictor. Then recheck when fitted with IMI in the model.
Step 3: Report final model
Note that effect measure modification (i.e. interactions) will not be investigated, because of the small number of exposed cows (most are from 30 to 60 cows). Stratification of analyses can lead to very few cases existing within the exposed group (i.e. issues with positivity).
Check proprtional hazards assumption for covariates
mm0 <- coxph(Surv(CM2TAR, CM2) ~ factor(ABX) + factor(Parity) + factor(Breed) +
factor(Dryfreestall) + factor(calvedalone) + factor(Summer) + factor(Bedtype) +
factor(East) + cluster(HerdID), data=postcalve %>% filter(!is.na(CM2)))
SR <- cox.zph(mm0)
SR## chisq df p
## factor(ABX) 1.336 1 0.25
## factor(Parity) 4.034 3 0.26
## factor(Breed) 0.656 1 0.42
## factor(Dryfreestall) 1.445 1 0.23
## factor(calvedalone) 0.523 1 0.47
## factor(Summer) 2.052 1 0.15
## factor(Bedtype) 2.456 3 0.48
## factor(East) 3.535 1 0.06
## GLOBAL 14.454 12 0.27ggcoxzph(SR,var = c("factor(Parity)","factor(Breed)"))
ggcoxzph(SR,var = c("factor(Dryfreestall)","factor(calvedalone)"))
ggcoxzph(SR,var = c("factor(Summer)","factor(East)"))
ggcoxzph(SR,var = c("factor(Bedtype)","factor(ABX)"))
Decision: Will not need to create seperate baseline hazards for any variables.
Create functions to simplify the modelling syntax
library(broom)
coxfunc <- function(expo,exposure,exponame,Other){
OtherCov <- "+ factor(ABX) + factor(Summer) + factor(Dryfreestall) + factor(calvedalone) + factor(Bedtype) + factor(Parity) + factor(Breed) + factor(East) + cluster(HerdID)"
OtherIMI <- c(exposure,Other)
Formula <- formula(paste("Surv(CM2TAR, CM2) ~ ",paste(OtherIMI, collapse=" + "),OtherCov))
mm0 <- coxph(Formula,data=postcalve) %>%
tidy(exp=T) %>% select(term,estimate,conf.low,conf.high) %>% slice(1)
mm0$term <- exponame
mm0$Exposure <- mm0$term
mm0$term <- NULL
a <- table(expo,postcalve$CM2) %>% 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$HR <- 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,HR,estimate,conf.low,conf.high)
}
shoen.resid <- function(expo){
coxph(Surv(CM2TAR, CM2) ~ expo + factor(ABX) + factor(Summer) + factor(Dryfreestall) + factor(calvedalone) +
factor(Bedtype) +
strata(Parity) + strata(Breed) + strata(East) +
cluster(HerdID), data=postcalve) %>% cox.zph
}
shoen.resid.plot <- function(expo){
coxph(Surv(CM2TAR, CM2) ~ expo + factor(ABX) + factor(Summer) + factor(Dryfreestall) + factor(calvedalone) +
factor(Bedtype) +
strata(Parity) + strata(Breed) + strata(East) +
cluster(HerdID), data=postcalve) %>% cox.zph %>% ggcoxzph(var = "expo")
}Exposure = Staphylococcus aureus
shoen.resid(postcalve$Staaur)## chisq df p
## expo 6.7309 1 0.0095
## factor(ABX) 1.5400 1 0.2146
## factor(Summer) 2.1765 1 0.1401
## factor(Dryfreestall) 0.0251 1 0.8741
## factor(calvedalone) 0.2666 1 0.6057
## factor(Bedtype) 3.1970 3 0.3622
## GLOBAL 14.5645 8 0.0682shoen.resid.plot(postcalve$Staaur)
#Staaur <- coxfunc(postcalve$Staaur,"Staphylococcus aureus")
Staaur <- coxfunc(postcalve$Staaur,
exposure="Staaur",
exponame="Staphylococcus aureus",
Other=c("Stachr","strata(Staepi)","Stahaem","Stasci","Stasim","Staxyl","Staspp",
"Aerococ","Entero","LactoALL","StrepALL",
"strata(Bac)","Coryn"))
StaaurDecision: A few outliers. Will keep model.
Exposure = Non aureus Staphylococcus sp. (whole group)
shoen.resid(postcalve$NAS)## chisq df p
## expo 0.00995 1 0.92
## factor(ABX) 1.55529 1 0.21
## factor(Summer) 2.29216 1 0.13
## factor(Dryfreestall) 0.02359 1 0.88
## factor(calvedalone) 0.27866 1 0.60
## factor(Bedtype) 3.19520 3 0.36
## GLOBAL 7.31138 8 0.50shoen.resid.plot(postcalve$NAS)
#NAS <- coxfunc(postcalve$NAS,"Non aureus Staphylococcus sp.")
NAS<- coxfunc(postcalve$NAS,
exposure="NAS",
exponame="Non aureus Staphylococcus sp.",
Other=c("strata(Staaur)","Aerococ","Entero","LactoALL","StrepALL",
"strata(Bac)","Coryn"))
NAS
Exposure = Staphylococcus chromogenes
shoen.resid(postcalve$Stachr)## chisq df p
## expo 2.899 1 0.089
## factor(ABX) 1.607 1 0.205
## factor(Summer) 2.322 1 0.128
## factor(Dryfreestall) 0.030 1 0.863
## factor(calvedalone) 0.271 1 0.602
## factor(Bedtype) 3.225 3 0.358
## GLOBAL 10.484 8 0.233shoen.resid.plot(postcalve$Stachr)
#Stachr <- coxfunc(postcalve$Stachr,"Staphylococcus chromogenes")
Stachr <- coxfunc(postcalve$Stachr,
exposure="Stachr",
exponame="Staphylococcus chromogenes",
Other=c("strata(Staaur)","strata(Staepi)","Stahaem","Stasci","Stasim","Staxyl","Staspp",
"Aerococ","Entero","LactoALL","StrepALL",
"strata(Bac)","Coryn"))
Stachr
Exposure = Staphylococcus epidermidis
shoen.resid(postcalve$Staepi)## chisq df p
## expo 1.9327 1 0.16
## factor(ABX) 1.5795 1 0.21
## factor(Summer) 2.3532 1 0.13
## factor(Dryfreestall) 0.0239 1 0.88
## factor(calvedalone) 0.2639 1 0.61
## factor(Bedtype) 3.2019 3 0.36
## GLOBAL 9.8692 8 0.27shoen.resid.plot(postcalve$Staepi)
#Staepi <- coxfunc(postcalve$Staepi,"Staphylococcus epidermidis")
Staepi <- coxfunc(postcalve$Staepi,
exposure="Staepi",
exponame="Staphylococcus epidermidis",
Other=c("strata(Staaur)","Stachr","Stahaem","Stasci","Stasim","Staxyl","Staspp",
"Aerococ","Entero","LactoALL","StrepALL",
"strata(Bac)","Coryn"))
Staepi
Exposure = Staphylococcus haemolyticus
shoen.resid(postcalve$Stahaem)## chisq df p
## expo 0.7186 1 0.40
## factor(ABX) 1.5498 1 0.21
## factor(Summer) 2.3018 1 0.13
## factor(Dryfreestall) 0.0337 1 0.85
## factor(calvedalone) 0.2662 1 0.61
## factor(Bedtype) 3.1934 3 0.36
## GLOBAL 8.4405 8 0.39shoen.resid.plot(postcalve$Stahaem)
#Stahaem <- coxfunc(postcalve$Stahaem,"Staphylococcus haemolyticus")
Stahaem <- coxfunc(postcalve$Stahaem,
exposure="Stahaem",
exponame="Staphylococcus haemolyticus",
Other=c("strata(Staaur)","Stachr","strata(Staepi)","Stasci","Stasim","Staxyl","Staspp",
"Aerococ","Entero","LactoALL","StrepALL",
"strata(Bac)","Coryn"))
Stahaem
Exposure = Staphylococcus sciuri
shoen.resid(postcalve$Stasci)## chisq df p
## expo 1.3312 1 0.25
## factor(ABX) 1.5786 1 0.21
## factor(Summer) 2.2951 1 0.13
## factor(Dryfreestall) 0.0268 1 0.87
## factor(calvedalone) 0.2726 1 0.60
## factor(Bedtype) 3.1953 3 0.36
## GLOBAL 8.3492 8 0.40shoen.resid.plot(postcalve$Stasci)
#Stasci <- coxfunc(postcalve$Stasci,"Staphylococcus sciuri")
Stasci <- coxfunc(postcalve$Stasci,
exposure="Stasci",
exponame="Staphylococcus sciuri",
Other=c("strata(Staaur)","Stachr","strata(Staepi)","Stahaem","Stasci","Stasim","Staxyl","Staspp",
"Aerococ","Entero","LactoALL","StrepALL",
"strata(Bac)","Coryn"))
Stasci
Exposure = Staphylococcus simulans
shoen.resid(postcalve$Stasim)## chisq df p
## expo 0.9104 1 0.34
## factor(ABX) 1.5749 1 0.21
## factor(Summer) 2.3184 1 0.13
## factor(Dryfreestall) 0.0281 1 0.87
## factor(calvedalone) 0.2678 1 0.60
## factor(Bedtype) 3.1657 3 0.37
## GLOBAL 8.2807 8 0.41shoen.resid.plot(postcalve$Stasim)
#Stasim <- coxfunc(postcalve$Stasim,"Staphylococcus simulans")
Stasim <- coxfunc(postcalve$Stasim,
exposure="Stasim",
exponame="Staphylococcus simulans",
Other=c("strata(Staaur)","Stachr","strata(Staepi)","Stahaem","Stasci","Stasim","Staxyl","Staspp",
"Aerococ","Entero","LactoALL","StrepALL",
"strata(Bac)","Coryn"))
Stasim
Exposure = Staphylococcus xylosus
shoen.resid(postcalve$Staxyl)## chisq df p
## expo 0.2499 1 0.62
## factor(ABX) 1.6056 1 0.21
## factor(Summer) 2.3830 1 0.12
## factor(Dryfreestall) 0.0173 1 0.90
## factor(calvedalone) 0.2627 1 0.61
## factor(Bedtype) 3.2699 3 0.35
## GLOBAL 8.3421 8 0.40shoen.resid.plot(postcalve$Staxyl)
#Staxyl <- coxfunc(postcalve$Staxyl,"Staphylococcus xylosus")
Staxyl <- coxfunc(postcalve$Staxyl,
exposure="Staxyl",
exponame="Staphylococcus xylosus",
Other=c("strata(Staaur)","Stachr","strata(Staepi)","Stahaem","Stasci","Stasim","Staxyl","Staspp",
"Aerococ","Entero","LactoALL","StrepALL",
"strata(Bac)","Coryn"))
Staxyl
Exposure = Staphylococcus species (i.e. species unable to be determined using MALDI-TOF)
shoen.resid(postcalve$Staspp)## chisq df p
## expo 0.8147 1 0.37
## factor(ABX) 1.5897 1 0.21
## factor(Summer) 2.2957 1 0.13
## factor(Dryfreestall) 0.0259 1 0.87
## factor(calvedalone) 0.2963 1 0.59
## factor(Bedtype) 3.2346 3 0.36
## GLOBAL 8.3831 8 0.40shoen.resid.plot(postcalve$Staspp)
#Staspp <- coxfunc(postcalve$Staspp,"Staphylococcus species")
Staspp <- coxfunc(postcalve$Staspp,
exposure="Staspp",
exponame="Staphylococcus sp.",
Other=c("strata(Staaur)","Stachr","strata(Staepi)","Stahaem","Stasci","Stasim","Staxyl","Staspp",
"Aerococ","Entero","LactoALL","StrepALL",
"strata(Bac)","Coryn"))
Staspp
### Exposure = SSLO
shoen.resid(postcalve$SSLO)## chisq df p
## expo 0.8438 1 0.36
## factor(ABX) 1.5996 1 0.21
## factor(Summer) 2.3006 1 0.13
## factor(Dryfreestall) 0.0393 1 0.84
## factor(calvedalone) 0.2813 1 0.60
## factor(Bedtype) 3.2948 3 0.35
## GLOBAL 7.9831 8 0.44shoen.resid.plot(postcalve$SSLO)
#SSLO <- coxfunc(postcalve$SSLO,"SSLO")
SSLO <- coxfunc(postcalve$SSLO,
exposure="SSLO",
exponame="SSLO",
Other=c("strata(Staaur)","Stachr","strata(Staepi)","Stahaem","Stasci","Stasim","Staxyl","Staspp",
"strata(Bac)","Coryn"))
SSLO
Exposure = Aerococcus species
shoen.resid(postcalve$Aerococ)## chisq df p
## expo 1.3420 1 0.25
## factor(ABX) 1.5975 1 0.21
## factor(Summer) 2.3255 1 0.13
## factor(Dryfreestall) 0.0293 1 0.86
## factor(calvedalone) 0.2740 1 0.60
## factor(Bedtype) 3.2147 3 0.36
## GLOBAL 8.5464 8 0.38shoen.resid.plot(postcalve$Aerococ)
#Aerococ <- coxfunc(postcalve$Aerococ,"Aerococcus sp.")
Aerococ <- coxfunc(postcalve$Aerococ,
exposure="Aerococ",
exponame="Aerococcus spp.",
Other=c("strata(Staaur)","Stachr","strata(Staepi)","Stahaem","Stasci","Stasim","Staxyl","Staspp",
"Aerococ","Entero","LactoALL","StrepALL",
"strata(Bac)","Coryn"))
Aerococ
Exposure = Enterococcus species
shoen.resid(postcalve$Entero)## chisq df p
## expo 0.0512 1 0.82
## factor(ABX) 1.5740 1 0.21
## factor(Summer) 2.3049 1 0.13
## factor(Dryfreestall) 0.0301 1 0.86
## factor(calvedalone) 0.2766 1 0.60
## factor(Bedtype) 3.1448 3 0.37
## GLOBAL 7.3866 8 0.50shoen.resid.plot(postcalve$Entero)
#Entero <- coxfunc(postcalve$Entero,"Enterococcus sp.")
Entero <- coxfunc(postcalve$Entero,
exposure="Entero",
exponame="Enterococcus spp.",
Other=c("strata(Staaur)","Stachr","strata(Staepi)","Stahaem","Stasci","Stasim","Staxyl","Staspp",
"Aerococ","Entero","LactoALL","StrepALL",
"strata(Bac)","Coryn"))
Entero
Exposure = Lactococcus (all species combined)
shoen.resid(postcalve$LactoALL)## chisq df p
## expo 0.1394 1 0.71
## factor(ABX) 1.5834 1 0.21
## factor(Summer) 2.3075 1 0.13
## factor(Dryfreestall) 0.0293 1 0.86
## factor(calvedalone) 0.2695 1 0.60
## factor(Bedtype) 3.2285 3 0.36
## GLOBAL 7.3857 8 0.50shoen.resid.plot(postcalve$LactoALL)
#LactoALL <- coxfunc(postcalve$LactoALL,"Lactococcus spp.")
LactoALL <- coxfunc(postcalve$LactoALL,
exposure="LactoALL",
exponame="Lactococcus spp.",
Other=c("strata(Staaur)","Stachr","strata(Staepi)","Stahaem","Stasci","Stasim","Staxyl","Staspp",
"Aerococ","Entero","LactoALL","StrepALL",
"strata(Bac)","Coryn"))
LactoALL
Exposure = Lactococcus garviae
shoen.resid(postcalve$Lactogar)## chisq df p
## expo 0.0477 1 0.83
## factor(ABX) 1.5837 1 0.21
## factor(Summer) 2.3061 1 0.13
## factor(Dryfreestall) 0.0293 1 0.86
## factor(calvedalone) 0.2687 1 0.60
## factor(Bedtype) 3.1986 3 0.36
## GLOBAL 7.3403 8 0.50shoen.resid.plot(postcalve$Lactogar)
#Lactogar <- coxfunc(postcalve$Lactogar,"Lactococcus garviae")
Lactogar <- coxfunc(postcalve$Lactogar,
exposure="Lactogar",
exponame="Lactococcus garviae",
Other=c("strata(Staaur)","Stachr","strata(Staepi)","Stahaem","Stasci","Stasim","Staxyl","Staspp",
"Aerococ","Entero","Lactolac","Lactosp","StrepALL",
"strata(Bac)","Coryn"))
Lactogar
Exposure = Lactococcus lactis
shoen.resid(postcalve$Lactolac)## chisq df p
## expo 0.1023 1 0.75
## factor(ABX) 1.5822 1 0.21
## factor(Summer) 2.3037 1 0.13
## factor(Dryfreestall) 0.0275 1 0.87
## factor(calvedalone) 0.2666 1 0.61
## factor(Bedtype) 3.2202 3 0.36
## GLOBAL 7.3760 8 0.50shoen.resid.plot(postcalve$Lactolac)
#Lactolac <- coxfunc(postcalve$Lactolac,"Lactococcus lactis")
Lactolac <- coxfunc(postcalve$Lactolac,
exposure="Lactolac",
exponame="Lactococcus lactis",
Other=c("strata(Staaur)","Stachr","strata(Staepi)","Stahaem","Stasci","Stasim","Staxyl","Staspp",
"Aerococ","Entero","Lactogar","StrepALL",
"strata(Bac)","Coryn"))
Lactolac
Exposure = Streptococcus
shoen.resid(postcalve$StrepALL)## chisq df p
## expo 0.4777 1 0.49
## factor(ABX) 1.5684 1 0.21
## factor(Summer) 2.2667 1 0.13
## factor(Dryfreestall) 0.0348 1 0.85
## factor(calvedalone) 0.2900 1 0.59
## factor(Bedtype) 3.2385 3 0.36
## GLOBAL 7.9249 8 0.44shoen.resid.plot(postcalve$StrepALL)
#StrepALL <- coxfunc(postcalve$StrepALL,"Streptococcus spp.")
StrepALL <- coxfunc(postcalve$StrepALL,
exposure="StrepALL",
exponame="Streptococcus spp.",
Other=c("strata(Staaur)","Stachr","strata(Staepi)","Stahaem","Stasci","Stasim","Staxyl","Staspp",
"Aerococ","Entero","LactoALL","StrepALL",
"strata(Bac)","Coryn"))
StrepALL
Exposure = Streptococcus dysgalactiae
shoen.resid(postcalve$Strepdys)## chisq df p
## expo 0.0545 1 0.82
## factor(ABX) 1.5763 1 0.21
## factor(Summer) 2.3214 1 0.13
## factor(Dryfreestall) 0.0285 1 0.87
## factor(calvedalone) 0.2803 1 0.60
## factor(Bedtype) 3.1723 3 0.37
## GLOBAL 7.4253 8 0.49shoen.resid.plot(postcalve$Strepdys)
#Strepdys <- coxfunc(postcalve$Strepdys,"Streptococcus dysgalactiae")
Strepdys <- coxfunc(postcalve$Strepdys,
exposure="Strepdys",
exponame="Streptococcus dysgalactiae",
Other=c("strata(Staaur)","Stachr","strata(Staepi)","Stahaem","Stasci","Stasim","Staxyl","Staspp",
"Aerococ","Entero","LactoALL","Strepub","Strepspp",
"strata(Bac)","Coryn"))
Strepdys
Exposure = Streptococcus uberis
shoen.resid(postcalve$Strepub)## chisq df p
## expo 0.0162 1 0.90
## factor(ABX) 1.5741 1 0.21
## factor(Summer) 2.2797 1 0.13
## factor(Dryfreestall) 0.0290 1 0.86
## factor(calvedalone) 0.2779 1 0.60
## factor(Bedtype) 3.1950 3 0.36
## GLOBAL 7.3346 8 0.50shoen.resid.plot(postcalve$Strepub)
#Strepub <- coxfunc(postcalve$Strepub,"Streptococcus uberis")
Strepub <- coxfunc(postcalve$Strepub,
exposure="Strepub",
exponame="Streptococcus uberis",
Other=c("strata(Staaur)","Stachr","strata(Staepi)","Stahaem","Stasci","Stasim","Staxyl","Staspp",
"Aerococ","Entero","LactoALL","Strepdys","Strepspp",
"strata(Bac)","Coryn"))
Strepub
Exposure = Strep species (unable to be determined using MALDI-TOF)
shoen.resid(postcalve$Strepspp)## chisq df p
## expo 0.7026 1 0.40
## factor(ABX) 1.5772 1 0.21
## factor(Summer) 2.2848 1 0.13
## factor(Dryfreestall) 0.0332 1 0.86
## factor(calvedalone) 0.2805 1 0.60
## factor(Bedtype) 3.2720 3 0.35
## GLOBAL 8.2264 8 0.41shoen.resid.plot(postcalve$Strepspp)
#Strepspp <- coxfunc(postcalve$Strepspp,"Streptococcus sp.")
Strepspp <- coxfunc(postcalve$Strepspp,
exposure="Strepspp",
exponame="Streptococcus sp.",
Other=c("strata(Staaur)","Stachr","strata(Staepi)","Stahaem","Stasci","Stasim","Staxyl","Staspp",
"Aerococ","Entero","LactoALL","Strepub","Strepdys",
"strata(Bac)","Coryn"))
Strepspp
Exposure = Bacillus
shoen.resid(postcalve$Bac)## chisq df p
## expo 2.5176 1 0.11
## factor(ABX) 1.5645 1 0.21
## factor(Summer) 2.3066 1 0.13
## factor(Dryfreestall) 0.0262 1 0.87
## factor(calvedalone) 0.2692 1 0.60
## factor(Bedtype) 3.2257 3 0.36
## GLOBAL 10.2679 8 0.25shoen.resid.plot(postcalve$Bac)
#Bac <- coxfunc(postcalve$Bac,"Bacillus sp.")
Bac <- coxfunc(postcalve$Bac,
exposure="Bac",
exponame="Bacillus spp.",
Other=c("strata(Staaur)","Stachr","strata(Staepi)","Stahaem","Stasci","Stasim","Staxyl","Staspp",
"Aerococ","Entero","LactoALL","Strepub","Strepdys",
"Coryn"))
Bac
Exposure = Corynebacterium
shoen.resid(postcalve$Coryn)## chisq df p
## expo 0.2440 1 0.62
## factor(ABX) 1.5806 1 0.21
## factor(Summer) 2.3000 1 0.13
## factor(Dryfreestall) 0.0271 1 0.87
## factor(calvedalone) 0.2787 1 0.60
## factor(Bedtype) 3.2094 3 0.36
## GLOBAL 7.8865 8 0.44shoen.resid.plot(postcalve$Coryn)
#Coryn <- coxfunc(postcalve$Coryn,"Corynebacterium sp.")
Coryn <- coxfunc(postcalve$Coryn,
exposure="Coryn",
exponame="Corynebacterium spp.",
Other=c("strata(Staaur)","Stachr","strata(Staepi)","Stahaem","Stasci","Stasim","Staxyl","Staspp",
"Aerococ","Entero","LactoALL","Strepub","Strepdys",
"strata(Bac)","Coryn"))
CorynExposure = All IMI
shoen.resid(postcalve$IMI)## chisq df p
## expo 1.5278 1 0.22
## factor(ABX) 1.5718 1 0.21
## factor(Summer) 2.2840 1 0.13
## factor(Dryfreestall) 0.0267 1 0.87
## factor(calvedalone) 0.2825 1 0.60
## factor(Bedtype) 3.2287 3 0.36
## GLOBAL 8.5862 8 0.38shoen.resid.plot(postcalve$IMI)
#IMIALL <- coxfunc(postcalve$IMI,"All pathogens")
IMIALL <- coxfunc(postcalve$IMI,
exposure="IMI",
exponame="All pathogens",
Other=c())
IMIALL
Exposure = Major pathogens
shoen.resid(postcalve$IMIMAJ)## chisq df p
## expo 2.3733 1 0.12
## factor(ABX) 1.5952 1 0.21
## factor(Summer) 2.2834 1 0.13
## factor(Dryfreestall) 0.0378 1 0.85
## factor(calvedalone) 0.2814 1 0.60
## factor(Bedtype) 3.3303 3 0.34
## GLOBAL 9.2661 8 0.32shoen.resid.plot(postcalve$IMIMAJ)
#IMIMAJ <- coxfunc(postcalve$IMIMAJ,"Major pathogens")
IMIMAJ <- coxfunc(postcalve$IMIMAJ,
exposure="IMIMAJ",
exponame="Major pathogens",
Other=c("IMIMIN"))
IMIMAJ
Exposure = Minor pathogens (NAS, Bacillus and Corynebacterium)
shoen.resid(postcalve$IMIMIN)## chisq df p
## expo 1.2962 1 0.25
## factor(ABX) 1.5680 1 0.21
## factor(Summer) 2.2961 1 0.13
## factor(Dryfreestall) 0.0238 1 0.88
## factor(calvedalone) 0.2807 1 0.60
## factor(Bedtype) 3.1758 3 0.37
## GLOBAL 8.3898 8 0.40shoen.resid.plot(postcalve$IMIMIN)
#IMIMIN <- coxfunc(postcalve$IMIMIN,"Minor pathogens")
IMIMIN <- coxfunc(postcalve$IMIMIN,
exposure="IMIMIN",
exponame="Minor pathogens",
Other=c("IMIMAJ"))
IMIMINOutcome 1: Final models
Blank <- Staaur
Blank$Exposure <- NA
Blank$RiskExp0 <- NA
Blank$RiskExp1 <- NA
Blank$HR <- 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
## 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("Hazard ratio","\n",FP$HR))
library(forestplot)
tiff("FP1.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= TRUE,
mean=c(NA,NA,FP$estimate),
lower=c(NA,NA,FP$conf.low), upper=c(NA,NA,FP$conf.high),
title="",
xlab="<---Negative association-- --Positive association--->",
#hrzl_lines=list("3" = gpar(lwd=1, col="#000000"),
# "7" = gpar(lwd=110, 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
## 2Outcome 2: Culling risk 0-120 DIM
Kaplan Meier curve & log rank test for culling risk
Kaplan Meier curve
library(ggplot2)
library(survminer)
library(survival)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ 1, data = postcalve)
summary(KM, times = c(1,30,60,90,120*(1:10)))## Call: survfit(formula = Surv(Rem2TAR, Rem2) ~ 1, data = postcalve)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 2578 4 0.998 0.000774 0.997 1.000
## 30 2446 137 0.945 0.004473 0.937 0.954
## 60 2375 68 0.919 0.005370 0.909 0.930
## 90 2324 50 0.900 0.005914 0.888 0.911
## 120 2271 51 0.880 0.006399 0.867 0.893ggsurvplot(KM, data = postcalve, title = "", pval = T, conf.int = F,risk.table = T, risk.table.y.text.col = TRUE , surv.plot.height = 5, tables.theme = theme_cleantable(), ggtheme = theme_bw(base_family = "Times"), palette = "Set1")
KMplot <- function(Survfuction){
ggsurvplot(KM, data = postcalve, title = "", pval = T, censor=F, conf.int = F,risk.table = T, surv.plot.height = 5, risk.table.pos="in", tables.theme = theme_cleantable(), ggtheme = theme_bw(base_family = "Times"), palette = "Set1",risk.table.fontsize=3)}
KM <- survfit(Surv(Rem2TAR, Rem2) ~ Staaur, data = postcalve)
Staaur <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ NAS, data = postcalve)
NAS <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ Stachr, data = postcalve)
Stachr <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ Staepi, data = postcalve)
Staepi <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ Stahaem, data = postcalve)
Stahaem <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ Stasci, data = postcalve)
Stasci <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ Stasim, data = postcalve)
Stasim <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ Staxyl, data = postcalve)
Staxyl <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ Staspp, data = postcalve)
Staspp <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ SSLO, data = postcalve)
SSLO <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ Aerococ, data = postcalve)
Aerococ <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ Entero, data = postcalve)
Entero <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ LactoALL, data = postcalve)
LactoALL <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ Lactolac, data = postcalve)
Lactolac <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ Lactogar, data = postcalve)
Lactogar <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ StrepALL, data = postcalve)
StrepALL <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ Strepdys, data = postcalve)
Strepdys <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ Strepub, data = postcalve)
Strepub <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ Strepspp, data = postcalve)
Strepspp <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ Bac, data = postcalve)
Bac <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ Coryn, data = postcalve)
Coryn <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ IMI, data = postcalve)
IMI <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ IMIMAJ, data = postcalve)
IMIMAJ <- KMplot(KM)
KM <- survfit(Surv(Rem2TAR, Rem2) ~ IMIMIN, data = postcalve)
IMIMIN <- KMplot(KM)
splots <- list(Staaur,NAS,Stachr,Staepi,Stahaem,Stasci,Stasim,Staxyl,Staspp,
Aerococ,Entero,LactoALL,Lactogar,Lactolac,
StrepALL,Strepdys,Strepub,Strepspp,
Bac,Coryn,IMI,IMIMAJ,IMIMIN)
arrange_ggsurvplots(splots,ncol=2)
Export KM for IMIMAJ and IMIMIN
KM <- survfit(Surv(Rem2TAR, Rem2) ~ IMIMAJ, data = postcalve)
IMIMAJ <- ggsurvplot(KM, data = postcalve, title = "Major IMI exposure", pval = T, censor=F, conf.int = F,risk.table = T, surv.plot.height = 5, risk.table.pos="in", tables.theme = theme_cleantable(), ggtheme = theme_bw(base_family = "Times"), xlab="Days in Milk",palette = "Set1",risk.table.fontsize=3,legend.title = "",legend.labs = c("Unexposed", "Exposed"),legend = c(0.82, 0.25))
KM <- survfit(Surv(Rem2TAR, Rem2) ~ IMIMIN, data = postcalve)
IMIMIN <- ggsurvplot(KM, data = postcalve, title = "Minor IMI exposure", pval = T, censor=F, conf.int = F,risk.table = T, surv.plot.height = 5, risk.table.pos="in", tables.theme = theme_cleantable(), ggtheme = theme_bw(base_family = "Times"), xlab="Days in Milk",palette = "Set1",risk.table.fontsize=3,legend.title = "",legend.labs = c("Unexposed", "Exposed"),legend = c(0.82, 0.25))
splots <- list(IMIMAJ,IMIMIN)
#tiff("KM1.tiff", units="in", width=8, height=5, res=300)
arrange_ggsurvplots(splots,ncol=2,nrow=1)
#dev.off()
tiff("KM2.tiff", units="in", width=4, height=9, res=300)
arrange_ggsurvplots(splots,ncol=1,nrow=2)
dev.off()## quartz_off_screen
## 2
Cox proportional hazards regression for culling or death (1-120 DIM)
Model building plan
Model type: Cox proportional hazards analysis with farm included as a cluster variable (robust sandwich standard error estimator) to account for lack of indepedence.
Step 1: Create model with confounders as covariates
Step 2: Investigate if covariates meet proportional hazards assumption for model without IMI as a predictor. Then recheck when fitted with IMI in the model.
Step 3: Report final model
Note that effect measure modification (i.e. interactions) will not be investigated, because of the small number of exposed cows (most are from 30 to 60 cows). Stratification of analyses can lead to very few cases existing within the exposed group (i.e. issues with positivity).
Check proprtional hazards assumption for covariates
mm0 <- coxph(Surv(Rem2TAR, Rem2) ~ factor(ABX) + factor(Parity) + factor(Breed) + factor(Bedtype) +
factor(Dryfreestall) + factor(calvedalone) + factor(Summer) +
factor(East) + cluster(HerdID), data=postcalve %>% filter(!is.na(Rem2)))
SR <- cox.zph(mm0)
SR## chisq df p
## factor(ABX) 2.05797 1 0.1514
## factor(Parity) 1.48998 3 0.6846
## factor(Breed) 0.00545 1 0.9411
## factor(Bedtype) 3.84025 3 0.2792
## factor(Dryfreestall) 0.43082 1 0.5116
## factor(calvedalone) 2.65207 1 0.1034
## factor(Summer) 8.37275 1 0.0038
## factor(East) 0.32917 1 0.5662
## GLOBAL 18.85054 12 0.0922ggcoxzph(SR,var = c("factor(Parity)","factor(Breed)"))
ggcoxzph(SR,var = c("factor(Dryfreestall)","factor(calvedalone)"))
ggcoxzph(SR,var = c("factor(Summer)","factor(East)"))
ggcoxzph(SR,var = c("factor(Bedtype)","factor(ABX)"))
Decision: Will create separate baselines for summer variable for all Rem120 cox models.
Create functions to simplify the modelling syntax
library(broom)
coxfunc <- function(expo,exposure,exponame,Other){
OtherCov <- "+ factor(ABX) + strata(Summer) + factor(Dryfreestall) + factor(calvedalone) + factor(Bedtype) + factor(Parity) + factor(Breed) + factor(East) + cluster(HerdID)"
OtherIMI <- c(exposure,Other)
Formula <- formula(paste("Surv(Rem2TAR, Rem2) ~ ",paste(OtherIMI, collapse=" + "),OtherCov))
mm0 <- coxph(Formula,data=postcalve) %>%
tidy(exp=T) %>% select(term,estimate,conf.low,conf.high) %>% slice(1)
mm0$term <- exponame
mm0$Exposure <- mm0$term
mm0$term <- NULL
a <- table(expo,postcalve$Rem2) %>% 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$HR <- 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,HR,estimate,conf.low,conf.high)
}
shoen.resid <- function(expo){
coxph(Surv(Rem2TAR, Rem2) ~ expo +
strata(ABX) + factor(Dryfreestall) + factor(East) + factor(Parity) + factor(Breed) + factor(Bedtype) +
strata(calvedalone) + strata(Summer) +
cluster(HerdID), data=postcalve) %>% cox.zph
}
shoen.resid.plot <- function(expo){
coxph(Surv(Rem2TAR, Rem2) ~ expo +
strata(ABX) + factor(Dryfreestall) + factor(East) + factor(Parity) + factor(Breed) + factor(Bedtype) +
strata(calvedalone) + strata(Summer) +
cluster(HerdID), data=postcalve) %>% cox.zph %>% ggcoxzph(var = "expo")
}Exposure = Staphylococcus aureus
shoen.resid(postcalve$Staaur)## chisq df p
## expo 3.7322 1 0.053
## factor(Dryfreestall) 0.1419 1 0.706
## factor(East) 0.0348 1 0.852
## factor(Parity) 1.5403 3 0.673
## factor(Breed) 0.0127 1 0.910
## factor(Bedtype) 3.7092 3 0.295
## GLOBAL 9.8421 10 0.454shoen.resid.plot(postcalve$Staaur)
#Staaur <- coxfunc(postcalve$Staaur,"Staphylococcus aureus")
Staaur <- coxfunc(postcalve$Staaur,
exposure="Staaur",
exponame="Staphylococcus aureus",
Other=c("Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","strata(Staspp)",
"Aerococ","strata(Entero)","LactoALL","StrepALL",
"Bac","strata(Coryn)"))
StaaurDecision: A few outliers. Will keep model.
Exposure = Non aureus Staphylococcus sp. (whole group)
shoen.resid(postcalve$NAS)## chisq df p
## expo 0.89185 1 0.34
## factor(Dryfreestall) 0.14454 1 0.70
## factor(East) 0.03326 1 0.86
## factor(Parity) 1.49917 3 0.68
## factor(Breed) 0.00975 1 0.92
## factor(Bedtype) 3.73018 3 0.29
## GLOBAL 6.76555 10 0.75shoen.resid.plot(postcalve$NAS)
#NAS <- coxfunc(postcalve$NAS,"Non aureus Staphylococcus sp.")
NAS<- coxfunc(postcalve$NAS,
exposure="NAS",
exponame="Non aureus Staphylococcus sp.",
Other=c("strata(Staaur)","factor(Staaur)","Aerococ","strata(Entero)","LactoALL","StrepALL",
"Bac","strata(Coryn)"))
NAS
Exposure = Staphylococcus chromogenes
shoen.resid(postcalve$Stachr)## chisq df p
## expo 0.9935 1 0.32
## factor(Dryfreestall) 0.1421 1 0.71
## factor(East) 0.0346 1 0.85
## factor(Parity) 1.5447 3 0.67
## factor(Breed) 0.0128 1 0.91
## factor(Bedtype) 3.7071 3 0.29
## GLOBAL 7.1078 10 0.72shoen.resid.plot(postcalve$Stachr)
#Stachr <- coxfunc(postcalve$Stachr,"Staphylococcus chromogenes")
Stachr <- coxfunc(postcalve$Stachr,
exposure="Stachr",
exponame="Staphylococcus chromogenes",
Other=c("strata(Staaur)","Staepi","Stahaem","Stasci","Stasim","Staxyl","strata(Staspp)",
"Aerococ","strata(Entero)","LactoALL","StrepALL",
"Bac","strata(Coryn)"))
Stachr
Exposure = Staphylococcus epidermidis
shoen.resid(postcalve$Staepi)## chisq df p
## expo 0.1446 1 0.70
## factor(Dryfreestall) 0.1448 1 0.70
## factor(East) 0.0321 1 0.86
## factor(Parity) 1.5323 3 0.67
## factor(Breed) 0.0118 1 0.91
## factor(Bedtype) 3.7247 3 0.29
## GLOBAL 6.1472 10 0.80shoen.resid.plot(postcalve$Staepi)
#Staepi <- coxfunc(postcalve$Staepi,"Staphylococcus epidermidis")
Staepi <- coxfunc(postcalve$Staepi,
exposure="Staepi",
exponame="Staphylococcus epidermidis",
Other=c("strata(Staaur)","Stachr","Stahaem","Stasci","Stasim","Staxyl","strata(Staspp)",
"Aerococ","strata(Entero)","LactoALL","StrepALL",
"Bac","strata(Coryn)"))
Staepi
Exposure = Staphylococcus haemolyticus
shoen.resid(postcalve$Stahaem)## chisq df p
## expo 0.8521 1 0.36
## factor(Dryfreestall) 0.1420 1 0.71
## factor(East) 0.0343 1 0.85
## factor(Parity) 1.5409 3 0.67
## factor(Breed) 0.0127 1 0.91
## factor(Bedtype) 3.7095 3 0.29
## GLOBAL 6.6988 10 0.75shoen.resid.plot(postcalve$Stahaem)
#Stahaem <- coxfunc(postcalve$Stahaem,"Staphylococcus haemolyticus")
Stahaem <- coxfunc(postcalve$Stahaem,
exposure="Stahaem",
exponame="Staphylococcus haemolyticus",
Other=c("strata(Staaur)","factor(Staaur)","Stachr","Staepi","Stasci","Stasim","Staxyl","strata(Staspp)",
"Aerococ","strata(Entero)","LactoALL","StrepALL",
"Bac","strata(Coryn)"))
Stahaem
Exposure = Staphylococcus sciuri
shoen.resid(postcalve$Stasci)## chisq df p
## expo 2.5031 1 0.11
## factor(Dryfreestall) 0.1405 1 0.71
## factor(East) 0.0359 1 0.85
## factor(Parity) 1.5132 3 0.68
## factor(Breed) 0.0122 1 0.91
## factor(Bedtype) 3.6826 3 0.30
## GLOBAL 8.5661 10 0.57shoen.resid.plot(postcalve$Stasci)
#Stasci <- coxfunc(postcalve$Stasci,"Staphylococcus sciuri")
Stasci <- coxfunc(postcalve$Stasci,
exposure="Stasci",
exponame="Staphylococcus sciuri",
Other=c("strata(Staaur)","strata(Staaur)","Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","strata(Staspp)",
"Aerococ","strata(Entero)","LactoALL","StrepALL",
"Bac","strata(Coryn)"))
Stasci
Exposure = Staphylococcus simulans
shoen.resid(postcalve$Stasim)## chisq df p
## expo 0.6549 1 0.42
## factor(Dryfreestall) 0.1429 1 0.71
## factor(East) 0.0342 1 0.85
## factor(Parity) 1.5351 3 0.67
## factor(Breed) 0.0129 1 0.91
## factor(Bedtype) 3.7047 3 0.30
## GLOBAL 6.3496 10 0.79shoen.resid.plot(postcalve$Stasim)
#Stasim <- coxfunc(postcalve$Stasim,"Staphylococcus simulans")
Stasim <- coxfunc(postcalve$Stasim,
exposure="Stasim",
exponame="Staphylococcus simulans",
Other=c("strata(Staaur)","Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","strata(Staspp)",
"Aerococ","strata(Entero)","LactoALL","StrepALL",
"Bac","strata(Coryn)"))
Stasim
Exposure = Staphylococcus xylosus
shoen.resid(postcalve$Staxyl)## chisq df p
## expo 0.1646 1 0.68
## factor(Dryfreestall) 0.1423 1 0.71
## factor(East) 0.0344 1 0.85
## factor(Parity) 1.5444 3 0.67
## factor(Breed) 0.0125 1 0.91
## factor(Bedtype) 3.7135 3 0.29
## GLOBAL 6.1294 10 0.80shoen.resid.plot(postcalve$Staxyl)
#Staxyl <- coxfunc(postcalve$Staxyl,"Staphylococcus xylosus")
Staxyl <- coxfunc(postcalve$Staxyl,
exposure="Staxyl",
exponame="Staphylococcus xylosus",
Other=c("strata(Staaur)","Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","strata(Staspp)",
"Aerococ","strata(Entero)","LactoALL","StrepALL",
"Bac","strata(Coryn)"))
Staxyl
Exposure = Staphylococcus species (i.e. species unable to be determined using MALDI-TOF)
shoen.resid(postcalve$Staspp)## chisq df p
## expo 2.56123 1 0.11
## factor(Dryfreestall) 0.15028 1 0.70
## factor(East) 0.03619 1 0.85
## factor(Parity) 1.51705 3 0.68
## factor(Breed) 0.00891 1 0.92
## factor(Bedtype) 3.77068 3 0.29
## GLOBAL 8.85043 10 0.55shoen.resid.plot(postcalve$Staspp)
#Staspp <- coxfunc(postcalve$Staspp,"Staphylococcus species")
Staspp <- coxfunc(postcalve$Staspp,
exposure="strata(Staspp)",
exponame="Staphylococcus sp.",
Other=c("strata(Staaur)","Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","strata(Staspp)",
"Aerococ","strata(Entero)","LactoALL","StrepALL",
"Bac","strata(Coryn)"))
Staspp
### Exposure = SSLO
shoen.resid(postcalve$SSLO)## chisq df p
## expo 1.4685 1 0.23
## factor(Dryfreestall) 0.1415 1 0.71
## factor(East) 0.0346 1 0.85
## factor(Parity) 1.5435 3 0.67
## factor(Breed) 0.0132 1 0.91
## factor(Bedtype) 3.6636 3 0.30
## GLOBAL 7.1714 10 0.71shoen.resid.plot(postcalve$SSLO)
#SSLO <- coxfunc(postcalve$SSLO,"SSLO")
SSLO <- coxfunc(postcalve$SSLO,
exposure="SSLO",
exponame="SSLO",
Other=c("strata(Staaur)","Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","strata(Staspp)",
"Bac","strata(Coryn)"))
SSLO
Exposure = Aerococcus species
shoen.resid(postcalve$Aerococ)## chisq df p
## expo 0.0898 1 0.76
## factor(Dryfreestall) 0.1431 1 0.71
## factor(East) 0.0331 1 0.86
## factor(Parity) 1.5616 3 0.67
## factor(Breed) 0.0128 1 0.91
## factor(Bedtype) 3.6529 3 0.30
## GLOBAL 6.0830 10 0.81shoen.resid.plot(postcalve$Aerococ)
#Aerococ <- coxfunc(postcalve$Aerococ,"Aerococcus sp.")
Aerococ <- coxfunc(postcalve$Aerococ,
exposure="Aerococ",
exponame="Aerococcus spp.",
Other=c("strata(Staaur)","Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","strata(Staspp)",
"Aerococ","strata(Entero)","LactoALL","StrepALL",
"Bac","strata(Coryn)"))
Aerococ
Exposure = Enterococcus species
shoen.resid(postcalve$Entero)## chisq df p
## expo 5.8941 1 0.015
## factor(Dryfreestall) 0.1394 1 0.709
## factor(East) 0.0364 1 0.849
## factor(Parity) 1.5214 3 0.677
## factor(Breed) 0.0131 1 0.909
## factor(Bedtype) 3.6945 3 0.296
## GLOBAL 11.7062 10 0.305shoen.resid.plot(postcalve$Entero)
#Entero <- coxfunc(postcalve$Entero,"Enterococcus sp.")
Entero <- coxfunc(postcalve$Entero,
exposure="strata(Entero)",
exponame="Enterococcus spp.",
Other=c("strata(Staaur)","Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","strata(Staspp)",
"Aerococ","strata(Entero)","LactoALL","StrepALL",
"Bac","strata(Coryn)"))
EnteroOutliers in early lactation.
Exposure = Lactococcus (all species combined)
shoen.resid(postcalve$LactoALL)## chisq df p
## expo 0.2006 1 0.65
## factor(Dryfreestall) 0.1442 1 0.70
## factor(East) 0.0343 1 0.85
## factor(Parity) 1.5379 3 0.67
## factor(Breed) 0.0134 1 0.91
## factor(Bedtype) 3.7050 3 0.30
## GLOBAL 5.9396 10 0.82shoen.resid.plot(postcalve$LactoALL)
#LactoALL <- coxfunc(postcalve$LactoALL,"Lactococcus spp.")
LactoALL <- coxfunc(postcalve$LactoALL,
exposure="LactoALL",
exponame="Lactococcus spp.",
Other=c("strata(Staaur)","Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","strata(Staspp)",
"Aerococ","strata(Entero)","LactoALL","StrepALL",
"Bac","strata(Coryn)"))
LactoALL
Exposure = Lactococcus garviae
shoen.resid(postcalve$Lactogar)## chisq df p
## expo 0.9657 1 0.33
## factor(Dryfreestall) 0.1417 1 0.71
## factor(East) 0.0338 1 0.85
## factor(Parity) 1.5504 3 0.67
## factor(Breed) 0.0125 1 0.91
## factor(Bedtype) 3.7070 3 0.29
## GLOBAL 6.4452 10 0.78shoen.resid.plot(postcalve$Lactogar)
#Lactogar <- coxfunc(postcalve$Lactogar,"Lactococcus garviae")
Lactogar <- coxfunc(postcalve$Lactogar,
exposure="Lactogar",
exponame="Lactococcus garviae",
Other=c("strata(Staaur)","Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","strata(Staspp)",
"Aerococ","strata(Entero)","Lactolac","Lactosp","StrepALL",
"Bac","strata(Coryn)"))
Lactogar
Exposure = Lactococcus lactis
shoen.resid(postcalve$Lactolac)## chisq df p
## expo 0.1032 1 0.75
## factor(Dryfreestall) 0.1467 1 0.70
## factor(East) 0.0312 1 0.86
## factor(Parity) 1.5352 3 0.67
## factor(Breed) 0.0135 1 0.91
## factor(Bedtype) 3.7028 3 0.30
## GLOBAL 6.1006 10 0.81shoen.resid.plot(postcalve$Lactolac)
#Lactolac <- coxfunc(postcalve$Lactolac,"Lactococcus lactis")
Lactolac <- coxfunc(postcalve$Lactolac,
exposure="Lactolac",
exponame="Lactococcus lactis",
Other=c("strata(Staaur)","Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","strata(Staspp)",
"Aerococ","strata(Entero)","Lactogar","StrepALL",
"Bac","strata(Coryn)"))
Lactolac
Exposure = Streptococcus
shoen.resid(postcalve$StrepALL)## chisq df p
## expo 0.5208 1 0.47
## factor(Dryfreestall) 0.1443 1 0.70
## factor(East) 0.0339 1 0.85
## factor(Parity) 1.5497 3 0.67
## factor(Breed) 0.0125 1 0.91
## factor(Bedtype) 3.7180 3 0.29
## GLOBAL 6.4715 10 0.77shoen.resid.plot(postcalve$StrepALL)
#StrepALL <- coxfunc(postcalve$StrepALL,"Streptococcus spp.")
StrepALL <- coxfunc(postcalve$StrepALL,
exposure="StrepALL",
exponame="Streptococcus spp.",
Other=c("strata(Staaur)","Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","strata(Staspp)",
"Aerococ","strata(Entero)","LactoALL","StrepALL",
"Bac","strata(Coryn)"))
StrepALL
Exposure = Streptococcus dysgalactiae
shoen.resid(postcalve$Strepdys)## chisq df p
## expo 0.1211 1 0.73
## factor(Dryfreestall) 0.1410 1 0.71
## factor(East) 0.0354 1 0.85
## factor(Parity) 1.5400 3 0.67
## factor(Breed) 0.0132 1 0.91
## factor(Bedtype) 3.6793 3 0.30
## GLOBAL 6.2013 10 0.80shoen.resid.plot(postcalve$Strepdys)
#Strepdys <- coxfunc(postcalve$Strepdys,"Streptococcus dysgalactiae")
Strepdys <- coxfunc(postcalve$Strepdys,
exposure="Strepdys",
exponame="Streptococcus dysgalactiae",
Other=c("strata(Staaur)","Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","strata(Staspp)",
"Aerococ","strata(Entero)","LactoALL","Strepub","Strepspp",
"Bac","strata(Coryn)"))
Strepdys
Exposure = Streptococcus uberis
shoen.resid(postcalve$Strepub)## chisq df p
## expo 2.6468 1 0.10
## factor(Dryfreestall) 0.1484 1 0.70
## factor(East) 0.0344 1 0.85
## factor(Parity) 1.5514 3 0.67
## factor(Breed) 0.0127 1 0.91
## factor(Bedtype) 3.7094 3 0.29
## GLOBAL 8.8461 10 0.55shoen.resid.plot(postcalve$Strepub)
#Strepub <- coxfunc(postcalve$Strepub,"Streptococcus uberis")
Strepub <- coxfunc(postcalve$Strepub,
exposure="Strepub",
exponame="Streptococcus uberis",
Other=c("strata(Staaur)","Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","strata(Staspp)",
"Aerococ","strata(Entero)","LactoALL","Strepdys","Strepspp",
"Bac","strata(Coryn)"))
Strepub
Exposure = Strep species (unable to be determined using MALDI-TOF)
shoen.resid(postcalve$Strepspp)## chisq df p
## expo 2.9264 1 0.087
## factor(Dryfreestall) 0.1419 1 0.706
## factor(East) 0.0346 1 0.852
## factor(Parity) 1.5425 3 0.672
## factor(Breed) 0.0128 1 0.910
## factor(Bedtype) 3.7056 3 0.295
## GLOBAL 8.6927 10 0.561shoen.resid.plot(postcalve$Strepspp)
#Strepspp <- coxfunc(postcalve$Strepspp,"Streptococcus sp.")
Strepspp <- coxfunc(postcalve$Strepspp,
exposure="Strepspp",
exponame="Streptococcus sp.",
Other=c("strata(Staaur)","Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","strata(Staspp)",
"Aerococ","strata(Entero)","LactoALL","Strepub","Strepdys",
"Bac","strata(Coryn)"))
Strepspp
Exposure = Bacillus
shoen.resid(postcalve$Bac)## chisq df p
## expo 0.3182 1 0.57
## factor(Dryfreestall) 0.1419 1 0.71
## factor(East) 0.0344 1 0.85
## factor(Parity) 1.5416 3 0.67
## factor(Breed) 0.0128 1 0.91
## factor(Bedtype) 3.7063 3 0.29
## GLOBAL 6.2644 10 0.79shoen.resid.plot(postcalve$Bac)
#Bac <- coxfunc(postcalve$Bac,"Bacillus sp.")
Bac <- coxfunc(postcalve$Bac,
exposure="Bac",
exponame="Bacillus spp.",
Other=c("strata(Staaur)","Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","strata(Staspp)",
"Aerococ","strata(Entero)","LactoALL","Strepub","Strepdys",
"Bac","strata(Coryn)"))
Bac
Exposure = Corynebacterium
shoen.resid(postcalve$Coryn)## chisq df p
## expo 5.6577 1 0.017
## factor(Dryfreestall) 0.1463 1 0.702
## factor(East) 0.0313 1 0.860
## factor(Parity) 1.5473 3 0.671
## factor(Breed) 0.0134 1 0.908
## factor(Bedtype) 3.7187 3 0.293
## GLOBAL 11.0959 10 0.350shoen.resid.plot(postcalve$Coryn)
#Coryn <- coxfunc(postcalve$Coryn,"Corynebacterium sp.")
Coryn <- coxfunc(postcalve$Coryn,
exposure="Coryn",
exponame="Corynebacterium spp.",
Other=c("strata(Staaur)","Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","strata(Staspp)",
"Aerococ","strata(Entero)","LactoALL","Strepub","Strepdys",
"Bac"))
CorynExposure = All IMI
shoen.resid(postcalve$IMI)## chisq df p
## expo 0.4081 1 0.52
## factor(Dryfreestall) 0.1414 1 0.71
## factor(East) 0.0347 1 0.85
## factor(Parity) 1.5544 3 0.67
## factor(Breed) 0.0140 1 0.91
## factor(Bedtype) 3.6720 3 0.30
## GLOBAL 6.3973 10 0.78shoen.resid.plot(postcalve$IMI)
#IMIALL <- coxfunc(postcalve$IMI,"All pathogens")
IMIALL <- coxfunc(postcalve$IMI,
exposure="IMI",
exponame="All pathogens",
Other=c())
IMIALL
Exposure = Major pathogens
shoen.resid(postcalve$IMIMAJ)## chisq df p
## expo 2.8381 1 0.092
## factor(Dryfreestall) 0.1412 1 0.707
## factor(East) 0.0335 1 0.855
## factor(Parity) 1.5489 3 0.671
## factor(Breed) 0.0133 1 0.908
## factor(Bedtype) 3.6540 3 0.301
## GLOBAL 8.5921 10 0.571shoen.resid.plot(postcalve$IMIMAJ)
#IMIMAJ <- coxfunc(postcalve$IMIMAJ,"Major pathogens")
IMIMAJ <- coxfunc(postcalve$IMIMAJ,
exposure="IMIMAJ",
exponame="Major pathogens",
Other=c("IMIMIN"))
IMIMAJ
Exposure = Minor pathogens (NAS, Bacillus and Corynebacterium)
shoen.resid(postcalve$IMIMIN)## chisq df p
## expo 0.8331 1 0.36
## factor(Dryfreestall) 0.1479 1 0.70
## factor(East) 0.0329 1 0.86
## factor(Parity) 1.5176 3 0.68
## factor(Breed) 0.0101 1 0.92
## factor(Bedtype) 3.7255 3 0.29
## GLOBAL 6.6532 10 0.76shoen.resid.plot(postcalve$IMIMIN)
#IMIMIN <- coxfunc(postcalve$IMIMIN,"Minor pathogens")
IMIMIN <- coxfunc(postcalve$IMIMIN,
exposure="IMIMIN",
exponame="Minor pathogens",
Other=c("IMIMAJ"))
IMIMINOutcome 2 Final models
Blank <- Staaur
Blank$Exposure <- NA
Blank$RiskExp0 <- NA
Blank$RiskExp1 <- NA
Blank$HR <- 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
## 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("Hazard ratio","\n",FP$HR))
tiff("FP2.tiff", units="in", width=8, height=6, res=300)
library(forestplot)
forestplot(labeltext=tabletext, graph.pos=4,graphwidth = unit(35,'mm'),
clip=c(0.2, 5),
xticks = c(.2, 1,5),
xlog= TRUE,
mean=c(NA,NA,FP$estimate),
lower=c(NA,NA,FP$conf.low), upper=c(NA,NA,FP$conf.high),
title="",
xlab="Hazard ratio for culling or death (Exposued/Unexposed)",
#hrzl_lines=list("3" = gpar(lwd=1, col="#000000"),
# "7" = gpar(lwd=110, 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