Notes on dataset

Datasets

1. calf.survival

The objective of this dataset was to evaluate if pre-weaning scours at weeks 1, 2, 3, and 4+ were associated with death, culling and first calving.

Exposures

scour.survival - In the calf.survival dataset, calves with scours were retrospectively enrolled at the time of their first case of diarrhoea. For each enrolled calf with diarrhea, two unexposed calves (i.e. calves without diarrhea) were enrolled, matched for age (±3d) and date of birth (±15d). If unexposed calves experienced a case of diarrhea between enrollment and weaning, then they were right censored and re-enrolled as an exposed calf, and two more unexposed animals were enrolled. We decided to match unexposed animals with cases in order to control for group-level confounders (e.g. weather, season).

agecat - The age of enrollment was classified into weeks (w1 = week 1, w2 = week 2, w3 = week 3, week4+ = weeks ≥4). For this analysis we will stratify scour.survival into weeks so that we can explore if age at scours affects calf survival/calving.

Outcomes

dead - Records of death from any cause were used to classify if subjects had died between enrollment and censoring. Subjects were right censored if they were sold, calved (1st calving), changed exposure status, or at their last health record entry (i.e. evidence that the calf was still in the herd). Use survival analysis for this outcome. Covariate = birth.season

Death

The “dead” outcome was originally split into two timer periods after finding that the proportional hazards assumption was violated in a number of models. The first time period (dead30) included the first 30 days immediately after enrollment, which was used to estimate the short-term impact of diarrhea on survival. The second time period (dead31_to_c1) was used to estimate the effect of diarrhea on survival from 30 days after enrollment until first calving.

The dead outcome was subsequently split into 3 time periods after finding that the 2-way split continued to violate the proportional hazards assumption. The 3 time periods were 0-20 days, 21-100days and 100+ days after enrollment.

c1 - Calving records were used to classify if subjects had calved between enrollment and censoring. Subjects were right censored if they were sold, died, or at at their last health record entry. Use survival analysis for this outcome. Covariate = birth.season

sold - Records were used to classify if subjects had been culled (i.e. “SOLD” events) from the herd between enrollment and censoring. Subjects were right censored if they died, calved (1st calving), or at their last health record entry. Use survival analysis for this outcome. Covariate = birth.season

2. calf.lact

The objective of the calf.lact dataset was to evaluate if pre-weaning scours at weeks 1, 2, 3 and 4+ were associated with health and productivity during the first lactation. In order to be included in the calf.lact dataset, subjects had to have a recorded calving date.

Exposures

scour - Subjects that were enrolled as part of the scour.survival dataset were classified as having had at least 1 case of scours between birth and weaning (scour==1) or having no cases of scours between birth and weaning (scours==0).

scour_w1 - subjects that had experience a case of diarrhea at days 1-7 were classified as scour_w1==1. scour_w2 - subjects that had experience a case of diarrhea at days 8-14 were classified as scour_w2==1. scour_w3 - subjects that had experience a case of diarrhea at days 15-21 were classified as scour_w3==1. scour_w4 - subjects that had experience a case of diarrhea between 22 days old and weaning were classified as scour_w4+==1.

scour_w2_prev - subjects that had experienced a case of diarrhea prior to day 8 were classified as scour_w2_prev==1 scour_w3_prev - subjects that had experienced a case of diarrhea prior to day 15 were classified as scour_w3_prev==1 scour_w4_prev - subjects that had experienced a case of diarrhea prior to day 22 were classified as scour_w4_prev==1

Outcomes

c1age - age of first calving for cows that calved. Use linear regression for this outcome, Covariate = enrollment.season. UPDATE: DON’T EVALUTE THIS OUTCOME

preg200 - subjects were classified as having experienced a pregnancy event between first calving and censoring. Subjects were right censored if they died, were sold or reached 200 days in milk. Use survival analysis for this outcome. Covariate = calving.season

rem300 - subjects were classified as having experienced a removal event (sold or died) between first calving and censoring. Subjects were right censored if they reached 300 dim without having experienced the event. Use survival analysis for this outcome. Covariate = calving.season

mast120 - subjects were classified as having experienced a mastitis event between calving and censoring. Subjects were censored if they died, were sold or reached 120 days in milk. Use survival analysis for this outcome. Covariate = calving.season

met10 - subjects were classified as having experienced a metritis event between calving and censoring. Subjects were censored if they died, were sold or reached 10 days in milk. Use survival analysis for this outcome. Covariate = calving.season

rfm10 - subjects were classified as having experienced a retained placenta event between calving and censoring. Subjects were censored if they died, were sold or reached 10 days in milk. Use survival analysis for this outcome. Covariate = calving.season

maxlscc - the maximum logSCC value from the first lactation of each subject was determined using test-day records. Use linear regression for this outcome. Covariates = calving.season

me305 - the estimated 305-d mature equivalent milk production Use linear regression for this outcome. Covariates = calving.season

Exposures

scour - Subjects that were enrolled as part of the scour.survival dataset were classified as having had at least 1 case of scours between birth and weaning (scours==1) or having no cases of scours between birth and weaning (scours==0).

scour_w1 - subjects that had experience a case of diarrhea at days 1-7 were classified as scour_w1==1. scour_w2 - subjects that had experience a case of diarrhea at days 8-14 were classified as scour_w2==1. scour_w3 - subjects that had experience a case of diarrhea at days 15-21 were classified as scour_w3==1. scour_w4 - subjects that had experience a case of diarrhea between 22 days old and weaning were classified as scour_w4+==1.

scour_w2_prev - subjects that had experienced a case of diarrhoea prior to day 8 were classified as scour_w2_prev==1 scour_w3_prev - subjects that had experienced a case of diarrhoea prior to day 15 were classified as scour_w3_prev==1 scour_w4_prev - subjects that had experienced a case of diarrhoea prior to day 22 were classified as scour_w4_prev==1

A guide to using RMarkdown: https://rmarkdown.rstudio.com/lesson-1.html

Load packages

#Load packages

library(tidyverse)
library(epiR)
library(broom)
library(survival)
library(survminer)
library(table1)
library(stats)
library(broom)
library(infer)
library(emmeans)
library(stringr)

Import data

igg <- readRDS("colostrum.rds")
calf.lact <- readRDS("calf.lact.update.rds") %>% distinct(id, .keep_all=T)
calf.test <- readRDS("calf.test.rds")
calf.survival <- readRDS("calf.survival.rds")

calf.lact$scours_cases <- calf.lact$scour_w1 + calf.lact$scour_w2 + calf.lact$scour_w3 + calf.lact$scour_w4

Restricting dataset to cows in the herd long enough to calve in

calf.survival <- calf.survival %>% filter(bdate < as.Date("2020-06-23"))

Checking out the datasets

calf.survival

print(summarytools::dfSummary(calf.survival,valid.col=FALSE, graph.magnif=0.8, style="grid"), method = "render")

Data Frame Summary

calf.survival

Dimensions: 11135 x 24
Duplicates: 2
No Variable Stats / Values Freqs (% of Valid) Graph Missing
1 id [integer]
Mean (sd) : 44077.2 (8353)
min ≤ med ≤ max:
29712 ≤ 43547 ≤ 60243
IQR (CV) : 11509.5 (0.2)
 9833 distinct values 0 (0.0%)
2 bdate [Date]
min : 2015-01-22
med : 2018-05-12
max : 2020-06-22
range : 5y 5m 0d
 1852 distinct values 0 (0.0%)
3 age_cont [numeric]
Mean (sd) : 11 (7.3)
min ≤ med ≤ max:
0 ≤ 10 ≤ 71
IQR (CV) : 6 (0.7)
 64 distinct values 0 (0.0%)
4 agecat [factor]
1. w1
2. w2
3. w3
4. w4+
3205(28.8%)
5990(53.8%)
1207(10.8%)
733(6.6%)
 0 (0.0%)
5 birth.season [character]
1. Autumn
2. Spring
3. Summer
4. Winter
3583(32.2%)
2499(22.4%)
2305(20.7%)
2748(24.7%)
 0 (0.0%)
6 exitage [difftime]
min : 1
med : 868
max : 2594
units : days
 1761 distinct values 4967 (44.6%)
7 weanage [numeric]
Mean (sd) : 66.8 (6.9)
min ≤ med ≤ max:
41 ≤ 66 ≤ 99
IQR (CV) : 9 (0.1)
 50 distinct values 304 (2.7%)
8 scour.survival [numeric]
Min : 0
Mean : 0.3
Max : 1
0:7449(66.9%)
1:3686(33.1%)
 0 (0.0%)
9 dead [numeric]
Min : 0
Mean : 0.1
Max : 1
0:10385(93.3%)
1:750(6.7%)
 0 (0.0%)
10 deadtar [numeric]
Mean (sd) : 550.4 (268.4)
min ≤ med ≤ max:
0 ≤ 669 ≤ 1373
IQR (CV) : 187.5 (0.5)
 858 distinct values 0 (0.0%)
11 sold [numeric]
Min : 0
Mean : 0.1
Max : 1
0:9906(89.0%)
1:1229(11.0%)
 0 (0.0%)
12 soldtar [numeric]
Mean (sd) : 550.4 (268.4)
min ≤ med ≤ max:
0 ≤ 669 ≤ 1373
IQR (CV) : 187.5 (0.5)
 858 distinct values 0 (0.0%)
13 c1 [numeric]
Min : 0
Mean : 0.6
Max : 1
0:4010(36.0%)
1:7125(64.0%)
 0 (0.0%)
14 c1tar [numeric]
Mean (sd) : 550.4 (268.4)
min ≤ med ≤ max:
0 ≤ 669 ≤ 1373
IQR (CV) : 187.5 (0.5)
 858 distinct values 0 (0.0%)
15 dead30 [numeric]
Min : 0
Mean : 0
Max : 1
0:10921(98.1%)
1:214(1.9%)
 0 (0.0%)
16 dead30tar [numeric]
Mean (sd) : 27.1 (8.1)
min ≤ med ≤ max:
0 ≤ 30 ≤ 30
IQR (CV) : 0 (0.3)
 31 distinct values 0 (0.0%)
17 dead31_to_c1 [numeric]
Min : 0
Mean : 0.1
Max : 1
0:9246(94.5%)
1:536(5.5%)
 1353 (12.2%)
18 dead31_to_c1tar [numeric]
Mean (sd) : 595.7 (187.8)
min ≤ med ≤ max:
0 ≤ 645 ≤ 1343
IQR (CV) : 64 (0.3)
 828 distinct values 1353 (12.2%)
19 dead20 [numeric]
Min : 0
Mean : 0
Max : 1
0:10937(98.2%)
1:198(1.8%)
 0 (0.0%)
20 dead20tar [numeric]
Mean (sd) : 18.3 (5)
min ≤ med ≤ max:
0 ≤ 20 ≤ 20
IQR (CV) : 0 (0.3)
 21 distinct values 0 (0.0%)
21 dead21_to_100 [numeric]
Min : 0
Mean : 0
Max : 1
0:9652(98.2%)
1:179(1.8%)
 1304 (11.7%)
22 dead21_to_100tar [numeric]
Mean (sd) : 78.8 (8.2)
min ≤ med ≤ max:
0 ≤ 80 ≤ 80
IQR (CV) : 0 (0.1)
 76 distinct values 1304 (11.7%)
23 dead101_to_c1 [numeric]
Min : 0
Mean : 0
Max : 1
0:9184(96.1%)
1:373(3.9%)
 1578 (14.2%)
24 dead101_to_c1tar [numeric]
Mean (sd) : 538.9 (168.8)
min ≤ med ≤ max:
1 ≤ 576 ≤ 1273
IQR (CV) : 61 (0.3)
 762 distinct values 1578 (14.2%)

Generated by summarytools 1.0.1 (R version 4.3.3)
2024-05-01

calf.lact

print(summarytools::dfSummary(calf.lact,valid.col=FALSE, graph.magnif=0.8, style="grid"), method = "render")

Data Frame Summary

calf.lact

Dimensions: 7047 x 26
Duplicates: 0
No Variable Stats / Values Freqs (% of Valid) Graph Missing
1 id [integer]
Mean (sd) : 41658.9 (6791.1)
min ≤ med ≤ max:
29712 ≤ 42296 ≤ 60225
IQR (CV) : 11573.5 (0.2)
 7047 distinct values 0 (0.0%)
2 bdate [Date]
min : 2015-01-22
med : 2018-01-01
max : 2020-06-23
range : 5y 5m 1d
 1758 distinct values 0 (0.0%)
3 c1age [numeric]
Mean (sd) : 714.6 (55.2)
min ≤ med ≤ max:
601 ≤ 696 ≤ 1104
IQR (CV) : 53 (0.1)
 334 distinct values 0 (0.0%)
4 enroll.season [character]
1. Autumn
2. Spring
3. Summer
4. Winter
1771(25.1%)
1930(27.4%)
1631(23.1%)
1715(24.3%)
 0 (0.0%)
5 calving.season [character]
1. Autumn
2. Spring
3. Summer
4. Winter
1627(23.1%)
1789(25.4%)
1793(25.4%)
1838(26.1%)
 0 (0.0%)
6 scour [numeric]
Min : 0
Mean : 0.4
Max : 1
0:4471(63.4%)
1:2576(36.6%)
 0 (0.0%)
7 scour_w1 [numeric]
Min : 0
Mean : 0.1
Max : 1
0:6539(92.8%)
1:508(7.2%)
 0 (0.0%)
8 scour_w2 [numeric]
Min : 0
Mean : 0.2
Max : 1
0:5527(78.4%)
1:1520(21.6%)
 0 (0.0%)
9 scour_w3 [numeric]
Min : 0
Mean : 0.1
Max : 1
0:6573(93.3%)
1:474(6.7%)
 0 (0.0%)
10 scour_w4 [numeric]
Min : 0
Mean : 0
Max : 1
0:6764(96.0%)
1:283(4.0%)
 0 (0.0%)
11 scour_w2_prev [numeric]
Min : 0
Mean : 0.1
Max : 1
0:6539(92.8%)
1:508(7.2%)
 0 (0.0%)
12 scour_w3_prev [numeric]
Min : 0
Mean : 0.3
Max : 1
0:5072(72.0%)
1:1975(28.0%)
 0 (0.0%)
13 scour_w4_prev [numeric]
Min : 0
Mean : 0.3
Max : 1
0:4680(66.4%)
1:2367(33.6%)
 0 (0.0%)
14 preg200 [numeric]
Min : 0
Mean : 0.7
Max : 1
0:2047(29.0%)
1:5000(71.0%)
 0 (0.0%)
15 preg200tar [numeric]
Mean (sd) : 124.3 (50.2)
min ≤ med ≤ max:
0 ≤ 116 ≤ 200
IQR (CV) : 67 (0.4)
 201 distinct values 0 (0.0%)
16 rem300 [numeric]
Min : 0
Mean : 0.2
Max : 1
0:5824(82.6%)
1:1223(17.4%)
 0 (0.0%)
17 rem300tar [numeric]
Mean (sd) : 243.4 (95)
min ≤ med ≤ max:
0 ≤ 300 ≤ 300
IQR (CV) : 96 (0.4)
 301 distinct values 0 (0.0%)
18 mast120 [numeric]
Min : 0
Mean : 0
Max : 1
0:6867(97.4%)
1:180(2.6%)
 0 (0.0%)
19 mast120tar [numeric]
Mean (sd) : 107.4 (29.7)
min ≤ med ≤ max:
0 ≤ 120 ≤ 120
IQR (CV) : 0 (0.3)
 121 distinct values 0 (0.0%)
20 met10 [numeric]
Min : 0
Mean : 0
Max : 1
0:6709(95.2%)
1:338(4.8%)
 0 (0.0%)
21 met10tar [numeric]
Mean (sd) : 9.7 (1.3)
min ≤ med ≤ max:
0 ≤ 10 ≤ 10
IQR (CV) : 0 (0.1)
 11 distinct values 0 (0.0%)
22 rfm10 [numeric]
Min : 0
Mean : 0
Max : 1
0:6874(97.5%)
1:173(2.5%)
 0 (0.0%)
23 rfm10tar [numeric]
Mean (sd) : 9.8 (1.4)
min ≤ med ≤ max:
0 ≤ 10 ≤ 10
IQR (CV) : 0 (0.1)
 11 distinct values 0 (0.0%)
24 maxlscc [numeric]
Mean (sd) : 5.5 (1.2)
min ≤ med ≤ max:
1.6 ≤ 5.4 ≤ 9.2
IQR (CV) : 1.6 (0.2)
 1557 distinct values 891 (12.6%)
25 me305 [integer]
Mean (sd) : 13451.8 (2888.2)
min ≤ med ≤ max:
0 ≤ 13700 ≤ 23560
IQR (CV) : 3612.5 (0.2)
 1303 distinct values 891 (12.6%)
26 scours_cases [numeric]
Mean (sd) : 0.4 (0.5)
min ≤ med ≤ max:
0 ≤ 0 ≤ 3
IQR (CV) : 1 (1.4)
0:4471(63.4%)
1:2372(33.7%)
2:199(2.8%)
3:5(0.1%)
 0 (0.0%)

Generated by summarytools 1.0.1 (R version 4.3.3)
2024-05-01

Descriptive results

#?summarytools::dfSummary
calf.survival <- calf.survival %>% distinct(id,scour.survival, .keep_all = TRUE)

calf.exposure <- calf.survival %>% group_by(id) %>% summarise(
  times_enrolled_scour = sum(scour.survival==1),
  times_enrolled_no_scour = sum(scour.survival==0)
)

calf.exposure %>% summarise(
  enrolled_scour = sum(times_enrolled_no_scour==0 & times_enrolled_scour==1),
  enrolled_no_scour = sum(times_enrolled_no_scour==1 & times_enrolled_scour==0),
  enrolled_switched = sum(times_enrolled_no_scour==1 & times_enrolled_scour==1)
)
calf.survival <- calf.survival %>% mutate(
  c1tar600 = case_when(
    c1tar < 590 ~ 0,
    c1tar >= 590 ~ c1tar - 590
  )
)

check <- calf.survival %>% filter(c1tar>500 & c1==1)
#hist(calf.survival$c1tar600,breaks="FD")

calf.desc <- calf.survival %>% arrange(scour.survival) %>% distinct(id, .keep_all = TRUE)

calf.lact <- calf.lact %>% mutate(
  multi_scours = case_when(
    scours_cases > 1 ~ 1,
    scours_cases <= 1 ~ 0
  )
) %>% left_join(calf.desc %>% select(id,weanage), by = "id")

