Case Control study investigating Coxiellosis and Abortion

Acknowledgements

Authors and affiliations

Thank you

Funding

Funding bodies include:

Study objectives

The primary objective of this study was to investigate if coxiella exposure status is associated with pregnancy loss in cattle. [reword as necessary]

Import data and start packages

Code

coxiella <- readRDS("~/Dropbox/R backup/2023/coxiella abortion/coxiella abortion case control_final.rds")
#coxiella_old <- readRDS("~/Dropbox/R backup/2023/coxiella abortion/coxiella abortion case control_updated.rds")

#write.csv(coxiella,'coxiella.csv')
#library(foreign)
#write.dta(data_out, "coxiella abortion data.dta")

library(tidymodels) ; library(tidyverse) ; library(janitor); library(survival); library(summarytools); library(DT); library(table1); library(epiR)

Dataframe summary

Code

print(summarytools::dfSummary(coxiella %>% select(-id) ,valid.col=FALSE, graph.magnif=0.8, style="grid"),method = "render")

Data Frame Summary

coxiella

Dimensions: 377 x 17
Duplicates: 0

lact [numeric]

Mean (sd) : 2.5 (1.2)
min ≤ med ≤ max:
1 ≤ 2 ≤ 7
IQR (CV) : 1 (0.5)

1	:	58	(	15.4%	)
2	:	189	(	50.3%	)
3	:	59	(	15.7%	)
4	:	35	(	9.3%	)
5	:	25	(	6.6%	)
6	:	9	(	2.4%	)
7	:	1	(	0.3%	)

1 (0.3%)

date_at_enrollment [Date]

1. 2022-07-04
2. 2022-07-19
3. 2022-08-02
4. 2022-10-13
5. 2023-08-22
6. 2023-08-23
7. 2023-11-20

26	(	6.9%	)
46	(	12.2%	)
45	(	11.9%	)
49	(	13.0%	)
65	(	17.2%	)
18	(	4.8%	)
128	(	34.0%	)

0 (0.0%)

conception_date [Date]

min : 2021-12-21
med : 2023-06-02
max : 2023-10-13
range : 1y 9m 22d

135 distinct values

0 (0.0%)

calving_to_enrollment_d [numeric]

Mean (sd) : 162 (44.4)
min ≤ med ≤ max:
106 ≤ 149 ≤ 321
IQR (CV) : 56 (0.3)

134 distinct values

1 (0.3%)

conception_to_enrollment_d [numeric]

Mean (sd) : 65.4 (26.8)
min ≤ med ≤ max:
36 ≤ 59 ≤ 253
IQR (CV) : 27 (0.4)

75 distinct values

0 (0.0%)

conception_to_preg_dx_1 [numeric]

Mean (sd) : 45.2 (25.4)
min ≤ med ≤ max:
30 ≤ 34 ≤ 209
IQR (CV) : 14 (0.6)

44 distinct values

0 (0.0%)

linked_sample [numeric]

Mean (sd) : 76.7 (45.5)
min ≤ med ≤ max:
1 ≤ 81 ≤ 199
IQR (CV) : 81 (0.6)

135 distinct values

0 (0.0%)

case [numeric]

Min : 0
Mean : 0.3
Max : 1

0	:	251	(	66.6%	)
1	:	126	(	33.4%	)

0 (0.0%)

coxiella_test_status [numeric]

Min : 0
Mean : 0.3
Max : 1

0	:	257	(	68.2%	)
1	:	120	(	31.8%	)

0 (0.0%)

elisa_value [numeric]

Mean (sd) : 53 (71)
min ≤ med ≤ max:
-15.6 ≤ 15.9 ≤ 250.2
IQR (CV) : 93.3 (1.3)

374 distinct values

0 (0.0%)

et [integer]

Min : 0
Mean : 0.5
Max : 1

0	:	184	(	48.8%	)
1	:	193	(	51.2%	)

0 (0.0%)

exclude_control [numeric]

Min : 0
Mean : 0.1
Max : 1

0	:	335	(	88.9%	)
1	:	42	(	11.1%	)

0 (0.0%)

exclude_case [numeric]

Min : 0
Mean : 0
Max : 1

0	:	375	(	99.5%	)
1	:	2	(	0.5%	)

0 (0.0%)

exclude_duplicate [numeric]

Min : 0
Mean : 0
Max : 1

