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
- Summary of dataset
- Outcome 3 - Milk yield (1-120 DIM)
- Model building plan
- 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 = SSLO
- 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 3: Final models
- Model building plan
Part 2: Milk yield analysis
library(knitr)
library(readr)
library(ggplot2)
library(survminer)## Loading required package: ggpubr## Loading required package: magrittrlibrary(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, unionlibrary(survival)
library(broom)
library(forestplot)## Loading required package: grid## Loading required package: checkmatelibrary(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"))
# Switching to DHI dataset
Summary of dataset
library(summarytools)
print(summarytools::dfSummary(dhi,valid.col=FALSE, graph.col=F, style="grid"))## Data Frame Summary
## dhi
## Dimensions: 7069 x 84
## Duplicates: 0
##
## +----+---------------+---------------------------------+----------------------+----------+
## | No | Variable | Stats / Values | Freqs (% of Valid) | Missing |
## +====+===============+=================================+======================+==========+
## | 1 | HerdID | Mean (sd) : 40 (23.3) | 69 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 1 < 40 < 80 | | |
## | | | IQR (CV) : 40 (0.6) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 2 | Cow | Mean (sd) : 12167.6 (12484) | 2076 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 4 < 8316 < 82342 | | |
## | | | IQR (CV) : 9087 (1) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 3 | ABX | Min : 0 | 0 : 542 ( 7.7%) | 0 |
## | | [numeric] | Mean : 0.9 | 1 : 6527 (92.3%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 4 | Calv2Date | min : 2017-08-28 | 240 distinct values | 0 |
## | | [Date] | med : 2017-12-07 | | (0%) |
## | | | max : 2018-07-15 | | |
## | | | range : 10m 17d | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 5 | Herdcow | Mean (sd) : 4346 (3068.8) | 2180 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 38 < 3897.9 < 17700.8 | | |
## | | | IQR (CV) : 4058.3 (0.7) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 6 | PeakSCC42D | Mean (sd) : 4.4 (1.4) | 465 distinct values | 403 |
## | | [numeric] | min < med < max: | | (5.7%) |
## | | | 0 < 4.4 < 9.2 | | |
## | | | IQR (CV) : 1.6 (0.3) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 7 | PeakSCC3D | Mean (sd) : 4.8 (1.4) | 596 distinct values | 201 |
## | | [numeric] | min < med < max: | | (2.84%) |
## | | | 0 < 4.8 < 9.2 | | |
## | | | IQR (CV) : 1.6 (0.3) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 8 | PeakSCCD | Mean (sd) : 5.4 (1.5) | 791 distinct values | 201 |
## | | [numeric] | min < med < max: | | (2.84%) |
## | | | 0 < 5.3 < 9.2 | | |
## | | | IQR (CV) : 1.8 (0.3) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 9 | SCCDDO | Mean (sd) : 4.4 (1.5) | 481 distinct values | 201 |
## | | [numeric] | min < med < max: | | (2.84%) |
## | | | 0 < 4.4 < 9.2 | | |
## | | | IQR (CV) : 1.6 (0.3) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 10 | MYDDO | Mean (sd) : 60.7 (21.5) | 119 distinct values | 201 |
## | | [numeric] | min < med < max: | | (2.84%) |
## | | | 0 < 61 < 179 | | |
## | | | IQR (CV) : 27 (0.4) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 11 | MYMEDDO | 1. 29270 | 23 ( 0.3%) | 201 |
## | | [character] | 2. 30520 | 23 ( 0.3%) | (2.84%) |
## | | | 3. 29160 | 21 ( 0.3%) | |
## | | | 4. 32280 | 20 ( 0.3%) | |
## | | | 5. 25850 | 19 ( 0.3%) | |
## | | | 6. 27690 | 19 ( 0.3%) | |
## | | | 7. 30970 | 19 ( 0.3%) | |
## | | | 8. 29570 | 18 ( 0.3%) | |
## | | | 9. 29500 | 17 ( 0.2%) | |
## | | | 10. 26870 | 16 ( 0.2%) | |
## | | | [ 1309 others ] | 6673 (97.2%) | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 12 | CM1Dcount | Mean (sd) : 0.3 (0.6) | 0 : 5551 (82.2%) | 315 |
## | | [numeric] | min < med < max: | 1 : 841 (12.4%) | (4.46%) |
## | | | 0 < 0 < 7 | 2 : 254 ( 3.8%) | |
## | | | IQR (CV) : 0 (2.5) | 3 : 74 ( 1.1%) | |
## | | | | 4 : 21 ( 0.3%) | |
## | | | | 5 : 7 ( 0.1%) | |
## | | | | 6 : 3 ( 0.0%) | |
## | | | | 7 : 3 ( 0.0%) | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 13 | PeakSCC42 | Mean (sd) : 4.3 (1.4) | 442 distinct values | 460 |
## | | [numeric] | min < med < max: | | (6.51%) |
## | | | 0 < 4.3 < 8.8 | | |
## | | | IQR (CV) : 1.6 (0.3) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 14 | PeakSCC3 | Mean (sd) : 4.7 (1.4) | 575 distinct values | 185 |
## | | [numeric] | min < med < max: | | (2.62%) |
## | | | 0 < 4.6 < 9.2 | | |
## | | | IQR (CV) : 1.7 (0.3) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 15 | PeakSCC | Mean (sd) : 5.3 (1.5) | 773 distinct values | 185 |
## | | [numeric] | min < med < max: | | (2.62%) |
## | | | 0 < 5.2 < 9.2 | | |
## | | | IQR (CV) : 1.9 (0.3) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 16 | SCCDO | Mean (sd) : 4.3 (1.4) | 459 distinct values | 185 |
## | | [numeric] | min < med < max: | | (2.62%) |
## | | | 0 < 4.3 < 8.8 | | |
## | | | IQR (CV) : 1.6 (0.3) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 17 | MYDO | Mean (sd) : 29.6 (10.1) | 121 distinct values | 185 |
## | | [numeric] | min < med < max: | | (2.62%) |
## | | | 0 < 29.9 < 81.2 | | |
## | | | IQR (CV) : 13.2 (0.3) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 18 | MYMEDO | Mean (sd) : 13671.5 (2566.3) | 1321 distinct values | 185 |
## | | [numeric] | min < med < max: | | (2.62%) |
## | | | 4018.8 < 13680.3 < 23863.5 | | |
## | | | IQR (CV) : 3501.7 (0.2) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 19 | CM1count | Mean (sd) : 0.2 (0.6) | 0 : 5573 (82.5%) | 315 |
## | | [numeric] | min < med < max: | 1 : 830 (12.3%) | (4.46%) |
## | | | 0 < 0 < 7 | 2 : 250 ( 3.7%) | |
## | | | IQR (CV) : 0 (2.6) | 3 : 67 ( 1.0%) | |
## | | | | 4 : 21 ( 0.3%) | |
## | | | | 5 : 10 ( 0.1%) | |
## | | | | 7 : 3 ( 0.0%) | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 20 | DryDate | min : 2017-08-14 | 146 distinct values | 93 |
## | | [Date] | med : 2017-10-12 | | (1.32%) |
## | | | max : 2018-05-02 | | |
## | | | range : 8m 18d | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 21 | Calv1Date | min : 2015-08-17 | 433 distinct values | 110 |
## | | [Date] | med : 2016-12-08 | | (1.56%) |
## | | | max : 2017-08-09 | | |
## | | | range : 1y 11m 23d | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 22 | Parity | Mean (sd) : 0.9 (0.9) | 0 : 3105 (43.9%) | 0 |
## | | [numeric] | min < med < max: | 1 : 2349 (33.2%) | (0%) |
## | | | 0 < 1 < 3 | 2 : 1026 (14.5%) | |
## | | | IQR (CV) : 1 (1.1) | 3 : 589 ( 8.3%) | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 23 | Breed | Min : 0 | 0 : 6080 (86.0%) | 0 |
## | | [numeric] | Mean : 0.1 | 1 : 989 (14.0%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 24 | Summer | Min : 0 | 0 : 3365 (47.6%) | 0 |
## | | [numeric] | Mean : 0.5 | 1 : 3704 (52.4%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 25 | Path1 | 1. Staphyloccocus chromogene | 498 (34.4%) | 5623 |
## | | [character] | 2. Bacillus sp. | 136 ( 9.4%) | (79.54%) |
## | | | 3. Aerococcus viridans | 113 ( 7.8%) | |
## | | | 4. Staphylococcus sp. | 113 ( 7.8%) | |
## | | | 5. Lactococcus garvieae | 87 ( 6.0%) | |
## | | | 6. Corynebacterium sp. | 71 ( 4.9%) | |
## | | | 7. Staphylococcus simulans | 65 ( 4.5%) | |
## | | | 8. Staphylococcus epidermidi | 33 ( 2.3%) | |
## | | | 9. Staphylococcus xylosus/sa | 33 ( 2.3%) | |
## | | | 10. Staphylococcus haemolytic | 29 ( 2.0%) | |
## | | | [ 21 others ] | 268 (18.5%) | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 26 | Path2 | 1. Staphyloccocus chromogene | 21 (19.1%) | 6959 |
## | | [character] | 2. Lactococcus garvieae | 12 (10.9%) | (98.44%) |
## | | | 3. Corynebacterium sp. | 11 (10.0%) | |
## | | | 4. Pseudomonas sp. | 10 ( 9.1%) | |
## | | | 5. Enterococcus sp. | 7 ( 6.4%) | |
## | | | 6. Staphylococcus xylosus/sa | 7 ( 6.4%) | |
## | | | 7. Gram Positive Cocci | 6 ( 5.5%) | |
## | | | 8. Bacillus sp. | 4 ( 3.6%) | |
## | | | 9. Lactococcus sp. | 4 ( 3.6%) | |
## | | | 10. Staphylococcus sciuri | 4 ( 3.6%) | |
## | | | [ 8 others ] | 24 (21.8%) | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 27 | Bedtype | Mean (sd) : 1.4 (1.2) | 0 : 2319 (32.8%) | 0 |
## | | [numeric] | min < med < max: | 1 : 1528 (21.6%) | (0%) |
## | | | 0 < 1 < 3 | 2 : 1393 (19.7%) | |
## | | | IQR (CV) : 3 (0.9) | 3 : 1829 (25.9%) | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 28 | East | Min : 0 | 0 : 2726 (38.6%) | 0 |
## | | [numeric] | Mean : 0.6 | 1 : 4343 (61.4%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 29 | Herdsize | Mean (sd) : 2228.3 (1870.1) | 60 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 235 < 1800 < 9650 | | |
## | | | IQR (CV) : 1462 (0.8) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 30 | Dryfreestall | Min : 0 | 0 : 2701 (38.2%) | 0 |
## | | [numeric] | Mean : 0.6 | 1 : 4368 (61.8%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 31 | ITS | Min : 0 | 0 : 1494 (21.6%) | 162 |
## | | [numeric] | Mean : 0.8 | 1 : 5413 (78.4%) | (2.29%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 32 | Vax | Min : 0 | 0 : 574 ( 8.1%) | 0 |
## | | [numeric] | Mean : 0.9 | 1 : 6495 (91.9%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 33 | calvedalone | Min : 0 | 0 : 4643 (65.7%) | 0 |
## | | [numeric] | Mean : 0.3 | 1 : 2426 (34.3%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 34 | Aerococ | Min : 0 | 0 : 6620 (93.7%) | 0 |
## | | [numeric] | Mean : 0.1 | 1 : 449 ( 6.3%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 35 | Bac | Min : 0 | 0 : 6650 (94.1%) | 0 |
## | | [numeric] | Mean : 0.1 | 1 : 419 ( 5.9%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 36 | Coryn | Min : 0 | 0 : 6799 (96.2%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 270 ( 3.8%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 37 | Entero | Min : 0 | 0 : 6942 (98.2%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 127 ( 1.8%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 38 | LactoALL | Min : 0 | 0 : 6713 (95.0%) | 0 |
## | | [numeric] | Mean : 0.1 | 1 : 356 ( 5.0%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 39 | Lactogar | Min : 0 | 0 : 6830 (96.6%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 239 ( 3.4%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 40 | Lactolac | Min : 0 | 0 : 6966 (98.5%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 103 ( 1.5%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 41 | Lactosp | Min : 0 | 0 : 7038 (99.6%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 31 ( 0.4%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 42 | Pseudo | Min : 0 | 0 : 7041 (99.6%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 28 ( 0.4%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 43 | Staaur | Min : 0 | 0 : 7000 (99.0%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 69 ( 1.0%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 44 | NAS | Min : 0 | 0 : 4967 (70.3%) | 0 |
## | | [numeric] | Mean : 0.3 | 1 : 2102 (29.7%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 45 | Stachr | Min : 0 | 0 : 5727 (81.0%) | 0 |
## | | [numeric] | Mean : 0.2 | 1 : 1342 (19.0%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 46 | Staepi | Min : 0 | 0 : 6938 (98.2%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 131 ( 1.8%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 47 | Stahaem | Min : 0 | 0 : 6951 (98.3%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 118 ( 1.7%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 48 | Stasci | Min : 0 | 0 : 6975 (98.7%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 94 ( 1.3%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 49 | Stasim | Min : 0 | 0 : 6930 (98.0%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 139 ( 2.0%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 50 | Staxyl | Min : 0 | 0 : 6957 (98.4%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 112 ( 1.6%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 51 | Staspp | Min : 0 | 0 : 6602 (93.4%) | 0 |
## | | [numeric] | Mean : 0.1 | 1 : 467 ( 6.6%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 52 | StrepALL | Min : 0 | 0 : 6815 (96.4%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 254 ( 3.6%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 53 | Strepdys | Min : 0 | 0 : 6996 (99.0%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 73 ( 1.0%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 54 | Strepub | Min : 0 | 0 : 6995 (99.0%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 74 ( 1.1%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 55 | Strepspp | Min : 0 | 0 : 6949 (98.3%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 120 ( 1.7%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 56 | Yeast | Min : 0 | 0 : 7022 (99.3%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 47 ( 0.7%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 57 | Coliform | Min : 0 | 0 : 7044 (99.7%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 25 ( 0.4%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 58 | MiscGP | Min : 0 | 0 : 6893 (97.5%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 176 ( 2.5%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 59 | MiscGN | Min : 0 | 0 : 6990 (98.9%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 79 ( 1.1%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 60 | Enrolldate | min : 2017-08-10 | 56 distinct values | 0 |
## | | [Date] | med : 2017-09-09 | | (0%) |
## | | | max : 2018-03-27 | | |
## | | | range : 7m 17d | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 61 | IMI | Min : 0 | 0 : 3631 (51.4%) | 0 |
## | | [numeric] | Mean : 0.5 | 1 : 3438 (48.6%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 62 | IMIMAJ | Min : 0 | 0 : 5865 (83.0%) | 0 |
## | | [numeric] | Mean : 0.2 | 1 : 1204 (17.0%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 63 | IMIMAJ2 | Min : 0 | 0 : 5643 (79.8%) | 0 |
## | | [numeric] | Mean : 0.2 | 1 : 1426 (20.2%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 64 | IMIMIN | Min : 0 | 0 : 4488 (63.5%) | 0 |
## | | [numeric] | Mean : 0.4 | 1 : 2581 (36.5%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 65 | DIMDO | Mean (sd) : 333.4 (59.4) | 276 distinct values | 160 |
## | | [numeric] | min < med < max: | | (2.26%) |
## | | | 212 < 313 < 823 | | |
## | | | IQR (CV) : 58 (0.2) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 66 | DIMEN | Mean (sd) : 319.8 (63.3) | 289 distinct values | 110 |
## | | [numeric] | min < med < max: | | (1.56%) |
## | | | 212 < 301 < 887 | | |
## | | | IQR (CV) : 61 (0.2) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 67 | RecCM3HT | Min : 0 | 0 : 6488 (96.1%) | 315 |
## | | [numeric] | Mean : 0 | 1 : 266 ( 3.9%) | (4.46%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 68 | RecCM14 | Min : 0 | 0 : 6731 (99.7%) | 315 |
## | | [numeric] | Mean : 0 | 1 : 23 ( 0.3%) | (4.46%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 69 | RecCM3HTD | Min : 0 | 0 : 6492 (96.1%) | 315 |
## | | [numeric] | Mean : 0 | 1 : 262 ( 3.9%) | (4.46%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 70 | RecCM14D | Min : 0 | 0 : 6723 (99.5%) | 315 |
## | | [numeric] | Mean : 0 | 1 : 31 ( 0.5%) | (4.46%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 71 | En2Calv | min : 3 | 103 distinct values | 0 |
## | | [difftime] | med : 69 | | (0%) |
## | | | max : 129 | | |
## | | | units : days | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 72 | Calvedlate | 1 distinct value | 0 : 7069 (100.0%) | 0 |
## | | [numeric] | | | (0%) |
## +----+---------------+---------------------------------+----------------------+----------+
## | 73 | DPlength | min : 3 | 99 distinct values | 93 |
## | | [difftime] | med : 56 | | (1.32%) |
## | | | max : 129 | | |
## | | | units : days | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 74 | CM1 | Min : 0 | 0 : 5573 (82.5%) | 315 |
## | | [numeric] | Mean : 0.2 | 1 : 1181 (17.5%) | (4.46%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 75 | CM1D | Min : 0 | 0 : 5551 (82.2%) | 315 |
## | | [numeric] | Mean : 0.2 | 1 : 1203 (17.8%) | (4.46%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 76 | RemTAR | Mean (sd) : 189.3 (19.5) | 155 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 44 < 188 < 249 | | |
## | | | IQR (CV) : 25 (0.1) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 77 | SCC | Mean (sd) : 3.8 (1.6) | 827 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 0 < 3.6 < 9.2 | | |
## | | | IQR (CV) : 2 (0.4) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 78 | MY | Mean (sd) : 45.4 (12.6) | 179 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 0 < 46.3 < 97.1 | | |
## | | | IQR (CV) : 15.4 (0.3) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 79 | TestDIM | Mean (sd) : 59.5 (33.9) | 119 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 1 < 59 < 119 | | |
## | | | IQR (CV) : 59 (0.6) | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 80 | TestDIMcat10 | 1. (0,10] | 522 ( 7.4%) | 0 |
## | | [factor] | 2. (10,20] | 627 ( 8.9%) | (0%) |
## | | | 3. (20,30] | 660 ( 9.3%) | |
## | | | 4. (30,40] | 593 ( 8.4%) | |
## | | | 5. (40,50] | 652 ( 9.2%) | |
## | | | 6. (50,60] | 598 ( 8.5%) | |
## | | | 7. (60,70] | 581 ( 8.2%) | |
## | | | 8. (70,80] | 592 ( 8.4%) | |
## | | | 9. (80,90] | 581 ( 8.2%) | |
## | | | 10. (90,100] | 582 ( 8.2%) | |
## | | | [ 2 others ] | 1081 (15.3%) | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 81 | TestDIMcat20 | 1. 0-20 | 1149 (16.2%) | 0 |
## | | [factor] | 2. 21-40 | 1253 (17.7%) | (0%) |
## | | | 3. 41-60 | 1250 (17.7%) | |
## | | | 4. 61-80 | 1173 (16.6%) | |
## | | | 5. 81-100 | 1163 (16.4%) | |
## | | | 6. 101-120 | 1081 (15.3%) | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 82 | TestNo | Mean (sd) : 2.3 (1.1) | 1 : 2181 (30.9%) | 0 |
## | | [integer] | min < med < max: | 2 : 1990 (28.1%) | (0%) |
## | | | 1 < 2 < 5 | 3 : 1718 (24.3%) | |
## | | | IQR (CV) : 2 (0.5) | 4 : 1144 (16.2%) | |
## | | | | 5 : 36 ( 0.5%) | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 83 | SCM | Min : 0 | 0 : 5972 (84.5%) | 0 |
## | | [numeric] | Mean : 0.2 | 1 : 1097 (15.5%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+
## | 84 | SSLO | Min : 0 | 0 : 5995 (84.8%) | 0 |
## | | [numeric] | Mean : 0.2 | 1 : 1074 (15.2%) | (0%) |
## | | | Max : 1 | | |
## +----+---------------+---------------------------------+----------------------+----------+Outcome 3 - Milk yield (1-120 DIM)
Model building plan
Model type: Linear mixed model. Random intercepts fitted for cow and herd (i.e. clustering of tests with cows, and cows within herds)
Step 1: Look for interactions with DIM by fitting a interaction model (Exposure*DIM) and plotting milk yield 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: Check residuals for normality and homoscedasticity.
Step 4: Report final model
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)
my.plot <- function(expo,exponame){
mm0 <- geeglm(MY ~ factor(expo)*dhi$TestDIMcat20, id=HerdID, 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(20,55))) + aes(x=TestDIMcat20, y=estimate, group=Exposure, colour=Exposure) +
labs(y = "Milk yield)", 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
}
lmmfunc <- 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) + MYDO"
OtherIMI <- c(expo,Other)
Formula <- formula(paste("MY ~ ",paste(OtherIMI, collapse=" + "),OtherCov))
mm1 <- geeglm(Formula, id=HerdID,data=dhi)
mm0 <- mm1 %>% tidy(conf.int=T) %>% select(estimate,conf.low,conf.high) %>% slice(2)
Formula2 <- formula(paste("~",exposure))
mm0$sign <- sign(mm0$estimate)
Exp0 <- paste0(dhi %>% filter(expo2==0) %>% summarize(MY = mean(MY, na.rm = TRUE)) %>% select(MY) %>% round(1) %>% sprintf(fmt ='%#.1f')," (",
dhi %>% filter(expo2==0) %>% summarize(SD = sd(MY, na.rm = TRUE)) %>% select(SD) %>% round(1) %>% sprintf(fmt ='%#.1f'),")")
Exp1 <- paste0(dhi %>% filter(expo2==1) %>% summarize(MY = mean(MY, na.rm = TRUE)) %>% select(MY) %>% round(1) %>% sprintf(fmt ='%#.1f')," (",
dhi %>% filter(expo2==1) %>% summarize(SD = sd(MY, na.rm = TRUE)) %>% select(SD) %>% round(1) %>% sprintf(fmt ='%#.1f'),")")
a <- NA %>% as.data.frame
a$MYEX0 <- Exp0
a$MYEX1 <- Exp1
a <- merge(a,mm0)
a$Exposure <- exponame
a
a$MYDiff[a$sign==-1] <- paste0("− ",substring(sprintf("%.1f",round(a$estimate[a$sign==-1],1)),2),
sep=""," (",sprintf("%.1f",round(a$conf.low[a$sign==-1],1)),", ",
sprintf("%.1f",round(a$conf.high[a$sign==-1],1)),")")
a$MYDiff[a$sign>=0] <- paste0("+ ",sprintf("%.1f",round(a$estimate[a$sign>=0],1)),
sep=""," (",sprintf("%.1f",round(a$conf.low[a$sign>=0],1)),", ",
sprintf("%.1f",round(a$conf.high[a$sign>=0],1)),")")
a %>% select(Exposure,MYEX0,MYEX1,MYDiff,estimate,conf.low,conf.high)
}
residplot <- function(expo){
mm0 <- geeglm(MY ~ factor(ABX) + factor(expo) + factor(TestDIMcat20) + factor(Dryfreestall) +
factor(East) + factor(Parity) + factor(Breed) + factor(Bedtype) +
factor(calvedalone) + factor(Summer), id=HerdID,data=dhi)
ggplot(data.frame(eta=predict(mm0,type="link"),pearson=residuals(mm0,type="pearson")),
aes(x=eta,y=pearson)) + geom_point() + theme_bw()
}
residplotqq <- function(expo){
mm0 <- geeglm(MY ~ factor(ABX) + factor(expo) + factor(TestDIMcat20) + factor(Dryfreestall) +
factor(East) + factor(Parity) + factor(Breed) + factor(Bedtype) +
factor(calvedalone) + factor(Summer), id=HerdID,data=dhi)
qqnorm(residuals(mm0))
}Exposure = Staphylococcus aureus
Run model and check residuals
my.plot(dhi$Staaur, "Staphylococcus aureus")
residplot(dhi$Staaur)
residplotqq(dhi$Staaur)
#Staaur <- lmmfunc(dhi$Staaur,"Staphylococcus aureus")
Staaur <- lmmfunc(expo="factor(Staaur)",
exposure="Staaur",
exponame="Staphylococcus aureus",
expo2 = dhi$Staaur,
Other=c("Stachr","Staepi","Stahaem","Stasci","Stasim","Staxyl","Staspp",
"Aerococ","Entero","LactoALL","StrepALL",
"Bac","Coryn"))
StaaurExposure = Non aureus Staphylococcus sp. (whole group)
Run model and check residuals
my.plot(dhi$NAS,"Non aureus Staphylococcus sp.")
residplot(dhi$NAS)
residplotqq(dhi$NAS)
#NAS <- lmmfunc(dhi$NAS,"Non aureus Staphylococcus sp.")
NAS<- lmmfunc("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
my.plot(dhi$Stachr,"Staphylococcus chromogenes")
residplot(dhi$Stachr)
residplotqq(dhi$Stachr)
Stachr <- lmmfunc(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 <- lmmfunc(dhi$Stachr,"Staphylococcus chromogenes")
Stachr
Exposure = Staphylococcus epidermidis
Run model and check residuals
my.plot(dhi$Staepi,"Staphylococcus epidermidis")
residplot(dhi$Staepi)
residplotqq(dhi$Staepi)
#Staepi <- lmmfunc(dhi$Staepi,"Staphylococcus epidermidis")
Staepi <- lmmfunc("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
my.plot(dhi$Stahaem,"Staphylococcus haemolyticus")
residplot(dhi$Stahaem)
residplotqq(dhi$Stahaem)
#Stahaem <- lmmfunc(dhi$Stahaem,"Staphylococcus haemolyticus")
Stahaem <- lmmfunc("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
my.plot(dhi$Stasci,"Staphylococcus sciuri")
residplot(dhi$Stasci)
residplotqq(dhi$Stasci)
#Stasci <- lmmfunc(dhi$Stasci,"Staphylococcus sciuri")
Stasci <- lmmfunc("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
my.plot(dhi$Stasim,"Staphylococcus simulans")
residplot(dhi$Stasim)
residplotqq(dhi$Stasim)
#Stasim <- lmmfunc(dhi$Stasim,"Staphylococcus simulans")
Stasim <- lmmfunc("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
Exposure = Staphylococcus xylosus
Run model and check residuals
my.plot(dhi$Staxyl,"Staphylococcus xylosus")
residplot(dhi$Staxyl)
residplotqq(dhi$Staxyl)
Staxyl <- lmmfunc("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 <- lmmfunc(dhi$Staxyl,"Staphylococcus xylosus")
Staxyl
Exposure = Staphylococcus species (i.e. species unable to be determined using MALDI-TOF)
Run model and check residuals
my.plot(dhi$Staspp,"Staphylococcus species")
residplot(dhi$Staspp)
residplotqq(dhi$Staspp)
#Staspp <- lmmfunc(dhi$Staspp,"Staphylococcus species")
Staspp <- lmmfunc("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
my.plot(dhi$SSLO,"SSLO")
residplot(dhi$SSLO)
residplotqq(dhi$SSLO)
#SSLO <- lmmfunc(dhi$SSLO,SSLO)
SSLO <- lmmfunc("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
my.plot(dhi$Aerococ,"Aerococcus sp.")
residplot(dhi$Aerococ)
residplotqq(dhi$Aerococ)
#Aerococ <- lmmfunc(dhi$Aerococ,"Aerococcus sp.")
Aerococ <- lmmfunc("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
my.plot(dhi$Entero,"Enterococcus sp.")
residplot(dhi$Entero)
residplotqq(dhi$Entero)
Entero <- lmmfunc("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 <- lmmfunc(dhi$Entero,"Enterococcus sp.")
Entero
Exposure = Lactococcus (all species combined)
Run model and check residuals
my.plot(dhi$LactoALL,"Lactococcus spp.")
residplot(dhi$LactoALL)
residplotqq(dhi$LactoALL)
#LactoALL <- lmmfunc(dhi$LactoALL,"Lactococcus spp.")
LactoALL <- lmmfunc("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
my.plot(dhi$Lactogar,"Lactococcus garviae")
residplot(dhi$Lactogar)
residplotqq(dhi$Lactogar)
#Lactogar <- lmmfunc(dhi$Lactogar,"Lactococcus garviae")
Lactogar <- lmmfunc("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
my.plot(dhi$Lactolac,"Lactococcus lactis")
residplot(dhi$Lactolac)
residplotqq(dhi$Lactolac)
#Lactolac <- lmmfunc(dhi$Lactolac,"Lactococcus lactis")
Lactolac <- lmmfunc("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
my.plot(dhi$StrepALL,"Streptococcus spp.")
residplot(dhi$StrepALL)
residplotqq(dhi$StrepALL)
#StrepALL <- lmmfunc(dhi$StrepALL,"Streptococcus spp.")
StrepALL <- lmmfunc("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
Exposure = Streptococcus dysgalactiae
Run model and check residuals
my.plot(dhi$Strepdys,"Streptococcus dysgalactiae")
residplot(dhi$Strepdys)
residplotqq(dhi$Strepdys)
Strepdys <- lmmfunc("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 <- lmmfunc(dhi$Strepdys,"Streptococcus dysgalactiae")
Strepdys
Exposure = Streptococcus uberis
Run model and check residuals
my.plot(dhi$Strepub,"Streptococcus uberis")
residplot(dhi$Strepub)
residplotqq(dhi$Strepub)
Strepub <- lmmfunc("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 <- lmmfunc(dhi$Strepub,"Streptococcus uberis")
Strepub
Exposure = Strep species (unable to be determined using MALDI-TOF)
Run model and check residuals
my.plot(dhi$Strepspp,"Streptococcus sp.")
residplot(dhi$Strepspp)
residplotqq(dhi$Strepspp)
#Strepspp <- lmmfunc(dhi$Strepspp,"Streptococcus sp.")
Strepspp <- lmmfunc("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
my.plot(dhi$Bac,"Bacillus sp.")
residplot(dhi$Bac)
residplotqq(dhi$Bac)
Bac <- lmmfunc("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 <- lmmfunc(dhi$Bac,"Bacillus sp.")
Bac
Exposure = Corynebacterium
Run model and check residuals
my.plot(dhi$Coryn,"Corynebacterium sp.")
residplot(dhi$Coryn)
residplotqq(dhi$Coryn)
Coryn <- lmmfunc("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 <- lmmfunc(dhi$Coryn,"Corynebacterium sp.")
Coryn
Exposure = All IMI
Run model and check residuals
dhi <- dhi %>% filter(!is.na(MYDO))
my.plot(dhi$IMI,"All pathogens")
residplot(dhi$IMI)
residplotqq(dhi$IMI)
#IMIALL <- lmmfunc(dhi$IMI,"All pathogens")
IMIALL <- lmmfunc("factor(IMI)",
exposure="IMI",
exponame="All pathogens",
expo2 = dhi$IMI,
Other=c("MYDO"))
IMIALL
Exposure = Major pathogens
Run model and check residuals
my.plot(dhi$IMIMAJ,"Major pathogens")
residplot(dhi$IMIMAJ)
residplotqq(dhi$IMIMAJ)
IMIMAJ <- lmmfunc("factor(IMIMAJ)",
exposure="IMIMAJ",
exponame="Major pathogens",
expo2 = dhi$IMIMAJ,
Other=c("IMIMIN","MYDO"))
#IMIMAJ <- lmmfunc(dhi$IMIMAJ,"Major pathogens")
IMIMAJExposure = Minor pathogens (NAS, Bacillus and Corynebacterium)
Run model and check residuals
my.plot(dhi$IMIMIN,"Minor pathogens")
residplot(dhi$IMIMIN)
residplotqq(dhi$IMIMIN)
IMIMIN <- lmmfunc("factor(IMIMIN)",
exposure="IMIMIN",
exponame="Minor pathogens",
expo2 = dhi$IMIMIN,
Other=c("IMIMAJ","MYDO"))
#IMIMIN <- lmmfunc(dhi$IMIMIN,"Minor pathogens")
IMIMIN
Outcome 3: Final models
library(Ecfun)
Blank <- Staaur
BlankBlank$Exposure <- NA
Blank$MYEX0 <- NA
Blank$MYEX1 <- NA
Blank$MYDiff <- 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$MYEX1),
c("unexposed","\n",FP$MYEX0),
c("Difference (95% CI)","\n",FP$MYDiff))
library(forestplot)
tiff("FP3.tiff", units="in", width=8, height=6, res=300)
forestplot(labeltext=tabletext, graph.pos=4,
graphwidth = unit(35,'mm'),
clip=c(-5, 0,5),
xticks = c(-5, 0,5),
mean=c(NA,NA,FP$estimate),
lower=c(NA,NA,FP$conf.low), upper=c(NA,NA,FP$conf.high),
title="",
xlab="Difference in milk yield (exposed minus 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=0, 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