Calves enrolled without diarrhea

summarytools::dfSummary(calf.desc %>% filter(scour.survival==0),valid.col=FALSE,
                        graph.magnif=0.8, style="grid",varnumbers=F,na.col=F,graph.col=F)

Calves enrolled with diarrhea

summarytools::dfSummary(calf.desc %>% filter(scour.survival==1),valid.col=FALSE,
                        graph.magnif=0.8, style="grid",varnumbers=F,na.col=F,graph.col=F)
print(summarytools::dfSummary(calf.lact %>% filter(scour==1),valid.col=FALSE, graph.magnif=0.8, style="grid"), method = "render")

Data Frame Summary

calf.lact

Dimensions: 2576 x 28
Duplicates: 0
No Variable Stats / Values Freqs (% of Valid) Graph Missing
1 id [integer]
Mean (sd) : 41563 (6722.2)
min ≤ med ≤ max:
29712 ≤ 42224 ≤ 60175
IQR (CV) : 10701 (0.2)
 2576 distinct values 0 (0.0%)
2 bdate [Date]
min : 2015-01-22
med : 2017-12-22
max : 2020-05-22
range : 5y 4m 0d
 1218 distinct values 0 (0.0%)
3 c1age [numeric]
Mean (sd) : 716.6 (56.2)
min ≤ med ≤ max:
601 ≤ 698 ≤ 1094
IQR (CV) : 55 (0.1)
 272 distinct values 0 (0.0%)
4 enroll.season [character]
1. Autumn
2. Spring
3. Summer
4. Winter
643(25.0%)
745(28.9%)
608(23.6%)
580(22.5%)
 0 (0.0%)
5 calving.season [character]
1. Autumn
2. Spring
3. Summer
4. Winter
575(22.3%)
687(26.7%)
669(26.0%)
645(25.0%)
 0 (0.0%)
6 scour [numeric] 1 distinct value
1:2576(100.0%)
 0 (0.0%)
7 scour_w1 [numeric]
Min : 0
Mean : 0.2
Max : 1
0:2068(80.3%)
1:508(19.7%)
 0 (0.0%)
8 scour_w2 [numeric]
Min : 0
Mean : 0.6
Max : 1
0:1056(41.0%)
1:1520(59.0%)
 0 (0.0%)
9 scour_w3 [numeric]
Min : 0
Mean : 0.2
Max : 1
0:2102(81.6%)
1:474(18.4%)
 0 (0.0%)
10 scour_w4 [numeric]
Min : 0
Mean : 0.1
Max : 1
0:2293(89.0%)
1:283(11.0%)
 0 (0.0%)
11 scour_w2_prev [numeric]
Min : 0
Mean : 0.2
Max : 1
0:2068(80.3%)
1:508(19.7%)
 0 (0.0%)
12 scour_w3_prev [numeric]
Min : 0
Mean : 0.8
Max : 1
0:601(23.3%)
1:1975(76.7%)
 0 (0.0%)
13 scour_w4_prev [numeric]
Min : 0
Mean : 0.9
Max : 1
0:209(8.1%)
1:2367(91.9%)
 0 (0.0%)
14 preg200 [numeric]
Min : 0
Mean : 0.7
Max : 1
0:782(30.4%)
1:1794(69.6%)
 0 (0.0%)
15 preg200tar [numeric]
Mean (sd) : 124.2 (50.3)
min ≤ med ≤ max:
1 ≤ 114 ≤ 200
IQR (CV) : 66 (0.4)
 199 distinct values 0 (0.0%)
16 rem300 [numeric]
Min : 0
Mean : 0.2
Max : 1
0:2098(81.4%)
1:478(18.6%)
 0 (0.0%)
17 rem300tar [numeric]
Mean (sd) : 241.4 (96.3)
min ≤ med ≤ max:
1 ≤ 300 ≤ 300
IQR (CV) : 104.2 (0.4)
 276 distinct values 0 (0.0%)
18 mast120 [numeric]
Min : 0
Mean : 0
Max : 1
0:2510(97.4%)
1:66(2.6%)
 0 (0.0%)
19 mast120tar [numeric]
Mean (sd) : 107 (29.8)
min ≤ med ≤ max:
0 ≤ 120 ≤ 120
IQR (CV) : 0 (0.3)
 118 distinct values 0 (0.0%)
20 met10 [numeric]
Min : 0
Mean : 0
Max : 1
0:2448(95.0%)
1:128(5.0%)
 0 (0.0%)
21 met10tar [numeric]
Mean (sd) : 9.7 (1.3)
min ≤ med ≤ max:
1 ≤ 10 ≤ 10
IQR (CV) : 0 (0.1)
1:5(0.2%)
2:13(0.5%)
3:13(0.5%)
4:52(2.0%)
5:15(0.6%)
6:13(0.5%)
7:15(0.6%)
8:7(0.3%)
9:16(0.6%)
10:2427(94.2%)
 0 (0.0%)
22 rfm10 [numeric]
Min : 0
Mean : 0
Max : 1
0:2517(97.7%)
1:59(2.3%)
 0 (0.0%)
23 rfm10tar [numeric]
Mean (sd) : 9.8 (1.3)
min ≤ med ≤ max:
0 ≤ 10 ≤ 10
IQR (CV) : 0 (0.1)
0:4(0.2%)
1:25(1.0%)
2:17(0.7%)
3:9(0.3%)
4:9(0.3%)
5:6(0.2%)
6:4(0.2%)
7:2(0.1%)
9:4(0.2%)
10:2496(96.9%)
 0 (0.0%)
24 maxlscc [numeric]
Mean (sd) : 5.5 (1.2)
min ≤ med ≤ max:
2.1 ≤ 5.4 ≤ 9.2
IQR (CV) : 1.6 (0.2)
 943 distinct values 340 (13.2%)
25 me305 [integer]
Mean (sd) : 13400.5 (2941.8)
min ≤ med ≤ max:
3400 ≤ 13680 ≤ 23560
IQR (CV) : 3672.5 (0.2)
 959 distinct values 340 (13.2%)
26 scours_cases [numeric]
Mean (sd) : 1.1 (0.3)
min ≤ med ≤ max:
1 ≤ 1 ≤ 3
IQR (CV) : 0 (0.3)
1:2372(92.1%)
2:199(7.7%)
3:5(0.2%)
 0 (0.0%)
27 multi_scours [numeric]
Min : 0
Mean : 0.1
Max : 1
0:2372(92.1%)
1:204(7.9%)
 0 (0.0%)
28 weanage [numeric]
Mean (sd) : 67.2 (7.2)
min ≤ med ≤ max:
42 ≤ 66.5 ≤ 99
IQR (CV) : 10 (0.1)
 48 distinct values 0 (0.0%)

Generated by summarytools 1.0.1 (R version 4.3.3)
2024-05-01

1. Survival analysis

Variables include: dead, dead20, dead21_to_100, dead101_to_c1, c1, mast120, rem300, met10, rfm10, preg200

Calculating rates

rate <- function(predictor,time.at.risk,outcome,group,dataframe){

  pred.label <- deparse(substitute(predictor))
  out.label <- deparse(substitute(outcome))
  predictor <-  enquo(predictor)
  outcome <-  enquo(outcome)
  tar <- enquo(time.at.risk)
  
  tab <- dataframe %>% filter(!is.na(!!outcome)) %>% group_by(!!predictor) %>% summarise(
  n.case = sum(!!outcome==1),
  time.at.risk = sum(!!tar),
  incidence.rate = n.case / time.at.risk,
  incidence.rate.10000 = round(incidence.rate*10000,2),
)

colnames(tab)[1] <- "exposed"
tab$exposed[tab$exposed==0] <- "Unexposed"
tab$exposed[tab$exposed==1] <- "Exposed"

tab$incidence.rate <- NULL
tab <- tab %>%
  pivot_wider(
    names_from = exposed,
    values_from = c(n.case, time.at.risk,incidence.rate.10000)
  )

tab$exposed.cases <- paste0(tab$n.case_Exposed, " cases in ",tab$time.at.risk_Exposed," days")
tab$unexposed.cases <- paste0(tab$n.case_Unexposed, " cases in ",tab$time.at.risk_Unexposed," days")
tab$exposed.incidence.rate <- paste0(tab$incidence.rate.10000_Exposed," cases per 10,000 days at risk")
tab$unexposed.incidence.rate <- paste0(tab$incidence.rate.10000_Unexposed," cases per 10,000 days at risk")
tab$incidence.rate.ratio <- round(tab$incidence.rate.10000_Exposed/tab$incidence.rate.10000_Unexposed,2)
tab$outcome <- out.label
tab$predictor <- pred.label
tab$week <- group
tab %>% select(outcome,predictor,week,exposed.cases,exposed.incidence.rate,unexposed.cases,unexposed.incidence.rate,incidence.rate.ratio,
               n.case_Exposed,time.at.risk_Exposed,incidence.rate.10000_Exposed,n.case_Unexposed,time.at.risk_Unexposed,incidence.rate.10000_Unexposed)
}

Rates in the “calf.survival” dataset

Example - death between enrollment and first calving

dead.scour.rate <- rate(predictor = scour.survival,
         outcome = dead,
         time.at.risk = deadtar,
         group = "All weeks",
         dataframe = calf.survival)

dead.scour.w1.rate <- rate(predictor = scour.survival,
         outcome = dead,
         time.at.risk = deadtar,
         group = "w1",
         dataframe = calf.survival %>% filter(agecat=="w1"))

dead.scour.w2.rate <- rate(predictor = scour.survival,
         outcome = dead,
         time.at.risk = deadtar,
         group = "w2",
         dataframe = calf.survival %>% filter(agecat=="w2"))

dead.scour.w3.rate <- rate(predictor = scour.survival,
         outcome = dead,
         time.at.risk = deadtar,
         group = "w3",
         dataframe = calf.survival %>% filter(agecat=="w3"))

dead.scour.w4.rate <- rate(predictor = scour.survival,
         outcome = dead,
         time.at.risk = deadtar,
         group = "w4+",
         dataframe = calf.survival %>% filter(agecat=="w4+"))

dead.scour.rates.all <- rbind(dead.scour.rate,
                              dead.scour.w1.rate,
                              dead.scour.w2.rate,
                              dead.scour.w3.rate,
                              dead.scour.w4.rate)

dead.scour.rates.all
write.csv(dead.scour.rates.all,'dead.scour.rates.all.csv')

Death enrollment - 20 days

dead20.scour.rate <- rate(predictor = scour.survival,
         outcome = dead20,
         time.at.risk = dead20tar,
         group = "All weeks",
         dataframe = calf.survival)

dead20.scour.w1.rate <- rate(predictor = scour.survival,
         outcome = dead20,
         time.at.risk = dead20tar,
         group = "w1",
         dataframe = calf.survival %>% filter(agecat=="w1"))

dead20.scour.w2.rate <- rate(predictor = scour.survival,
         outcome = dead20,
         time.at.risk = dead20tar,
         group = "w2",
         dataframe = calf.survival %>% filter(agecat=="w2"))

dead20.scour.w3.rate <- rate(predictor = scour.survival,
         outcome = dead20,
         time.at.risk = dead20tar,
         group = "w3",
         dataframe = calf.survival %>% filter(agecat=="w3"))

dead20.scour.w4.rate <- rate(predictor = scour.survival,
         outcome = dead20,
         time.at.risk = dead20tar,
         group = "w4+",
         dataframe = calf.survival %>% filter(agecat=="w4+"))

dead20.scour.rates.all <- rbind(dead20.scour.rate,
                              dead20.scour.w1.rate,
                              dead20.scour.w2.rate,
                              dead20.scour.w3.rate,
                              dead20.scour.w4.rate)

dead20.scour.rates.all

Death 21 days - 100

dead21_to_100.scour.rate <- rate(predictor = scour.survival,
         outcome = dead21_to_100,
         time.at.risk = dead21_to_100tar,
         group = "All weeks",
         dataframe = calf.survival)

dead21_to_100.scour.w1.rate <- rate(predictor = scour.survival,
         outcome = dead21_to_100,
         time.at.risk = dead21_to_100tar,
         group = "w1",
         dataframe = calf.survival %>% filter(agecat=="w1"))

dead21_to_100.scour.w2.rate <- rate(predictor = scour.survival,
         outcome = dead21_to_100,
         time.at.risk = dead21_to_100tar,
         group = "w2",
         dataframe = calf.survival %>% filter(agecat=="w2"))

dead21_to_100.scour.w3.rate <- rate(predictor = scour.survival,
         outcome = dead21_to_100,
         time.at.risk = dead21_to_100tar,
         group = "w3",
         dataframe = calf.survival %>% filter(agecat=="w3"))

dead21_to_100.scour.w4.rate <- rate(predictor = scour.survival,
         outcome = dead21_to_100,
         time.at.risk = dead21_to_100tar,
         group = "w4+",
         dataframe = calf.survival %>% filter(agecat=="w4+"))

dead21_to_100.scour.rates.all <- rbind(dead21_to_100.scour.rate,
                              dead21_to_100.scour.w1.rate,
                              dead21_to_100.scour.w2.rate,
                              dead21_to_100.scour.w3.rate,
                              dead21_to_100.scour.w4.rate)

dead21_to_100.scour.rates.all

Death 101 days - calving

dead101_to_c1.scour.rate <- rate(predictor = scour.survival,
         outcome = dead101_to_c1,
         time.at.risk = dead101_to_c1tar,
         group = "All weeks",
         dataframe = calf.survival)

dead101_to_c1.scour.w1.rate <- rate(predictor = scour.survival,
         outcome = dead101_to_c1,
         time.at.risk = dead101_to_c1tar,
         group = "w1",
         dataframe = calf.survival %>% filter(agecat=="w1"))

dead101_to_c1.scour.w2.rate <- rate(predictor = scour.survival,
         outcome = dead101_to_c1,
         time.at.risk = dead101_to_c1tar,
         group = "w2",
         dataframe = calf.survival %>% filter(agecat=="w2"))

dead101_to_c1.scour.w3.rate <- rate(predictor = scour.survival,
         outcome = dead101_to_c1,
         time.at.risk = dead101_to_c1tar,
         group = "w3",
         dataframe = calf.survival %>% filter(agecat=="w3"))

dead101_to_c1.scour.w4.rate <- rate(predictor = scour.survival,
         outcome = dead101_to_c1,
         time.at.risk = dead101_to_c1tar,
         group = "w4+",
         dataframe = calf.survival %>% filter(agecat=="w4+"))

dead101_to_c1.scour.rates.all <- rbind(dead101_to_c1.scour.rate,
                              dead101_to_c1.scour.w1.rate,
                              dead101_to_c1.scour.w2.rate,
                              dead101_to_c1.scour.w3.rate,
                              dead101_to_c1.scour.w4.rate)

dead101_to_c1.scour.rates.all

c1

c1.scour.rate <- rate(predictor = scour.survival,
         outcome = c1,
         time.at.risk = c1tar600,
         group = "All weeks",
         dataframe = calf.survival)

c1.scour.w1.rate <- rate(predictor = scour.survival,
         outcome = c1,
         time.at.risk = c1tar600,
         group = "w1",
         dataframe = calf.survival %>% filter(agecat=="w1"))

c1.scour.w2.rate <- rate(predictor = scour.survival,
         outcome = c1,
         time.at.risk = c1tar600,
         group = "w2",
         dataframe = calf.survival %>% filter(agecat=="w2"))

c1.scour.w3.rate <- rate(predictor = scour.survival,
         outcome = c1,
         time.at.risk = c1tar600,
         group = "w3",
         dataframe = calf.survival %>% filter(agecat=="w3"))

c1.scour.w4.rate <- rate(predictor = scour.survival,
         outcome = c1,
         time.at.risk = c1tar600,
         group = "w4+",
         dataframe = calf.survival %>% filter(agecat=="w4+"))

c1.scour.rates.all <- rbind(c1.scour.rate,
                            c1.scour.w1.rate,
                            c1.scour.w2.rate,
                            c1.scour.w3.rate,
                            c1.scour.w4.rate)
c1.scour.rates.all

Sold

sold.scour.rate <- rate(predictor = scour.survival,
         outcome = sold,
         time.at.risk = soldtar,
         group = "All weeks",
         dataframe = calf.survival)

sold.scour.w1.rate <- rate(predictor = scour.survival,
         outcome = sold,
         time.at.risk = soldtar,
         group = "w1",
         dataframe = calf.survival %>% filter(agecat=="w1"))

sold.scour.w2.rate <- rate(predictor = scour.survival,
         outcome = sold,
         time.at.risk = soldtar,
         group = "w2",
         dataframe = calf.survival %>% filter(agecat=="w2"))

sold.scour.w3.rate <- rate(predictor = scour.survival,
         outcome = sold,
         time.at.risk = soldtar,
         group = "w3",
         dataframe = calf.survival %>% filter(agecat=="w3"))

sold.scour.w4.rate <- rate(predictor = scour.survival,
         outcome = sold,
         time.at.risk = soldtar,
         group = "w4+",
         dataframe = calf.survival %>% filter(agecat=="w4+"))

sold.scour.rates.all <- rbind(sold.scour.rate,
                              sold.scour.w1.rate,
                              sold.scour.w2.rate,
                              sold.scour.w3.rate,
                              sold.scour.w4.rate)

sold.scour.rates.all

####Table combining all rates in calf.survival dataset

calf.survival.rates.all <- rbind(dead.scour.rates.all,
                                 dead20.scour.rates.all,
                                 dead21_to_100.scour.rates.all,
                                 dead101_to_c1.scour.rates.all,
                                 c1.scour.rates.all,
                                 sold.scour.rates.all)

calf.survival.rates.all
write.csv(calf.survival.rates.all,'calf.survival.rates.all.csv')

Rates in the “calf.lact” dataset

Eg. Removal from the herd between 1 and 300 days in milk.

rem300.scour.rate <- rate(predictor = scour,
         group = "All weeks",
         outcome = rem300,
         time.at.risk = rem300tar,
         dataframe = calf.lact)

rem300.scour.w1.rate <- rate(predictor = scour_w1,
         group = "w1",
         outcome = rem300,
         time.at.risk = rem300tar,
         dataframe = calf.lact)

rem300.scour.w2.rate <- rate(predictor = scour_w2,
         group = "w2",
         outcome = rem300,
         time.at.risk = rem300tar,
         dataframe = calf.lact)

