Analysis log - Evaluation of a soluble milk LPS/IgG ELISA for identification of Gram-negative mastitis pathogens

Author

Gabriella Harb, Sam Rowe, Ruth Zadoks

Published

July 27, 2026

Study Objectives

  1. Evaluate the diagnostic performance of the soluble milk lipopolysaccharide/immunoglobulin G (sLPS/IgG) complex ELISA in differentiating mastitis cases associated with Gram-negative (GN) pathogens from non-Gram-negative (non-GN) cases, and compare its diagnostic performance to that of N-acetyl-β-D-glucosaminidase (NAGase) and lactate dehydrogenase (LDH).
  2. Investigate how microbial culture classification and intramammary inflammation relate to sLPS/IgG complexes by evaluating associations between pathogen groups, sLPS/IgG concentrations, and established inflammatory markers (NAGase and LDH).

Import data and packages

Code 
library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)
library(summarytools)
library(pROC)
library(tidyverse)
library(lme4)
library(lmerTest)
library(performance)
library(brms)
library(patchwork)
library(loo)
library(tidybayes)
library(rstatix)
library(knitr)
library(kableExtra)
library(broom.mixed)

koru_df <- read_excel("20260528.Koru.df.xlsx", sheet = "combined")

Objective 1

For this analysis, all ‘Gram positive’, ‘Contaminated’, and ‘No Growth’ samples will be grouped as non-GN.

1.1 Data exploration

Code 
koru_df_1 <- koru_df %>%
  mutate(pathogen = as.character(Diagnosis)) %>%
  filter(!is.na(pathogen) & pathogen != "N/A")
Code 
print(dfSummary(koru_df_1[, 1:8],
            valid.col = FALSE,
            graph.magnif = 0.8,
            style = "multiline"),
  method = "render"
)

Data Frame Summary

koru_df_1

Dimensions: 622 x 8
Duplicates: 0
No Variable Stats / Values Freqs (% of Valid) Graph Missing
1 Date [character]
1. 2022
2. 2024
3. Unknown
322 ( 51.8% )
299 ( 48.1% )
1 ( 0.2% )
0 (0.0%)
2 Sample number [numeric]
Mean (sd) : 311.5 (179.7)
min ≤ med ≤ max:
1 ≤ 311.5 ≤ 622
IQR (CV) : 310.5 (0.6)
622 distinct values 0 (0.0%)
3 sLPS_IgG [numeric]
Mean (sd) : 0.1 (0.2)
min ≤ med ≤ max:
0 ≤ 0 ≤ 1.7
IQR (CV) : 0 (2.6)
241 distinct values 0 (0.0%)
4 LDH_activity [numeric]
Mean (sd) : 0.7 (1)
min ≤ med ≤ max:
0 ≤ 0.4 ≤ 6.9
IQR (CV) : 1 (1.4)
499 distinct values 0 (0.0%)
5 NAGase_activity [numeric]
Mean (sd) : 4448.4 (5166.2)
min ≤ med ≤ max:
0 ≤ 2658.9 ≤ 32853.4
IQR (CV) : 5285.7 (1.2)
588 distinct values 0 (0.0%)
6 Diagnosis_old [character]
1. Contamination
2. E. coli
3. Klebsiella pneumoniae
4. No growth
5. Pasteurella multocida
6. Serratia liquefaciens
7. Serratia marcescens
8. Staphylococcus aureus
9. Staphylococcus sp. (non a
10. Streptococcus dysgalactia
11. Streptococcus uberis
54 ( 13.2% )
52 ( 12.7% )
79 ( 19.3% )
51 ( 12.4% )
5 ( 1.2% )
1 ( 0.2% )
4 ( 1.0% )
11 ( 2.7% )
20 ( 4.9% )
20 ( 4.9% )
113 ( 27.6% )
212 (34.1%)
7 Diagnosis [character]
1. Contamination
2. E. coli
3. Klebsiella pneumoniae
4. No growth
5. Pasteurella multocida
6. Serratia liquefaciens
7. Serratia marcescens
8. Staphylococcus aureus
9. Staphylococcus sp. (non a
10. Streptococcus dysgalactia
11. Streptococcus uberis
89 ( 14.3% )
52 ( 8.4% )
79 ( 12.7% )
228 ( 36.7% )
5 ( 0.8% )
1 ( 0.2% )
4 ( 0.6% )
11 ( 1.8% )
20 ( 3.2% )
20 ( 3.2% )
113 ( 18.2% )
0 (0.0%)
8 Gram_status [character]
1. Gram negative
2. Gram positive
3. N/A
141 ( 22.7% )
164 ( 26.4% )
317 ( 51.0% )
0 (0.0%)

Generated by summarytools 1.1.5 (R version 4.6.0)
2026-07-27

Code 
koru_df_1$pathogen <- factor(
  koru_df_1$pathogen,
  levels = c(
    "No growth",
    "Contamination",
    "E. coli",
    "Klebsiella pneumoniae",
    "Pasteurella multocida",
    "Serratia marcescens",
    "Serratia liquefaciens",
    "Staphylococcus aureus",
    "Staphylococcus sp. (non aureus)",
    "Streptococcus uberis",
    "Streptococcus dysgalactiae"
  )
)
Code 
koru_df_gram_group <- koru_df_1 %>%
  filter(!is.na(`Gram_status`)) %>%  
  mutate(
    GN_status = ifelse(`Gram_status` == "Gram negative", 1, 0),
    GN_label  = ifelse(GN_status == 1, "Gram-negative", "non-Gram-negative")
  )

There is a greater number of cases in the non-Gram-negative category (77.4%), particularly due to a large number of No Growth samples (36.7%) (Table 1). Pasteurella and Serratia groups have <10 observations each so, where relevant, they will be combined as ‘Other Gram-negative (GN)’ in plots and removed from statistical tests.

Code 
pathogen_summary <- koru_df_1 %>%
  group_by(pathogen) %>%
  summarise(
    n = n(),
    median_LDH = median(LDH_activity, na.rm = TRUE),
    Gram_status = first(Gram_status),
    IQR_LDH = IQR(LDH_activity, na.rm = TRUE),
    median_NAG = median(NAGase_activity, na.rm = TRUE),
    IQR_NAG = IQR(NAGase_activity, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(desc(median_LDH)) %>%
  mutate(
    Percent = round(100 * n / sum(n), 1),
    `n (%)` = paste0(n, " (", Percent, "%)"), 
    `Median LDH activity (IQR)` =    
      paste0(round(median_LDH, 2), " (", round(IQR_LDH, 2), ")"),  
    `Median NAGase activity (IQR)` =    
      paste0(round(median_NAG, 0), " (", round(IQR_NAG, 0), ")")
  ) %>%
  select(
    Pathogen = pathogen,
    `n (%)`,
    `Median LDH activity (IQR)`,
    `Median NAGase activity (IQR)`,
    `Gram status` = Gram_status
  )

kable(
  pathogen_summary,
  caption = "Table 1. Median (Interquartile range (IQR)) LDH and NAGase activities by pathogen group.",
  booktabs = TRUE,
  align = c("l", "c", "c", "c")
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center"
  )
Table 1. Median (Interquartile range (IQR)) LDH and NAGase activities by pathogen group.
Pathogen n (%) Median LDH activity (IQR) Median NAGase activity (IQR) Gram status
Pasteurella multocida 5 (0.8%) 2.39 (2.04) 18748 (5665) Gram negative
Streptococcus uberis 113 (18.2%) 1.07 (1.18) 4155 (7175) Gram positive
E. coli 52 (8.4%) 0.82 (1.45) 3669 (7803) Gram negative
Klebsiella pneumoniae 79 (12.7%) 0.77 (0.78) 4358 (4544) Gram negative
Serratia marcescens 4 (0.6%) 0.6 (0.36) 3794 (1180) Gram negative
Streptococcus dysgalactiae 20 (3.2%) 0.53 (0.86) 4482 (5121) Gram positive
Serratia liquefaciens 1 (0.2%) 0.41 (0) 2456 (0) Gram negative
Staphylococcus aureus 11 (1.8%) 0.34 (1.11) 4202 (3425) Gram positive
Contamination 89 (14.3%) 0.21 (0.74) 2396 (5221) N/A
Staphylococcus sp. (non aureus) 20 (3.2%) 0.15 (0.32) 1618 (2628) Gram positive
No growth 228 (36.7%) 0.07 (0.29) 1370 (3620) N/A
Code 
valid_pathogens <- koru_df_1 %>%
  group_by(pathogen) %>%
  summarise(n = n()) %>%
  filter(n >= 10) %>%
  pull(pathogen)

koru_df_1_filtered <- koru_df_1 %>%
  filter(pathogen %in% valid_pathogens)

The majority of the data sits near zero but with a large number of high value outliers for both groups. Gram-negative samples exhibit moderately higher sLPS/IgG values than non-Gram-negative samples (Wilcoxon rank-sum test, p < 0.001; effect size r = 0.38; Figure 1).

Code 
ggplot(
  koru_df_gram_group, aes(x = GN_label, y = sLPS_IgG, fill = GN_label)
) +
  geom_boxplot(
    width = 0.6,
    alpha = 0.8,
    outlier.shape = NA
  ) +
  geom_jitter(
    width = 0.15,
    alpha = 0.3,
    size = 1,
    colour = "grey30"
  ) +
  scale_fill_manual(
    values = c(
      "Gram-negative" = "grey60",
      "non-Gram-negative" = "white"
    )
  ) +
  labs(
    x = "Sample gram status",
    y = "sLPS/IgG (OD)"
  ) +
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.grid.major = element_line(colour = "grey85"),
    panel.grid.minor = element_line(colour = "grey92"),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 11)
  )

Figure 1. Distribution of sLPS/IgG ELISA results in Gram-negative and non-Gram-negative samples. Boxes represent the median and interquartile range, whiskers extend to 1.5 × the interquartile range, and points represent individual observations.

Code 
# Difference in distributions?
wilcox.test(sLPS_IgG ~ GN_label, data = koru_df_gram_group)

    Wilcoxon rank sum test with continuity correction

data:  sLPS_IgG by GN_label
W = 51924, p-value < 2.2e-16
alternative hypothesis: true location shift is not equal to 0
Code 
# Magnitude of difference between groups (calculates r and uses conventional Cohen thresholds where 0.3-0.5 = moderate effect)
koru_df_gram_group %>%
  wilcox_effsize(sLPS_IgG ~ GN_label,
    ci = TRUE,
    conf.level = 0.95,
    nboot = 2000
  )
# A tibble: 1 × 9
  .y.      group1        group2 effsize    n1    n2 conf.low conf.high magnitude
* <chr>    <chr>         <chr>    <dbl> <int> <int>    <dbl>     <dbl> <ord>    
1 sLPS_IgG Gram-negative non-G…   0.385   141   481     0.32      0.45 moderate

1.2 ROC Analysis - sLPS/IgG ELISA

Code 
roc_obj <- roc(koru_df_gram_group$GN_status, koru_df_gram_group$sLPS_IgG)

roc_df <- data.frame(specificity = roc_obj$specificities, sensitivity = roc_obj$sensitivities)

best <- coords(
  roc_obj,
  x = "best",
  best.method = "youden",
  ret = c("threshold", "sensitivity", "specificity"),
  transpose = FALSE
)

best_thresh <- best$threshold
best_sens   <- best$sensitivity
best_spec   <- best$specificity

ggplot(roc_df, aes(x = 1 - specificity, y = sensitivity)) +
  geom_line(
    colour = "blue",
    linewidth = 1.2
  ) +
  geom_abline(
    linetype = "dashed",
    colour = "grey60"
  ) +
  geom_point(
    aes(
      x = 1 - best_spec,
      y = best_sens
    ),
    colour = "red",
    size = 3
  ) +
  labs(
    x = "False positive rate",
    y = "True positive rate",
    subtitle = paste(
      "AUC =", round(auc(roc_obj), 3),
      "| Cut-point =", round(best_thresh, 4),
      "| Se =", round(100 * best_sens, 0), "%",
      "| Sp =", round(100 * best_spec, 0), "%"
    )
  ) +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(colour = "grey85"),
    panel.grid.minor = element_line(colour = "grey92")
  )

Figure 2. Receiver Operating Characteristic (ROC) curve describing the ability of the sLPS/IgG ELISA to discriminate Gram-negative from non-Gram-negative mastitis cases. The curve plots the True Positive Rate against the False Positive Rate at various probability thresholds. The dashed diagonal line represents a baseline classifier performing no better than random chance. The Area Under the Curve (AUC) achieved by the model is 0.766, demonstrating a moderate degree of discriminatory performance for detecting Gram-negative cases. The optimal cut-point determined using Youden’s index was 0.0247, corresponding to a sensitivity of 75% and specificity of 68%.

Code 
ci.auc(roc_obj)
95% CI: 0.7207-0.8106 (DeLong)

Moving the cut-point from 0.03 to 0.0247 (the optimal cut-point) gains about 13 percentage points of sensitivity while losing only about 9 percentage points of specificity, resulting in a higher Youden index (0.427) (calculated as Se+Sp-1) (Table 2). The optimal cut-points are chosen at highest Youden index. Increasing efficiency regions means that sensitivity rises quickly relative to the false positives. When efficiency = 0, false positives are increasing with no gain in detecting true positives. Efficiency at different cut-points for this test is highly variable, which reflects the skewed data distribution.

Code 
thresholds <- unique(c(
  best$threshold,
  seq(0, 0.45, by = 0.015)
))

roc_mid <- coords(
  roc_obj,
  x = thresholds,
  input = "threshold",
  ret = c("threshold", "sensitivity", "specificity"),
  transpose = FALSE
) %>%
  mutate(
    FPR = 1 - specificity,
    Youden = sensitivity + specificity - 1
  ) %>%
  arrange(threshold)

roc_mid <- roc_mid %>%
  arrange(FPR) %>%
  mutate(
    delta_sens = sensitivity - lag(sensitivity),
    delta_fpr  = FPR - lag(FPR),
    efficiency = delta_sens / delta_fpr
  )
Code 
roc_table <- roc_mid %>%
  distinct(threshold, .keep_all = TRUE) %>%
  mutate(
    Threshold = round(threshold, 4),
    Sensitivity = paste0(round(100 * sensitivity, 0), "%"),
    Specificity = paste0(round(100 * specificity, 0), "%"),
    `False positive rate` = paste0(round(100 * FPR, 0), "%"),
    `Youden index` = round(Youden, 3),
    Efficiency = round(efficiency, 2)
  ) %>%
  select(
    Threshold,
    Sensitivity,
    Specificity,
    `False positive rate`,
    `Youden index`,
    Efficiency
  )

opt_row <- which(abs(roc_table$Threshold - 0.0247) < 0.0001)

kable(
  roc_table,
  caption = "Table 2. Diagnostic performance at selected sLPS/IgG thresholds for discrimination of Gram-negative and non-Gram-negative cases. The optimal cut-point is highlighted in grey and bolded. For efficiency: NA = no previous cut-point available for comparison; NaN = no change between successive cut-points; -Inf = performance changed in an unfavourable direction.",
  booktabs = TRUE,
  align = "cccccc"
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center"
  ) %>%
  row_spec(
    opt_row,
    bold = TRUE,
    background = "grey90"
  )
Table 2. Diagnostic performance at selected sLPS/IgG thresholds for discrimination of Gram-negative and non-Gram-negative cases. The optimal cut-point is highlighted in grey and bolded. For efficiency: NA = no previous cut-point available for comparison; NaN = no change between successive cut-points; -Inf = performance changed in an unfavourable direction.
Threshold Sensitivity Specificity False positive rate Youden index Efficiency
0.4350 12% 98% 2% 0.096 NA
0.4500 12% 98% 2% 0.096 NaN
0.4050 13% 97% 3% 0.101 3.41
0.4200 13% 97% 3% 0.101 NaN
0.3750 13% 97% 3% 0.106 3.41
0.3900 13% 97% 3% 0.106 NaN
0.3300 16% 97% 3% 0.130 6.82
0.3450 14% 97% 3% 0.109 -Inf
0.3600 13% 97% 3% 0.101 -Inf
0.3000 17% 96% 4% 0.135 17.06
0.3150 17% 96% 4% 0.135 NaN
0.2850 17% 96% 4% 0.133 0.00
0.2550 18% 96% 4% 0.138 3.41
0.2700 18% 96% 4% 0.138 NaN
0.2400 19% 96% 4% 0.150 6.82
0.2250 19% 95% 5% 0.146 0.00
0.1650 23% 95% 5% 0.182 6.82
0.1800 23% 95% 5% 0.175 -Inf
0.1950 21% 95% 5% 0.161 -Inf
0.2100 20% 95% 5% 0.147 -Inf
0.1500 24% 95% 5% 0.187 20.47
0.1200 28% 94% 6% 0.225 10.23
0.1350 28% 94% 6% 0.218 -Inf
0.1050 30% 94% 6% 0.238 10.23
0.0900 34% 94% 6% 0.278 20.47
0.0750 35% 93% 7% 0.286 2.27
0.0600 40% 92% 8% 0.325 4.78
0.0450 45% 88% 12% 0.331 1.14
0.0300 61% 77% 23% 0.377 1.42
0.0247 75% 68% 32% 0.434 1.66
0.0150 93% 27% 73% 0.201 0.43
0.0000 100% 0% 100% 0.000 0.26

1.3 Distribution of sLPS/IgG concentrations relative to the ROC-derived diagnostic threshold

The following density plot is used to visualise the overlap responsible for the observed false-positive and false-negative classifications. Gram-negative cases exhibited somewhat higher sLPS/IgG concentrations than non-Gram-negative cases. However, there is considerable overlap between the distributions (Figure 3). The ROC-derived optimal cut-point occurred within this overlap region, so a proportion of Gram-negative and non-Gram-negative cases could not be fully distinguished on the basis of sLPS/IgG concentration alone. This overlap is consistent with the moderate discriminatory ability observed in ROC analysis (AUC = 0.766).

Code 
ggplot(
  koru_df_gram_group,
  aes(
    x = sLPS_IgG,
    fill = GN_label
  )
) +
  geom_density(
    alpha = 0.4
  ) +
  geom_vline(
    xintercept = best_thresh,
    linetype = "dashed",
    colour = "black",
    linewidth = 0.8
  ) +
  scale_fill_manual(
    values = c(
      "Gram-negative" = "grey60",
      "non-Gram-negative" = "white"
    )
  ) + 
  scale_x_log10() +
  labs(
    x = "sLPS/IgG OD (log tranformed axis)",
    y = "Density", 
    fill = "Gram status"
  ) +
  theme_bw()

Figure 3.Distribution of sLPS/IgG optical density (OD) values among samples with Gram-negative and non-Gram negative culture infections. Density plots are shown on a log-transformed x-axis. The dashed, vertical line indicates the optimal sLPS/IgG threshold used to differentiate between groups.

1.4 NAGase/LDH activity vs infection type + ROC curve comparisons

This section aims to determine whether the ELISA can provide information beyond inflammatory markers already being explored for mastitis identification (NAGase and LDH).

Code 
ggplot(
  koru_df_gram_group, aes(x = GN_label, y = NAGase_activity, fill = GN_label)
) +
  geom_boxplot(
    width = 0.6,
    alpha = 0.8,
    outlier.shape = NA
  ) +
  geom_jitter(
    width = 0.15,
    alpha = 0.3,
    size = 1,
    colour = "grey30"
  ) +
  scale_fill_manual(
    values = c(
      "Gram-negative" = "grey60",
      "non-Gram-negative" = "white"
    )
  ) +
  labs(
    x = "Sample gram status",
    y = "NAGase activity (ΔOD/min/uL)"
  ) +
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.grid.major = element_line(colour = "grey85"),
    panel.grid.minor = element_line(colour = "grey92"),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 11)
  )