0	:	369	(	97.9%	)
1	:	8	(	2.1%	)

0 (0.0%)

exclude_incomplete_pair [numeric]

Min : 0
Mean : 0.1
Max : 1

0	:	336	(	89.1%	)
1	:	41	(	10.9%	)

0 (0.0%)

reason_exclude_control [character]

1. aborted prior to current
2. abortion event within 50d
3. pregnancy not confirmed

19	(	45.2%	)
15	(	35.7%	)
8	(	19.0%	)

335 (88.9%)

reason_exclude_case [character]

1. no abortion event in reco

2	(	100.0%	)

375 (99.5%)

Generated by summarytools 1.0.1 (R version 4.3.2)
2024-06-17

Reasons for exclusion

Show duplicates

Code

get_dupes(coxiella,id) %>% select(date_at_enrollment,case,linked_sample,exclude_duplicate) %>% datatable(filter = "top")

Cases

Code

coxiella %>% tabyl(reason_exclude_case)

reason_exclude_case	n	percent	valid_percent
no abortion event in records	2	0.005305	1
NA	375	0.994695	NA

Controls

Code

coxiella %>% tabyl(reason_exclude_control)

reason_exclude_control	n	percent	valid_percent
aborted prior to current conception (this lactation)	19	0.0503979	0.4523810
abortion event within 50d of enrollment	15	0.0397878	0.3571429
pregnancy not confirmed	8	0.0212202	0.1904762
NA	335	0.8885942	NA

Incomplete pairs

These are situations when exclusion of 1 case or both controls leaves an incomplete matched set. These have been excluded.

Code

coxiella %>% filter(exclude_control == 0,exclude_case == 0,exclude_duplicate == 0) %>% tabyl(exclude_incomplete_pair)

exclude_incomplete_pair	n	percent
0	302	0.9292308
1	23	0.0707692

Code

coxiella %>% filter(exclude_incomplete_pair == 1) %>% select(date_at_enrollment,case,linked_sample,exclude_control,exclude_case,exclude_duplicate) %>% datatable(filter = "top")

Code

coxiella %>% filter(exclude_incomplete_pair == 1,
                    exclude_control == 0,exclude_case == 0,exclude_duplicate == 0) %>% tabyl(case)

case	n	percent
0	14	0.6086957
1	9	0.3913043

Remove exclusions from dataset

Code

coxiella_filtered <- coxiella %>% filter(exclude_control == 0,
                                         exclude_case == 0,
                                         exclude_duplicate == 0,
                                         exclude_incomplete_pair == 0)

Dataframe summary of new dataset

Code

print(summarytools::dfSummary(coxiella_filtered %>% select(-id),valid.col=FALSE, graph.magnif=0.8, style="grid"),method = "render")

Data Frame Summary

coxiella_filtered

Dimensions: 302 x 17
Duplicates: 0

lact [numeric]

Mean (sd) : 2.5 (1.3)
min ≤ med ≤ max:
1 ≤ 2 ≤ 7
IQR (CV) : 1 (0.5)

1	:	55	(	18.3%	)
2	:	144	(	47.8%	)
3	:	45	(	15.0%	)
4	:	28	(	9.3%	)
5	:	19	(	6.3%	)
6	:	9	(	3.0%	)
7	:	1	(	0.3%	)

1 (0.3%)

date_at_enrollment [Date]

1. 2022-07-04
2. 2022-07-19
3. 2022-08-02
4. 2022-10-13
5. 2023-08-22
6. 2023-08-23
7. 2023-11-20

25	(	8.3%	)
44	(	14.6%	)
38	(	12.6%	)
47	(	15.6%	)
52	(	17.2%	)
15	(	5.0%	)
81	(	26.8%	)

0 (0.0%)

conception_date [Date]

min : 2022-04-08
med : 2022-09-02
max : 2023-10-13
range : 1y 6m 5d

108 distinct values

0 (0.0%)

calving_to_enrollment_d [numeric]

Mean (sd) : 156 (39.6)
min ≤ med ≤ max:
106 ≤ 147 ≤ 315
IQR (CV) : 47 (0.3)

109 distinct values

1 (0.3%)

conception_to_enrollment_d [numeric]

Mean (sd) : 60.6 (15.1)
min ≤ med ≤ max:
36 ≤ 58 ≤ 111
IQR (CV) : 23.8 (0.2)

