Negatively controlled clinical trial evaluating the efficiacy of dry cow therapy
Analysis by Sam Rowe (samrowe101@gmail.com)
11 of November, 2021
Acknowledgements
Author list for this study:
Author affiliations:
pre code {
white-space:
}library(knitr)
library(readr)
library(tidyverse)
library(ggplot2)
library(janitor)
load("data.Rdata")Inspect Data
print(summarytools::dfSummary(data,valid.col=FALSE, graph.magnif=0.8, style="grid"), method = "render")Data Frame Summary
data
Dimensions: 2431 x 33Duplicates: 1
| No | Variable | Stats / Values | Freqs (% of Valid) | Graph | Missing | ||||||||||||||||||||||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | id [integer] | Mean (sd) : 192158.3 (78733.7) min < med < max: 1002 < 202751 < 800096 IQR (CV) : 113662.5 (0.4) | 2430 distinct values |  |
0 (0%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 2 | calve1date [Date] | min : 2018-04-01 med : 2019-08-26 max : 2019-11-09 range : 1y 7m 8d | 245 distinct values |  |
0 (0%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 3 | tx [factor] | 1. ITS 2. ABXITS |
|
 |
0 (0%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 4 | parity [numeric] | Mean (sd) : 3.2 (1.3) min < med < max: 2 < 3 < 9 IQR (CV) : 2 (0.4) |
|
 |
0 (0%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 5 | drydate [Date] | min : 2020-06-03 med : 2020-07-03 max : 2020-07-31 range : 1m 28d | 43 distinct values |  |
0 (0%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 6 | BIRTH [Date] | min : 2009-10-06 med : 2016-01-26 max : 2017-06-06 range : 7y 8m 0d | 1116 distinct values |  |
0 (0%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 7 | calve2date [Date] | min : 2020-06-16 med : 2020-08-18 max : 2020-10-27 range : 4m 11d | 100 distinct values |  |
44 (1.81%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 8 | remdp [numeric] | Min : 0 Mean : 0 Max : 1 |
|
 |
0 (0%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 9 | remdptype [character] | 1. DIED 2. SOLD |
|
 |
2387 (98.19%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 10 | remdpdate [Date] | min : 2020-06-14 med : 2020-07-22 max : 2020-10-28 range : 4m 14d | 31 distinct values |  |
2387 (98.19%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 11 | remdpremark [character] | 1. ABORTION 2. DIGESTIV 3. DISTOCIA 4. INFECTIO 5. INJURY 6. METABOLI 7. UNKOW |
|
 |
2387 (98.19%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 12 | rem2 [numeric] | Min : 0 Mean : 0.2 Max : 1 |
|
 |
44 (1.81%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 13 | rem2type [character] | 1. DIED 2. SOLD |
|
 |
2011 (82.72%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 14 | rem2date [Date] | min : 2020-07-02 med : 2020-11-18 max : 2021-06-09 range : 11m 7d | 167 distinct values |  |
2011 (82.72%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 15 | rem2remark [character] | 1. ACCIDEN 2. ANYSICK 3. DIGESTIV 4. INFECTIO 5. INJURY 6. LAMES 7. LOWPR 8. MASTITIS 9. METABOLI 10. RPLP 11. UNKOW |
|
 |
2011 (82.72%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 16 | cm1 [numeric] | Min : 0 Mean : 0.1 Max : 1 |
|
 |
0 (0%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 17 | cm1last [Date] | min : 2018-08-28 med : 2020-01-04 max : 2020-07-23 range : 1y 10m 25d | 167 distinct values |  |
2179 (89.63%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 18 | cm1count [integer] | Mean (sd) : 1.2 (0.6) min < med < max: 1 < 1 < 5 IQR (CV) : 0 (0.5) |
|
 |
2179 (89.63%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 19 | cm2 [numeric] | Min : 0 Mean : 0.1 Max : 1 |
|
|
44 (1.81%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 20 | cm2date [Date] | min : 2020-07-23 med : 2020-10-09 max : 2021-06-05 range : 10m 13d | 127 distinct values |  |
2196 (90.33%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 21 | cm2remark [character] | 1. N016Q3.1 2. N052Q2.1 3. N074Q1.1 4. N016Q4.1 5. N050Q4.1 6. N077Q1.1 7. N087Q1.1 8. N092Q1.1 9. N014Q2.1 10. N014Q3.1 [ 152 others ] |
|
 |
2196 (90.33%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 22 | dimdo [numeric] | Mean (sd) : 320.3 (49.4) min < med < max: 259 < 307 < 838 IQR (CV) : 59 (0.2) | 217 distinct values |  |
0 (0%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 23 | paritycat [character] | 1. 2 2. 3 3. 4 4. 5+ |
|
 |
0 (0%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 24 | postcalve [difftime] | min : 225 med : 295 max : 358 units : days | 100 distinct values |  |
44 (1.81%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 25 | remdpdateout [Date] | min : 2020-06-14 med : 2020-08-18 max : 2020-10-28 range : 4m 14d | 107 distinct values |  |
0 (0%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 26 | remdptar [numeric] | Mean (sd) : 46.7 (8.6) min < med < max: 1 < 47 < 127 IQR (CV) : 8 (0.2) | 77 distinct values |  |
0 (0%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 27 | rem2dateout [Date] | min : 2020-07-02 med : 2021-06-09 max : 2021-06-09 range : 11m 7d | 167 distinct values |  |
44 (1.81%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 28 | rem2tar [numeric] | Mean (sd) : 262.3 (82.4) min < med < max: 1 < 287 < 358 IQR (CV) : 38 (0.3) | 270 distinct values |  |
44 (1.81%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 29 | cm2dateout [Date] | min : 2020-07-02 med : 2021-06-09 max : 2021-06-09 range : 11m 7d | 225 distinct values |  |
44 (1.81%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 30 | cm2tar [numeric] | Mean (sd) : 244.8 (98.4) min < med < max: 1 < 284 < 358 IQR (CV) : 43 (0.4) | 293 distinct values |  |
44 (1.81%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 31 | age [numeric] | Mean (sd) : 4.7 (1.4) min < med < max: 3.1 < 4.4 < 10.8 IQR (CV) : 2.1 (0.3) | 1120 distinct values |  |
0 (0%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 32 | drylength [numeric] | Mean (sd) : 46.9 (7.6) min < med < max: 1 < 47 < 127 IQR (CV) : 8 (0.2) | 69 distinct values |  |
44 (1.81%) | ||||||||||||||||||||||||||||||||||||||||||||||||
| 33 | drylengthcat [factor] | 1. (0,30] 2. (30,60] 3. (60,90] 4. (90,150] |
|
 |
44 (1.81%) |
Generated by summarytools 0.9.6 (R version 4.0.0)
2021-11-13
Compare treatment groups at baseline
library(table1)
data$cm1 <- as.factor(data$cm1)
label(data$paritycat) <- "Parity"
label(data$age) <- "Age at dry-off"
label(data$dimdo) <- "DIM at dry-off"
label(data$cm1) <- "Previous clinical mastitis"
label(data$drylength) <- "Dry period length"
table1::table1(~ age + paritycat + dimdo + cm1 + drylength | tx, data=data)| ITS (N=1104) |
ABXITS (N=1327) |
Overall (N=2431) |
|
|---|---|---|---|
| Age at dry-off | |||
| Mean (SD) | 4.82 (1.44) | 4.64 (1.38) | 4.72 (1.41) |
| Median [Min, Max] | 4.58 [3.13, 10.8] | 4.32 [3.14, 10.1] | 4.43 [3.13, 10.8] |
| Parity | |||
| 2 | 434 (39.3%) | 615 (46.3%) | 1049 (43.2%) |
| 3 | 240 (21.7%) | 247 (18.6%) | 487 (20.0%) |
| 4 | 264 (23.9%) | 283 (21.3%) | 547 (22.5%) |
| 5+ | 166 (15.0%) | 182 (13.7%) | 348 (14.3%) |
| DIM at dry-off | |||
| Mean (SD) | 325 (50.3) | 317 (48.4) | 320 (49.4) |
| Median [Min, Max] | 310 [260, 671] | 304 [259, 838] | 307 [259, 838] |
| Previous clinical mastitis | |||
| 0 | 975 (88.3%) | 1204 (90.7%) | 2179 (89.6%) |
| 1 | 129 (11.7%) | 123 (9.3%) | 252 (10.4%) |
| Dry period length | |||
| Mean (SD) | 47.0 (7.78) | 46.9 (7.41) | 46.9 (7.58) |
| Median [Min, Max] | 48.0 [1.00, 127] | 47.0 [5.00, 83.0] | 47.0 [1.00, 127] |
| Missing | 30 (2.7%) | 14 (1.1%) | 44 (1.8%) |
More data
data <- data %>% mutate(
dimdo3 = ntile(dimdo,3) %>% as.factor,
dplength3 = ntile(drylength,3) %>% as.factor
)
data$dimdo3 <- data$dimdo3 %>% recode(
"1" = "Low DIM",
"2" = "Mid DIM",
"3" = "High DIM")
data$dplength3 <- data$dplength3 %>% recode(
"1" = "Low dplength",
"2" = "Mid dplength",
"3" = "High dplength")
data %>% group_by(dplength3) %>% summarise(
dplength = mean(drylength,na.rm=T)
)Compare treatment groups (all variables)
table1::table1(~ factor(paritycat) + dimdo + factor(cm1) + drylength + factor(cm2) + cm2tar + factor(remdp) + remdptar + factor(rem2) + rem2tar | tx, data=data)| ITS (N=1104) |
ABXITS (N=1327) |
Overall (N=2431) |
|
|---|---|---|---|
| factor(paritycat) | |||
| 2 | 434 (39.3%) | 615 (46.3%) | 1049 (43.2%) |
| 3 | 240 (21.7%) | 247 (18.6%) | 487 (20.0%) |
| 4 | 264 (23.9%) | 283 (21.3%) | 547 (22.5%) |
| 5+ | 166 (15.0%) | 182 (13.7%) | 348 (14.3%) |
| DIM at dry-off | |||
| Mean (SD) | 325 (50.3) | 317 (48.4) | 320 (49.4) |
| Median [Min, Max] | 310 [260, 671] | 304 [259, 838] | 307 [259, 838] |
| factor(cm1) | |||
| 0 | 975 (88.3%) | 1204 (90.7%) | 2179 (89.6%) |
| 1 | 129 (11.7%) | 123 (9.3%) | 252 (10.4%) |
| Dry period length | |||
| Mean (SD) | 47.0 (7.78) | 46.9 (7.41) | 46.9 (7.58) |
| Median [Min, Max] | 48.0 [1.00, 127] | 47.0 [5.00, 83.0] | 47.0 [1.00, 127] |
| Missing | 30 (2.7%) | 14 (1.1%) | 44 (1.8%) |
| factor(cm2) | |||
| 0 | 930 (84.2%) | 1222 (92.1%) | 2152 (88.5%) |
| 1 | 144 (13.0%) | 91 (6.9%) | 235 (9.7%) |
| Missing | 30 (2.7%) | 14 (1.1%) | 44 (1.8%) |
| cm2tar | |||
| Mean (SD) | 238 (103) | 250 (93.8) | 245 (98.4) |
| Median [Min, Max] | 282 [1.00, 358] | 286 [1.00, 356] | 284 [1.00, 358] |
| Missing | 30 (2.7%) | 14 (1.1%) | 44 (1.8%) |
| factor(remdp) | |||
| 0 | 1074 (97.3%) | 1313 (98.9%) | 2387 (98.2%) |
| 1 | 30 (2.7%) | 14 (1.1%) | 44 (1.8%) |
| remdptar | |||
| Mean (SD) | 46.6 (9.25) | 46.7 (7.94) | 46.7 (8.56) |
| Median [Min, Max] | 47.0 [1.00, 127] | 47.0 [4.00, 83.0] | 47.0 [1.00, 127] |
| factor(rem2) | |||
| 0 | 881 (79.8%) | 1086 (81.8%) | 1967 (80.9%) |
| 1 | 193 (17.5%) | 227 (17.1%) | 420 (17.3%) |
| Missing | 30 (2.7%) | 14 (1.1%) | 44 (1.8%) |
| rem2tar | |||
| Mean (SD) | 262 (83.1) | 263 (81.9) | 262 (82.4) |
| Median [Min, Max] | 287 [1.00, 358] | 287 [2.00, 356] | 287 [1.00, 358] |
| Missing | 30 (2.7%) | 14 (1.1%) | 44 (1.8%) |
data %>% tabyl(drydate,tx)tab <- data %>% group_by(tx,drydate) %>%
summarise(n = n())Analysis
Extract crude frequencies
cm.freq.tx <- data %>% group_by(tx) %>% summarise(
n.cm.0 = sum(cm2==0, na.rm=T),
n.cm.1 = sum(cm2==1, na.rm=T),
p.cm.1 = n.cm.1/(n.cm.1+n.cm.0),
cm = paste0(n.cm.1, "/" ,n.cm.0 + n.cm.1, " (",100*round(p.cm.1,4),"%)")
)
cm.freq.tx$stratatype <- "All cows"
cm.freq.tx$strata <- NA
cm.freq.parity <- data %>% group_by(paritycat,tx) %>% summarise(
n.cm.0 = sum(cm2==0, na.rm=T),
n.cm.1 = sum(cm2==1, na.rm=T),
p.cm.1 = n.cm.1/(n.cm.1+n.cm.0),
cm = paste0(n.cm.1, "/" ,n.cm.0 + n.cm.1, " (",100*round(p.cm.1,4),"%)")
)
cm.freq.parity$strata <- cm.freq.parity$paritycat
cm.freq.parity$paritycat <- NULL
cm.freq.parity$stratatype <- "Parity"
cm.freq.cm1 <- data %>% group_by(cm1,tx) %>% summarise(
n.cm.0 = sum(cm2==0, na.rm=T),
n.cm.1 = sum(cm2==1, na.rm=T),
p.cm.1 = n.cm.1/(n.cm.1+n.cm.0),
cm = paste0(n.cm.1, "/" ,n.cm.0 + n.cm.1, " (",100*round(p.cm.1,4),"%)")
)
cm.freq.cm1$strata <- cm.freq.cm1$cm1
cm.freq.cm1$cm1 <- NULL
cm.freq.cm1$stratatype <- "Clinical mastitis history"
cm.freq.dp <- data %>% group_by(dplength3,tx) %>% summarise(
n.cm.0 = sum(cm2==0, na.rm=T),
n.cm.1 = sum(cm2==1, na.rm=T),
p.cm.1 = n.cm.1/(n.cm.1+n.cm.0),
cm = paste0(n.cm.1, "/" ,n.cm.0 + n.cm.1, " (",100*round(p.cm.1,4),"%)")
)
cm.freq.dp$strata <- cm.freq.dp$dplength3
cm.freq.dp$dplength3 <- NULL
cm.freq.dp$stratatype <- "Dry period length"
cm.freq <- rbind(cm.freq.tx,
cm.freq.parity,
cm.freq.cm1,
cm.freq.dp) %>% select(stratatype,strata,tx,cm)
cm.freq$risk <- cm.freq$cm
cm.freq$outcome <- "Clinical mastitis"
cm.freq$cm <- NULL
cm.freq$outcome## [1] "Clinical mastitis" "Clinical mastitis" "Clinical mastitis"
## [4] "Clinical mastitis" "Clinical mastitis" "Clinical mastitis"
## [7] "Clinical mastitis" "Clinical mastitis" "Clinical mastitis"
## [10] "Clinical mastitis" "Clinical mastitis" "Clinical mastitis"
## [13] "Clinical mastitis" "Clinical mastitis" "Clinical mastitis"
## [16] "Clinical mastitis" "Clinical mastitis" "Clinical mastitis"
## [19] "Clinical mastitis" "Clinical mastitis" "Clinical mastitis"
## [22] "Clinical mastitis"#Culling during dry period
cull.dp.tx <- data %>% group_by(tx) %>% summarise(
n.cull.0 = sum(remdp==0, na.rm=T),
n.cull.1 = sum(remdp==1, na.rm=T),
p.cull.1 = n.cull.1/(n.cull.1+n.cull.0),
cull = paste0(n.cull.1, "/" ,n.cull.0 + n.cull.1, " (",100*round(p.cull.1,4),"%)")
)
cull.dp.tx$stratatype <- "All cows"
cull.dp.tx$strata <- NA
cull.dp.parity <- data %>% group_by(tx,paritycat) %>% summarise(
n.cull.0 = sum(remdp==0, na.rm=T),
n.cull.1 = sum(remdp==1, na.rm=T),
p.cull.1 = n.cull.1/(n.cull.1+n.cull.0),
cull = paste0(n.cull.1, "/" ,n.cull.0 + n.cull.1, " (",100*round(p.cull.1,4),"%)")
)
cull.dp.parity$stratatype <- "Parity"
cull.dp.parity$strata <- cull.dp.parity$paritycat
cull.dp.parity$paritycat <- NULL
cull.dp.cm1 <- data %>% group_by(tx,cm1) %>% summarise(
n.cull.0 = sum(remdp==0, na.rm=T),
n.cull.1 = sum(remdp==1, na.rm=T),
p.cull.1 = n.cull.1/(n.cull.1+n.cull.0),
cull = paste0(n.cull.1, "/" ,n.cull.0 + n.cull.1, " (",100*round(p.cull.1,4),"%)")
)
cull.dp.cm1$stratatype <- "Clinical mastitis history"
cull.dp.cm1$strata <- cull.dp.cm1$cm1
cull.dp.cm1$cm1 <- NULL
cull.dp <- rbind(cull.dp.tx,
cull.dp.parity,
cull.dp.cm1) %>% select(stratatype,strata,tx,cull)
cull.dp$outcome <- "Removal during the dry period"
cull.dp$risk <- cull.dp$cull
cull.dp$cull <- NULL
#Culling during 1-200 DIM
cull.lact.tx <- data %>% group_by(tx) %>% summarise(
n.cull.0 = sum(rem2==0, na.rm=T),
n.cull.1 = sum(rem2==1, na.rm=T),
p.cull.1 = n.cull.1/(n.cull.1+n.cull.0),
cull = paste0(n.cull.1, "/" ,n.cull.0 + n.cull.1, " (",100*round(p.cull.1,4),"%)")
)
cull.lact.tx$stratatype <- "All cows"
cull.lact.tx$strata <- NA
#Parity strata
cull.lact.parity <- data %>% group_by(tx,paritycat) %>% summarise(
n.cull.0 = sum(rem2==0, na.rm=T),
n.cull.1 = sum(rem2==1, na.rm=T),
p.cull.1 = n.cull.1/(n.cull.1+n.cull.0),
cull = paste0(n.cull.1, "/" ,n.cull.0 + n.cull.1, " (",100*round(p.cull.1,4),"%)")
)
cull.lact.parity$stratatype <- "Parity"
cull.lact.parity$strata <- cull.lact.parity$paritycat
cull.lact.parity$paritycat <- NULL
#clinical mastitis history
cull.lact.cm1 <- data %>% group_by(tx,cm1) %>% summarise(
n.cull.0 = sum(rem2==0, na.rm=T),
n.cull.1 = sum(rem2==1, na.rm=T),
p.cull.1 = n.cull.1/(n.cull.1+n.cull.0),
cull = paste0(n.cull.1, "/" ,n.cull.0 + n.cull.1, " (",100*round(p.cull.1,4),"%)")
)
cull.lact.cm1$stratatype <- "Clinical mastitis history"
cull.lact.cm1$strata <- cull.lact.cm1$cm1
cull.lact.cm1$cm1 <- NULL
#dry period length
cull.lact.dp <- data %>% group_by(tx,dplength3) %>% summarise(
n.cull.0 = sum(rem2==0, na.rm=T),
n.cull.1 = sum(rem2==1, na.rm=T),
p.cull.1 = n.cull.1/(n.cull.1+n.cull.0),
cull = paste0(n.cull.1, "/" ,n.cull.0 + n.cull.1, " (",100*round(p.cull.1,4),"%)")
)
cull.lact.dp$stratatype <- "Dry period length"
cull.lact.dp$strata <- cull.lact.dp$dplength3
cull.lact.dp$dplength3 <- NULL
cull.lact <- rbind(cull.lact.tx,
cull.lact.parity,
cull.lact.cm1,
cull.lact.dp) %>% select(stratatype,strata,tx,cull)
cull.lact$outcome <- "Removal during 1-200 DIM"
cull.lact$risk <- cull.lact$cull
cull.lact$cull <- NULL
freq.tab <- rbind(cm.freq,
cull.lact,
cull.dp) %>% arrange(outcome,stratatype,strata)Clinical mastitis
Raw numbers
data %>% filter(!is.na(cm2)) %>%
tabyl(tx,cm2) %>%
adorn_totals("row") %>%
adorn_percentages("row") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns %>%
adorn_titleSurvival analysis
#Will limit to 200 DIM
data %>% tabyl(cm2)data$cm2[data$cm2tar >= 200] <- 0
data$cm2tar[data$cm2tar > 200] <- 200library(survival)
library(survminer)
KM <- survfit(Surv(cm2tar, cm2) ~ tx, data = data)
summary(KM)## Call: survfit(formula = Surv(cm2tar, cm2) ~ tx, data = data)
##
## 44 observations deleted due to missingness
## tx=ITS
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 1074 2 0.998 0.00132 0.996 1.000
## 2 1071 8 0.991 0.00293 0.985 0.996
## 3 1059 5 0.986 0.00359 0.979 0.993
## 4 1052 3 0.983 0.00393 0.976 0.991
## 5 1046 1 0.982 0.00404 0.974 0.990
## 6 1040 2 0.980 0.00424 0.972 0.989
## 7 1035 4 0.977 0.00463 0.968 0.986
## 8 1029 3 0.974 0.00490 0.964 0.983
## 9 1023 1 0.973 0.00499 0.963 0.983
## 10 1022 2 0.971 0.00516 0.961 0.981
## 11 1018 1 0.970 0.00524 0.960 0.980
## 13 1011 1 0.969 0.00532 0.959 0.979
## 14 1009 1 0.968 0.00540 0.957 0.979
## 15 1006 2 0.966 0.00556 0.955 0.977
## 17 1002 2 0.964 0.00571 0.953 0.975
## 18 1000 1 0.963 0.00579 0.952 0.975
## 19 997 1 0.962 0.00586 0.951 0.974
## 21 993 1 0.961 0.00594 0.950 0.973
## 22 991 1 0.960 0.00601 0.949 0.972
## 23 989 2 0.958 0.00615 0.946 0.970
## 25 983 1 0.957 0.00622 0.945 0.970
## 26 982 2 0.955 0.00636 0.943 0.968
## 27 979 1 0.954 0.00643 0.942 0.967
## 28 977 1 0.953 0.00649 0.941 0.966
## 30 974 2 0.951 0.00663 0.939 0.965
## 31 972 3 0.949 0.00682 0.935 0.962
## 32 968 2 0.947 0.00695 0.933 0.960
## 33 964 1 0.946 0.00701 0.932 0.959
## 34 961 1 0.945 0.00707 0.931 0.959
## 35 960 3 0.942 0.00725 0.928 0.956
## 37 955 1 0.941 0.00731 0.926 0.955
## 38 954 2 0.939 0.00742 0.924 0.953
## 39 952 1 0.938 0.00748 0.923 0.953
## 40 951 1 0.937 0.00754 0.922 0.952
## 41 950 2 0.935 0.00765 0.920 0.950
## 44 946 1 0.934 0.00771 0.919 0.949
## 47 943 3 0.931 0.00787 0.916 0.946
## 48 938 1 0.930 0.00792 0.914 0.945
## 50 937 2 0.928 0.00803 0.912 0.944
## 52 934 2 0.926 0.00814 0.910 0.942
## 56 928 1 0.925 0.00819 0.909 0.941
## 58 927 1 0.924 0.00824 0.908 0.940
## 62 925 1 0.923 0.00829 0.907 0.939
## 64 923 1 0.922 0.00834 0.906 0.938
## 66 922 3 0.919 0.00849 0.902 0.936
## 68 919 2 0.917 0.00859 0.900 0.934
## 73 915 1 0.916 0.00864 0.899 0.933
## 74 914 1 0.915 0.00869 0.898 0.932
## 77 912 2 0.913 0.00878 0.896 0.930
## 79 910 2 0.911 0.00888 0.894 0.928
## 80 907 1 0.910 0.00893 0.893 0.928
## 81 905 2 0.908 0.00902 0.890 0.926
## 83 902 1 0.907 0.00906 0.889 0.925
## 88 900 1 0.906 0.00911 0.888 0.924
## 89 898 1 0.905 0.00916 0.887 0.923
## 90 897 1 0.904 0.00920 0.886 0.922
## 91 896 1 0.903 0.00925 0.885 0.921
## 92 894 1 0.902 0.00929 0.884 0.920
## 95 891 1 0.901 0.00934 0.883 0.919
## 96 889 1 0.900 0.00938 0.882 0.918
## 101 886 1 0.899 0.00942 0.880 0.917
## 102 885 1 0.898 0.00947 0.879 0.916
## 105 882 1 0.897 0.00951 0.878 0.916
## 106 881 1 0.896 0.00955 0.877 0.915
## 113 876 1 0.895 0.00960 0.876 0.914
## 114 874 1 0.894 0.00964 0.875 0.913
## 115 873 1 0.893 0.00969 0.874 0.912
## 116 872 1 0.892 0.00973 0.873 0.911
## 120 867 1 0.891 0.00977 0.872 0.910
## 121 866 1 0.890 0.00981 0.871 0.909
## 134 858 1 0.889 0.00986 0.869 0.908
## 141 854 1 0.887 0.00990 0.868 0.907
## 142 853 1 0.886 0.00994 0.867 0.906
## 149 851 1 0.885 0.00999 0.866 0.905
## 161 847 1 0.884 0.01003 0.865 0.904
## 164 844 1 0.883 0.01007 0.864 0.903
## 169 843 1 0.882 0.01011 0.863 0.902
## 170 841 1 0.881 0.01016 0.862 0.901
## 172 839 1 0.880 0.01020 0.860 0.900
## 175 838 1 0.879 0.01024 0.859 0.899
## 180 834 1 0.878 0.01028 0.858 0.898
## 182 831 1 0.877 0.01032 0.857 0.897
## 191 827 1 0.876 0.01037 0.856 0.896
## 194 825 1 0.875 0.01041 0.855 0.895
## 199 823 1 0.874 0.01045 0.854 0.895
##
## tx=ABXITS
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 1313 3 0.998 0.00132 0.995 1.000
## 2 1310 3 0.995 0.00186 0.992 0.999
## 3 1301 3 0.993 0.00228 0.989 0.998
## 4 1291 3 0.991 0.00264 0.986 0.996
## 5 1283 2 0.989 0.00285 0.984 0.995
## 6 1278 3 0.987 0.00314 0.981 0.993
## 8 1266 1 0.986 0.00323 0.980 0.993
## 9 1259 1 0.985 0.00333 0.979 0.992
## 10 1252 1 0.985 0.00341 0.978 0.991
## 13 1243 1 0.984 0.00350 0.977 0.991
## 14 1241 1 0.983 0.00359 0.976 0.990
## 21 1230 1 0.982 0.00367 0.975 0.989
## 22 1228 1 0.981 0.00376 0.974 0.989
## 23 1227 1 0.981 0.00384 0.973 0.988
## 24 1225 2 0.979 0.00400 0.971 0.987
## 26 1221 1 0.978 0.00407 0.970 0.986
## 27 1219 2 0.977 0.00422 0.968 0.985
## 28 1215 1 0.976 0.00429 0.967 0.984
## 29 1214 1 0.975 0.00436 0.966 0.984
## 31 1213 1 0.974 0.00443 0.966 0.983
## 32 1212 1 0.973 0.00450 0.965 0.982
## 37 1206 1 0.973 0.00457 0.964 0.982
## 40 1202 3 0.970 0.00477 0.961 0.980
## 41 1198 3 0.968 0.00496 0.958 0.978
## 46 1193 1 0.967 0.00502 0.957 0.977
## 47 1191 1 0.966 0.00508 0.956 0.976
## 48 1188 2 0.964 0.00520 0.954 0.975
## 49 1186 1 0.964 0.00526 0.953 0.974
## 51 1185 1 0.963 0.00532 0.952 0.973
## 52 1182 1 0.962 0.00538 0.952 0.973
## 61 1177 1 0.961 0.00543 0.951 0.972
## 62 1176 1 0.960 0.00549 0.950 0.971
## 63 1175 2 0.959 0.00560 0.948 0.970
## 65 1170 1 0.958 0.00566 0.947 0.969
## 66 1169 2 0.956 0.00576 0.945 0.968
## 69 1166 1 0.956 0.00582 0.944 0.967
## 70 1165 1 0.955 0.00587 0.943 0.966
## 71 1164 1 0.954 0.00592 0.942 0.966
## 74 1163 1 0.953 0.00597 0.941 0.965
## 76 1162 1 0.952 0.00602 0.940 0.964
## 78 1161 1 0.951 0.00607 0.940 0.963
## 86 1156 1 0.951 0.00612 0.939 0.963
## 87 1153 1 0.950 0.00617 0.938 0.962
## 94 1147 1 0.949 0.00622 0.937 0.961
## 101 1140 1 0.948 0.00627 0.936 0.960
## 102 1139 1 0.947 0.00632 0.935 0.960
## 107 1136 1 0.946 0.00637 0.934 0.959
## 110 1134 1 0.946 0.00642 0.933 0.958
## 111 1133 1 0.945 0.00647 0.932 0.958
## 118 1128 2 0.943 0.00657 0.930 0.956
## 120 1125 1 0.942 0.00661 0.929 0.955
## 124 1119 1 0.941 0.00666 0.928 0.955
## 128 1117 1 0.941 0.00671 0.928 0.954
## 129 1116 1 0.940 0.00675 0.927 0.953
## 139 1111 1 0.939 0.00680 0.926 0.952
## 146 1107 1 0.938 0.00685 0.925 0.952
## 150 1105 1 0.937 0.00689 0.924 0.951
## 154 1101 1 0.936 0.00694 0.923 0.950
## 155 1100 1 0.935 0.00699 0.922 0.949
## 167 1094 1 0.935 0.00703 0.921 0.948
## 187 1084 1 0.934 0.00708 0.920 0.948
## 191 1082 1 0.933 0.00712 0.919 0.947
## 192 1081 1 0.932 0.00717 0.918 0.946
## 193 1080 1 0.931 0.00721 0.917 0.945
## 197 1077 1 0.930 0.00726 0.916 0.945ggsurvplot(KM, data = data,
pval.coord = c(0.03, 0.1),
fun = "cumhaz",xlim=c(0,200),legend.labs = c("ITS","ITS+ABX"),
title = "", pval = T, censor=F, conf.int = F,risk.table = T,
surv.plot.height = 5, tables.theme = theme_cleantable(),
ggtheme = theme_bw(base_family = "Times"), xlab="Days in milk",
palette = "Set1",risk.table.fontsize=3,break.time.by = 10)
tiff("KM.tiff", units="in", width=4, height=9, res=300)
KM## Call: survfit(formula = Surv(cm2tar, cm2) ~ tx, data = data)
##
## 44 observations deleted due to missingness
## n events median 0.95LCL 0.95UCL
## tx=ITS 1074 128 NA NA NA
## tx=ABXITS 1313 86 NA NA NAdev.off()## quartz_off_screen
## 2Cox models
library(broom)
SR <- coxph(Surv(cm2tar, cm2) ~ tx + dimdo + factor(paritycat) + cm1 + drylengthcat + cluster(drydate), data=data)
car::Anova(SR)output <- SR %>% tidy(exp=TRUE,conf.int=T) %>% select(term,estimate,conf.low,conf.high,p.value)
output$estimate <- output$estimate %>% round(2)
output$conf.low <- output$conf.low %>% round(2)
output$conf.high <- output$conf.high %>% round(2)
output$p.value <- output$p.value %>% round(4)
outputcm.main <- output %>% filter(term=="txABXITS") %>% mutate(
estimate = paste0(round(estimate,2),
" (",
round(conf.low,2),
" to ",
round(conf.high,2),
")"),
conf.low = NULL,
conf.high = NULL,
p.value = NULL)
cm.main$contrast <- "ABXITS / ITS"
cm.main$strata <- NA
cm.main$outcome <- "Clinical mastitis"
cm.main$stratatype <- "All cows"
cm.main$term <- NULL
tab.cm <- cm.mainSR <- cox.zph(SR)
SR## chisq df p
## tx 0.196 1 0.66
## dimdo 1.459 1 0.23
## factor(paritycat) 1.816 3 0.61
## cm1 0.315 1 0.57
## drylengthcat 4.069 3 0.25
## GLOBAL 8.199 9 0.51ggcoxzph(SR)
Stratified analysis
Parity
data %>% filter(!is.na(cm2)) %>%
tabyl(tx,cm2,paritycat) %>%
adorn_totals("row") %>%
adorn_percentages("row") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns %>%
adorn_title## $`2`
## cm2
## tx 0 1
## ITS 90.6% (384) 9.4% (40)
## ABXITS 94.6% (577) 5.4% (33)
## Total 92.9% (961) 7.1% (73)
##
## $`3`
## cm2
## tx 0 1
## ITS 90.2% (211) 9.8% (23)
## ABXITS 93.9% (230) 6.1% (15)
## Total 92.1% (441) 7.9% (38)
##
## $`4`
## cm2
## tx 0 1
## ITS 86.1% (223) 13.9% (36)
## ABXITS 92.5% (258) 7.5% (21)
## Total 89.4% (481) 10.6% (57)
##
## $`5+`
## cm2
## tx 0 1
## ITS 81.5% (128) 18.5% (29)
## ABXITS 90.5% (162) 9.5% (17)
## Total 86.3% (290) 13.7% (46)KM <- survfit(Surv(cm2tar, cm2) ~ tx, data = data %>% filter(paritycat=="2"))
summary(KM)## Call: survfit(formula = Surv(cm2tar, cm2) ~ tx, data = data %>% filter(paritycat ==
## "2"))
##
## 15 observations deleted due to missingness
## tx=ITS
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 424 1 0.998 0.00236 0.993 1.000
## 2 423 1 0.995 0.00333 0.989 1.000
## 3 422 3 0.988 0.00524 0.978 0.999
## 4 419 2 0.983 0.00619 0.971 0.996
## 7 415 1 0.981 0.00661 0.968 0.994
## 9 414 1 0.979 0.00701 0.965 0.993
## 13 412 1 0.976 0.00738 0.962 0.991
## 14 410 1 0.974 0.00774 0.959 0.989
## 17 407 1 0.972 0.00808 0.956 0.988
## 18 406 1 0.969 0.00841 0.953 0.986
## 21 401 1 0.967 0.00873 0.950 0.984
## 22 399 1 0.964 0.00904 0.947 0.982
## 23 398 1 0.962 0.00933 0.944 0.980
## 26 396 1 0.960 0.00962 0.941 0.979
## 30 394 2 0.955 0.01017 0.935 0.975
## 31 392 1 0.952 0.01043 0.932 0.973
## 32 390 1 0.950 0.01069 0.929 0.971
## 35 389 1 0.947 0.01093 0.926 0.969
## 40 388 1 0.945 0.01117 0.923 0.967
## 41 387 2 0.940 0.01164 0.917 0.963
## 47 384 1 0.938 0.01186 0.915 0.961
## 58 382 1 0.935 0.01208 0.912 0.959
## 62 381 1 0.933 0.01230 0.909 0.957
## 66 380 1 0.930 0.01251 0.906 0.955
## 68 379 2 0.925 0.01291 0.900 0.951
## 90 375 1 0.923 0.01311 0.897 0.949
## 102 374 1 0.920 0.01331 0.895 0.947
## 113 371 1 0.918 0.01350 0.892 0.945
## 115 370 1 0.915 0.01369 0.889 0.943
## 121 368 1 0.913 0.01388 0.886 0.941
## 141 364 1 0.910 0.01407 0.883 0.938
## 170 360 1 0.908 0.01425 0.880 0.936
## 180 358 1 0.905 0.01444 0.877 0.934
## 191 356 1 0.903 0.01462 0.875 0.932
##
## tx=ABXITS
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 2 610 1 0.998 0.00164 0.995 1.000
## 3 607 1 0.997 0.00232 0.992 1.000
## 4 604 2 0.993 0.00328 0.987 1.000
## 5 600 1 0.992 0.00367 0.985 0.999
## 6 598 2 0.988 0.00434 0.980 0.997
## 9 591 1 0.987 0.00465 0.978 0.996
## 14 583 1 0.985 0.00494 0.975 0.995
## 24 578 1 0.983 0.00521 0.973 0.994
## 28 577 1 0.982 0.00548 0.971 0.992
## 29 576 1 0.980 0.00573 0.969 0.991
## 47 573 1 0.978 0.00597 0.967 0.990
## 48 571 2 0.975 0.00642 0.962 0.987
## 52 568 1 0.973 0.00663 0.960 0.986
## 61 567 1 0.971 0.00684 0.958 0.985
## 62 566 1 0.970 0.00704 0.956 0.984
## 63 565 1 0.968 0.00723 0.954 0.982
## 71 563 1 0.966 0.00742 0.952 0.981
## 86 560 1 0.965 0.00761 0.950 0.980
## 87 558 1 0.963 0.00779 0.948 0.978
## 94 557 1 0.961 0.00796 0.946 0.977
## 102 555 1 0.959 0.00813 0.944 0.975
## 107 554 1 0.958 0.00830 0.941 0.974
## 110 553 1 0.956 0.00847 0.939 0.973
## 120 549 1 0.954 0.00863 0.937 0.971
## 124 547 1 0.952 0.00879 0.935 0.970
## 155 540 1 0.951 0.00894 0.933 0.968
## 167 537 1 0.949 0.00910 0.931 0.967
## 191 534 1 0.947 0.00926 0.929 0.965
## 193 533 1 0.945 0.00941 0.927 0.964
## 197 532 1 0.944 0.00956 0.925 0.962KM1 <- ggsurvplot(KM, data = data,
pval.coord = c(0.03, 0.1),
fun = "cumhaz",xlim=c(0,200),ylim=c(0,0.25),legend.labs = c("ITS","ITS+ABX"),
title = "Parity 2", pval = T, censor=F, conf.int = F,risk.table = F,
surv.plot.height = 5, tables.theme = theme_cleantable(),
ggtheme = theme_bw(base_family = "Times"), xlab="Days in milk",
palette = "Set1",risk.table.fontsize=3,break.time.by = 10)
KM <- survfit(Surv(cm2tar, cm2) ~ tx, data = data %>% filter(paritycat=="3"))
KM2 <- ggsurvplot(KM, data = data,
pval.coord = c(0.03, 0.1),
fun = "cumhaz",xlim=c(0,200),ylim=c(0,0.25),legend.labs = c("ITS","ITS+ABX"),
title = "Parity 3", pval = T, censor=F, conf.int = F,risk.table = F,
surv.plot.height = 5, tables.theme = theme_cleantable(),
ggtheme = theme_bw(base_family = "Times"), xlab="Days in milk",
palette = "Set1",risk.table.fontsize=3,break.time.by = 10)
KM <- survfit(Surv(cm2tar, cm2) ~ tx, data = data %>% filter(paritycat=="4"))
KM3 <- ggsurvplot(KM, data = data,
pval.coord = c(0.03, 0.1),
fun = "cumhaz",xlim=c(0,200),ylim=c(0,0.25),legend.labs = c("ITS","ITS+ABX"),
title = "Parity 4", pval = T, censor=F, conf.int = F,risk.table = F,
surv.plot.height = 5, tables.theme = theme_cleantable(),
ggtheme = theme_bw(base_family = "Times"), xlab="Days in milk",
palette = "Set1",risk.table.fontsize=3,break.time.by = 10)
KM <- survfit(Surv(cm2tar, cm2) ~ tx, data = data %>% filter(paritycat=="5+"))
KM4 <- ggsurvplot(KM, data = data,
pval.coord = c(0.03, 0.1),
fun = "cumhaz",xlim=c(0,200),ylim=c(0,0.25),legend.labs = c("ITS","ITS+ABX"),
title = "Parity 5+", pval = T, censor=F, conf.int = F,risk.table = F,
surv.plot.height = 5, tables.theme = theme_cleantable(),
ggtheme = theme_bw(base_family = "Times"), xlab="Days in milk",
palette = "Set1",risk.table.fontsize=3,break.time.by = 10)
splots <- list(KM1,KM2,KM3,KM4)
arrange_ggsurvplots(splots,ncol=2,nrow=2)
library(broom)
library(emmeans)
data$tx <- fct_relevel(data$tx,"ABXITS","ITS")
CPH <- coxph(Surv(cm2tar, cm2) ~ tx*factor(paritycat) + dimdo + cm1 + drylengthcat + cluster(drydate), data=data)
car::Anova(CPH)cm.parity.p <- car::Anova(CPH) %>% tail(1)
out <- emmeans(CPH,pairwise ~ tx | factor(paritycat),adjust = "none",type = "response") %>% confint
out## $emmeans
## paritycat = 2:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 0.115 0.0764 Inf 0.0311 0.423
## ITS 0.205 0.1535 Inf 0.0475 0.888
##
## paritycat = 3:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 0.125 0.0846 Inf 0.0335 0.470
## ITS 0.206 0.1515 Inf 0.0486 0.871
##
## paritycat = 4:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 0.157 0.1073 Inf 0.0409 0.600
## ITS 0.283 0.2041 Inf 0.0686 1.164
##
## paritycat = 5+:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 0.188 0.1295 Inf 0.0485 0.726
## ITS 0.367 0.2483 Inf 0.0977 1.381
##
## Results are averaged over the levels of: cm1, drylengthcat
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## paritycat = 2:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 0.559 0.120 Inf 0.366 0.852
##
## paritycat = 3:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 0.610 0.231 Inf 0.291 1.280
##
## paritycat = 4:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 0.554 0.188 Inf 0.285 1.076
##
## paritycat = 5+:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 0.511 0.124 Inf 0.318 0.820
##
## Results are averaged over the levels of: cm1, drylengthcat
## Confidence level used: 0.95
## Intervals are back-transformed from the log scaletab <- data.frame(
outcome = "Clinical mastitis",
stratatype = "Parity",
strata = out[["contrasts"]][["paritycat"]],
contrast = out[["contrasts"]][["contrast"]],
est = out[["contrasts"]][["ratio"]],
lcl = out[["contrasts"]][["asymp.LCL"]],
ucl = out[["contrasts"]][["asymp.UCL"]]
)
tab <- tab %>% mutate(
estimate = paste0(round(est,2),
" (",
round(lcl,2),
" to ",
round(ucl,2),
")"),
est = NULL,
lcl = NULL,
ucl = NULL
)
tab.cm.parity <- tab
tab.cm.paritydata$tx <- fct_relevel(data$tx,"ITS","ABXITS")Clinical mastitis history
data %>% filter(!is.na(cm2)) %>%
tabyl(tx,cm2,cm1) %>%
adorn_totals("row") %>%
adorn_percentages("row") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns %>%
adorn_title## $`0`
## cm2
## tx 0 1
## ITS 89.9% (855) 10.1% (96)
## ABXITS 94.2% (1121) 5.8% (69)
## Total 92.3% (1976) 7.7% (165)
##
## $`1`
## cm2
## tx 0 1
## ITS 74.0% (91) 26.0% (32)
## ABXITS 86.2% (106) 13.8% (17)
## Total 80.1% (197) 19.9% (49)KM <- survfit(Surv(cm2tar, cm2) ~ tx, data = data %>% filter(cm1==0))
summary(KM)## Call: survfit(formula = Surv(cm2tar, cm2) ~ tx, data = data %>% filter(cm1 ==
## 0))
##
## 38 observations deleted due to missingness
## tx=ITS
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 2 950 4 0.996 0.00210 0.992 1.000
## 3 945 3 0.993 0.00278 0.987 0.998
## 4 940 2 0.991 0.00315 0.984 0.997
## 5 937 1 0.989 0.00332 0.983 0.996
## 6 932 2 0.987 0.00363 0.980 0.994
## 7 928 4 0.983 0.00419 0.975 0.991
## 8 922 2 0.981 0.00445 0.972 0.990
## 9 918 1 0.980 0.00457 0.971 0.989
## 10 917 1 0.979 0.00469 0.970 0.988
## 11 915 1 0.978 0.00480 0.968 0.987
## 13 909 1 0.977 0.00492 0.967 0.986
## 14 907 1 0.976 0.00503 0.966 0.985
## 15 905 2 0.973 0.00524 0.963 0.984
## 17 901 2 0.971 0.00545 0.961 0.982
## 18 899 1 0.970 0.00555 0.959 0.981
## 19 897 1 0.969 0.00565 0.958 0.980
## 21 894 1 0.968 0.00574 0.957 0.979
## 22 892 1 0.967 0.00584 0.956 0.978
## 23 891 2 0.965 0.00602 0.953 0.977
## 25 885 1 0.964 0.00612 0.952 0.976
## 26 884 1 0.963 0.00621 0.951 0.975
## 28 881 1 0.961 0.00629 0.949 0.974
## 30 878 2 0.959 0.00647 0.947 0.972
## 31 876 3 0.956 0.00672 0.943 0.969
## 32 872 2 0.954 0.00688 0.940 0.967
## 34 867 1 0.953 0.00696 0.939 0.966
## 35 866 3 0.949 0.00719 0.935 0.964
## 37 862 1 0.948 0.00727 0.934 0.963
## 38 861 2 0.946 0.00741 0.932 0.961
## 39 859 1 0.945 0.00749 0.930 0.960
## 40 858 1 0.944 0.00756 0.929 0.959
## 41 857 2 0.942 0.00770 0.927 0.957
## 44 853 1 0.941 0.00777 0.926 0.956
## 47 851 3 0.937 0.00797 0.922 0.953
## 48 846 1 0.936 0.00804 0.921 0.952
## 52 844 1 0.935 0.00811 0.919 0.951
## 56 839 1 0.934 0.00817 0.918 0.950
## 58 838 1 0.933 0.00824 0.917 0.949
## 62 836 1 0.932 0.00831 0.916 0.948
## 64 835 1 0.931 0.00837 0.914 0.947
## 66 834 2 0.928 0.00850 0.912 0.945
## 68 832 1 0.927 0.00856 0.911 0.944
## 73 830 1 0.926 0.00862 0.909 0.943
## 74 829 1 0.925 0.00868 0.908 0.942
## 77 827 1 0.924 0.00875 0.907 0.941
## 79 826 1 0.923 0.00881 0.906 0.940
## 81 823 2 0.921 0.00893 0.903 0.938
## 83 820 1 0.919 0.00899 0.902 0.937
## 88 818 1 0.918 0.00905 0.901 0.936
## 89 816 1 0.917 0.00910 0.900 0.935
## 90 815 1 0.916 0.00916 0.898 0.934
## 91 814 1 0.915 0.00922 0.897 0.933
## 95 810 1 0.914 0.00928 0.896 0.932
## 101 806 1 0.913 0.00933 0.895 0.931
## 105 803 1 0.912 0.00939 0.893 0.930
## 113 798 1 0.910 0.00945 0.892 0.929
## 116 796 1 0.909 0.00951 0.891 0.928
## 120 791 1 0.908 0.00956 0.890 0.927
## 121 790 1 0.907 0.00962 0.888 0.926
## 141 781 1 0.906 0.00968 0.887 0.925
## 142 780 1 0.905 0.00974 0.886 0.924
## 149 778 1 0.903 0.00979 0.884 0.923
## 161 775 1 0.902 0.00985 0.883 0.922
## 170 771 1 0.901 0.00990 0.882 0.921
## 172 769 1 0.900 0.00996 0.881 0.920
## 175 768 1 0.899 0.01002 0.879 0.919
## 180 764 1 0.898 0.01007 0.878 0.918
## 182 761 1 0.896 0.01013 0.877 0.917
## 191 757 1 0.895 0.01018 0.876 0.915
## 199 754 1 0.894 0.01024 0.874 0.914
##
## tx=ABXITS
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 1190 2 0.998 0.00119 0.996 1.000
## 2 1188 3 0.996 0.00188 0.992 0.999
## 3 1181 1 0.995 0.00205 0.991 0.999
## 4 1175 3 0.992 0.00252 0.987 0.997
## 5 1167 2 0.991 0.00279 0.985 0.996
## 6 1162 2 0.989 0.00303 0.983 0.995
## 9 1148 1 0.988 0.00315 0.982 0.994
## 10 1143 1 0.987 0.00326 0.981 0.994
## 13 1135 1 0.986 0.00337 0.980 0.993
## 14 1133 1 0.986 0.00348 0.979 0.992
## 21 1122 1 0.985 0.00359 0.978 0.992
## 22 1120 1 0.984 0.00369 0.977 0.991
## 23 1119 1 0.983 0.00379 0.976 0.990
## 24 1117 2 0.981 0.00398 0.973 0.989
## 26 1113 1 0.980 0.00408 0.972 0.988
## 27 1111 2 0.978 0.00425 0.970 0.987
## 29 1107 1 0.978 0.00434 0.969 0.986
## 31 1106 1 0.977 0.00443 0.968 0.985
## 37 1102 1 0.976 0.00451 0.967 0.985
## 40 1098 3 0.973 0.00475 0.964 0.983
## 41 1094 2 0.971 0.00491 0.962 0.981
## 46 1090 1 0.971 0.00498 0.961 0.980
## 47 1088 1 0.970 0.00506 0.960 0.980
## 48 1085 2 0.968 0.00521 0.958 0.978
## 49 1083 1 0.967 0.00528 0.957 0.977
## 51 1082 1 0.966 0.00535 0.956 0.977
## 52 1079 1 0.965 0.00542 0.955 0.976
## 61 1075 1 0.964 0.00549 0.954 0.975
## 62 1074 1 0.963 0.00555 0.953 0.974
## 63 1073 2 0.962 0.00569 0.950 0.973
## 65 1068 1 0.961 0.00575 0.949 0.972
## 66 1067 1 0.960 0.00582 0.948 0.971
## 70 1065 1 0.959 0.00588 0.947 0.970
## 74 1064 1 0.958 0.00594 0.946 0.970
## 76 1063 1 0.957 0.00601 0.945 0.969
## 86 1058 1 0.956 0.00607 0.944 0.968
## 87 1056 1 0.955 0.00613 0.943 0.967
## 94 1051 1 0.954 0.00619 0.942 0.967
## 101 1045 1 0.953 0.00625 0.941 0.966
## 102 1044 1 0.953 0.00631 0.940 0.965
## 107 1042 1 0.952 0.00637 0.939 0.964
## 111 1040 1 0.951 0.00643 0.938 0.963
## 118 1035 1 0.950 0.00649 0.937 0.963
## 124 1030 1 0.949 0.00655 0.936 0.962
## 128 1028 1 0.948 0.00661 0.935 0.961
## 146 1021 1 0.947 0.00667 0.934 0.960
## 150 1019 1 0.946 0.00672 0.933 0.959
## 154 1016 1 0.945 0.00678 0.932 0.959
## 155 1015 1 0.944 0.00684 0.931 0.958
## 167 1009 1 0.943 0.00689 0.930 0.957
## 187 1000 1 0.942 0.00695 0.929 0.956
## 191 998 1 0.941 0.00701 0.928 0.955
## 192 997 1 0.940 0.00707 0.927 0.954
## 193 996 1 0.939 0.00712 0.926 0.954
## 197 993 1 0.939 0.00718 0.925 0.953KM1 <- ggsurvplot(KM, data = data,
pval.coord = c(0.03, 0.1),
fun = "cumhaz",xlim=c(0,200),ylim=c(0,0.25),legend.labs = c("ITS","ITS+ABX"),
title = "Cows with no history of clinical mastitis", pval = T, censor=F, conf.int = F,risk.table = F,
surv.plot.height = 5, tables.theme = theme_cleantable(),
ggtheme = theme_bw(base_family = "Times"), xlab="Days in milk",
palette = "Set1",risk.table.fontsize=3,break.time.by = 10)
KM <- survfit(Surv(cm2tar, cm2) ~ tx, data = data %>% filter(cm1==1))
summary(KM)## Call: survfit(formula = Surv(cm2tar, cm2) ~ tx, data = data %>% filter(cm1 ==
## 1))
##
## 6 observations deleted due to missingness
## tx=ITS
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 123 2 0.984 0.0114 0.962 1.000
## 2 121 4 0.951 0.0194 0.914 0.990
## 3 114 2 0.935 0.0224 0.892 0.979
## 4 112 1 0.926 0.0237 0.881 0.974
## 8 107 1 0.918 0.0250 0.870 0.968
## 10 105 1 0.909 0.0262 0.859 0.962
## 26 98 1 0.900 0.0276 0.847 0.955
## 27 97 1 0.890 0.0288 0.836 0.949
## 33 95 1 0.881 0.0300 0.824 0.942
## 50 92 2 0.862 0.0322 0.801 0.927
## 52 90 1 0.852 0.0333 0.789 0.920
## 66 88 1 0.842 0.0343 0.778 0.912
## 68 87 1 0.833 0.0352 0.767 0.905
## 77 85 1 0.823 0.0361 0.755 0.897
## 79 84 1 0.813 0.0370 0.744 0.889
## 80 83 1 0.803 0.0378 0.733 0.881
## 92 82 1 0.794 0.0386 0.721 0.873
## 96 81 1 0.784 0.0394 0.710 0.865
## 102 80 1 0.774 0.0401 0.699 0.857
## 106 79 1 0.764 0.0408 0.688 0.848
## 114 78 1 0.754 0.0414 0.677 0.840
## 115 77 1 0.745 0.0420 0.667 0.832
## 134 74 1 0.735 0.0426 0.656 0.823
## 164 72 1 0.724 0.0432 0.644 0.814
## 169 71 1 0.714 0.0438 0.633 0.805
## 194 70 1 0.704 0.0444 0.622 0.796
##
## tx=ABXITS
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 123 1 0.992 0.0081 0.976 1.000
## 3 120 2 0.975 0.0141 0.948 1.000
## 6 116 1 0.967 0.0163 0.936 0.999
## 8 114 1 0.958 0.0182 0.923 0.995
## 28 108 1 0.950 0.0201 0.911 0.990
## 32 107 1 0.941 0.0218 0.899 0.984
## 41 104 1 0.932 0.0234 0.887 0.979
## 66 102 1 0.923 0.0248 0.875 0.973
## 69 101 1 0.913 0.0262 0.863 0.966
## 71 100 1 0.904 0.0275 0.852 0.960
## 78 99 1 0.895 0.0287 0.841 0.953
## 110 94 1 0.886 0.0299 0.829 0.946
## 118 93 1 0.876 0.0311 0.817 0.939
## 120 92 1 0.867 0.0322 0.806 0.932
## 129 89 1 0.857 0.0333 0.794 0.925
## 139 87 1 0.847 0.0343 0.782 0.917KM2 <- ggsurvplot(KM, data = data,
pval.coord = c(0.03, 0.1),
fun = "cumhaz",xlim=c(0,200),ylim=c(0,0.25),legend.labs = c("ITS","ITS+ABX"),
title = "Cows with at least 1 previous case of clinical mastitis", pval = T, censor=F, conf.int = F,risk.table = F,
surv.plot.height = 5, tables.theme = theme_cleantable(),
ggtheme = theme_bw(base_family = "Times"), xlab="Days in milk",
palette = "Set1",risk.table.fontsize=3,break.time.by = 10)
splots <- list(KM1,KM2)
arrange_ggsurvplots(splots,ncol=2,nrow=1)
data$tx <- fct_relevel(data$tx,"ABXITS","ITS")
CPH <- coxph(Surv(cm2tar, cm2) ~ tx*cm1 + factor(paritycat) + dimdo + drylengthcat + cluster(drydate), data=data)
car::Anova(CPH)cm.cm1.p <- car::Anova(CPH) %>% tail(1)
out <- emmeans(CPH,pairwise ~ tx | cm1,adjust = "none",type = "response") %>% confint
out## $emmeans
## cm1 = 0:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 0.0865 0.0555 Inf 0.0246 0.304
## ITS 0.1497 0.1075 Inf 0.0366 0.612
##
## cm1 = 1:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 0.1973 0.1371 Inf 0.0506 0.770
## ITS 0.4088 0.2842 Inf 0.1046 1.597
##
## Results are averaged over the levels of: paritycat, drylengthcat
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## cm1 = 0:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 0.578 0.104 Inf 0.406 0.823
##
## cm1 = 1:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 0.483 0.127 Inf 0.288 0.808
##
## Results are averaged over the levels of: paritycat, drylengthcat
## Confidence level used: 0.95
## Intervals are back-transformed from the log scaletab <- data.frame(
outcome = "Clinical mastitis",
stratatype = "Clinical mastitis history",
strata = out[["contrasts"]][["cm1"]],
contrast = out[["contrasts"]][["contrast"]],
est = out[["contrasts"]][["ratio"]],
lcl = out[["contrasts"]][["asymp.LCL"]],
ucl = out[["contrasts"]][["asymp.UCL"]]
)
tab <- tab %>% mutate(
estimate = paste0(round(est,2),
" (",
round(lcl,2),
" to ",
round(ucl,2),
")"),
est = NULL,
lcl = NULL,
ucl = NULL
)
tab.cm.cm1 <- tab
tab.cm.cm1data$tx <- fct_relevel(data$tx,"ITS","ABXITS")Dry period length
data %>% filter(!is.na(cm2)) %>%
tabyl(tx,cm2,drylengthcat) %>%
adorn_totals("row") %>%
adorn_percentages("row") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns %>%
adorn_title## $`(0,30]`
## cm2
## tx 0 1
## ITS 91.3% (21) 8.7% (2)
## ABXITS 86.7% (26) 13.3% (4)
## Total 88.7% (47) 11.3% (6)
##
## $`(30,60]`
## cm2
## tx 0 1
## ITS 88.3% (910) 11.7% (121)
## ABXITS 93.6% (1184) 6.4% (81)
## Total 91.2% (2094) 8.8% (202)
##
## $`(60,90]`
## cm2
## tx 0 1
## ITS 73.7% (14) 26.3% (5)
## ABXITS 94.4% (17) 5.6% (1)
## Total 83.8% (31) 16.2% (6)
##
## $`(90,150]`
## cm2
## tx 0 1
## ITS 100.0% (1) 0.0% (0)
## ABXITS - (0) - (0)
## Total 100.0% (1) 0.0% (0)data %>% filter(!is.na(cm2)) %>%
tabyl(tx,cm2,dplength3) %>%
adorn_totals("row") %>%
adorn_percentages("row") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns %>%
adorn_title## $`Mid dplength`
## cm2
## tx 0 1
## ITS 89.1% (319) 10.9% (39)
## ABXITS 96.1% (421) 3.9% (17)
## Total 93.0% (740) 7.0% (56)
##
## $`High dplength`
## cm2
## tx 0 1
## ITS 87.5% (316) 12.5% (45)
## ABXITS 91.2% (396) 8.8% (38)
## Total 89.6% (712) 10.4% (83)
##
## $`Low dplength`
## cm2
## tx 0 1
## ITS 87.6% (311) 12.4% (44)
## ABXITS 93.0% (410) 7.0% (31)
## Total 90.6% (721) 9.4% (75)KM <- survfit(Surv(cm2tar, cm2) ~ tx, data = data %>% filter(dplength3=="Low dplength"))
summary(KM)## Call: survfit(formula = Surv(cm2tar, cm2) ~ tx, data = data %>% filter(dplength3 ==
## "Low dplength"))
##
## tx=ITS
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 2 355 3 0.992 0.00486 0.982 1.000
## 3 351 3 0.983 0.00685 0.970 0.997
## 7 342 1 0.980 0.00741 0.966 0.995
## 8 340 1 0.977 0.00793 0.962 0.993
## 14 334 1 0.974 0.00843 0.958 0.991
## 15 332 1 0.971 0.00890 0.954 0.989
## 17 329 1 0.969 0.00935 0.950 0.987
## 18 328 1 0.966 0.00978 0.947 0.985
## 19 327 1 0.963 0.01018 0.943 0.983
## 21 324 1 0.960 0.01057 0.939 0.981
## 26 321 1 0.957 0.01096 0.935 0.978
## 31 319 1 0.954 0.01133 0.932 0.976
## 32 318 2 0.948 0.01202 0.924 0.971
## 33 316 1 0.945 0.01235 0.921 0.969
## 35 314 2 0.939 0.01299 0.914 0.964
## 37 311 1 0.936 0.01329 0.910 0.962
## 38 310 1 0.933 0.01359 0.906 0.960
## 39 309 1 0.930 0.01387 0.903 0.957
## 40 308 1 0.927 0.01415 0.899 0.955
## 41 307 2 0.921 0.01469 0.892 0.950
## 47 303 1 0.917 0.01495 0.889 0.947
## 52 301 1 0.914 0.01521 0.885 0.945
## 62 300 1 0.911 0.01546 0.882 0.942
## 66 298 1 0.908 0.01571 0.878 0.940
## 73 296 1 0.905 0.01595 0.875 0.937
## 77 294 1 0.902 0.01619 0.871 0.934
## 79 293 2 0.896 0.01666 0.864 0.929
## 81 290 1 0.893 0.01688 0.860 0.927
## 83 288 1 0.890 0.01711 0.857 0.924
## 91 286 1 0.887 0.01733 0.853 0.921
## 92 285 1 0.884 0.01755 0.850 0.919
## 102 282 1 0.880 0.01776 0.846 0.916
## 142 278 1 0.877 0.01798 0.843 0.913
## 169 275 1 0.874 0.01819 0.839 0.911
## 175 273 1 0.871 0.01841 0.836 0.908
## 182 272 1 0.868 0.01861 0.832 0.905
##
## tx=ABXITS
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 441 1 0.998 0.00227 0.993 1.000
## 4 431 1 0.995 0.00323 0.989 1.000
## 6 428 1 0.993 0.00397 0.985 1.000
## 10 416 1 0.991 0.00463 0.982 1.000
## 21 409 1 0.988 0.00521 0.978 0.999
## 22 408 1 0.986 0.00573 0.975 0.997
## 26 406 1 0.983 0.00621 0.971 0.996
## 28 405 1 0.981 0.00666 0.968 0.994
## 40 402 3 0.974 0.00783 0.958 0.989
## 41 399 1 0.971 0.00819 0.955 0.987
## 47 397 1 0.969 0.00852 0.952 0.986
## 48 394 1 0.966 0.00885 0.949 0.984
## 49 393 1 0.964 0.00916 0.946 0.982
## 51 392 1 0.961 0.00946 0.943 0.980
## 52 389 1 0.959 0.00976 0.940 0.978
## 63 386 1 0.956 0.01004 0.937 0.976
## 65 384 1 0.954 0.01032 0.934 0.974
## 71 383 1 0.951 0.01059 0.931 0.972
## 74 382 1 0.949 0.01085 0.928 0.970
## 87 378 1 0.946 0.01111 0.925 0.969
## 120 367 1 0.944 0.01137 0.922 0.966
## 129 363 1 0.941 0.01164 0.919 0.964
## 146 360 1 0.939 0.01189 0.916 0.962
## 154 358 1 0.936 0.01215 0.913 0.960
## 155 357 1 0.933 0.01239 0.909 0.958
## 167 355 1 0.931 0.01263 0.906 0.956
## 187 351 1 0.928 0.01287 0.903 0.954
## 191 350 1 0.926 0.01311 0.900 0.952
## 192 349 1 0.923 0.01333 0.897 0.949KM1 <- ggsurvplot(KM, data = data,
pval.coord = c(0.03, 0.1),
fun = "cumhaz",xlim=c(0,200),ylim=c(0,0.25),legend.labs = c("ITS","ITS+ABX"),
title = "Cows with shorter dry period length", pval = T, censor=F, conf.int = F,risk.table = F,
surv.plot.height = 5, tables.theme = theme_cleantable(),
ggtheme = theme_bw(base_family = "Times"), xlab="Days in milk",
palette = "Set1",risk.table.fontsize=3,break.time.by = 10)
#KM1
KM <- survfit(Surv(cm2tar, cm2) ~ tx, data = data %>% filter(dplength3=="Mid dplength"))
summary(KM)## Call: survfit(formula = Surv(cm2tar, cm2) ~ tx, data = data %>% filter(dplength3 ==
## "Mid dplength"))
##
## tx=ITS
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 358 1 0.997 0.00279 0.992 1.000
## 2 356 2 0.992 0.00483 0.982 1.000
## 3 351 1 0.989 0.00558 0.978 1.000
## 4 349 1 0.986 0.00624 0.974 0.998
## 10 344 2 0.980 0.00740 0.966 0.995
## 11 341 1 0.977 0.00792 0.962 0.993
## 22 338 1 0.974 0.00841 0.958 0.991
## 23 337 1 0.972 0.00887 0.954 0.989
## 26 336 1 0.969 0.00930 0.951 0.987
## 27 335 1 0.966 0.00971 0.947 0.985
## 30 334 1 0.963 0.01010 0.943 0.983
## 34 332 1 0.960 0.01048 0.940 0.981
## 38 330 1 0.957 0.01085 0.936 0.979
## 44 328 1 0.954 0.01120 0.932 0.976
## 47 327 1 0.951 0.01154 0.929 0.974
## 52 324 1 0.948 0.01187 0.925 0.972
## 58 319 1 0.945 0.01220 0.922 0.970
## 66 318 1 0.942 0.01252 0.918 0.967
## 68 317 1 0.939 0.01283 0.915 0.965
## 80 315 1 0.936 0.01313 0.911 0.962
## 81 314 1 0.933 0.01342 0.907 0.960
## 88 313 1 0.930 0.01371 0.904 0.958
## 89 312 1 0.927 0.01398 0.900 0.955
## 101 308 1 0.924 0.01426 0.897 0.953
## 105 307 1 0.921 0.01453 0.893 0.950
## 106 306 1 0.918 0.01479 0.890 0.948
## 113 302 1 0.915 0.01505 0.886 0.945
## 115 300 1 0.912 0.01530 0.883 0.943
## 116 299 1 0.909 0.01555 0.879 0.940
## 134 294 1 0.906 0.01580 0.876 0.938
## 141 291 1 0.903 0.01605 0.872 0.935
## 149 289 1 0.900 0.01630 0.869 0.932
## 170 286 1 0.897 0.01654 0.865 0.930
## 180 283 1 0.894 0.01679 0.861 0.927
## 191 279 1 0.890 0.01703 0.858 0.924
## 194 278 1 0.887 0.01727 0.854 0.922
## 199 277 1 0.884 0.01750 0.850 0.919
##
## tx=ABXITS
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 2 438 1 0.998 0.00228 0.993 1.000
## 4 434 1 0.995 0.00323 0.989 1.000
## 6 430 2 0.991 0.00458 0.982 1.000
## 23 415 1 0.988 0.00516 0.978 0.999
## 27 414 2 0.984 0.00614 0.972 0.996
## 32 411 1 0.981 0.00657 0.968 0.994
## 37 408 1 0.979 0.00698 0.965 0.993
## 41 405 1 0.976 0.00737 0.962 0.991
## 61 401 1 0.974 0.00775 0.959 0.989
## 62 400 1 0.972 0.00810 0.956 0.988
## 76 397 1 0.969 0.00844 0.953 0.986
## 94 393 1 0.967 0.00877 0.950 0.984
## 101 390 1 0.964 0.00909 0.946 0.982
## 102 389 1 0.962 0.00940 0.943 0.980
## 139 383 1 0.959 0.00971 0.940 0.978KM2 <- ggsurvplot(KM, data = data,
pval.coord = c(0.03, 0.1),
fun = "cumhaz",xlim=c(0,200),ylim=c(0,0.25),legend.labs = c("ITS","ITS+ABX"),
title = "Cows with middle dry period length", pval = T, censor=F, conf.int = F,risk.table = F,
surv.plot.height = 5, tables.theme = theme_cleantable(),
ggtheme = theme_bw(base_family = "Times"), xlab="Days in milk",
palette = "Set1",risk.table.fontsize=3,break.time.by = 10)
#KM2
KM <- survfit(Surv(cm2tar, cm2) ~ tx, data = data %>% filter(dplength3=="High dplength"))
summary(KM)## Call: survfit(formula = Surv(cm2tar, cm2) ~ tx, data = data %>% filter(dplength3 ==
## "High dplength"))
##
## tx=ITS
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 361 1 0.997 0.00277 0.992 1.000
## 2 360 3 0.989 0.00551 0.978 1.000
## 3 357 1 0.986 0.00615 0.974 0.998
## 4 356 2 0.981 0.00726 0.966 0.995
## 5 353 1 0.978 0.00775 0.963 0.993
## 6 351 2 0.972 0.00865 0.955 0.989
## 7 349 3 0.964 0.00983 0.945 0.983
## 8 345 2 0.958 0.01054 0.938 0.979
## 9 341 1 0.956 0.01087 0.934 0.977
## 13 337 1 0.953 0.01121 0.931 0.975
## 15 335 1 0.950 0.01153 0.927 0.973
## 17 334 1 0.947 0.01184 0.924 0.970
## 23 329 1 0.944 0.01215 0.921 0.968
## 25 326 1 0.941 0.01245 0.917 0.966
## 28 324 1 0.938 0.01275 0.914 0.964
## 30 321 1 0.935 0.01304 0.910 0.961
## 31 320 2 0.930 0.01360 0.903 0.957
## 35 315 1 0.927 0.01387 0.900 0.954
## 47 313 1 0.924 0.01414 0.896 0.952
## 48 312 1 0.921 0.01440 0.893 0.949
## 50 311 2 0.915 0.01490 0.886 0.944
## 56 309 1 0.912 0.01515 0.883 0.942
## 64 307 1 0.909 0.01538 0.879 0.939
## 66 306 1 0.906 0.01562 0.876 0.937
## 68 305 1 0.903 0.01585 0.872 0.934
## 74 304 1 0.900 0.01607 0.869 0.932
## 77 303 1 0.897 0.01629 0.866 0.929
## 90 300 1 0.894 0.01651 0.862 0.927
## 95 298 1 0.891 0.01672 0.859 0.924
## 96 297 1 0.888 0.01693 0.855 0.922
## 114 293 1 0.885 0.01714 0.852 0.919
## 120 291 1 0.882 0.01735 0.849 0.917
## 121 290 1 0.879 0.01756 0.845 0.914
## 161 285 1 0.876 0.01776 0.842 0.911
## 164 283 1 0.873 0.01797 0.838 0.909
## 172 281 1 0.870 0.01817 0.835 0.906
##
## tx=ABXITS
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 434 2 0.995 0.00325 0.989 1.000
## 2 432 2 0.991 0.00459 0.982 1.000
## 3 429 3 0.984 0.00605 0.972 0.996
## 4 426 1 0.982 0.00646 0.969 0.994
## 5 423 2 0.977 0.00722 0.963 0.991
## 8 417 1 0.975 0.00757 0.960 0.990
## 9 414 1 0.972 0.00791 0.957 0.988
## 13 409 1 0.970 0.00824 0.954 0.986
## 14 408 1 0.967 0.00856 0.951 0.984
## 24 404 2 0.963 0.00916 0.945 0.981
## 29 399 1 0.960 0.00945 0.942 0.979
## 31 398 1 0.958 0.00973 0.939 0.977
## 41 394 1 0.955 0.01000 0.936 0.975
## 46 392 1 0.953 0.01027 0.933 0.973
## 48 391 1 0.951 0.01053 0.930 0.971
## 63 390 1 0.948 0.01078 0.927 0.969
## 66 388 2 0.943 0.01127 0.921 0.966
## 69 386 1 0.941 0.01150 0.918 0.964
## 70 385 1 0.938 0.01173 0.916 0.962
## 78 384 1 0.936 0.01195 0.913 0.960
## 86 382 1 0.933 0.01216 0.910 0.958
## 107 379 1 0.931 0.01238 0.907 0.956
## 110 377 1 0.928 0.01259 0.904 0.954
## 111 376 1 0.926 0.01280 0.901 0.951
## 118 375 2 0.921 0.01320 0.896 0.947
## 124 371 1 0.919 0.01339 0.893 0.945
## 128 370 1 0.916 0.01358 0.890 0.943
## 150 364 1 0.914 0.01378 0.887 0.941
## 193 355 1 0.911 0.01398 0.884 0.939
## 197 354 1 0.908 0.01417 0.881 0.937KM3 <- ggsurvplot(KM, data = data,
pval.coord = c(0.03, 0.1),
fun = "cumhaz",xlim=c(0,200),ylim=c(0,0.25),legend.labs = c("ITS","ITS+ABX"),
title = "Cows with longer dry period length", pval = T, censor=F, conf.int = F,risk.table = F,
surv.plot.height = 5, tables.theme = theme_cleantable(),
ggtheme = theme_bw(base_family = "Times"), xlab="Days in milk",
palette = "Set1",risk.table.fontsize=3,break.time.by = 10)
KM3
splots <- list(KM1,KM2,KM3)
arrange_ggsurvplots(splots,ncol=2,nrow=2)
data$tx <- fct_relevel(data$tx,"ABXITS","ITS")
CPH <- coxph(Surv(cm2tar, cm2) ~ tx*dplength3 + cm1 + factor(paritycat) + dimdo + cluster(drydate), data=data)
car::Anova(CPH)cm.dp.p <- car::Anova(CPH) %>% tail(1)
out <- emmeans(CPH,pairwise ~ tx | dplength3,adjust = "none",type = "response") %>% confint
out## $emmeans
## dplength3 = Low dplength:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 3.39 1.658 Inf 1.299 8.84
## ITS 5.70 3.243 Inf 1.868 17.39
##
## dplength3 = Mid dplength:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 1.67 0.846 Inf 0.619 4.51
## ITS 4.68 2.721 Inf 1.498 14.63
##
## dplength3 = High dplength:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 4.00 2.129 Inf 1.407 11.35
## ITS 5.75 3.313 Inf 1.859 17.78
##
## Results are averaged over the levels of: cm1, paritycat
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## dplength3 = Low dplength:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 0.595 0.1711 Inf 0.338 1.045
##
## dplength3 = Mid dplength:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 0.357 0.0876 Inf 0.221 0.577
##
## dplength3 = High dplength:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 0.695 0.1316 Inf 0.480 1.007
##
## Results are averaged over the levels of: cm1, paritycat
## Confidence level used: 0.95
## Intervals are back-transformed from the log scaletab <- data.frame(
outcome = "Clinical mastitis",
stratatype = "Dry period length",
strata = out[["contrasts"]][["dplength3"]],
contrast = out[["contrasts"]][["contrast"]],
est = out[["contrasts"]][["ratio"]],
lcl = out[["contrasts"]][["asymp.LCL"]],
ucl = out[["contrasts"]][["asymp.UCL"]]
)
tab <- tab %>% mutate(
estimate = paste0(round(est,2),
" (",
round(lcl,2),
" to ",
round(ucl,2),
")"),
est = NULL,
lcl = NULL,
ucl = NULL
)
tab.cm.dp <- tab
tab.cm.dpdata$tx <- fct_relevel(data$tx,"ITS","ABXITS")cm.hr <- rbind(tab.cm,
tab.cm.cm1,
tab.cm.dp,
tab.cm.parity)
cm.p <- rbind(cm.cm1.p,
cm.parity.p,
cm.dp.p)
cm.p$outcome <- "Clinical mastitis"Culling during the dry period
Raw numbers
data %>%
tabyl(tx,remdp) %>%
adorn_totals("row") %>%
adorn_percentages("row") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns %>%
adorn_titleSurvival analysis
#Will limit dplength to 70 as most cows calved within 70d
#data$remdp[data$remdptar > 60] <- 0
#data$remdptar[data$remdptar > 60] <- 60KM <- survfit(Surv(remdptar, remdp) ~ tx, data = data)
summary(KM)## Call: survfit(formula = Surv(remdptar, remdp) ~ tx, data = data)
##
## tx=ITS
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 2 1103 1 0.999 0.000906 0.9973 1.000
## 3 1102 1 0.998 0.001281 0.9957 1.000
## 4 1101 1 0.997 0.001568 0.9942 1.000
## 5 1099 1 0.996 0.001810 0.9928 1.000
## 6 1098 2 0.995 0.002216 0.9902 0.999
## 8 1096 1 0.994 0.002392 0.9890 0.998
## 9 1094 1 0.993 0.002557 0.9877 0.998
## 11 1093 1 0.992 0.002711 0.9865 0.997
## 12 1092 3 0.989 0.003127 0.9830 0.995
## 15 1084 1 0.988 0.003254 0.9818 0.995
## 16 1083 1 0.987 0.003377 0.9807 0.994
## 19 1079 1 0.986 0.003495 0.9795 0.993
## 29 1069 3 0.984 0.003834 0.9761 0.991
## 36 1038 1 0.983 0.003945 0.9750 0.990
## 38 1014 1 0.982 0.004059 0.9738 0.990
## 41 936 1 0.981 0.004188 0.9725 0.989
## 49 475 1 0.979 0.004660 0.9695 0.988
## 60 35 1 0.951 0.027926 0.8974 1.000
## 66 12 1 0.871 0.080049 0.7278 1.000
## 69 9 1 0.775 0.115740 0.5779 1.000
## 70 8 1 0.678 0.135863 0.4575 1.000
## 75 6 1 0.565 0.153139 0.3320 0.961
## 76 5 1 0.452 0.158798 0.2269 0.900
## 77 4 1 0.339 0.154124 0.1390 0.826
## 100 2 1 0.169 0.142454 0.0326 0.880
##
## tx=ABXITS
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 4 1327 2 0.998 0.00106 0.996 1.000
## 9 1321 1 0.998 0.00131 0.995 1.000
## 10 1320 1 0.997 0.00151 0.994 1.000
## 12 1318 1 0.996 0.00169 0.993 1.000
## 13 1317 1 0.995 0.00185 0.992 0.999
## 19 1310 1 0.995 0.00199 0.991 0.999
## 24 1303 1 0.994 0.00213 0.990 0.998
## 35 1266 1 0.993 0.00227 0.989 0.998
## 44 967 1 0.992 0.00249 0.987 0.997
## 49 545 1 0.990 0.00308 0.984 0.996
## 60 32 1 0.959 0.03061 0.901 1.000
## 70 7 1 0.822 0.12957 0.604 1.000
## 71 5 1 0.658 0.17995 0.385 1.000ggsurvplot(KM, data = data,ylim = c(0,0.05),xlim = c(0,55),
pval.coord = c(0.03, 0.02),
fun = "cumhaz",legend.labs = c("ITS","ITS+ABX"),
title = "", pval = T, censor=F, conf.int = F,risk.table = T,
surv.plot.height = 5, tables.theme = theme_cleantable(),
ggtheme = theme_bw(base_family = "Times"), xlab="Days post dry-off",
palette = "Set1",risk.table.fontsize=3,break.time.by = 10)
Cox models
library(broom)
SR <- coxph(Surv(remdptar, remdp) ~ tx + dimdo + factor(paritycat) + cm1 + cluster(drydate), data=data)
car::Anova(SR)output <- SR %>% tidy(exp=TRUE,conf.int=T) %>% select(term,estimate,conf.low,conf.high,p.value)
output$estimate <- output$estimate %>% round(2)
output$conf.low <- output$conf.low %>% round(2)
output$conf.high <- output$conf.high %>% round(2)
output$p.value <- output$p.value %>% round(4)
outputculldp.main <- output %>% filter(term=="txABXITS") %>% mutate(
estimate = paste0(round(estimate,2),
" (",
round(conf.low,2),
" to ",
round(conf.high,2),
")"),
conf.low = NULL,
conf.high = NULL,
p.value = NULL)
culldp.main$contrast <- "ABXITS / ITS"
culldp.main$strata <- NA
culldp.main$outcome <- "Removal during the dry period"
culldp.main$stratatype <- "All cows"
culldp.main$term <- NULL
tab.cull.dp <- culldp.mainSR <- cox.zph(SR)
SR## chisq df p
## tx 0.079916 1 0.78
## dimdo 0.279254 1 0.60
## factor(paritycat) 0.868549 3 0.83
## cm1 0.000838 1 0.98
## GLOBAL 1.667169 6 0.95ggcoxzph(SR)
Stratified analysis
Parity
data %>%
tabyl(tx,remdp,paritycat) %>%
adorn_totals("row") %>%
adorn_percentages("row") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns %>%
adorn_title## $`2`
## remdp
## tx 0 1
## ITS 97.7% (424) 2.3% (10)
## ABXITS 99.2% (610) 0.8% (5)
## Total 98.6% (1034) 1.4% (15)
##
## $`3`
## remdp
## tx 0 1
## ITS 97.5% (234) 2.5% (6)
## ABXITS 99.2% (245) 0.8% (2)
## Total 98.4% (479) 1.6% (8)
##
## $`4`
## remdp
## tx 0 1
## ITS 98.1% (259) 1.9% (5)
## ABXITS 98.6% (279) 1.4% (4)
## Total 98.4% (538) 1.6% (9)
##
## $`5+`
## remdp
## tx 0 1
## ITS 94.6% (157) 5.4% (9)
## ABXITS 98.4% (179) 1.6% (3)
## Total 96.6% (336) 3.4% (12)data$tx <- fct_relevel(data$tx,"ABXITS","ITS")
CPH <- coxph(Surv(remdptar, remdp) ~ tx*factor(paritycat) + dimdo + cm1 + cluster(drydate), data=data)
car::Anova(CPH)cull.dp.parity.p <- car::Anova(CPH) %>% tail(1)
out <- emmeans(CPH,pairwise ~ tx | factor(paritycat),adjust = "none",type = "response") %>% confint
out## $emmeans
## paritycat = 2:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 0.603 0.508 Inf 0.1159 3.14
## ITS 1.612 1.758 Inf 0.1899 13.68
##
## paritycat = 3:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 0.603 0.700 Inf 0.0621 5.86
## ITS 1.291 1.403 Inf 0.1534 10.87
##
## paritycat = 4:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 0.885 1.048 Inf 0.0869 9.01
## ITS 1.478 1.310 Inf 0.2603 8.40
##
## paritycat = 5+:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 1.154 1.377 Inf 0.1111 11.98
## ITS 2.569 2.645 Inf 0.3414 19.33
##
## Results are averaged over the levels of: cm1
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## paritycat = 2:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 0.374 0.208 Inf 0.126 1.11
##
## paritycat = 3:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 0.467 0.343 Inf 0.111 1.97
##
## paritycat = 4:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 0.598 0.415 Inf 0.154 2.33
##
## paritycat = 5+:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 0.449 0.249 Inf 0.152 1.33
##
## Results are averaged over the levels of: cm1
## Confidence level used: 0.95
## Intervals are back-transformed from the log scaletab <- data.frame(
outcome = "Removal during the dry period",
stratatype = "Parity",
strata = out[["contrasts"]][["paritycat"]],
contrast = out[["contrasts"]][["contrast"]],
est = out[["contrasts"]][["ratio"]],
lcl = out[["contrasts"]][["asymp.LCL"]],
ucl = out[["contrasts"]][["asymp.UCL"]]
)
tab <- tab %>% mutate(
estimate = paste0(round(est,2),
" (",
round(lcl,2),
" to ",
round(ucl,2),
")"),
est = NULL,
lcl = NULL,
ucl = NULL
)
tab.cull.dp.parity <- tab
tab.cull.dp.paritydata$tx <- fct_relevel(data$tx,"ITS","ABXITS")Clinical mastitis history
data %>%
tabyl(tx,remdp,cm1) %>%
adorn_totals("row") %>%
adorn_percentages("row") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns %>%
adorn_title## $`0`
## remdp
## tx 0 1
## ITS 97.5% (951) 2.5% (24)
## ABXITS 98.8% (1190) 1.2% (14)
## Total 98.3% (2141) 1.7% (38)
##
## $`1`
## remdp
## tx 0 1
## ITS 95.3% (123) 4.7% (6)
## ABXITS 100.0% (123) 0.0% (0)
## Total 97.6% (246) 2.4% (6)data$tx <- fct_relevel(data$tx,"ABXITS","ITS")
CPH <- coxph(Surv(remdptar, remdp) ~ tx*cm1 + factor(paritycat) + dimdo + cluster(drydate), data=data)
car::Anova(CPH)cull.dp.cm1.p <- car::Anova(CPH) %>% tail(1)
out <- emmeans(CPH,pairwise ~ tx | cm1,adjust = "none",type = "response") %>% confint
out## $emmeans
## cm1 = 0:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 0.738 0.708 Inf 1.13e-01 4.83
## ITS 1.426 1.428 Inf 2.00e-01 10.15
##
## cm1 = 1:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 0.000 0.000 Inf 1.00e-08 0.00
## ITS 1.328 1.209 Inf 2.23e-01 7.91
##
## Results are averaged over the levels of: paritycat
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## cm1 = 0:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 0.518 1.58e-01 Inf 2.85e-01 9.4e-01
##
## cm1 = 1:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 0.000 2.00e-08 Inf 1.00e-08 1.0e-07
##
## Results are averaged over the levels of: paritycat
## Confidence level used: 0.95
## Intervals are back-transformed from the log scaletab <- data.frame(
outcome = "Removal during the dry period",
stratatype = "Clinical mastitis history",
strata = out[["contrasts"]][["cm1"]],
contrast = out[["contrasts"]][["contrast"]],
est = out[["contrasts"]][["ratio"]],
lcl = out[["contrasts"]][["asymp.LCL"]],
ucl = out[["contrasts"]][["asymp.UCL"]]
)
tab <- tab %>% mutate(
estimate = paste0(round(est,2),
" (",
round(lcl,2),
" to ",
round(ucl,2),
")"),
est = NULL,
lcl = NULL,
ucl = NULL
)
tab.cull.dp.cm1 <- tab
tab.cull.dp.cm1data$tx <- fct_relevel(data$tx,"ITS","ABXITS")cull.dp.hr <- rbind(tab.cull.dp,
tab.cull.dp.cm1,
tab.cull.dp.parity)
cull.dp.p <- rbind(cull.dp.cm1.p,cull.dp.parity.p)
cull.dp.p$outcome <- "Removal during the dry period"Culling during the subsequent lactation
Raw numbers
data %>% filter(!is.na(rem2)) %>%
tabyl(tx,rem2) %>%
adorn_totals("row") %>%
adorn_percentages("row") %>%
adorn_pct_formatting(digits = 1) %>%
adorn_ns %>%
adorn_titleSurvival analysis
#Will limit to 200 DIM
data$rem2[data$rem2tar >= 200] <- 0
data$rem2tar[data$rem2tar > 200] <- 200KM <- survfit(Surv(rem2tar, rem2) ~ tx, data = data)
summary(KM)## Call: survfit(formula = Surv(rem2tar, rem2) ~ tx, data = data)
##
## 44 observations deleted due to missingness
## tx=ITS
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 1074 1 0.999 0.000931 0.997 1.000
## 2 1073 4 0.995 0.002077 0.991 0.999
## 3 1069 2 0.993 0.002455 0.989 0.998
## 4 1067 4 0.990 0.003072 0.984 0.996
## 5 1063 5 0.985 0.003697 0.978 0.992
## 6 1058 3 0.982 0.004023 0.974 0.990
## 7 1055 2 0.980 0.004225 0.972 0.989
## 8 1053 3 0.978 0.004510 0.969 0.987
## 10 1050 5 0.973 0.004946 0.963 0.983
## 11 1045 4 0.969 0.005266 0.959 0.980
## 12 1041 3 0.966 0.005492 0.956 0.977
## 13 1038 2 0.965 0.005637 0.954 0.976
## 14 1036 2 0.963 0.005778 0.951 0.974
## 15 1034 2 0.961 0.005915 0.949 0.973
## 16 1032 2 0.959 0.006048 0.947 0.971
## 18 1030 2 0.957 0.006178 0.945 0.969
## 20 1028 3 0.954 0.006367 0.942 0.967
## 21 1025 1 0.953 0.006429 0.941 0.966
## 22 1024 1 0.953 0.006490 0.940 0.965
## 23 1023 2 0.951 0.006609 0.938 0.964
## 24 1021 2 0.949 0.006726 0.936 0.962
## 26 1019 1 0.948 0.006784 0.935 0.961
## 27 1018 1 0.947 0.006841 0.934 0.960
## 29 1017 2 0.945 0.006953 0.932 0.959
## 30 1015 1 0.944 0.007008 0.930 0.958
## 31 1014 1 0.943 0.007063 0.929 0.957
## 32 1013 2 0.941 0.007170 0.927 0.955
## 33 1011 2 0.939 0.007276 0.925 0.954
## 36 1009 2 0.938 0.007380 0.923 0.952
## 41 1007 1 0.937 0.007431 0.922 0.951
## 43 1006 1 0.936 0.007482 0.921 0.951
## 44 1005 2 0.934 0.007582 0.919 0.949
## 47 1003 2 0.932 0.007680 0.917 0.947
## 50 1001 1 0.931 0.007729 0.916 0.946
## 52 1000 2 0.929 0.007825 0.914 0.945
## 53 998 1 0.928 0.007872 0.913 0.944
## 54 997 2 0.926 0.007966 0.911 0.942
## 56 995 2 0.925 0.008058 0.909 0.941
## 60 993 1 0.924 0.008103 0.908 0.940
## 63 992 1 0.923 0.008148 0.907 0.939
## 66 991 1 0.922 0.008193 0.906 0.938
## 69 990 1 0.921 0.008238 0.905 0.937
## 72 989 1 0.920 0.008282 0.904 0.936
## 75 988 1 0.919 0.008326 0.903 0.935
## 76 987 1 0.918 0.008369 0.902 0.935
## 79 986 1 0.917 0.008412 0.901 0.934
## 80 985 1 0.916 0.008455 0.900 0.933
## 81 984 1 0.915 0.008497 0.899 0.932
## 82 983 1 0.914 0.008540 0.898 0.931
## 87 982 1 0.913 0.008582 0.897 0.930
## 88 981 1 0.912 0.008623 0.896 0.930
## 91 980 1 0.912 0.008665 0.895 0.929
## 94 979 2 0.910 0.008746 0.893 0.927
## 95 977 1 0.909 0.008787 0.892 0.926
## 96 976 2 0.907 0.008867 0.890 0.924
## 104 974 2 0.905 0.008946 0.888 0.923
## 106 972 1 0.904 0.008985 0.887 0.922
## 107 971 1 0.903 0.009024 0.886 0.921
## 108 970 1 0.902 0.009063 0.885 0.920
## 109 969 1 0.901 0.009101 0.884 0.919
## 111 968 1 0.900 0.009139 0.883 0.918
## 113 967 1 0.899 0.009177 0.882 0.918
## 116 966 1 0.899 0.009214 0.881 0.917
## 117 965 1 0.898 0.009252 0.880 0.916
## 118 964 1 0.897 0.009289 0.879 0.915
## 119 963 1 0.896 0.009326 0.878 0.914
## 122 962 1 0.895 0.009363 0.877 0.913
## 123 961 1 0.894 0.009399 0.876 0.912
## 126 960 1 0.893 0.009435 0.875 0.912
## 127 959 1 0.892 0.009471 0.874 0.911
## 128 958 1 0.891 0.009507 0.873 0.910
## 129 957 1 0.890 0.009543 0.872 0.909
## 130 956 2 0.888 0.009613 0.870 0.907
## 131 954 1 0.887 0.009648 0.869 0.906
## 135 953 1 0.886 0.009683 0.868 0.906
## 136 952 1 0.885 0.009717 0.867 0.905
## 137 951 2 0.884 0.009785 0.865 0.903
## 148 949 1 0.883 0.009819 0.864 0.902
## 149 948 2 0.881 0.009887 0.862 0.900
## 153 946 1 0.880 0.009920 0.861 0.900
## 155 945 1 0.879 0.009953 0.860 0.899
## 157 944 1 0.878 0.009986 0.859 0.898
## 161 943 1 0.877 0.010019 0.858 0.897
## 163 942 1 0.876 0.010051 0.857 0.896
## 167 941 1 0.875 0.010083 0.856 0.895
## 169 940 1 0.874 0.010116 0.855 0.894
## 171 939 1 0.873 0.010148 0.854 0.893
## 172 938 1 0.872 0.010179 0.853 0.893
## 175 937 1 0.872 0.010211 0.852 0.892
## 178 936 1 0.871 0.010243 0.851 0.891
## 179 935 1 0.870 0.010274 0.850 0.890
## 180 934 1 0.869 0.010305 0.849 0.889
## 181 933 1 0.868 0.010336 0.848 0.888
## 184 932 3 0.865 0.010428 0.845 0.886
## 189 929 1 0.864 0.010458 0.844 0.885
## 192 928 1 0.863 0.010488 0.843 0.884
## 193 927 1 0.862 0.010518 0.842 0.883
## 196 926 1 0.861 0.010548 0.841 0.882
## 197 925 2 0.859 0.010607 0.839 0.880
##
## tx=ABXITS
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 2 1313 6 0.995 0.00186 0.992 0.999
## 3 1307 7 0.990 0.00273 0.985 0.995
## 4 1300 5 0.986 0.00321 0.980 0.993
## 5 1295 3 0.984 0.00346 0.977 0.991
## 6 1292 6 0.979 0.00392 0.972 0.987
## 7 1286 3 0.977 0.00412 0.969 0.985
## 8 1283 7 0.972 0.00457 0.963 0.981
## 9 1276 6 0.967 0.00491 0.958 0.977
## 10 1270 4 0.964 0.00513 0.954 0.974
## 11 1266 2 0.963 0.00523 0.952 0.973
## 12 1264 3 0.960 0.00538 0.950 0.971
## 13 1261 1 0.960 0.00543 0.949 0.970
## 14 1260 2 0.958 0.00553 0.947 0.969
## 15 1258 3 0.956 0.00567 0.945 0.967
## 17 1255 1 0.955 0.00572 0.944 0.966
## 18 1254 2 0.954 0.00581 0.942 0.965
## 19 1252 3 0.951 0.00594 0.940 0.963
## 21 1249 1 0.950 0.00599 0.939 0.962
## 23 1248 1 0.950 0.00603 0.938 0.962
## 25 1247 2 0.948 0.00612 0.936 0.960
## 26 1245 1 0.947 0.00616 0.935 0.960
## 27 1244 2 0.946 0.00624 0.934 0.958
## 32 1242 1 0.945 0.00628 0.933 0.958
## 33 1241 1 0.944 0.00632 0.932 0.957
## 35 1240 3 0.942 0.00644 0.930 0.955
## 36 1237 3 0.940 0.00656 0.927 0.953
## 37 1234 1 0.939 0.00660 0.926 0.952
## 38 1233 2 0.938 0.00668 0.925 0.951
## 39 1231 1 0.937 0.00672 0.924 0.950
## 40 1230 1 0.936 0.00675 0.923 0.949
## 41 1229 3 0.934 0.00686 0.920 0.947
## 46 1226 1 0.933 0.00690 0.920 0.947
## 47 1225 2 0.931 0.00697 0.918 0.945
## 51 1223 2 0.930 0.00704 0.916 0.944
## 53 1221 1 0.929 0.00708 0.915 0.943
## 56 1220 1 0.928 0.00711 0.915 0.942
## 57 1219 1 0.928 0.00715 0.914 0.942
## 59 1218 2 0.926 0.00722 0.912 0.940
## 63 1216 2 0.925 0.00729 0.910 0.939
## 64 1214 1 0.924 0.00732 0.910 0.938
## 66 1213 1 0.923 0.00735 0.909 0.938
## 67 1212 2 0.922 0.00742 0.907 0.936
## 79 1210 1 0.921 0.00745 0.906 0.936
## 80 1209 1 0.920 0.00749 0.905 0.935
## 84 1208 2 0.919 0.00755 0.904 0.933
## 86 1206 2 0.917 0.00761 0.902 0.932
## 88 1204 3 0.915 0.00771 0.900 0.930
## 89 1201 1 0.914 0.00774 0.899 0.929
## 92 1200 1 0.913 0.00777 0.898 0.929
## 93 1199 1 0.912 0.00780 0.897 0.928
## 94 1198 1 0.912 0.00783 0.896 0.927
## 96 1197 1 0.911 0.00786 0.896 0.926
## 97 1196 1 0.910 0.00789 0.895 0.926
## 98 1195 1 0.909 0.00792 0.894 0.925
## 99 1194 2 0.908 0.00798 0.892 0.924
## 102 1192 1 0.907 0.00801 0.892 0.923
## 104 1191 1 0.906 0.00804 0.891 0.922
## 108 1190 1 0.906 0.00807 0.890 0.922
## 112 1189 2 0.904 0.00813 0.888 0.920
## 115 1187 2 0.903 0.00819 0.887 0.919
## 119 1185 1 0.902 0.00821 0.886 0.918
## 120 1184 1 0.901 0.00824 0.885 0.917
## 122 1183 2 0.899 0.00830 0.883 0.916
## 123 1181 2 0.898 0.00835 0.882 0.914
## 124 1179 1 0.897 0.00838 0.881 0.914
## 127 1178 1 0.896 0.00841 0.880 0.913
## 131 1177 2 0.895 0.00846 0.878 0.912
## 134 1175 1 0.894 0.00849 0.878 0.911
## 137 1174 1 0.893 0.00852 0.877 0.910
## 139 1173 1 0.893 0.00854 0.876 0.910
## 141 1172 1 0.892 0.00857 0.875 0.909
## 144 1171 1 0.891 0.00860 0.874 0.908
## 147 1170 1 0.890 0.00862 0.874 0.907
## 148 1169 1 0.890 0.00865 0.873 0.907
## 150 1168 1 0.889 0.00868 0.872 0.906
## 152 1167 1 0.888 0.00870 0.871 0.905
## 153 1166 1 0.887 0.00873 0.870 0.905
## 156 1165 1 0.887 0.00875 0.870 0.904
## 157 1164 1 0.886 0.00878 0.869 0.903
## 158 1163 1 0.885 0.00880 0.868 0.902
## 159 1162 1 0.884 0.00883 0.867 0.902
## 161 1161 2 0.883 0.00888 0.865 0.900
## 167 1159 1 0.882 0.00890 0.865 0.900
## 168 1158 2 0.880 0.00895 0.863 0.898
## 169 1156 1 0.880 0.00898 0.862 0.897
## 170 1155 1 0.879 0.00900 0.861 0.897
## 172 1154 1 0.878 0.00903 0.861 0.896
## 173 1153 1 0.877 0.00905 0.860 0.895
## 175 1152 1 0.877 0.00908 0.859 0.895
## 180 1151 2 0.875 0.00912 0.857 0.893
## 188 1149 1 0.874 0.00915 0.857 0.892
## 189 1148 1 0.874 0.00917 0.856 0.892
## 193 1147 1 0.873 0.00920 0.855 0.891
## 195 1146 1 0.872 0.00922 0.854 0.890
## 196 1145 1 0.871 0.00924 0.853 0.890ggsurvplot(KM, data = data,ylim = c(0,0.2),xlim = c(0,200),
pval.coord = c(0.03, 0.1),
fun = "cumhaz",legend.labs = c("ITS","ITS+ABX"),
title = "", pval = T, censor=F, conf.int = F,risk.table = T,
surv.plot.height = 5, tables.theme = theme_cleantable(),
ggtheme = theme_bw(base_family = "Times"), xlab="Days in milk",
palette = "Set1",risk.table.fontsize=3,break.time.by = 10)
Cox models
library(broom)
SR <- coxph(Surv(rem2tar, rem2) ~ tx + dimdo + factor(paritycat) + cm1 + drylength + cluster(drydate), data=data)
car::Anova(SR)output <- SR %>% tidy(exp=TRUE,conf.int=T) %>% select(term,estimate,conf.low,conf.high,p.value)
output$estimate <- output$estimate %>% round(2)
output$conf.low <- output$conf.low %>% round(2)
output$conf.high <- output$conf.high %>% round(2)
output$p.value <- output$p.value %>% round(4)
outputSR <- cox.zph(SR)
SR## chisq df p
## tx 1.559 1 0.212
## dimdo 0.198 1 0.656
## factor(paritycat) 2.791 3 0.425
## cm1 5.424 1 0.020
## drylength 6.459 1 0.011
## GLOBAL 15.850 7 0.027ggcoxzph(SR)
Redo model with stratified baseline hazards for CM1
library(broom)
SR <- coxph(Surv(rem2tar, rem2) ~ tx + dimdo + factor(paritycat) + drylength + strata(cm1) + cluster(drydate), data=data)
car::Anova(SR)output <- SR %>% tidy(exp=TRUE,conf.int=T) %>% select(term,estimate,conf.low,conf.high,p.value)
output$estimate <- output$estimate %>% round(2)
output$conf.low <- output$conf.low %>% round(2)
output$conf.high <- output$conf.high %>% round(2)
output$p.value <- output$p.value %>% round(4)
outputcull.lact.main <- output %>% filter(term=="txABXITS") %>% mutate(
estimate = paste0(round(estimate,2),
" (",
round(conf.low,2),
" to ",
round(conf.high,2),
")"),
conf.low = NULL,
conf.high = NULL,
p.value = NULL)
cull.lact.main$contrast <- "ABXITS / ITS"
cull.lact.main$strata <- NA
cull.lact.main$outcome <- "Removal during 1-200 DIM"
cull.lact.main$stratatype <- "All cows"
cull.lact.main$term <- NULL
tab.cull.lact <- cull.lact.mainSR <- cox.zph(SR)
SR## chisq df p
## tx 1.7884 1 0.1811
## dimdo 0.0174 1 0.8950
## factor(paritycat) 1.7177 3 0.6330
## drylength 6.7163 1 0.0096
## GLOBAL 10.4524 6 0.1069ggcoxzph(SR)
Stratified models
Parity
data$tx <- fct_relevel(data$tx,"ABXITS","ITS")
CPH <- coxph(Surv(rem2tar, rem2) ~ tx*factor(paritycat) + dimdo + drylength + strata(cm1) + cluster(drydate), data=data)
car::Anova(CPH)cull.lact.parity.p <- car::Anova(CPH) %>% tail(1)
out <- emmeans(CPH,pairwise ~ tx | paritycat,adjust = "none",type = "response") %>% confint
out## $emmeans
## paritycat = 2:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 1.35 0.651 Inf 0.522 3.48
## ITS 1.38 0.726 Inf 0.491 3.87
##
## paritycat = 3:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 2.60 1.221 Inf 1.037 6.53
## ITS 2.43 1.210 Inf 0.917 6.45
##
## paritycat = 4:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 2.40 1.208 Inf 0.898 6.44
## ITS 2.26 1.149 Inf 0.834 6.12
##
## paritycat = 5+:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 3.51 1.875 Inf 1.231 10.00
## ITS 4.09 2.070 Inf 1.519 11.03
##
## Results are averaged over the levels of: cm1
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## paritycat = 2:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 0.978 0.166 Inf 0.701 1.36
##
## paritycat = 3:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 1.070 0.250 Inf 0.676 1.69
##
## paritycat = 4:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 1.064 0.209 Inf 0.724 1.56
##
## paritycat = 5+:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 0.857 0.180 Inf 0.568 1.29
##
## Results are averaged over the levels of: cm1
## Confidence level used: 0.95
## Intervals are back-transformed from the log scaletab <- data.frame(
outcome = "Removal during 1-200 DIM",
stratatype = "Parity",
strata = out[["contrasts"]][["paritycat"]],
contrast = out[["contrasts"]][["contrast"]],
est = out[["contrasts"]][["ratio"]],
lcl = out[["contrasts"]][["asymp.LCL"]],
ucl = out[["contrasts"]][["asymp.UCL"]]
)
tab <- tab %>% mutate(
estimate = paste0(round(est,2),
" (",
round(lcl,2),
" to ",
round(ucl,2),
")"),
est = NULL,
lcl = NULL,
ucl = NULL
)
tab.cull.lact.parity <- tab
tab.cull.lact.paritydata$tx <- fct_relevel(data$tx,"ITS","ABXITS")Clinical mastitis history
data$tx <- fct_relevel(data$tx,"ABXITS","ITS")
CPH <- coxph(Surv(rem2tar, rem2) ~ tx*cm1 + factor(paritycat) + dimdo + drylength + cluster(drydate), data=data)
car::Anova(CPH)cull.lact.cm1.p <- car::Anova(CPH) %>% tail(1)
out <- emmeans(CPH,pairwise ~ tx | cm1,adjust = "none",type = "response") %>% confint
out## $emmeans
## cm1 = 0:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 2.38 1.11 Inf 0.956 5.93
## ITS 2.32 1.16 Inf 0.871 6.18
##
## cm1 = 1:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 3.53 2.08 Inf 1.116 11.18
## ITS 4.31 2.04 Inf 1.698 10.92
##
## Results are averaged over the levels of: paritycat
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## cm1 = 0:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 1.03 0.124 Inf 0.809 1.30
##
## cm1 = 1:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 0.82 0.256 Inf 0.445 1.51
##
## Results are averaged over the levels of: paritycat
## Confidence level used: 0.95
## Intervals are back-transformed from the log scaletab <- data.frame(
outcome = "Removal during 1-200 DIM",
stratatype = "Clinical mastitis history",
strata = out[["contrasts"]][["cm1"]],
contrast = out[["contrasts"]][["contrast"]],
est = out[["contrasts"]][["ratio"]],
lcl = out[["contrasts"]][["asymp.LCL"]],
ucl = out[["contrasts"]][["asymp.UCL"]]
)
tab <- tab %>% mutate(
estimate = paste0(round(est,2),
" (",
round(lcl,2),
" to ",
round(ucl,2),
")"),
est = NULL,
lcl = NULL,
ucl = NULL
)
tab.cull.lact.cm1 <- tab
tab.cull.lact.cm1data$tx <- fct_relevel(data$tx,"ITS","ABXITS")Dry period length
data$tx <- fct_relevel(data$tx,"ABXITS","ITS")
CPH <- coxph(Surv(rem2tar, rem2) ~ tx*dplength3 + cm1 + factor(paritycat) + dimdo + cluster(drydate), data=data)
car::Anova(CPH)cull.lact.dplength.p <- car::Anova(CPH) %>% tail(1)
out <- emmeans(CPH,pairwise ~ tx | dplength3,adjust = "none",type = "response") %>% confint
out## $emmeans
## dplength3 = Low dplength:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 6.30 1.89 Inf 3.49 11.36
## ITS 5.20 2.04 Inf 2.41 11.21
##
## dplength3 = Mid dplength:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 4.18 1.59 Inf 1.98 8.80
## ITS 4.86 1.45 Inf 2.71 8.71
##
## dplength3 = High dplength:
## tx response SE df asymp.LCL asymp.UCL
## ABXITS 4.26 1.47 Inf 2.17 8.37
## ITS 4.77 1.67 Inf 2.41 9.46
##
## Results are averaged over the levels of: cm1, paritycat
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## dplength3 = Low dplength:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 1.211 0.261 Inf 0.793 1.85
##
## dplength3 = Mid dplength:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 0.859 0.171 Inf 0.582 1.27
##
## dplength3 = High dplength:
## contrast ratio SE df asymp.LCL asymp.UCL
## ABXITS / ITS 0.894 0.158 Inf 0.632 1.27
##
## Results are averaged over the levels of: cm1, paritycat
## Confidence level used: 0.95
## Intervals are back-transformed from the log scaletab <- data.frame(
outcome = "Removal during 1-200 DIM",
stratatype = "Dry period length",
strata = out[["contrasts"]][["dplength3"]],
contrast = out[["contrasts"]][["contrast"]],
est = out[["contrasts"]][["ratio"]],
lcl = out[["contrasts"]][["asymp.LCL"]],
ucl = out[["contrasts"]][["asymp.UCL"]]
)
tab <- tab %>% mutate(
estimate = paste0(round(est,2),
" (",
round(lcl,2),
" to ",
round(ucl,2),
")"),
est = NULL,
lcl = NULL,
ucl = NULL
)
tab.cull.lact.dplength <- tab
tab.cull.lact.dplengthdata$tx <- fct_relevel(data$tx,"ITS","ABXITS")cull.lact.hr <- rbind(tab.cull.lact,
tab.cull.lact.parity,
tab.cull.lact.dplength,
tab.cull.lact.cm1)
cull.lact.p <- rbind(
cull.lact.parity.p,
cull.lact.cm1.p,
cull.lact.dplength.p
)
cull.lact.p$outcome <- "Removal during 1-200 DIM"p <- rbind(cm.p,
cull.dp.p,
cull.lact.p) %>% as.data.frame()
p$`Pr(>Chisq)` <- round(p$`Pr(>Chisq)`,3)hr <- rbind(cm.hr,
cull.dp.hr,
cull.lact.hr)
freq <- merge(freq.tab,
hr,
all.x=T,
all.y=T,
by=c("outcome","stratatype","strata"))
p$stratatype <- c("Clinical mastitis history",
"Parity",
"Dry period length",
"Clinical mastitis history",
"Parity",
"Parity",
"Clinical mastitis history",
"Dry period length")
p$p <- p$`Pr(>Chisq)`
p <- p %>% select(outcome,stratatype,p)
freq2 <- merge(freq,
p,
all.x=T,
by=c("outcome","stratatype"))
freq2 <- freq2 %>% arrange(outcome,stratatype,strata,tx)
write.csv(freq2,"freq.table.csv")