Figure 4. Distribution of NAGase activity results in Gram-negative and non-Gram-negative samples. Boxes represent the median and interquartile range, whiskers extend to 1.5 × the interquartile range, and points represent individual observations.

Code 
ggplot(
  koru_df_gram_group, aes(x = GN_label, y = LDH_activity, fill = GN_label)
) +
  geom_boxplot(
    width = 0.6,
    alpha = 0.8,
    outlier.shape = NA
  ) +
  geom_jitter(
    width = 0.15,
    alpha = 0.3,
    size = 1,
    colour = "grey30"
  ) +
  scale_fill_manual(
    values = c(
      "Gram-negative" = "grey60",
      "non-Gram-negative" = "white"
    )
  ) +
  labs(
    x = "Sample gram status",
    y = "LDH activity (µmoles/min/L))"
  ) +
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.grid.major = element_line(colour = "grey85"),
    panel.grid.minor = element_line(colour = "grey92"),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 11)
  )

Figure 5. Distribution of LDH activity results in Gram-negative and non-Gram-negative samples. Boxes represent the median and interquartile range, whiskers extend to 1.5 × the interquartile range, and points represent individual observations.

Code 
# Extract optimal values
elisa_best <- coords(
  roc_obj_grouped,
  x = "best",
  best.method = "youden",
  ret = c("threshold", "sensitivity", "specificity"),
  transpose = FALSE
)

ldh_best <- coords(
  roc_ldh_grouped,
  x = "best",
  best.method = "youden",
  ret = c("threshold", "sensitivity", "specificity"),
  transpose = FALSE
)

nag_best <- coords(
  roc_nag_grouped,
  x = "best",
  best.method = "youden",
  ret = c("threshold", "sensitivity", "specificity"),
  transpose = FALSE
)

ggplot() +
  geom_line(
    data = data.frame(
      fpr = 1 - roc_obj_grouped$specificities,
      tpr = roc_obj_grouped$sensitivities
    ),
    aes(fpr, tpr),
    colour = "blue",
    linewidth = 1.2
  ) +

  geom_line(
    data = data.frame(
      fpr = 1 - roc_ldh_grouped$specificities,
      tpr = roc_ldh_grouped$sensitivities
    ),
    aes(fpr, tpr),
    colour = "green3",
    linewidth = 1.2
  ) +

  geom_line(
    data = data.frame(
      fpr = 1 - roc_nag_grouped$specificities,
      tpr = roc_nag_grouped$sensitivities
    ),
    aes(fpr, tpr),
    colour = "red",
    linewidth = 1.2
  ) +

  geom_abline(
    linetype = "dashed",
    colour = "grey60"
  ) +

  labs(
    x = "False positive rate",
    y = "True positive rate"
  ) +

  annotate(
    "text",
    x = 0.68,
    y = 0.30,
    hjust = 0,
    size = 3.5,
    label = paste0(
      "ELISA: AUC = ", round(auc(roc_obj_grouped), 3), "\n",
      "LDH: AUC = ", round(auc(roc_ldh_grouped), 3), "\n",
      "NAGase: AUC = ", round(auc(roc_nag_grouped), 3)
    )
  ) +

  theme_minimal() +
  theme(
    panel.grid.major = element_line(colour = "grey85"),
    panel.grid.minor = element_line(colour = "grey92")
  )

Figure 6. Receiver Operating Characteristic (ROC) curves describing the ability of the sLPS/IgG ELISA, LDH activity, and NAGase activity to discriminate Gram-negative from non-Gram-negative mastitis cases. The curve plots the True Positive Rate against the False Positive Rate at various probability thresholds. The dashed diagonal line represents a baseline classifier performing no better than random chance. The Area Under the Curve (AUC) achieved by LDH is 0.708 and by NAGase is 0.642, demonstrating a limited degree of discriminatory performance for detecting Gram-negative cases.

Code 
knitr::kable(
  roc_summary,
  caption = "Table 3. Optimal cut-points and diagnostic performance for ELISA, LDH and NAGase."
)
Table 3. Optimal cut-points and diagnostic performance for ELISA, LDH and NAGase.
Test AUC Optimal cut-point Sensitivity (%) Specificity (%)
ELISA 0.766 0.0247 75.2 68.2
LDH 0.708 0.3300 80.1 57.0
NAGase 0.642 1515.3407 83.0 42.8
Code 
roc.test(roc_obj_grouped, roc_ldh_grouped)

    DeLong's test for two correlated ROC curves

data:  roc_obj_grouped and roc_ldh_grouped
Z = 2.2153, p-value = 0.02674
alternative hypothesis: true difference in AUC is not equal to 0
95 percent confidence interval:
 0.006599103 0.107908203
sample estimates:
AUC of roc1 AUC of roc2 
  0.7656109   0.7083573
Code 
roc.test(roc_obj_grouped, roc_nag_grouped)

    DeLong's test for two correlated ROC curves

data:  roc_obj_grouped and roc_nag_grouped
Z = 5.2606, p-value = 1.436e-07
alternative hypothesis: true difference in AUC is not equal to 0
95 percent confidence interval:
 0.07748821 0.16951493
sample estimates:
AUC of roc1 AUC of roc2 
  0.7656109   0.6421094

The ELISA is significantly better at discriminating GN from non-GN cases than LDH activity (AUC = 0.708; p = 0.02674) and NAGase activity (AUC = 0.642; p < 0.001). Although performance exceeded that of LDH and NAGase, discrimination remains insufficient for use as a standalone diagnostic test.

1.5 Sensitivity analysis (Gram-status)

Instead of grouping all non-Gram-negative cases (i.e Gram positive + Contaminated + No Growth), this section tests whether excluding ‘Contaminated’ and ‘No Growth’ samples from the non-GN category changes the outcome (i.e. comparing only Gram-positive (GP) vs Gram-negative).

Code 
# Exclude N/As 
koru_df_gram_excl <- koru_df %>%
  filter(!is.na(`Gram_status`) & `Gram_status` != "N/A") %>%
  mutate(GN_status = ifelse(`Gram_status` == "Gram negative", 1, 0),
         GN_label = ifelse(GN_status == 1, "Gram-negative", "Gram-positive"))

GN samples exhibited moderately higher sLPS/IgG values than GP samples (Wilcoxon rank-sum test, p < 0.001; effect size r = 0.32). The majority of the data sit near zero but with a large number of high value outliers for both groups (Figure 7).

Code 
ggplot(
  koru_df_gram_excl, aes(x = GN_label, y = sLPS_IgG, fill = GN_label)
) +
  geom_boxplot(
    width = 0.6,
    alpha = 0.8,
    outlier.shape = NA
  ) +
  geom_jitter(
    width = 0.15,
    alpha = 0.3,
    size = 1,
    colour = "grey30"
  ) +
  scale_fill_manual(
    values = c(
      "Gram-negative" = "grey60",
      "non-Gram-negative" = "white"
    )
  ) +
  labs(
    x = "Sample gram status",
    y = "sLPS/IgG (OD)"
  ) +
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.grid.major = element_line(colour = "grey85"),
    panel.grid.minor = element_line(colour = "grey92"),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 11)
  )

Figure 7. Distribution of sLPS/IgG ELISA results in Gram-negative and Gram-positive samples. Boxes represent the median and interquartile range, whiskers extend to 1.5 × the interquartile range, and points represent individual observations.

Code 
# Difference in distributions?
wilcox.test(sLPS_IgG ~ GN_label, data = koru_df_gram_excl)

    Wilcoxon rank sum test with continuity correction

data:  sLPS_IgG by GN_label
W = 15804, p-value = 3.32e-08
alternative hypothesis: true location shift is not equal to 0
Code 
# Magnitude of difference between groups
koru_df_gram_excl %>%
  wilcox_effsize(sLPS_IgG ~ GN_label, 
  ci = TRUE,
    conf.level = 0.95,
    nboot = 2000
  )
# A tibble: 1 × 9
  .y.      group1        group2 effsize    n1    n2 conf.low conf.high magnitude
* <chr>    <chr>         <chr>    <dbl> <int> <int>    <dbl>     <dbl> <ord>    
1 sLPS_IgG Gram-negative Gram-…   0.316   141   164     0.21      0.42 moderate
Code 
roc_obj2 <- roc(koru_df_gram_excl$GN_status, koru_df_gram_excl$sLPS_IgG)

roc_df2 <- data.frame(specificity = roc_obj2$specificities, sensitivity = roc_obj2$sensitivities)

best2 <- coords(
  roc_obj2,
  x = "best",
  best.method = "youden",
  ret = c("threshold", "sensitivity", "specificity"),
  transpose = FALSE
)

best_thresh2 <- best2$threshold
best_sens2   <- best2$sensitivity
best_spec2   <- best2$specificity

ggplot(roc_df2, aes(x = 1 - specificity, y = sensitivity)) +
  geom_line(
    colour = "blue",
    linewidth = 1.2
  ) +
  geom_abline(
    linetype = "dashed",
    colour = "grey60"
  ) +
  geom_point(
    aes(
      x = 1 - best_spec2,
      y = best_sens2
    ),
    colour = "red",
    size = 3
  ) +
  labs(
    x = "False positive rate",
    y = "True positive rate",
    subtitle = paste(
      "AUC =", round(auc(roc_obj2), 3),
      "| Cut-point =", round(best_thresh2, 4),
      "| Se =", round(100 * best_sens2, 1), "%",
      "| Sp =", round(100 * best_spec2, 1), "%"
    )
  ) +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(colour = "grey85"),
    panel.grid.minor = element_line(colour = "grey92")
  )

Figure 8. Receiver Operating Characteristic (ROC) curve describing the ability of the sLPS/IgG ELISA to discriminate Gram-negative from GP mastitis cases. The curve plots the True Positive Rate against the False Positive Rate at various probability thresholds. The dashed diagonal line represents a baseline classifier performing no better than random chance. The Area Under the Curve (AUC) achieved by the model is 0.683, demonstrating a limited degree of discriminatory performance for detecting Gram-negative cases. The optimal cut-point determined using Youden’s index was 0.0247, corresponding to a sensitivity of 75% and specificity of 55%.

Code 
roc_ldh_gram <- roc(koru_df_gram_excl$GN_status, koru_df_gram_excl$LDH_activity)
roc_nag_gram <- roc(koru_df_gram_excl$GN_status, koru_df_gram_excl$NAGase_activity)

auc(roc_obj2)
Area under the curve: 0.6834
Code 
auc(roc_ldh_gram)
Area under the curve: 0.4855
Code 
auc(roc_nag_gram)
Area under the curve: 0.5525
Code 
roc.test(roc_obj2, roc_ldh_gram)

    Bootstrap test for two correlated ROC curves

data:  roc_obj2 and roc_ldh_gram
D = 3.8523, boot.n = 2000, boot.stratified = 1, p-value = 0.000117
alternative hypothesis: true difference in AUC is not equal to 0
sample estimates:
AUC of roc1 AUC of roc2 
  0.6834458   0.4854696
Code 
roc.test(roc_obj2, roc_nag_gram)

    DeLong's test for two correlated ROC curves

data:  roc_obj2 and roc_nag_gram
Z = 4.2528, p-value = 2.111e-05
alternative hypothesis: true difference in AUC is not equal to 0
95 percent confidence interval:
 0.07058635 0.19126281
sample estimates:
AUC of roc1 AUC of roc2 
  0.6834458   0.5525212

When the non-GN group is limited to GP samples, the ELISA is significantly better at discriminating GN from non-GN cases than LDH activity (AUC = 0.4855; p < 0.001) and NAGase activity (AUC = 0.5525; p < 0.001).

1.6 Sensitivity analysis (Sample type - enrolment only)

This section tests whether excluding follow-up samples changes the performance (i.e. only using enrolment samples).

