Investigation of intramammary infections in primiparous cows during early lactation on organic dairy farms
Felipe Pena Mosca
6/16/2023
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 59prop.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 104prop.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 25prop.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 46prop.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 24prop.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 12prop.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 39prop.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.1372.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 ' ' 1emmeans(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 scale2.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 scaleemmeans_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")