External validation of four predictive algorithms used to guide selective dry cow therapy
Sam Rowe (samrowe101@gmail.com)
June 26th, 2020
- Create algorithm groups
- Section 1: Descriptive analyses
- Cow demographics
- Herd demographics
- Intramammary infections
- Kappa matrix for ALG agreement with each other
- Proportion of cows testing positive across herds
- Proportion of cows testing positive within herds
- 2 x 2 tables for each algorithm with and without CM criteria
- Kappas for ALG with and without CM criteria
- Section 2: Agreement between algorithms and lab culture
library(knitr)
library(dplyr)##
## 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(broom)
#Set WD ----
wd <- function(){setwd("~/Dropbox/R backup/2020/IMI cohort")}
wd1 <- function(){setwd("~/Dropbox/R backup/2020/Algorithm80")}
wd1()
#dataset = cow
load(file="Alg80Final1.Rdata")
load(file="Alg80Final2.Rdata")
load(file="Alg80Final3.Rdata")
knitr::opts_chunk$set(warning=FALSE, message=FALSE)
Describe dataset
print(summarytools::dfSummary(COW,valid.col=FALSE, graph.col=F, style="grid"))## Data Frame Summary
## COW
## Dimensions: 1594 x 51
## Duplicates: 0
##
## +----+-------------+--------------------------------+----------------------+---------+
## | No | Variable | Stats / Values | Freqs (% of Valid) | Missing |
## +====+=============+================================+======================+=========+
## | 1 | Herdcow | 1. 1 297.892598095354 | 1 ( 0.1%) | 0 |
## | | [character] | 2. 1 299.819945967576 | 1 ( 0.1%) | (0%) |
## | | | 3. 1 301.496268633627 | 1 ( 0.1%) | |
## | | | 4. 1 302.866307138975 | 1 ( 0.1%) | |
## | | | 5. 1 465.725240887801 | 1 ( 0.1%) | |
## | | | 6. 1 469.497603827751 | 1 ( 0.1%) | |
## | | | 7. 1 469.804214540483 | 1 ( 0.1%) | |
## | | | 8. 1 470.416836433391 | 1 ( 0.1%) | |
## | | | 9. 1 471.181493694309 | 1 ( 0.1%) | |
## | | | 10. 1 475.743628438679 | 1 ( 0.1%) | |
## | | | [ 1584 others ] | 1584 (99.4%) | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 2 | IMI | Min : 0 | 0 : 832 (52.2%) | 0 |
## | | [numeric] | Mean : 0.5 | 1 : 762 (47.8%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 3 | IMIMAJ | Min : 0 | 0 : 1335 (83.8%) | 0 |
## | | [numeric] | Mean : 0.2 | 1 : 259 (16.2%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 4 | IMIMAJ2 | Min : 0 | 0 : 1217 (76.3%) | 0 |
## | | [numeric] | Mean : 0.2 | 1 : 377 (23.6%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 5 | IMIMIN | Min : 0 | 0 : 1025 (64.3%) | 0 |
## | | [numeric] | Mean : 0.4 | 1 : 569 (35.7%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 6 | Aerococ | Min : 0 | 0 : 1498 (94.0%) | 0 |
## | | [numeric] | Mean : 0.1 | 1 : 96 ( 6.0%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 7 | Bac | Min : 0 | 0 : 1507 (94.5%) | 0 |
## | | [numeric] | Mean : 0.1 | 1 : 87 ( 5.5%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 8 | Coryn | Min : 0 | 0 : 1536 (96.4%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 58 ( 3.6%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 9 | Entero | Min : 0 | 0 : 1572 (98.6%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 22 ( 1.4%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 10 | LactoALL | Min : 0 | 0 : 1509 (94.7%) | 0 |
## | | [numeric] | Mean : 0.1 | 1 : 85 ( 5.3%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 11 | Lactogar | Min : 0 | 0 : 1535 (96.3%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 59 ( 3.7%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 12 | Lactolac | Min : 0 | 0 : 1575 (98.8%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 19 ( 1.2%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 13 | Lactosp | Min : 0 | 0 : 1584 (99.4%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 10 ( 0.6%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 14 | Pseudo | Min : 0 | 0 : 1589 (99.7%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 5 ( 0.3%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 15 | Staaur | Min : 0 | 0 : 1574 (98.8%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 20 ( 1.2%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 16 | NAS | Min : 0 | 0 : 1127 (70.7%) | 0 |
## | | [numeric] | Mean : 0.3 | 1 : 467 (29.3%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 17 | Stachr | Min : 0 | 0 : 1286 (80.7%) | 0 |
## | | [numeric] | Mean : 0.2 | 1 : 308 (19.3%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 18 | Staepi | Min : 0 | 0 : 1571 (98.6%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 23 ( 1.4%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 19 | Stahaem | Min : 0 | 0 : 1571 (98.6%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 23 ( 1.4%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 20 | Stasci | Min : 0 | 0 : 1576 (98.9%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 18 ( 1.1%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 21 | Stasim | Min : 0 | 0 : 1567 (98.3%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 27 ( 1.7%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 22 | Staxyl | Min : 0 | 0 : 1570 (98.5%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 24 ( 1.5%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 23 | Staspp | Min : 0 | 0 : 1485 (93.2%) | 0 |
## | | [numeric] | Mean : 0.1 | 1 : 109 ( 6.8%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 24 | SSLO | Min : 0 | 0 : 1360 (85.3%) | 0 |
## | | [numeric] | Mean : 0.1 | 1 : 234 (14.7%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 25 | StrepALL | Min : 0 | 0 : 1542 (96.7%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 52 ( 3.3%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 26 | Strepdys | Min : 0 | 0 : 1577 (98.9%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 17 ( 1.1%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 27 | Strepub | Min : 0 | 0 : 1580 (99.1%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 14 ( 0.9%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 28 | Strepspp | Min : 0 | 0 : 1571 (98.6%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 23 ( 1.4%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 29 | Yeast | Min : 0 | 0 : 1586 (99.5%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 8 ( 0.5%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 30 | Coliform | Min : 0 | 0 : 1591 (99.8%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 3 ( 0.2%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 31 | MiscGP | Min : 0 | 0 : 1561 (97.9%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 33 ( 2.1%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 32 | MiscGN | Min : 0 | 0 : 1576 (98.9%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 18 ( 1.1%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 33 | Cow | Mean (sd) : 13741.3 (13826.3) | 1553 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 19 < 8839.5 < 82342 | | |
## | | | IQR (CV) : 12136.8 (1) | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 34 | HerdID | Mean (sd) : 40.6 (23.2) | 56 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 1 < 40 < 80 | | |
## | | | IQR (CV) : 43 (0.6) | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 35 | Parity | 1. 1 | 672 (42.2%) | 0 |
## | | [character] | 2. 2 | 524 (32.9%) | (0%) |
## | | | 3. 3 | 253 (15.9%) | |
## | | | 4. 4 | 145 ( 9.1%) | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 36 | PeakSCC42 | Mean (sd) : 4.4 (1.3) | 411 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 0 < 4.3 < 8.6 | | |
## | | | IQR (CV) : 1.5 (0.3) | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 37 | PeakSCC | Mean (sd) : 5.4 (1.3) | 670 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 2.7 < 5.3 < 9.2 | | |
## | | | IQR (CV) : 1.8 (0.2) | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 38 | PeakSCC3 | Mean (sd) : 4.8 (1.2) | 503 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 0 < 4.7 < 9.2 | | |
## | | | IQR (CV) : 1.6 (0.3) | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 39 | SCCDO | Mean (sd) : 4.4 (1.3) | 411 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 0 < 4.3 < 8.6 | | |
## | | | IQR (CV) : 1.5 (0.3) | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 40 | MYDO | Mean (sd) : 65.4 (22.1) | 122 distinct values | 0 |
## | | [numeric] | min < med < max: | | (0%) |
## | | | 0 < 65 < 179 | | |
## | | | IQR (CV) : 29 (0.3) | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 41 | CM1count | Mean (sd) : 0.3 (0.6) | 0 : 1305 (81.9%) | 0 |
## | | [numeric] | min < med < max: | 1 : 211 (13.2%) | (0%) |
## | | | 0 < 0 < 5 | 2 : 55 ( 3.5%) | |
## | | | IQR (CV) : 0 (2.5) | 3 : 14 ( 0.9%) | |
## | | | | 4 : 6 ( 0.4%) | |
## | | | | 5 : 3 ( 0.2%) | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 42 | CM1 | Min : 0 | 0 : 1305 (81.9%) | 0 |
## | | [numeric] | Mean : 0.2 | 1 : 289 (18.1%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 43 | RecCM3HT | Min : 0 | 0 : 1534 (96.2%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 60 ( 3.8%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 44 | RecCM14 | Min : 0 | 0 : 1589 (99.7%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 5 ( 0.3%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 45 | RemALL | Min : 0 | 0 : 1338 (83.9%) | 0 |
## | | [numeric] | Mean : 0.2 | 1 : 256 (16.1%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 46 | RemLate | Min : 0 | 0 : 1580 (99.1%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 14 ( 0.9%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 47 | RemDP | Min : 0 | 0 : 1556 (98.6%) | 16 |
## | | [numeric] | Mean : 0 | 1 : 22 ( 1.4%) | (1%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 48 | Calvedlate | Min : 0 | 0 : 1585 (99.4%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 9 ( 0.6%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 49 | nocalf | Min : 0 | 0 : 1548 (97.1%) | 0 |
## | | [numeric] | Mean : 0 | 1 : 46 ( 2.9%) | (0%) |
## | | | Max : 1 | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 50 | PeakSCC42D | Mean (sd) : 4.4 (1.4) | 426 distinct values | 28 |
## | | [numeric] | min < med < max: | | (1.76%) |
## | | | 0 < 4.5 < 9.2 | | |
## | | | IQR (CV) : 1.6 (0.3) | | |
## +----+-------------+--------------------------------+----------------------+---------+
## | 51 | PeakSCC3D | Mean (sd) : 4.9 (1.2) | 520 distinct values | 16 |
## | | [numeric] | min < med < max: | | (1%) |
## | | | 0 < 4.8 < 9.2 | | |
## | | | IQR (CV) : 1.5 (0.2) | | |
## +----+-------------+--------------------------------+----------------------+---------+Create algorithm groups
New Zealand
SCC > 120 for Parity = 1 or any clinical mastitis SCC > 150 for Parity > 1 or any clinical mastitis
library(janitor)
COW$ALGNZ <- 0
COW$ALGNZ[COW$PeakSCC>=log(150) & COW$Parity>=2] <- 1
COW$ALGNZ[COW$PeakSCC>=log(120) & COW$Parity==1] <- 1
COW$ALGNZ[COW$CM1>0] <- 1
COW %>% tabyl(Parity,ALGNZ) %>%
adorn_percentages("row") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns()NZ minus clinical mastitis component
COW$ALG2NZ <- 0
COW$ALG2NZ[COW$PeakSCC>=log(120) & COW$Parity==1] <- 1
COW$ALG2NZ[COW$PeakSCC>=log(150) & COW$Parity>=2] <- 1
COW$ALG2NZ %>% tabylUnited States
≥2 cases of clinical mastitis during lactation Any SCC > 200,000 cell/ml during lactation
COW$ALGUS <- 0
COW$ALGUS[COW$PeakSCC>log(200)] <- 1
#COW$ALGUS[COW$RecCM14>=1] <- 1
COW$ALGUS[COW$CM1count>=2] <- 1
COW %>% tabyl(Parity,ALGUS) %>%
adorn_percentages("row") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns()COW$ALG2US <- 0
COW$ALG2US[COW$PeakSCC>log(200)] <- 1
#COW$ALG2US[COW$RecCM14>=1] <- 1
#COW$ALG2US[COW$CM1count>=2] <- 1
COW$ALG2US %>% tabylUnited Kingdom
≥1 cases of clinical mastitis during the period of the last 3 herd tests (~ 3 months) Any SCC > 200,000 cell/ml during the last 3 herd tests
COW$ALGUK <- 0
COW$ALGUK[COW$PeakSCC3>log(200)] <- 1
COW$ALGUK[COW$RecCM3HT>=1] <- 1
COW %>% tabyl(Parity,ALGUK) %>%
adorn_percentages("row") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns()Without CM criteria
COW$ALG2UK <- 0
COW$ALG2UK[COW$PeakSCC3>log(200)] <- 1
COW$ALG2UK %>% tabylNetherlands
Primiparous cows: > 150.000 cells/mL at last test <6 weeks prior to dry-off Multiparous cows: > 50.000 cells/mL at last test <6 weeks prior to dry-off
COW$ALGNE <- 0
COW$ALGNE[COW$Parity==1 & COW$PeakSCC42 >log(150)] <- 1
COW$ALGNE[COW$Parity>=2 & COW$PeakSCC42 >log(50)] <- 1
COW %>% tabyl(Parity,ALGNE) %>%
adorn_percentages("row") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns()ALG2 for Dutch criteria is the same (clinical mastitis is not included in algorithm)
COW$ALG2NE <- 0
COW$ALG2NE[COW$Parity==1 & COW$PeakSCC42 >log(150)] <- 1
COW$ALG2NE[COW$Parity>=2 & COW$PeakSCC42 >log(50)] <- 1
COW$ALG2NE %>% tabylSection 1: Descriptive analyses
Cow demographics
dem$Herdcow <- paste0(dem$HerdID," ",dem$Herdcow)
enlist <- COW$Herdcow
dem <- dem[dem$Herdcow %in% enlist,]
table1::table1(~ factor(Parity) + MYDO + SCCDO + factor(SCCDO > log(200)) + PeakSCC + factor(CM1) + factor(Breed), data=dem)| Overall (N=1594) |
|
|---|---|
| factor(Parity) | |
| 1 | 672 (42.2%) |
| 2 | 524 (32.9%) |
| 3 | 253 (15.9%) |
| 4 | 145 (9.1%) |
| MYDO | |
| Mean (SD) | 65.4 (22.1) |
| Median [Min, Max] | 65.0 [0, 179] |
| SCCDO | |
| Mean (SD) | 4.38 (1.31) |
| Median [Min, Max] | 4.34 [0, 8.60] |
| factor(SCCDO > log(200)) | |
| FALSE | 1258 (78.9%) |
| TRUE | 336 (21.1%) |
| PeakSCC | |
| Mean (SD) | 5.45 (1.32) |
| Median [Min, Max] | 5.28 [2.71, 9.21] |
| factor(CM1) | |
| 0 | 1305 (81.9%) |
| 1 | 289 (18.1%) |
| factor(Breed) | |
| 0 | 1388 (87.1%) |
| 1 | 206 (12.9%) |
Herd demographics
a <- dem %>% select(HerdID,calvedalone)
b <- merge(Herd,a,by="HerdID") %>% filter(Season==1)
c <- b %>% distinct(HerdID, .keep_all= TRUE) %>% select(State,Breed,ProdL,Herdsize,Lacthouse,DryHouse,Bedtype,MilkSchedual,Predip,Postdip,Forestripping,TS,Housedrysum,Housedrywin,DryPaddock,DCT,ITS,ETS,Vax,Calvingarea,calvedalone)
c$ProdL <- c$ProdL*0.453592
table1::table1(~factor(State) + factor(Breed) + ProdL + Herdsize + factor(Lacthouse) + factor(DryHouse) + factor(Bedtype) + factor(MilkSchedual) + factor(Predip) + factor(Postdip) + factor(Forestripping) + factor(TS) + factor(DryPaddock) + factor(DCT) + factor(ITS) + factor(ETS) +factor(Vax) + factor(Calvingarea) + factor(calvedalone), data=c)| Overall (N=56) |
|
|---|---|
| factor(State) | |
| California | 11 (19.6%) |
| Idaho | 5 (8.9%) |
| indiana | 1 (1.8%) |
| Michigan | 5 (8.9%) |
| minnesota | 8 (14.3%) |
| New York | 6 (10.7%) |
| oregon | 1 (1.8%) |
| texas | 2 (3.6%) |
| Wisconsin | 17 (30.4%) |
| factor(Breed) | |
| h | 48 (85.7%) |
| j | 3 (5.4%) |
| MIX | 5 (8.9%) |
| ProdL | |
| Mean (SD) | 38.2 (5.14) |
| Median [Min, Max] | 38.6 [22.7, 47.2] |
| Herdsize | |
| Mean (SD) | 2410 (2040) |
| Median [Min, Max] | 1820 [235, 9650] |
| factor(Lacthouse) | |
| Dry lot | 1 (1.8%) |
| Free stall | 55 (98.2%) |
| factor(DryHouse) | |
| Free stall | 38 (67.9%) |
| Loose - bedded pack | 2 (3.6%) |
| Loose - dry lot | 14 (25.0%) |
| Mixed | 2 (3.6%) |
| factor(Bedtype) | |
| 0 | 17 (30.4%) |
| 1 | 15 (26.8%) |
| 2 | 8 (14.3%) |
| 3 | 16 (28.6%) |
| factor(MilkSchedual) | |
| 2x | 16 (28.6%) |
| 3x | 40 (71.4%) |
| factor(Predip) | |
| n | 2 (3.6%) |
| y | 54 (96.4%) |
| factor(Postdip) | |
| y | 56 (100%) |
| factor(Forestripping) | |
| n | 7 (12.5%) |
| y | 49 (87.5%) |
| factor(TS) | |
| n | 46 (82.1%) |
| y | 10 (17.9%) |
| factor(DryPaddock) | |
| No | 54 (96.4%) |
| Yes - summer | 2 (3.6%) |
| factor(DCT) | |
| Never | 1 (1.8%) |
| Some cows --> | 6 (10.7%) |
| Yes, for all cows---> | 49 (87.5%) |
| factor(ITS) | |
| 0 | 1 (1.8%) |
| n | 14 (25.0%) |
| y | 41 (73.2%) |
| factor(ETS) | |
| 0 | 1 (1.8%) |
| n | 49 (87.5%) |
| y | 6 (10.7%) |
| factor(Vax) | |
| NO | 3 (5.4%) |
| YES | 53 (94.6%) |
| factor(Calvingarea) | |
| In a dedicated calving area (separated from all other cows) | 19 (33.9%) |
| In a dedicated calving area (with other cows that are close to calving) | 35 (62.5%) |
| With the other dry cows in the dry cow area | 2 (3.6%) |
| factor(calvedalone) | |
| 0 | 37 (66.1%) |
| 1 | 19 (33.9%) |
Intramammary infections
countexp <- function(Expo,Exponame){
a <- NA %>% as.data.frame
a$Exposure <- Exponame
a$Count <- table(COW$HerdID,Expo) %>% as.data.frame %>% filter(Expo==1) %>% select(Freq) %>% filter(Freq>0) %>% nrow
a$Herd <- a$Count
b <- a %>% select(Exposure,Herd)
a$Exposure <- Exponame
a$Count <- table(COW$HerdID,Expo) %>% as.data.frame %>% filter(Expo==1) %>% select(Freq) %>% sum
a$QTR <- a$Count
a <- a %>% select(Exposure,QTR)
a }
countexp(COW$Lactolac,"Lactococcus lactis")Staaur <- countexp(COW$Staaur,"Staphylococcus aureus")
NAS <- countexp(COW$NAS,"Non aureus Staphylococcus")
Stachr <- countexp(COW$Stachr,"Staphylococcus chromogenes")
Staepi <- countexp(COW$Staepi,"Staphylococcus epidermidis")
Stahaem <- countexp(COW$Stahaem,"Staphylococcus haemolyticus")
Stasci <- countexp(COW$Stasci,"Staphylococcus sciuri")
Stasim <- countexp(COW$Stasim,"Staphylococcus simulans")
Staxyl <- countexp(COW$Staxyl,"Staphylococcus xylosus")
Staspp <- countexp(COW$Staspp,"Staphylococcus sp.")
SSLO <- countexp(COW$SSLO,"Strep and Strep-like organisms.")
Aerococ <- countexp(COW$Aerococ,"Aerococcus spp.")
Entero <- countexp(COW$Entero,"Enterococcus spp.")
LactoALL <- countexp(COW$LactoALL,"Lactococcus spp.")
Lactolac <- countexp(COW$Lactolac,"Lactococcus lactis")
Lactogar <- countexp(COW$Lactogar,"Lactococcus garviae")
StrepALL <- countexp(COW$StrepALL,"Streptococcus spp.")
Strepdys <- countexp(COW$Strepdys,"Streptococcus dysgalactiae")
Strepub <- countexp(COW$Strepub,"Streptococcus uberis")
Strepspp <- countexp(COW$Strepspp,"Streptococcus sp.")
Bac <- countexp(COW$Bac,"Bacillus spp.")
Coryn <- countexp(COW$Coryn,"Corynebacterium spp.")
IMI <- countexp(COW$IMI,"All pathogens")
IMIMAJ <- countexp(COW$IMIMAJ,"Major pathogens")
IMIMIN <- countexp(COW$IMIMIN,"Minor pathogens")
Yeast <- countexp(COW$Yeast,"Yeast")
coli <- countexp(COW$Coliform,"Coliform")
Pseudo <- countexp(COW$Pseudo,"Pseudo")
miscgn <- countexp(COW$MiscGN,"Misc GN")
miscgp <- countexp(COW$MiscGP,"Misc GP")
table <- rbind(NAS,Stachr,Staepi,Stahaem,Stasci,Stasim,Staxyl,Staspp,
SSLO,Aerococ,Entero,LactoALL,Lactogar,Lactolac,
StrepALL,Strepdys,Strepub,Strepspp,Staaur,
Bac,Coryn,IMI,IMIMAJ,IMIMIN,coli,Pseudo,miscgn,miscgp,Yeast)
QTRatrisk <- table(COW$IMI) %>% as.data.frame %>% select(Freq) %>% sum
table$QTRtotal <- QTRatrisk
table$QTRPer <- table$QTR/table$QTRtotal
table$QTRExp <- paste0(table$QTR," (",sprintf("%.1f",round(100*table$QTRPer, 1)),")")
table <- table %>% select(Exposure,QTRExp)
tableKappa matrix for ALG agreement with each other
library(psych)
a <- cohen.kappa(table(COW$ALGNZ,COW$ALGNE)) %>% tidy %>% select(estimate) %>% slice(n())
b <- cohen.kappa(table(COW$ALGNZ,COW$ALGUK)) %>% tidy %>% select(estimate) %>% slice(n())
c <- cohen.kappa(table(COW$ALGNZ,COW$ALGUS)) %>% tidy %>% select(estimate) %>% slice(n())
d <- cohen.kappa(table(COW$ALGNE,COW$ALGUK)) %>% tidy %>% select(estimate) %>% slice(n())
e <- cohen.kappa(table(COW$ALGNE,COW$ALGUS)) %>% tidy %>% select(estimate) %>% slice(n())
f <- cohen.kappa(table(COW$ALGUK,COW$ALGUS)) %>% tidy %>% select(estimate) %>% slice(n())
Kappamatrix <- matrix(c(1,a,b,c,
a,1,d,e,
b,d,1,f,
c,e,f,1),nrow=4,ncol=4)
is.num <- sapply(Kappamatrix, is.numeric)
Kappamatrix[is.num] <- lapply(Kappamatrix[is.num], round, 2)
colnames(Kappamatrix) <- c("NZ","NETH","UK","US")
rownames(Kappamatrix) <- c("NZ","NETH","UK","US")
write.csv(Kappamatrix,"kappamatrix.csv")Proportion of cows testing positive across herds
table1::table1(~ factor(ALGNE) + factor(ALGNZ) + factor(ALGUK) + factor(ALGUS), data=COW)| Overall (N=1594) |
|
|---|---|
| factor(ALGNE) | |
| 0 | 793 (49.7%) |
| 1 | 801 (50.3%) |
| factor(ALGNZ) | |
| 0 | 598 (37.5%) |
| 1 | 996 (62.5%) |
| factor(ALGUK) | |
| 0 | 1099 (68.9%) |
| 1 | 495 (31.1%) |
| factor(ALGUS) | |
| 0 | 801 (50.3%) |
| 1 | 793 (49.7%) |
Proportion of cows testing positive within herds
tab0 <- COW %>% group_by(HerdID) %>% summarise(NE = mean(is.na(ALGNE)),
n = n())
tab1 <- COW %>% group_by(COW$HerdID) %>% summarise(NE = mean(ALGNE),
NZ = mean(ALGNZ),
UK = mean(ALGUK),
US = mean(ALGUS),
n = n())
table1::table1(~ NE + NZ + UK + US + n , data=tab1)| Overall (N=56) |
|
|---|---|
| NE | |
| Mean (SD) | 0.506 (0.168) |
| Median [Min, Max] | 0.517 [0, 0.839] |
| NZ | |
| Mean (SD) | 0.620 (0.188) |
| Median [Min, Max] | 0.611 [0.0909, 1.00] |
| UK | |
| Mean (SD) | 0.313 (0.167) |
| Median [Min, Max] | 0.316 [0, 0.677] |
| US | |
| Mean (SD) | 0.494 (0.199) |
| Median [Min, Max] | 0.500 [0, 0.921] |
| n | |
| Mean (SD) | 28.5 (7.21) |
| Median [Min, Max] | 29.5 [11.0, 38.0] |
library(tidyr)
library(ggplot2)
library(ggridges)
library(wesanderson)
theme_set(theme_ridges())
tab1$Herd <- tab1$`COW$HerdID`
tab1$`COW$HerdID` <- NULL
tab1$n <- NULL
tablong <- tab1 %>%
pivot_longer(-Herd, names_to = "Algorithm", values_to = "proportion")
tablong$Algorithm <- recode(tablong$Algorithm,
UK = "United Kingdom",
US = "United States",
NE = "Netherlands",
NZ = "New Zealand")
hist <- ggplot(tablong, aes(x = proportion, y = Algorithm)) +
xlab("Proportion of herd test-positive") +
geom_density_ridges(aes(fill = Algorithm),scale=.75) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_discrete(expand = c(0, 0)) +
coord_cartesian(clip = "off") +
theme_ridges(grid = FALSE, center_axis_labels = TRUE) +
theme(legend.position = "none",text=element_text(family="Times New Roman"),
panel.border = element_rect(colour = "black",fill=NA, size=1),
axis.text=element_text(size=12,face="bold"),axis.title=element_text(size=12,face="bold")) +
xlim(0, 1) +
scale_fill_manual(values = wes_palette("Moonrise3", n = 4))
tiff("Hist.tiff", units="in", width=8*.75, height=6*.75, res=300)
hist
dev.off()## quartz_off_screen
## 2hist2 <- ggplot(tablong, aes(x = proportion, y = Algorithm)) +
xlab("Proportion of herd test-positive") +
geom_density_ridges(aes(fill = Algorithm),scale=1.4) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_discrete(expand = c(0, 0)) +
coord_cartesian(clip = "off") +
theme_ridges(grid = FALSE, center_axis_labels = TRUE) +
theme(legend.position = "none",text=element_text(family="Times New Roman"),
panel.border = element_rect(colour = "black",fill=NA, size=1),
axis.text=element_text(size=12,face="bold"),axis.title=element_text(size=12,face="bold")) +
xlim(0, 1) +
scale_fill_manual(values = wes_palette("Moonrise3", n = 4))
tiff("Hist2.tiff", units="in", width=8*.75, height=6*.75, res=300)
hist2
dev.off()## quartz_off_screen
## 2hist
hist
2 x 2 tables for each algorithm with and without CM criteria
NZ Algorithm
COW %>% tabyl(ALGNZ,ALG2NZ)COW %>% filter(ALGNZ==1,ALG2NZ==0) %>% select(IMI) %>% table## .
## 0 1
## 21 15NZ: 36 cows with CM but low SCC. 15 / 21 had IMI at enrollment.
Netherlands Algorithm - not conducted as algorithms are identical
UK Algorithm
COW %>% tabyl(ALGUK,ALG2UK)COW %>% filter(ALGUK==1,ALG2UK==0) %>% select(IMI) %>% table## .
## 0 1
## 8 6UK: 14 cows with CM but low SCC. 6/14 had IMI at enrollment.
US Algorithm
COW %>% tabyl(ALGUS,ALG2US)COW %>% filter(ALGUS==1,ALG2US==0) %>% select(IMI) %>% table## .
## 0
## 4US: 4 cows with CM but low SCC. 0/4 had IMI at enrollment.
Kappas for ALG with and without CM criteria
library(psych)
NZ <- cohen.kappa(table(COW$ALGNZ,COW$ALG2NZ)) %>% tidy %>% select(estimate) %>% slice(n())
US <- cohen.kappa(table(COW$ALGUS,COW$ALG2US)) %>% tidy %>% select(estimate) %>% slice(n())
UK <- cohen.kappa(table(COW$ALGUK,COW$ALG2UK)) %>% tidy %>% select(estimate) %>% slice(n())
NE <- cohen.kappa(table(COW$ALGNE,COW$ALG2NE)) %>% tidy %>% select(estimate) %>% slice(n())
rbind(NE,NZ,UK,US)Section 2: Agreement between algorithms and lab culture
Calculate Kappa for each algorithm vs lab culture
All pathogens
library(psych)
library(broom)
UK <- cohen.kappa(table(COW$IMI,COW$ALGUK)) %>% tidy %>% filter(type=="weighted")
UK$type <- "UK"
US <- cohen.kappa(table(COW$IMI,COW$ALGUS)) %>% tidy %>% filter(type=="weighted")
US$type <- "US"
NE <- cohen.kappa(table(COW$IMI,COW$ALGNE)) %>% tidy %>% filter(type=="weighted")
NE$type <- "Netherlands"
NZ <- cohen.kappa(table(COW$IMI,COW$ALGNZ)) %>% tidy %>% filter(type=="weighted")
NZ$type <- "NZ"
rbind(NZ,NE,UK,US)Major pathogens
UK <- cohen.kappa(table(COW$IMIMAJ,COW$ALGUK)) %>% tidy %>% filter(type=="weighted")
UK$type <- "UK"
US <- cohen.kappa(table(COW$IMIMAJ,COW$ALGUS)) %>% tidy %>% filter(type=="weighted")
US$type <- "US"
NE <- cohen.kappa(table(COW$IMIMAJ,COW$ALGNE)) %>% tidy %>% filter(type=="weighted")
NE$type <- "Netherlands"
NZ <- cohen.kappa(table(COW$IMIMAJ,COW$ALGNZ)) %>% tidy %>% filter(type=="weighted")
NZ$type <- "NZ"
rbind(NZ,NE,UK,US)Calculating SE and SP for detecting IMI in late lactation
TestSESPgee <- function(Test,Test1,IMItype){
mm0 <- geeglm(Test ~ factor(IMItype), data = COW, id=HerdID, family = "binomial")
mm1 <- geeglm(Test ~ IMItype,data = COW, id=HerdID, family = "binomial")
SESP <- emmeans(mm0,~IMItype,type="response") %>% tidy()
SESP$IMI <- SESP$IMItype
SESP$IMItype <- NULL
SESP$SP[SESP$IMI==0] <- 1-SESP$prob[SESP$IMI==0]
SESP$SPLO[SESP$IMI==0] <- 1-SESP$asymp.UCL[SESP$IMI==0]
SESP$SPHI[SESP$IMI==0] <- 1-SESP$asymp.LCL[SESP$IMI==0]
SESP$SE[SESP$IMI==1] <- SESP$prob[SESP$IMI==1]
SESP$SELO[SESP$IMI==1] <- SESP$asymp.LCL[SESP$IMI==1]
SESP$SEHI[SESP$IMI==1] <- SESP$asymp.UCL[SESP$IMI==1]
SE <- SESP[SESP$IMI==1,] %>% select("SE","SELO","SEHI")
SP <- SESP[SESP$IMI==0,] %>% select("SP","SPLO","SPHI")
# SESP1 <- emmeans(mm1,~IMItype,type="response") %>% tidy() %>% select(prob,asymp.LCL,asymp.UCL)
# SESP1$AP <- SESP1$prob
# SESP1$APLCL <- SESP1$asymp.LCL
# SESP1$APUCL <- SESP1$asymp.UCL
# SESP1 <- SESP1 %>% select(AP,APLCL,APUCL)
# Tab <- cbind(Test1,SESP1,SE,SP)
Tab <- cbind(Test1,SE,SP)
Tab
}
TestPVgee <- function(Test,Test1,IMItype,IMIname){
mm0 <- geeglm(IMItype ~ factor(Test),data = COW, family = binomial,id=HerdID,)
mm1 <- geeglm(IMItype ~ Test,data = COW, family = binomial,id=HerdID)
SESP <- emmeans(mm0,~Test,type="response") %>% tidy()
SESP$IMI <- SESP$Test
SESP$NPV[SESP$IMI==0] <- 1-SESP$prob[SESP$IMI==0]
SESP$NPVLO[SESP$IMI==0] <- 1-SESP$asymp.UCL[SESP$IMI==0]
SESP$NPVHI[SESP$IMI==0] <- 1-SESP$asymp.LCL[SESP$IMI==0]
SESP$PPV[SESP$IMI==1] <- SESP$prob[SESP$IMI==1]
SESP$PPVLO[SESP$IMI==1] <- SESP$asymp.LCL[SESP$IMI==1]
SESP$PPVHI[SESP$IMI==1] <- SESP$asymp.UCL[SESP$IMI==1]
PPV <- SESP[SESP$IMI==1,] %>% select("PPV","PPVLO","PPVHI")
NPV <- SESP[SESP$IMI==0,] %>% select("NPV","NPVLO","NPVHI")
# SESP1 <- emmeans(mm1,~Test,type="response") %>% tidy() %>% select(prob,asymp.LCL,asymp.UCL)
# SESP1$Prev <- SESP1$prob
# SESP1 <- SESP1 %>% select(Prev)
# Tab <- cbind(Test1,SESP1,PPV,NPV)
Tab <- cbind(Test1,PPV,NPV)
Tab
}
runtestsgee <- function(IMItype,IMIname){
a <- merge(TestSESPgee(COW$ALGUS," United States",IMItype=IMItype),TestPVgee(COW$ALGUS," United States",IMItype=IMItype))
b <- merge(TestSESPgee(COW$ALGNE," Netherlands",IMItype=IMItype),TestPVgee(COW$ALGNE," Netherlands",IMItype=IMItype))
c <- merge(TestSESPgee(COW$ALGNZ," New Zealand",IMItype=IMItype),TestPVgee(COW$ALGNZ," New Zealand",IMItype=IMItype))
d <- merge(TestSESPgee(COW$ALGUK," United Kingdom",IMItype=IMItype),TestPVgee(COW$ALGUK," United Kingdom",IMItype=IMItype))
Tab1 <- rbind(b,c,d,a)
Tab1$YI <- Tab1$SE + Tab1$SP - 1
is.num <- sapply(Tab1, is.numeric)
Tab1[is.num] <- lapply(Tab1[is.num], round, 2)
Tab1$SE <- sprintf("%.2f",round(Tab1$SE, 2))
Tab1$SELO <- sprintf("%.2f",round(Tab1$SELO, 2))
Tab1$SEHI <- sprintf("%.2f",round(Tab1$SEHI, 2))
#Concatenate to put CI in the same cell
# Tab1 <- Tab1 %>% mutate(APcl = paste0(AP," (", APLCL, " - ", APUCL, ")"))
Tab1 <- Tab1 %>% mutate(SEcl = paste0(SE," (", SELO, " - ", SEHI, ")"))
Tab1 <- Tab1 %>% mutate(SPcl = paste0(SP," (", SPLO, " - ", SPHI, ")"))
Tab1 <- Tab1 %>% mutate(PPVcl = paste0(PPV," (", PPVLO, " - ", PPVHI, ")"))
Tab1 <- Tab1 %>% mutate(NPVcl = paste0(NPV," (", NPVLO, " - ", NPVHI, ")"))
# TabGEE <- Tab1 %>% select(Test1,APcl,SEcl,SPcl,YI,PPVcl,NPVcl,Prev,SE,SP)
TabGEE <- Tab1 %>% select(Test1,SEcl,SPcl,YI,PPVcl,NPVcl,SE,SP)
Tab1 <- NULL
TabGEE$Model <- "QTRGEE"
TabGEE$Dx <- IMIname
TabGEE
}
library(geepack)
library(emmeans)
library(tidyr)
NAS <- runtestsgee(COW$NAS,"Non aureus Staphylococcus spp.")
Stachr <- runtestsgee(COW$Stachr,"Staphylococcus chromogenes")
Staepi <- runtestsgee(COW$Staepi,"Staphylococcus epidermidis")
Stahaem <- runtestsgee(COW$Stahaem,"Staphylococcus haemolyticus")
Stasci <- runtestsgee(COW$Stasci,"Staphylococcus sciuri")
Stasim <- runtestsgee(COW$Stasim,"Staphylococcus simulans")
Staxylo <- runtestsgee(COW$Staxyl,"Staphylococcus xylosus")
Stasp <- runtestsgee(COW$Staspp,"Staphylococcus sp.")
SSLO<- runtestsgee(COW$SSLO,"Streptococcus and Strep-like organisms")
Aerococ <- runtestsgee(COW$Aerococ,"Aerococcus spp.")
Entero <- runtestsgee(COW$Entero,"Enterococcus spp.")
LactoALL <- runtestsgee(COW$LactoALL,"Lactococcus spp.")
Lactogar <- runtestsgee(COW$Lactogar,"Lactococcus garviae")
Lactolac <- runtestsgee(COW$Lactolac,"Lactococcus lactis")
StrepALL <- runtestsgee(COW$StrepALL,"Streptococcus spp.")
Strepub <- runtestsgee(COW$Strepub,"Streptococcus uberis")
Strepdys <- runtestsgee(COW$Strepdys,"Streptococcus dysgalactiae")
Staur <- runtestsgee(COW$Staaur,"Staphylococcus aureus")
Bac <- runtestsgee(COW$Bac,"Bacillus spp.")
Coryn <- runtestsgee(COW$Coryn,"Corynebacterium spp.")
IMI <- runtestsgee(COW$IMI,"All pathogens")
IMIMAJ <- runtestsgee(COW$IMIMAJ,"Major pathogens")
IMIMIN <- runtestsgee(COW$IMIMIN,"Minor pathogens")
table0 <- rbind(NAS,Stachr,Staepi,Stahaem,Stasci,Stasim,Staxylo,Stasp,
SSLO,Aerococ,Entero,LactoALL,Lactogar,Lactolac,StrepALL,Strepub,Strepdys,Staur,Bac,Coryn,IMI,IMIMAJ,IMIMIN)
table <- table0
table <- table %>% select(Test1,SEcl,Dx)
table <- spread(table, Test1, SEcl)
pathsetable <- table
pathsetabletabledet <- rbind(IMI,IMIMAJ)
tabledet <- tabledet %>% select(Test1,SEcl,SPcl,PPVcl,NPVcl,Dx)
tabledet <- gather(tabledet, condition, measurement, SEcl:NPVcl, factor_key=TRUE)
tabledet <- spread(tabledet, Test1, measurement)
tabledetcountexp <- function(Expo,Exponame){
a <- NA %>% as.data.frame
a$Exposure <- Exponame
a$Count <- table(COW$HerdID,Expo) %>% as.data.frame %>% filter(Expo==1) %>% select(Freq) %>% sum
a$Cow <- a$Count
a <- a %>% select(Exposure,Cow)
a
}
Staaur <- countexp(COW$Staaur,"Staphylococcus aureus")
NAS <- countexp(COW$NAS,"Non aureus Staphylococcus")
Stachr <- countexp(COW$Stachr,"Staphylococcus chromogenes")
Staepi <- countexp(COW$Staepi,"Staphylococcus epidermidis")
Stahaem <- countexp(COW$Stahaem,"Staphylococcus haemolyticus")
Stasci <- countexp(COW$Stasci,"Staphylococcus sciuri")
Stasim <- countexp(COW$Stasim,"Staphylococcus simulans")
Staxyl <- countexp(COW$Staxyl,"Staphylococcus xylosus")
Staspp <- countexp(COW$Staspp,"Staphylococcus sp.")
SSLO <- countexp(COW$SSLO,"Strep and Strep-like organisms.")
Aerococ <- countexp(COW$Aerococ,"Aerococcus spp.")
Entero <- countexp(COW$Entero,"Enterococcus spp.")
LactoALL <- countexp(COW$LactoALL,"Lactococcus spp.")
Lactolac <- countexp(COW$Lactolac,"Lactococcus lactis")
Lactogar <- countexp(COW$Lactogar,"Lactococcus garviae")
StrepALL <- countexp(COW$StrepALL,"Streptococcus spp.")
Strepdys <- countexp(COW$Strepdys,"Streptococcus dysgalactiae")
Strepub <- countexp(COW$Strepub,"Streptococcus uberis")
Bac <- countexp(COW$Bac,"Bacillus spp.")
Coryn <- countexp(COW$Coryn,"Corynebacterium spp.")
IMI <- countexp(COW$IMI,"All pathogens")
IMIMAJ <- countexp(COW$IMIMAJ,"Major pathogens")
IMIMIN <- countexp(COW$IMIMIN,"Minor pathogens")
tablecount <- rbind(NAS,Stachr,Staepi,Stahaem,Stasci,Stasim,Staxyl,Staspp,
SSLO,Aerococ,Entero,LactoALL,Lactogar,Lactolac,
StrepALL,Strepdys,Strepub,Staaur,
Bac,Coryn,IMI,IMIMAJ,IMIMIN)
cowsatrisk <- table(COW$IMI) %>% as.data.frame %>% select(Freq) %>% sum
table <- tablecount
table$Cowstotal <- cowsatrisk
table$CowPer <- table$Cow/table$Cowstotal
table$CowExp <- paste0(table$Cow," (",sprintf("%.1f",round(100*table$CowPer, 1)),")")
table <- table %>% select(Exposure,CowExp)
table