Code 
koru_enrol <- koru_df %>%
  filter(Sample_Type == "Enrolment")
Code 
koru_enrol_gram <- koru_enrol %>%
  filter(!is.na(Gram_status)) %>%
  mutate(
    GN_status = ifelse(Gram_status == "Gram negative", 1, 0),
    GN_label  = ifelse(GN_status == 1, "Gram-negative", "non-Gram-negative")
  )
Code 
koru_enrol_gram %>%
  mutate(
    Gram_status = ifelse(
      Gram_status == "N/A",
      "Non-Gram-negative",
      Gram_status
    )
  ) %>%
  count(Diagnosis, Gram_status) %>%
  mutate(
    Percent = round(100 * n / sum(n), 1),
    `n (%)` = paste0(n, " (", Percent, "%)")
  ) %>%
  select(
    'Culture result' = Diagnosis,
    `Gram status` = Gram_status,
    `n (%)`
  ) %>%
  kable(caption = "Table 4. Gram status of each culture result and the number (n) of each result obtained within the reduced dataset (enrolment samples only).") %>%
  kable_styling(position = "left")
Table 4. Gram status of each culture result and the number (n) of each result obtained within the reduced dataset (enrolment samples only).
Culture result Gram status n (%)
Contamination Non-Gram-negative 73 (15.5%)
E. coli Gram negative 52 (11%)
Klebsiella pneumoniae Gram negative 79 (16.7%)
No growth Non-Gram-negative 103 (21.8%)
Pasteurella multocida Gram negative 5 (1.1%)
Serratia liquefaciens Gram negative 1 (0.2%)
Serratia marcescens Gram negative 4 (0.8%)
Staphylococcus aureus Gram positive 11 (2.3%)
Staphylococcus sp. (non aureus) Gram positive 20 (4.2%)
Streptococcus dysgalactiae Gram positive 20 (4.2%)
Streptococcus uberis Gram positive 104 (22%)

Gram-negative enrolment samples exhibited moderately higher sLPS/IgG values than non-Gram-negative enrolment samples (Wilcoxon rank-sum test, p < 0.001; effect size r = 0.31; Figure 9). The majority of the data sits near zero but with a large number of high value outliers for both groups.

Code 
ggplot(
  koru_enrol_gram, aes(x = GN_label, y = sLPS_IgG, fill = GN_label)
) +
  geom_boxplot(
    width = 0.6,
    alpha = 0.8,
    outlier.shape = NA
  ) +
  geom_jitter(
    width = 0.15,
    alpha = 0.3,
    size = 1,
    colour = "grey30"
  ) +
  scale_fill_manual(
    values = c(
      "Gram-negative" = "grey60",
      "non-Gram-negative" = "white"
    )
  ) +
  labs(
    x = "Sample gram status",
    y = "sLPS/IgG (OD)"
  ) +
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.grid.major = element_line(colour = "grey85"),
    panel.grid.minor = element_line(colour = "grey92"),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 11)
  )

Figure 9. Distribution of sLPS/IgG ELISA results in Gram-negative and non-GN enrolment samples. Boxes represent the median and interquartile range, whiskers extend to 1.5 × the interquartile range, and points represent individual observations.

Code 
# Difference in distributions?
wilcox.test(sLPS_IgG ~ GN_label, data = koru_enrol_gram)

    Wilcoxon rank sum test with continuity correction

data:  sLPS_IgG by GN_label
W = 32604, p-value = 8.273e-12
alternative hypothesis: true location shift is not equal to 0
Code 
# Magnitude of difference between groups 
koru_enrol_gram %>%
  wilcox_effsize(sLPS_IgG ~ GN_label)
# A tibble: 1 × 7
  .y.      group1        group2            effsize    n1    n2 magnitude
* <chr>    <chr>         <chr>               <dbl> <int> <int> <ord>    
1 sLPS_IgG Gram-negative non-Gram-negative   0.315   141   331 moderate
Code 
roc_enrol <- roc(koru_enrol_gram$GN_status, koru_enrol_gram$sLPS_IgG)
ci.auc(roc_enrol)
95% CI: 0.6449-0.7523 (DeLong)
Code 
roc_obj_enrol <- roc(
  koru_enrol_gram$GN_status,
  koru_enrol_gram$sLPS_IgG
)

roc_df_enrol <- data.frame(
  specificity = roc_obj_enrol$specificities,
  sensitivity = roc_obj_enrol$sensitivities
)

best_enrol <- coords(
  roc_obj_enrol,
  x = "best",
  best.method = "youden",
  ret = c("threshold", "sensitivity", "specificity"),
  transpose = FALSE
)

best_thresh_enrol <- best_enrol$threshold
best_sens_enrol   <- best_enrol$sensitivity
best_spec_enrol   <- best_enrol$specificity

ggplot(roc_df_enrol, aes(x = 1 - specificity, y = sensitivity)) +
  geom_line(
    colour = "blue",
    linewidth = 1.2
  ) +
  geom_abline(
    linetype = "dashed",
    colour = "grey60"
  ) +
  geom_point(
    aes(
      x = 1 - best_spec_enrol,
      y = best_sens_enrol
    ),
    colour = "red",
    size = 3
  ) +
  labs(
    x = "False positive rate",
    y = "True positive rate",
    subtitle = paste(
      "AUC =", round(auc(roc_obj_enrol), 3),
      "| Cut-point =", round(best_thresh_enrol, 4),
      "| Se =", round(100 * best_sens_enrol, 1), "%",
      "| Sp =", round(100 * best_spec_enrol, 1), "%"
    )
  ) +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(colour = "grey85"),
    panel.grid.minor = element_line(colour = "grey92")
  )

Figure 10. Receiver Operating Characteristic (ROC) curve describing the ability of the sLPS/IgG ELISA to discriminate Gram-negative from GP mastitis cases in enrolment samples only. The curve plots the True Positive Rate against the False Positive Rate at various probability thresholds. The dashed diagonal line represents a baseline classifier performing no better than random chance. The Area Under the Curve (AUC) achieved by the model is 0.699, demonstrating a limited degree of discriminatory performance for detecting Gram-negative cases. The optimal cut-point determined using Youden’s index was 0.0247, corresponding to a sensitivity of 75% and specificity of 58%.

Code 
roc_ldh_enrol <- roc(koru_enrol_gram$GN_status, koru_enrol_gram$LDH_activity)
roc_nag_enrol <- roc(koru_enrol_gram$GN_status, koru_enrol_gram$NAGase_activity)

auc(roc_enrol)
Area under the curve: 0.6986
Code 
auc(roc_ldh_enrol)
Area under the curve: 0.6064
Code 
auc(roc_nag_enrol)
Area under the curve: 0.5887
Code 
roc.test(roc_enrol, roc_ldh_enrol)

    DeLong's test for two correlated ROC curves

data:  roc_enrol and roc_ldh_enrol
Z = 3.0317, p-value = 0.002432
alternative hypothesis: true difference in AUC is not equal to 0
95 percent confidence interval:
 0.03259751 0.15182108
sample estimates:
AUC of roc1 AUC of roc2 
  0.6986030   0.6063937
Code 
roc.test(roc_enrol, roc_nag_enrol)

    DeLong's test for two correlated ROC curves

data:  roc_enrol and roc_nag_enrol
Z = 4.0823, p-value = 4.458e-05
alternative hypothesis: true difference in AUC is not equal to 0
95 percent confidence interval:
 0.05712904 0.16264341
sample estimates:
AUC of roc1 AUC of roc2 
  0.6986030   0.5887168

When the dataset is limited to enrolment samples, the ELISA is significantly better at discriminating GN from non-GN cases than LDH activity (AUC = 0.6064; p = 0.02432) and NAGase activity (AUC = 0.5887; p < 0.001).

1.7 Sensitivity analysis (Sample type - follow up only)

This section tests whether excluding enrolment samples changes the performance (i.e. only using follow up samples).

Code 
koru_follow <- koru_df %>%
  filter(Sample_Type == "Follow-up")
Code 
koru_follow_gram <- koru_follow %>%
  filter(!is.na(Gram_status)) %>%
  mutate(
    GN_status = ifelse(Gram_status == "Gram negative", 1, 0),
    GN_label  = ifelse(GN_status == 1, "Gram-negative", "non-Gram-negative")
  )
Code 
ggplot(
  koru_follow_gram, aes(x = GN_label, y = sLPS_IgG, fill = GN_label)
) +
  geom_boxplot(
    width = 0.6,
    alpha = 0.8,
    outlier.shape = NA
  ) +
  geom_jitter(
    width = 0.15,
    alpha = 0.3,
    size = 1,
    colour = "grey30"
  ) +
  scale_fill_manual(
    values = c(
      "Gram-negative" = "grey60",
      "non-Gram-negative" = "white"
    )
  ) +
  labs(
    x = "Sample gram status",
    y = "sLPS/IgG (OD)"
  ) +
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.grid.major = element_line(colour = "grey85"),
    panel.grid.minor = element_line(colour = "grey92"),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 11)
  )

Figure. 11. Distribution of sLPS/IgG ELISA results in Gram-negative and non-GN follow-up samples. Boxes represent the median and interquartile range, whiskers extend to 1.5 × the interquartile range, and points represent individual observations.

There are no GN test results in follow up samples. Performance cannot be compared between GN and non-GN. This is expected given follow-up samples are mostly non-clinical. Table 5 shows a breakdown of culture results within follow up samples.

Code 
koru_follow_gram %>%
  mutate(
    Gram_status = ifelse(
      Gram_status == "N/A",
      "Non-Gram-negative",
      Gram_status
    )
  ) %>%
  count(Diagnosis, Gram_status) %>%
  mutate(
    Percent = round(100 * n / sum(n), 1),
    `n (%)` = paste0(n, " (", Percent, "%)")
  ) %>%
  select(
    'Culture result' = Diagnosis,
    `Gram status` = Gram_status,
    `n (%)`
  ) %>%
  kable(caption = "Table 5. Gram status of each culture result and the number (n) of each result obtained within the reduced dataset (follow-up samples only)") %>%
  kable_styling(position = "left")
Table 5. Gram status of each culture result and the number (n) of each result obtained within the reduced dataset (follow-up samples only)
Culture result Gram status n (%)
Contamination Non-Gram-negative 16 (10.7%)
No growth Non-Gram-negative 125 (83.3%)
Streptococcus uberis Gram positive 9 (6%)

1.8 Objective 1 conclusions

Code 
obj1_summary <- tribble(
  ~Metric, ~`All samples`, ~`Clinical samples`,
  "Overall mean sLPS/IgG", "0.081", "0.102",
  "Overall median sLPS/IgG", "0.023", "0.026",
  "GN median sLPS/IgG", "0.040", "0.040",
  "Non-GN median sLPS/IgG", "0.021", "0.023",
  "Group difference (95% CI)", "0.38 (0.31-0.45)", "0.31 (0.21-0.41)",
  "P-value (distribution)", "<0.001", "<0.001",
  "AUC", "0.766", "0.699",
  "Optimal cut-point", "0.0247", "0.0247",
  "Se @ optimal cut-point", "75%", "75%",
  "Sp @ optimal cut-point", "68%", "58%"
)

kable(
  obj1_summary,
  caption = "Table 6. Comparison of sLPS/IgG and ROC-derived diagnostic performance metrics for identification of Gram-negative mastitis cases using all samples and clinical (enrolment) samples only.",
  align = c("l", "c", "c", "c")
)
Table 6. Comparison of sLPS/IgG and ROC-derived diagnostic performance metrics for identification of Gram-negative mastitis cases using all samples and clinical (enrolment) samples only.
Metric All samples Clinical samples
Overall mean sLPS/IgG 0.081 0.102
Overall median sLPS/IgG 0.023 0.026
GN median sLPS/IgG 0.040 0.040
Non-GN median sLPS/IgG 0.021 0.023
Group difference (95% CI) 0.38 (0.31-0.45) 0.31 (0.21-0.41)
P-value (distribution) <0.001 <0.001
AUC 0.766 0.699
Optimal cut-point 0.0247 0.0247
Se @ optimal cut-point 75% 75%
Sp @ optimal cut-point 68% 58%

The ELISA has moderate accuracy in discriminating between GN and non-GN cases of clinical mastitis when all samples are included in analysis (AUC = 0.766). In practice, this is unlikely to be acceptable as a standalone diagnostic test. The ELISA performs better when grouping all non-GN together compared to GN vs GP only and when retaining follow-up samples in the data set. Grouping non-GN and retaining non-clinical follow-up samples likely improve the ELISAs ability to distinguish between groups because of the wider spread in low values contributed by Contamination/No Growth categories. Grouping all non-GN results better reflect samples obtained in practice as Contamination/No Growth outcomes are common in standard laboratory milk culturing. The ELISA also performed better than LDH and NAGase in discriminating between GN and non-GN cases of clinical mastitis. While these two bio markers are not yet widely used in mastitis diagnostics, there has been research into their usefulness in identifying inflammation during mastitis. These results indicate that further work is warranted into how these three bio markers (sLPS/IgG as measured by the ELISA, NAGase, and LDH) may be used appropriately.

Rowe et al. (2024) discussed whether the performance of evaluated point-of-care tests in distinguishing GP from non-GP cases was sufficient to support treatment decisions. Caution is required when making these comparisons, as the current study assessed a different diagnostic target (GN vs non-GN cases). Despite this limitation and the moderate discriminatory ability of the ELISA, its sensitivity (75%) was comparable to several of the currently available point-of-care tests. However, specificity was lower (68%) than most published alternatives (Se/Sp: 0.84/0.75 for Biplate, 0.76/0.90 for Accumast, 0.89/0.79 for Check-Up, 0.67/0.83 for Petrifilm, and 0.55/0.81 for Mastatest). While the ELISA would correctly identify a reasonable proportion of true Gram-negative cases, the lower specificity indicates a greater number of false-positive results. In this context, false positives represent non-GN cases incorrectly classified as GN. As the purpose of the test is to aid treatment selection, this could result in approximately one third of the Gram-positive cases being incorrectly considered less likely to benefit from antimicrobial therapy, potentially reducing appropriate treatment allocation. Reducing antimicrobial use overall is beneficial in a ONE Health context, but can also reduce cow health and welfare when its application is warranted.

Objective 2

2.1 Scatterplots by Gram group

Code 
df_pathogen_plot <- koru_df_1 %>%
  mutate(
    pathogen_plot = case_when(
      pathogen %in% c(
        "Pasteurella multocida",
        "Serratia marcescens",
        "Serratia liquefaciens"
      ) ~ "Other GN",
      TRUE ~ as.character(pathogen)
    )
  )

df_pathogen_plot$pathogen_plot <- factor(
  df_pathogen_plot$pathogen_plot,
  levels = c(
    "No growth",
    "Contamination",
    "E. coli",
    "Klebsiella pneumoniae",
    "Other GN",
    "Staphylococcus aureus",
    "Staphylococcus sp. (non aureus)",
    "Streptococcus uberis",
    "Streptococcus dysgalactiae"
  )
)

Positive associations between sLPS/IgG concentration and both NAGase and LDH activity were observed within Gram-negative, Gram-positive, no-growth and contamination groups (Figure 12: Figure 13). Although sLPS/IgG concentrations were somewhat higher among Gram-negative cases, the persistence of these associations within non-Gram-negative groups suggests that sLPS/IgG may reflect not only exposure to Gram-negative bacterial LPS but also the severity of the host inflammatory response. NAGase and LDH appear to be correlated with each other to a similar level, regardless of Gram status (Figure 14).

Code 
df_pathogen_plot2 <- df_pathogen_plot %>%
  mutate(
    pathogen_group = case_when(
      Diagnosis == "No growth" ~ "No growth",
      Diagnosis == "Contamination" ~ "Contamination",
      Gram_status == "Gram negative" ~ "Gram-negative",
      Gram_status == "Gram positive" ~ "Gram-positive"
    )
  )