55 distinct values

0 (0.0%)

conception_to_preg_dx_1 [numeric]

Mean (sd) : 40.6 (15.3)
min ≤ med ≤ max:
30 ≤ 34 ≤ 82
IQR (CV) : 7 (0.4)

29 distinct values

0 (0.0%)

linked_sample [numeric]

Mean (sd) : 69.9 (44.4)
min ≤ med ≤ max:
1 ≤ 63 ≤ 150
IQR (CV) : 77 (0.6)

115 distinct values

0 (0.0%)

case [numeric]

Min : 0
Mean : 0.4
Max : 1

0	:	187	(	61.9%	)
1	:	115	(	38.1%	)

0 (0.0%)

coxiella_test_status [numeric]

Min : 0
Mean : 0.3
Max : 1

0	:	203	(	67.2%	)
1	:	99	(	32.8%	)

0 (0.0%)

elisa_value [numeric]

Mean (sd) : 53.8 (70.5)
min ≤ med ≤ max:
-15.6 ≤ 17.1 ≤ 250.2
IQR (CV) : 101.6 (1.3)

300 distinct values

0 (0.0%)

et [integer]

Min : 0
Mean : 0.5
Max : 1

0	:	145	(	48.0%	)
1	:	157	(	52.0%	)

0 (0.0%)

exclude_control [numeric]

1 distinct value

0	:	302	(	100.0%	)

0 (0.0%)

exclude_case [numeric]

1 distinct value

0	:	302	(	100.0%	)

0 (0.0%)

exclude_duplicate [numeric]

1 distinct value

0	:	302	(	100.0%	)

0 (0.0%)

exclude_incomplete_pair [numeric]

1 distinct value

0	:	302	(	100.0%	)

0 (0.0%)

reason_exclude_control [character]

All NA's

302 (100.0%)

reason_exclude_case [character]

All NA's

302 (100.0%)

Generated by summarytools 1.0.1 (R version 4.3.2)
2024-06-17

Compare cases and controls

Code

coxiella_filtered <- coxiella_filtered %>% mutate(
  conception_to_enrollment = as.numeric(date_at_enrollment - conception_date),
  exposure_label = case_when(
    coxiella_test_status == 1 ~ "Coxiella seropositive",
    coxiella_test_status == 0 ~ "Coxiella seronegative"
  ),
  case_label = case_when(
    case == 1 ~ "Case",
    case == 0 ~ "Control"
  ))
  

table1(~ factor(exposure_label) + elisa_value + factor(date_at_enrollment) + factor(et) + factor(lact) + 
         conception_to_enrollment_d + calving_to_enrollment_d + conception_to_preg_dx_1 | factor(case_label), data=coxiella_filtered)

	Case (N=115)	Control (N=187)	Overall (N=302)
factor(exposure_label)
Coxiella seronegative	77 (67.0%)	126 (67.4%)	203 (67.2%)
Coxiella seropositive	38 (33.0%)	61 (32.6%)	99 (32.8%)
elisa_value
Mean (SD)	49.4 (66.6)	56.5 (72.8)	53.8 (70.5)
Median [Min, Max]	15.5 [-15.6, 236]	17.8 [-14.3, 250]	17.1 [-15.6, 250]
factor(date_at_enrollment)
2022-07-04	9 (7.8%)	16 (8.6%)	25 (8.3%)
2022-07-19	16 (13.9%)	28 (15.0%)	44 (14.6%)
2022-08-02	14 (12.2%)	24 (12.8%)	38 (12.6%)
2022-10-13	18 (15.7%)	29 (15.5%)	47 (15.6%)
2023-08-22	20 (17.4%)	32 (17.1%)	52 (17.2%)
2023-08-23	6 (5.2%)	9 (4.8%)	15 (5.0%)
2023-11-20	32 (27.8%)	49 (26.2%)	81 (26.8%)
factor(et)
0	44 (38.3%)	101 (54.0%)	145 (48.0%)
1	71 (61.7%)	86 (46.0%)	157 (52.0%)
factor(lact)
1	20 (17.4%)	35 (18.7%)	55 (18.2%)
2	55 (47.8%)	89 (47.6%)	144 (47.7%)
3	18 (15.7%)	27 (14.4%)	45 (14.9%)
4	10 (8.7%)	18 (9.6%)	28 (9.3%)
5	8 (7.0%)	11 (5.9%)	19 (6.3%)
6	3 (2.6%)	6 (3.2%)	9 (3.0%)
7	1 (0.9%)	0 (0%)	1 (0.3%)
Missing	0 (0%)	1 (0.5%)	1 (0.3%)
conception_to_enrollment_d
Mean (SD)	61.3 (15.2)	60.2 (15.1)	60.6 (15.1)
Median [Min, Max]	59.0 [38.0, 109]	58.0 [36.0, 111]	58.0 [36.0, 111]
calving_to_enrollment_d
Mean (SD)	154 (36.3)	157 (41.6)	156 (39.6)
Median [Min, Max]	148 [111, 302]	145 [106, 315]	147 [106, 315]
Missing	0 (0%)	1 (0.5%)	1 (0.3%)
conception_to_preg_dx_1
Mean (SD)	36.5 (10.9)	43.2 (17.1)	40.6 (15.3)
Median [Min, Max]	33.0 [30.0, 82.0]	34.0 [30.0, 80.0]	34.0 [30.0, 82.0]