rem300.scour.w3.rate <- rate(predictor = scour_w3,
         group = "w3",
         outcome = rem300,
         time.at.risk = rem300tar,
         dataframe = calf.lact)

rem300.scour.w4.rate <- rate(predictor = scour_w4,
         group = "w4+",
         outcome = rem300,
         time.at.risk = rem300tar,
         dataframe = calf.lact)

rem300.scours.rates.all <- rbind(rem300.scour.rate,
                                 rem300.scour.w1.rate,
                                 rem300.scour.w2.rate,
                                 rem300.scour.w3.rate,
                                 rem300.scour.w4.rate)

rem300.dead.scours.rates <- rbind(dead.scour.rates.all,
      rem300.scours.rates.all)

rem300.dead.scours.rates

Mastitis within first 120 DIM

mast120.scour.rate <- rate(predictor = scour,
         group = "All weeks",
         outcome = mast120,
         time.at.risk = mast120tar,
         dataframe = calf.lact)

mast120.scour.w1.rate <- rate(predictor = scour_w1,
         group = "w1",
         outcome = mast120,
         time.at.risk = mast120tar,
         dataframe = calf.lact)

mast120.scour.w2.rate <- rate(predictor = scour_w2,
         group = "w2",
         outcome = mast120,
         time.at.risk = mast120tar,
         dataframe = calf.lact)

mast120.scour.w3.rate <- rate(predictor = scour_w3,
         group = "w3",
         outcome = mast120,
         time.at.risk = mast120tar,
         dataframe = calf.lact)

mast120.scour.w4.rate <- rate(predictor = scour_w4,
         group = "w4+",
         outcome = mast120,
         time.at.risk = mast120tar,
         dataframe = calf.lact)

mast120.scours.rates.all <- rbind(mast120.scour.rate,
                                  mast120.scour.w1.rate,
                                  mast120.scour.w2.rate,
                                  mast120.scour.w3.rate,
                                  mast120.scour.w4.rate)
mast120.scours.rates.all

Metritis within first 10 DIM

met10.scour.rate <- rate(predictor = scour,
         group = "All weeks",
         outcome = met10,
         time.at.risk = met10tar,
         dataframe = calf.lact)

met10.scour.w1.rate <- rate(predictor = scour_w1,
         group = "w1",
         outcome = met10,
         time.at.risk = met10tar,
         dataframe = calf.lact)

met10.scour.w2.rate <- rate(predictor = scour_w2,
         group = "w2",
         outcome = met10,
         time.at.risk = met10tar,
         dataframe = calf.lact)

met10.scour.w3.rate <- rate(predictor = scour_w3,
         group = "w3",
         outcome = met10,
         time.at.risk = met10tar,
         dataframe = calf.lact)

met10.scour.w4.rate <- rate(predictor = scour_w4,
         group = "w4+",
         outcome = met10,
         time.at.risk = met10tar,
         dataframe = calf.lact)

met10.scours.rates.all <- rbind(met10.scour.rate,
                                met10.scour.w1.rate,
                                met10.scour.w2.rate,
                                met10.scour.w3.rate,
                                met10.scour.w4.rate)
met10.scours.rates.all

Retained foetal membranes during first 10 DIM

rfm10.scour.rate <- rate(predictor = scour,
         group = "All weeks",
         outcome = rfm10,
         time.at.risk = rfm10tar,
         dataframe = calf.lact)

rfm10.scour.w1.rate <- rate(predictor = scour_w1,
         group = "w1",
         outcome = rfm10,
         time.at.risk = rfm10tar,
         dataframe = calf.lact)

rfm10.scour.w2.rate <- rate(predictor = scour_w2,
         group = "w2",
         outcome = rfm10,
         time.at.risk = rfm10tar,
         dataframe = calf.lact)

rfm10.scour.w3.rate <- rate(predictor = scour_w3,
         group = "w3",
         outcome = rfm10,
         time.at.risk = rfm10tar,
         dataframe = calf.lact)

rfm10.scour.w4.rate <- rate(predictor = scour_w4,
         group = "w4+",
         outcome = rfm10,
         time.at.risk = rfm10tar,
         dataframe = calf.lact)

rfm10.scours.rates.all <- rbind(rfm10.scour.rate,
                                rfm10.scour.w1.rate,
                                rfm10.scour.w2.rate,
                                rfm10.scour.w3.rate,
                                rfm10.scour.w4.rate)

Pregnancy within 200 DIM

preg200.scour.rate <- rate(predictor = scour,
         group = "All weeks",
         outcome = preg200,
         time.at.risk = preg200tar,
         dataframe = calf.lact)

preg200.scour.w1.rate <- rate(predictor = scour_w1,
         group = "w1",
         outcome = preg200,
         time.at.risk = preg200tar,
         dataframe = calf.lact)

preg200.scour.w2.rate <- rate(predictor = scour_w2,
         group = "w2",
         outcome = preg200,
         time.at.risk = preg200tar,
         dataframe = calf.lact)

preg200.scour.w3.rate <- rate(predictor = scour_w3,
         group = "w3",
         outcome = preg200,
         time.at.risk = preg200tar,
         dataframe = calf.lact)

preg200.scour.w4.rate <- rate(predictor = scour_w4,
         group = "w4+",
         outcome = preg200,
         time.at.risk = preg200tar,
         dataframe = calf.lact)

preg200.scours.rates.all <- rbind(preg200.scour.rate,
                                  preg200.scour.w1.rate,
                                  preg200.scour.w2.rate,
                                  preg200.scour.w3.rate,
                                  preg200.scour.w4.rate)

####Table combining calf.lact rates

calf.lact.rates.all <- rbind(rem300.scours.rates.all,
                             mast120.scours.rates.all,
                             met10.scours.rates.all,
                             rfm10.scours.rates.all,
                             preg200.scours.rates.all)

calf.lact.rates.all
write.csv(calf.lact.rates.all,'calf.lact.rates.all.csv')

###Table combining survival analysis rates

survival.analysis.rates.all <- rbind(calf.survival.rates.all,
                                     calf.lact.rates.all)

survival.analysis.rates.all

Kaplan meier plots

These functions will create Kaplan Meier curves and estimate the incidence of the outcome at a particular point in time. For example, you may want to estimate the proportion of cows dead by 800 days.

KMplot <- function(survival.function,dataframe,xmin,xmax,ymax,x.label,y.label,title.label,time.break){
ggsurvplot(survival.function, data = dataframe,
           fun = "event",xlim=c(xmin,xmax),
           ylim=c(0,ymax),
           title = title.label, pval = F, censor=F, conf.int = F,risk.table = F,
           surv.plot.height = 5, tables.theme = theme_cleantable(),
           #risk.table.pos="in",
           ggtheme = theme_bw(base_family = "Times"), 
           xlab=x.label,
           ylab=y.label,
           palette = "Set1",risk.table.fontsize=3,break.time.by = time.break)
}



#Note that this function is not being used currently
risk.estimate <- function(survival.function,rcensor,group){
cmtable <- summary(survival.function, times = rcensor)
keep <- cmtable
cmtable <- data.frame(
  group = cmtable[["strata"]],
  days = cmtable[["time"]],
  incidence = round((1 - cmtable[["surv"]])*100,1)
)


cmtable
cmtable$predictor <- km[["call"]][["formula"]][[3]] %>% as.character
cmtable$label <- paste0(km[["call"]][["formula"]][[2]][[3]],".",cmtable$days[1],".",cmtable$group)
days <- cmtable$days[1]
cmtable <- cmtable %>% select(label,incidence) %>% pivot_wider(
   names_from = label,
   values_from = c(incidence))
colnames(cmtable) <- c("risk.unexposed","risk.exposed")
cmtable$predictor <- km[["call"]][["formula"]][[3]] %>% as.character
cmtable$week <- group
cmtable$outcome <- km[["call"]][["formula"]][[2]][[3]] %>% as.character
cmtable$risk.day <- days
cmtable %>% select(outcome,predictor,week,risk.day,risk.exposed,risk.unexposed)
}

KM plots using the calf.survival dataset

Death from enrollment until calving

km <- survfit(Surv(deadtar, dead) ~ scour.survival, data = calf.survival)
dead.scour <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.2,
       xmin = 0, xmax = 800,time.break=200,
       title.label = "Diarrhoea at all ages",
       x.label = "Days post enrollment",
       y.label = "Proportion of subjects that have died")

km <- survfit(Surv(deadtar, dead) ~ scour.survival, data = calf.survival %>% filter(agecat=="w1"))
dead.scour.w1 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.2,
       xmin = 0, xmax = 800, time.break=200,
       title.label = "Diarrhoea during first week of life",
       x.label = "Days post enrollment",
       y.label = "Proportion of subjects that have died")

km <- survfit(Surv(deadtar, dead) ~ scour.survival, data = calf.survival %>% filter(agecat=="w2"))
dead.scour.w2 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.2,xmin = 0, xmax = 800, time.break=200,
       title.label = "Diarrhoea during second week of life",
       x.label = "Days post enrollment",
       y.label = "Proportion of subjects that have died")

km <- survfit(Surv(deadtar, dead) ~ scour.survival, data = calf.survival %>% filter(agecat=="w3"))
dead.scour.w3 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.2,xmin = 0, xmax = 800, time.break=200,
       title.label = "Diarrhoea during third week of life",
       x.label = "Days post enrollment",
       y.label = "Proportion of subjects that have died")

km <- survfit(Surv(deadtar, dead) ~ scour.survival, data = calf.survival %>% filter(agecat=="w4+"))
dead.scour.w4 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.2,xmin = 0, xmax = 800, time.break=200,
       title.label = "Diarrhoea between 28d and weaning",
       x.label = "Days post enrollment",
       y.label = "Proportion of subjects that have died")

dead.scour

splots <- list(dead.scour.w1,dead.scour.w2,dead.scour.w3,dead.scour.w4)
arrange_ggsurvplots(splots,ncol=2)

arrange_ggsurvplots(splots,ncol=2,nrow=2)

km <- survfit(Surv(deadtar, dead) ~ scour.survival, data = calf.survival)


tab <- summary(km, times = c(20,100,700))
tibble(
  time = tab$time,
  survival = tab$surv,
  death = round((1 - survival) * 100,1)
)
plot_out <- ggsurvplot(km, data = calf.survival,
           pval.coord = c(0.03, 0.20),
           fun = "event",xlim=c(0,800),ylim=c(0,0.2),
           legend.labs = c("No Diarrhea","Diarrhea"),legend = c(0.2,0.7),legend.title = "",
           title = "", 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 since enrollment",
           ylab="Cumulative incidence of death",
           palette = "Set1",break.time.by = 100)

plot_out

ggsave("output.png",width = 5, height = 5)
dead.scour

splots <- list(dead.scour.w1,dead.scour.w2,dead.scour.w3,dead.scour.w4)
arrange_ggsurvplots(splots,ncol=2,nrow=2)

tiff("KM.tiff", units="in", width=13, height=6, res=100)
dev.off()
## quartz_off_screen 
##                 2

Death enrollment - 20 days

km <- survfit(Surv(dead20tar, dead20) ~ scour.survival, data = calf.survival)
dead20.scour <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.05,
       xmin = 0, xmax = 20,time.break=10,
       title.label = "Diarrhoea at all ages",
       x.label = "Days post enrollment",
       y.label = "Proportion of subjects that have died within 20 days")

km <- survfit(Surv(dead20tar, dead20) ~ scour.survival, data = calf.survival %>% filter(agecat=="w1"))
dead20.scour.w1 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.05,
       xmin = 0, xmax = 20,time.break=10,
       title.label = "Diarrhoea during first week of life",
       x.label = "Days post enrollment",
       y.label = "Proportion of subjects that have died within 20 days")

km <- survfit(Surv(dead20tar, dead20) ~ scour.survival, data = calf.survival %>% filter(agecat=="w2"))
dead20.scour.w2 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.05,
       xmin = 0, xmax = 20,time.break=10,
       title.label = "Diarrhoea during second week of life",
       x.label = "Days post enrollment",
       y.label = "Proportion of subjects that have died within 20 days")

km <- survfit(Surv(dead20tar, dead20) ~ scour.survival, data = calf.survival %>% filter(agecat=="w3"))
dead20.scour.w3 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.05,
       xmin = 0, xmax = 20,time.break=10,
       title.label = "Diarrhoea during third week of life",
       x.label = "Days post enrollment",
       y.label = "Proportion of subjects that have died within 20 days")

km <- survfit(Surv(dead20tar, dead20) ~ scour.survival, data = calf.survival %>% filter(agecat=="w4+"))
dead20.scour.w4 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.05,
       xmin = 0, xmax = 20,time.break=10,
       title.label = "Diarrhoea between 28d and weaning",
       x.label = "Days post enrollment",
       y.label = "Proportion of subjects that have died within 20 days")

splots <- list(dead20.scour.w1,dead20.scour.w2,dead20.scour.w3,dead20.scour.w4)
dead20.scour

arrange_ggsurvplots(splots,ncol=2)

Death 21 days - 100

km <- survfit(Surv(dead21_to_100tar, dead21_to_100) ~ scour.survival, data = calf.survival)
dead21_to_100.scour <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.04,
       xmin = 0, xmax = 80,time.break = 10,
       title.label = "Diarrhoea at all ages",
       x.label = "Days post 21d",
       y.label = "Proportion of subjects that have died between 21-100d")

km <- survfit(Surv(dead21_to_100tar, dead21_to_100) ~ scour.survival, data = calf.survival %>% filter(agecat=="w1"))
dead21_to_100.scour.w1 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.04,
       xmin = 0, xmax = 80,time.break = 10,
       title.label = "Diarrhoea during first week of life",
       x.label = "Days post 21d",
       y.label = "Proportion of subjects that have died between 21-100d")

km <- survfit(Surv(dead21_to_100tar, dead21_to_100) ~ scour.survival, data = calf.survival %>% filter(agecat=="w2"))
dead21_to_100.scour.w2 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.04,
       xmin = 0, xmax = 80,time.break = 10,
       title.label = "Diarrhoea during second week of life",
       x.label = "Days post 21d",
       y.label = "Proportion of subjects that have died between 21-100d")

km <- survfit(Surv(dead21_to_100tar, dead21_to_100) ~ scour.survival, data = calf.survival %>% filter(agecat=="w3"))
dead21_to_100.scour.w3 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.04,
       xmin = 0, xmax = 80,time.break = 10,
       title.label = "Diarrhoea during third week of life",
       x.label = "Days post 21d",
       y.label = "Proportion of subjects that have died between 21-100d")

km <- survfit(Surv(dead21_to_100tar, dead21_to_100) ~ scour.survival, data = calf.survival %>% filter(agecat=="w4+"))
dead21_to_100.scour.w4 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.04,
       xmin = 0, xmax = 80,time.break = 10,
       title.label = "Diarrhoea between 28d and weaning",
       x.label = "Days post 21d",
       y.label = "Proportion of subjects that have died between 21-100d")

splots <- list(dead21_to_100.scour.w1,dead21_to_100.scour.w2,dead21_to_100.scour.w3,dead21_to_100.scour.w4)
dead21_to_100.scour

arrange_ggsurvplots(splots,ncol=2)

Death 101 days - calving

km <- survfit(Surv(dead101_to_c1tar, dead101_to_c1) ~ scour.survival, data = calf.survival)
dead101_to_c1.scour <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.1,
       xmin = 0, xmax = 700,time.break = 200,
       title.label = "Diarrhoea at all ages",
       x.label = "Days post 101d",
       y.label = "Proportion of subjects that died between 101d and calving")

km <- survfit(Surv(dead101_to_c1tar, dead101_to_c1) ~ scour.survival, data = calf.survival %>% filter(agecat=="w1"))
dead101_to_c1.scour.w1 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.1,
       xmin = 0, xmax = 700,time.break = 200,
       title.label = "Diarrhoea during first week of life",
       x.label = "Days post 101d",
       y.label = "Proportion of subjects that died between 101d and calving")

km <- survfit(Surv(dead101_to_c1tar, dead101_to_c1) ~ scour.survival, data = calf.survival %>% filter(agecat=="w2"))
dead101_to_c1.scour.w2 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.1,
       xmin = 0, xmax = 700,time.break = 200,
       title.label = "Diarrhoea during second week of life",
       x.label = "Days post 101d",
       y.label = "Proportion of subjects that died between 101d and calving")

km <- survfit(Surv(dead101_to_c1tar, dead101_to_c1) ~ scour.survival, data = calf.survival %>% filter(agecat=="w3"))
dead101_to_c1.scour.w3 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.1,
       xmin = 0, xmax = 700,time.break = 200,
       title.label = "Diarrhoea during third week of life",
       x.label = "Days post 101d",
       y.label = "Proportion of subjects that died between 101d and calving")

km <- survfit(Surv(dead101_to_c1tar, dead101_to_c1) ~ scour.survival, data = calf.survival %>% filter(agecat=="w4+"))
dead101_to_c1.scour.w4 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.1,
       xmin = 0, xmax = 700,time.break = 200,
       title.label = "Diarrhoea between 28d and weaning",
       x.label = "Days post 101d",
       y.label = "Proportion of subjects that died between 101d and calving")

splots <- list(dead101_to_c1.scour.w1,dead101_to_c1.scour.w2,dead101_to_c1.scour.w3,dead101_to_c1.scour.w4)
dead101_to_c1.scour

arrange_ggsurvplots(splots,ncol=2)

arrange_ggsurvplots(splots,ncol=2, nrow=2)

c1

km <- survfit(Surv(c1tar, c1) ~ scour.survival, data = calf.survival)
c1.scour <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 1.0,
       xmin = 600, xmax = 1000,time.break= 100,
       title.label = "Diarrhoea at all ages",
       x.label = "Days post enrollment",
       y.label = "Proportion of subjects that have calved")

km <- survfit(Surv(c1tar, c1) ~ scour.survival, data = calf.survival %>% filter(agecat=="w1"))
c1.scour.w1 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 1.0,
       xmin = 600, xmax = 1000,time.break= 100,
       title.label = "Diarrhoea during first week of life",
       x.label = "Days post enrollment",
       y.label = "Proportion of subjects that have calved")

km <- survfit(Surv(c1tar, c1) ~ scour.survival, data = calf.survival %>% filter(agecat=="w2"))
c1.scour.w2 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 1.0,
       xmin = 600, xmax = 1000,time.break= 100,
       title.label = "Diarrhoea during second week of life",
       x.label = "Days post enrollment",
       y.label = "Proportion of subjects that have calved")