df_pathogen_plot2$pathogen_group <- factor(
  df_pathogen_plot2$pathogen_group,
  levels = c(
    "Gram-negative",
    "Gram-positive",
    "No growth",
    "Contamination"
  )
)
Code 
ggplot(
  df_pathogen_plot2,
  aes(
    x = NAGase_activity,
    y = sLPS_IgG
  )
) +
  geom_point(
    colour = "grey50",
    alpha = 0.5,
    size = 1
  ) +
  geom_smooth(
    method = "lm",
    colour = "black",
    se = FALSE,
    linewidth = 0.8
  ) +
  scale_x_log10() +
  scale_y_log10() +
  facet_wrap(
    ~ pathogen_group,
    ncol = 2
  ) +
  labs(
    x = "NAGase activity",
    y = "sLPS/IgG (OD)"
  ) +

  theme_bw() +
  theme(
    strip.background = element_blank(),
    strip.text = element_text(face = "bold")
  )

Figure 12. Relationship between NAGase activity and sLPS/IgG optical density (OD) within Gram groups. Points represent individual samples and solid lines represent fitted trends. Both axes are presented on the log10 scale to aid visualisation.

Code 
ggplot(
  df_pathogen_plot2,
  aes(
    x = LDH_activity,
    y = sLPS_IgG
  )
) +
  geom_point(
    colour = "grey50",
    alpha = 0.5,
    size = 1
  ) +
  geom_smooth(
    method = "lm",
    colour = "black",
    se = FALSE,
    linewidth = 0.8
  ) +
  scale_x_log10() +
  scale_y_log10() +
  facet_wrap(
    ~ pathogen_group,
    ncol = 2
  ) +
  labs(
    x = "LDH activity",
    y = "sLPS/IgG (OD)"
  ) +

  theme_bw() +
  theme(
    strip.background = element_blank(),
    strip.text = element_text(face = "bold")
  )

Figure 13. Relationship between LDH activity and sLPS/IgG optical density (OD) within Gram groups. Points represent individual samples and solid lines represent fitted trends. Both axes are presented on the log10 scale to aid visualisation.

Code 
ggplot(
  df_pathogen_plot2,
  aes(
    x = LDH_activity,
    y = NAGase_activity
  )
) +
  geom_point(
    colour = "grey50",
    alpha = 0.5,
    size = 1
  ) +
  geom_smooth(
    method = "lm",
    colour = "black",
    se = FALSE,
    linewidth = 0.8
  ) +
  scale_x_log10() +
  scale_y_log10() +
  facet_wrap(
    ~ pathogen_group,
    ncol = 2
  ) +
  labs(
    x = "LDH activity",
    y = "NAGase activity"
  ) +

  theme_bw() +
  theme(
    strip.background = element_blank(),
    strip.text = element_text(face = "bold")
  )

Figure 14. Relationship between NAGase activity and LDH within Gram groups. Points represent individual samples and solid lines represent fitted trends. Both axes are presented on the log10 scale to aid visualisation.

2.2 Scatterplots by pathogen group

There are positive associations between LDH activity, NAGase activity, and sLPS/IgG optical density across most pathogen groups (Figure 15; Figure 16; Figure 17). As expected, due to the LPS in their outer cell membranes, these relationships were evident across all GN pathogens. However, they were also observed for the GP pathogens, S. aureus and the Streptococcus spp. and the No growth and Contamination categories. Relationships were weaker for No growth and Contamination, although the presence of a slight positive relationship may suggest that poor sample quality prevented culturing of the causative pathogen, the infection was associated with low bacterial numbers or the inflammatory response was still evident even after bacterial clearance. The flat regression line for non-aureus Staphylococcus species may suggest that these cases are generally less severe. The consistency of these positive trends indicate that sLPS/IgG concentration reflects not only pathogen classification but also the severity of the host inflammatory response. There appears to be a generic inflammatory response to these major pathogens. Localised inflammation in the mammary gland may trigger a systemic response causing leaky gut (pro inflammatory cytokines released into the bloodstream degrade tight junctions), leading to the detection of LPS originating from the GN bacteria in the gut microbiome rather than the mastitis pathogen itself (references? other explanations? come back to this…).

NAGase and LDH appear to be correlated with each other to a similar level, regardless of pathogen (Figure 17). Reflections of SCC? (references? other explanations? come back to this…).

Code 
df_pathogen_plot$pathogen_plot <- factor(
  df_pathogen_plot$pathogen_plot,
  levels = c(
    # Gram-negative pathogens
    "E. coli",
    "Klebsiella pneumoniae",
    "Other GN",

    # Major Gram-positive pathogens
    "Streptococcus uberis",
    "Streptococcus dysgalactiae",
    "Staphylococcus aureus",

    # Non-pathogen/minor categories
    "No growth",
    "Contamination",
    "Staphylococcus sp. (non aureus)"
  )
)
ggplot(
  df_pathogen_plot,
  aes(
    x = NAGase_activity,
    y = sLPS_IgG
  )
) +
  geom_point(
    colour = "grey50",
    alpha = 0.5,
    size = 1
  ) +
  geom_smooth(
    method = "lm",
    colour = "black",
    se = FALSE,
    linewidth = 0.8
  ) +
  scale_x_log10() +
  scale_y_log10() +
  facet_wrap(
    ~ pathogen_plot,
    ncol = 3
  ) +
  labs(
    x = "NAGase activity",
    y = "sLPS/IgG OD"
  ) +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(face = "bold")
  )

Figure 15. Relationship between NAGase activity and sLPS/IgG optical density (OD) within pathogen groups. Points represent individual samples and solid lines represent fitted trends. Both axes are presented on the log10 scale to aid visualisation. Rare Gram-negative pathogens are combined in the ‘Other GN’ category (Pasteurella multocida and Serratia spp.).

Code 
df_pathogen_plot$pathogen_plot <- factor(
  df_pathogen_plot$pathogen_plot,
  levels = c(
    # Gram-negative pathogens
    "E. coli",
    "Klebsiella pneumoniae",
    "Other GN",

    # Major Gram-positive pathogens
    "Streptococcus uberis",
    "Streptococcus dysgalactiae",
    "Staphylococcus aureus",

    # Non-pathogen/minor categories
    "No growth",
    "Contamination",
    "Staphylococcus sp. (non aureus)"
  )
)
ggplot(
  df_pathogen_plot,
  aes(
    x = LDH_activity,
    y = sLPS_IgG
  )
) +
  geom_point(
    colour = "grey50",
    alpha = 0.5,
    size = 1
  ) +
  geom_smooth(
    method = "lm",
    colour = "black",
    se = FALSE,
    linewidth = 0.8
  ) +
  scale_x_log10() +
  scale_y_log10() +
  facet_wrap(
    ~ pathogen_plot,
    ncol = 3
  ) +
  labs(
    x = "LDH activity",
    y = "sLPS/IgG OD"
  ) +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(face = "bold")
  )

Figure 16. Relationship between LDH activity and sLPS/IgG optical density (OD) within pathogen groups. Points represent individual samples and solid lines represent fitted trends. Both axes are presented on the log10 scale to aid visualisation. Rare Gram-negative pathogens are combined in the ‘Other GN’ category (Pasteurella multocida and Serratia spp.).

Code 
df_pathogen_plot$pathogen_plot <- factor(
  df_pathogen_plot$pathogen_plot,
  levels = c(
    # Gram-negative pathogens
    "E. coli",
    "Klebsiella pneumoniae",
    "Other GN",

    # Major Gram-positive pathogens
    "Streptococcus uberis",
    "Streptococcus dysgalactiae",
    "Staphylococcus aureus",

    # Non-pathogen/minor categories
    "No growth",
    "Contamination",
    "Staphylococcus sp. (non aureus)"
  )
)
ggplot(
  df_pathogen_plot,
  aes(
    x = LDH_activity,
    y = NAGase_activity
  )
) +
  geom_point(
    colour = "grey50",
    alpha = 0.5,
    size = 1
  ) +
  geom_smooth(
    method = "lm",
    colour = "black",
    se = FALSE,
    linewidth = 0.8
  ) +
  scale_x_log10() +
  scale_y_log10() +
  facet_wrap(
    ~ pathogen_plot,
    ncol = 3
  ) +
  labs(
    x = "LDG activity",
    y = "NAGase activity"
  ) +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(face = "bold")
  )

Figure 17. Relationship between NAGase activity and LDH activity within pathogen groups. Points represent individual samples and solid lines represent fitted trends. Both axes are presented on the log10 scale to aid visualisation. Rare Gram-negative pathogens are combined in the ‘Other GN’ category (Pasteurella multocida and Serratia spp.).

2.3 Scatterplots by year

Plots in section 2.2 show some evidence of clustering at high sLPS/IgG values. To explore whether the high values are associated with severe clinical mastitis, samples will be plotted by pathogen per year. All severe clinical mastitis cases in this dataset appeared in 2022 and were associated with K. pneuomiae and S. uberis.

Samples culture-positive for Gram-negative pathogens, particularly E. coli and K. pneumoniae, tended to exhibit higher sLPS/IgG OD values at increasing NAGase activity, suggesting a positive relationship between inflammatory activity and anti-LPS antibody response. In contrast, no-growth, contamination, and most Gram-positive pathogen groups generally showed lower sLPS/IgG OD values across the observed range of NAGase activity, although a small number of high-value outliers were present. Similar patterns were observed across collection years.

Code 
koru_df <- koru_df %>%
  mutate(
    Date = na_if(Date, "???")
  )
Code 
table(koru_df$Date, useNA = "ifany")

   2022    2024 Unknown 
    322     299       1
Code 
df_pathogen_plot <- koru_df %>%
  mutate(
    pathogen_plot = case_when(
      Diagnosis %in% c(
        "Pasteurella multocida",
        "Serratia marcescens",
        "Serratia liquefaciens"
      ) ~ "Other GN",
      TRUE ~ as.character(Diagnosis)
    )
  )
Code 
ggplot(
  df_pathogen_plot,
  aes(
    x = NAGase_activity,
    y = sLPS_IgG,
    colour = factor(Date)
  )
) +
  geom_point(
    alpha = 0.6,
    size = 2
  ) +
  scale_x_log10() +
  scale_y_log10() +
  facet_wrap(~ pathogen_plot, ncol = 3) +
  labs(
    x = "NAGase activity",
    y = "sLPS/IgG OD",
    colour = "Collection year"
  ) +
  theme_bw()

Figure 18. Relationship between NAGase activity and sLPS/IgG optical density (OD) across pathogen groups. Each point represents an individual milk sample, coloured by collection year. NAGase activity and sLPS/IgG OD are shown on log-transformed scales.

Code 
ggplot(
  df_pathogen_plot,
  aes(
    x = LDH_activity,
    y = sLPS_IgG,
    colour = factor(Date)
  )
) +
  geom_point(
    alpha = 0.6,
    size = 2
  ) +
  scale_x_log10() +
  scale_y_log10() +
  facet_wrap(~ pathogen_plot, ncol = 3) +
  labs(
    x = "LDH activity",
    y = "sLPS/IgG OD",
    colour = "Collection year"
  ) +
  theme_bw()

Figure 19. Relationship between LDH activity and sLPS/IgG optical density (OD) across pathogen groups. Each point represents an individual milk sample, coloured by collection year. LDH activity and sLPS/IgG OD are shown on log-transformed scales.

Code 
table(koru_df$Date, koru_df$Sample_Type, useNA = "ifany")
         
          Enrolment Follow-up
  2022          322         0
  2024          150       149
  Unknown         0         1
Code 
koru_enrol_year <- koru_df %>%
  filter(Sample_Type == "Enrolment")
Code 
df_year_compare_enrol <- koru_enrol_year %>%
  filter(
    Diagnosis %in% c(
      "Klebsiella pneumoniae",
      "Streptococcus uberis"
    )
  )
Code 
ggplot(
  df_year_compare_enrol,
  aes(
    x = NAGase_activity,
    y = sLPS_IgG,
    colour = factor(Date)
  )
) +
  geom_point(size = 2, alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE) +
  scale_x_log10() +
  scale_y_log10() +
  facet_wrap(~ Diagnosis) +
  labs(
    x = "NAGase activity",
    y = "sLPS/IgG OD",
    colour = "Collection year"
  ) +
  theme_bw()

Figure 20. Relationship between NAGase activity and sLPS/IgG optical density (OD) for Klebsiella pneumoniae and Streptococcus uberis. Each point represents an individual milk sample, coloured by collection year. NAGase activity and sLPS/IgG OD are shown on log-transformed scales.

Code 
ggplot(
  df_year_compare_enrol,
  aes(
    x = LDH_activity,
    y = sLPS_IgG,
    colour = factor(Date)
  )
) +
  geom_point(size = 2, alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE) +
  scale_x_log10() +
  scale_y_log10() +
  facet_wrap(~ Diagnosis) +
  labs(
    x = "LDH activity",
    y = "sLPS/IgG OD",
    colour = "Collection year"
  ) +
  theme_bw()

Figure 21. Relationship between LDH activity and sLPS/IgG optical density (OD) for Klebsiella pneumoniae and Streptococcus uberis. Each point represents an individual milk sample, coloured by collection year. LDH activity and sLPS/IgG OD are shown on log-transformed scales.

To explore whether the high values are associated with severe clinical mastitis, samples will be plotted by pathogen per year. All severe clinical mastitis cases in this dataset appeared in 2022 and were associated with K. pneuomiae and S. uberis.

While exact severity classification of cases is unavailable, high biomarker values do appear to be associated with the more severe clinical mastitis cases sampled in 2022 (Figures 18-21).

2.4 Boxplots by pathogen group

Contamination and no growth categories appear to have lowest inflammatory marker activities and sLPS/IgG, GN pathogens tend to have slightly higher median levels, and levels vary greatly within the same pathogen group (Figure 22; Figure 23; Figure 24).

Code 
ggplot(
  koru_df_1_filtered,
  aes(
    x = pathogen,
    y = sLPS_IgG,
    fill = Gram_status
  )
) +
  geom_boxplot(
    outlier.alpha = 0.3,
    colour = "grey30"
  ) +
  geom_hline(
    yintercept = best_thresh,
    linetype = "dashed",
    colour = "red",
    linewidth = 0.8
  ) +
  scale_fill_manual(
    values = c(
      "Gram positive" = "white",
      "Gram negative" = "grey60",
      "N/A" = "grey85"
    )
  ) +
  scale_y_log10() +
  labs(
    x = "Pathogen",
    y = "sLPS/IgG (log scale)",
    fill = "Gram status"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "right",
    panel.grid.major = element_line(colour = "grey85"),
    panel.grid.minor = element_line(colour = "grey92")
  )

Figure 22. Distribution of sLPS/IgG optical density among pathogen groups isolated from mastitis samples. Boxplots display the median, interquartile range, and observations within 1.5 × the interquartile range, with individual outliers shown as points. sLPS/IgG is presented on a log10 scale to improve visualisation of the highly skewed distribution. Boxes are shaded according to Gram status. The dashed horizontal line indicates the optimal sLPS/IgG cut-point identified from ROC analysis for discriminating Gram-negative from non-Gram-negative infections (0.0247).

Code 
ggplot(
  koru_df_1_filtered,
  aes(
    x = pathogen,
    y = NAGase_activity,
    fill = Gram_status
  )
) +
  geom_boxplot(
    outlier.alpha = 0.3,
    colour = "grey30"
    ) +
  geom_hline(
    yintercept = nag_best$threshold,
    linetype = "dashed",
    colour = "red",
    linewidth = 0.8
    ) +
  scale_fill_manual(
    values = c(
      "Gram positive" = "white",
      "Gram negative" = "grey60",
      "N/A" = "grey85"
    )
  ) +
  scale_y_log10() +
  labs(
    x = "Pathogen",
    y = "NAGase activity (log scale)", 
    fill = "Gram status"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "right",
    panel.grid.major = element_line(colour = "grey85"),
    panel.grid.minor = element_line(colour = "grey92")
  )