Comparison of elisa values between cases and controls

Code

my_plot <- ggplot(coxiella_filtered, aes(x = elisa_value, fill = factor(case_label))) +
  geom_density(alpha = 0.5) +
  geom_vline(xintercept = 40, linetype = "dashed", color = "black") +
  labs(title = "",
       x = "Antibody ELISA units",
       y = "Density") +
  theme_minimal() +
  theme(legend.position = c(0.85, 0.85),  # Move legend to top right inside the plot
        legend.title = element_blank(),  # Remove title from legend
        panel.border = element_rect(color = "black", fill = NA),  # Add a black border around the plot
        text = element_text(family = "Times New Roman"))  # Set font family to Times New Roman


ggsave("elisa result.tiff", plot = my_plot, device = "tiff", dpi = 200, width = 5, height = 5)

my_plot

Compare seropositive and seronegative

Code

table1(~ factor(case_label) + factor(date_at_enrollment) + factor(et) + factor(lact) + 
         conception_to_enrollment_d + calving_to_enrollment_d + conception_to_preg_dx_1 | factor(exposure_label), data=coxiella_filtered)

	Coxiella seronegative (N=203)	Coxiella seropositive (N=99)	Overall (N=302)
factor(case_label)
Case	77 (37.9%)	38 (38.4%)	115 (38.1%)
Control	126 (62.1%)	61 (61.6%)	187 (61.9%)
factor(date_at_enrollment)
2022-07-04	17 (8.4%)	8 (8.1%)	25 (8.3%)
2022-07-19	28 (13.8%)	16 (16.2%)	44 (14.6%)
2022-08-02	20 (9.9%)	18 (18.2%)	38 (12.6%)
2022-10-13	34 (16.7%)	13 (13.1%)	47 (15.6%)
2023-08-22	37 (18.2%)	15 (15.2%)	52 (17.2%)
2023-08-23	12 (5.9%)	3 (3.0%)	15 (5.0%)
2023-11-20	55 (27.1%)	26 (26.3%)	81 (26.8%)
factor(et)
0	100 (49.3%)	45 (45.5%)	145 (48.0%)
1	103 (50.7%)	54 (54.5%)	157 (52.0%)
factor(lact)
1	40 (19.7%)	15 (15.2%)	55 (18.2%)
2	94 (46.3%)	50 (50.5%)	144 (47.7%)
3	30 (14.8%)	15 (15.2%)	45 (14.9%)
4	19 (9.4%)	9 (9.1%)	28 (9.3%)
5	12 (5.9%)	7 (7.1%)	19 (6.3%)
6	6 (3.0%)	3 (3.0%)	9 (3.0%)
7	1 (0.5%)	0 (0%)	1 (0.3%)
Missing	1 (0.5%)	0 (0%)	1 (0.3%)
conception_to_enrollment_d
Mean (SD)	60.6 (16.1)	60.7 (13.0)	60.6 (15.1)
Median [Min, Max]	57.0 [36.0, 111]	60.0 [38.0, 95.0]	58.0 [36.0, 111]
calving_to_enrollment_d
Mean (SD)	157 (39.7)	154 (39.6)	156 (39.6)
Median [Min, Max]	148 [106, 315]	144 [111, 315]	147 [106, 315]
Missing	1 (0.5%)	0 (0%)	1 (0.3%)
conception_to_preg_dx_1
Mean (SD)	41.4 (16.0)	39.1 (13.8)	40.6 (15.3)
Median [Min, Max]	34.0 [30.0, 82.0]	34.0 [30.0, 75.0]	34.0 [30.0, 82.0]