km <- survfit(Surv(c1tar, c1) ~ scour.survival, data = calf.survival %>% filter(agecat=="w3"))
c1.scour.w3 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 1.0,
       xmin = 600, xmax = 1000,time.break= 100,
       title.label = "Diarrhoea during third week of life",
       x.label = "Days post enrollment",
       y.label = "Proportion of subjects that have calved")

km <- survfit(Surv(c1tar, c1) ~ scour.survival, data = calf.survival %>% filter(agecat=="w4+"))
c1.scour.w4 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 1.0,
       xmin = 600, xmax = 1000,time.break= 100,
       title.label = "Diarrhoea between 28d and calving",
       x.label = "Days post enrollment",
       y.label = "Proportion of subjects that have calved")

splots <- list(c1.scour.w1,c1.scour.w2,c1.scour.w3,c1.scour.w4)
c1.scour

arrange_ggsurvplots(splots,ncol=2)

Median age to c1

Example interpretation of the first output: The median days to calving was 4 days later for calves with scours (691d) than calves without scours (687d).

Calculate median days to first calving

km <- survfit(Surv(c1tar, c1) ~ scour.survival, data = calf.survival)
km
## Call: survfit(formula = Surv(c1tar, c1) ~ scour.survival, data = calf.survival)
## 
##                     n events median 0.95LCL 0.95UCL
## scour.survival=0 7330   4470    687     686     689
## scour.survival=1 3643   2576    691     689     693
km.w1 <- survfit(Surv(c1tar, c1) ~ scour.survival, data = calf.survival %>% filter(agecat=="w1"))
km.w1
## Call: survfit(formula = Surv(c1tar, c1) ~ scour.survival, data = calf.survival %>% 
##     filter(agecat == "w1"))
## 
##                     n events median 0.95LCL 0.95UCL
## scour.survival=0 2118   1112    694     692     698
## scour.survival=1 1059    770    694     690     702
km.w2 <- survfit(Surv(c1tar, c1) ~ scour.survival, data = calf.survival %>% filter(agecat=="w2"))
km.w2
## Call: survfit(formula = Surv(c1tar, c1) ~ scour.survival, data = calf.survival %>% 
##     filter(agecat == "w2"))
## 
##                     n events median 0.95LCL 0.95UCL
## scour.survival=0 3933   2380    688     686     692
## scour.survival=1 1955   1317    693     689     697
km.w3 <- survfit(Surv(c1tar, c1) ~ scour.survival, data = calf.survival %>% filter(agecat=="w3"))
km.w3
## Call: survfit(formula = Surv(c1tar, c1) ~ scour.survival, data = calf.survival %>% 
##     filter(agecat == "w3"))
## 
##                    n events median 0.95LCL 0.95UCL
## scour.survival=0 795    604    683     678     689
## scour.survival=1 389    292    687     680     692
km.w4 <- survfit(Surv(c1tar, c1) ~ scour.survival, data = calf.survival %>% filter(agecat=="w4+"))
km.w4
## Call: survfit(formula = Surv(c1tar, c1) ~ scour.survival, data = calf.survival %>% 
##     filter(agecat == "w4+"))
## 
##                    n events median 0.95LCL 0.95UCL
## scour.survival=0 484    374    666     663     673
## scour.survival=1 240    197    667     661     673

Sold

km <- survfit(Surv(soldtar, sold) ~ scour.survival, data = calf.survival)
sold.scour <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.4,
       xmin = 0, xmax = 1000,time.break= 200,
       title.label = "Diarrhoea at all ages",
       x.label = "Days post enrollment",
       y.label = "Proportion of subjects that were sold")

km <- survfit(Surv(soldtar, sold) ~ scour.survival, data = calf.survival %>% filter(agecat=="w1"))
sold.scour.w1 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.4,
       xmin = 0, xmax = 1000,time.break= 200,
       title.label = "Diarrhoea during first week of life",
       x.label = "Days post enrollment",
       y.label = "Proportion of subjects that were sold")

km <- survfit(Surv(soldtar, sold) ~ scour.survival, data = calf.survival %>% filter(agecat=="w2"))
sold.scour.w2 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.4,
       xmin = 0, xmax = 1000,time.break= 200,
       title.label = "Diarrhoea during second week of life",
       x.label = "Days post enrollment",
       y.label = "Proportion of subjects that were sold")

km <- survfit(Surv(soldtar, sold) ~ scour.survival, data = calf.survival %>% filter(agecat=="w3"))
sold.scour.w3 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.4,
       xmin = 0, xmax = 1000,time.break= 200,
       title.label = "Diarrhoea during third week of life",
       x.label = "Days post enrollment",
       y.label = "Proportion of subjects that were sold")

km <- survfit(Surv(soldtar, sold) ~ scour.survival, data = calf.survival %>% filter(agecat=="w4+"))
sold.scour.w4 <- KMplot(dataframe = calf.survival,survival.function=km,
       ymax = 0.4,
       xmin = 0, xmax = 1000,time.break= 200,
       title.label = "Diarrhoea between 28d and calving",
       x.label = "Days post enrollment",
       y.label = "Proportion of subjects that were sold")

splots <- list(sold.scour.w1,sold.scour.w2,sold.scour.w3,sold.scour.w4)
sold.scour

arrange_ggsurvplots(splots,ncol=2)

KM plots using the calf.lact dataset

Eg. Removal (culled/died) from 1-300 DIM

km <- survfit(Surv(rem300tar, rem300) ~ scour, data = calf.lact)
rem300.scour <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 0.3,
       xmin = 0, xmax = 300,time.break=50,
       title.label = "Scours at all ages",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that were removed from the herd by 300 DIM")

km <- survfit(Surv(rem300tar, rem300) ~ scour_w1, data = calf.lact)
rem300.scour.w1 <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 0.3,
       xmin = 0, xmax = 300,time.break=50,
       title.label = "Diarrhoea during first week of life",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that were removed from the herd by 300 DIM")

km <- survfit(Surv(rem300tar, rem300) ~ scour_w2, data = calf.lact)
rem300.scour.w2 <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 0.3,xmin = 0, xmax = 300,time.break=50,
       title.label = "Diarrhoea during second week of life",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that were removed from the herd by 300 DIM")

km <- survfit(Surv(rem300tar, rem300) ~ scour_w3, data = calf.lact)
rem300.scour.w3 <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 0.3,xmin = 0, xmax = 300,time.break=50,
       title.label = "Diarrhoea during third week of life",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that were removed from the herd by 300 DIM")

km <- survfit(Surv(rem300tar, rem300) ~ scour_w4, data = calf.lact)
rem300.scour.w4 <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 0.3,xmin = 0, xmax = 300,time.break=50,
       title.label = "Diarrhoea between 28d and weaning",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that were removed from the herd by 300 DIM")

splots <- list(rem300.scour.w1,rem300.scour.w2,rem300.scour.w3,rem300.scour.w4)
rem300.scour

arrange_ggsurvplots(splots,ncol=2)

Mastitis within first 120 DIM

km <- survfit(Surv(mast120tar, mast120) ~ scour, data = calf.lact)
mast120.scour <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 0.03,
       xmin = 0, xmax = 120,time.break=30,
       title.label = "Scours at all ages",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that contracted mastitis within 120 DIM")

km <- survfit(Surv(mast120tar, mast120) ~ scour_w1, data = calf.lact)
mast120.scour.w1 <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 0.03,
       xmin = 0, xmax = 120,time.break=30,
       title.label = "Scours during first week of life",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that contracted mastitis within 120 DIM")

km <- survfit(Surv(mast120tar, mast120) ~ scour_w2, data = calf.lact)
mast120.scour.w2 <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 0.03,
       xmin = 0, xmax = 120,time.break=30,
       title.label = "Scours during second week of life",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that contracted mastitis within 120 DIM")

km <- survfit(Surv(mast120tar, mast120) ~ scour_w3, data = calf.lact)
mast120.scour.w3 <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 0.03,
       xmin = 0, xmax = 120,time.break=30,
       title.label = "Scours during third week of life",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that contracted mastitis within 120 DIM")

km <- survfit(Surv(mast120tar, mast120) ~ scour_w4, data = calf.lact)
mast120.scour.w4 <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 0.03,
       xmin = 0, xmax = 120,time.break=30,
       title.label = "Scours between 28d and weaning",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that contracted mastitis within 120 DIM")

splots <- list(mast120.scour.w1,mast120.scour.w2,mast120.scour.w3,mast120.scour.w4)
mast120.scour

arrange_ggsurvplots(splots,ncol=2)

Metritis within first 10 DIM

km <- survfit(Surv(met10tar, met10) ~ scour, data = calf.lact)
met10.scour <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 0.06,
       xmin = 0, xmax = 10,time.break=2,
       title.label = "Scours at all ages",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that contracted metritis within 10 DIM")

km <- survfit(Surv(met10tar, met10) ~ scour_w1, data = calf.lact)
met10.scour.w1 <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 0.06,
       xmin = 0, xmax = 10,time.break=2,
       title.label = "Scours during the first week of life",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that contracted metritis within 10 DIM")

km <- survfit(Surv(met10tar, met10) ~ scour_w2, data = calf.lact)
met10.scour.w2 <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 0.06,
       xmin = 0, xmax = 10,time.break=2,
       title.label = "Scours during the second week of life",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that contracted metritis within 10 DIM")

km <- survfit(Surv(met10tar, met10) ~ scour_w3, data = calf.lact)
met10.scour.w3 <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 0.06,
       xmin = 0, xmax = 10,time.break=2,
       title.label = "Scours during the third week of life",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that contracted metritis within 10 DIM")

km <- survfit(Surv(met10tar, met10) ~ scour_w4, data = calf.lact)
met10.scour.w4 <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 0.06,
       xmin = 0, xmax = 10,time.break=2,
       title.label = "Scours between 28d and calving",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that contracted metritis within 10 DIM")

met10.scour

splots <- list(met10.scour.w1,met10.scour.w2,met10.scour.w3,met10.scour.w4)
arrange_ggsurvplots(splots,ncol=2)

Retained foetal membranes during first 10 DIM

km <- survfit(Surv(rfm10tar, rfm10) ~ scour, data = calf.lact)
rfm10.scour <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 0.04,
       xmin = 0, xmax = 10,time.break=2,
       title.label = "Scours at all ages",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that had retained foetal membranes within 10 DIM")

km <- survfit(Surv(rfm10tar, rfm10) ~ scour_w1, data = calf.lact)
rfm10.scour.w1 <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 0.04,
       xmin = 0, xmax = 10,time.break=2,
       title.label = "Scours during first week of life",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that had retained foetal membranes within 10 DIM")

km <- survfit(Surv(rfm10tar, rfm10) ~ scour_w2, data = calf.lact)
rfm10.scour.w2 <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 0.04,
       xmin = 0, xmax = 10,time.break=2,
       title.label = "Scours during second week of life",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that had retained foetal membranes within 10 DIM")

km <- survfit(Surv(rfm10tar, rfm10) ~ scour_w3, data = calf.lact)
rfm10.scour.w3 <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 0.04,
       xmin = 0, xmax = 10,time.break=2,
       title.label = "Scours during third week of life",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that had retained foetal membranes within 10 DIM")

km <- survfit(Surv(rfm10tar, rfm10) ~ scour_w4, data = calf.lact)
rfm10.scour.w4 <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 0.04,
       xmin = 0, xmax = 10,time.break=2,
       title.label = "Scours between 28d and calving",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that had retained foetal membranes within 10 DIM")

rfm10.scour

splots <- list(rfm10.scour.w1,rfm10.scour.w2,rfm10.scour.w3,rfm10.scour.w4)
arrange_ggsurvplots(splots,ncol=2)

Pregnancy within 200 DIM

km <- survfit(Surv(preg200tar, preg200) ~ scour, data = calf.lact)
preg200.scour <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 1.0,
       xmin = 50, xmax = 200,time.break=50,
       title.label = "Scours at all ages",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that are pregnant")

km <- survfit(Surv(preg200tar, preg200) ~ scour_w1, data = calf.lact)
preg200.scour.w1 <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 1.0,
       xmin = 50, xmax = 200,time.break=50,
       title.label = "Scours during first week of life",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that are pregnant")

km <- survfit(Surv(preg200tar, preg200) ~ scour_w2, data = calf.lact)
preg200.scour.w2 <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 1.0,
       xmin = 50, xmax = 200,time.break=50,
       title.label = "Scours during second week of life",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that are pregnant")

km <- survfit(Surv(preg200tar, preg200) ~ scour_w3, data = calf.lact)
preg200.scour.w3 <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 1.0,
       xmin = 50, xmax = 200,time.break=50,
       title.label = "Scours during third week of life",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that are pregnant")

km <- survfit(Surv(preg200tar, preg200) ~ scour_w4, data = calf.lact)
preg200.scour.w4 <- KMplot(dataframe = calf.lact,survival.function=km,
       ymax = 1.0,
       xmin = 50, xmax = 200,time.break=50,
       title.label = "Scours between 28d and weaning",
       x.label = "Days after first calving",
       y.label = "Proportion of subjects that are pregnant")

splots <- list(preg200.scour.w1,preg200.scour.w2,preg200.scour.w3,preg200.scour.w4)
preg200.scour

arrange_ggsurvplots(splots,ncol=2)

Cox regression

cox.all.extract <- function(model,group){
cph <- model %>% tidy(exp=T,conf.int=T) %>% select(term,estimate,conf.low,conf.high,p.value) %>% slice(1)
cph$predictor <- cph$term
mm0 <- cph
mm0$hazard.ratio <- paste0(sprintf("%.2f",round(mm0$estimate,2))," (",sprintf("%.2f",round(mm0$conf.low,2)),
                   ", ",sprintf("%.2f",round(mm0$conf.high,2)),")")
mm0$p.value.hr <- mm0$p.value
mm0$outcome <- as.character(model[["call"]][["formula"]][[2]][[3]])
mm0$week <- group
mm0 %>% select(outcome,predictor,week,hazard.ratio)
}


cox.weeks.extract <- function(model){
  cph <- model
em <- emmeans(cph,pairwise ~ scour.survival | agecat,adjust = "none",type = "response",reverse = TRUE) %>% confint
tab1 <- data.frame(
  contrast = em[["contrasts"]][["agecat"]],
  hr = paste0(round(1/em[["contrasts"]][["ratio"]],2)," (",
              round(1/em[["contrasts"]][["asymp.UCL"]],2), " to ",
              round(1/em[["contrasts"]][["asymp.LCL"]],2),")")
)
tab1$outcome <- as.character(model[["call"]][["formula"]][[2]][[3]])
tab1$predictor <- as.character(model[["terms"]][[3]][[2]][[2]])
tab1$week <- tab1$contrast
tab1$hazard.ratio <- tab1$hr
tab1 <- tab1 %>% select(outcome,predictor,week,hazard.ratio)
tab1
}

Cox regression in the calf.survival dataset

Eg. Death between enrollment and calving

Step 1: Run the model. You need to enter the event (e.g. ‘dead), the time at risk (’deadtar’), the predictor/exposure (scour.survival), covariates (birth.season) and the dataframe/database (calf.survival). After running the model, look at the summary.

cox <- coxph(Surv(deadtar,dead) ~ scour.survival + strata(birth.season),data = calf.survival)
summary(cox)
## Call:
## coxph(formula = Surv(deadtar, dead) ~ scour.survival + strata(birth.season), 
##     data = calf.survival)
## 
##   n= 10973, number of events= 734 
## 
##                   coef exp(coef) se(coef)     z Pr(>|z|)    
## scour.survival 0.39504   1.48445  0.07428 5.318 1.05e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                exp(coef) exp(-coef) lower .95 upper .95
## scour.survival     1.484     0.6737     1.283     1.717
## 
## Concordance= 0.559  (se = 0.01 )
## Likelihood ratio test= 27.79  on 1 df,   p=1e-07
## Wald test            = 28.29  on 1 df,   p=1e-07
## Score (logrank) test = 28.65  on 1 df,   p=9e-08

p-value << 0.05; significant

Step 2: Check the assumption of proportional hazards. The p-value should be greater than 0.05. if it is less than 0.05, then it is likely that the effect of your predictor (e.g. scours) impacts risk of outcome (e.g. death) differently over time. In this case, the p-value for scour.survival is very small, which suggests that pre-weaning scours (any week) doesn’t have a uniform impact on death risk between the enrollment and calving. This is because, most of the death due to scours happens in the short-term and less in the long-term.

cox.zph(cox)
##                chisq df       p
## scour.survival  28.9  1 7.5e-08
## GLOBAL          28.9  1 7.5e-08
cox.zph(cox) %>% ggcoxzph

p << 0.05; assumption violated

Step 3: Extract the hazard ratio from the model

dead.scour.cox <- cox.all.extract(model = cox, group = "All weeks")
dead.scour.cox

Step 4: Get hazard ratios for individual weeks by creating a model with an interaction term (don’t worry too much about what this means). Note that you can only do this for the calf.survival dataset. I.e. don’t use the cox.weeks.extract function on the other dataset.

cox.weeks <- coxph(Surv(deadtar,dead) ~ scour.survival*agecat + birth.season,data = calf.survival)
dead.scour.cox.weeks <- cox.weeks.extract(model=cox.weeks)
dead.scour.cox.weeks

Step 5: Check assumptions for the new model. Again, this model isn’t great because the p-values are all < 0.05.

cox.zph(cox.weeks)
##                       chisq df       p
## scour.survival        29.28  1 6.3e-08
## agecat                11.97  3  0.0075
## birth.season           2.36  3  0.5003
## scour.survival:agecat 13.02  3  0.0046
## GLOBAL                52.72 10 8.4e-08
cox.zph(cox.weeks) %>% ggcoxzph

dead.scour.hr <- rbind(dead.scour.cox,
      dead.scour.cox.weeks)

dead.scour.hr

p << 0.05; assumption violated

Death enrollment - 20 days

step 1 - initial model

cox <- coxph(Surv(dead20tar,dead20) ~ scour.survival + strata(birth.season),data = calf.survival)
summary(cox)
## Call:
## coxph(formula = Surv(dead20tar, dead20) ~ scour.survival + strata(birth.season), 
##     data = calf.survival)
## 
##   n= 10973, number of events= 197 
## 
##                  coef exp(coef) se(coef)     z Pr(>|z|)    
## scour.survival 0.9077    2.4786   0.1438 6.312 2.75e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                exp(coef) exp(-coef) lower .95 upper .95
## scour.survival     2.479     0.4034      1.87     3.286
## 
## Concordance= 0.616  (se = 0.018 )
## Likelihood ratio test= 40.02  on 1 df,   p=3e-10
## Wald test            = 39.84  on 1 df,   p=3e-10
## Score (logrank) test = 42.64  on 1 df,   p=7e-11