Figure 23. Distribution of NAGase activity among pathogen groups isolated from mastitis samples. Boxplots display the median, interquartile range, and observations within 1.5 × the interquartile range, with individual outliers shown as points. NAGase activity is presented on a log10 scale to improve visualisation of the highly skewed distribution. Boxes are shaded according to Gram status. The dashed horizontal line indicates the optimal NAGase cut-point identified from ROC analysis for discriminating Gram-negative from non-Gram-negative infections (1515.34).

Code 
ggplot(
  koru_df_1_filtered,
  aes(
    x = pathogen,
    y = LDH_activity,
    fill = Gram_status
  )
) +
  geom_boxplot(
    outlier.alpha = 0.3,
    colour = "grey30"
  ) +
  geom_hline(
    yintercept = ldh_best$threshold,
    linetype = "dashed",
    colour = "red",
    linewidth = 0.8
    ) +
  scale_fill_manual(
    values = c(
      "Gram positive" = "white",
      "Gram negative" = "grey60",
      "N/A" = "grey85"
    )
  ) +
  scale_y_log10() +
  labs(
    x = "Pathogen",
    y = "LDH activity (log scale)", 
    fill = "Gram status"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "right",
    panel.grid.major = element_line(colour = "grey85"),
    panel.grid.minor = element_line(colour = "grey92")
  )

Figure 24. Distribution of LDH among pathogen groups isolated from mastitis samples. Boxplots display the median, interquartile range, and observations within 1.5 × the interquartile range, with individual outliers shown as points. LDH activity is presented on a log10 scale to improve visualisation of the highly skewed distribution. Boxes are shaded according to Gram status. The dashed horizontal line indicates the optimal LDH cut-point identified from ROC analysis for discriminating Gram-negative from non-Gram-negative infections (0.3300).

2.5 Statistical comparisons - pathogen groups within biomarkers

Code 
kruskal.test(sLPS_IgG ~ pathogen, data = koru_df_1_filtered)

    Kruskal-Wallis rank sum test

data:  sLPS_IgG by pathogen
Kruskal-Wallis chi-squared = 133.33, df = 7, p-value < 2.2e-16
Code 
kruskal.test(NAGase_activity ~ pathogen, data = koru_df_1_filtered)

    Kruskal-Wallis rank sum test

data:  NAGase_activity by pathogen
Kruskal-Wallis chi-squared = 58.176, df = 7, p-value = 3.488e-10
Code 
kruskal.test(LDH_activity ~ pathogen, data = koru_df_1_filtered)

    Kruskal-Wallis rank sum test

data:  LDH_activity by pathogen
Kruskal-Wallis chi-squared = 193, df = 7, p-value < 2.2e-16
  • Estimated biomarker levels do differ across pathogen groups [sLPS/IgG OD (p < 0.001); NAGase activity (p < 0.001); LDH activity (p < 0.001)].
Code 
pairwise.wilcox.test(
  koru_df_1_filtered$sLPS_IgG,
  koru_df_1_filtered$pathogen,
  p.adjust.method = "BH"
)

    Pairwise comparisons using Wilcoxon rank sum test with continuity correction 

data:  koru_df_1_filtered$sLPS_IgG and koru_df_1_filtered$pathogen 

                                No growth Contamination E. coli
Contamination                   6.6e-06   -             -      
E. coli                         < 2e-16   5.7e-07       -      
Klebsiella pneumoniae           6.2e-15   9.9e-05       0.10847
Staphylococcus aureus           0.35148   0.26280       0.00092
Staphylococcus sp. (non aureus) 0.11351   0.24628       7.5e-06
Streptococcus uberis            4.8e-07   0.39689       4.6e-05
Streptococcus dysgalactiae      0.00029   0.24944       0.02522
                                Klebsiella pneumoniae Staphylococcus aureus
Contamination                   -                     -                    
E. coli                         -                     -                    
Klebsiella pneumoniae           -                     -                    
Staphylococcus aureus           0.02657               -                    
Staphylococcus sp. (non aureus) 0.00059               0.86303              
Streptococcus uberis            0.00982               0.50555              
Streptococcus dysgalactiae      0.26280               0.10665              
                                Staphylococcus sp. (non aureus)
Contamination                   -                              
E. coli                         -                              
Klebsiella pneumoniae           -                              
Staphylococcus aureus           -                              
Staphylococcus sp. (non aureus) -                              
Streptococcus uberis            0.24628                        
Streptococcus dysgalactiae      0.04993                        
                                Streptococcus uberis
Contamination                   -                   
E. coli                         -                   
Klebsiella pneumoniae           -                   
Staphylococcus aureus           -                   
Staphylococcus sp. (non aureus) -                   
Streptococcus uberis            -                   
Streptococcus dysgalactiae      0.42207             

P value adjustment method: BH
  • Compared to No growth, all pathogens other than Staphylococcus species (S. aureus (p = 0.3511); non aureus (p = 0.1135)) result in significantly higher inflammation when measured by sLPS/IgG OD (p < 0.05). There is no significant difference in inflammation between pathogens in the GN group (p = 0.10847). Within the GP group, most pathogens do not differ significantly in inflammation level. There is weak evidence that Staphylococcus spp. (non aureus) and Streptococcus dysgalactiae differ in NAGase levels (p = 0.04993) (where inflammation associated with S. dysgalactiae tends to be higher).
Code 
pairwise.wilcox.test(
  koru_df_1_filtered$NAGase_activity,
  koru_df_1_filtered$pathogen,
  p.adjust.method = "BH"
)

    Pairwise comparisons using Wilcoxon rank sum test with continuity correction 

data:  koru_df_1_filtered$NAGase_activity and koru_df_1_filtered$pathogen 

                                No growth Contamination E. coli
Contamination                   0.0271    -             -      
E. coli                         2.9e-05   0.0722        -      
Klebsiella pneumoniae           1.5e-06   0.0577        0.9264 
Staphylococcus aureus           0.0474    0.4570        0.9146 
Staphylococcus sp. (non aureus) 0.6570    0.6061        0.0474 
Streptococcus uberis            1.9e-05   0.1274        0.6374 
Streptococcus dysgalactiae      0.0019    0.0974        0.8733 
                                Klebsiella pneumoniae Staphylococcus aureus
Contamination                   -                     -                    
E. coli                         -                     -                    
Klebsiella pneumoniae           -                     -                    
Staphylococcus aureus           0.9146                -                    
Staphylococcus sp. (non aureus) 0.0235                0.0974               
Streptococcus uberis            0.8167                0.9720               
Streptococcus dysgalactiae      0.8079                0.8079               
                                Staphylococcus sp. (non aureus)
Contamination                   -                              
E. coli                         -                              
Klebsiella pneumoniae           -                              
Staphylococcus aureus           -                              
Staphylococcus sp. (non aureus) -                              
Streptococcus uberis            0.0974                         
Streptococcus dysgalactiae      0.0271                         
                                Streptococcus uberis
Contamination                   -                   
E. coli                         -                   
Klebsiella pneumoniae           -                   
Staphylococcus aureus           -                   
Staphylococcus sp. (non aureus) -                   
Streptococcus uberis            -                   
Streptococcus dysgalactiae      0.6594              

P value adjustment method: BH
  • Compared to No growth, most pathogens result in significantly higher inflammation when measured by NAGase activity (p < 0.05). There is no significant difference in inflammation between No Growth and Staphylococcus spp. (non aureus) (p = 0.6570). There is no significant difference in inflammation between pathogens in the GN group (p = 0.9264). Within the GP group, most pathogens do not differ significantly in inflammation level. The only GP pathogens that differ significantly are Staphylococcus spp. (non aureus) and Streptococcus dysgalactiae (p = 0.0271) (where inflammation with S. dysgalactiae is higher).
Code 
pairwise.wilcox.test(
  koru_df_1_filtered$LDH_activity,
  koru_df_1_filtered$pathogen,
  p.adjust.method = "BH"
)

    Pairwise comparisons using Wilcoxon rank sum test with continuity correction 

data:  koru_df_1_filtered$LDH_activity and koru_df_1_filtered$pathogen 

                                No growth Contamination E. coli
Contamination                   0.0064    -             -      
E. coli                         4.3e-12   3.4e-05       -      
Klebsiella pneumoniae           < 2e-16   3.0e-06       0.5429 
Staphylococcus aureus           0.0269    0.3382        0.2174 
Staphylococcus sp. (non aureus) 0.0442    0.9781        0.0013 
Streptococcus uberis            < 2e-16   1.1e-10       0.4967 
Streptococcus dysgalactiae      8.8e-05   0.0440        0.2969 
                                Klebsiella pneumoniae Staphylococcus aureus
Contamination                   -                     -                    
E. coli                         -                     -                    
Klebsiella pneumoniae           -                     -                    
Staphylococcus aureus           0.2567                -                    
Staphylococcus sp. (non aureus) 0.0001                0.3623               
Streptococcus uberis            0.0440                0.0811               
Streptococcus dysgalactiae      0.3641                0.6339               
                                Staphylococcus sp. (non aureus)
Contamination                   -                              
E. coli                         -                              
Klebsiella pneumoniae           -                              
Staphylococcus aureus           -                              
Staphylococcus sp. (non aureus) -                              
Streptococcus uberis            6.0e-06                        
Streptococcus dysgalactiae      0.0338                         
                                Streptococcus uberis
Contamination                   -                   
E. coli                         -                   
Klebsiella pneumoniae           -                   
Staphylococcus aureus           -                   
Staphylococcus sp. (non aureus) -                   
Streptococcus uberis            -                   
Streptococcus dysgalactiae      0.0484              

P value adjustment method: BH
  • Compared to No growth, all pathogens result in significantly higher inflammation when measured by LDH activity (p < 0.05). There is weak evidence that there is a difference in inflammation between pathogens in the GN group (p = 0.5429). Within the GP group, most pathogens do not differ significantly in inflammation level. The only GP pathogens that differ significantly are Staphylococcus spp. (non aureus) and S. ubris (p < 0.001) (where inflammation with S. ubris is higher) and Streptococcus uberis and Streptococcus dysgalactiae (p = 0.0484) (where inflammation with S. dysgalactiae is higher).

2.6 Distribution of raw data

From previous sections, it is clear that there are some differences in inflammatory response between pathogen groups but biomarkers overlap in distibution substantially. Therefore, diagnostic performance is moderate. Here we attempt to understand why this may occur by investigating the relationships between biomarkers.

Code 
h1 <- ggplot(koru_df, aes(x = LDH_activity)) + geom_histogram(bins = 30)
h2 <- ggplot(koru_df, aes(x = NAGase_activity)) + geom_histogram(bins = 30)
h3 <- ggplot(koru_df, aes(x = sLPS_IgG)) + geom_histogram(bins = 30)

h1 + h2 + h3

  • LDH activity: strong right skew, clear outliers.
  • NAGase activity: strong right skew, clear outliers.
  • sLPS/IgG: values concentrated near 0.
Code 
df_model <- koru_df %>%
  mutate(
    log_LDH = log1p(LDH_activity),
    log_NAG = log1p(NAGase_activity),
    log_LPS = log(sLPS_IgG)
  )

2.8 Group pathogens for linear models

Pathogen was identified as a possible effect modifier/confounder (COME BACK TO THIS) so will be added as a covariate in certain models to assess performance.

Code 
df_model_pathgroup <- df_model %>%
  mutate(
    diagnosis_group = case_when(
      Diagnosis %in% c("Staphylococcus aureus", "Staphylococcus sp. (non aureus)") ~ "Staphylococcus",
      Diagnosis %in% c("Streptococcus dysgalactiae", "Streptococcus uberis") ~ "Streptococcus",
      Diagnosis %in% c("Serratia liquefaciens", "Serratia marcescens") ~ "Serratia",
      Diagnosis == "Klebsiella pneumoniae" ~ "Klebsiella",
      Diagnosis == "E. coli" ~ "E. coli",
      Diagnosis == "Pasteurella multocida" ~ "Pasteurella",
      Diagnosis == "No growth" ~ "No growth",
      Diagnosis == "Contamination" ~ "Contamination",
      TRUE ~ "Other"
    )
  )

df_model_pathgroup$diagnosis_group <- factor(df_model_pathgroup$diagnosis_group)

df_model_pathgroup$diagnosis_group <- relevel(df_model_pathgroup$diagnosis_group, ref = "E. coli")

2.9 Finding the optimal model

Q1. Does the choice of model family matter?
Q2. Does controlling for pathogen matter?
Q3. Does clustering in Cow ID and/or Farm matter?
Q4. Do log transformations make a difference to model fit?

Code 
all_models <- list(
  lm_nag = lm(sLPS_IgG ~ NAGase_activity, data = df_model),
  lm_nag_path = lm(sLPS_IgG ~ NAGase_activity + diagnosis_group, data = df_model_pathgroup),
  lmer_nag_cow = lmer(sLPS_IgG ~ NAGase_activity + (1|Cow_ID), data = df_model),
  lmer_nag_farm = lmer(sLPS_IgG ~ NAGase_activity + (1|Farm), data = df_model),
  lmer_nag_both = lmer(sLPS_IgG ~ NAGase_activity + (1|Cow_ID) + (1|Farm), data = df_model),
  lmer_nag_cowpath = lmer(sLPS_IgG ~ NAGase_activity + diagnosis_group + (1|Cow_ID), data = df_model_pathgroup),
  lmer_nag_farmpath = lmer(sLPS_IgG ~ NAGase_activity + diagnosis_group + (1|Farm), data = df_model_pathgroup),
  lmer_nag_bothpath = lmer(sLPS_IgG ~ NAGase_activity + diagnosis_group + (1|Cow_ID) + (1|Farm), data = df_model_pathgroup),
  glm_nag = glm(sLPS_IgG ~ NAGase_activity, family = Gamma(link = "log"), data = df_model),
  glm_nag_path = glm(sLPS_IgG ~ NAGase_activity + diagnosis_group, family = Gamma(link = "log"), data = df_model_pathgroup),
  glm_nag_identity = glm(sLPS_IgG ~ NAGase_activity, family = Gamma(link = "identity"), data = df_model),
  glm_nag_pathidentity = glm(sLPS_IgG ~ NAGase_activity + diagnosis_group, family = Gamma(link = "identity"), data = df_model_pathgroup),
  
  mod_nag_logY = lm(log_LPS ~ NAGase_activity, data = df_model),
  mod_nag_logXY = lm(sLPS_IgG ~ log_NAG, data = df_model),
  mod_nag_logY_path = lm(log_LPS ~ NAGase_activity + diagnosis_group, data = df_model_pathgroup),
  mod_nag_logXY_path = lm(log_LPS ~ log_NAG + diagnosis_group, data = df_model_pathgroup),
  
  
  
  lm_ldh = lm(sLPS_IgG ~ LDH_activity, data = df_model),
  lm_ldh_path = lm(sLPS_IgG ~ LDH_activity + diagnosis_group, data = df_model_pathgroup),
  lmer_ldh_cow = lmer(sLPS_IgG ~ LDH_activity + (1|Cow_ID), data = df_model),
  lmer_ldh_farm = lmer(sLPS_IgG ~ LDH_activity + (1|Farm), data = df_model),
  lmer_ldh_both = lmer(sLPS_IgG ~ LDH_activity + (1|Cow_ID) + (1|Farm), data = df_model),
  lmer_ldh_cowpath = lmer(sLPS_IgG ~ LDH_activity + diagnosis_group + (1|Cow_ID), data = df_model_pathgroup),
  lmer_ldh_farmpath = lmer(sLPS_IgG ~ LDH_activity + diagnosis_group + (1|Farm), data = df_model_pathgroup),
  lmer_ldh_bothpath = lmer(sLPS_IgG ~ LDH_activity + diagnosis_group + (1|Cow_ID) + (1|Farm), data = df_model_pathgroup),
  glm_ldh = glm(sLPS_IgG ~ LDH_activity, family = Gamma(link = "log"), data = df_model),
  glm_ldh_path = glm(sLPS_IgG ~ LDH_activity + diagnosis_group, family = Gamma(link = "log"), data = df_model_pathgroup),
  glm_ldh_identity = glm(sLPS_IgG ~ LDH_activity, family = Gamma(link = "identity"), data = df_model),
  glm_ldh_pathidentity = glm(sLPS_IgG ~ LDH_activity + diagnosis_group, family = Gamma(link = "identity"), data = df_model_pathgroup),
  
  mod_ldh_logY = lm(log_LPS ~ LDH_activity, data = df_model),
  mod_ldh_logXY = lm(sLPS_IgG ~ log_LDH, data = df_model),
  mod_ldh_logY_path = lm(log_LPS ~ LDH_activity + diagnosis_group, data = df_model_pathgroup),
  mod_ldh_logXY_path = lm(log_LPS ~ log_LDH + diagnosis_group, data = df_model_pathgroup),
  
  null_cow = lmer(sLPS_IgG ~ (1 | Cow_ID), data = df_model),
  null_farm = lmer(sLPS_IgG ~ (1 | Farm), data = df_model)
)