Analysis using tabular methods

All data (including excluded cows)

Code

data_tbl <- coxiella %>%
  mutate(
    coxiella_test_status = factor(coxiella_test_status,levels = c(1,0)),
    case = factor(case,levels = c(1,0))) %>%
  group_by(coxiella_test_status,case) %>% 
  summarise(n = n())

data_out <- epi.2by2(dat = data_tbl, method = "case.control", digits = 2, 
   conf.level = 0.95, units = 100, interpret = FALSE, outcome = "as.columns")

data_out

             Outcome +    Outcome -      Total                       Odds
Exposed +           40           80        120        0.50 (0.33 to 0.71)
Exposed -           86          171        257        0.50 (0.38 to 0.65)
Total              126          251        377        0.50 (0.40 to 0.62)

Point estimates and 95% CIs:
-------------------------------------------------------------------
Exposure odds ratio                            0.99 (0.63, 1.57)
Attrib fraction (est) in the exposed (%)      -0.58 (-64.20, 37.89)
Attrib fraction (est) in the population (%)   -0.19 (-15.94, 13.43)
-------------------------------------------------------------------
Uncorrected chi2 test that OR = 1: chi2(1) = 0.001 Pr>chi2 = 0.980
Fisher exact test that OR = 1: Pr>chi2 = 1.000
 Wald confidence limits
 CI: confidence interval

Code

tabular_all_data <- tibble(
    association = "All enrolled cows - Tabular",
  odds_ratio = paste0(
    round(data_out$massoc.summary$est[1],2)," (95% CI: ",
    round(data_out$massoc.summary$lower[1],2)," to ",
    round(data_out$massoc.summary$upper[1],2), ")"
  ),
  p_value = round(data_out$massoc.detail$chi2.strata.uncor$p.value.2s,3)
)

All data (final data)

Code

data_tbl <- coxiella_filtered %>%
  mutate(
    coxiella_test_status = factor(coxiella_test_status,levels = c(1,0)),
    case = factor(case,levels = c(1,0))) %>%
  group_by(coxiella_test_status,case) %>% 
  summarise(n = n())

data_out <- epi.2by2(dat = data_tbl, method = "case.control", digits = 2, 
   conf.level = 0.95, units = 100, interpret = FALSE, outcome = "as.columns")

data_out

             Outcome +    Outcome -      Total                       Odds
Exposed +           38           61         99        0.62 (0.41 to 0.94)
Exposed -           77          126        203        0.61 (0.46 to 0.81)
Total              115          187        302        0.61 (0.49 to 0.78)

Point estimates and 95% CIs:
-------------------------------------------------------------------
Exposure odds ratio                            1.02 (0.62, 1.67)
Attrib fraction (est) in the exposed (%)      1.90 (-66.33, 41.78)
Attrib fraction (est) in the population (%)   0.63 (-16.91, 15.54)
-------------------------------------------------------------------
Uncorrected chi2 test that OR = 1: chi2(1) = 0.006 Pr>chi2 = 0.939
Fisher exact test that OR = 1: Pr>chi2 = 1.000
 Wald confidence limits
 CI: confidence interval

Code

tabular_filtered <- tibble(
    association = "Final dataset (removed ineligible cows) - Tabular",
  odds_ratio = paste0(
    round(data_out$massoc.summary$est[1],2)," (95% CI: ",
    round(data_out$massoc.summary$lower[1],2)," to ",
    round(data_out$massoc.summary$upper[1],2), ")"
  ),
  p_value = round(data_out$massoc.detail$chi2.strata.uncor$p.value.2s,3)
)

AI pregnancies only

Code

data_tbl <- coxiella_filtered %>% filter(et == 0) %>%
  mutate(
    coxiella_test_status = factor(coxiella_test_status,levels = c(1,0)),
    case = factor(case,levels = c(1,0))) %>%
  group_by(coxiella_test_status,case) %>% 
  summarise(n = n())

data_out <- epi.2by2(dat = data_tbl, method = "case.control", digits = 2, 
   conf.level = 0.95, units = 100, interpret = FALSE, outcome = "as.columns")