p-value << 0.05; significant

step 2 - check assumption of proportional hazards in initial model

cox.zph(cox)
##                chisq df   p
## scour.survival 0.452  1 0.5
## GLOBAL         0.452  1 0.5
cox.zph(cox) %>% ggcoxzph

p = 0.6 > 0.05; assumption met

step 3 - hazard ratios

dead20.scour.cox <- cox.all.extract(model = cox, group = "All weeks")
dead20.scour.cox

step 4 - hazard ratios for individual weeks

cox.weeks <- coxph(Surv(dead20tar,dead20) ~ scour.survival*agecat + strata(birth.season),data = calf.survival)
dead20.scour.cox.weeks <- cox.weeks.extract(model=cox.weeks)
dead20.scour.cox.weeks

step 5 - check assumption of proportional hazards in individual week model

cox.zph(cox.weeks)
##                       chisq df    p
## scour.survival        0.293  1 0.59
## agecat                0.949  3 0.81
## scour.survival:agecat 2.081  3 0.56
## GLOBAL                7.232  7 0.41
cox.zph(cox.weeks) %>% ggcoxzph

p = 0.69 > 0.05; assumption met

step 6 - combining all weeks hazard ratio

dead20.scour.hr <- rbind(dead20.scour.cox,
      dead20.scour.cox.weeks)

dead20.scour.hr

Death 21 days - 100

step 1 - initial model

cox <- coxph(Surv(dead21_to_100tar,dead21_to_100) ~ scour.survival + strata(birth.season),data = calf.survival)
summary(cox)
## Call:
## coxph(formula = Surv(dead21_to_100tar, dead21_to_100) ~ scour.survival + 
##     strata(birth.season), data = calf.survival)
## 
##   n= 9690, number of events= 175 
##    (1283 observations deleted due to missingness)
## 
##                  coef exp(coef) se(coef)     z Pr(>|z|)    
## scour.survival 0.6386    1.8939   0.1516 4.213 2.52e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                exp(coef) exp(-coef) lower .95 upper .95
## scour.survival     1.894      0.528     1.407     2.549
## 
## Concordance= 0.58  (se = 0.019 )
## Likelihood ratio test= 17.65  on 1 df,   p=3e-05
## Wald test            = 17.75  on 1 df,   p=3e-05
## Score (logrank) test = 18.36  on 1 df,   p=2e-05

p-value < 0.05; significant

step 2 - check assumption of proportional hazards in initial model

cox.zph(cox)
##                chisq df    p
## scour.survival 0.227  1 0.63
## GLOBAL         0.227  1 0.63
cox.zph(cox) %>% ggcoxzph

p = 0.63 > 0.05; assumption met

step 3 - hazard ratios

dead21_to_100.scour.cox <- cox.all.extract(model = cox, group = "All weeks")
dead21_to_100.scour.cox

step 4 - hazard ratios for individual weeks

cox.weeks <- coxph(Surv(dead21_to_100tar,dead21_to_100) ~ scour.survival*agecat + birth.season,data = calf.survival)
dead21_to_100.scour.cox.weeks <- cox.weeks.extract(model=cox.weeks)
dead21_to_100.scour.cox.weeks

step 5 -check assumption of proportional hazards in individual week model

cox.zph(cox.weeks)
##                       chisq df    p
## scour.survival        0.318  1 0.57
## agecat                1.835  3 0.61
## birth.season          3.961  3 0.27
## scour.survival:agecat 3.475  3 0.32
## GLOBAL                8.682 10 0.56
cox.zph(cox.weeks) %>% ggcoxzph

p = 0.57 > 0.05; assumption met

step 6 - combining all weeks hazard ratio

dead21_to_100.scour.hr <- rbind(dead21_to_100.scour.cox,
      dead21_to_100.scour.cox.weeks)

dead21_to_100.scour.hr

Death 101 days - calving

step 1 - initial model

cox <- coxph(Surv(dead101_to_c1tar,dead101_to_c1) ~ scour.survival + strata(birth.season),data = calf.survival)
summary(cox)
## Call:
## coxph(formula = Surv(dead101_to_c1tar, dead101_to_c1) ~ scour.survival + 
##     strata(birth.season), data = calf.survival)
## 
##   n= 9421, number of events= 362 
##    (1552 observations deleted due to missingness)
## 
##                    coef exp(coef) se(coef)      z Pr(>|z|)
## scour.survival -0.01541   0.98470  0.10960 -0.141    0.888
## 
##                exp(coef) exp(-coef) lower .95 upper .95
## scour.survival    0.9847      1.016    0.7944     1.221
## 
## Concordance= 0.495  (se = 0.013 )
## Likelihood ratio test= 0.02  on 1 df,   p=0.9
## Wald test            = 0.02  on 1 df,   p=0.9
## Score (logrank) test = 0.02  on 1 df,   p=0.9

p-value = 0.891 > 0.05; insignificant

step 2 - check assumption of proportional hazards in initial model

cox.zph(cox)
##                chisq df    p
## scour.survival  2.59  1 0.11
## GLOBAL          2.59  1 0.11
cox.zph(cox) %>% ggcoxzph

p = 0.11 > 0.05; assumption met

step 3 - hazard ratios

dead101_to_c1.scour.cox <- cox.all.extract(model = cox, group = "All weeks")
dead101_to_c1.scour.cox

step 4 - hazard ratios for individual weeks

cox.weeks <- coxph(Surv(dead101_to_c1tar,dead101_to_c1) ~ scour.survival*agecat + birth.season,data = calf.survival)
dead101_to_c1.scour.cox.weeks <- cox.weeks.extract(model=cox.weeks)
dead101_to_c1.scour.cox.weeks

step 5 -check assumption of proportional hazards in individual week model

cox.zph(cox.weeks)
##                       chisq df     p
## scour.survival         2.97  1 0.085
## agecat                 2.13  3 0.546
## birth.season           3.42  3 0.332
## scour.survival:agecat  2.96  3 0.398
## GLOBAL                12.32 10 0.264
cox.zph(cox.weeks) %>% ggcoxzph

p = 0.085 > 0.05; assumption met

step 6 - combining all weeks hazard ratio

dead101_to_c1.scour.hr <- rbind(dead101_to_c1.scour.cox,
      dead101_to_c1.scour.cox.weeks)

dead101_to_c1.scour.hr

c1

step 1 - initial model

cox <- coxph(Surv(c1tar,c1) ~ scour.survival + strata(birth.season),data = calf.survival)
summary(cox)
## Call:
## coxph(formula = Surv(c1tar, c1) ~ scour.survival + strata(birth.season), 
##     data = calf.survival)
## 
##   n= 10973, number of events= 7046 
## 
##                    coef exp(coef) se(coef)      z Pr(>|z|)  
## scour.survival -0.04946   0.95175  0.02479 -1.995    0.046 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                exp(coef) exp(-coef) lower .95 upper .95
## scour.survival    0.9517      1.051    0.9066    0.9991
## 
## Concordance= 0.509  (se = 0.003 )
## Likelihood ratio test= 4  on 1 df,   p=0.05
## Wald test            = 3.98  on 1 df,   p=0.05
## Score (logrank) test = 3.98  on 1 df,   p=0.05

p-value = 0.0422 < 0.05; significant

step 2 - check assumption of proportional hazards in initial model

cox.zph(cox)
##                chisq df     p
## scour.survival  4.56  1 0.033
## GLOBAL          4.56  1 0.033
cox.zph(cox) %>% ggcoxzph

p = 0.028 < 0.05; assumption violated

step 3 - hazard ratios

c1.scour.cox <- cox.all.extract(model = cox, group = "All weeks")
c1.scour.cox

step 4 - hazard ratios for individual weeks

cox.weeks <- coxph(Surv(c1tar,c1) ~ scour.survival*agecat + birth.season,data = calf.survival)
c1.scour.cox.weeks <- cox.weeks.extract(model=cox.weeks)
c1.scour.cox.weeks

step 5 - check assumption of proportional hazards in individual week model

cox.zph(cox.weeks)
##                        chisq df      p
## scour.survival          4.81  1  0.028
## agecat                312.33  3 <2e-16
## birth.season           10.54  3  0.015
## scour.survival:agecat  77.50  3 <2e-16
## GLOBAL                331.29 10 <2e-16
cox.zph(cox.weeks) %>% ggcoxzph

p = 0.024 < 0.05; assumption violated

step 6- combining all weeks hazard ratio

c1.scour.hr <- rbind(c1.scour.cox,
      c1.scour.cox.weeks)

c1.scour.hr

Sold

step 1 - initial model

cox <- coxph(Surv(soldtar,sold) ~ scour.survival + strata(birth.season),data = calf.survival)
summary(cox)
## Call:
## coxph(formula = Surv(soldtar, sold) ~ scour.survival + strata(birth.season), 
##     data = calf.survival)
## 
##   n= 10973, number of events= 1199 
## 
##                    coef exp(coef) se(coef)      z Pr(>|z|)  
## scour.survival -0.12637   0.88129  0.06125 -2.063   0.0391 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                exp(coef) exp(-coef) lower .95 upper .95
## scour.survival    0.8813      1.135    0.7816    0.9937
## 
## Concordance= 0.513  (se = 0.007 )
## Likelihood ratio test= 4.31  on 1 df,   p=0.04
## Wald test            = 4.26  on 1 df,   p=0.04
## Score (logrank) test = 4.26  on 1 df,   p=0.04

p-value = 0.00978 < 0.05; significant

step 2 - check assumption of proportional hazards in initial model

cox.zph(cox)
##                chisq df   p
## scour.survival  0.15  1 0.7
## GLOBAL          0.15  1 0.7
cox.zph(cox) %>% ggcoxzph

p = 0.95 > 0.05; assumption met

step 3 - hazard ratios

sold.scour.cox <- cox.all.extract(model = cox, group = "All weeks")
sold.scour.cox

step 4 - hazard ratios for individual weeks

cox.weeks <- coxph(Surv(soldtar,sold) ~ scour.survival*agecat + birth.season,data = calf.survival)
sold.scour.cox.weeks <- cox.weeks.extract(model=cox.weeks)
sold.scour.cox.weeks

step 5 - check assumption of proportional hazards in individual week model

cox.zph(cox.weeks)
##                         chisq df       p
## scour.survival          0.738  1    0.39
## agecat                 21.935  3 6.7e-05
## birth.season          138.680  3 < 2e-16
## scour.survival:agecat   2.596  3    0.46
## GLOBAL                181.301 10 < 2e-16
cox.zph(cox.weeks) %>% ggcoxzph

p = 0.60 > 0.05; assumption met

step 6- combining all weeks hazard ratio

sold.scour.hr <- rbind(sold.scour.cox,
      sold.scour.cox.weeks)

sold.scour.hr

###Table combining calf.survival hazard ratios

calf.survival.hr.all <- rbind(dead.scour.hr,
                              dead20.scour.hr,
                              dead21_to_100.scour.hr,
                              dead101_to_c1.scour.hr,
                              c1.scour.hr,
                              sold.scour.hr)

calf.survival.hr.all
write.csv(calf.survival.hr.all,'calf.survival.hr.all.csv')

Interpretation on significant hazard ratios: 1. Preweaning scours increase chance of death by 1.47x as compared to without preweaning scours, when stratified into weeks, weeks 1-3 are significant, but week 4+ isnt

  1. Preweaning scours increase chance of death by 20 days of age by 2.46x as compared to without preweaning scours, when stratified into weeks, weeks 1-3 are significant but week 4 isn’t. However, week 3 seems to have a really high hazard ratio; considering the range of the confidence interval, this may be unreliable due to small sample size.

  2. Preweaning scours increase chance of death between 21-100d by 1.86x, when stratified into weeks, weeks 2 and 3 are significant but weeks 1 and 4 are not.

  3. Preweaning scours does not seem to have an effect on death between 101d and calving (HR = 0.99) i.e. does not seem to have a long term impact

  4. Considering the CI and HR, preweaning scours during week 2 seem to have the greatest impact short and medium term (<100d), although I am unsure if this is a right interpretation. Dead20_week2 HR 2.76 (1.82, 4.18) vs dead20_week1 HR 2.02 (1.33, 3.07). Although not significantly different, both the upper and lower range of the CI in week 2 is higher than week 1, with the overall value of the range being similar. Same can be said for dead21_to_100.

  5. Hazard ratios in week 3 have limited reliability due to CI having a large range; possibly due to small sample size and having insufficient power. Also unsure about this interpretation.

Cox regression in the calf.lact dataset

E.g. Removal (death/culling) during 1-300 DIM

Step 1: Run the model. You need to enter the event (e.g. ‘rem300’), the time at risk (‘rem300tar’), the predictor/exposure (‘scour’), covariates (‘calving.season’ - note that it should have “strata” in front of it) and the dataframe/database (‘calf.lact’). After running the model, look at the summary.