NAGase model estimates

Code 
nag_compare <- purrr::imap_dfr(
  all_models[grep("nag", names(all_models), ignore.case = TRUE)],
  ~ broom.mixed::tidy(.x, conf.int = TRUE) %>%
      filter(term %in% c("NAGase_activity", "log_NAG")) %>%
      mutate(Model = .y) 
%>%
  select(
    Model,
    Estimate = estimate,
    LowerCI = conf.low,
    UpperCI = conf.high,
    p.value
  ))
Code 
nag_compare_table <- nag_compare %>%
  mutate(
    Estimate = sprintf("%.7f", Estimate),
    LowerCI  = sprintf("%.7f", LowerCI),
    UpperCI  = sprintf("%.7f", UpperCI),
    `95% CI` = paste0(LowerCI, " to ", UpperCI),
    `P-value` = format.pval(p.value, digits = 3, eps = 0.001)
  ) %>%
  select(
    Model,
    Estimate,
    `95% CI`,
    `P-value`
  )

kable(
  nag_compare_table,
  caption = "Table 7. Comparison of NAGase model estimates and 95% confidence intervals across linear, mixed-effects, and generalized linear modelling approaches.",
  booktabs = TRUE,
  align = c("l", "c", "c", "c")
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center"
  )
Table 7. Comparison of NAGase model estimates and 95% confidence intervals across linear, mixed-effects, and generalized linear modelling approaches.
Model Estimate 95% CI P-value
lm_nag 0.0000218 0.0000191 to 0.0000245 <0.001
lm_nag_path 0.0000195 0.0000166 to 0.0000223 <0.001
lmer_nag_cow 0.0000218 0.0000191 to 0.0000245 <0.001
lmer_nag_farm 0.0000210 0.0000183 to 0.0000237 <0.001
lmer_nag_both 0.0000205 0.0000178 to 0.0000232 <0.001
lmer_nag_cowpath 0.0000195 0.0000166 to 0.0000223 <0.001
lmer_nag_farmpath 0.0000190 0.0000161 to 0.0000218 <0.001
lmer_nag_bothpath 0.0000186 0.0000158 to 0.0000214 <0.001
glm_nag 0.0001690 0.0001368 to 0.0002040 <0.001
glm_nag_path 0.0001426 0.0001177 to 0.0001691 <0.001
glm_nag_identity 0.0000123 0.0000098 to 0.0000153 <0.001
glm_nag_pathidentity 0.0000085 0.0000065 to 0.0000109 <0.001
mod_nag_logY 0.0001186 0.0001049 to 0.0001322 <0.001
mod_nag_logXY 0.0291226 0.0212099 to 0.0370353 <0.001
mod_nag_logY_path 0.0001016 0.0000881 to 0.0001152 <0.001
mod_nag_logXY_path 0.1718888 0.1351443 to 0.2086333 <0.001

LDH model estimates

Code 
ldh_compare <- purrr::imap_dfr(
  all_models[grep("ldh", names(all_models), ignore.case = TRUE)],
  ~ tryCatch(
    broom.mixed::tidy(.x, conf.int = TRUE) %>%
      filter(term %in% c("LDH_activity", "log_LDH")) %>%
      mutate(Model = .y),
    
    error = function(e) {
      broom.mixed::tidy(.x) %>%
        filter(term %in% c("LDH_activity", "log_LDH")) %>%
        mutate(
          conf.low = NA,
          conf.high = NA,
          Model = .y
        )
    }
  )
) %>%
  select(
    Model,
    Estimate = estimate,
    LowerCI = conf.low,
    UpperCI = conf.high,
    p.value
  )
Code 
ldh_compare_table <- ldh_compare %>%
  mutate(
    Estimate = sprintf("%.7f", Estimate),
    LowerCI  = sprintf("%.7f", LowerCI),
    UpperCI  = sprintf("%.7f", UpperCI),
    `95% CI` = paste0(LowerCI, " to ", UpperCI),
    `P-value` = format.pval(p.value, digits = 3, eps = 0.001)
  ) %>%
  select(
    Model,
    Estimate,
    `95% CI`,
    `P-value`
  )

kable(
  ldh_compare_table,
  caption = "Table 8. Comparison of LDH model estimates and 95% confidence intervals across linear, mixed-effects, and generalized linear modelling approaches.",
  booktabs = TRUE,
  align = c("l", "c", "c", "c")
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center"
  )
Table 8. Comparison of LDH model estimates and 95% confidence intervals across linear, mixed-effects, and generalized linear modelling approaches.
Model Estimate 95% CI P-value
lm_ldh 0.1437476 0.1315262 to 0.1559690 <0.001
lm_ldh_path 0.1479150 0.1340708 to 0.1617592 <0.001
lmer_ldh_cow 0.1426066 0.1304114 to 0.1548018 <0.001
lmer_ldh_farm 0.1411458 0.1287845 to 0.1535070 <0.001
lmer_ldh_both 0.1402691 0.1279304 to 0.1526077 <0.001
lmer_ldh_cowpath 0.1479154 0.1340594 to 0.1617713 <0.001
lmer_ldh_farmpath 0.1458462 0.1319836 to 0.1597089 <0.001
lmer_ldh_bothpath 0.1450147 0.1311380 to 0.1588913 <0.001
glm_ldh 0.9747019 0.8112846 to 1.1524772 <0.001
glm_ldh_path 0.8721448 0.7336041 to 1.0203876 <0.001
glm_ldh_identity 0.0734226 0.0619486 to 0.0869013 <0.001
glm_ldh_pathidentity 0.0596506 NA to NA <0.001
mod_ldh_logY 0.7392484 0.6752899 to 0.8032069 <0.001
mod_ldh_logXY 0.2768743 0.2454682 to 0.3082803 <0.001
mod_ldh_logY_path 0.6840188 0.6137678 to 0.7542698 <0.001
mod_ldh_logXY_path 1.5370789 1.3643791 to 1.7097788 <0.001

For all tested models (raw-scale, partially log-transformed, and log-log linear), the relationship between the inflammatory biomarker and sLPS/IgG is positive and significant. Confidence intervals do not include zero for any model. The gamma GLM using an identity link and including pathogen groups for LDH failed to calculate confidence intervals, likely because the link allowed the algorithm to test impossible negative values. Log transformation and Gamma GLMs did not alter the final conclusions and are not preferred because of their added complexity (model fit compared below). Retaining pathogen group in the LM NAGase model results in a change in estimate of ~12%, and so this model is preferred. For consistency and given pathogen group was identified as a possible confounder, the same model structure will be adopted for LDH, despite a change in estimate of <10%.

Model assumtions

Raw scale vs log transformations in LMs (pathogen adjusted)

Raw-scale NAGase

Code 
check_model(all_models$lm_nag_path)

LogY NAGase

Code 
check_model(all_models$mod_nag_logY_path)

Log-Log NAGase

Code 
check_model(all_models$mod_nag_logXY_path)

Raw-scale LDH

Code 
check_model(all_models$lm_ldh)

LogY LDH

Code 
check_model(all_models$mod_ldh_logY)

Log-log LDH

Code 
check_model(all_models$mod_ldh_logXY)

Raw scale LM vs gamma GLMs (pathogen adjusted)

Raw-scale NAGase

Code 
check_model(all_models$lm_nag_path)

Gamma GLM NAGase

Code 
check_model(all_models$glm_nag_path)

Raw-scale LDH

Code 
check_model(all_models$lm_ldh_path)

Gamma GLM NAGase

Code 
check_model(all_models$glm_ldh_path)

Differences in model fit are modest.

Clustering

Code 
performance::icc(all_models$null_cow)
# Intraclass Correlation Coefficient

    Adjusted ICC: 0.697
  Unadjusted ICC: 0.697
Code 
performance::icc(all_models$null_farm)
# Intraclass Correlation Coefficient

    Adjusted ICC: 0.082
  Unadjusted ICC: 0.082

Weak farm clustering and substantial cow clustering appears to exist prior to the inclusion of explanatory variables. Check when variables are added to the model:

Code 
vc <- VarCorr(all_models$lmer_nag_cowpath)

var_random <- as.numeric(vc$Cow_ID)

var_resid <- sigma(all_models$lmer_nag_cowpath)^2

icc <- var_random / (var_random + var_resid)

icc
[1] 0
Code 
vc <- VarCorr(all_models$lmer_ldh_cowpath)

var_random <- as.numeric(vc$Cow_ID)

var_resid <- sigma(all_models$lmer_nag_cowpath)^2

icc <- var_random / (var_random + var_resid)

icc
[1] 0
Code 
vc <- VarCorr(all_models$lmer_nag_farmpath)

var_random <- as.numeric(vc$Farm)

var_resid <- sigma(all_models$lmer_nag_farmpath)^2

icc <- var_random / (var_random + var_resid)

icc
[1] 0.0267942
Code 
vc <- VarCorr(all_models$lmer_ldh_farmpath)

var_random <- as.numeric(vc$Farm)

var_resid <- sigma(all_models$lmer_ldh_farmpath)^2

icc <- var_random / (var_random + var_resid)

icc
[1] 0.01990326

After adjustment for biomarker activity (NAGase or LDH) and pathogen group, the estimated cow-level ICC fell to 0 for both NAGase and LDH models, while farm-level ICCs were low (NAGase: 0.027; LDH: 0.020). These findings indicate that the apparent clustering observed in the null models was largely explained by measured covariates rather than dependence among observations from the same cow or farm. Therefore, inclusion of random intercepts for cow or farm provided little additional explanatory value, supporting the use of simpler fixed-effects models for primary inference.

NAGase regression lines

Code 
nag_models <- c(
  "lm_nag",
  "lm_nag_path",
  "lmer_nag_cow",
  "lmer_nag_farm",
  "lmer_nag_both",
  "lmer_nag_cowpath",
  "lmer_nag_farmpath",
  "lmer_nag_bothpath",
  "glm_nag",
  "glm_nag_path",
  "glm_nag_identity",
  "glm_nag_pathidentity",
  "mod_nag_logY",
  "mod_nag_logY_path"
)

pred_df_nag <- data.frame(
  NAGase_activity = seq(
    min(df_model$NAGase_activity),
    max(df_model$NAGase_activity),
    length.out = 200
  ),
  diagnosis_group = "E. coli"
)


pred_df_nag <- pred_df_nag %>%
  mutate(
    log_NAG = log1p(NAGase_activity))


plot_df_nag <- purrr::map_dfr(
  nag_models,
  function(mod_name){

    mod <- all_models[[mod_name]]

    pred <- predict(
      mod,
      newdata = pred_df_nag,
      re.form = NA,
      type = "response"
    )

    # Back-transform log(Y) models
    if(mod_name %in% c(
      "mod_nag_logY",
      "mod_nag_logY_path",
      "mod_nag_logXY_path"
    )){
      pred <- exp(pred)
    }

    pred_df_nag %>%
      mutate(
        Model = mod_name,
        Predicted_sLPS = pred
      )
  }
)

ggplot() +
  geom_point(
    data = df_model,
    aes(
      NAGase_activity,
      sLPS_IgG
    ),
    alpha = 0.15
  ) +
  geom_line(
    data = plot_df_nag,
    aes(
      NAGase_activity,
      Predicted_sLPS,
      colour = Model
    ),
    linewidth = 1
  ) +
  labs(
    x = "NAGase activity (ΔOD/min/uL)",
    y = "sLPS/IgG (OD)"
  ) +
  coord_cartesian(ylim = c(0, 3)) +
  theme_bw()

Figure 27. Predicted relationships between NAGase activity and sLPS/IgG OD obtained from candidate linear, mixed-effects, and generalised linear models. Models included unadjusted models, pathogen-adjusted models, and models incorporating cow and/or farm as random effects. The close agreement among fitted regression lines indicates that the estimated association between NAGase activity and sLPS/IgG was largely robust to model specification.

LDH regression lines

Code 
ldh_models <- c(
  "lm_ldh",
  "lm_ldh_path",
  "lmer_ldh_cow",
  "lmer_ldh_farm",
  "lmer_ldh_both",
  "lmer_ldh_cowpath",
  "lmer_ldh_farmpath",
  "lmer_ldh_bothpath",
  "glm_ldh",
  "glm_ldh_path",
  "glm_ldh_identity",
  "glm_ldh_pathidentity",
  "mod_ldh_logY",
  "mod_ldh_logY_path"
)

pred_df_ldh <- data.frame(
  LDH_activity = seq(
    min(df_model$LDH_activity),
    max(df_model$LDH_activity),
    length.out = 200
  ),
  diagnosis_group = "E. coli"
)


pred_df_ldh <- pred_df_ldh %>%
  mutate(
    log_LDH = log1p(LDH_activity))


plot_df_ldh <- purrr::map_dfr(
  ldh_models,
  function(mod_name){

    mod <- all_models[[mod_name]]

    pred <- predict(
      mod,
      newdata = pred_df_ldh,
      re.form = NA,
      type = "response"
    )

    # Back-transform log(Y) models
    if(mod_name %in% c(
      "mod_ldh_logY",
      "mod_ldh_logY_path",
      "mod_ldh_logXY_path"
    )){
      pred <- exp(pred)
    }

    pred_df_ldh %>%
      mutate(
        Model = mod_name,
        Predicted_sLPS = pred
      )
  }
)

ggplot() +
  geom_point(
    data = df_model,
    aes(
      LDH_activity,
      sLPS_IgG
    ),
    alpha = 0.15
  ) +
  geom_line(
    data = plot_df_ldh,
    aes(
      LDH_activity,
      Predicted_sLPS,
      colour = Model
    ),
    linewidth = 1
  ) +
   labs(
    x = "LDH activity (umoles/min/L)",
    y = "sLPS/IgG (OD)"
  ) +
  coord_cartesian(ylim = c(0, 3)) +
  theme_bw()

Figure 28. Predicted relationships between LDH activity and sLPS/IgG OD obtained from candidate linear, mixed-effects, and generalised linear models. Models included unadjusted models, pathogen-adjusted models, and models incorporating cow and/or farm as random effects. The close agreement among fitted regression lines indicates that the estimated association between LDH activity and sLPS/IgG was largely robust to model specification.