data_out

             Outcome +    Outcome -      Total                       Odds
Exposed +           14           31         45        0.45 (0.22 to 0.80)
Exposed -           30           70        100        0.43 (0.27 to 0.64)
Total               44          101        145        0.44 (0.29 to 0.61)

Point estimates and 95% CIs:
-------------------------------------------------------------------
Exposure odds ratio                            1.05 (0.49, 2.26)
Attrib fraction (est) in the exposed (%)      5.07 (-121.83, 58.20)
Attrib fraction (est) in the population (%)   1.62 (-25.06, 22.61)
-------------------------------------------------------------------
Uncorrected chi2 test that OR = 1: chi2(1) = 0.018 Pr>chi2 = 0.893
Fisher exact test that OR = 1: Pr>chi2 = 1.000
 Wald confidence limits
 CI: confidence interval

Code

tabular_ai <- tibble(
  association = "AI Pregnancies only - Tabular",
  odds_ratio = paste0(
    round(data_out$massoc.summary$est[1],2)," (95% CI: ",
    round(data_out$massoc.summary$lower[1],2)," to ",
    round(data_out$massoc.summary$upper[1],2), ")"
  ),
  p_value = round(data_out$massoc.detail$chi2.strata.uncor$p.value.2s,3)
)

ET pregnancies only

Code

data_tbl <- coxiella_filtered %>% filter(et == 1) %>%
  mutate(
    coxiella_test_status = factor(coxiella_test_status,levels = c(1,0)),
    case = factor(case,levels = c(1,0))) %>%
  group_by(coxiella_test_status,case) %>% 
  summarise(n = n())

data_out <- epi.2by2(dat = data_tbl, method = "case.control", digits = 2, 
   conf.level = 0.95, units = 100, interpret = FALSE, outcome = "as.columns")

data_out

             Outcome +    Outcome -      Total                       Odds
Exposed +           24           30         54        0.80 (0.46 to 1.35)
Exposed -           47           56        103        0.84 (0.56 to 1.24)
Total               71           86        157        0.83 (0.60 to 1.12)

Point estimates and 95% CIs:
-------------------------------------------------------------------
Exposure odds ratio                            0.95 (0.49, 1.85)
Attrib fraction (est) in the exposed (%)      -4.88 (-115.12, 48.58)
Attrib fraction (est) in the population (%)   -1.66 (-27.57, 18.99)
-------------------------------------------------------------------
Uncorrected chi2 test that OR = 1: chi2(1) = 0.020 Pr>chi2 = 0.887
Fisher exact test that OR = 1: Pr>chi2 = 1.000
 Wald confidence limits
 CI: confidence interval

Code

tabular_et <- tibble(
  association = "ET Pregnancies only - Tabular",
  odds_ratio = paste0(
    round(data_out$massoc.summary$est[1],2)," (95% CI: ",
    round(data_out$massoc.summary$lower[1],2)," to ",
    round(data_out$massoc.summary$upper[1],2), ")"
  ),
  p_value = round(data_out$massoc.detail$chi2.strata.uncor$p.value.2s,3)
)

Summary of tabular methods

Code

rbind(tabular_all_data,
      tabular_filtered,
      tabular_ai,
      tabular_et) %>% select(association,odds_ratio,p_value)

association	odds_ratio	p_value
All enrolled cows - Tabular	0.99 (95% CI: 0.63 to 1.57)	0.980
Final dataset (removed ineligible cows) - Tabular	1.02 (95% CI: 0.62 to 1.67)	0.939
AI Pregnancies only - Tabular	1.05 (95% CI: 0.49 to 2.26)	0.893
ET Pregnancies only - Tabular	0.95 (95% CI: 0.49 to 1.85)	0.887

Models

Univariable conditional logistic regression. No covariates.

Code

cond_logit <- clogit(case ~ coxiella_test_status + strata(linked_sample),
                     data = coxiella_filtered)

summary(cond_logit)

Call:
coxph(formula = Surv(rep(1, 302L), case) ~ coxiella_test_status + 
    strata(linked_sample), data = coxiella_filtered, method = "exact")

  n= 302, number of events= 115 

                        coef exp(coef) se(coef)     z Pr(>|z|)
coxiella_test_status 0.03394   1.03453  0.24106 0.141    0.888

                     exp(coef) exp(-coef) lower .95 upper .95
