1 Exploratory analysis

1.1 Samples by postpartum sample, DIM and farm

# Investigating association between number of samples, postpartum sample and farm

samples_farm <- ggplot(df, aes(x = SampleOrder.overall, y = number,fill = Site)) +
  geom_bar(position = "stack", stat = "identity") +
  labs(fill = "Sample", y = "Count", x = "Postpartum sample") +
  scale_fill_viridis(discrete = TRUE)

ggarrange(samples_farm,common.legend = TRUE)

ggsave(plot = last_plot(),"./figures/exploratory_sample_farm.png",width = 20, height = 10, units = "cm") 

# Investigating association between days in milk and postpartum sample

sample_dim <- ggplot(df, aes(x = as.numeric(DIM), fill = as.factor(SampleOrder.overall))) +
  geom_histogram(binwidth = .5, alpha = .5, position = "identity") +
  labs(fill = "Postpartum sample", y = "Number of samples", x = "DIM") + 
  scale_fill_viridis(discrete = TRUE)

ggarrange(sample_dim,common.legend = TRUE)

ggsave(plot = last_plot(),"./figures/exploratory_dim_sample.png",width = 20, height = 10, units = "cm") 

1.2 Prevalence by postpartum sample

# Prevalence of all microorganisms in dataset

table1(~`Staphylococcus aureus` + nas + `Staphylococcus chromogenes` + nas.non.chrom + `Staphylococcus epidermidis` + `Staphylococcus haemolyticus` + `Staphylococcus haemolyticus` + `Staphylococcus hominis` + `Staphylococcus sciuri` + `Staphylococcus xylosus/saprophyticus` + `Staphylococcus sp.` + strep + strep.genus + `Streptococcus uberis` + `Streptococcus dysgalactiae` + `Streptococcus sp.` + strep.like + `Aerococcus sp.` + `Aerococcus viridans` + `Enterococcus casseliflavus` + `Enterococcus hirae` + `Enterococcus mundtii` + `Enterococcus sp.` + `Enterococcus faecalis` + `Enterococcus saccarolyticus` + `Lactococcus garvieae` + `Lactococcus lactis` + `Lactococcus sp.` + `Micrococcus sp.` + gram.neg +`Escherichia coli` + `Escherichia sp.` + `Klebsiella oxytoca` + `Pseudomonas sp.` + `Citrobacter freundii` + `Enterobacter cloacae` + `Enterobacter sp.` + `Gram Negative Organism` + `Pantoea sp.` + `Serratia sp.` +  other + `Arthrobacter gandavensis` + `Arthrobacter sp.` + `Bacillus sp.` + `Corynebacterium sp.` + `Gram Positive Cocci` + `Gram Positive Rod` + `Helcococcus ovis` + `Trueperella pyogenes` + Yeast + `No Growth` | SampleOrder.overall  , data = df)
1
(N=424)
2
(N=404)
3
(N=373)
4
(N=320)
5
(N=152)
Overall
(N=1673)
Staphylococcus aureus
0 362 (85.4%) 349 (86.4%) 321 (86.1%) 277 (86.6%) 133 (87.5%) 1442 (86.2%)
1 62 (14.6%) 55 (13.6%) 52 (13.9%) 43 (13.4%) 19 (12.5%) 231 (13.8%)
nas
0 292 (68.9%) 277 (68.6%) 257 (68.9%) 241 (75.3%) 108 (71.1%) 1175 (70.2%)
1 132 (31.1%) 127 (31.4%) 116 (31.1%) 79 (24.7%) 44 (28.9%) 498 (29.8%)
Staphylococcus chromogenes
0 316 (74.5%) 313 (77.5%) 281 (75.3%) 257 (80.3%) 119 (78.3%) 1286 (76.9%)
1 108 (25.5%) 91 (22.5%) 92 (24.7%) 63 (19.7%) 33 (21.7%) 387 (23.1%)
nas.non.chrom
0 394 (92.9%) 362 (89.6%) 342 (91.7%) 299 (93.4%) 139 (91.4%) 1536 (91.8%)
1 30 (7.1%) 42 (10.4%) 31 (8.3%) 21 (6.6%) 13 (8.6%) 137 (8.2%)
Staphylococcus epidermidis
0 424 (100%) 404 (100%) 372 (99.7%) 320 (100%) 152 (100%) 1672 (99.9%)
1 0 (0%) 0 (0%) 1 (0.3%) 0 (0%) 0 (0%) 1 (0.1%)
Staphylococcus haemolyticus
0 424 (100%) 401 (99.3%) 371 (99.5%) 317 (99.1%) 151 (99.3%) 1664 (99.5%)
1 0 (0%) 3 (0.7%) 2 (0.5%) 3 (0.9%) 1 (0.7%) 9 (0.5%)
Staphylococcus hominis
0 424 (100%) 403 (99.8%) 371 (99.5%) 320 (100%) 152 (100%) 1670 (99.8%)
1 0 (0%) 1 (0.2%) 2 (0.5%) 0 (0%) 0 (0%) 3 (0.2%)
Staphylococcus sciuri
0 423 (99.8%) 403 (99.8%) 372 (99.7%) 319 (99.7%) 152 (100%) 1669 (99.8%)
1 1 (0.2%) 1 (0.2%) 1 (0.3%) 1 (0.3%) 0 (0%) 4 (0.2%)
Staphylococcus xylosus/saprophyticus
0 421 (99.3%) 400 (99.0%) 371 (99.5%) 319 (99.7%) 152 (100%) 1663 (99.4%)
1 3 (0.7%) 4 (1.0%) 2 (0.5%) 1 (0.3%) 0 (0%) 10 (0.6%)
Staphylococcus sp.
0 398 (93.9%) 369 (91.3%) 349 (93.6%) 304 (95.0%) 140 (92.1%) 1560 (93.2%)
1 26 (6.1%) 35 (8.7%) 24 (6.4%) 16 (5.0%) 12 (7.9%) 113 (6.8%)
strep
0 350 (82.5%) 346 (85.6%) 320 (85.8%) 283 (88.4%) 138 (90.8%) 1437 (85.9%)
1 74 (17.5%) 58 (14.4%) 53 (14.2%) 37 (11.6%) 14 (9.2%) 236 (14.1%)
strep.genus
0 375 (88.4%) 368 (91.1%) 342 (91.7%) 296 (92.5%) 144 (94.7%) 1525 (91.2%)
1 49 (11.6%) 36 (8.9%) 31 (8.3%) 24 (7.5%) 8 (5.3%) 148 (8.8%)
Streptococcus uberis
0 418 (98.6%) 402 (99.5%) 369 (98.9%) 317 (99.1%) 151 (99.3%) 1657 (99.0%)
1 6 (1.4%) 2 (0.5%) 4 (1.1%) 3 (0.9%) 1 (0.7%) 16 (1.0%)
Streptococcus dysgalactiae
0 391 (92.2%) 387 (95.8%) 356 (95.4%) 303 (94.7%) 147 (96.7%) 1584 (94.7%)
1 33 (7.8%) 17 (4.2%) 17 (4.6%) 17 (5.3%) 5 (3.3%) 89 (5.3%)
Streptococcus sp.
0 414 (97.6%) 387 (95.8%) 363 (97.3%) 315 (98.4%) 150 (98.7%) 1629 (97.4%)
1 10 (2.4%) 17 (4.2%) 10 (2.7%) 5 (1.6%) 2 (1.3%) 44 (2.6%)
strep.like
0 399 (94.1%) 381 (94.3%) 350 (93.8%) 307 (95.9%) 146 (96.1%) 1583 (94.6%)
1 25 (5.9%) 23 (5.7%) 23 (6.2%) 13 (4.1%) 6 (3.9%) 90 (5.4%)
Aerococcus sp.
0 424 (100%) 401 (99.3%) 372 (99.7%) 320 (100%) 152 (100%) 1669 (99.8%)
1 0 (0%) 3 (0.7%) 1 (0.3%) 0 (0%) 0 (0%) 4 (0.2%)
Aerococcus viridans
0 414 (97.6%) 399 (98.8%) 367 (98.4%) 316 (98.8%) 152 (100%) 1648 (98.5%)
1 10 (2.4%) 5 (1.2%) 6 (1.6%) 4 (1.3%) 0 (0%) 25 (1.5%)
Enterococcus casseliflavus
0 418 (98.6%) 401 (99.3%) 371 (99.5%) 318 (99.4%) 151 (99.3%) 1659 (99.2%)
1 6 (1.4%) 3 (0.7%) 2 (0.5%) 2 (0.6%) 1 (0.7%) 14 (0.8%)
Enterococcus hirae
0 421 (99.3%) 402 (99.5%) 371 (99.5%) 319 (99.7%) 151 (99.3%) 1664 (99.5%)
1 3 (0.7%) 2 (0.5%) 2 (0.5%) 1 (0.3%) 1 (0.7%) 9 (0.5%)
Enterococcus mundtii
0 422 (99.5%) 401 (99.3%) 371 (99.5%) 317 (99.1%) 152 (100%) 1663 (99.4%)
1 2 (0.5%) 3 (0.7%) 2 (0.5%) 3 (0.9%) 0 (0%) 10 (0.6%)
Enterococcus sp.
0 419 (98.8%) 398 (98.5%) 366 (98.1%) 319 (99.7%) 149 (98.0%) 1651 (98.7%)
1 5 (1.2%) 6 (1.5%) 7 (1.9%) 1 (0.3%) 3 (2.0%) 22 (1.3%)
Enterococcus faecalis
0 424 (100%) 403 (99.8%) 373 (100%) 320 (100%) 152 (100%) 1672 (99.9%)
1 0 (0%) 1 (0.2%) 0 (0%) 0 (0%) 0 (0%) 1 (0.1%)
Enterococcus saccarolyticus
0 424 (100%) 404 (100%) 373 (100%) 320 (100%) 151 (99.3%) 1672 (99.9%)
1 0 (0%) 0 (0%) 0 (0%) 0 (0%) 1 (0.7%) 1 (0.1%)
Lactococcus garvieae
0 423 (99.8%) 403 (99.8%) 373 (100%) 319 (99.7%) 152 (100%) 1670 (99.8%)
1 1 (0.2%) 1 (0.2%) 0 (0%) 1 (0.3%) 0 (0%) 3 (0.2%)
Lactococcus lactis
0 424 (100%) 404 (100%) 371 (99.5%) 319 (99.7%) 152 (100%) 1670 (99.8%)
1 0 (0%) 0 (0%) 2 (0.5%) 1 (0.3%) 0 (0%) 3 (0.2%)
Lactococcus sp.
0 424 (100%) 403 (99.8%) 373 (100%) 320 (100%) 152 (100%) 1672 (99.9%)
1 0 (0%) 1 (0.2%) 0 (0%) 0 (0%) 0 (0%) 1 (0.1%)
Micrococcus sp.
0 423 (99.8%) 403 (99.8%) 372 (99.7%) 320 (100%) 152 (100%) 1670 (99.8%)
1 1 (0.2%) 1 (0.2%) 1 (0.3%) 0 (0%) 0 (0%) 3 (0.2%)
gram.neg
0 412 (97.2%) 394 (97.5%) 363 (97.3%) 315 (98.4%) 150 (98.7%) 1634 (97.7%)
1 12 (2.8%) 10 (2.5%) 10 (2.7%) 5 (1.6%) 2 (1.3%) 39 (2.3%)
Escherichia coli
0 418 (98.6%) 399 (98.8%) 371 (99.5%) 320 (100%) 152 (100%) 1660 (99.2%)
1 6 (1.4%) 5 (1.2%) 2 (0.5%) 0 (0%) 0 (0%) 13 (0.8%)
Escherichia sp.
0 422 (99.5%) 404 (100%) 372 (99.7%) 320 (100%) 152 (100%) 1670 (99.8%)
1 2 (0.5%) 0 (0%) 1 (0.3%) 0 (0%) 0 (0%) 3 (0.2%)
Klebsiella oxytoca
0 424 (100%) 404 (100%) 372 (99.7%) 318 (99.4%) 152 (100%) 1670 (99.8%)
1 0 (0%) 0 (0%) 1 (0.3%) 2 (0.6%) 0 (0%) 3 (0.2%)
Pseudomonas sp.
0 423 (99.8%) 401 (99.3%) 371 (99.5%) 318 (99.4%) 151 (99.3%) 1664 (99.5%)
1 1 (0.2%) 3 (0.7%) 2 (0.5%) 2 (0.6%) 1 (0.7%) 9 (0.5%)
Citrobacter freundii
0 424 (100%) 404 (100%) 373 (100%) 319 (99.7%) 152 (100%) 1672 (99.9%)
1 0 (0%) 0 (0%) 0 (0%) 1 (0.3%) 0 (0%) 1 (0.1%)
Enterobacter cloacae
0 424 (100%) 404 (100%) 372 (99.7%) 320 (100%) 152 (100%) 1672 (99.9%)
1 0 (0%) 0 (0%) 1 (0.3%) 0 (0%) 0 (0%) 1 (0.1%)
Enterobacter sp.
0 423 (99.8%) 404 (100%) 373 (100%) 320 (100%) 152 (100%) 1672 (99.9%)
1 1 (0.2%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 1 (0.1%)
Gram Negative Organism
0 422 (99.5%) 404 (100%) 372 (99.7%) 319 (99.7%) 151 (99.3%) 1668 (99.7%)
1 2 (0.5%) 0 (0%) 1 (0.3%) 1 (0.3%) 1 (0.7%) 5 (0.3%)
Pantoea sp.
0 423 (99.8%) 403 (99.8%) 372 (99.7%) 320 (100%) 152 (100%) 1670 (99.8%)
1 1 (0.2%) 1 (0.2%) 1 (0.3%) 0 (0%) 0 (0%) 3 (0.2%)
Serratia sp.
0 424 (100%) 403 (99.8%) 372 (99.7%) 320 (100%) 151 (99.3%) 1670 (99.8%)
1 0 (0%) 1 (0.2%) 1 (0.3%) 0 (0%) 1 (0.7%) 3 (0.2%)
other
0 384 (90.6%) 358 (88.6%) 319 (85.5%) 277 (86.6%) 130 (85.5%) 1468 (87.7%)
1 40 (9.4%) 46 (11.4%) 54 (14.5%) 43 (13.4%) 22 (14.5%) 205 (12.3%)
Arthrobacter gandavensis
0 423 (99.8%) 403 (99.8%) 373 (100%) 320 (100%) 152 (100%) 1671 (99.9%)
1 1 (0.2%) 1 (0.2%) 0 (0%) 0 (0%) 0 (0%) 2 (0.1%)
Arthrobacter sp.
0 424 (100%) 403 (99.8%) 373 (100%) 320 (100%) 152 (100%) 1672 (99.9%)
1 0 (0%) 1 (0.2%) 0 (0%) 0 (0%) 0 (0%) 1 (0.1%)
Bacillus sp.
0 406 (95.8%) 376 (93.1%) 346 (92.8%) 295 (92.2%) 142 (93.4%) 1565 (93.5%)
1 18 (4.2%) 28 (6.9%) 27 (7.2%) 25 (7.8%) 10 (6.6%) 108 (6.5%)
Corynebacterium sp.
0 412 (97.2%) 389 (96.3%) 353 (94.6%) 303 (94.7%) 145 (95.4%) 1602 (95.8%)
1 12 (2.8%) 15 (3.7%) 20 (5.4%) 17 (5.3%) 7 (4.6%) 71 (4.2%)
Gram Positive Cocci
0 423 (99.8%) 400 (99.0%) 371 (99.5%) 316 (98.8%) 149 (98.0%) 1659 (99.2%)
1 1 (0.2%) 4 (1.0%) 2 (0.5%) 4 (1.3%) 3 (2.0%) 14 (0.8%)
Gram Positive Rod
0 418 (98.6%) 403 (99.8%) 366 (98.1%) 316 (98.8%) 150 (98.7%) 1653 (98.8%)
1 6 (1.4%) 1 (0.2%) 7 (1.9%) 4 (1.3%) 2 (1.3%) 20 (1.2%)
Helcococcus ovis
0 424 (100%) 404 (100%) 372 (99.7%) 320 (100%) 152 (100%) 1672 (99.9%)
1 0 (0%) 0 (0%) 1 (0.3%) 0 (0%) 0 (0%) 1 (0.1%)
Trueperella pyogenes
0 423 (99.8%) 403 (99.8%) 372 (99.7%) 320 (100%) 152 (100%) 1670 (99.8%)
1 1 (0.2%) 1 (0.2%) 1 (0.3%) 0 (0%) 0 (0%) 3 (0.2%)
Yeast
0 423 (99.8%) 404 (100%) 372 (99.7%) 320 (100%) 152 (100%) 1671 (99.9%)
1 1 (0.2%) 0 (0%) 1 (0.3%) 0 (0%) 0 (0%) 2 (0.1%)
No Growth
0 251 (59.2%) 228 (56.4%) 227 (60.9%) 172 (53.8%) 85 (55.9%) 963 (57.6%)
1 173 (40.8%) 176 (43.6%) 146 (39.1%) 148 (46.3%) 67 (44.1%) 710 (42.4%)

1.3 Period prevalence by farm

# ---- Period prevalence by farm (table) ------

table1(~IMI.cow + staph.aureus.IMI.cow + nas.IMI.cow + sch.IMI.cow + nas.non.chrom.IMI.cow + staph.epidermidis.IMI.cow + staph.haemolyticus.IMI.cow + staph.hominis.IMI.cow + staph.sciuri.IMI.cow + staph.xylosus.saprophyticus.IMI.cow + staph.spp.IMI.cow + strep.IMI.cow + strep.genus.IMI.cow + strep.uberis.IMI.cow + strep.dys.IMI.cow + strep.sp.IMI.cow + strep.like.IMI.cow + aerococcus.spp.IMI.cow + aero.vir.IMI.cow + aero.spp.IMI.cow + enterococcus.spp.IMI.cow + entero.cass.IMI.cow + entero.hyr.IMI.cow + entero.mun.IMI.cow + entero.spp.IMI.cow + entero.fae.IMI.cow + entero.sac.IMI.cow + lactococcus.spp.IMI.cow + lact.garv.IMI.cow + lact.lact.IMI.cow + lact.spp.IMI.cow + micro.spp.IMI.cow + gram.neg.IMI.cow + escherichia.spp.IMI.cow + kleb.IMI.cow + pseud.IMI.cow + citr.IMI.cow + enterobacter.genus.IMI.cow + gram.neg.org.IMI.cow + pantoea.IMI.cow + serratia.IMI.cow + other.IMI.cow + arthro.gand.IMI.cow + arthro.spp.IMI.cow + bacillus.IMI.cow + coryne.IMI.cow + gram.pos.cocci.IMI.cow + gram.pos.rod.IMI.cow + helcoccus.IMI.cow + trueperella.IMI.cow + yeast.IMI.cow |Site , data = df.per.cow)
A
(N=123)
B
(N=117)
C
(N=111)
D
(N=21)
E
(N=52)
Overall
(N=424)
IMI.cow
0 47 (38.2%) 11 (9.4%) 6 (5.4%) 1 (4.8%) 0 (0%) 65 (15.3%)
1 76 (61.8%) 106 (90.6%) 105 (94.6%) 20 (95.2%) 52 (100%) 359 (84.7%)
staph.aureus.IMI.cow
0 108 (87.8%) 86 (73.5%) 101 (91.0%) 21 (100%) 28 (53.8%) 344 (81.1%)
1 15 (12.2%) 31 (26.5%) 10 (9.0%) 0 (0%) 24 (46.2%) 80 (18.9%)
nas.IMI.cow
0 72 (58.5%) 46 (39.3%) 65 (58.6%) 3 (14.3%) 17 (32.7%) 203 (47.9%)
1 51 (41.5%) 71 (60.7%) 46 (41.4%) 18 (85.7%) 35 (67.3%) 221 (52.1%)
sch.IMI.cow
0 81 (65.9%) 68 (58.1%) 76 (68.5%) 11 (52.4%) 27 (51.9%) 263 (62.0%)
1 42 (34.1%) 49 (41.9%) 35 (31.5%) 10 (47.6%) 25 (48.1%) 161 (38.0%)
nas.non.chrom.IMI.cow
0 109 (88.6%) 82 (70.1%) 98 (88.3%) 8 (38.1%) 31 (59.6%) 328 (77.4%)
1 14 (11.4%) 35 (29.9%) 13 (11.7%) 13 (61.9%) 21 (40.4%) 96 (22.6%)
staph.epidermidis.IMI.cow
0 122 (99.2%) 117 (100%) 111 (100%) 21 (100%) 52 (100%) 423 (99.8%)
1 1 (0.8%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 1 (0.2%)
staph.haemolyticus.IMI.cow
0 121 (98.4%) 115 (98.3%) 110 (99.1%) 19 (90.5%) 50 (96.2%) 415 (97.9%)
1 2 (1.6%) 2 (1.7%) 1 (0.9%) 2 (9.5%) 2 (3.8%) 9 (2.1%)
staph.hominis.IMI.cow
0 123 (100%) 117 (100%) 111 (100%) 19 (90.5%) 51 (98.1%) 421 (99.3%)
1 0 (0%) 0 (0%) 0 (0%) 2 (9.5%) 1 (1.9%) 3 (0.7%)
staph.sciuri.IMI.cow
0 123 (100%) 117 (100%) 111 (100%) 21 (100%) 48 (92.3%) 420 (99.1%)
1 0 (0%) 0 (0%) 0 (0%) 0 (0%) 4 (7.7%) 4 (0.9%)
staph.xylosus.saprophyticus.IMI.cow
0 122 (99.2%) 117 (100%) 110 (99.1%) 13 (61.9%) 52 (100%) 414 (97.6%)
1 1 (0.8%) 0 (0%) 1 (0.9%) 8 (38.1%) 0 (0%) 10 (2.4%)
staph.spp.IMI.cow
0 113 (91.9%) 83 (70.9%) 100 (90.1%) 17 (81.0%) 35 (67.3%) 348 (82.1%)
1 10 (8.1%) 34 (29.1%) 11 (9.9%) 4 (19.0%) 17 (32.7%) 76 (17.9%)
strep.IMI.cow
0 112 (91.1%) 82 (70.1%) 66 (59.5%) 9 (42.9%) 19 (36.5%) 288 (67.9%)
1 11 (8.9%) 35 (29.9%) 45 (40.5%) 12 (57.1%) 33 (63.5%) 136 (32.1%)
strep.genus.IMI.cow
0 118 (95.9%) 88 (75.2%) 102 (91.9%) 19 (90.5%) 26 (50.0%) 353 (83.3%)
1 5 (4.1%) 29 (24.8%) 9 (8.1%) 2 (9.5%) 26 (50.0%) 71 (16.7%)
strep.uberis.IMI.cow
0 123 (100%) 117 (100%) 111 (100%) 21 (100%) 43 (82.7%) 415 (97.9%)
1 0 (0%) 0 (0%) 0 (0%) 0 (0%) 9 (17.3%) 9 (2.1%)
strep.dys.IMI.cow
0 119 (96.7%) 92 (78.6%) 109 (98.2%) 19 (90.5%) 38 (73.1%) 377 (88.9%)
1 4 (3.3%) 25 (21.4%) 2 (1.8%) 2 (9.5%) 14 (26.9%) 47 (11.1%)
strep.sp.IMI.cow
0 121 (98.4%) 104 (88.9%) 103 (92.8%) 21 (100%) 43 (82.7%) 392 (92.5%)
1 2 (1.6%) 13 (11.1%) 8 (7.2%) 0 (0%) 9 (17.3%) 32 (7.5%)
strep.like.IMI.cow
0 116 (94.3%) 111 (94.9%) 71 (64.0%) 11 (52.4%) 43 (82.7%) 352 (83.0%)
1 7 (5.7%) 6 (5.1%) 40 (36.0%) 10 (47.6%) 9 (17.3%) 72 (17.0%)
aerococcus.spp.IMI.cow
0 120 (97.6%) 117 (100%) 108 (97.3%) 11 (52.4%) 45 (86.5%) 401 (94.6%)
1 3 (2.4%) 0 (0%) 3 (2.7%) 10 (47.6%) 7 (13.5%) 23 (5.4%)
aero.vir.IMI.cow
0 120 (97.6%) 117 (100%) 108 (97.3%) 11 (52.4%) 47 (90.4%) 403 (95.0%)
1 3 (2.4%) 0 (0%) 3 (2.7%) 10 (47.6%) 5 (9.6%) 21 (5.0%)
aero.spp.IMI.cow
0 123 (100%) 117 (100%) 111 (100%) 20 (95.2%) 49 (94.2%) 420 (99.1%)
1 0 (0%) 0 (0%) 0 (0%) 1 (4.8%) 3 (5.8%) 4 (0.9%)
enterococcus.spp.IMI.cow
0 120 (97.6%) 112 (95.7%) 77 (69.4%) 21 (100%) 51 (98.1%) 381 (89.9%)
1 3 (2.4%) 5 (4.3%) 34 (30.6%) 0 (0%) 1 (1.9%) 43 (10.1%)
entero.cass.IMI.cow
0 122 (99.2%) 116 (99.1%) 100 (90.1%) 21 (100%) 52 (100%) 411 (96.9%)
1 1 (0.8%) 1 (0.9%) 11 (9.9%) 0 (0%) 0 (0%) 13 (3.1%)
entero.hyr.IMI.cow
0 123 (100%) 117 (100%) 102 (91.9%) 21 (100%) 52 (100%) 415 (97.9%)
1 0 (0%) 0 (0%) 9 (8.1%) 0 (0%) 0 (0%) 9 (2.1%)
entero.mun.IMI.cow
0 123 (100%) 116 (99.1%) 102 (91.9%) 21 (100%) 52 (100%) 414 (97.6%)
1 0 (0%) 1 (0.9%) 9 (8.1%) 0 (0%) 0 (0%) 10 (2.4%)
entero.spp.IMI.cow
0 121 (98.4%) 114 (97.4%) 100 (90.1%) 21 (100%) 51 (98.1%) 407 (96.0%)
1 2 (1.6%) 3 (2.6%) 11 (9.9%) 0 (0%) 1 (1.9%) 17 (4.0%)
entero.fae.IMI.cow
0 123 (100%) 117 (100%) 110 (99.1%) 21 (100%) 52 (100%) 423 (99.8%)
1 0 (0%) 0 (0%) 1 (0.9%) 0 (0%) 0 (0%) 1 (0.2%)
entero.sac.IMI.cow
0 123 (100%) 117 (100%) 110 (99.1%) 21 (100%) 52 (100%) 423 (99.8%)
1 0 (0%) 0 (0%) 1 (0.9%) 0 (0%) 0 (0%) 1 (0.2%)
lactococcus.spp.IMI.cow
0 122 (99.2%) 117 (100%) 107 (96.4%) 20 (95.2%) 52 (100%) 418 (98.6%)
1 1 (0.8%) 0 (0%) 4 (3.6%) 1 (4.8%) 0 (0%) 6 (1.4%)
lact.garv.IMI.cow
0 122 (99.2%) 117 (100%) 109 (98.2%) 21 (100%) 52 (100%) 421 (99.3%)
1 1 (0.8%) 0 (0%) 2 (1.8%) 0 (0%) 0 (0%) 3 (0.7%)
lact.lact.IMI.cow
0 123 (100%) 117 (100%) 110 (99.1%) 20 (95.2%) 52 (100%) 422 (99.5%)
1 0 (0%) 0 (0%) 1 (0.9%) 1 (4.8%) 0 (0%) 2 (0.5%)
lact.spp.IMI.cow
0 123 (100%) 117 (100%) 110 (99.1%) 21 (100%) 52 (100%) 423 (99.8%)
1 0 (0%) 0 (0%) 1 (0.9%) 0 (0%) 0 (0%) 1 (0.2%)
micro.spp.IMI.cow
0 123 (100%) 116 (99.1%) 110 (99.1%) 21 (100%) 51 (98.1%) 421 (99.3%)
1 0 (0%) 1 (0.9%) 1 (0.9%) 0 (0%) 1 (1.9%) 3 (0.7%)
gram.neg.IMI.cow
0 119 (96.7%) 111 (94.9%) 88 (79.3%) 20 (95.2%) 49 (94.2%) 387 (91.3%)
1 4 (3.3%) 6 (5.1%) 23 (20.7%) 1 (4.8%) 3 (5.8%) 37 (8.7%)
escherichia.spp.IMI.cow
0 120 (97.6%) 113 (96.6%) 107 (96.4%) 20 (95.2%) 50 (96.2%) 410 (96.7%)
1 3 (2.4%) 4 (3.4%) 4 (3.6%) 1 (4.8%) 2 (3.8%) 14 (3.3%)
kleb.IMI.cow
0 123 (100%) 117 (100%) 108 (97.3%) 21 (100%) 52 (100%) 421 (99.3%)
1 0 (0%) 0 (0%) 3 (2.7%) 0 (0%) 0 (0%) 3 (0.7%)
pseud.IMI.cow
0 123 (100%) 117 (100%) 108 (97.3%) 21 (100%) 52 (100%) 421 (99.3%)
1 0 (0%) 0 (0%) 3 (2.7%) 0 (0%) 0 (0%) 3 (0.7%)
citr.IMI.cow
0 123 (100%) 117 (100%) 110 (99.1%) 21 (100%) 52 (100%) 423 (99.8%)
1 0 (0%) 0 (0%) 1 (0.9%) 0 (0%) 0 (0%) 1 (0.2%)
enterobacter.genus.IMI.cow
0 123 (100%) 117 (100%) 110 (99.1%) 21 (100%) 51 (98.1%) 422 (99.5%)
1 0 (0%) 0 (0%) 1 (0.9%) 0 (0%) 1 (1.9%) 2 (0.5%)
gram.neg.org.IMI.cow
0 122 (99.2%) 116 (99.1%) 108 (97.3%) 21 (100%) 52 (100%) 419 (98.8%)
1 1 (0.8%) 1 (0.9%) 3 (2.7%) 0 (0%) 0 (0%) 5 (1.2%)
pantoea.IMI.cow
0 123 (100%) 117 (100%) 108 (97.3%) 21 (100%) 52 (100%) 421 (99.3%)
1 0 (0%) 0 (0%) 3 (2.7%) 0 (0%) 0 (0%) 3 (0.7%)
serratia.IMI.cow
0 123 (100%) 117 (100%) 108 (97.3%) 21 (100%) 52 (100%) 421 (99.3%)
1 0 (0%) 0 (0%) 3 (2.7%) 0 (0%) 0 (0%) 3 (0.7%)
other.IMI.cow
0 108 (87.8%) 84 (71.8%) 32 (28.8%) 19 (90.5%) 43 (82.7%) 286 (67.5%)
1 15 (12.2%) 33 (28.2%) 79 (71.2%) 2 (9.5%) 9 (17.3%) 138 (32.5%)
arthro.gand.IMI.cow
0 123 (100%) 116 (99.1%) 110 (99.1%) 21 (100%) 52 (100%) 422 (99.5%)
1 0 (0%) 1 (0.9%) 1 (0.9%) 0 (0%) 0 (0%) 2 (0.5%)
arthro.spp.IMI.cow
0 123 (100%) 116 (99.1%) 111 (100%) 21 (100%) 52 (100%) 423 (99.8%)
1 0 (0%) 1 (0.9%) 0 (0%) 0 (0%) 0 (0%) 1 (0.2%)
bacillus.IMI.cow
0 117 (95.1%) 99 (84.6%) 56 (50.5%) 20 (95.2%) 45 (86.5%) 337 (79.5%)
1 6 (4.9%) 18 (15.4%) 55 (49.5%) 1 (4.8%) 7 (13.5%) 87 (20.5%)
coryne.IMI.cow
0 117 (95.1%) 108 (92.3%) 69 (62.2%) 19 (90.5%) 52 (100%) 365 (86.1%)
1 6 (4.9%) 9 (7.7%) 42 (37.8%) 2 (9.5%) 0 (0%) 59 (13.9%)
gram.pos.cocci.IMI.cow
0 123 (100%) 114 (97.4%) 100 (90.1%) 21 (100%) 52 (100%) 410 (96.7%)
1 0 (0%) 3 (2.6%) 11 (9.9%) 0 (0%) 0 (0%) 14 (3.3%)
gram.pos.rod.IMI.cow
0 123 (100%) 116 (99.1%) 111 (100%) 21 (100%) 52 (100%) 423 (99.8%)
1 0 (0%) 1 (0.9%) 0 (0%) 0 (0%) 0 (0%) 1 (0.2%)
helcoccus.IMI.cow
0 123 (100%) 117 (100%) 111 (100%) 21 (100%) 51 (98.1%) 423 (99.8%)
1 0 (0%) 0 (0%) 0 (0%) 0 (0%) 1 (1.9%) 1 (0.2%)
trueperella.IMI.cow
0 121 (98.4%) 116 (99.1%) 111 (100%) 21 (100%) 52 (100%) 421 (99.3%)
1 2 (1.6%) 1 (0.9%) 0 (0%) 0 (0%) 0 (0%) 3 (0.7%)
yeast.IMI.cow
0 123 (100%) 115 (98.3%) 111 (100%) 21 (100%) 52 (100%) 422 (99.5%)
1 0 (0%) 2 (1.7%) 0 (0%) 0 (0%) 0 (0%) 2 (0.5%)
# ---- Period prevalence by farm (visualization) ------

# ---- Staph aureus ----

sau.table.farm <- table(df.per.cow$staph.aureus.IMI.cow,df.per.cow$Site) %>%
  as.data.frame() %>% 
  group_by(Var2) %>%
  mutate(Total = sum(Freq)) %>%
  filter(Var1 == 1) %>%
  mutate(Perc = round((Freq/Total)*100,2))

group.1.farm <- sau.table.farm %>%
  ggplot(aes(x = Var2,y = Perc,fill = Var2)) +
  geom_bar(stat = "identity") +
  labs(y = "CI(%)",x = element_blank(),fill = "Farm",title = "SAU") +
  ylim(0,70) + scale_fill_viridis(option = "cividis",discrete = TRUE)

# ----- Staph chromogenes -------

sch.table.farm <- table(df.per.cow$sch.IMI.cow,df.per.cow$Site) %>%
  as.data.frame() %>% 
  group_by(Var2) %>%
  mutate(Total = sum(Freq)) %>%
  filter(Var1 == 1) %>%
  mutate(Perc = round((Freq/Total)*100,2))

group.2.farm <- sch.table.farm %>%
  ggplot(aes(x = Var2,y = Perc,fill = Var2)) +
  geom_bar(stat = "identity") +
  labs(y = "CI(%)",x = element_blank(),fill = "Farm",title = "SCH") +
  ylim(0,70) + scale_fill_viridis(option = "cividis",discrete = TRUE)

# ---- Staph chromogenes ----

nas.non.chr.table.farm <- table(df.per.cow$nas.non.chrom.IMI.cow,df.per.cow$Site) %>%
  as.data.frame() %>% 
  group_by(Var2) %>%
  mutate(Total = sum(Freq)) %>%
  filter(Var1 == 1) %>%
  mutate(Perc = round((Freq/Total)*100,2))

group.3.farm <- nas.non.chr.table.farm %>%
  ggplot(aes(x = Var2,y = Perc,fill = Var2)) +
  geom_bar(stat = "identity") +
  labs(y = "CI(%)",x = element_blank(),fill = "Farm",title = "NASM non-SCH") +
  ylim(0,70) + scale_fill_viridis(option = "cividis",discrete = TRUE)

# Strep genus 

strep.genus.table.farm <- table(df.per.cow$strep.genus.IMI.cow,df.per.cow$Site) %>%
  as.data.frame() %>% 
  group_by(Var2) %>%
  mutate(Total = sum(Freq)) %>%
  filter(Var1 == 1) %>%
  mutate(Perc = round((Freq/Total)*100,2))

group.4.farm <- strep.genus.table.farm %>%
  ggplot(aes(x = Var2,y = Perc,fill = Var2)) +
  geom_bar(stat = "identity") +
  labs(y = "CI(%)",x = element_blank(),fill = "Farm",title = "Strep") +
  ylim(0,70) + scale_fill_viridis(option = "cividis",discrete = TRUE)

# ---- Strep-like organisms ----

strep.like.org.table.farm <- table(df.per.cow$strep.like.IMI.cow,df.per.cow$Site) %>%
  as.data.frame() %>% 
  group_by(Var2) %>%
  mutate(Total = sum(Freq)) %>%
  filter(Var1 == 1) %>%
  mutate(Perc = round((Freq/Total)*100,2))

group.5.farm <- strep.like.org.table.farm %>%
  ggplot(aes(x = Var2,y = Perc,fill = Var2)) +
  geom_bar(stat = "identity") +
  labs(y = "CI(%)",x = element_blank(),fill = "Farm",title = "SLO") +
  ylim(0,70) + scale_fill_viridis(option = "cividis",discrete = TRUE)

# ---- Gram-negative ----

gram.neg.table.farm <- table(df.per.cow$gram.neg.IMI.cow,df.per.cow$Site) %>%
  as.data.frame() %>% 
  group_by(Var2) %>%
  mutate(Total = sum(Freq)) %>%
  filter(Var1 == 1) %>%
  mutate(Perc = round((Freq/Total)*100,2))

group.6.farm <- gram.neg.table.farm %>%
  ggplot(aes(x = Var2,y = Perc,fill = Var2)) +
  geom_bar(stat = "identity") +
  labs(y = "CI(%)",x = element_blank(),fill = "Farm",title = "GN") +
  ylim(0,70) + scale_fill_viridis(option = "cividis",discrete = TRUE)

# ---- Other ----

other.table.farm <- table(df.per.cow$other.IMI.cow,df.per.cow$Site) %>%
  as.data.frame() %>% 
  group_by(Var2) %>%
  mutate(Total = sum(Freq)) %>%
  filter(Var1 == 1) %>%
  mutate(Perc = round((Freq/Total)*100,2))

group.7.farm <- other.table.farm %>%
  ggplot(aes(x = Var2,y = Perc,fill = Var2)) +
  geom_bar(stat = "identity") +
  labs(y = "CI(%)",x = element_blank(),fill = "Farm",title = "Other") +
  ylim(0,70) + scale_fill_viridis(option = "cividis",discrete = TRUE)

groups.figure.farm <- ggarrange(group.1.farm + rremove("ylab") + rremove("xlab"), group.2.farm + rremove("ylab") + rremove("xlab"), group.3.farm + rremove("ylab") + rremove("xlab"),group.4.farm + rremove("ylab") + rremove("xlab"), group.5.farm + rremove("ylab") + rremove("xlab"), group.6.farm + rremove("ylab") + rremove("xlab"),group.7.farm + rremove("ylab") + rremove("xlab"),
                          ncol = 4, nrow = 2, common.legend = TRUE,align = "h",labels = c("A", "B", "C","D","E","F","G"))
## Warning: Removed 1 rows containing missing values (`position_stack()`).
plot_2 <- annotate_figure(groups.figure.farm,
                bottom = text_grob("Farm", color = "black",
                                             hjust = 0.75, x = 0.5),
                left = text_grob("Period prevalence (%)", color = "black", rot = 90))

plot_2

ggsave(plot = last_plot(),"./figures/plot_2_farm.png",width = 20, height = 10, units = "cm") 

1.4 Prevalence at calving and cumulative incidence

# ---- Prevalence at calving and Cumulative incidence  -----

# ---- Staph aureus -----

# Create table

sau.table <- table(df.per.cow$staph.aureus.IMI.cow,df.per.cow$calving.sau) %>%
as.data.frame() %>% 
group_by(Var2) %>%
mutate(Total = sum(Freq)) %>%
filter(Var1 == 1) %>%
mutate(Total = ifelse(Var2 == 1,424,Total)) %>%
mutate(Total = ifelse(Var2 == 0,(Total - 20),Total)) %>%
mutate(Perc = round((Freq/Total)*100,2))

sau.table$Var2 <- fct_recode(sau.table$Var2,"CI" = "0","PC" = "1")

sau.table$Var2 <- relevel(sau.table$Var2, ref = "PC")

# Show table

sau.table
## # A tibble: 2 × 5
## # Groups:   Var2 [2]
##   Var1  Var2   Freq Total  Perc
##   <fct> <fct> <int> <dbl> <dbl>
## 1 1     CI       18   342  5.26
## 2 1     PC       62   424 14.6
# Create figure

group.1 <- sau.table %>%
ggplot(aes(x = Var2,y = Perc,fill = Var2)) +
geom_bar(stat = "identity") +
labs(y = "%",x = element_blank(),fill = "Period",title = "SAU") +
ylim(0,40) + scale_fill_viridis(discrete = TRUE) 

# Proportion IMI in first sample

table(subset(df.per.cow,df.per.cow$numberpercow > 1)$staph.aureus.IMI.cow,subset(df.per.cow,df.per.cow$numberpercow > 1)$calving.sau)
##    
##       0   1
##   0 327   0
##   1  18  59
prop.table(table(subset(df.per.cow,df.per.cow$numberpercow > 1)$staph.aureus.IMI.cow,subset(df.per.cow,df.per.cow$numberpercow > 1)$calving.sau),1)
##    
##             0         1
##   0 1.0000000 0.0000000
##   1 0.2337662 0.7662338
# ----- Staph chromogenes -------

# Create table

sch.table <- table(df.per.cow$sch.IMI.cow,df.per.cow$calving.sch) %>%
  as.data.frame() %>% 
  group_by(Var2) %>%
  mutate(Total = sum(Freq)) %>%
  filter(Var1 == 1) %>%
  mutate(Total = ifelse(Var2 == 1,424,Total)) %>%
  mutate(Total = ifelse(Var2 == 0,(Total - 20),Total)) %>%
  mutate(Perc = round((Freq/Total)*100,2)) 

sch.table$Var2 <- fct_recode(sch.table$Var2,"CI" = "0","PC" = "1")

sch.table$Var2 <- relevel(sch.table$Var2, ref = "PC")

# Show table

sch.table
## # A tibble: 2 × 5
## # Groups:   Var2 [2]
##   Var1  Var2   Freq Total  Perc
##   <fct> <fct> <int> <dbl> <dbl>
## 1 1     CI       53   296  17.9
## 2 1     PC      108   424  25.5
# Create figure

group.2 <- sch.table %>%
  ggplot(aes(x = Var2,y = Perc,fill = Var2)) +
  geom_bar(stat = "identity") +
  labs(y = "%",x = element_blank(),fill = "Period",title = "SCH") +
  ylim(0,40) + scale_fill_viridis(discrete = TRUE) 

# Proportion IMI in first sample

table(subset(df.per.cow,df.per.cow$numberpercow > 1)$sch.IMI.cow,subset(df.per.cow,df.per.cow$numberpercow > 1)$calving.sch)
##    
##       0   1
##   0 247   0
##   1  53 104
prop.table(table(subset(df.per.cow,df.per.cow$numberpercow > 1)$sch.IMI.cow,subset(df.per.cow,df.per.cow$numberpercow > 1)$calving.sch),1)
##    
##             0         1
##   0 1.0000000 0.0000000
##   1 0.3375796 0.6624204
# ------ NAS-non chromogenes -------

# Create table

nas.non.ch.table <- table(df.per.cow$nas.non.chrom.IMI.cow,df.per.cow$calving.nas.nonch) %>%
  as.data.frame() %>% 
  group_by(Var2) %>%
  mutate(Total = sum(Freq)) %>%
  filter(Var1 == 1) %>%
  mutate(Total = ifelse(Var2 == 1,424,Total)) %>%
  mutate(Total = ifelse(Var2 == 0,(Total - 20),Total)) %>%
  mutate(Perc = round((Freq/Total)*100,2)) 

nas.non.ch.table$Var2 <- fct_recode(nas.non.ch.table$Var2,"CI" = "0","PC" = "1")

nas.non.ch.table$Var2 <- relevel(nas.non.ch.table$Var2, ref = "PC")

# Show table

nas.non.ch.table
## # A tibble: 2 × 5
## # Groups:   Var2 [2]
##   Var1  Var2   Freq Total  Perc
##   <fct> <fct> <int> <dbl> <dbl>
## 1 1     CI       66   374 17.6 
## 2 1     PC       30   424  7.08
# Create figure

group.3 <- nas.non.ch.table %>%
  ggplot(aes(x = Var2,y = Perc,fill = Var2)) +
  geom_bar(stat = "identity") +
  labs(y = "%",x = element_blank(),fill = "Period",title = "NASM non-SCH") +
  ylim(0,40) + scale_fill_viridis(discrete = TRUE)

# Proportion IMI in first sample

table(subset(df.per.cow,df.per.cow$numberpercow > 1)$nas.non.chrom.IMI.cow,subset(df.per.cow,df.per.cow$numberpercow > 1)$calving.nas.nonch)
##    
##       0   1
##   0 313   0
##   1  66  25
prop.table(table(subset(df.per.cow,df.per.cow$numberpercow > 1)$nas.non.chrom.IMI.cow,subset(df.per.cow,df.per.cow$numberpercow > 1)$calving.nas.nonch),1)
##    
##             0         1
##   0 1.0000000 0.0000000
##   1 0.7252747 0.2747253
# ---- Strep genus ---- 

# Create table

strep.genus.table <- table(df.per.cow$strep.genus.IMI.cow,df.per.cow$calving.strep.genus) %>%
  as.data.frame() %>% 
  group_by(Var2) %>%
  mutate(Total = sum(Freq)) %>%
  filter(Var1 == 1) %>%
  mutate(Total = ifelse(Var2 == 1,424,Total)) %>%
  mutate(Total = ifelse(Var2 == 0,(Total - 20),Total)) %>%
  mutate(Perc = round((Freq/Total)*100,2))

strep.genus.table$Var2 <- fct_recode(strep.genus.table$Var2,"CI" = "0","PC" = "1")

strep.genus.table$Var2 <- relevel(strep.genus.table$Var2, ref = "PC")

# Show table

strep.genus.table
## # A tibble: 2 × 5
## # Groups:   Var2 [2]
##   Var1  Var2   Freq Total  Perc
##   <fct> <fct> <int> <dbl> <dbl>
## 1 1     CI       22   355   6.2
## 2 1     PC       49   424  11.6
# Create figure

group.4 <- strep.genus.table %>%
  ggplot(aes(x = Var2,y = Perc,fill = Var2)) +
  geom_bar(stat="identity") +
  labs(y="%",x = element_blank(),fill = "Period",title = "Strep") +
  ylim(0,40) + scale_fill_viridis(discrete = TRUE)

# Proportion IMI in first sample

table(subset(df.per.cow,df.per.cow$numberpercow > 1)$strep.genus.IMI.cow,subset(df.per.cow,df.per.cow$numberpercow > 1)$calving.strep.genus)
##    
##       0   1
##   0 336   0
##   1  22  46
prop.table(table(subset(df.per.cow,df.per.cow$numberpercow > 1)$strep.genus.IMI.cow,subset(df.per.cow,df.per.cow$numberpercow > 1)$calving.strep.genus),1)
##    
##             0         1
##   0 1.0000000 0.0000000
##   1 0.3235294 0.6764706
# ---- Strep-like org ----

# Prevalence at calving and cumulative incidence

# Create table

strep.like.org.table <- table(df.per.cow$strep.like.IMI.cow,df.per.cow$calving.strep.like.org) %>%
  as.data.frame() %>% 
  group_by(Var2) %>%
  mutate(Total = sum(Freq)) %>%
  filter(Var1 == 1) %>%
  mutate(Total = ifelse(Var2 == 1,424,Total)) %>%
  mutate(Total = ifelse(Var2 == 0,(Total - 20),Total)) %>%
  mutate(Perc = round((Freq/Total)*100,2))

strep.like.org.table$Var2 <- fct_recode(strep.like.org.table$Var2,"CI" = "0","PC" = "1")

strep.like.org.table$Var2 <- relevel(strep.like.org.table$Var2, ref = "PC")

# Show table

strep.like.org.table
## # A tibble: 2 × 5
## # Groups:   Var2 [2]
##   Var1  Var2   Freq Total  Perc
##   <fct> <fct> <int> <dbl> <dbl>
## 1 1     CI       47   379  12.4
## 2 1     PC       25   424   5.9
# Create figure

 group.5 <- strep.like.org.table %>%
  ggplot(aes(x = Var2,y = Perc,fill = Var2)) +
  geom_bar(stat = "identity") +
  labs(y = "%",x = element_blank(),fill = "Period",title = "SLO") +
   ylim(0,40) + scale_fill_viridis(discrete = TRUE) 

# Proportion IMI in first sample
 
 table(subset(df.per.cow,df.per.cow$numberpercow > 1)$strep.like.IMI.cow,subset(df.per.cow,df.per.cow$numberpercow > 1)$calving.strep.like.org)
##    
##       0   1
##   0 333   0
##   1  47  24
prop.table(table(subset(df.per.cow,df.per.cow$numberpercow > 1)$strep.like.IMI.cow,subset(df.per.cow,df.per.cow$numberpercow > 1)$calving.strep.like.org),1)
##    
##             0         1
##   0 1.0000000 0.0000000
##   1 0.6619718 0.3380282
# ---- gram negative ----
 
# Create table

gram.negative.table <- table(df.per.cow$gram.neg.IMI.cow,df.per.cow$calving.gram.neg) %>%
  as.data.frame() %>% 
  group_by(Var2) %>%
  mutate(Total = sum(Freq)) %>%
  filter(Var1 == 1) %>%
  mutate(Total = ifelse(Var2 == 1,424,Total)) %>%
  mutate(Total = ifelse(Var2 == 0,(Total - 20),Total)) %>%
  mutate(Perc = round((Freq/Total)*100,2))

gram.negative.table$Var2 <- fct_recode(gram.negative.table$Var2,"CI" = "0","PC" = "1")

gram.negative.table$Var2 <- relevel(gram.negative.table$Var2, ref = "PC")

# Show table

gram.negative.table
## # A tibble: 2 × 5
## # Groups:   Var2 [2]
##   Var1  Var2   Freq Total  Perc
##   <fct> <fct> <int> <dbl> <dbl>
## 1 1     CI       25   392  6.38
## 2 1     PC       12   424  2.83
# Show figure

group.6 <- gram.negative.table %>%
  ggplot(aes(x = Var2,y = Perc,fill = Var2)) +
  geom_bar(stat = "identity") +
  labs(y = "%",x = element_blank(),fill = "Period",title = "GN") +
  ylim(0,40) + scale_fill_viridis(discrete = TRUE)

# Proportion IMI in first sample

table(subset(df.per.cow,df.per.cow$numberpercow > 1)$gram.neg.IMI.cow,subset(df.per.cow,df.per.cow$numberpercow > 1)$calving.gram.neg)
##    
##       0   1
##   0 367   0
##   1  25  12
prop.table(table(subset(df.per.cow,df.per.cow$numberpercow > 1)$gram.neg.IMI.cow,subset(df.per.cow,df.per.cow$numberpercow > 1)$calving.gram.neg),1)
##    
##             0         1
##   0 1.0000000 0.0000000
##   1 0.6756757 0.3243243
# ---- Other microorganisms ----

# Create table

others.table <- table(df.per.cow$other.IMI.cow,df.per.cow$calving.others) %>%
  as.data.frame() %>% 
  group_by(Var2) %>%
  mutate(Total = sum(Freq)) %>%
  filter(Var1 == 1) %>%
  mutate(Total = ifelse(Var2 == 1,424,Total)) %>%
  mutate(Total = ifelse(Var2 == 0,(Total - 20),Total)) %>%
  mutate(Perc = round((Freq/Total)*100,2))

others.table$Var2 <- fct_recode(others.table$Var2,"CI" = "0","PC" = "1")

others.table$Var2 <- relevel(others.table$Var2, ref = "PC")

# Show table

others.table
## # A tibble: 2 × 5
## # Groups:   Var2 [2]
##   Var1  Var2   Freq Total  Perc
##   <fct> <fct> <int> <dbl> <dbl>
## 1 1     CI       98   364 26.9 
## 2 1     PC       40   424  9.43
# Create figure

group.7 <- others.table %>%
  ggplot(aes(x = Var2,y = Perc,fill = Var2)) +
  geom_bar(stat = "identity") +
  labs(y = "%",x = element_blank(),fill = "Period",title = "Other") +
  ylim(0,40) + scale_fill_viridis(discrete = TRUE)

# Proportion IMI in first sample

table(subset(df.per.cow,df.per.cow$numberpercow > 1)$other.IMI.cow,subset(df.per.cow,df.per.cow$numberpercow > 1)$calving.others)
##    
##       0   1
##   0 267   0
##   1  98  39
prop.table(table(subset(df.per.cow,df.per.cow$numberpercow > 1)$other.IMI.cow,subset(df.per.cow,df.per.cow$numberpercow > 1)$calving.others),1)
##    
##             0         1
##   0 1.0000000 0.0000000
##   1 0.7153285 0.2846715
# Merge figures

groups.figure <- ggarrange(group.1 + rremove("ylab") + rremove("xlab"), group.2 + rremove("ylab") + rremove("xlab"), group.3 + rremove("ylab") + rremove("xlab"),group.4 + rremove("ylab") + rremove("xlab"), group.5 + rremove("ylab") + rremove("xlab"), group.6 + rremove("ylab") + rremove("xlab"), group.7 + rremove("ylab") + rremove("xlab"),
                  ncol = 4, nrow = 2, common.legend = TRUE,align = "h",labels = c("A", "B", "C","D","E","F","G"))

plot_1 <- annotate_figure(groups.figure,
                left = text_grob("%", color = "black", rot = 90))

plot_1

ggsave(plot = last_plot(),"./figures/plot_1_cav_ci.png",width = 20, height = 10, units = "cm") 

1.5 Persistence

# ---- Persistence by bacterial group -----

# Number of observations per farm

# All microorganisms

table(Persistence.analysis.long.4$Site)
## 
##   A   B   C   D   E 
##  99 204 256  43 119
# All except for gram-negative

table(subset(Persistence.analysis.long.4,Persistence.analysis.long.4$bacterial.group != "Gram-negative")$Site)
## 
##   A   B   C   D   E 
##  95 198 235  42 115
# Persistence per organism

table1(~number.weeks|organism , data = Persistence.analysis.long.4)
## Warning in .table1.internal(x = x, labels = labels, groupspan = groupspan, :
## Table has 43 columns. Are you sure this is what you want?
Aerococcus sp., n = 4
(N=4)
Aerococcus viridans n = 20
(N=20)
Arthrobacter gandavensis, n = 2
(N=2)
Arthrobacter sp., n = 1
(N=1)
Bacillus sp. n = 87
(N=87)
Citrobacter freundii, n = 1
(N=1)
Corynebacterium sp., n = 59
(N=59)
Enterobacter cloacae, n = 1
(N=1)
Enterobacter sp., n = 1
(N=1)
Enterococcus casseliflavus n = 13
(N=13)
Enterococcus faecalis n = 1
(N=1)
Enterococcus hirae, n = 9
(N=9)
Enterococcus mundtii, n = 10
(N=10)
Enterococcus saccarolyticus, n = 1
(N=1)
Enterococcus sp., n = 17
(N=17)
Escherichia coli, n = 13
(N=13)
Escherichia sp., n = 3
(N=3)
Helcococcus ovis, n = 1
(N=1)
Klebsiella oxytoca, n = 3
(N=3)
Lactococcus garvieae, n = 3
(N=3)
Lactococcus lactis, n = 2
(N=2)
Lactococcus sp. n = 1
(N=1)
Micrococcus sp., n = 3
(N=3)
Other gram-negative, n = 5
(N=5)
Other gram-positive coccci n = 14
(N=14)
Other gram-positive rod, n = 15
(N=15)
Pantoea sp., n = 3
(N=3)
Pseudomonas sp., n = 3
(N=3)
Serratia spp. n = 3
(N=3)
Staphylococcus aureus, n = 77
(N=77)
Staphylococcus chromogenes, n = 157
(N=157)
Staphylococcus epidermidis. n = 1
(N=1)
Staphylococcus haemolyticus, n = 9
(N=9)
Staphylococcus hominis, n = 3
(N=3)
Staphylococcus sciuri n = 3
(N=3)
Staphylococcus sp., n = 73
(N=73)
Staphylococcus xylosus, n = 9
(N=9)
Streptococcus dysgalactiae, n = 45
(N=45)
Streptococcus sp. n = 32
(N=32)
Streptococcus uberis, n = 8
(N=8)
Trueperella pyogenes, n = 3
(N=3)
Yeast, n = 2
(N=2)
Overall
(N=721)
number.weeks
1 4 (100%) 17 (85.0%) 2 (100%) 1 (100%) 67 (77.0%) 1 (100%) 47 (79.7%) 1 (100%) 1 (100%) 12 (92.3%) 1 (100%) 9 (100%) 10 (100%) 1 (100%) 14 (82.4%) 13 (100%) 3 (100%) 1 (100%) 3 (100%) 3 (100%) 1 (50.0%) 1 (100%) 3 (100%) 5 (100%) 14 (100%) 13 (86.7%) 3 (100%) 3 (100%) 3 (100%) 10 (13.0%) 49 (31.2%) 1 (100%) 9 (100%) 3 (100%) 3 (100%) 49 (67.1%) 9 (100%) 21 (46.7%) 22 (68.8%) 4 (50.0%) 3 (100%) 2 (100%) 442 (61.3%)
2 0 (0%) 2 (10.0%) 0 (0%) 0 (0%) 19 (21.8%) 0 (0%) 12 (20.3%) 0 (0%) 0 (0%) 1 (7.7%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 2 (11.8%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 1 (50.0%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 1 (6.7%) 0 (0%) 0 (0%) 0 (0%) 21 (27.3%) 37 (23.6%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 15 (20.5%) 0 (0%) 10 (22.2%) 9 (28.1%) 3 (37.5%) 0 (0%) 0 (0%) 133 (18.4%)
3 or + 0 (0%) 1 (5.0%) 0 (0%) 0 (0%) 1 (1.1%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 1 (5.9%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 1 (6.7%) 0 (0%) 0 (0%) 0 (0%) 46 (59.7%) 71 (45.2%) 0 (0%) 0 (0%) 0 (0%) 0 (0%) 9 (12.3%) 0 (0%) 14 (31.1%) 1 (3.1%) 1 (12.5%) 0 (0%) 0 (0%) 146 (20.2%)
# Persistence per bacterial group

table1(~number.weeks|bacterial.group , data = Persistence.analysis.long.4)
Staphylococcus aureus
(N=77)
Staphylococcus chromogenes
(N=157)
NASM non-chromogenes
(N=98)
Streptococcus spp.
(N=85)
Streptococcus-like org.
(N=84)
Gram-negative
(N=36)
Other
(N=184)
Overall
(N=721)
number.weeks
1 10 (13.0%) 49 (31.2%) 74 (75.5%) 47 (55.3%) 76 (90.5%) 36 (100%) 150 (81.5%) 442 (61.3%)
2 21 (27.3%) 37 (23.6%) 15 (15.3%) 22 (25.9%) 6 (7.1%) 0 (0%) 32 (17.4%) 133 (18.4%)
3 or + 46 (59.7%) 71 (45.2%) 9 (9.2%) 16 (18.8%) 2 (2.4%) 0 (0%) 2 (1.1%) 146 (20.2%)
# Persistence Staph aureus

Persistence.1 <- Persistence.Staph.aureus %>%
  ggplot(aes(fill = as.factor(number.weeks),y = n_repeated,x = organism)) + 
  geom_bar(position = "fill", stat = "identity") +
  labs(fill = "Number of samples",x = "",y = "Proportion") +
  coord_flip() +
  ggtitle("Staphylococcus aureus") +
  scale_fill_viridis(discrete = TRUE) 

# Persistence Staph chromogenes

Persistence.2 <- Persistence.Staph.chromogenes %>%
  ggplot(aes(fill = as.factor(number.weeks),y = n_repeated,x = organism)) + 
  geom_bar(position = "fill", stat = "identity") +
  labs(fill = "Number of samples",x = "",y = "Proportion") +
  coord_flip() +
  ggtitle("Staphylococcus chromogenes") + 
  scale_fill_viridis(discrete = TRUE) 

# Persistence NAS non chromogenes

Persistence.3 <- Persistence.NAS.non.chromogenes %>%
  arrange(variable) %>%
  ggplot(aes(fill = as.factor(number.weeks),y = n_repeated,x = organism)) + 
  geom_bar(position = "fill", stat = "identity") +
  labs(fill = "Number of samples",x = "",y = "Proportion") +
  coord_flip() + scale_x_discrete(limits = rev) +
  ggtitle("NASM non-chromogenes") + 
  scale_fill_viridis(discrete = TRUE) 

# Persistence Strep spp.

Persistence.4 <- Persistence.Strep %>%
  arrange(variable) %>%
  ggplot(aes(fill = as.factor(number.weeks),y = n_repeated,x = organism)) + 
  geom_bar(position = "fill", stat = "identity") +
  labs(fill = "Number of samples",x = "",y = "Proportion") +
  coord_flip() + scale_x_discrete(limits = rev) +
  ggtitle("Streptococcus spp.") + 
  scale_fill_viridis(discrete = TRUE) 

# Persistence Strep-like organisms

Persistence.5 <- Persistence.Strep.like %>%
  arrange(variable) %>%
  ggplot(aes(fill = as.factor(number.weeks),y = n_repeated,x = organism)) + 
  geom_bar(position = "fill", stat = "identity") +
  labs(fill = "Number of samples",x = "",y = "Proportion") +
  coord_flip() + scale_x_discrete(limits = rev) +
  ggtitle("Streptococcus-like organisms") + 
  scale_fill_viridis(discrete = TRUE) 

# Persistence gram-negative organisms

Persistence.6 <- Persistence.gram.neg %>%
  arrange(variable) %>%
  ggplot(aes(fill = as.factor(number.weeks),y = n_repeated,x = organism)) + 
  geom_bar(position = "fill", stat = "identity") +
  labs(fill = "Number of samples",x = "",y = "Proportion") +
  coord_flip() + scale_x_discrete(limits = rev) +
  ggtitle("Gram-negative") +
  scale_fill_viridis(discrete = TRUE)

# Persistence others

Persistence.7 <- Persistence.other %>%
  arrange(variable) %>%
  ggplot(aes(fill = as.factor(number.weeks),y = n_repeated,x = organism)) + 
  geom_bar(position = "fill", stat = "identity") +
  labs(fill = "Number of samples",x = "",y = "Proportion") +
  coord_flip() + scale_x_discrete(limits = rev) +
  ggtitle("Other") +
  scale_fill_viridis(discrete = TRUE)

# Create plot for all groups

persistence.plot <- ggarrange(Persistence.1 + rremove("xlab"),Persistence.2 + rremove("xlab"),Persistence.3 + rremove("xlab"),Persistence.4 + rremove("xlab"),Persistence.5 + rremove("xlab"),Persistence.6 + rremove("xlab"),Persistence.7 + rremove("xlab"),  ncol = 1, nrow = 7, common.legend = TRUE,align ="v",labels = c("A", "B", "C","D","E","F","H"),heights = c(1.5,1.5,3,2,5,4.5,4))

annotate_figure(persistence.plot,
                bottom = text_grob("Proportion", color = "black",
                                   hjust = 0.75, x = 0.5),
                left = text_grob("Taxonomic group", color = "black", rot = 90, size = 12))

ggsave(plot = last_plot(),"./figures/plot_persistence.png",width = 15, height = 30, units = "cm") 

2 Confirmatory Analysis

2.1 Prevalence by postpartum sample

# ---- Any IMI ----

# Built model

IMI.week <- glmer(data = df, IMI ~ as.factor(SampleOrder.overall) + Site + (1|CowId) ,family = binomial(link = "logit"))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0229252 (tol = 0.002, component 1)
# Summary model

summary(IMI.week)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: IMI ~ as.factor(SampleOrder.overall) + Site + (1 | CowId)
##    Data: df
## 
##      AIC      BIC   logLik deviance df.resid 
##   1998.0   2052.3   -989.0   1978.0     1663 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9758 -0.6551  0.1929  0.5965  1.7289 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  CowId  (Intercept) 2.051    1.432   
## Number of obs: 1673, groups:  CowId, 424
## 
## Fixed effects:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      -0.5528     0.2091  -2.644 0.008193 ** 
## as.factor(SampleOrder.overall)2  -0.1790     0.1760  -1.017 0.309239    
## as.factor(SampleOrder.overall)3   0.1288     0.1805   0.714 0.475514    
## as.factor(SampleOrder.overall)4  -0.2490     0.1870  -1.332 0.183000    
## as.factor(SampleOrder.overall)5  -0.1974     0.2419  -0.816 0.414369    
## SiteB                             1.3236     0.2585   5.121 3.04e-07 ***
## SiteC                             1.0257     0.2574   3.984 6.77e-05 ***
## SiteD                             1.6090     0.4862   3.309 0.000936 ***
## SiteE                             3.8877     0.4645   8.369  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) a.(SO.)2 a.(SO.)3 a.(SO.)4 a.(SO.)5 SiteB  SiteC  SiteD 
## as.fc(SO.)2 -0.394                                                         
## as.fc(SO.)3 -0.387  0.482                                                  
## as.fc(SO.)4 -0.349  0.468    0.461                                         
## as.fc(SO.)5 -0.224  0.363    0.360    0.362                                
## SiteB       -0.621 -0.024   -0.020   -0.058   -0.119                       
## SiteC       -0.620 -0.022   -0.022   -0.053   -0.092    0.544              
## SiteD       -0.343 -0.006    0.007    0.000    0.006    0.288  0.284       
## SiteE       -0.380 -0.020    0.019   -0.002    0.002    0.339  0.327  0.186
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.0229252 (tol = 0.002, component 1)
# Type III p-values

Anova(IMI.week,type = "III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: IMI
##                                  Chisq Df Pr(>Chisq)    
## (Intercept)                     6.9907  1   0.008193 ** 
## as.factor(SampleOrder.overall)  5.3571  4   0.252576    
## Site                           77.4508  4  6.037e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ICC

performance::icc(IMI.week,by_group = TRUE)
## # ICC by Group
## 
## Group |   ICC
## -------------
## CowId | 0.384
# ---- Staph aureus ----

# Built model

aureus.week <- glmer(data = subset(df,df$Site!="D"), `Staphylococcus aureus` ~ as.factor(SampleOrder.overall) + Site + (1|CowId) ,family = binomial(link = "logit"))
## Warning in (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf, :
## failure to converge in 10000 evaluations
## Warning in optwrap(optimizer, devfun, start, rho$lower, control = control, :
## convergence code 4 from Nelder_Mead: failure to converge in 10000 evaluations
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.016246 (tol = 0.002, component 1)
# Summary model

summary(aureus.week)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: `Staphylococcus aureus` ~ as.factor(SampleOrder.overall) + Site +  
##     (1 | CowId)
##    Data: subset(df, df$Site != "D")
## 
##      AIC      BIC   logLik deviance df.resid 
##    615.9    664.4   -299.0    597.9     1597 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.06260 -0.00818 -0.00579 -0.00413  2.44491 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  CowId  (Intercept) 184      13.56   
## Number of obs: 1606, groups:  CowId, 403
## 
## Fixed effects:
##                                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                     -10.08442    1.09250  -9.231   <2e-16 ***
## as.factor(SampleOrder.overall)2  -0.46772    0.47480  -0.985   0.3246    
## as.factor(SampleOrder.overall)3  -0.06941    0.47890  -0.145   0.8848    
## as.factor(SampleOrder.overall)4  -0.19501    0.49980  -0.390   0.6964    
## as.factor(SampleOrder.overall)5   0.06957    0.66614   0.104   0.9168    
## SiteB                             0.71867    1.17690   0.611   0.5414    
## SiteC                            -0.68906    1.49690  -0.460   0.6453    
## SiteE                             2.61488    1.33851   1.954   0.0508 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) a.(SO.)2 a.(SO.)3 a.(SO.)4 a.(SO.)5 SiteB  SiteC 
## as.fc(SO.)2 -0.162                                                  
## as.fc(SO.)3 -0.192  0.478                                           
## as.fc(SO.)4 -0.168  0.462    0.466                                  
## as.fc(SO.)5 -0.118  0.342    0.350    0.348                         
## SiteB       -0.680 -0.010   -0.014   -0.024   -0.081                
## SiteC       -0.535 -0.007   -0.011   -0.019   -0.051    0.511       
## SiteE       -0.588  0.014    0.027    0.034    0.015    0.567  0.446
## optimizer (Nelder_Mead) convergence code: 4 (failure to converge in 10000 evaluations)
## Model failed to converge with max|grad| = 0.016246 (tol = 0.002, component 1)
## failure to converge in 10000 evaluations
# Type III p-values

Anova(aureus.week,type = "III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Staphylococcus aureus
##                                  Chisq Df Pr(>Chisq)    
## (Intercept)                    85.2040  1     <2e-16 ***
## as.factor(SampleOrder.overall)  1.2791  4     0.8649    
## Site                            6.0331  3     0.1100    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ICC

performance::icc(aureus.week,by_group = TRUE)
## # ICC by Group
## 
## Group |   ICC
## -------------
## CowId | 0.982
# ---- Staph chromogenes ----

# Built model

s.chrom.week<-glmer(data = df,`Staphylococcus chromogenes`~as.factor(SampleOrder.overall) + Site + (1|CowId),family = binomial(link = "logit"))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.323579 (tol = 0.002, component 1)
# Summary model

summary(s.chrom.week)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: `Staphylococcus chromogenes` ~ as.factor(SampleOrder.overall) +  
##     Site + (1 | CowId)
##    Data: df
## 
##      AIC      BIC   logLik deviance df.resid 
##   1358.1   1412.3   -669.0   1338.1     1663 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.04003 -0.11295 -0.08140 -0.05292  2.71813 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  CowId  (Intercept) 19.06    4.366   
## Number of obs: 1673, groups:  CowId, 424
## 
## Fixed effects:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                     -3.94971    0.78352  -5.041 4.63e-07 ***
## as.factor(SampleOrder.overall)2 -0.45969    0.25870  -1.777   0.0756 .  
## as.factor(SampleOrder.overall)3 -0.08717    0.25895  -0.337   0.7364    
## as.factor(SampleOrder.overall)4 -0.66564    0.28965  -2.298   0.0216 *  
## as.factor(SampleOrder.overall)5 -0.38327    0.37496  -1.022   0.3067    
## SiteB                            0.48954    0.66556   0.736   0.4620    
## SiteC                           -0.87540    0.63016  -1.389   0.1648    
## SiteD                            0.97444    1.23823   0.787   0.4313    
## SiteE                            1.95455    0.93403   2.093   0.0364 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) a.(SO.)2 a.(SO.)3 a.(SO.)4 a.(SO.)5 SiteB  SiteC  SiteD 
## as.fc(SO.)2 -0.036                                                         
## as.fc(SO.)3 -0.105  0.462                                                  
## as.fc(SO.)4 -0.007  0.431    0.432                                         
## as.fc(SO.)5 -0.018  0.328    0.334    0.328                                
## SiteB       -0.347 -0.015   -0.028   -0.036   -0.071                       
## SiteC       -0.423 -0.017   -0.030   -0.035   -0.059    0.481              
## SiteD       -0.253 -0.011   -0.001   -0.007   -0.001    0.241  0.258       
## SiteE       -0.442 -0.036   -0.003   -0.022   -0.013    0.314  0.346  0.186
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.323579 (tol = 0.002, component 1)
# Type III p-values

Anova(s.chrom.week,type = "III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Staphylococcus chromogenes
##                                  Chisq Df Pr(>Chisq)    
## (Intercept)                    25.4116  1  4.631e-07 ***
## as.factor(SampleOrder.overall)  7.2981  4    0.12095    
## Site                           11.3115  4    0.02328 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ICC

performance::icc(s.chrom.week,by_group = TRUE)
## # ICC by Group
## 
## Group |   ICC
## -------------
## CowId | 0.853
# ---- NAS non-chromogenes ----

# Built model

nas.non.chrom.week <- glmer(data = df,nas.non.chrom~as.factor(SampleOrder.overall) + Site +  (1|CowId),family = binomial(link = "logit"))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.711704 (tol = 0.002, component 1)
# Summary model

summary(nas.non.chrom.week)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: nas.non.chrom ~ as.factor(SampleOrder.overall) + Site + (1 |  
##     CowId)
##    Data: df
## 
##      AIC      BIC   logLik deviance df.resid 
##    834.3    888.5   -407.1    814.3     1663 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1402 -0.1933 -0.1157 -0.0956  4.0843 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  CowId  (Intercept) 3.053    1.747   
## Number of obs: 1673, groups:  CowId, 424
## 
## Fixed effects:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      -4.5208     0.4635  -9.753  < 2e-16 ***
## as.factor(SampleOrder.overall)2   0.4578     0.2900   1.579  0.11439    
## as.factor(SampleOrder.overall)3   0.2352     0.3080   0.764  0.44504    
## as.factor(SampleOrder.overall)4   0.0197     0.3413   0.058  0.95398    
## as.factor(SampleOrder.overall)5   0.3811     0.4178   0.912  0.36172    
## SiteB                             1.2195     0.4193   2.909  0.00363 ** 
## SiteC                            -0.3113     0.4813  -0.647  0.51772    
## SiteD                             2.5833     0.6419   4.024 5.72e-05 ***
## SiteE                             2.1940     0.5017   4.373 1.23e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) a.(SO.)2 a.(SO.)3 a.(SO.)4 a.(SO.)5 SiteB  SiteC  SiteD 
## as.fc(SO.)2 -0.401                                                         
## as.fc(SO.)3 -0.361  0.521                                                  
## as.fc(SO.)4 -0.317  0.472    0.456                                         
## as.fc(SO.)5 -0.236  0.391    0.379    0.362                                
## SiteB       -0.628  0.012   -0.008   -0.032   -0.112                       
## SiteC       -0.466 -0.010   -0.021   -0.034   -0.089    0.556              
## SiteD       -0.504  0.040    0.030    0.032    0.035    0.434  0.349       
## SiteE       -0.601  0.044    0.050    0.055    0.039    0.546  0.449  0.385
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.711704 (tol = 0.002, component 1)
# Type III p-values

Anova(nas.non.chrom.week,type = "III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: nas.non.chrom
##                                  Chisq Df Pr(>Chisq)    
## (Intercept)                    95.1257  1  < 2.2e-16 ***
## as.factor(SampleOrder.overall)  3.3572  4     0.4999    
## Site                           39.8504  4  4.648e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ICC

performance::icc(nas.non.chrom.week,by_group = TRUE)
## # ICC by Group
## 
## Group |   ICC
## -------------
## CowId | 0.481
# ---- Strep spp. ----

# Built model

strep.genus.week <- glmer(data = df,strep.genus~as.factor(SampleOrder.overall) + Site +  (1|CowId),family = binomial(link = "logit"))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 1.47339 (tol = 0.002, component 1)
# Summary model

summary(strep.genus.week)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: strep.genus ~ as.factor(SampleOrder.overall) + Site + (1 | CowId)
##    Data: df
## 
##      AIC      BIC   logLik deviance df.resid 
##    619.1    673.3   -299.6    599.1     1663 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1998 -0.0224 -0.0133 -0.0089  4.2979 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  CowId  (Intercept) 44.45    6.667   
## Number of obs: 1673, groups:  CowId, 424
## 
## Fixed effects:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      -7.5534     0.8139  -9.281  < 2e-16 ***
## as.factor(SampleOrder.overall)2  -0.7724     0.4103  -1.882  0.05977 .  
## as.factor(SampleOrder.overall)3  -1.0355     0.4367  -2.371  0.01773 *  
## as.factor(SampleOrder.overall)4  -1.2636     0.4827  -2.617  0.00886 ** 
## as.factor(SampleOrder.overall)5  -2.1812     0.6891  -3.165  0.00155 ** 
## SiteB                             0.9966     0.8790   1.134  0.25687    
## SiteC                            -1.1117     1.2633  -0.880  0.37886    
## SiteD                            -0.5782     2.1328  -0.271  0.78632    
## SiteE                             3.6527     1.3246   2.758  0.00582 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) a.(SO.)2 a.(SO.)3 a.(SO.)4 a.(SO.)5 SiteB  SiteC  SiteD 
## as.fc(SO.)2 -0.105                                                         
## as.fc(SO.)3 -0.067  0.453                                                  
## as.fc(SO.)4 -0.043  0.415    0.422                                         
## as.fc(SO.)5  0.013  0.307    0.317    0.322                                
## SiteB       -0.491  0.012    0.014    0.008   -0.001                       
## SiteC       -0.290  0.024    0.029    0.028    0.026    0.481              
## SiteD       -0.183  0.015    0.024    0.025    0.025    0.280  0.210       
## SiteE       -0.192  0.040    0.061    0.070    0.066    0.500  0.391  0.224
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 1.47339 (tol = 0.002, component 1)
# Type III p-values

Anova(strep.genus.week,type = "III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: strep.genus
##                                 Chisq Df Pr(>Chisq)    
## (Intercept)                    86.134  1  < 2.2e-16 ***
## as.factor(SampleOrder.overall) 13.690  4   0.008354 ** 
## Site                           12.915  4   0.011698 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Adjusted means and multiple comparisons

emmeans(strep.genus.week,pairwise~SampleOrder.overall,type="response")$contrasts
##  contrast                                    odds.ratio    SE  df null z.ratio
##  SampleOrder.overall1 / SampleOrder.overall2       2.16 0.888 Inf    1   1.882
##  SampleOrder.overall1 / SampleOrder.overall3       2.82 1.230 Inf    1   2.371
##  SampleOrder.overall1 / SampleOrder.overall4       3.54 1.708 Inf    1   2.617
##  SampleOrder.overall1 / SampleOrder.overall5       8.86 6.103 Inf    1   3.165
##  SampleOrder.overall2 / SampleOrder.overall3       1.30 0.577 Inf    1   0.593
##  SampleOrder.overall2 / SampleOrder.overall4       1.63 0.795 Inf    1   1.009
##  SampleOrder.overall2 / SampleOrder.overall5       4.09 2.803 Inf    1   2.056
##  SampleOrder.overall3 / SampleOrder.overall4       1.26 0.623 Inf    1   0.460
##  SampleOrder.overall3 / SampleOrder.overall5       3.14 2.166 Inf    1   1.663
##  SampleOrder.overall4 / SampleOrder.overall5       2.50 1.759 Inf    1   1.306
##  p.value
##   0.3267
##   0.1232
##   0.0672
##   0.0134
##   0.9761
##   0.8512
##   0.2394
##   0.9908
##   0.4569
##   0.6874
## 
## Results are averaged over the levels of: Site 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio scale
# ICC

performance::icc(strep.genus.week,by_group = TRUE)
## # ICC by Group
## 
## Group |   ICC
## -------------
## CowId | 0.931
# ---- Strep.-like organisms ----

# Built model

strep.like.week <- glmer(data=df,strep.like~as.factor(SampleOrder.overall) + Site + (1|CowId),family = binomial(link = "logit"))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.336029 (tol = 0.002, component 1)
# Summary model

summary(strep.like.week)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: strep.like ~ as.factor(SampleOrder.overall) + Site + (1 | CowId)
##    Data: df
## 
##      AIC      BIC   logLik deviance df.resid 
##    623.9    678.1   -301.9    603.9     1663 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.0612 -0.2019 -0.1019 -0.0847  7.7828 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  CowId  (Intercept) 1.392    1.18    
## Number of obs: 1673, groups:  CowId, 424
## 
## Fixed effects:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                     -4.51728    0.49793  -9.072  < 2e-16 ***
## as.factor(SampleOrder.overall)2 -0.07596    0.32356  -0.235   0.8144    
## as.factor(SampleOrder.overall)3  0.04776    0.32510   0.147   0.8832    
## as.factor(SampleOrder.overall)4 -0.51915    0.38667  -1.343   0.1794    
## as.factor(SampleOrder.overall)5 -0.32586    0.50899  -0.640   0.5220    
## SiteB                           -0.29782    0.57616  -0.517   0.6052    
## SiteC                            1.96426    0.45335   4.333 1.47e-05 ***
## SiteD                            3.22040    0.59096   5.449 5.05e-08 ***
## SiteE                            1.35096    0.55934   2.415   0.0157 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) a.(SO.)2 a.(SO.)3 a.(SO.)4 a.(SO.)5 SiteB  SiteC  SiteD 
## as.fc(SO.)2 -0.298                                                         
## as.fc(SO.)3 -0.307  0.475                                                  
## as.fc(SO.)4 -0.214  0.401    0.405                                         
## as.fc(SO.)5 -0.148  0.306    0.313    0.274                                
## SiteB       -0.529 -0.007   -0.019   -0.029   -0.096                       
## SiteC       -0.746 -0.012   -0.020   -0.051   -0.096    0.619              
## SiteD       -0.638 -0.002    0.013   -0.013    0.011    0.461  0.616       
## SiteE       -0.587  0.004    0.020    0.021    0.012    0.495  0.631  0.487
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.336029 (tol = 0.002, component 1)
# Type III p-values

Anova(strep.like.week,type = "III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: strep.like
##                                  Chisq Df Pr(>Chisq)    
## (Intercept)                    82.3023  1  < 2.2e-16 ***
## as.factor(SampleOrder.overall)  2.6102  4      0.625    
## Site                           52.1242  4  1.299e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ICC

performance::icc(strep.like.week,by_group = TRUE)
## # ICC by Group
## 
## Group |   ICC
## -------------
## CowId | 0.297
# ---- Gram-negative organisms ----

# Built model

gram.neg.week <- glmer(data=df,gram.neg~as.factor(SampleOrder.overall) + Site +(1|CowId),family = binomial(link = "logit"))
## boundary (singular) fit: see help('isSingular')
# Summary model

summary(gram.neg.week)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: gram.neg ~ as.factor(SampleOrder.overall) + Site + (1 | CowId)
##    Data: df
## 
##      AIC      BIC   logLik deviance df.resid 
##    370.0    424.3   -175.0    350.0     1663 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.2532 -0.1742 -0.1259 -0.0995 10.0477 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  CowId  (Intercept) 4e-14    2e-07   
## Number of obs: 1673, groups:  CowId, 424
## 
## Fixed effects:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      -4.4588     0.5535  -8.055 7.95e-16 ***
## as.factor(SampleOrder.overall)2  -0.1558     0.4369  -0.357  0.72129    
## as.factor(SampleOrder.overall)3  -0.0967     0.4378  -0.221  0.82518    
## as.factor(SampleOrder.overall)4  -0.6865     0.5420  -1.267  0.20525    
## as.factor(SampleOrder.overall)5  -0.9803     0.7787  -1.259  0.20809    
## SiteB                             0.4108     0.6328   0.649  0.51629    
## SiteC                             1.7116     0.5482   3.123  0.00179 ** 
## SiteD                             0.4258     1.1262   0.378  0.70537    
## SiteE                             0.9643     0.7139   1.351  0.17678    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) a.(SO.)2 a.(SO.)3 a.(SO.)4 a.(SO.)5 SiteB  SiteC  SiteD 
## as.fc(SO.)2 -0.352                                                         
## as.fc(SO.)3 -0.344  0.456                                                  
## as.fc(SO.)4 -0.265  0.369    0.371                                         
## as.fc(SO.)5 -0.144  0.258    0.262    0.216                                
## SiteB       -0.710 -0.009   -0.022   -0.034   -0.084                       
## SiteC       -0.819 -0.012   -0.026   -0.041   -0.079    0.735              
## SiteD       -0.406  0.000    0.003    0.006    0.007    0.354  0.408       
## SiteE       -0.645  0.001    0.015    0.023    0.008    0.558  0.644  0.314
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# Type III p-values

Anova(gram.neg.week,type = "III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: gram.neg
##                                 Chisq Df Pr(>Chisq)    
## (Intercept)                    64.882  1  7.951e-16 ***
## as.factor(SampleOrder.overall)  2.881  4   0.577936    
## Site                           16.440  4   0.002482 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ICC

performance::icc(gram.neg.week,by_group = TRUE)
## Warning: Can't compute random effect variances. Some variance components equal
##   zero. Your model may suffer from singularity (see `?lme4::isSingular`
##   and `?performance::check_singularity`).
##   Solution: Respecify random structure! You may also decrease the
##   `tolerance` level to enforce the calculation of random effect variances.
## [1] NA
# ---- Other microorganisms ----

# Built model

other.week <- glmer(data = df,other~as.factor(SampleOrder.overall) + Site + (1|CowId),family = binomial(link = "logit"))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.15341 (tol = 0.002, component 1)
# Summary model

summary(other.week)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: other ~ as.factor(SampleOrder.overall) + Site + (1 | CowId)
##    Data: df
## 
##      AIC      BIC   logLik deviance df.resid 
##   1106.8   1161.0   -543.4   1086.8     1663 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.0296 -0.3023 -0.2222 -0.1672  5.1850 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  CowId  (Intercept) 0.5239   0.7238  
## Number of obs: 1673, groups:  CowId, 424
## 
## Fixed effects:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      -3.7207     0.3313 -11.231   <2e-16 ***
## as.factor(SampleOrder.overall)2   0.1877     0.2458   0.764   0.4450    
## as.factor(SampleOrder.overall)3   0.4492     0.2409   1.865   0.0622 .  
## as.factor(SampleOrder.overall)4   0.2829     0.2529   1.119   0.2633    
## as.factor(SampleOrder.overall)5   0.1555     0.3116   0.499   0.6177    
## SiteB                             0.8004     0.3288   2.434   0.0149 *  
## SiteC                             2.4299     0.3079   7.891    3e-15 ***
## SiteD                             0.1680     0.6871   0.244   0.8068    
## SiteE                             0.4934     0.4571   1.079   0.2804    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) a.(SO.)2 a.(SO.)3 a.(SO.)4 a.(SO.)5 SiteB  SiteC  SiteD 
## as.fc(SO.)2 -0.393                                                         
## as.fc(SO.)3 -0.411  0.533                                                  
## as.fc(SO.)4 -0.373  0.508    0.521                                         
## as.fc(SO.)5 -0.258  0.413    0.425    0.412                                
## SiteB       -0.672 -0.006   -0.015   -0.033   -0.098                       
## SiteC       -0.787  0.004    0.011   -0.015   -0.070    0.745              
## SiteD       -0.321  0.000    0.001    0.004    0.006    0.330  0.347       
## SiteE       -0.503  0.002    0.016    0.025    0.012    0.495  0.531  0.238
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.15341 (tol = 0.002, component 1)
# Type III p-values

Anova(other.week,type="III")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: other
##                                   Chisq Df Pr(>Chisq)    
## (Intercept)                    126.1261  1     <2e-16 ***
## as.factor(SampleOrder.overall)   3.7299  4     0.4438    
## Site                           100.6818  4     <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ICC

performance::icc(other.week,by_group = TRUE)
## # ICC by Group
## 
## Group |   ICC
## -------------
## CowId | 0.137

2.2 Period prevalence by farm

# ---- Any IMI ----

# Built model

IMI.farm <- glm(data = df.per.cow, IMI.cow ~Site  ,family = binomial(link = "logit"))

# Summary model

summary(IMI.farm)
## 
## Call:
## glm(formula = IMI.cow ~ Site, family = binomial(link = "logit"), 
##     data = df.per.cow)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.46760   0.00013   0.33338   0.44438   0.98128  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.4806     0.1856   2.590   0.0096 ** 
## SiteB         1.7850     0.3671   4.862 1.16e-06 ***
## SiteC         2.3816     0.4589   5.189 2.11e-07 ***
## SiteD         2.5151     1.0414   2.415   0.0157 *  
## SiteE        18.0855   904.5272   0.020   0.9840    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 363.28  on 423  degrees of freedom
## Residual deviance: 291.28  on 419  degrees of freedom
## AIC: 301.28
## 
## Number of Fisher Scoring iterations: 17
# Type III p-values

Anova(IMI.farm,type="III")
## Analysis of Deviance Table (Type III tests)
## 
## Response: IMI.cow
##      LR Chisq Df Pr(>Chisq)    
## Site   71.997  4  8.596e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Adjusted risk by farm and multiple comparisons

emmeans(IMI.farm,pairwise~Site,type = "response")
## $emmeans
##  Site     prob         SE  df asymp.LCL asymp.UCL
##  A    0.617886 0.04381250 Inf  0.529188  0.699367
##  B    0.905983 0.02698177 Inf  0.838172  0.947170
##  C    0.945946 0.02146277 Inf  0.884883  0.975515
##  D    0.952381 0.04647143 Inf  0.728568  0.993334
##  E    1.000000 0.00000782 Inf  0.000000  1.000000
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast odds.ratio       SE  df null z.ratio p.value
##  A / B      1.68e-01 6.16e-02 Inf    1  -4.862  <.0001
##  A / C      9.24e-02 4.24e-02 Inf    1  -5.189  <.0001
##  A / D      8.09e-02 8.42e-02 Inf    1  -2.415  0.1112
##  A / E      0.00e+00 1.26e-05 Inf    1  -0.020  1.0000
##  B / C      5.51e-01 2.90e-01 Inf    1  -1.135  0.7882
##  B / D      4.82e-01 5.17e-01 Inf    1  -0.681  0.9607
##  B / E      1.00e-07 7.54e-05 Inf    1  -0.018  1.0000
##  C / D      8.75e-01 9.69e-01 Inf    1  -0.121  1.0000
##  C / E      2.00e-07 1.37e-04 Inf    1  -0.017  1.0000
##  D / E      2.00e-07 1.56e-04 Inf    1  -0.017  1.0000
## 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio scale
# ---- Staph aureus ----

# Built model

aureus.farm <- glm(data=subset(df.per.cow,df.per.cow$Site!="D"), staph.aureus.IMI.cow~Site  ,family = binomial(link = "logit"))

# Summary model

summary(aureus.farm)
## 
## Call:
## glm(formula = staph.aureus.IMI.cow ~ Site, family = binomial(link = "logit"), 
##     data = subset(df.per.cow, df.per.cow$Site != "D"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1127  -0.7846  -0.5100  -0.4345   2.1941  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.9741     0.2755  -7.164 7.82e-13 ***
## SiteB         0.9537     0.3461   2.755  0.00586 ** 
## SiteC        -0.3385     0.4310  -0.785  0.43231    
## SiteE         1.8199     0.3915   4.648 3.35e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 401.66  on 402  degrees of freedom
## Residual deviance: 365.50  on 399  degrees of freedom
## AIC: 373.5
## 
## Number of Fisher Scoring iterations: 4
# Type III p-values

Anova(aureus.farm,type = "III")
## Analysis of Deviance Table (Type III tests)
## 
## Response: staph.aureus.IMI.cow
##      LR Chisq Df Pr(>Chisq)    
## Site   36.157  3  6.937e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Adjusted risk by farm and multiple comparisons

emmeans(aureus.farm,pairwise~Site,type = "response")
## $emmeans
##  Site   prob     SE  df asymp.LCL asymp.UCL
##  A    0.1220 0.0295 Inf    0.0749     0.192
##  B    0.2650 0.0408 Inf    0.1930     0.352
##  C    0.0901 0.0272 Inf    0.0492     0.159
##  E    0.4615 0.0691 Inf    0.3320     0.597
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast odds.ratio     SE  df null z.ratio p.value
##  A / B         0.385 0.1334 Inf    1  -2.755  0.0299
##  A / C         1.403 0.6046 Inf    1   0.785  0.8612
##  A / E         0.162 0.0634 Inf    1  -4.648  <.0001
##  B / C         3.641 1.4275 Inf    1   3.296  0.0054
##  B / E         0.421 0.1464 Inf    1  -2.487  0.0619
##  C / E         0.116 0.0500 Inf    1  -4.988  <.0001
## 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## Tests are performed on the log odds ratio scale
# ---- Staph chromogenes ----

# Built model

s.chrom.farm <- glm(data = df.per.cow,sch.IMI.cow~Site,family = binomial(link = "logit"))

# Summary model

summary(s.chrom.farm)
## 
## Call:
## glm(formula = sch.IMI.cow ~ Site, family = binomial(link = "logit"), 
##     data = df.per.cow)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1449  -0.9460  -0.8704   1.3194   1.5193  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.6568     0.1901  -3.454 0.000552 ***
## SiteB         0.3291     0.2670   1.233 0.217678    
## SiteC        -0.1186     0.2791  -0.425 0.670843    
## SiteD         0.5615     0.4765   1.178 0.238682    
## SiteE         0.5798     0.3364   1.723 0.084818 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 563.01  on 423  degrees of freedom
## Residual deviance: 556.47  on 419  degrees of freedom
## AIC: 566.47
## 
## Number of Fisher Scoring iterations: 4
# Type III p-values

Anova(s.chrom.farm,type = "III")
## Analysis of Deviance Table (Type III tests)
## 
## Response: sch.IMI.cow
##      LR Chisq Df Pr(>Chisq)
## Site   6.5344  4     0.1626
# Adjusted risk by farm and multiple comparisons

emmeans(s.chrom.farm,pairwise~Site,type = "response")
## $emmeans
##  Site  prob     SE  df asymp.LCL asymp.UCL
##  A    0.341 0.0428 Inf     0.263     0.429
##  B    0.419 0.0456 Inf     0.333     0.510
##  C    0.315 0.0441 Inf     0.236     0.407
##  D    0.476 0.1090 Inf     0.279     0.682
##  E    0.481 0.0693 Inf     0.350     0.615
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast odds.ratio    SE  df null z.ratio p.value
##  A / B         0.720 0.192 Inf    1  -1.233  0.7322
##  A / C         1.126 0.314 Inf    1   0.425  0.9932
##  A / D         0.570 0.272 Inf    1  -1.178  0.7639
##  A / E         0.560 0.188 Inf    1  -1.723  0.4194
##  B / C         1.565 0.434 Inf    1   1.615  0.4876
##  B / D         0.793 0.377 Inf    1  -0.489  0.9884
##  B / E         0.778 0.261 Inf    1  -0.749  0.9449
##  C / D         0.507 0.244 Inf    1  -1.410  0.6212
##  C / E         0.497 0.171 Inf    1  -2.027  0.2531
##  D / E         0.982 0.508 Inf    1  -0.035  1.0000
## 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio scale
# ---- NAS-non chromogenes -----

# Built model

nas.non.chrom.farm <- glm(data = df.per.cow,nas.non.chrom.IMI.cow~Site,family = binomial(link = "logit"))

# Summary model

summary(nas.non.chrom.farm)
## 
## Call:
## glm(formula = nas.non.chrom.IMI.cow ~ Site, family = binomial(link = "logit"), 
##     data = df.per.cow)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3893  -0.8431  -0.4991  -0.4916   2.0848  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.05229    0.28390  -7.229 4.87e-13 ***
## SiteB        1.20092    0.34838   3.447 0.000566 ***
## SiteC        0.03227    0.40954   0.079 0.937191    
## SiteD        2.53780    0.53153   4.775 1.80e-06 ***
## SiteE        1.66283    0.40060   4.151 3.31e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 453.6  on 423  degrees of freedom
## Residual deviance: 408.2  on 419  degrees of freedom
## AIC: 418.2
## 
## Number of Fisher Scoring iterations: 4
# Type III p-values

Anova(nas.non.chrom.farm,type = "III")
## Analysis of Deviance Table (Type III tests)
## 
## Response: nas.non.chrom.IMI.cow
##      LR Chisq Df Pr(>Chisq)    
## Site   45.404  4  3.276e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Adjusted risk by farm and multiple comparisons

emmeans(nas.non.chrom.farm,pairwise~Site,type = "response")
## $emmeans
##  Site  prob     SE  df asymp.LCL asymp.UCL
##  A    0.114 0.0286 Inf    0.0686     0.183
##  B    0.299 0.0423 Inf    0.2232     0.388
##  C    0.117 0.0305 Inf    0.0692     0.191
##  D    0.619 0.1060 Inf    0.4025     0.797
##  E    0.404 0.0680 Inf    0.2802     0.541
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast odds.ratio     SE  df null z.ratio p.value
##  A / B        0.3009 0.1048 Inf    1  -3.447  0.0051
##  A / C        0.9682 0.3965 Inf    1  -0.079  1.0000
##  A / D        0.0790 0.0420 Inf    1  -4.775  <.0001
##  A / E        0.1896 0.0760 Inf    1  -4.151  0.0003
##  B / C        3.2176 1.1507 Inf    1   3.268  0.0096
##  B / D        0.2627 0.1294 Inf    1  -2.714  0.0520
##  B / E        0.6301 0.2189 Inf    1  -1.330  0.6725
##  C / D        0.0816 0.0439 Inf    1  -4.660  <.0001
##  C / E        0.1958 0.0800 Inf    1  -3.990  0.0006
##  D / E        2.3988 1.2734 Inf    1   1.648  0.4664
## 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio scale
# ---- Strep genus ----

# Built model

strep.genus.farm <- glm(data = df.per.cow,strep.genus.IMI.cow~Site,family = binomial(link = "logit"))

# Summary model

summary(strep.genus.farm)
## 
## Call:
## glm(formula = strep.genus.IMI.cow ~ Site, family = binomial(link = "logit"), 
##     data = df.per.cow)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1774  -0.7548  -0.4112  -0.2881   2.5309  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -3.1612     0.4565  -6.925 4.35e-12 ***
## SiteB         2.0512     0.5042   4.068 4.74e-05 ***
## SiteC         0.7335     0.5738   1.278    0.201    
## SiteD         0.9100     0.8724   1.043    0.297    
## SiteE         3.1612     0.5341   5.918 3.25e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 383.15  on 423  degrees of freedom
## Residual deviance: 320.62  on 419  degrees of freedom
## AIC: 330.62
## 
## Number of Fisher Scoring iterations: 5
# Type III p-values

Anova(strep.genus.farm,type = "III")
## Analysis of Deviance Table (Type III tests)
## 
## Response: strep.genus.IMI.cow
##      LR Chisq Df Pr(>Chisq)    
## Site   62.524  4  8.546e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Adjusted risk by farm and multiple comparisons

emmeans(strep.genus.farm,pairwise~Site,type = "response")
## $emmeans
##  Site   prob     SE  df asymp.LCL asymp.UCL
##  A    0.0407 0.0178 Inf    0.0170    0.0939
##  B    0.2479 0.0399 Inf    0.1780    0.3339
##  C    0.0811 0.0259 Inf    0.0427    0.1485
##  D    0.0952 0.0641 Inf    0.0239    0.3113
##  E    0.5000 0.0693 Inf    0.3674    0.6326
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast odds.ratio     SE  df null z.ratio p.value
##  A / B        0.1286 0.0648 Inf    1  -4.068  0.0005
##  A / C        0.4802 0.2756 Inf    1  -1.278  0.7047
##  A / D        0.4025 0.3512 Inf    1  -1.043  0.8353
##  A / E        0.0424 0.0226 Inf    1  -5.918  <.0001
##  B / C        3.7348 1.5252 Inf    1   3.227  0.0110
##  B / D        3.1307 2.4219 Inf    1   1.475  0.5787
##  B / E        0.3295 0.1155 Inf    1  -3.168  0.0133
##  C / D        0.8382 0.6879 Inf    1  -0.215  0.9995
##  C / E        0.0882 0.0392 Inf    1  -5.458  <.0001
##  D / E        0.1053 0.0835 Inf    1  -2.837  0.0367
## 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio scale
# ---- Strep.-like organisms ----

# Built model

strep.like.farm <- glm(data=df.per.cow,strep.like.IMI.cow~Site,family = binomial(link = "logit"))

# Summary model

summary(strep.like.farm)
## 
## Call:
## glm(formula = strep.like.IMI.cow ~ Site, family = binomial(link = "logit"), 
##     data = df.per.cow)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1372  -0.6165  -0.3423  -0.3245   2.4374  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.8077     0.3892  -7.214 5.43e-13 ***
## SiteB        -0.1101     0.5720  -0.192    0.847    
## SiteC         2.2339     0.4365   5.117 3.10e-07 ***
## SiteD         2.7124     0.5851   4.635 3.56e-06 ***
## SiteE         1.2437     0.5346   2.326    0.020 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 386.34  on 423  degrees of freedom
## Residual deviance: 323.14  on 419  degrees of freedom
## AIC: 333.14
## 
## Number of Fisher Scoring iterations: 5
# Type III p-values

Anova(strep.like.farm,type = "III")
## Analysis of Deviance Table (Type III tests)
## 
## Response: strep.like.IMI.cow
##      LR Chisq Df Pr(>Chisq)    
## Site   63.199  4  6.163e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Adjusted risk by farm and multiple comparisons

emmeans(strep.like.farm,pairwise~Site,type = "response")
## $emmeans
##  Site   prob     SE  df asymp.LCL asymp.UCL
##  A    0.0569 0.0209 Inf    0.0274     0.115
##  B    0.0513 0.0204 Inf    0.0232     0.109
##  C    0.3604 0.0456 Inf    0.2766     0.454
##  D    0.4762 0.1090 Inf    0.2785     0.682
##  E    0.1731 0.0525 Inf    0.0926     0.300
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast odds.ratio     SE  df null z.ratio p.value
##  A / B        1.1164 0.6385 Inf    1   0.192  0.9997
##  A / C        0.1071 0.0468 Inf    1  -5.117  <.0001
##  A / D        0.0664 0.0388 Inf    1  -4.635  <.0001
##  A / E        0.2883 0.1541 Inf    1  -2.326  0.1365
##  B / C        0.0959 0.0445 Inf    1  -5.058  <.0001
##  B / D        0.0595 0.0360 Inf    1  -4.662  <.0001
##  B / E        0.2583 0.1438 Inf    1  -2.431  0.1070
##  C / D        0.6197 0.2972 Inf    1  -0.998  0.8565
##  C / E        2.6917 1.1210 Inf    1   2.378  0.1214
##  D / E        4.3434 2.4772 Inf    1   2.575  0.0750
## 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio scale
# ---- Gram-negative organisms ----

gram.neg.farm <- glm(data=df.per.cow,gram.neg.IMI.cow~Site,family = binomial(link = "logit"))

# Summary model

summary(gram.neg.farm)
## 
## Call:
## glm(formula = gram.neg.IMI.cow ~ Site, family = binomial(link = "logit"), 
##     data = df.per.cow)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.6815  -0.3447  -0.3245  -0.2571   2.6176  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -3.3928     0.5083  -6.674 2.48e-11 ***
## SiteB         0.4751     0.6588   0.721 0.470881    
## SiteC         2.0510     0.5597   3.665 0.000248 ***
## SiteD         0.3971     1.1439   0.347 0.728474    
## SiteE         0.5996     0.7824   0.766 0.443444    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 251.15  on 423  degrees of freedom
## Residual deviance: 226.86  on 419  degrees of freedom
## AIC: 236.86
## 
## Number of Fisher Scoring iterations: 6
# Type III p-values

Anova(gram.neg.farm,type = "III")
## Analysis of Deviance Table (Type III tests)
## 
## Response: gram.neg.IMI.cow
##      LR Chisq Df Pr(>Chisq)    
## Site   24.286  4  6.999e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(gram.neg.farm,pairwise~Site,type = "response")
## $emmeans
##  Site   prob     SE  df asymp.LCL asymp.UCL
##  A    0.0325 0.0160 Inf   0.01226    0.0834
##  B    0.0513 0.0204 Inf   0.02322    0.1095
##  C    0.2072 0.0385 Inf   0.14175    0.2926
##  D    0.0476 0.0465 Inf   0.00667    0.2714
##  E    0.0577 0.0323 Inf   0.01873    0.1642
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast odds.ratio     SE  df null z.ratio p.value
##  A / B         0.622 0.4097 Inf    1  -0.721  0.9517
##  A / C         0.129 0.0720 Inf    1  -3.665  0.0023
##  A / D         0.672 0.7690 Inf    1  -0.347  0.9969
##  A / E         0.549 0.4296 Inf    1  -0.766  0.9402
##  B / C         0.207 0.0993 Inf    1  -3.282  0.0091
##  B / D         1.081 1.1969 Inf    1   0.070  1.0000
##  B / E         0.883 0.6424 Inf    1  -0.171  0.9998
##  C / D         5.227 5.4945 Inf    1   1.573  0.5145
##  C / E         4.269 2.7287 Inf    1   2.271  0.1545
##  D / E         0.817 0.9676 Inf    1  -0.171  0.9998
## 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio scale

2.3 Persistence

# ---- Persistence by microorganism ----

# Built model

Model1 <- glmer(data = subset(Persistence.analysis.long.4,Persistence.analysis.long.4$bacterial.group != "Gram-negative"), persistent~bacterial.group + numberpercow + Site + (1|CowId) ,family = binomial(link = "logit"))
## boundary (singular) fit: see help('isSingular')
# Summary model

summary(Model1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: persistent ~ bacterial.group + numberpercow + Site + (1 | CowId)
##    Data: 
## subset(Persistence.analysis.long.4, Persistence.analysis.long.4$bacterial.group !=  
##     "Gram-negative")
## 
##      AIC      BIC   logLik deviance df.resid 
##    713.4    767.8   -344.7    689.4      673 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2809 -0.5355 -0.3537  0.5716  5.0880 
## 
## Random effects:
##  Groups Name        Variance  Std.Dev. 
##  CowId  (Intercept) 1.882e-14 1.372e-07
## Number of obs: 685, groups:  CowId, 344
## 
## Fixed effects:
##                                           Estimate Std. Error z value Pr(>|z|)
## (Intercept)                                -0.3887     0.6423  -0.605  0.54514
## bacterial.groupStaphylococcus chromogenes  -1.0887     0.3895  -2.795  0.00519
## bacterial.groupNASM non-chromogenes        -3.0983     0.4257  -7.278 3.38e-13
## bacterial.groupStreptococcus spp.          -2.2959     0.4152  -5.530 3.21e-08
## bacterial.groupStreptococcus-like org.     -4.1933     0.5267  -7.961 1.71e-15
## bacterial.groupOther                       -3.4564     0.4150  -8.328  < 2e-16
## numberpercow                                0.5874     0.1426   4.119 3.81e-05
## SiteB                                      -0.2065     0.3347  -0.617  0.53720
## SiteC                                      -0.3412     0.3408  -1.001  0.31670
## SiteD                                       0.1533     0.4713   0.325  0.74492
## SiteE                                       0.4151     0.3485   1.191  0.23368
##                                              
## (Intercept)                                  
## bacterial.groupStaphylococcus chromogenes ** 
## bacterial.groupNASM non-chromogenes       ***
## bacterial.groupStreptococcus spp.         ***
## bacterial.groupStreptococcus-like org.    ***
## bacterial.groupOther                      ***
## numberpercow                              ***
## SiteB                                        
## SiteC                                        
## SiteD                                        
## SiteE                                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) bct.Sc b.NASn bc.Ss. b.S-o. bctr.O nmbrpr SiteB  SiteC 
## bctrl.grpSc -0.469                                                        
## bctr.NASMn- -0.336  0.726                                                 
## bctrl.grSs. -0.305  0.725  0.687                                          
## bctrl.gS-o. -0.254  0.601  0.557  0.552                                   
## bctrl.grpOt -0.334  0.757  0.692  0.697  0.617                            
## numberpercw -0.772 -0.031 -0.108 -0.129 -0.107 -0.119                     
## SiteB        0.020  0.055 -0.003 -0.045  0.033  0.011 -0.398              
## SiteC       -0.001  0.005 -0.001 -0.021 -0.119 -0.167 -0.326  0.696       
## SiteD       -0.230 -0.044 -0.149 -0.078 -0.139 -0.073  0.083  0.372  0.374
## SiteE       -0.436  0.079 -0.022 -0.119 -0.014  0.013  0.177  0.502  0.470
##             SiteD 
## bctrl.grpSc       
## bctr.NASMn-       
## bctrl.grSs.       
## bctrl.gS-o.       
## bctrl.grpOt       
## numberpercw       
## SiteB             
## SiteC             
## SiteD             
## SiteE        0.402
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
# Type III p-values anova

Anova(Model1,type="III")
## Warning in printHypothesis(L, rhs, names(b)): one or more coefficients in the hypothesis include
##      arithmetic operators in their names;
##   the printed representation of the hypothesis will be omitted
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: persistent
##                    Chisq Df Pr(>Chisq)    
## (Intercept)       0.3661  1     0.5451    
## bacterial.group 142.6896  5  < 2.2e-16 ***
## numberpercow     16.9629  1  3.812e-05 ***
## Site              4.8118  4     0.3072    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# OR and 95%CI 

exp(Confint(Model1))
##                                             Estimate       2.5 %     97.5 %
## (Intercept)                               0.67796969 0.192507697 2.38765987
## bacterial.groupStaphylococcus chromogenes 0.33665644 0.156911500 0.72230241
## bacterial.groupNASM non-chromogenes       0.04512540 0.019591560 0.10393772
## bacterial.groupStreptococcus spp.         0.10067551 0.044617840 0.22716379
## bacterial.groupStreptococcus-like org.    0.01509583 0.005376378 0.04238619
## bacterial.groupOther                      0.03154232 0.013983793 0.07114792
## numberpercow                              1.79939213 1.360555572 2.37977198
## SiteB                                     0.81339239 0.422071490 1.56752398
## SiteC                                     0.71092958 0.364567171 1.38635868
## SiteD                                     1.16570191 0.462853207 2.93583566
## SiteE                                     1.51448251 0.764880569 2.99871296
# Adjusted risk by microorganism and multiple comparisons

emmeans(Model1,pairwise~bacterial.group,type = "response")
## $emmeans
##  bacterial.group             prob     SE  df asymp.LCL asymp.UCL
##  Staphylococcus aureus      0.885 0.0367 Inf    0.7915     0.940
##  Staphylococcus chromogenes 0.721 0.0388 Inf    0.6395     0.791
##  NASM non-chromogenes       0.258 0.0467 Inf    0.1770     0.359
##  Streptococcus spp.         0.436 0.0602 Inf    0.3239     0.556
##  Streptococcus-like org.    0.104 0.0359 Inf    0.0518     0.198
##  Other                      0.195 0.0357 Inf    0.1344     0.275
## 
## Results are averaged over the levels of: Site 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast                                               odds.ratio     SE  df
##  Staphylococcus aureus / Staphylococcus chromogenes          2.970  1.157 Inf
##  Staphylococcus aureus / (NASM non-chromogenes)             22.160  9.434 Inf
##  Staphylococcus aureus / Streptococcus spp.                  9.933  4.124 Inf
##  Staphylococcus aureus / (Streptococcus-like org.)          66.243 34.893 Inf
##  Staphylococcus aureus / Other                              31.703 13.158 Inf
##  Staphylococcus chromogenes / (NASM non-chromogenes)         7.460  2.267 Inf
##  Staphylococcus chromogenes / Streptococcus spp.             3.344  1.001 Inf
##  Staphylococcus chromogenes / (Streptococcus-like org.)     22.301  9.527 Inf
##  Staphylococcus chromogenes / Other                         10.673  3.003 Inf
##  (NASM non-chromogenes) / Streptococcus spp.                 0.448  0.149 Inf
##  (NASM non-chromogenes) / (Streptococcus-like org.)          2.989  1.366 Inf
##  (NASM non-chromogenes) / Other                              1.431  0.472 Inf
##  Streptococcus spp. / (Streptococcus-like org.)              6.669  3.044 Inf
##  Streptococcus spp. / Other                                  3.192  1.032 Inf
##  (Streptococcus-like org.) / Other                           0.479  0.203 Inf
##  null z.ratio p.value
##     1   2.795  0.0582
##     1   7.278  <.0001
##     1   5.530  <.0001
##     1   7.961  <.0001
##     1   8.328  <.0001
##     1   6.614  <.0001
##     1   4.034  0.0008
##     1   7.268  <.0001
##     1   8.416  <.0001
##     1  -2.410  0.1525
##     1   2.397  0.1571
##     1   1.086  0.8872
##     1   4.157  0.0005
##     1   3.589  0.0045
##     1  -1.738  0.5066
## 
## Results are averaged over the levels of: Site 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## Tests are performed on the log odds ratio scale
emmeans_2 <- emmeans(Model1,pairwise~bacterial.group,type = "response")
confint(emmeans_2)
## $emmeans
##  bacterial.group             prob     SE  df asymp.LCL asymp.UCL
##  Staphylococcus aureus      0.885 0.0367 Inf    0.7915     0.940
##  Staphylococcus chromogenes 0.721 0.0388 Inf    0.6395     0.791
##  NASM non-chromogenes       0.258 0.0467 Inf    0.1770     0.359
##  Streptococcus spp.         0.436 0.0602 Inf    0.3239     0.556
##  Streptococcus-like org.    0.104 0.0359 Inf    0.0518     0.198
##  Other                      0.195 0.0357 Inf    0.1344     0.275
## 
## Results are averaged over the levels of: Site 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast                                               odds.ratio     SE  df
##  Staphylococcus aureus / Staphylococcus chromogenes          2.970  1.157 Inf
##  Staphylococcus aureus / (NASM non-chromogenes)             22.160  9.434 Inf
##  Staphylococcus aureus / Streptococcus spp.                  9.933  4.124 Inf
##  Staphylococcus aureus / (Streptococcus-like org.)          66.243 34.893 Inf
##  Staphylococcus aureus / Other                              31.703 13.158 Inf
##  Staphylococcus chromogenes / (NASM non-chromogenes)         7.460  2.267 Inf
##  Staphylococcus chromogenes / Streptococcus spp.             3.344  1.001 Inf
##  Staphylococcus chromogenes / (Streptococcus-like org.)     22.301  9.527 Inf
##  Staphylococcus chromogenes / Other                         10.673  3.003 Inf
##  (NASM non-chromogenes) / Streptococcus spp.                 0.448  0.149 Inf
##  (NASM non-chromogenes) / (Streptococcus-like org.)          2.989  1.366 Inf
##  (NASM non-chromogenes) / Other                              1.431  0.472 Inf
##  Streptococcus spp. / (Streptococcus-like org.)              6.669  3.044 Inf
##  Streptococcus spp. / Other                                  3.192  1.032 Inf
##  (Streptococcus-like org.) / Other                           0.479  0.203 Inf
##  asymp.LCL asymp.UCL
##      0.979      9.01
##      6.588     74.55
##      3.042     32.43
##     14.765    297.20
##      9.716    103.45
##      3.138     17.73
##      1.425      7.85
##      6.602     75.34
##      4.788     23.79
##      0.174      1.16
##      0.813     10.99
##      0.559      3.66
##      1.816     24.49
##      1.270      8.02
##      0.143      1.60
## 
## Results are averaged over the levels of: Site 
## Confidence level used: 0.95 
## Conf-level adjustment: tukey method for comparing a family of 6 estimates 
## Intervals are back-transformed from the log odds ratio scale
# Adjusted risk by farm and multiple comparisons

emmeans(Model1,pairwise~Site,type = "response")
## $emmeans
##  Site  prob     SE  df asymp.LCL asymp.UCL
##  A    0.421 0.0656 Inf     0.300     0.552
##  B    0.371 0.0461 Inf     0.286     0.465
##  C    0.341 0.0458 Inf     0.257     0.435
##  D    0.459 0.1017 Inf     0.275     0.654
##  E    0.524 0.0651 Inf     0.397     0.647
## 
## Results are averaged over the levels of: bacterial.group 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale 
## 
## $contrasts
##  contrast odds.ratio    SE  df null z.ratio p.value
##  A / B         1.229 0.412 Inf    1   0.617  0.9724
##  A / C         1.407 0.479 Inf    1   1.001  0.8549
##  A / D         0.858 0.404 Inf    1  -0.325  0.9976
##  A / E         0.660 0.230 Inf    1  -1.191  0.7567
##  B / C         1.144 0.302 Inf    1   0.511  0.9863
##  B / D         0.698 0.325 Inf    1  -0.773  0.9385
##  B / E         0.537 0.183 Inf    1  -1.822  0.3609
##  C / D         0.610 0.285 Inf    1  -1.059  0.8274
##  C / E         0.469 0.167 Inf    1  -2.132  0.2066
##  D / E         0.770 0.354 Inf    1  -0.569  0.9795
## 
## Results are averaged over the levels of: bacterial.group 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## Tests are performed on the log odds ratio scale
# ICC

performance::icc(Model1,by_group = TRUE)
## Warning: Can't compute random effect variances. Some variance components equal
##   zero. Your model may suffer from singularity (see `?lme4::isSingular`
##   and `?performance::check_singularity`).
##   Solution: Respecify random structure! You may also decrease the
##   `tolerance` level to enforce the calculation of random effect variances.
## [1] NA
# Visualization Adj risks and 95%CI

# Extract information from emmeans

Adj.plot <- emmeans(Model1,pairwise~bacterial.group,type = "response")

# Create dataframe to use in ggplot

Adj.plot.emmeans <- as.data.frame(summary(Adj.plot)$emmean)
Adj.plot.emmeans
##  bacterial.group                 prob         SE  df asymp.LCL asymp.UCL
##  Staphylococcus aureus      0.8849219 0.03668296 Inf 0.7914803 0.9396816
##  Staphylococcus chromogenes 0.7213556 0.03877999 Inf 0.6394673 0.7907312
##  NASM non-chromogenes       0.2576111 0.04666799 Inf 0.1770170 0.3588971
##  Streptococcus spp.         0.4363559 0.06021807 Inf 0.3239143 0.5557451
##  Streptococcus-like org.    0.1040094 0.03585974 Inf 0.0517774 0.1979333
##  Other                      0.1952050 0.03574321 Inf 0.1344159 0.2747592
## 
## Results are averaged over the levels of: Site 
## Confidence level used: 0.95 
## Intervals are back-transformed from the logit scale
# Recode levels of factor for visualization

Adj.plot.emmeans$bacterial.group  <- fct_recode(Adj.plot.emmeans$bacterial.group,"SAU" = "Staphylococcus aureus","SCH" = "Staphylococcus chromogenes","Strep" = "Streptococcus spp.","SLO" = "Streptococcus-like org.","NASM non-SCH" = "NASM non-chromogenes")

# Create plot

ggplot(Adj.plot.emmeans, aes(x = bacterial.group , y = (prob)*100 ,fill = bacterial.group)) +
  geom_bar(position = "dodge", stat = "identity") + ylab("% Persistent (95%CI)") + 
  xlab("Bacterial group") +
  geom_errorbar(aes(ymin =  (asymp.LCL)*100, ymax =  (asymp.UCL)*100), width = 0.2) +
  theme_classic() +
  theme(legend.position = "none") +
  labs(fill = "Bacterial group") +
  scale_fill_viridis(discrete = TRUE) +
  ylim(0,100)

ggsave(plot = last_plot(),"./figures/plot_persistence_model.png",width = 20, height = 10, units = "cm")