For both biomarkers, LM and LMMs are almost identical, and cow ID and farm contribute essentially no variation. Gamma and log transformed models show modest improvements in diagnostic performance but the underlying relationship is not altered in terms of direction, statistical significance, or biological interpretation. Over the majority of the observed range, the Gamma, log-transformed LMs and raw-scale LM fits are fairly similar (Figure 27; Figure 28). However, Gamma and log-transformed models tend to be influenced more strongly by the skewed distribution, resulting in a sharp uptick at high values. As inference is unchanged and raw-scale models are more readily interpretable, untransformed LM models will be retained for primary reporting.

2.10 Chosen models

Code 
summary(all_models$lm_nag_path)

Call:
lm(formula = sLPS_IgG ~ NAGase_activity + diagnosis_group, data = df_model_pathgroup)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.57397 -0.06435 -0.00006  0.03387  1.46698 

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)    
(Intercept)                    8.425e-02  2.568e-02   3.281 0.001094 ** 
NAGase_activity                1.946e-05  1.442e-06  13.492  < 2e-16 ***
diagnosis_groupContamination  -1.008e-01  3.034e-02  -3.321 0.000951 ***
diagnosis_groupKlebsiella     -5.176e-02  3.090e-02  -1.675 0.094438 .  
diagnosis_groupNo growth      -1.142e-01  2.704e-02  -4.222 2.79e-05 ***
diagnosis_groupPasteurella     1.525e-01  8.310e-02   1.835 0.067057 .  
diagnosis_groupSerratia       -1.233e-01  8.105e-02  -1.522 0.128585    
diagnosis_groupStaphylococcus -1.128e-01  3.946e-02  -2.858 0.004412 ** 
diagnosis_groupStreptococcus  -9.888e-02  2.829e-02  -3.495 0.000508 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1729 on 613 degrees of freedom
Multiple R-squared:  0.3241,    Adjusted R-squared:  0.3153 
F-statistic: 36.74 on 8 and 613 DF,  p-value: < 2.2e-16

After controlling for pathogen group, NAGase activity remained positively associated with sLPS/IgG OD (β = 1.95 × 10-5, p < 0.001). This equates to an increase in sLPS/IgG of approximately 0.02 units for every 1000-unit increase in NAGase activity. This model explains approximately 32% of variation in sLPS/IgG OD.

Code 
summary(all_models$lm_ldh_path)

Call:
lm(formula = sLPS_IgG ~ LDH_activity + diagnosis_group, data = df_model_pathgroup)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.53863 -0.05446  0.00877  0.02781  1.31470 

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)    
(Intercept)                    0.0027151  0.0230208   0.118  0.90615    
LDH_activity                   0.1479150  0.0070495  20.982  < 2e-16 ***
diagnosis_groupContamination  -0.0140905  0.0269338  -0.523  0.60106    
diagnosis_groupKlebsiella     -0.0009254  0.0270199  -0.034  0.97269    
diagnosis_groupNo growth      -0.0145285  0.0244283  -0.595  0.55224    
diagnosis_groupPasteurella     0.1091353  0.0717371   1.521  0.12869    
diagnosis_groupSerratia       -0.0563197  0.0705704  -0.798  0.42514    
diagnosis_groupStaphylococcus -0.0287657  0.0347404  -0.828  0.40798    
diagnosis_groupStreptococcus  -0.0889902  0.0245910  -3.619  0.00032 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1502 on 613 degrees of freedom
Multiple R-squared:  0.4898,    Adjusted R-squared:  0.4831 
F-statistic: 73.56 on 8 and 613 DF,  p-value: < 2.2e-16

After controlling for pathogen group, LDH activity remained positively associated with sLPS/IgG OD (β = 0.148, p < 0.001). This model explains approximately 48% of the variation in sLPS/IgG. LDH explained substantially more variation in sLPS/IgG OD than NAGase.

2.11 Sensitivity analysis (enrolment samples only)

Code 
lm_enrol <- koru_df %>%
  filter(Sample_Type == "Enrolment")
Code 
lm_enrol_gram <- lm_enrol %>%
  filter(!is.na(Gram_status)) %>%
  mutate(
    GN_status = ifelse(Gram_status == "Gram negative", 1, 0),
    GN_label  = ifelse(GN_status == 1, "Gram-negative", "non-Gram-negative")
  )
Code 
df_lm_pathgroup_enrol <- lm_enrol_gram %>%
  mutate(
    diagnosis_group = case_when(
      Diagnosis %in% c("Staphylococcus aureus", "Staphylococcus sp. (non aureus)") ~ "Staphylococcus",
      Diagnosis %in% c("Streptococcus dysgalactiae", "Streptococcus uberis") ~ "Streptococcus",
      Diagnosis %in% c("Serratia liquefaciens", "Serratia marcescens") ~ "Serratia",
      Diagnosis == "Klebsiella pneumoniae" ~ "Klebsiella",
      Diagnosis == "E. coli" ~ "E. coli",
      Diagnosis == "Pasteurella multocida" ~ "Pasteurella",
      Diagnosis == "No growth" ~ "No growth",
      Diagnosis == "Contamination" ~ "Contamination",
      TRUE ~ "Other"
    )
  )

df_lm_pathgroup_enrol$diagnosis_group <- factor(df_lm_pathgroup_enrol$diagnosis_group)

df_lm_pathgroup_enrol$diagnosis_group <- relevel(df_lm_pathgroup_enrol$diagnosis_group, ref = "E. coli")
Code 
mod_nag_lm_path_enrol <- lm(sLPS_IgG ~ NAGase_activity + diagnosis_group, data = df_lm_pathgroup_enrol)

There is no difference in the conclusions drawn from these models. Increasing NAGase activity is associated with increasing sLPS/IgG OD. There is a slight increase in NAGase estimate when removing follow-up samples but confidence intervals overlap and the direction and magnitude of effect remained similar.

Code 
nag_enrol_compare <- bind_rows(

  tidy(all_models$lm_nag_path, conf.int = TRUE) %>%
    mutate(Model = "LM + pathogen"),
    
  tidy(mod_nag_lm_path_enrol, conf.int = TRUE) %>%
    mutate(Model = "LM + pathogen (enrolment only)")
  
) %>%
  filter(term == "NAGase_activity") %>%
  select(
    Model,
    Estimate = estimate,
    LowerCI = conf.low,
    UpperCI = conf.high,
    p.value
  )

nag_compare_table <- nag_enrol_compare %>%
  mutate(
    p.value = ifelse(p.value < 0.001, "<0.001", as.character(signif(p.value, 3)))
  )

kable(
  nag_compare_table,
  caption = "Table 9. Comparison of NAGase model estimates and 95% confidence intervals for the full dataset and enrolment samples only dataset.",
  booktabs = TRUE,
  align = c("l", "c", "c", "c", "c")
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center"
  )
Table 9. Comparison of NAGase model estimates and 95% confidence intervals for the full dataset and enrolment samples only dataset.
Model Estimate LowerCI UpperCI p.value
LM + pathogen 1.95e-05 1.66e-05 2.23e-05 <0.001
LM + pathogen (enrolment only) 2.38e-05 2.02e-05 2.73e-05 <0.001
Code 
mod_ldh_lm_path_enrol <- lm(sLPS_IgG ~ LDH_activity + diagnosis_group, data = df_lm_pathgroup_enrol)
Code 
ldh_enrol_compare <- bind_rows(

  tidy(all_models$lm_ldh_path, conf.int = TRUE) %>%
    mutate(Model = "LM + pathogen"),
    
  tidy(mod_ldh_lm_path_enrol, conf.int = TRUE) %>%
    mutate(Model = "LM + pathogen (enrolment only)")
  
) %>%
  filter(term == "LDH_activity") %>%
  select(
    Model,
    Estimate = estimate,
    LowerCI = conf.low,
    UpperCI = conf.high,
    p.value
  )

ldh_compare_table <- ldh_enrol_compare %>%
  mutate(
    p.value = ifelse(p.value < 0.001, "<0.001", as.character(signif(p.value, 3)))
  )

kable(
  ldh_compare_table,
  caption = "Table 10. Comparison of LDH model estimates and 95% confidence intervals for the full dataset and enrolment samples only dataset.",
  booktabs = TRUE,
  align = c("l", "c", "c", "c", "c")
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center"
  )
Table 10. Comparison of LDH model estimates and 95% confidence intervals for the full dataset and enrolment samples only dataset.
Model Estimate LowerCI UpperCI p.value
LM + pathogen 0.1479150 0.1340708 0.1617592 <0.001
LM + pathogen (enrolment only) 0.1527455 0.1364615 0.1690296 <0.001

There is no difference in the conclusions drawn from these models. Increasing LDH activity is associated with increasing sLPS/IgG OD. There is a slight increase in LDH estimate when removing follow-up samples but confidence intervals overlap substantially and direction and magnitude of effect remained similar.

2.12 Objective 2 conclusions

Inflammatory marker activities and sLPS/IgG OD differed across some pathogen groups, indicating that the severity and nature of inflammation may vary slightly depending on the causative organism. However, descriptive and non-parametric analyses demonstrated variability in biomarker activity between pathogens and overlap in inflammatory responses across pathogen categories. Overall, these findings suggest that while certain pathogens may be associated with differing levels of inflammation, the relationship is not strongly discriminatory, and substantial variability exists within groups. This supports the interpretation that pathogen identity alone does not fully explain variation in inflammatory marker activities or sLPS/IgG OD, and that other host or contextual factors (e.g. immune competence, age, extent of tissue damage, pathogen load) likely contribute to observed responses. The substantial overlap among pathogen groups for all three biomarkers supports the results of Objective 1, where pathogen identification was not achievable solely on the basis of inflammatory biomarker concentrations.

Modelling results indicated that both LDH and NAGase activity were positively associated with sLPS/IgG OD, suggesting that increasing inflammatory activity corresponds to increasing ELISA signal. After adjustment for pathogen group, both NAGase activity and LDH activity remained positively associated with sLPS/IgG OD (p < 0.001 for both models), demonstrating that sLPS/IgG reflects inflammation independently of the causative organism. The estimated NAGase coefficient was modestly reduced (~11%) following pathogen adjustment, whereas the LDH coefficient remained largely unchanged, indicating minimal confounding. The LDH model explained a greater proportion of the variation in sLPS/IgG OD (adjusted R² = 0.483) than the NAGase model (adjusted R² = 0.315), suggesting a stronger relationship between LDH activity and ELISA response.

Appendix A: Alternative models considered

Code 
final_models <- c(
  "lm_nag_path",
  "lm_ldh_path"
)

rejected_models <- all_models[
  !(names(all_models) %in% final_models)
]
Code 
summarise_model <- function(model, model_name){

  cat("\n\n# ", model_name, "\n\n")

  fit <- tryCatch(
    broom::glance(model),
    error = function(e) NULL
  )

  if(!is.null(fit)){

    if("AIC" %in% names(fit))
      cat("**AIC:**", round(fit$AIC, 2), "\n\n")

    if("r.squared" %in% names(fit))
      cat("**R²:**", round(fit$r.squared, 3), "\n\n")

    if("adj.r.squared" %in% names(fit))
      cat("**Adjusted R²:**", round(fit$adj.r.squared, 3), "\n\n")
  }

  model_table <- tryCatch(
    broom.mixed::tidy(
      model,
      conf.int = FALSE
    ),
    error = function(e) broom::tidy(model)
  )

if("p.value" %in% names(model_table)) {

  model_table$p.value <- format.pval(
    model_table$p.value,
    digits = 3,
    eps = 1e-50
  )

}

print(
  kable(
    model_table,
    digits = 4
  )
)

  cat("\n\n")
}
Code 
for(nm in names(rejected_models)){

  summarise_model(
    rejected_models[[nm]],
    nm
  )

}

lm_nag

AIC: -391.04

R²: 0.29

Adjusted R²: 0.289

term estimate std.error statistic p.value
(Intercept) -0.0155 0.0093 -1.6619 0.097
NAGase_activity 0.0000 0.0000 15.9254 4.02e-48

lmer_nag_cow

AIC: -354.3

effect group term estimate std.error statistic df p.value
fixed NA (Intercept) -0.0154 0.0093 -1.6488 619 0.0997
fixed NA NAGase_activity 0.0000 0.0000 15.9140 619 4.7e-48
ran_pars Cow_ID sd__(Intercept) 0.0000 NA NA NA NA
ran_pars Residual sd__Observation 0.1763 NA NA NA NA

lmer_nag_farm

AIC: -364.33

effect group term estimate std.error statistic df p.value
fixed NA (Intercept) -0.0012 0.0128 -0.0951 40.6132 0.925
fixed NA NAGase_activity 0.0000 0.0000 15.1404 592.2162 5.18e-44
ran_pars Farm sd__(Intercept) 0.0301 NA NA NA NA
ran_pars Residual sd__Observation 0.1735 NA NA NA NA

lmer_nag_both

AIC: -360.55

effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 0.0027 0.013 0.2066 38.6455 0.837
fixed NA NAGase_activity 0.0000 0.000 14.9059 488.9279 1.08e-41
ran_pars Cow_ID sd__(Intercept) 0.1081 NA NA NA NA
ran_pars Farm sd__(Intercept) 0.0311 NA NA NA NA
ran_pars Residual sd__Observation 0.1376 NA NA NA NA

lmer_nag_cowpath

AIC: -335.51

effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 0.0842 0.0257 3.2774 612 0.001107
fixed NA NAGase_activity 0.0000 0.0000 13.4835 612 1.82e-36
fixed NA diagnosis_groupContamination -0.1008 0.0304 -3.3180 612 0.000961
fixed NA diagnosis_groupKlebsiella -0.0518 0.0309 -1.6737 612 0.094708
fixed NA diagnosis_groupNo growth -0.1142 0.0271 -4.2179 612 2.84e-05
fixed NA diagnosis_groupPasteurella 0.1524 0.0832 1.8327 612 0.067331
fixed NA diagnosis_groupSerratia -0.1233 0.0811 -1.5205 612 0.128904
fixed NA diagnosis_groupStaphylococcus -0.1114 0.0399 -2.7921 612 0.005400
fixed NA diagnosis_groupStreptococcus -0.0989 0.0283 -3.4923 612 0.000513
ran_pars Cow_ID sd__(Intercept) 0.0000 NA NA NA NA
ran_pars Residual sd__Observation 0.1730 NA NA NA NA

lmer_nag_farmpath

AIC: -342.41

effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 0.0861 0.0262 3.2806 403.2332 0.00113
fixed NA NAGase_activity 0.0000 0.0000 13.1014 604.8627 1.06e-34
fixed NA diagnosis_groupContamination -0.0945 0.0305 -3.0995 600.2024 0.00203
fixed NA diagnosis_groupKlebsiella -0.0460 0.0322 -1.4318 470.5635 0.15287
fixed NA diagnosis_groupNo growth -0.1003 0.0280 -3.5853 498.7049 0.00037
fixed NA diagnosis_groupPasteurella 0.1710 0.0835 2.0488 589.8487 0.04093
fixed NA diagnosis_groupSerratia -0.1271 0.0809 -1.5722 611.2451 0.11643
fixed NA diagnosis_groupStaphylococcus -0.1169 0.0396 -2.9514 596.7209 0.00329
fixed NA diagnosis_groupStreptococcus -0.0879 0.0286 -3.0693 594.6387 0.00224
ran_pars Farm sd__(Intercept) 0.0284 NA NA NA NA
ran_pars Residual sd__Observation 0.1709 NA NA NA NA

lmer_nag_bothpath

AIC: -338.58

effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 0.0911 0.0263 3.4612 370.6281 0.000600
fixed NA NAGase_activity 0.0000 0.0000 12.8954 518.8422 3.36e-33
fixed NA diagnosis_groupContamination -0.0961 0.0307 -3.1316 572.9838 0.001827
fixed NA diagnosis_groupKlebsiella -0.0505 0.0322 -1.5678 474.9888 0.117605
fixed NA diagnosis_groupNo growth -0.1032 0.0280 -3.6851 509.6495 0.000253
fixed NA diagnosis_groupPasteurella 0.1630 0.0830 1.9638 471.7033 0.050146
fixed NA diagnosis_groupSerratia -0.1277 0.0796 -1.6038 207.5922 0.110275
fixed NA diagnosis_groupStaphylococcus -0.1168 0.0404 -2.8915 587.7185 0.003975
fixed NA diagnosis_groupStreptococcus -0.0912 0.0287 -3.1816 596.5758 0.001541
ran_pars Cow_ID sd__(Intercept) 0.0990 NA NA NA NA
ran_pars Farm sd__(Intercept) 0.0295 NA NA NA NA
ran_pars Residual sd__Observation 0.1407 NA NA NA NA

glm_nag

AIC: -2449.35

term estimate std.error statistic p.value
(Intercept) -3.7295 0.1008 -37.0032 < 1e-50
NAGase_activity 0.0002 0.0000 11.4262 1.43e-27

glm_nag_path

AIC: -2610.28

term estimate std.error statistic p.value
(Intercept) -2.7376 0.2152 -12.7200 4.63e-33
NAGase_activity 0.0001 0.0000 11.8023 4.12e-29
diagnosis_groupContamination -1.0797 0.2543 -4.2460 2.51e-05
diagnosis_groupKlebsiella -0.3276 0.2590 -1.2652 0.20628
diagnosis_groupNo growth -1.4298 0.2266 -6.3084 5.40e-10
diagnosis_groupPasteurella -0.3869 0.6965 -0.5556 0.57871
diagnosis_groupSerratia -1.3105 0.6793 -1.9292 0.05417
diagnosis_groupStaphylococcus -1.3255 0.3307 -4.0082 6.87e-05
diagnosis_groupStreptococcus -0.8903 0.2371 -3.7546 0.00019

glm_nag_identity

AIC: -2442.02

term estimate std.error statistic p.value
(Intercept) 0.0161 0.0026 6.1597 1.31e-09
NAGase_activity 0.0000 0.0000 7.9914 6.53e-15

glm_nag_pathidentity

AIC: -2526.33

term estimate std.error statistic p.value
(Intercept) 0.1049 0.0286 3.6724 0.000261
NAGase_activity 0.0000 0.0000 7.6639 7.09e-14
diagnosis_groupContamination -0.0888 0.0288 -3.0864 0.002118
diagnosis_groupKlebsiella -0.0541 0.0311 -1.7362 0.083037
diagnosis_groupNo growth -0.0931 0.0285 -3.2637 0.001161
diagnosis_groupPasteurella 0.2993 0.3645 0.8210 0.411973
diagnosis_groupSerratia -0.1015 0.0337 -3.0124 0.002698
diagnosis_groupStaphylococcus -0.0938 0.0291 -3.2231 0.001336
diagnosis_groupStreptococcus -0.0792 0.0289 -2.7374 0.006372

mod_nag_logY

AIC: 1632.22

R²: 0.319

Adjusted R²: 0.318

term estimate std.error statistic p.value
(Intercept) -4.0370 0.0474 -85.1461 <1e-50
NAGase_activity 0.0001 0.0000 17.0419 <1e-50

mod_nag_logXY

AIC: -228.05

R²: 0.078

Adjusted R²: 0.076

term estimate std.error statistic p.value
(Intercept) -0.1349 0.031 -4.3520 1.58e-05
log_NAG 0.0291 0.004 7.2277 1.45e-12

mod_nag_logY_path

AIC: 1540.95

R²: 0.425

Adjusted R²: 0.418

term estimate std.error statistic p.value
(Intercept) -3.1365 0.1230 -25.5075 < 1e-50
NAGase_activity 0.0001 0.0000 14.7181 3.29e-42
diagnosis_groupContamination -0.8921 0.1453 -6.1400 1.48e-09
diagnosis_groupKlebsiella -0.3399 0.1480 -2.2971 0.0220
diagnosis_groupNo growth -1.1648 0.1295 -8.9952 2.92e-18
diagnosis_groupPasteurella -0.4991 0.3979 -1.2542 0.2103
diagnosis_groupSerratia -0.7963 0.3881 -2.0518 0.0406
diagnosis_groupStaphylococcus -0.9817 0.1889 -5.1960 2.78e-07
diagnosis_groupStreptococcus -0.7860 0.1355 -5.8015 1.05e-08

mod_nag_logXY_path

AIC: 1648.94

R²: 0.316

Adjusted R²: 0.307

term estimate std.error statistic p.value
(Intercept) -3.8714 0.1959 -19.7663 < 1e-50
log_NAG 0.1719 0.0187 9.1867 6.19e-19
diagnosis_groupContamination -1.0042 0.1581 -6.3533 4.11e-10
diagnosis_groupKlebsiella -0.4128 0.1612 -2.5599 0.01071
diagnosis_groupNo growth -1.3172 0.1405 -9.3748 1.32e-19
diagnosis_groupPasteurella 0.5235 0.4240 1.2347 0.21742
diagnosis_groupSerratia -1.0985 0.4227 -2.5986 0.00959
diagnosis_groupStaphylococcus -1.2147 0.2050 -5.9254 5.20e-09
diagnosis_groupStreptococcus -0.7984 0.1478 -5.4016 9.46e-08

lm_ldh

AIC: -563.91

R²: 0.463

Adjusted R²: 0.462

term estimate std.error statistic p.value
(Intercept) -0.0218 0.0076 -2.8674 0.00428
LDH_activity 0.1437 0.0062 23.0982 < 1e-50

lmer_ldh_cow

AIC: -543.56

effect group term estimate std.error statistic df p.value
fixed NA (Intercept) -0.0199 0.0077 -2.5687 419.7140 0.0106
fixed NA LDH_activity 0.1426 0.0062 22.9664 589.0199 <1e-50
ran_pars Cow_ID sd__(Intercept) 0.0786 NA NA NA NA
ran_pars Residual sd__Observation 0.1323 NA NA NA NA

lmer_ldh_farm

AIC: -550.07

effect group term estimate std.error statistic df p.value
fixed NA (Intercept) -0.0123 0.0098 -1.2512 39.3037 0.218
fixed NA LDH_activity 0.1411 0.0063 22.4257 589.9302 <1e-50
ran_pars Farm sd__(Intercept) 0.0194 NA NA NA NA
ran_pars Residual sd__Observation 0.1519 NA NA NA NA

lmer_ldh_both

AIC: -546.19

effect group term estimate std.error statistic df p.value
fixed NA (Intercept) -0.0108 0.0099 -1.0884 38.3164 0.283
fixed NA LDH_activity 0.1403 0.0063 22.3311 543.3744 <1e-50
ran_pars Cow_ID sd__(Intercept) 0.0734 NA NA NA NA
ran_pars Farm sd__(Intercept) 0.0195 NA NA NA NA
ran_pars Residual sd__Observation 0.1337 NA NA NA NA

lmer_ldh_cowpath

AIC: -524.84

effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 0.0027 0.0230 0.1178 612 0.906246
fixed NA LDH_activity 0.1479 0.0071 20.9646 612 < 1e-50
fixed NA diagnosis_groupContamination -0.0141 0.0270 -0.5227 612 0.601362
fixed NA diagnosis_groupKlebsiella -0.0009 0.0270 -0.0342 612 0.972717
fixed NA diagnosis_groupNo growth -0.0145 0.0244 -0.5942 612 0.552572
fixed NA diagnosis_groupPasteurella 0.1091 0.0718 1.5201 612 0.129010
fixed NA diagnosis_groupSerratia -0.0563 0.0706 -0.7974 612 0.425523
fixed NA diagnosis_groupStaphylococcus -0.0288 0.0351 -0.8202 612 0.412421
fixed NA diagnosis_groupStreptococcus -0.0890 0.0246 -3.6159 612 0.000324
ran_pars Cow_ID sd__(Intercept) 0.0000 NA NA NA NA
ran_pars Residual sd__Observation 0.1503 NA NA NA NA

lmer_ldh_farmpath

AIC: -531.14

effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 0.0057 0.0233 0.2422 480.6195 0.80870
fixed NA LDH_activity 0.1458 0.0071 20.6614 610.3359 < 1e-50
fixed NA diagnosis_groupContamination -0.0084 0.0270 -0.3108 602.2124 0.75602
fixed NA diagnosis_groupKlebsiella 0.0082 0.0280 0.2935 455.3441 0.76928
fixed NA diagnosis_groupNo growth -0.0015 0.0251 -0.0585 502.8359 0.95337
fixed NA diagnosis_groupPasteurella 0.1209 0.0720 1.6791 591.3249 0.09367
fixed NA diagnosis_groupSerratia -0.0603 0.0704 -0.8567 611.6803 0.39195
fixed NA diagnosis_groupStaphylococcus -0.0339 0.0348 -0.9724 600.0989 0.33124
fixed NA diagnosis_groupStreptococcus -0.0796 0.0249 -3.2029 592.1319 0.00143
ran_pars Farm sd__(Intercept) 0.0212 NA NA NA NA
ran_pars Residual sd__Observation 0.1488 NA NA NA NA

lmer_ldh_bothpath

AIC: -527.22

effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 0.0085 0.0234 0.3615 460.8171 0.71788
fixed NA LDH_activity 0.1450 0.0071 20.5230 607.2846 < 1e-50
fixed NA diagnosis_groupContamination -0.0098 0.0271 -0.3608 581.4613 0.71839
fixed NA diagnosis_groupKlebsiella 0.0047 0.0280 0.1686 458.2161 0.86617
fixed NA diagnosis_groupNo growth -0.0021 0.0251 -0.0831 508.6241 0.93384
fixed NA diagnosis_groupPasteurella 0.1125 0.0717 1.5689 497.3873 0.11730
fixed NA diagnosis_groupSerratia -0.0654 0.0696 -0.9401 252.3770 0.34807
fixed NA diagnosis_groupStaphylococcus -0.0339 0.0355 -0.9549 580.1363 0.34004
fixed NA diagnosis_groupStreptococcus -0.0815 0.0249 -3.2752 593.9839 0.00112
ran_pars Cow_ID sd__(Intercept) 0.0807 NA NA NA NA
ran_pars Farm sd__(Intercept) 0.0217 NA NA NA NA
ran_pars Residual sd__Observation 0.1260 NA NA NA NA

glm_ldh

AIC: -2609.2

term estimate std.error statistic p.value
(Intercept) -3.7843 0.0837 -45.2179 < 1e-50
LDH_activity 0.9747 0.0685 14.2211 6.31e-40

glm_ldh_path

AIC: -2728.24

term estimate std.error statistic p.value
(Intercept) -3.2125 0.2062 -15.5819 2.32e-46
LDH_activity 0.8721 0.0631 13.8144 5.50e-38
diagnosis_groupContamination -0.5294 0.2412 -2.1947 0.02856
diagnosis_groupKlebsiella 0.1697 0.2420 0.7014 0.48334
diagnosis_groupNo growth -0.8677 0.2188 -3.9664 8.16e-05
diagnosis_groupPasteurella -1.1486 0.6425 -1.7878 0.07430
diagnosis_groupSerratia -0.8237 0.6320 -1.3033 0.19296
diagnosis_groupStaphylococcus -0.7085 0.3111 -2.2774 0.02311
diagnosis_groupStreptococcus -0.6944 0.2202 -3.1530 0.00169

glm_ldh_identity

AIC: -2625.79

term estimate std.error statistic p.value
(Intercept) 0.0152 0.0016 9.6662 1.11e-20
LDH_activity 0.0734 0.0068 10.7922 5.32e-25

glm_ldh_pathidentity

AIC: -2655.25

term estimate std.error statistic p.value
(Intercept) 0.0472 0.0144 3.2794 0.0011
LDH_activity 0.0597 0.0064 9.3574 1.52e-19
diagnosis_groupContamination -0.0301 0.0147 -2.0481 0.0410
diagnosis_groupKlebsiella 0.0035 0.0187 0.1869 0.8518
diagnosis_groupNo growth -0.0331 0.0144 -2.2919 0.0223
diagnosis_groupPasteurella 0.1497 0.2098 0.7135 0.4758
diagnosis_groupSerratia -0.0416 0.0212 -1.9680 0.0495
diagnosis_groupStaphylococcus -0.0339 0.0153 -2.2196 0.0268
diagnosis_groupStreptococcus -0.0352 0.0147 -2.3934 0.0170

mod_ldh_logY

AIC: 1494.97

R²: 0.454

Adjusted R²: 0.453

term estimate std.error statistic p.value
(Intercept) -4.0404 0.0398 -101.5988 <1e-50
LDH_activity 0.7392 0.0326 22.6981 <1e-50

mod_ldh_logXY

AIC: -423.03

R²: 0.326

Adjusted R²: 0.325

term estimate std.error statistic p.value
(Intercept) -0.0383 0.0098 -3.9231 9.72e-05
log_LDH 0.2769 0.0160 17.3128 < 1e-50

mod_ldh_logY_path

AIC: 1438.21

R²: 0.513

Adjusted R²: 0.506

term estimate std.error statistic p.value
(Intercept) -3.4391 0.1168 -29.4402 < 1e-50
LDH_activity 0.6840 0.0358 19.1215 < 1e-50
diagnosis_groupContamination -0.5170 0.1367 -3.7827 0.000170
diagnosis_groupKlebsiella -0.1155 0.1371 -0.8422 0.400013
diagnosis_groupNo growth -0.7449 0.1240 -6.0092 3.2e-09
diagnosis_groupPasteurella -0.5473 0.3640 -1.5035 0.133229
diagnosis_groupSerratia -0.5202 0.3581 -1.4527 0.146816
diagnosis_groupStaphylococcus -0.6276 0.1763 -3.5603 0.000399
diagnosis_groupStreptococcus -0.7485 0.1248 -5.9985 3.4e-09

mod_ldh_logXY_path

AIC: 1477.65

R²: 0.481

Adjusted R²: 0.474

term estimate std.error statistic p.value
(Intercept) -3.5814 0.1258 -28.4768 < 1e-50
log_LDH 1.5371 0.0879 17.4788 < 1e-50
diagnosis_groupContamination -0.5361 0.1413 -3.7948 0.000162
diagnosis_groupKlebsiella -0.2491 0.1409 -1.7684 0.077487
diagnosis_groupNo growth -0.7190 0.1293 -5.5600 4.03e-08
diagnosis_groupPasteurella -0.2704 0.3737 -0.7237 0.469514
diagnosis_groupSerratia -0.6537 0.3692 -1.7706 0.077128
diagnosis_groupStaphylococcus -0.6756 0.1818 -3.7153 0.000222
diagnosis_groupStreptococcus -0.8403 0.1287 -6.5307 1.38e-10

null_cow

AIC: -180.86

effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 0.0853 0.0089 9.538 533.261 5.15e-20
ran_pars Cow_ID sd__(Intercept) 0.1786 NA NA NA NA
ran_pars Residual sd__Observation 0.1176 NA NA NA NA

null_farm

AIC: -200.2

effect group term estimate std.error statistic df p.value
fixed NA (Intercept) 0.1068 0.016 6.6874 28.3607 2.76e-07
ran_pars Farm sd__(Intercept) 0.0600 NA NA NA NA
ran_pars Residual sd__Observation 0.2004 NA NA NA NA

References

Rowe, S., House, J. K., Pooley, H., Bullen, S., Humphris, M., Ingenhoff, L., Norris, J. M., & Zadoks, R. N. (2024). Evaluation of point-of-care tests for identification of pathogens to inform clinical mastitis treatment decisions in pasture- and confinement-managed dairy cows in Australia. Journal of Dairy Science, 107(10), 8271–8285. https://doi.org/10.3168/jds.2023-24612