coxiella_test_status     1.035     0.9666     0.645     1.659

Concordance= 0.508  (se = 0.037 )
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

Code

mod1 <- tidy(cond_logit, exp = T, conf.int = T) %>% 
  mutate(odds_ratio = paste0(round(estimate,2)," (95% CI: ",round(conf.low,2)," to ",round(conf.high,2),")"),
         model = "Univariable Model") %>%
  filter(term == "coxiella_test_status") %>%
  select(model,term,odds_ratio,p.value)

Multivariable conditional logistic regression. ET status as covariate.

Code

cond_logit <- clogit(case ~ coxiella_test_status + et + strata(linked_sample),
                     data = coxiella_filtered)

summary(cond_logit)

Call:
coxph(formula = Surv(rep(1, 302L), case) ~ coxiella_test_status + 
    et + strata(linked_sample), data = coxiella_filtered, method = "exact")

  n= 302, number of events= 115 

                        coef exp(coef) se(coef)      z Pr(>|z|)    
coxiella_test_status -0.1216    0.8855   0.2548 -0.477 0.633194    
et                    1.3617    3.9028   0.4066  3.349 0.000812 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                     exp(coef) exp(-coef) lower .95 upper .95
coxiella_test_status    0.8855     1.1293    0.5374     1.459
et                      3.9028     0.2562    1.7589     8.660

Concordance= 0.583  (se = 0.042 )
Likelihood ratio test= 13.95  on 2 df,   p=9e-04
Wald test            = 11.22  on 2 df,   p=0.004
Score (logrank) test = 12.69  on 2 df,   p=0.002

Code

mod2 <- tidy(cond_logit, exp = T, conf.int = T) %>% 
  mutate(odds_ratio = paste0(round(estimate,2)," (95% CI: ",round(conf.low,2)," to ",round(conf.high,2),")"),
         model = "Adjusted for ET status") %>%
  filter(term == "coxiella_test_status") %>%
  select(model,term,odds_ratio,p.value)

Multivariable conditional logistic regression. Covariates = Et status, parity, and calving to enrollment, conception to enrollment

Code

cond_logit <- clogit(case ~ coxiella_test_status + et + conception_to_enrollment_d + lact + strata(linked_sample),data = coxiella_filtered)


summary(cond_logit)

Call:
coxph(formula = Surv(rep(1, 302L), case) ~ coxiella_test_status + 
    et + conception_to_enrollment_d + lact + strata(linked_sample), 
    data = coxiella_filtered, method = "exact")

  n= 301, number of events= 115 
   (1 observation deleted due to missingness)

                               coef exp(coef) se(coef)      z Pr(>|z|)    
coxiella_test_status       -0.22052   0.80210  0.26281 -0.839 0.401412    
et                          1.36530   3.91689  0.41212  3.313 0.000923 ***
conception_to_enrollment_d  0.08439   1.08805  0.05777  1.461 0.144075    
lact                        1.67076   5.31619  1.03064  1.621 0.104998    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                           exp(coef) exp(-coef) lower .95 upper .95
coxiella_test_status          0.8021     1.2467    0.4792     1.343
et                            3.9169     0.2553    1.7464     8.785
conception_to_enrollment_d    1.0881     0.9191    0.9716     1.218
lact                          5.3162     0.1881    0.7052    40.076

Concordance= 0.583  (se = 0.047 )
Likelihood ratio test= 20.45  on 4 df,   p=4e-04
Wald test            = 15.01  on 4 df,   p=0.005
Score (logrank) test = 17.78  on 4 df,   p=0.001

Code

mod3 <- tidy(cond_logit, exp = T, conf.int = T) %>% 
  mutate(odds_ratio = paste0(round(estimate,2)," (95% CI: ",round(conf.low,2)," to ",round(conf.high,2),")"),
         model = "Adjusted for ET status and matching variables") %>%
  filter(term == "coxiella_test_status") %>%
  select(model,term,odds_ratio,p.value)

Interaction term between ET and coxiella status (ELISA)

Code

cond_logit <- clogit(case ~ coxiella_test_status*et + strata(linked_sample),
                     data = coxiella_filtered)

summary(cond_logit)