cox <- coxph(Surv(rem300tar,rem300) ~ scour + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(rem300tar, rem300) ~ scour + strata(calving.season), 
##     data = calf.lact)
## 
##   n= 7047, number of events= 1223 
## 
##          coef exp(coef) se(coef)     z Pr(>|z|)  
## scour 0.12035   1.12789  0.05864 2.052   0.0401 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## scour     1.128     0.8866     1.005     1.265
## 
## Concordance= 0.515  (se = 0.007 )
## Likelihood ratio test= 4.17  on 1 df,   p=0.04
## Wald test            = 4.21  on 1 df,   p=0.04
## Score (logrank) test = 4.22  on 1 df,   p=0.04

p-value = 0.0758 > 0.05; insignificant

Step 2: Check the assumption of proportional hazards.

cox.zph(cox)
##        chisq df    p
## scour  0.176  1 0.67
## GLOBAL 0.176  1 0.67
cox.zph(cox) %>% ggcoxzph

p = 0.85 > 0.05; assumption met

Step 3: Extract the hazard ratio from the model

rem300.scour.cox <- cox.all.extract(model = cox, group = "All weeks")
rem300.scour.cox

Steps 4+: Repeat the process for individual weeks. You can combine the code into a single chunk.

cox <- coxph(Surv(rem300tar,rem300) ~ scour_w1 + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(rem300tar, rem300) ~ scour_w1 + strata(calving.season), 
##     data = calf.lact)
## 
##   n= 7047, number of events= 1223 
## 
##            coef exp(coef) se(coef)     z Pr(>|z|)
## scour_w1 0.0510    1.0523   0.1069 0.477    0.633
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## scour_w1     1.052     0.9503    0.8534     1.298
## 
## Concordance= 0.502  (se = 0.004 )
## Likelihood ratio test= 0.22  on 1 df,   p=0.6
## Wald test            = 0.23  on 1 df,   p=0.6
## Score (logrank) test = 0.23  on 1 df,   p=0.6
cox.zph(cox)
##          chisq df    p
## scour_w1   0.2  1 0.65
## GLOBAL     0.2  1 0.65
cox.zph(cox) %>% ggcoxzph

rem300.scour.w1.cox <- cox.all.extract(model = cox, group = "w1")
rem300.scour.w1.cox
cox <- coxph(Surv(rem300tar,rem300) ~ scour_w2 + scour_w2_prev + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(rem300tar, rem300) ~ scour_w2 + scour_w2_prev + 
##     strata(calving.season), data = calf.lact)
## 
##   n= 7047, number of events= 1223 
## 
##                  coef exp(coef) se(coef)     z Pr(>|z|)  
## scour_w2      0.14258   1.15325  0.06760 2.109   0.0349 *
## scour_w2_prev 0.06856   1.07096  0.10729 0.639   0.5228  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## scour_w2          1.153     0.8671    1.0101     1.317
## scour_w2_prev     1.071     0.9337    0.8679     1.322
## 
## Concordance= 0.515  (se = 0.007 )
## Likelihood ratio test= 4.56  on 2 df,   p=0.1
## Wald test            = 4.67  on 2 df,   p=0.1
## Score (logrank) test = 4.68  on 2 df,   p=0.1
cox.zph(cox)
##               chisq df    p
## scour_w2       1.79  1 0.18
## scour_w2_prev  0.20  1 0.65
## GLOBAL         1.91  2 0.39
cox.zph(cox) %>% ggcoxzph

rem300.scour.w2.cox <- cox.all.extract(model = cox, group = "w2")
rem300.scour.w2.cox
cox <- coxph(Surv(rem300tar,rem300) ~ scour_w3 + scour_w3_prev + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(rem300tar, rem300) ~ scour_w3 + scour_w3_prev + 
##     strata(calving.season), data = calf.lact)
## 
##   n= 7047, number of events= 1223 
## 
##                   coef exp(coef) se(coef)      z Pr(>|z|)  
## scour_w3      -0.03679   0.96388  0.11800 -0.312   0.7552  
## scour_w3_prev  0.12755   1.13604  0.06229  2.048   0.0406 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## scour_w3         0.9639     1.0375    0.7649     1.215
## scour_w3_prev    1.1360     0.8803    1.0055     1.284
## 
## Concordance= 0.514  (se = 0.007 )
## Likelihood ratio test= 4.31  on 2 df,   p=0.1
## Wald test            = 4.38  on 2 df,   p=0.1
## Score (logrank) test = 4.39  on 2 df,   p=0.1
cox.zph(cox)
##               chisq df    p
## scour_w3      0.652  1 0.42
## scour_w3_prev 1.062  1 0.30
## GLOBAL        1.624  2 0.44
cox.zph(cox) %>% ggcoxzph

rem300.scour.w3.cox <- cox.all.extract(model = cox, group = "w3")
rem300.scour.w3.cox
cox <- coxph(Surv(rem300tar,rem300) ~ scour_w4 + scour_w4_prev + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(rem300tar, rem300) ~ scour_w4 + scour_w4_prev + 
##     strata(calving.season), data = calf.lact)
## 
##   n= 7047, number of events= 1223 
## 
##                  coef exp(coef) se(coef)     z Pr(>|z|)  
## scour_w4      0.08944   1.09356  0.13937 0.642   0.5210  
## scour_w4_prev 0.11479   1.12164  0.05974 1.922   0.0547 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## scour_w4          1.094     0.9144    0.8322     1.437
## scour_w4_prev     1.122     0.8916    0.9977     1.261
## 
## Concordance= 0.515  (se = 0.007 )
## Likelihood ratio test= 3.98  on 2 df,   p=0.1
## Wald test            = 4.03  on 2 df,   p=0.1
## Score (logrank) test = 4.04  on 2 df,   p=0.1
cox.zph(cox)
##                chisq df    p
## scour_w4      0.0863  1 0.77
## scour_w4_prev 0.1502  1 0.70
## GLOBAL        0.2300  2 0.89
cox.zph(cox) %>% ggcoxzph

rem300.scour.w4.cox <- cox.all.extract(model = cox, group = "w4+")
rem300.scour.w4.cox

Final step is to combine all the outcomes into a single table

rem300.cox <- rbind(rem300.scour.cox,
                    rem300.scour.w1.cox,
                    rem300.scour.w2.cox,
                    rem300.scour.w3.cox,
                    rem300.scour.w4.cox)
rem300.cox

Mastitis within first 120 DIM

step 1 - initial model

cox <- coxph(Surv(mast120tar,mast120) ~ scour + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(mast120tar, mast120) ~ scour + strata(calving.season), 
##     data = calf.lact)
## 
##   n= 7047, number of events= 180 
## 
##          coef exp(coef) se(coef)     z Pr(>|z|)
## scour 0.01015   1.01021  0.15474 0.066    0.948
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## scour      1.01     0.9899    0.7459     1.368
## 
## Concordance= 0.5  (se = 0.018 )
## Likelihood ratio test= 0  on 1 df,   p=0.9
## Wald test            = 0  on 1 df,   p=0.9
## Score (logrank) test = 0  on 1 df,   p=0.9

p-value = 0.983 > 0.05; insignificant

step 2 - check assumption of proportional hazards in initial model

cox.zph(cox)
##        chisq df    p
## scour   0.63  1 0.43
## GLOBAL  0.63  1 0.43
cox.zph(cox) %>% ggcoxzph

p = 0.41 > 0.05; assumption met

step 3 - hazard ratios

mast120.scour.cox <- cox.all.extract(model = cox, group = "All weeks")
mast120.scour.cox

step 4 - individual weeks

cox <- coxph(Surv(mast120tar,mast120) ~ scour_w1 + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(mast120tar, mast120) ~ scour_w1 + strata(calving.season), 
##     data = calf.lact)
## 
##   n= 7047, number of events= 180 
## 
##             coef exp(coef) se(coef)      z Pr(>|z|)
## scour_w1 -0.4047    0.6672   0.3421 -1.183    0.237
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## scour_w1    0.6672      1.499    0.3412     1.305
## 
## Concordance= 0.512  (se = 0.008 )
## Likelihood ratio test= 1.58  on 1 df,   p=0.2
## Wald test            = 1.4  on 1 df,   p=0.2
## Score (logrank) test = 1.42  on 1 df,   p=0.2
cox.zph(cox)
##          chisq df    p
## scour_w1  1.38  1 0.24
## GLOBAL    1.38  1 0.24
cox.zph(cox) %>% ggcoxzph

mast120.scour.w1.cox <- cox.all.extract(model = cox, group = "w1")
mast120.scour.w1.cox
cox <- coxph(Surv(mast120tar,mast120) ~ scour_w2 + scour_w2_prev + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(mast120tar, mast120) ~ scour_w2 + scour_w2_prev + 
##     strata(calving.season), data = calf.lact)
## 
##   n= 7047, number of events= 180 
## 
##                  coef exp(coef) se(coef)      z Pr(>|z|)
## scour_w2       0.1604    1.1740   0.1739  0.922    0.356
## scour_w2_prev -0.3847    0.6806   0.3429 -1.122    0.262
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## scour_w2         1.1740     0.8518    0.8349     1.651
## scour_w2_prev    0.6806     1.4692    0.3476     1.333
## 
## Concordance= 0.522  (se = 0.017 )
## Likelihood ratio test= 2.41  on 2 df,   p=0.3
## Wald test            = 2.26  on 2 df,   p=0.3
## Score (logrank) test = 2.28  on 2 df,   p=0.3
cox.zph(cox)
##               chisq df    p
## scour_w2      0.148  1 0.70
## scour_w2_prev 1.371  1 0.24
## GLOBAL        1.467  2 0.48
cox.zph(cox) %>% ggcoxzph

mast120.scour.w2.cox <- cox.all.extract(model = cox, group = "w2")
mast120.scour.w2.cox
cox <- coxph(Surv(mast120tar,mast120) ~ scour_w3 + scour_w3_prev + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(mast120tar, mast120) ~ scour_w3 + scour_w3_prev + 
##     strata(calving.season), data = calf.lact)
## 
##   n= 7047, number of events= 180 
## 
##                   coef exp(coef) se(coef)      z Pr(>|z|)
## scour_w3      -0.21331   0.80790  0.32598 -0.654    0.513
## scour_w3_prev  0.01559   1.01571  0.16579  0.094    0.925
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## scour_w3         0.8079     1.2378    0.4265     1.531
## scour_w3_prev    1.0157     0.9845    0.7339     1.406
## 
## Concordance= 0.509  (se = 0.018 )
## Likelihood ratio test= 0.47  on 2 df,   p=0.8
## Wald test            = 0.45  on 2 df,   p=0.8
## Score (logrank) test = 0.45  on 2 df,   p=0.8
cox.zph(cox)
##                chisq df    p
## scour_w3      0.2738  1 0.60
## scour_w3_prev 0.0413  1 0.84
## GLOBAL        0.3277  2 0.85
cox.zph(cox) %>% ggcoxzph

mast120.scour.w3.cox <- cox.all.extract(model = cox, group = "w3")
mast120.scour.w3.cox
cox <- coxph(Surv(mast120tar,mast120) ~ scour_w4 + scour_w4_prev + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(mast120tar, mast120) ~ scour_w4 + scour_w4_prev + 
##     strata(calving.season), data = calf.lact)
## 
##   n= 7047, number of events= 180 
## 
##                    coef exp(coef)  se(coef)      z Pr(>|z|)
## scour_w4      -0.038280  0.962444  0.385911 -0.099    0.921
## scour_w4_prev -0.006499  0.993522  0.158301 -0.041    0.967
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## scour_w4         0.9624      1.039    0.4517     2.051
## scour_w4_prev    0.9935      1.007    0.7285     1.355
## 
## Concordance= 0.5  (se = 0.018 )
## Likelihood ratio test= 0.01  on 2 df,   p=1
## Wald test            = 0.01  on 2 df,   p=1
## Score (logrank) test = 0.01  on 2 df,   p=1
cox.zph(cox)
##               chisq df    p
## scour_w4      1.594  1 0.21
## scour_w4_prev 0.285  1 0.59
## GLOBAL        1.919  2 0.38
cox.zph(cox) %>% ggcoxzph

mast120.scour.w4.cox <- cox.all.extract(model = cox, group = "w4+")
mast120.scour.w4.cox

step 5 - combining all weeks

mast120.cox <- rbind(mast120.scour.cox,
                    mast120.scour.w1.cox,
                    mast120.scour.w2.cox,
                    mast120.scour.w3.cox,
                    mast120.scour.w4.cox)
mast120.cox

Metritis within first 10 DIM

step 1 - initial model

cox <- coxph(Surv(met10tar,met10) ~ scour + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(met10tar, met10) ~ scour + strata(calving.season), 
##     data = calf.lact)
## 
##   n= 7047, number of events= 338 
## 
##          coef exp(coef) se(coef)     z Pr(>|z|)
## scour 0.05457   1.05608  0.11218 0.486    0.627
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## scour     1.056     0.9469    0.8476     1.316
## 
## Concordance= 0.506  (se = 0.013 )
## Likelihood ratio test= 0.24  on 1 df,   p=0.6
## Wald test            = 0.24  on 1 df,   p=0.6
## Score (logrank) test = 0.24  on 1 df,   p=0.6

p-value = 0.665 > 0.05; insignificant

step 2 - check assumption of proportional hazards in initial model

cox.zph(cox)
##        chisq df    p
## scour   0.93  1 0.33
## GLOBAL  0.93  1 0.33
cox.zph(cox) %>% ggcoxzph

p = 0.35 > 0.05; assumption met

step 3 - hazard ratios

met10.scour.cox <- cox.all.extract(model = cox, group = "All weeks")
met10.scour.cox

step 4 - individual weeks

cox <- coxph(Surv(met10tar,met10) ~ scour_w1 + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(met10tar, met10) ~ scour_w1 + strata(calving.season), 
##     data = calf.lact)
## 
##   n= 7047, number of events= 338 
## 
##             coef exp(coef) se(coef)      z Pr(>|z|)
## scour_w1 -0.1030    0.9021   0.2206 -0.467     0.64
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## scour_w1    0.9021      1.109    0.5855      1.39
## 
## Concordance= 0.503  (se = 0.007 )
## Likelihood ratio test= 0.22  on 1 df,   p=0.6
## Wald test            = 0.22  on 1 df,   p=0.6
## Score (logrank) test = 0.22  on 1 df,   p=0.6
cox.zph(cox)
##           chisq df    p
## scour_w1 0.0877  1 0.77
## GLOBAL   0.0877  1 0.77
cox.zph(cox) %>% ggcoxzph

met10.scour.w1.cox <- cox.all.extract(model = cox, group = "w1")
met10.scour.w1.cox
cox <- coxph(Surv(met10tar,met10) ~ scour_w2 + scour_w2_prev + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(met10tar, met10) ~ scour_w2 + scour_w2_prev + 
##     strata(calving.season), data = calf.lact)
## 
##   n= 7047, number of events= 338 
## 
##                   coef exp(coef) se(coef)      z Pr(>|z|)
## scour_w2       0.19280   1.21264  0.12583  1.532    0.125
## scour_w2_prev -0.07887   0.92416  0.22122 -0.356    0.721
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## scour_w2         1.2126     0.8246    0.9476     1.552
## scour_w2_prev    0.9242     1.0821    0.5990     1.426
## 
## Concordance= 0.52  (se = 0.013 )
## Likelihood ratio test= 2.5  on 2 df,   p=0.3
## Wald test            = 2.57  on 2 df,   p=0.3
## Score (logrank) test = 2.58  on 2 df,   p=0.3
cox.zph(cox)
##                chisq df    p
## scour_w2      0.0852  1 0.77
## scour_w2_prev 0.0872  1 0.77
## GLOBAL        0.1867  2 0.91
cox.zph(cox) %>% ggcoxzph

met10.scour.w2.cox <- cox.all.extract(model = cox, group = "w2")
met10.scour.w2.cox
cox <- coxph(Surv(met10tar,met10) ~ scour_w3 + scour_w3_prev + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(met10tar, met10) ~ scour_w3 + scour_w3_prev + 
##     strata(calving.season), data = calf.lact)
## 
##   n= 7047, number of events= 338 
## 
##                  coef exp(coef) se(coef)     z Pr(>|z|)
## scour_w3      0.07382   1.07662  0.21236 0.348    0.728
## scour_w3_prev 0.10717   1.11312  0.11889 0.901    0.367
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## scour_w3          1.077     0.9288    0.7101     1.632
## scour_w3_prev     1.113     0.8984    0.8817     1.405
## 
## Concordance= 0.512  (se = 0.013 )
## Likelihood ratio test= 0.88  on 2 df,   p=0.6
## Wald test            = 0.89  on 2 df,   p=0.6
## Score (logrank) test = 0.89  on 2 df,   p=0.6
cox.zph(cox)
##                chisq df     p
## scour_w3      3.4651  1 0.063
## scour_w3_prev 0.0225  1 0.881
## GLOBAL        3.5454  2 0.170
cox.zph(cox) %>% ggcoxzph

met10.scour.w3.cox <- cox.all.extract(model = cox, group = "w3")
met10.scour.w3.cox
cox <- coxph(Surv(met10tar,met10) ~ scour_w4 + scour_w4_prev + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(met10tar, met10) ~ scour_w4 + scour_w4_prev + 
##     strata(calving.season), data = calf.lact)
## 
##   n= 7047, number of events= 338 
## 
##                  coef exp(coef) se(coef)      z Pr(>|z|)
## scour_w4      -0.2084    0.8119   0.3069 -0.679    0.497
## scour_w4_prev  0.1070    1.1129   0.1134  0.944    0.345
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## scour_w4         0.8119     1.2317    0.4449     1.482
## scour_w4_prev    1.1129     0.8985    0.8911     1.390
## 
## Concordance= 0.513  (se = 0.013 )
## Likelihood ratio test= 1.42  on 2 df,   p=0.5
## Wald test            = 1.39  on 2 df,   p=0.5
## Score (logrank) test = 1.4  on 2 df,   p=0.5
cox.zph(cox)
##               chisq df    p
## scour_w4      0.902  1 0.34
## scour_w4_prev 0.752  1 0.39
## GLOBAL        1.707  2 0.43
cox.zph(cox) %>% ggcoxzph

met10.scour.w4.cox <- cox.all.extract(model = cox, group = "w4+")
met10.scour.w4.cox

step 5 - combining all weeks

met10.cox <- rbind(met10.scour.cox,
                    met10.scour.w1.cox,
                    met10.scour.w2.cox,
                    met10.scour.w3.cox,
                    met10.scour.w4.cox)
met10.cox

Retained foetal membranes within first 10 DIM

step 1 - initial model

cox <- coxph(Surv(rfm10tar,rfm10) ~ scour + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(rfm10tar, rfm10) ~ scour + strata(calving.season), 
##     data = calf.lact)
## 
##   n= 7047, number of events= 173 
## 
##          coef exp(coef) se(coef)      z Pr(>|z|)
## scour -0.1092    0.8966   0.1604 -0.681    0.496
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## scour    0.8966      1.115    0.6546     1.228
## 
## Concordance= 0.512  (se = 0.018 )
## Likelihood ratio test= 0.47  on 1 df,   p=0.5
## Wald test            = 0.46  on 1 df,   p=0.5
## Score (logrank) test = 0.46  on 1 df,   p=0.5

p-value = 0.557 > 0.05; insignificant

step 2 - check assumption of proportional hazards in initial model

cox.zph(cox)
##        chisq df    p
## scour  0.389  1 0.53
## GLOBAL 0.389  1 0.53
cox.zph(cox) %>% ggcoxzph

p = 0.61 > 0.05; assumption met

step 3 - hazard ratios

rfm10.scour.cox <- cox.all.extract(model = cox, group = "All weeks")
rfm10.scour.cox

step 4 - individual weeks

cox <- coxph(Surv(rfm10tar,rfm10) ~ scour_w1 + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(rfm10tar, rfm10) ~ scour_w1 + strata(calving.season), 
##     data = calf.lact)
## 
##   n= 7047, number of events= 173 
## 
##            coef exp(coef) se(coef)     z Pr(>|z|)
## scour_w1 0.3599    1.4332   0.2555 1.408    0.159
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## scour_w1     1.433     0.6977    0.8685     2.365
## 
## Concordance= 0.515  (se = 0.011 )
## Likelihood ratio test= 1.8  on 1 df,   p=0.2
## Wald test            = 1.98  on 1 df,   p=0.2
## Score (logrank) test = 2  on 1 df,   p=0.2
cox.zph(cox)
##          chisq df     p
## scour_w1  3.74  1 0.053
## GLOBAL    3.74  1 0.053
cox.zph(cox) %>% ggcoxzph

rfm10.scour.w1.cox <- cox.all.extract(model = cox, group = "w1")
rfm10.scour.w1.cox
cox <- coxph(Surv(rfm10tar,rfm10) ~ scour_w2 + scour_w2_prev + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(rfm10tar, rfm10) ~ scour_w2 + scour_w2_prev + 
##     strata(calving.season), data = calf.lact)
## 
##   n= 7047, number of events= 173 
## 
##                  coef exp(coef) se(coef)      z Pr(>|z|)
## scour_w2      -0.1384    0.8707   0.1943 -0.712    0.476
## scour_w2_prev  0.3442    1.4109   0.2564  1.342    0.179
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## scour_w2         0.8707     1.1485    0.5950     1.274
## scour_w2_prev    1.4109     0.7088    0.8536     2.332
## 
## Concordance= 0.525  (se = 0.018 )
## Likelihood ratio test= 2.32  on 2 df,   p=0.3
## Wald test            = 2.48  on 2 df,   p=0.3
## Score (logrank) test = 2.51  on 2 df,   p=0.3
cox.zph(cox)
##               chisq df     p
## scour_w2       2.16  1 0.142
## scour_w2_prev  3.72  1 0.054
## GLOBAL         5.45  2 0.066
cox.zph(cox) %>% ggcoxzph

rfm10.scour.w2.cox <- cox.all.extract(model = cox, group = "w2")
rfm10.scour.w2.cox
cox <- coxph(Surv(rfm10tar,rfm10) ~ scour_w3 + scour_w3_prev + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(rfm10tar, rfm10) ~ scour_w3 + scour_w3_prev + 
##     strata(calving.season), data = calf.lact)
## 
##   n= 7047, number of events= 173 
## 
##                   coef exp(coef) se(coef)      z Pr(>|z|)
## scour_w3      -0.07145   0.93105  0.31230 -0.229    0.819
## scour_w3_prev -0.01151   0.98855  0.17031 -0.068    0.946
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## scour_w3         0.9310      1.074    0.5048     1.717
## scour_w3_prev    0.9886      1.012    0.7080     1.380
## 
## Concordance= 0.505  (se = 0.018 )
## Likelihood ratio test= 0.06  on 2 df,   p=1
## Wald test            = 0.06  on 2 df,   p=1
## Score (logrank) test = 0.06  on 2 df,   p=1
cox.zph(cox)
##                chisq df     p
## scour_w3      3.3503  1 0.067
## scour_w3_prev 0.0198  1 0.888
## GLOBAL        3.4163  2 0.181
cox.zph(cox) %>% ggcoxzph

rfm10.scour.w3.cox <- cox.all.extract(model = cox, group = "w3")
rfm10.scour.w3.cox
cox <- coxph(Surv(rfm10tar,rfm10) ~ scour_w4 + scour_w4_prev + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(rfm10tar, rfm10) ~ scour_w4 + scour_w4_prev + 
##     strata(calving.season), data = calf.lact)
## 
##   n= 7047, number of events= 173 
## 
##                  coef exp(coef) se(coef)      z Pr(>|z|)
## scour_w4      -0.5895    0.5546   0.5062 -1.165    0.244
## scour_w4_prev -0.0609    0.9409   0.1627 -0.374    0.708
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## scour_w4         0.5546      1.803    0.2056     1.496
## scour_w4_prev    0.9409      1.063    0.6840     1.294
## 
## Concordance= 0.515  (se = 0.018 )
## Likelihood ratio test= 1.76  on 2 df,   p=0.4
## Wald test            = 1.48  on 2 df,   p=0.5
## Score (logrank) test = 1.51  on 2 df,   p=0.5
cox.zph(cox)
##               chisq df    p
## scour_w4      1.084  1 0.30
## scour_w4_prev 0.969  1 0.32
## GLOBAL        2.006  2 0.37
cox.zph(cox) %>% ggcoxzph

rfm10.scour.w4.cox <- cox.all.extract(model = cox, group = "w4+")
rfm10.scour.w4.cox

step 5 - combining all weeks

rfm10.cox <- rbind(rfm10.scour.cox,
                    rfm10.scour.w1.cox,
                    rfm10.scour.w2.cox,
                    rfm10.scour.w3.cox,
                    rfm10.scour.w4.cox)
rfm10.cox

Pregnancy within 200 DIM

step 1 - initial model

cox <- coxph(Surv(preg200tar,preg200) ~ scour + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(preg200tar, preg200) ~ scour + strata(calving.season), 
##     data = calf.lact)
## 
##   n= 7047, number of events= 5000 
## 
##           coef exp(coef) se(coef)     z Pr(>|z|)
## scour -0.03069   0.96977  0.02951 -1.04    0.298
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## scour    0.9698      1.031    0.9153     1.028
## 
## Concordance= 0.505  (se = 0.004 )
## Likelihood ratio test= 1.08  on 1 df,   p=0.3
## Wald test            = 1.08  on 1 df,   p=0.3
## Score (logrank) test = 1.08  on 1 df,   p=0.3

p-value = 0.436 > 0.05; insignificant

step 2 - check assumption of proportional hazards in initial model

cox.zph(cox)
##        chisq df    p
## scour   0.39  1 0.53
## GLOBAL  0.39  1 0.53
cox.zph(cox) %>% ggcoxzph

p = 0.61 > 0.05; assumption met

step 3 - hazard ratios

preg200.scour.cox <- cox.all.extract(model = cox, group = "All weeks")
preg200.scour.cox

step 4 - individual weeks

cox <- coxph(Surv(preg200tar,preg200) ~ scour_w1 + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(preg200tar, preg200) ~ scour_w1 + strata(calving.season), 
##     data = calf.lact)
## 
##   n= 7047, number of events= 5000 
## 
##              coef exp(coef) se(coef)      z Pr(>|z|)  
## scour_w1 -0.11394   0.89231  0.05526 -2.062   0.0392 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## scour_w1    0.8923      1.121    0.8007    0.9944
## 
## Concordance= 0.507  (se = 0.002 )
## Likelihood ratio test= 4.39  on 1 df,   p=0.04
## Wald test            = 4.25  on 1 df,   p=0.04
## Score (logrank) test = 4.26  on 1 df,   p=0.04
cox.zph(cox)
##          chisq df       p
## scour_w1  11.5  1 0.00071
## GLOBAL    11.5  1 0.00071
cox.zph(cox) %>% ggcoxzph

preg200.scour.w1.cox <- cox.all.extract(model = cox, group = "w1")
preg200.scour.w1.cox
cox <- coxph(Surv(preg200tar,preg200) ~ scour_w2 + scour_w2_prev + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(preg200tar, preg200) ~ scour_w2 + scour_w2_prev + 
##     strata(calving.season), data = calf.lact)
## 
##   n= 7047, number of events= 5000 
## 
##                   coef exp(coef) se(coef)      z Pr(>|z|)  
## scour_w2      -0.02332   0.97695  0.03484 -0.669   0.5032  
## scour_w2_prev -0.11673   0.88983  0.05542 -2.106   0.0352 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## scour_w2         0.9769      1.024    0.9125    1.0460
## scour_w2_prev    0.8898      1.124    0.7982    0.9919
## 
## Concordance= 0.509  (se = 0.003 )
## Likelihood ratio test= 4.84  on 2 df,   p=0.09
## Wald test            = 4.7  on 2 df,   p=0.1
## Score (logrank) test = 4.7  on 2 df,   p=0.1
cox.zph(cox)
##                 chisq df       p
## scour_w2       0.0118  1 0.91333
## scour_w2_prev 11.4485  1 0.00072
## GLOBAL        11.4743  2 0.00322
cox.zph(cox) %>% ggcoxzph

preg200.scour.w2.cox <- cox.all.extract(model = cox, group = "w2")
preg200.scour.w2.cox
cox <- coxph(Surv(preg200tar,preg200) ~ scour_w3 + scour_w3_prev + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(preg200tar, preg200) ~ scour_w3 + scour_w3_prev + 
##     strata(calving.season), data = calf.lact)
## 
##   n= 7047, number of events= 5000 
## 
##                     coef  exp(coef)   se(coef)      z Pr(>|z|)
## scour_w3       0.0004126  1.0004127  0.0566074  0.007    0.994
## scour_w3_prev -0.0509608  0.9503159  0.0318255 -1.601    0.109
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## scour_w3         1.0004     0.9996    0.8954     1.118
## scour_w3_prev    0.9503     1.0523    0.8928     1.011
## 
## Concordance= 0.508  (se = 0.004 )
## Likelihood ratio test= 2.59  on 2 df,   p=0.3
## Wald test            = 2.57  on 2 df,   p=0.3
## Score (logrank) test = 2.57  on 2 df,   p=0.3
cox.zph(cox)
##               chisq df     p
## scour_w3       5.07  1 0.024
## scour_w3_prev  1.97  1 0.161
## GLOBAL         6.70  2 0.035
cox.zph(cox) %>% ggcoxzph

preg200.scour.w3.cox <- cox.all.extract(model = cox, group = "w3")
preg200.scour.w3.cox
cox <- coxph(Surv(preg200tar,preg200) ~ scour_w4 + scour_w4_prev + strata(calving.season),data = calf.lact)
summary(cox)
## Call:
## coxph(formula = Surv(preg200tar, preg200) ~ scour_w4 + scour_w4_prev + 
##     strata(calving.season), data = calf.lact)
## 
##   n= 7047, number of events= 5000 
## 
##                   coef exp(coef) se(coef)      z Pr(>|z|)
## scour_w4       0.03780   1.03852  0.07179  0.526    0.599
## scour_w4_prev -0.03860   0.96213  0.03016 -1.280    0.201
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## scour_w4         1.0385     0.9629    0.9022     1.195
## scour_w4_prev    0.9621     1.0394    0.9069     1.021
## 
## Concordance= 0.504  (se = 0.004 )
## Likelihood ratio test= 1.96  on 2 df,   p=0.4
## Wald test            = 1.95  on 2 df,   p=0.4
## Score (logrank) test = 1.95  on 2 df,   p=0.4
cox.zph(cox)
##                 chisq df     p
## scour_w4      4.97170  1 0.026
## scour_w4_prev 0.00834  1 0.927
## GLOBAL        4.99605  2 0.082
cox.zph(cox) %>% ggcoxzph

preg200.scour.w4.cox <- cox.all.extract(model = cox, group = "w4+")
preg200.scour.w4.cox

step 5 - combining all weeks

preg200.cox <- rbind(preg200.scour.cox,
                    preg200.scour.w1.cox,
                    preg200.scour.w2.cox,
                    preg200.scour.w3.cox,
                    preg200.scour.w4.cox)
preg200.cox

###Table combining calf.lact hazard ratios

calf.lact.hr.all <- rbind(rem300.cox,
                          mast120.cox,
                          met10.cox,
                          rfm10.cox,
                          preg200.cox)

calf.lact.hr.all
write.csv(calf.lact.hr.all,'calf.lact.hr.all.csv')

None significant

2. Analysis for continuous outcomes in calf.lact dataset = Linear regression

Variables = me305, maxlscc, c1age

continuous.table <- function(predictor,outcome){
  predictor <-  enquo(predictor)
  outcome <- enquo(outcome)

  calf.lact %>% group_by(!!predictor) %>% summarise(
    mean = mean(!!outcome, na.rm=T),
    sd = sd(!!outcome, na.rm=T))

}


linear.regression.extract <- function(predictor,outcome,covariate,dataframe){

lm <- lm(formula(paste0(substitute(outcome)," ~ ",substitute(predictor),"+", substitute(covariate))), data=dataframe)
aic <- glance(lm) %>% select(AIC)
aic <- aic$AIC %>% round(2)
tab <- tidy(lm,conf.int=T)
tab <- tab %>% slice_head(n=2) %>% filter(term != "(Intercept)")
tab$difference <- paste0(round(tab$estimate,2),
                             " (",
                             round(tab$conf.low,2),
                             " to ",
                             round(tab$conf.high,2),
                             ")")
tab$exposure <- deparse(substitute(predictor))
tab$outcome <- deparse(substitute(outcome))
tab$covariate <- deparse(substitute(covariate))
tab$p.value <- format.pval(tab$p.value,nsmall=4,digits=3)
tab$predictor <- tab$exposure
tab$aic <- aic
tab
tab1 <- emmeans(lm,formula(paste0("~",substitute(predictor)))) %>% as_tibble
colnames(tab1) <- c("predictor","emmean","se","df","lower.CL","upper.CL")
tab1$estimate <- paste0(round(tab1$emmean,2)," (",round(tab1$lower.CL,2)," to ",round(tab1$upper.CL,2),")")
adjusted.mean.exp1 <- tab1$estimate[tab1$predictor==1]
adjusted.mean.exp0 <- tab1$estimate[tab1$predictor==0]

tab$unexposed <- adjusted.mean.exp0
tab$exposed <- adjusted.mean.exp1

tab %>% select(outcome,predictor,covariate,exposed,unexposed,difference,p.value)
}

Use this code to look at the descriptive statistics (i.e. crude mean for each group)

continuous.table(predictor=scour,
                 outcome = maxlscc)
continuous.table(predictor=scour,
                 outcome = me305)
continuous.table(predictor=scour,
                 outcome = c1age)

##Linear regression with a dichotomous predictor (e.g. scour), and no covariates.

model_maxlscc <- lm(maxlscc ~ scour, data = calf.lact) 
summary(model_maxlscc)
## 
## Call:
## lm(formula = maxlscc ~ scour, data = calf.lact)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8236 -0.8179 -0.0811  0.7562  3.7760 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.43300    0.01944 279.473   <2e-16 ***
## scour        0.07615    0.03226   2.361   0.0183 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.217 on 6154 degrees of freedom
##   (891 observations deleted due to missingness)
## Multiple R-squared:  0.0009048,  Adjusted R-squared:  0.0007424 
## F-statistic: 5.573 on 1 and 6154 DF,  p-value: 0.01827
model_me305 <- lm(me305 ~ scour, data = calf.lact) 
summary(model_me305)
## 
## Call:
## lm(formula = me305 ~ scour, data = calf.lact)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13481.1  -1660.7    249.2   1958.9  10159.5 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 13481.14      46.13 292.245   <2e-16 ***
## scour         -80.63      76.54  -1.053    0.292    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2888 on 6154 degrees of freedom
##   (891 observations deleted due to missingness)
## Multiple R-squared:  0.0001803,  Adjusted R-squared:  1.784e-05 
## F-statistic:  1.11 on 1 and 6154 DF,  p-value: 0.2922
model_c1age <- lm(c1age ~ scour, data = calf.lact) 
summary(model_c1age)
## 
## Call:
## lm(formula = c1age ~ scour, data = calf.lact)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -115.63  -36.44  -18.63   16.56  390.56 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 713.4413     0.8248 865.002   <2e-16 ***
## scour         3.1891     1.3642   2.338   0.0194 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 55.15 on 7045 degrees of freedom
## Multiple R-squared:  0.0007752,  Adjusted R-squared:  0.0006333 
## F-statistic: 5.465 on 1 and 7045 DF,  p-value: 0.01943

###Checking model assumptions of homoscedasticity and normality of residuals with no covariates

model.diag.metrics <- augment(model_maxlscc)

par(mfrow = c(2, 2))
plot(model_maxlscc)

#maxlscc residuals have even scatter, normally distributed and no fanning; assumptions met

model.diag.metrics <- augment(model_me305)

par(mfrow = c(2, 2))
plot(model_me305)

#me305 residuals have even scatter, normally distributed and no fanning; assumptions met

model.diag.metrics <- augment(model_c1age)

par(mfrow = c(2, 2))
plot(model_c1age)

#c1age residuals no fanning but have uneven scatter, not normally distributed on Q-Q plot; assumptions violated

###Transforming c1age model and rechecking assumptions

model_c1age <- lm(log(c1age) ~ scour, data = calf.lact)
model_c1age <- lm(log10(c1age) ~ scour, data = calf.lact)
model_c1age <- lm(c1age ~ scour, data = calf.lact) 
summary(model_c1age)
## 
## Call:
## lm(formula = c1age ~ scour, data = calf.lact)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -115.63  -36.44  -18.63   16.56  390.56 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 713.4413     0.8248 865.002   <2e-16 ***
## scour         3.1891     1.3642   2.338   0.0194 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 55.15 on 7045 degrees of freedom
## Multiple R-squared:  0.0007752,  Adjusted R-squared:  0.0006333 
## F-statistic: 5.465 on 1 and 7045 DF,  p-value: 0.01943
#For log model
exp(0.004347 + 6.567389) - exp(6.567389)
## [1] 3.099665
#For log10 model
10^(2.8521808 + 0.0018877) - 10^(2.8521808)
## [1] 3.099372
model.diag.metrics <- augment(model_c1age)

par(mfrow = c(2, 2))
plot(model_c1age)

#still violated
hist(calf.lact$c1age %>% log,breaks="FD")

##Linear regression with a dichotomous predictor (e.g. scour) with calving.season as a covariate

model_maxlscc <- lm(maxlscc ~ scour + calving.season, data = calf.lact) 
summary(model_maxlscc)
## 
## Call:
## lm(formula = maxlscc ~ scour + calving.season, data = calf.lact)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6878 -0.8286 -0.0778  0.7575  3.8987 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           5.545313   0.034419 161.110  < 2e-16 ***
## scour                 0.075219   0.032155   2.339 0.019353 *  
## calving.seasonSpring -0.162882   0.044009  -3.701 0.000217 ***
## calving.seasonSummer -0.003565   0.045640  -0.078 0.937748    
## calving.seasonWinter -0.248026   0.044019  -5.635 1.83e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.213 on 6151 degrees of freedom
##   (891 observations deleted due to missingness)
## Multiple R-squared:  0.008512,   Adjusted R-squared:  0.007867 
## F-statistic:  13.2 on 4 and 6151 DF,  p-value: 1.038e-10
model_me305 <- lm(me305 ~ scour + calving.season, data = calf.lact) 
summary(model_me305)
## 
## Call:
## lm(formula = me305 ~ scour + calving.season, data = calf.lact)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13284.5  -1634.3    245.5   1956.0  10614.3 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          13284.46      81.48 163.039  < 2e-16 ***
## scour                  -78.31      76.12  -1.029 0.303621    
## calving.seasonSpring   388.83     104.18   3.732 0.000192 ***
## calving.seasonSummer  -260.43     108.04  -2.410 0.015962 *  
## calving.seasonWinter   558.62     104.21   5.361 8.59e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2871 on 6151 degrees of freedom
##   (891 observations deleted due to missingness)
## Multiple R-squared:  0.01249,    Adjusted R-squared:  0.01185 
## F-statistic: 19.45 on 4 and 6151 DF,  p-value: 6.388e-16
model_c1age <- lm(c1age ~ scour + calving.season, data = calf.lact) 
summary(model_c1age)
## 
## Call:
## lm(formula = c1age ~ scour + calving.season, data = calf.lact)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -120.19  -36.09  -18.09   16.91  385.91 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          711.0865     1.4483 490.984  < 2e-16 ***
## scour                  3.0997     1.3633   2.274 0.023014 *  
## calving.seasonSpring   1.9183     1.8878   1.016 0.309582    
## calving.seasonSummer   7.0018     1.8865   3.712 0.000208 ***
## calving.seasonWinter   0.4564     1.8753   0.243 0.807740    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 55.09 on 7042 degrees of freedom
## Multiple R-squared:  0.003325,   Adjusted R-squared:  0.002759 
## F-statistic: 5.874 on 4 and 7042 DF,  p-value: 0.0001025

###Checking model assumptions of homoscedasticity and normality of residuals with calving.season as a covariate

model.diag.metrics <- augment(model_maxlscc)

par(mfrow = c(2, 2))
plot(model_maxlscc)

#maxlscc residuals have even scatter, normally distributed and no fanning; assumptions met

model.diag.metrics <- augment(model_me305)

par(mfrow = c(2, 2))
plot(model_me305)

#me305 residuals have even scatter, normally distributed and no fanning; assumptions met

model.diag.metrics <- augment(model_c1age)

par(mfrow = c(2, 2))
plot(model_c1age)

#c1age residuals no fanning but have uneven scatter, not normally distributed on Q-Q plot; assumptions violated

##Linear regression with a dichotomous predictor (e.g. scour) with enroll.season as a covariate

model_maxlscc <- lm(maxlscc ~ scour + enroll.season, data = calf.lact) 
summary(model_maxlscc)
## 
## Call:
## lm(formula = maxlscc ~ scour + enroll.season, data = calf.lact)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7023 -0.8138 -0.0820  0.7636  3.8098 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          5.54164    0.03395 163.238  < 2e-16 ***
## scour                0.07841    0.03221   2.434  0.01495 *  
## enroll.seasonSpring -0.22985    0.04309  -5.334 9.93e-08 ***
## enroll.seasonSummer -0.05682    0.04537  -1.252  0.21052    
## enroll.seasonWinter -0.12421    0.04442  -2.796  0.00519 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.214 on 6151 degrees of freedom
##   (891 observations deleted due to missingness)
## Multiple R-squared:  0.00606,    Adjusted R-squared:  0.005414 
## F-statistic: 9.375 on 4 and 6151 DF,  p-value: 1.494e-07
model_me305 <- lm(me305 ~ scour + enroll.season, data = calf.lact) 
summary(model_me305)
## 
## Call:
## lm(formula = me305 ~ scour + enroll.season, data = calf.lact)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13226.0  -1648.9    271.3   1942.2  10423.4 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         13226.03      80.55 164.198  < 2e-16 ***
## scour                 -89.41      76.42  -1.170   0.2421    
## enroll.seasonSpring   566.18     102.24   5.538 3.19e-08 ***
## enroll.seasonSummer   182.70     107.65   1.697   0.0897 .  
## enroll.seasonWinter   221.38     105.40   2.100   0.0357 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2881 on 6151 degrees of freedom
##   (891 observations deleted due to missingness)
## Multiple R-squared:  0.005498,   Adjusted R-squared:  0.004852 
## F-statistic: 8.502 on 4 and 6151 DF,  p-value: 7.744e-07
model_c1age <- lm(c1age ~ scour + enroll.season, data = calf.lact) 
summary(model_c1age)
## 
## Call:
## lm(formula = c1age ~ scour + enroll.season, data = calf.lact)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -112.85  -36.62  -18.68   16.66  394.89 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          710.617      1.399 508.080  < 2e-16 ***
## scour                  3.230      1.363   2.370 0.017802 *  
## enroll.seasonSpring    6.065      1.812   3.348 0.000819 ***
## enroll.seasonSummer   -1.507      1.889  -0.798 0.425009    
## enroll.seasonWinter    6.150      1.865   3.297 0.000982 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 55.05 on 7042 degrees of freedom
## Multiple R-squared:  0.004691,   Adjusted R-squared:  0.004126 
## F-statistic: 8.298 on 4 and 7042 DF,  p-value: 1.13e-06

###Checking model assumptions of homoscedasticity and normality of residuals with enroll.season as a covariate

model.diag.metrics <- augment(model_maxlscc)

par(mfrow = c(2, 2))
plot(model_maxlscc)

#maxlscc residuals have even scatter, normally distributed and no fanning; assumptions met

model.diag.metrics <- augment(model_me305)

par(mfrow = c(2, 2))
plot(model_me305)

#me305 residuals have even scatter, normally distributed and no fanning; assumptions met

model.diag.metrics <- augment(model_c1age)

par(mfrow = c(2, 2))
plot(model_c1age)

#c1age residuals no fanning but have uneven scatter, not normally distributed on Q-Q plot; assumptions violated

##Linear regression with a dichotomous predictor (e.g. scour) with enroll.season and calving.season as covariates

model_maxlscc <- lm(maxlscc ~ scour + enroll.season + calving.season, data = calf.lact)
summary(model_maxlscc)
## 
## Call:
## lm(formula = maxlscc ~ scour + enroll.season + calving.season, 
##     data = calf.lact)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6486 -0.8266 -0.0835  0.7664  3.8707 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           5.563686   0.042202 131.834  < 2e-16 ***
## scour                 0.076842   0.032170   2.389   0.0169 *  
## enroll.seasonSpring  -0.103915   0.052852  -1.966   0.0493 *  
## enroll.seasonSummer  -0.009346   0.051034  -0.183   0.8547    
## enroll.seasonWinter  -0.036745   0.049621  -0.741   0.4590    
## calving.seasonSpring -0.125199   0.053971  -2.320   0.0204 *  
## calving.seasonSummer -0.007024   0.050343  -0.140   0.8890    
## calving.seasonWinter -0.201695   0.049453  -4.079 4.59e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.213 on 6148 degrees of freedom
##   (891 observations deleted due to missingness)
## Multiple R-squared:  0.009369,   Adjusted R-squared:  0.008241 
## F-statistic: 8.306 on 7 and 6148 DF,  p-value: 3.968e-10
model_me305 <- lm(me305 ~ scour + enroll.season + calving.season, data = calf.lact) 
summary(model_me305)
## 
## Call:
## lm(formula = me305 ~ scour + enroll.season + calving.season, 
##     data = calf.lact)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13346.5  -1629.0    247.9   1964.5  10670.0 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          13346.46      99.89 133.617  < 2e-16 ***
## scour                  -85.48      76.14  -1.123  0.26160    
## enroll.seasonSpring    149.48     125.09   1.195  0.23213    
## enroll.seasonSummer     94.36     120.79   0.781  0.43473    
## enroll.seasonWinter   -152.81     117.44  -1.301  0.19327    
## calving.seasonSpring   230.75     127.74   1.806  0.07091 .  
## calving.seasonSummer  -371.00     119.15  -3.114  0.00186 ** 
## calving.seasonWinter   496.16     117.05   4.239 2.28e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2870 on 6148 degrees of freedom
##   (891 observations deleted due to missingness)
## Multiple R-squared:  0.01371,    Adjusted R-squared:  0.01259 
## F-statistic: 12.21 on 7 and 6148 DF,  p-value: 1.366e-15
model_c1age <- lm(c1age ~ scour + enroll.season + calving.season, data = calf.lact) 
summary(model_c1age)
## 
## Call:
## lm(formula = c1age ~ scour + enroll.season + calving.season, 
##     data = calf.lact)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -119.78  -35.94  -16.92   17.20  389.12 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           705.421      1.737 406.130  < 2e-16 ***
## scour                   3.089      1.357   2.277   0.0228 *  
## enroll.seasonSpring    11.322      2.211   5.120 3.14e-07 ***
## enroll.seasonSummer    -2.804      2.108  -1.330   0.1834    
## enroll.seasonWinter    13.380      2.083   6.423 1.42e-10 ***
## calving.seasonSpring    2.217      2.287   0.970   0.3323    
## calving.seasonSummer   12.268      2.068   5.933 3.12e-09 ***
## calving.seasonWinter   -5.119      2.099  -2.438   0.0148 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 54.8 on 7039 degrees of freedom
## Multiple R-squared:  0.01408,    Adjusted R-squared:  0.0131 
## F-statistic: 14.36 on 7 and 7039 DF,  p-value: < 2.2e-16

###Checking model assumptions of homoscedasticity and normality of residuals with enroll.season and calving.season as covariates

model.diag.metrics <- augment(model_maxlscc)

par(mfrow = c(2, 2))
plot(model_maxlscc)

#maxlscc residuals have even scatter, normally distributed and no fanning; assumptions met

model.diag.metrics <- augment(model_me305)

par(mfrow = c(2, 2))
plot(model_me305)

#me305 residuals have even scatter, normally distributed and no fanning; assumptions met

model.diag.metrics <- augment(model_c1age)

par(mfrow = c(2, 2))
plot(model_c1age)

#c1age residuals no fanning but have uneven scatter, not normally distributed on Q-Q plot; assumptions violated

###**Transforming c1age model and rechecking assumptions

model_c1age <- lm(log(c1age) ~ scour + enroll.season + calving.season, data = calf.lact) 
summary(model_c1age)
## 
## Call:
## lm(formula = log(c1age) ~ scour + enroll.season + calving.season, 
##     data = calf.lact)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.17806 -0.04944 -0.02176  0.02640  0.43770 
## 
## Coefficients:
##                       Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)           6.557005   0.002276 2880.800  < 2e-16 ***
## scour                 0.004183   0.001778    2.353   0.0187 *  
## enroll.seasonSpring   0.015512   0.002898    5.353 8.93e-08 ***
## enroll.seasonSummer  -0.003482   0.002762   -1.261   0.2074    
## enroll.seasonWinter   0.017561   0.002730    6.433 1.33e-10 ***
## calving.seasonSpring  0.002437   0.002997    0.813   0.4162    
## calving.seasonSummer  0.015472   0.002710    5.710 1.18e-08 ***
## calving.seasonWinter -0.006895   0.002751   -2.506   0.0122 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07182 on 7039 degrees of freedom
## Multiple R-squared:  0.01409,    Adjusted R-squared:  0.01311 
## F-statistic: 14.37 on 7 and 7039 DF,  p-value: < 2.2e-16
model.diag.metrics <- augment(model_c1age)

par(mfrow = c(2, 2))
plot(model_c1age)

#still violated

**Functions to extract multivariable linear regression output

Eg. maxlscc

maxlscc.scour <- linear.regression.extract(predictor = scour,
                          outcome = maxlscc,
                          covariate = enroll.season,
                          dataframe = calf.lact)
maxlscc.scour
#difference = 0.07* [0.01, 0.14]

maxlscc.scour_w1 <- linear.regression.extract(predictor = scour_w1,
                          outcome = maxlscc,
                          covariate = enroll.season,
                          dataframe = calf.lact)
maxlscc.scour_w1
maxlscc.scour_w2 <- linear.regression.extract(predictor = scour_w2,
                          outcome = maxlscc,
                          covariate = enroll.season,
                          dataframe = calf.lact)
maxlscc.scour_w2
maxlscc.scour_w3 <- linear.regression.extract(predictor = scour_w3,
                          outcome = maxlscc,
                          covariate = enroll.season,
                          dataframe = calf.lact)
maxlscc.scour_w3
maxlscc.scour_w4 <- linear.regression.extract(predictor = scour_w4,
                          outcome = maxlscc,
                          covariate = enroll.season,
                          dataframe = calf.lact)
maxlscc.scour_w4
maxlscc.scour.all <- rbind(maxlscc.scour,
                           maxlscc.scour_w1,
                           maxlscc.scour_w2,
                           maxlscc.scour_w3,
                           maxlscc.scour_w4)
maxlscc.scour.all

Eg. c1age

c1age.scour <- linear.regression.extract(predictor = scour,
                          outcome = c1age,
                          covariate = enroll.season,
                          dataframe = calf.lact)

c1age.scour
#difference = 3.24* [0.59, 5.9]

c1age.scour_w1 <- linear.regression.extract(predictor = scour_w1,
                          outcome = c1age,
                          covariate = enroll.season,
                          dataframe = calf.lact)
c1age.scour_w1
c1age.scour_w2 <- linear.regression.extract(predictor = scour_w2,
                          outcome = c1age,
                          covariate = enroll.season,
                          dataframe = calf.lact)
c1age.scour_w2
c1age.scour_w3 <- linear.regression.extract(predictor = scour_w3,
                          outcome = c1age,
                          covariate = enroll.season,
                          dataframe = calf.lact)
c1age.scour_w3
c1age.scour_w4 <- linear.regression.extract(predictor = scour_w4,
                          outcome = c1age,
                          covariate = enroll.season,
                          dataframe = calf.lact)
c1age.scour_w4
c1age.scour.all <- rbind(c1age.scour,
                         c1age.scour_w1,
                         c1age.scour_w2,
                         c1age.scour_w3,
                         c1age.scour_w4)
c1age.scour.all

me305

me305.scour <- linear.regression.extract(predictor = scour,
                          outcome = me305,
                          covariate = enroll.season,
                          dataframe = calf.lact)

me305.scour
#difference = -70.28 [-218.51, 77.95]

me305.scour_w1 <- linear.regression.extract(predictor = scour_w1,
                          outcome = me305,
                          covariate = enroll.season,
                          dataframe = calf.lact)
me305.scour_w1
me305.scour_w2 <- linear.regression.extract(predictor = scour_w2,
                          outcome = me305,
                          covariate = enroll.season,
                          dataframe = calf.lact)
me305.scour_w2
me305.scour_w3 <- linear.regression.extract(predictor = scour_w3,
                          outcome = me305,
                          covariate = enroll.season,
                          dataframe = calf.lact)
me305.scour_w3
me305.scour_w4 <- linear.regression.extract(predictor = scour_w4,
                          outcome = me305,
                          covariate = enroll.season,
                          dataframe = calf.lact)
me305.scour_w4
me305.scour.all <- rbind(me305.scour,
                         me305.scour_w1,
                         me305.scour_w2,
                         me305.scour_w3,
                         me305.scour_w4)
me305.scour.all

Combining calf.lact tables

calf.lact.continuous <- rbind(maxlscc.scour.all,
                              c1age.scour.all,
                              me305.scour.all)
calf.lact.continuous
write.csv(calf.lact.continuous,'calf.lact.continuous.csv')

Post Hoc considerations for confounding

Below is a directed acyclic graph that considers potential causal paths that confound the relationship between diarrhea (exposure) and survival (one outcome evaluated in this study).


Our analytical approach did not control for all possible confounding. These include backdoor paths via colostral IgG, birthweight, and previous diseases. Therefore, we recommend that readers use caution when generalising these findings.

# DAG code for DAGGITY
# dag {
# bb="0,0,1,1"
# "Colostral IgG Transfer" [latent,pos="0.214,0.441"]
# "Delayed weaning" [pos="0.577,0.601"]
# "Disease events after diarrhea" [pos="0.532,0.291"]
# "Disease events before diarrhea" [pos="0.267,0.298"]
# "Growth rates" [pos="0.480,0.466"]
# "Immune fitness" [latent,pos="0.184,0.146"]
# Birthweight [latent,pos="0.176,0.595"]
# Diarrhea [exposure,pos="0.390,0.436"]
# Survival [outcome,pos="0.816,0.448"]
# "Colostral IgG Transfer" -> "Disease events after diarrhea"
# "Colostral IgG Transfer" -> "Disease events before diarrhea"
# "Colostral IgG Transfer" -> Diarrhea
# "Disease events after diarrhea" -> "Delayed weaning"
# "Disease events after diarrhea" -> "Growth rates"
# "Disease events after diarrhea" -> Survival
# "Disease events before diarrhea" -> "Disease events after diarrhea"
# "Disease events before diarrhea" -> "Growth rates"
# "Disease events before diarrhea" -> Diarrhea
# "Disease events before diarrhea" -> Survival
# "Growth rates" -> "Delayed weaning"
# "Immune fitness" -> "Disease events after diarrhea"
# "Immune fitness" -> "Disease events before diarrhea"
# "Immune fitness" -> Diarrhea
# Birthweight -> "Colostral IgG Transfer"
# Birthweight -> "Growth rates"
# Birthweight -> Diarrhea
# Birthweight -> Survival
# Diarrhea -> "Delayed weaning"
# Diarrhea -> "Disease events after diarrhea"
# Diarrhea -> "Growth rates"
# Diarrhea -> Survival
# }

Description of transfer of passive immunity in this herd

check <- calf.lact %>% left_join(
  igg %>% mutate(id = cow_id) %>% select(id,sts)
)
## Joining with `by = join_by(id)`
igg_filt <- igg %>% filter(cow_id %in% calf.survival$id) %>% distinct

hist(igg_filt$sts, breaks = "FD")

# igg_filt %>% tabyl(class) %>%
#   adorn_totals("row") %>%
#   adorn_pct_formatting(digits = 1) %>%
#   adorn_ns

igg_import <- igg %>% mutate(id = cow_id) %>% select(id,sts,class)

calf.survival.sens <- calf.survival %>% left_join(igg_import, by = "id") %>%
  filter(!is.na(sts))



ids <- calf.survival.sens %>% filter(dead21_to_100 == 1)


calf.lact.sens <- calf.lact %>% left_join(igg_import, by = "id") %>%
  filter(!is.na(sts))
# cox <- coxph(Surv(dead20tar,dead20) ~ scour.survival + strata(birth.season),data = calf.survival)
# summary(cox)
# 
# cox <- coxph(Surv(dead20tar,dead20) ~ scour.survival + strata(birth.season),data = calf.survival.sens)
# summary(cox)
# cox <- coxph(Surv(c1tar,c1) ~ scour.survival + strata(birth.season),data = calf.survival)
# summary(cox)
# 
# cox <- coxph(Surv(c1tar,c1) ~ scour.survival + strata(birth.season), data = calf.survival.sens)
# summary(cox)
# 
# cox <- coxph(Surv(c1tar,c1) ~ scour.survival + strata(birth.season) + factor(class), data = calf.survival.sens)
# summary(cox)
# 
# cox <- coxph(Surv(c1tar,c1) ~ scour.survival + strata(birth.season) + sts + I(sts^2), data = calf.survival.sens)
# summary(cox)
# cox <- coxph(Surv(rem300tar,rem300) ~ scour + strata(calving.season),data = calf.lact)
# summary(cox)
# 
# cox <- coxph(Surv(rem300tar,rem300) ~ scour + strata(calving.season),data = calf.lact.sens)
# summary(cox)
# 
# cox <- coxph(Surv(rem300tar,rem300) ~ scour + factor(class) + strata(calving.season),data = calf.lact.sens)
# summary(cox)
# 
# cox <- coxph(Surv(rem300tar,rem300) ~ sts + I(sts^2),data = calf.lact.sens)
# summary(cox)

Possible confounding by weaning age (requested by reviewers)

Relationship between scour (exposure) and weaning age (outcome). Calves with diarrhea were weaned 0.3 days later (-0.04 to 0.6). This indicates a weak relationship.

lm <- lm(weanage ~ scour, calf.lact)
lm %>% tidy(conf.int = T)

Time to first calving

Relationship between weaning age (exposure) and calving age. HR = 0.998 (95% CI: 0.995 to 1.001). No strong association observed. Effect estimate for scour did not change when including weaning age in the model.

cox <- coxph(Surv(c1tar,c1) ~ weanage + strata(birth.season),data = calf.survival)
cox %>% tidy(exp = T, conf.int =T)
cox <- coxph(Surv(c1tar,c1) ~ scour.survival + strata(birth.season),data = calf.survival)
cox %>% tidy(exp = T, conf.int =T)
cox <- coxph(Surv(c1tar,c1) ~ weanage + scour.survival + strata(birth.season),data = calf.survival)
cox %>% tidy(exp = T, conf.int =T)

Time to pregnancy during first lactation

cox <- coxph(Surv(preg200tar,preg200) ~ weanage + strata(calving.season),data = calf.lact)
cox %>% tidy(conf.int =T)
cox <- coxph(Surv(preg200tar,preg200) ~ scour + strata(calving.season),data = calf.lact)
cox %>% tidy(conf.int =T)
cox <- coxph(Surv(preg200tar,preg200) ~ weanage + scour + strata(calving.season),data = calf.lact)
cox %>% tidy(conf.int =T)
cox %>% tidy(exp = T, conf.int =T)

Time to culling/death during first lactation

No evidence of confounding

cox <- coxph(Surv(rem300tar,rem300) ~ scour + strata(calving.season),data = calf.lact)
cox %>% tidy(conf.int =T)
cox <- coxph(Surv(rem300tar,rem300) ~ weanage + strata(calving.season),data = calf.lact)
cox %>% tidy(conf.int =T)
cox <- coxph(Surv(rem300tar,rem300) ~ scour + weanage + strata(calving.season),data = calf.lact)
cox %>% tidy(conf.int =T)