Call:
coxph(formula = Surv(rep(1, 302L), case) ~ coxiella_test_status * 
    et + strata(linked_sample), data = coxiella_filtered, method = "exact")

  n= 302, number of events= 115 

                            coef exp(coef) se(coef)      z Pr(>|z|)   
coxiella_test_status    -0.15394   0.85732  0.36415 -0.423  0.67248   
et                       1.34227   3.82772  0.43509  3.085  0.00204 **
coxiella_test_status:et  0.06423   1.06633  0.51468  0.125  0.90069   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                        exp(coef) exp(-coef) lower .95 upper .95
coxiella_test_status       0.8573     1.1664    0.4199     1.750
et                         3.8277     0.2613    1.6315     8.980
coxiella_test_status:et    1.0663     0.9378    0.3889     2.924

Concordance= 0.583  (se = 0.042 )
Likelihood ratio test= 13.97  on 3 df,   p=0.003
Wald test            = 11.23  on 3 df,   p=0.01
Score (logrank) test = 12.7  on 3 df,   p=0.005

Code

mod4 <- tidy(cond_logit, exp = T, conf.int = T) %>% 
  mutate(odds_ratio = paste0(round(estimate,2)," (95% CI: ",round(conf.low,2)," to ",round(conf.high,2),")"),
         model = "Interaction Model") %>%
  select(model,term,odds_ratio,p.value)

Stratified model (equivalent to interaction model)

AI pregnacies only

Code

cond_logit <- clogit(case ~ coxiella_test_status + strata(linked_sample),
                     data = coxiella_filtered %>% filter(et == 0))

summary(cond_logit)

Call:
coxph(formula = Surv(rep(1, 145L), case) ~ coxiella_test_status + 
    strata(linked_sample), data = coxiella_filtered %>% filter(et == 
    0), method = "exact")

  n= 145, number of events= 44 

                        coef exp(coef) se(coef)      z Pr(>|z|)
coxiella_test_status -0.1693    0.8442   0.3836 -0.441    0.659

                     exp(coef) exp(-coef) lower .95 upper .95
coxiella_test_status    0.8442      1.185    0.3981     1.791

Concordance= 0.508  (se = 0.073 )
Likelihood ratio test= 0.2  on 1 df,   p=0.7
Wald test            = 0.19  on 1 df,   p=0.7
Score (logrank) test = 0.2  on 1 df,   p=0.7

Code

mod5 <- tidy(cond_logit, exp = T, conf.int = T) %>% 
  mutate(odds_ratio = paste0(round(estimate,2)," (95% CI: ",round(conf.low,2)," to ",round(conf.high,2),")"),
         model = "AI Pregnancies only") %>%
  select(model,term,odds_ratio,p.value)

ET pregnacies only

Code

cond_logit <- clogit(case ~ coxiella_test_status + strata(linked_sample),
                     data = coxiella_filtered %>% filter(et == 1))

summary(cond_logit)

Call:
coxph(formula = Surv(rep(1, 157L), case) ~ coxiella_test_status + 
    strata(linked_sample), data = coxiella_filtered %>% filter(et == 
    1), method = "exact")

  n= 157, number of events= 71 

                        coef exp(coef) se(coef)     z Pr(>|z|)
coxiella_test_status 0.05122   1.05255  0.39178 0.131    0.896

                     exp(coef) exp(-coef) lower .95 upper .95
coxiella_test_status     1.053     0.9501    0.4884     2.268

Concordance= 0.5  (se = 0.05 )
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

Code

mod6 <- tidy(cond_logit, exp = T, conf.int = T) %>% 
  mutate(odds_ratio = paste0(round(estimate,2)," (95% CI: ",round(conf.low,2)," to ",round(conf.high,2),")"),
         model = "ET Pregnancies only") %>%
  select(model,term,odds_ratio,p.value)

Summary of models

Code

rbind(mod1,mod2,mod3,mod5,mod6) %>% select(model,odds_ratio,p.value)

model	odds_ratio	p.value
Univariable Model	1.03 (95% CI: 0.64 to 1.66)	0.8880172
Adjusted for ET status	0.89 (95% CI: 0.54 to 1.46)	0.6331938
Adjusted for ET status and matching variables	0.8 (95% CI: 0.48 to 1.34)	0.4014121
AI Pregnancies only	0.84 (95% CI: 0.4 to 1.79)	0.6589281
ET Pregnancies only	1.05 (95% CI: 0.49 to 2.27)	0.8959880