Prediction of transfer of passive immunity in Holstein and Holstein x Angus calves

Author

Dini Hapukotuwa and Sam Rowe

Published

July 15, 2025

Authors and affiliations

Dini Hapukotuwa1, John House1, Katie Denholm2, and Sam Rowe1

1Sydney School of Veterinary Science, University of Sydney, Camden 2750, NSW, Australia

2School of Biodiversity, One Health and Veterinary Medicine, University of Glasgow, Glasgow, Scotland

Acknowledgements

We would like to thank the participating farm, Anna Waldron, and Jennie Mohler

Study objectives

To determine if there is a difference in the relationship between serum total solids (STS) and serum IgG for Holstein and Holstein X Angus calves.

  1. Describe effect of breed on STS and serum IgG

  2. Compare the correlation between STS and IgG between breeds

  3. Compare the use of STS to predict ‘excellent’ transfer of passive immunity between breeds

Import data

Show the code 
library(tidyverse); library(janitor); library(epiR); library(readxl); library(here); library(emmeans); library(table1); library(pROC); library(performance); library(plotROC);library(ggbeeswarm); library(ggpubr);library(DiagrammeR)


data <- readRDS(here("Data in/igg_trial.rds")) %>%
  select(calf_id:sts,igg,tx_group,breed:lost_to_follow_up_reason,time_to_feed,vol,
         igg_gt_10:sts_grams_l)

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

Data Frame Summary

data

Dimensions: 400 x 21
Duplicates: 0
No Variable Stats / Values Freqs (% of Valid) Graph Missing
1 calf_id [numeric]
Mean (sd) : 426256.8 (426039.1)
min ≤ med ≤ max:
108720 ≤ 108936.5 ≤ 996765
IQR (CV) : 887831.5 (1)
400 distinct values 0 (0.0%)
2 enrol_date [POSIXct, POSIXt]
min : 2024-10-24
med : 2024-11-03
max : 2024-11-10
range : 17d
18 distinct values 0 (0.0%)
3 time_born [POSIXct, POSIXt]
min : 2024-10-24 11:00:00
med : 2024-11-03 12:15:00
max : 2024-11-10 23:45:00
range : 17d 12H 45M 0S
383 distinct values 0 (0.0%)
4 time_fed [POSIXct, POSIXt]
min : 2024-10-24 11:40:00
med : 2024-11-03 13:17:30
max : 2024-11-11 00:30:00
range : 17d 12H 50M 0S
386 distinct values 0 (0.0%)
5 sts [numeric]
Mean (sd) : 6 (0.6)
min ≤ med ≤ max:
4 ≤ 6 ≤ 8
IQR (CV) : 0.6 (0.1)
32 distinct values 2 (0.5%)
6 igg [numeric]
Mean (sd) : 31.5 (7.8)
min ≤ med ≤ max:
7.6 ≤ 32 ≤ 53.9
IQR (CV) : 10.6 (0.2)
367 distinct values 8 (2.0%)
7 tx_group [factor]
1. 200g
2. 250g
3. 300g
4. 350g
45 ( 11.2% )
125 ( 31.2% )
123 ( 30.8% )
107 ( 26.8% )
0 (0.0%)
8 breed [character]
1. Angus x Holstein
2. Holstein
126 ( 31.5% )
274 ( 68.5% )
0 (0.0%)
9 sex [character]
1. F
2. M
319 ( 79.8% )
81 ( 20.2% )
0 (0.0%)
10 excluded [numeric] 1 distinct value
0 : 400 ( 100.0% )
0 (0.0%)
11 excluded_reason [character]
All NA's
400 (100.0%)
12 lost_to_follow_up [numeric]
Min : 0
Mean : 0
Max : 1
0 : 399 ( 99.8% )
1 : 1 ( 0.2% )
0 (0.0%)
13 lost_to_follow_up_reason [character] 1. Died 3/11/24
1 ( 100.0% )
399 (99.8%)
14 time_to_feed [numeric]
Mean (sd) : 59.7 (29.4)
min ≤ med ≤ max:
0 ≤ 50 ≤ 195
IQR (CV) : 35 (0.5)
37 distinct values 0 (0.0%)
15 vol [numeric]
Mean (sd) : 3.1 (0.6)
min ≤ med ≤ max:
1.8 ≤ 3.1 ≤ 3.9
IQR (CV) : 0.8 (0.2)
22 distinct values 0 (0.0%)
16 igg_gt_10 [numeric]
Min : 0
Mean : 1
Max : 1
0 : 2 ( 0.5% )
1 : 390 ( 99.5% )
8 (2.0%)
17 igg_gt_25 [factor]
1. IgG >= 25
2. IgG < 25
306 ( 78.1% )
86 ( 21.9% )
8 (2.0%)
18 sts_gte_6.2 [factor]
1. STS >= 6.2
2. STS < 6.2
141 ( 35.4% )
257 ( 64.6% )
2 (0.5%)
19 igg_category [factor]
1. Excellent (>= 25g/L)
2. Good (18 to <25g/L)
3. Fair (10 to <18g/L)
4. Poor (< 10g/L)
306 ( 78.1% )
69 ( 17.6% )
15 ( 3.8% )
2 ( 0.5% )
8 (2.0%)
20 sts_category [factor]
1. Excellent (>= 6.2g/dL)
2. Good (5.8 to 6.1g/dL)
3. Fair (5.1 to <5.7g/dL)
4. Poor (< 5.1g/dL)
141 ( 37.9% )
133 ( 35.8% )
98 ( 26.3% )
0 ( 0.0% )
28 (7.0%)
21 sts_grams_l [numeric]
Mean (sd) : 59.8 (6.2)
min ≤ med ≤ max:
40 ≤ 60 ≤ 80
IQR (CV) : 6 (0.1)
32 distinct values 2 (0.5%)

Generated by summarytools 1.1.4 (R version 4.5.1)
2025-07-15

Causal diagram (directed acyclic graph)

This diagram demonstrates that adjusting for IgG dose will correct for the unequal IgG doses administered to Holstein and Angus X calves, due to having 2 dose randomization periods in the clinical trial that these data are collected from.

Adjusting for sex will remove the indirect effect of breed on serum IgG mediated via the higher prevalence of male calves in the Angus X group.

The hashed lines indicate the causal effects of interest (breed impacting IgG, and breed impacting STS).

Show the code 
DiagrammeR::grViz("
digraph {
  graph [ranksep = 0.2]

  node [shape = plaintext]
    breed  [label = 'Breed']
    random    [label = 'Randomization period']
    dose   [label = 'Colostrum IgG Dose', shape = box]
    vol   [label = 'Colostrum volume']
    time  [label = 'Time to feeding']
    sex [label = 'Sex', shape = box]
    igg  [label = 'Serum IgG']
    sts [label = 'Serum Total Solids']
    non_igg [label = 'Non-IgG proteins']

  edge [minlen = 2]
    breed->sex
    random->breed
    random->dose
    dose->vol
    dose->igg
    vol->igg
    time->igg
    sex->igg
    igg->sts
    non_igg->sts

  edge [style = dashed, minlen = 5]
    breed->non_igg
    breed->sts
    breed->igg

  { rank = same }
}
")

Descriptive stats

Show the code 
data <- data %>% filter(
  !is.na(igg) & !is.na(sts_grams_l)) %>%
  mutate(
    dose = case_when(
      tx_group == "200g" ~ 200,
      tx_group == "250g" ~ 250,
      tx_group == "300g" ~ 300,
      tx_group == "350g" ~ 350
    )
  )

table1(~ tx_group + dose + time_to_feed + vol + sex | factor(breed), data=data)
Angus x Holstein
(N=119)		 Holstein
(N=272)		 Overall
(N=391)
tx_group
200g 21 (17.6%) 24 (8.8%) 45 (11.5%)
250g 35 (29.4%) 88 (32.4%) 123 (31.5%)
300g 39 (32.8%) 79 (29.0%) 118 (30.2%)
350g 24 (20.2%) 81 (29.8%) 105 (26.9%)
dose
Mean (SD) 278 (50.3) 290 (48.4) 286 (49.2)
Median [Min, Max] 300 [200, 350] 300 [200, 350] 300 [200, 350]
time_to_feed
Mean (SD) 57.3 (28.9) 60.7 (29.9) 59.7 (29.6)
Median [Min, Max] 50.0 [0, 190] 50.0 [0, 195] 50.0 [0, 195]
vol
Mean (SD) 3.02 (0.591) 3.10 (0.538) 3.08 (0.555)
Median [Min, Max] 3.10 [1.80, 3.90] 3.10 [1.80, 3.90] 3.10 [1.80, 3.90]
sex
F 60 (50.4%) 255 (93.8%) 315 (80.6%)
M 59 (49.6%) 17 (6.3%) 76 (19.4%)

Is IgG concentration in serum affected by breed?

Show the code 
table1(~ igg + sts_grams_l | factor(breed), data=data)
Angus x Holstein
(N=119)		 Holstein
(N=272)		 Overall
(N=391)
igg
Mean (SD) 30.2 (7.30) 32.1 (8.02) 31.5 (7.84)
Median [Min, Max] 30.8 [7.60, 46.4] 32.8 [8.79, 53.9] 32.1 [7.60, 53.9]
sts_grams_l
Mean (SD) 56.9 (5.47) 61.0 (6.08) 59.8 (6.19)
Median [Min, Max] 58.0 [40.0, 72.0] 60.0 [45.0, 80.0] 60.0 [40.0, 80.0]
Show the code 
figure <- ggplot(data,aes(x = breed, y = igg)) +
  # Add boxplot with fixed fill color
  geom_boxplot(fill = "white",  # Set all boxplots to light blue
               outlier.shape = 16, outlier.size = 2, alpha = 1) +
  
  # Labels
  labs(
    title = "",
    x = "Breed",
    y = "Serum IgG (g/L)",
    fill = "Treatment"
  ) +
  
  # Custom theme
  theme_minimal() +
  coord_cartesian(ylim = c(0, 60)) +
  theme(
    legend.position = "none",
    panel.border = element_rect(
      color = "black", fill = NA, linewidth = 1
    ),
    text = element_text(family = "Times New Roman")
  )

figure

Model based on DAG. Adjusts for oral IgG dose and sex.

Show the code 
lm_adjust_dose_sex <- lm(igg ~ breed + tx_group + sex, data = data)
summary(lm_adjust_dose_sex)

Call:
lm(formula = igg ~ breed + tx_group + sex, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-23.3755  -4.6829   0.1765   4.8955  19.0287 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    25.0439     1.2301  20.360  < 2e-16 ***
breedHolstein   0.3418     0.9123   0.375 0.708128    
tx_group250g    4.6900     1.2395   3.784 0.000179 ***
tx_group300g    6.7804     1.2406   5.465 8.32e-08 ***
tx_group350g   11.1605     1.2738   8.762  < 2e-16 ***
sexM           -1.3979     1.0502  -1.331 0.183962    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.055 on 385 degrees of freedom
Multiple R-squared:  0.2017,    Adjusted R-squared:  0.1913 
F-statistic: 19.45 on 5 and 385 DF,  p-value: < 2.2e-16
Show the code 
emmeans(lm_adjust_dose_sex,~breed)
 breed            emmean    SE  df lower.CL upper.CL
 Angus x Holstein   30.0 0.655 385     28.7     31.3
 Holstein           30.3 0.661 385     29.0     31.6

Results are averaged over the levels of: tx_group, sex 
Confidence level used: 0.95

Univariable model

Show the code 
lm <- lm(igg ~ breed, data = data)
summary(lm)

Call:
lm(formula = igg ~ breed, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-23.3179  -5.7253   0.6473   4.9716  21.8270 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    30.2033     0.7155  42.211   <2e-16 ***
breedHolstein   1.9052     0.8579   2.221   0.0269 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.806 on 389 degrees of freedom
Multiple R-squared:  0.01252,   Adjusted R-squared:  0.009982 
F-statistic: 4.932 on 1 and 389 DF,  p-value: 0.02694

Maximal model

Show the code 
lm <- lm(igg ~ breed + tx_group + sex + time_to_feed + vol, data = data)
summary(lm)

Call:
lm(formula = igg ~ breed + tx_group + sex + time_to_feed + vol, 
    data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-21.7424  -4.7066   0.1031   4.7301  18.6934 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   29.22072    2.68176  10.896  < 2e-16 ***
breedHolstein  0.42735    0.90861   0.470   0.6384    
tx_group250g   5.50600    1.36260   4.041 6.44e-05 ***
tx_group300g   8.12175    1.67123   4.860 1.72e-06 ***
tx_group350g  12.94647    1.88861   6.855 2.86e-11 ***
sexM          -1.30078    1.04553  -1.244   0.2142    
time_to_feed  -0.02580    0.01208  -2.137   0.0332 *  
vol           -1.25352    1.04297  -1.202   0.2302    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.017 on 383 degrees of freedom
Multiple R-squared:  0.2143,    Adjusted R-squared:  0.1999 
F-statistic: 14.92 on 7 and 383 DF,  p-value: < 2.2e-16

Comparison of effect estimates with different models

Show the code 
lm_univariable <- lm(igg ~ breed, data = data)

lm_adjust_all <- lm(igg ~ breed + tx_group + sex + time_to_feed + vol, data = data)

lm_adjust_dose_sex <- lm(igg ~ breed + tx_group + sex, data = data)

lm_adjust_dose <- lm(igg ~ breed + tx_group, data = data)

rbind(
    broom::tidy(lm_adjust_dose_sex, conf.int=T) %>% 
    select(term,estimate,conf.low,conf.high) %>%
    filter(term == "breedHolstein") %>%
    mutate(model = "DAG model - Adjusted for dose, sex"),
    
    broom::tidy(lm_univariable, conf.int=T) %>% 
    select(term,estimate,conf.low,conf.high) %>%
    filter(term == "breedHolstein") %>%
    mutate(model = "No adjustments"),

  broom::tidy(lm_adjust_all, conf.int=T) %>% 
    select(term,estimate,conf.low,conf.high) %>%
    filter(term == "breedHolstein") %>%
    mutate(model = "Adjusted for dose, sex, time to feed, volume"),
  
  broom::tidy(lm_adjust_dose, conf.int=T) %>% 
    select(term,estimate,conf.low,conf.high) %>%
    filter(term == "breedHolstein") %>%
    mutate(model = "Adjusted for dose")
)
term estimate conf.low conf.high model
breedHolstein 0.3417896 -1.4519115 2.135491 DAG model - Adjusted for dose, sex
breedHolstein 1.9052491 0.2185495 3.591949 No adjustments
breedHolstein 0.4273528 -1.3591309 2.213836 Adjusted for dose, sex, time to feed, volume
breedHolstein 0.9615939 -0.5823834 2.505571 Adjusted for dose

Is total protein concentration in serum affected by breed?

Show the code 
figure <- ggplot(data,aes(x = breed, y = sts_grams_l)) +
  # Add boxplot with fixed fill color
  geom_boxplot(fill = "white",  # Set all boxplots to light blue
               outlier.shape = 16, outlier.size = 2, alpha = 1) +
  
  # Labels
  labs(
    title = "",
    x = "Breed",
    y = "Serum Total Solids (g/L)"
  ) +
  
  # Custom theme
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.border = element_rect(
      color = "black", fill = NA, linewidth = 1
    ),
    text = element_text(family = "Times New Roman")
  )

figure

Show the code 
# Reshape the data to long format
data_long <- data %>%
  select(breed, igg, sts_grams_l) %>%
  pivot_longer(cols = c(igg, sts_grams_l), names_to = "measurement", values_to = "value") %>%
  mutate(
    measurement = recode(measurement,
                         igg = "Serum IgG",
                         sts_grams_l = "Serum Total Solids")
  )

# Plot both panels
figure <- ggplot(data_long, aes(x = breed, y = value)) +
  geom_jitter() +
  # geom_boxplot(fill = "white",
  #              outlier.shape = 16, outlier.size = 2, alpha = 1) +
  labs(
    x = NULL,
    y = "Grams per L"
  ) +
  facet_wrap(~ measurement, ncol = 2, scales = "free_y") +  # One column, different y-axes
  coord_cartesian(ylim = c(0, 80)) +  # This will apply only to the IgG plot
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1),
    text = element_text(family = "Times New Roman"),
    strip.text = element_text(size = 12, face = "bold")
  )

print(figure)

Show the code 
figure <- ggplot(data_long, aes(x = breed, y = value)) +
  geom_quasirandom(alpha=.4,colour = "blue") +
  labs(
    x = NULL,
    y = "Concentration in serum (G/L)"
  ) +
  facet_wrap(~ measurement, ncol = 2) +
  coord_cartesian(ylim = c(0, 100)) +
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1),
    text = element_text(family = "Times New Roman"),
    strip.text = element_text(size = 12, face = "bold")
  )

figure

Show the code 
ggsave(figure,filename = "short communication figure 1.tiff",dpi = 300,width = 7, height = 4)
Show the code 
# Create the background band data
background_df <- data.frame(
  measurement = rep(c("Serum IgG", "Serum Total Solids"), times = c(4, 4)),
  ymin = c(0, 10, 18, 25, 0, 51, 58, 62),
  ymax = c(10, 18, 25, 100, 51, 58, 62, 100),
  fill = c(
    "Poor (<10 g/L)", "Fair (10–17.9 g/L)", "Good (18–24.9 g/L)", "Excellent (≥25 g/L)",
    "Poor (<51 g/dL)", "Fair (51–58 g/dL)", "Good (58–62 g/dL)", "Excellent (≥62 g/dL)"
  )
)

    
# Assign colors (same as your previous figure)
fill_colors <- c(
  "Poor (<10 g/L)" = "pink",
  "Fair (10–17.9 g/L)" = "#FFF2CC",
  "Good (18–24.9 g/L)" = "#C3D69B",
  "Excellent (≥25 g/L)" = "#A2D9CE",
  "Poor (<51 g/dL)" = "pink",
  "Fair (51–58 g/dL)" = "#FFF2CC",
  "Good (58–62 g/dL)" = "#C3D69B",
  "Excellent (≥62 g/dL)" = "#A2D9CE"
)

figure <- ggplot(data_long, aes(x = breed, y = value)) +
  # background bands
  geom_rect(
    data = background_df,
    aes(xmin = -Inf, xmax = Inf, ymin = ymin, ymax = ymax, fill = fill),
    inherit.aes = FALSE,
    alpha = 1
  ) +
  # horizontal lines every 5
  geom_hline(
  yintercept = seq(0, 100, by = 5),
  color = "black",
  size = 0.3,
  alpha = 0.5,
  linetype = "dashed") +
  # your points
  geom_quasirandom(alpha = 0.4, colour = "blue") +
  scale_fill_manual(values = fill_colors) +
  # set y-axis breaks/labels at every 5
  scale_y_continuous(
    breaks = seq(0, 100, by = 5),
    limits = c(0, 100),
    expand = expansion(mult = c(0, 0))
  ) +
  labs(
    x = NULL,
    y = "Concentration in serum (G/L)"
  ) +
  facet_wrap(~ measurement, ncol = 2, scales = "fixed") +
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1),
    panel.grid.major = element_blank(),   # turn off default grid
    panel.grid.minor = element_blank(),
    text = element_text(family = "Times New Roman"),
    strip.text = element_text(size = 12, face = "bold")
  )

figure

Show the code 
ggsave(figure,filename = "short communication figure 1a.tiff",dpi = 300,width = 7, height = 4)

Model based on DAG. Adjusts for dose and sex.

Show the code 
lm_adjust_dose_sex <- lm(sts_grams_l ~ breed + tx_group + sex, data = data)
lm_adjust_dose_sex %>% summary

Call:
lm(formula = sts_grams_l ~ breed + tx_group + sex, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-16.724  -3.631  -0.344   3.701  18.462 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    55.0689     0.9787  56.265  < 2e-16 ***
breedHolstein   3.1184     0.7259   4.296 2.21e-05 ***
tx_group250g    1.3509     0.9862   1.370  0.17158    
tx_group300g    2.8492     0.9871   2.886  0.00412 ** 
tx_group350g    5.5358     1.0135   5.462 8.47e-08 ***
sexM           -1.1942     0.8356  -1.429  0.15379    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.613 on 385 degrees of freedom
Multiple R-squared:  0.188, Adjusted R-squared:  0.1775 
F-statistic: 17.83 on 5 and 385 DF,  p-value: 6.612e-16
Show the code 
emmeans(lm_adjust_dose_sex,~breed)
 breed            emmean    SE  df lower.CL upper.CL
 Angus x Holstein   56.9 0.521 385     55.9     57.9
 Holstein           60.0 0.526 385     59.0     61.1

Results are averaged over the levels of: tx_group, sex 
Confidence level used: 0.95

Univariable model

Show the code 
lm <- lm(sts_grams_l ~ breed, data = data)
summary(lm)

Call:
lm(formula = sts_grams_l ~ breed, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.9244  -3.0257  -0.9244   3.0756  18.9743 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    56.9244     0.5410 105.226  < 2e-16 ***
breedHolstein   4.1014     0.6486   6.323 7.03e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.901 on 389 degrees of freedom
Multiple R-squared:  0.09321,   Adjusted R-squared:  0.09088 
F-statistic: 39.99 on 1 and 389 DF,  p-value: 7.03e-10

Maximal model

Show the code 
lm <- lm(sts_grams_l ~ tx_group + breed + igg + time_to_feed + vol + sex, data = data)
summary(lm)

Call:
lm(formula = sts_grams_l ~ tx_group + breed + igg + time_to_feed + 
    vol + sex, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.5728  -2.7550  -0.4452   2.4897  14.2778 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   41.580742   1.732703  23.998  < 2e-16 ***
tx_group250g  -1.172431   0.785426  -1.493    0.136    
tx_group300g  -0.778562   0.972079  -0.801    0.424    
tx_group350g  -0.493685   1.129652  -0.437    0.662    
breedHolstein  2.939544   0.513065   5.729 2.05e-08 ***
igg            0.563175   0.028845  19.524  < 2e-16 ***
time_to_feed  -0.004065   0.006857  -0.593    0.554    
vol           -0.178485   0.589877  -0.303    0.762    
sexM          -0.392308   0.591400  -0.663    0.508    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.961 on 382 degrees of freedom
Multiple R-squared:  0.5988,    Adjusted R-squared:  0.5904 
F-statistic: 71.26 on 8 and 382 DF,  p-value: < 2.2e-16

Comparing models

Show the code 
lm_univariable <- lm(sts_grams_l ~ breed, data = data)

lm_adjust_all <- lm(sts_grams_l ~ tx_group + breed + igg + time_to_feed + vol + sex, data = data)

lm_adjust_dose_sex <- lm(sts_grams_l ~ breed + tx_group + sex, data = data)

lm_adjust_dose <- lm(sts_grams_l ~ breed + tx_group, data = data)

rbind(
    broom::tidy(lm_adjust_dose_sex, conf.int=T) %>% 
    select(term,estimate,conf.low,conf.high) %>%
    filter(term == "breedHolstein") %>%
    mutate(model = "DAG-based model - Adjusted for dose, sex"),
    
    broom::tidy(lm_univariable, conf.int=T) %>% 
    select(term,estimate,conf.low,conf.high) %>%
    filter(term == "breedHolstein") %>%
    mutate(model = "No adjustments"),
    
    broom::tidy(lm_adjust_all, conf.int=T) %>% 
    select(term,estimate,conf.low,conf.high) %>%
    filter(term == "breedHolstein") %>%
    mutate(model = "Adjusted for serum IgG, dose, sex, time to feed, volume"),
  
    broom::tidy(lm_adjust_dose, conf.int=T) %>% 
    select(term,estimate,conf.low,conf.high) %>%
    filter(term == "breedHolstein") %>%
    mutate(model = "Adjusted for dose")
)
term estimate conf.low conf.high model
breedHolstein 3.118389 1.691185 4.545594 DAG-based model - Adjusted for dose, sex
breedHolstein 4.101365 2.826158 5.376573 No adjustments
breedHolstein 2.939544 1.930758 3.948330 Adjusted for serum IgG, dose, sex, time to feed, volume
breedHolstein 3.647880 2.418945 4.876815 Adjusted for dose

Does breed impact ability of STS to predict serum IgG?

STS vs IgG stratified by breed

Show the code 
# data_test <- data %>% mutate(
#   sts_grams_l = case_when(
#     breed == "Holstein" & sts_grams_l > 70 ~ 50,
#     T ~ sts_grams_l
#   )
# )

scatter_figure <- ggplot(data, aes(x = igg, y = sts_grams_l, color = breed)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, aes(group = breed), formula = y ~ x, size = 0.5) +
  stat_cor(
    aes(label = paste(..r.label..)), 
    method = "pearson", 
    label.x.npc = "left", 
    label.y.npc = "top", 
    show.legend = FALSE
  ) +
  labs(
    x = "IgG (g/L)",
    y = "Serum total solids (g/L)",
    color = NULL
  ) +
  theme_minimal() +
  theme(
    legend.position = c(0.98, 0.02),
    legend.justification = c("right", "bottom"),
    legend.background = element_rect(fill = "white", color = "black"),
    legend.title = element_blank(),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1),
    text = element_text(family = "Times New Roman"),
    strip.text = element_text(size = 12, face = "bold")
  )

scatter_figure

Show the code 
ggsave(scatter_figure,filename = "short communication figure 2.tiff",dpi = 300,width = 4, height = 4)
Show the code 
ggplot(data, aes(y = sts_grams_l, x = igg)) +
  geom_point() +                          # Add points
  geom_smooth(method = "lm", se = FALSE, formula = y ~ x,linetype = 1,size = 0.5) +
  labs(
    title = "",
    x = "IgG (g/L)",
    y = "Serum total solids (g/L)",
    color = "Breed"
  ) +
  facet_wrap(~ breed, ncol = 2) +
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1),
    text = element_text(family = "Times New Roman"),
    strip.text = element_text(size = 12, face = "bold")
  )

Compare slopes by breed

Show the code 
model <- lm(sts_grams_l ~ igg * breed, data = data)
summary(model)

Call:
lm(formula = sts_grams_l ~ igg * breed, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.9088  -2.7851  -0.3085   2.5322  14.4353 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)       40.042367   1.551137  25.815   <2e-16 ***
igg                0.558946   0.049931  11.194   <2e-16 ***
breedHolstein      2.913811   1.841759   1.582    0.114    
igg:breedHolstein  0.003819   0.058254   0.066    0.948    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.96 on 387 degrees of freedom
Multiple R-squared:  0.5938,    Adjusted R-squared:  0.5907 
F-statistic: 188.6 on 3 and 387 DF,  p-value: < 2.2e-16
Show the code 
emtrends(model, ~ breed, var = "igg")
 breed            igg.trend     SE  df lower.CL upper.CL
 Angus x Holstein     0.559 0.0499 387    0.461    0.657
 Holstein             0.563 0.0300 387    0.504    0.622

Confidence level used: 0.95

Comparing Area Under the Curve (AUC) for STS as a predictor of IgG >25g/L for angus and holstein calves.

Holstein

Show the code 
roc_holstein <- roc(data$igg[data$breed=="Holstein"] >= 25,
                    data$sts_grams_l[data$breed=="Holstein"])
roc_holstein

Call:
roc.default(response = data$igg[data$breed == "Holstein"] >=     25, predictor = data$sts_grams_l[data$breed == "Holstein"])

Data: data$sts_grams_l[data$breed == "Holstein"] in 54 controls (data$igg[data$breed == "Holstein"] >= 25 FALSE) < 218 cases (data$igg[data$breed == "Holstein"] >= 25 TRUE).
Area under the curve: 0.8626

Angus x Holstein

Show the code 
roc_angus <- roc(data$igg[data$breed=="Angus x Holstein"] >= 25,
                    data$sts_grams_l[data$breed=="Angus x Holstein"])
roc_angus

Call:
roc.default(response = data$igg[data$breed == "Angus x Holstein"] >=     25, predictor = data$sts_grams_l[data$breed == "Angus x Holstein"])

Data: data$sts_grams_l[data$breed == "Angus x Holstein"] in 31 controls (data$igg[data$breed == "Angus x Holstein"] >= 25 FALSE) < 88 cases (data$igg[data$breed == "Angus x Holstein"] >= 25 TRUE).
Area under the curve: 0.8605

Generating a ROC curve to compare performance of STS to predict IgG > 25 in holstein and angus cross calves

Show the code 
roc_figure <- ggplot(data %>% filter(!is.na(breed)), 
       aes(d = igg >= 25, m = sts_grams_l, color = breed)) +
  
  geom_roc(cutoffs.at = c(59,62),
           labels = TRUE,
           size = 0.6,
           pointsize = 0.3,
           labelsize = 3,
           label.family = "Times New Roman",  # set font
           label.fontface = "plain") +
  
  style_roc(ylab = "Sensitivity", xlab = "1 - Specificity") +
  
  facet_wrap(~ breed) +
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1),
    text = element_text(family = "Times New Roman"),
    strip.text = element_text(size = 12, face = "bold")
  )

ggsave(roc_figure,filename = "short communication figure 3.tiff",dpi = 300,width = 7, height = 4)
roc_figure

Holsteins

Show the code 
tibble(
  threshold = roc_holstein$thresholds,
  sensitivity = roc_holstein$sensitivities,
  specificity = roc_holstein$specificities,
  sum_se_sp = sensitivity + specificity
)
threshold sensitivity specificity sum_se_sp
-Inf 1.0000000 0.0000000 1.000000
46.0 1.0000000 0.0185185 1.018519
47.5 1.0000000 0.0370370 1.037037
49.0 1.0000000 0.1296296 1.129630
51.0 1.0000000 0.1666667 1.166667
52.5 0.9770642 0.3888889 1.365953
53.5 0.9770642 0.4259259 1.402990
54.5 0.9495413 0.5185185 1.468060
55.5 0.9403670 0.5555556 1.495923
56.5 0.8899083 0.6851852 1.575093
57.5 0.8807339 0.7222222 1.602956
58.5 0.7844037 0.7777778 1.562181
59.5 0.7568807 0.8148148 1.571695
60.5 0.5366972 0.9074074 1.444105
61.5 0.5137615 0.9074074 1.421169
62.5 0.3807339 0.9629630 1.343697
63.5 0.3715596 0.9629630 1.334523
65.0 0.3027523 0.9629630 1.265715
66.5 0.2293578 0.9629630 1.192321
67.5 0.2201835 0.9629630 1.183146
68.5 0.1238532 0.9814815 1.105335
69.5 0.1146789 0.9814815 1.096160
70.5 0.0688073 1.0000000 1.068807
71.5 0.0596330 1.0000000 1.059633
73.0 0.0458716 1.0000000 1.045872
74.5 0.0275229 1.0000000 1.027523
75.5 0.0183486 1.0000000 1.018349
77.0 0.0091743 1.0000000 1.009174
79.0 0.0045872 1.0000000 1.004587
Inf 0.0000000 1.0000000 1.000000

  • This receiver operating characteristic curve demonstrates the performance of the serum total solids test across multiple thresholds.

  • ROC curve for Holsteins to assess the performance of STS to predict IgG >25 g/L.

    • 6.2 g/dl STS appears to have a sensitivity of 51% and specificity of 89%.

    • 5.9 g/dl STS has sensitivity of 77% and specificity of 75%.

    • AUC is 0.85.

Holstein x Angus

Show the code 
tibble(
  threshold = roc_angus$thresholds,
  sensitivity = roc_angus$sensitivities,
  specificity = roc_angus$specificities,
  sum_se_sp = sensitivity + specificity
)
threshold sensitivity specificity sum_se_sp
-Inf 1.0000000 0.0000000 1.000000
41.0 1.0000000 0.0322581 1.032258
43.5 1.0000000 0.0645161 1.064516
46.5 1.0000000 0.0967742 1.096774
49.0 1.0000000 0.2580645 1.258064
50.5 0.9659091 0.4516129 1.417522
51.5 0.9545455 0.5483871 1.502933
53.0 0.9090909 0.6129032 1.521994
54.5 0.8068182 0.6774194 1.484238
55.5 0.7613636 0.7096774 1.471041
56.5 0.6704545 0.8709677 1.541422
57.5 0.6590909 0.8709677 1.530059
58.5 0.5454545 0.9354839 1.480938
59.5 0.4772727 0.9354839 1.412757
60.5 0.2727273 1.0000000 1.272727
61.5 0.2500000 1.0000000 1.250000
62.5 0.1363636 1.0000000 1.136364
63.5 0.1250000 1.0000000 1.125000
65.0 0.0795455 1.0000000 1.079546
67.0 0.0568182 1.0000000 1.056818
70.0 0.0113636 1.0000000 1.011364
Inf 0.0000000 1.0000000 1.000000

  • ROC curve for Angus X to assess the performance of STS to predict IgG >25 g/L.

    • 5.9 g/dl STS has a sensitivity of 58% and specificity of 95%.

      • This cut-off shows similar performance to the current established threshold of 6.2g/dl for Holsteins. However, this is an interim recommendation only as there is not enough data to make this conclusion and more research required to verify this.

    • 5.6 g/dl STS has a sensitivity of 77% and specificity of 67%.

    • AUC is 0.85.

  • Conclusion:

    • STS is a good test for predicting serum IgG of >25g/L as the area under the curve is high 0.85.

      • STS test for estimating serum IgG perform equivalently for both breeds as the AUC is the same.

    • A lower cutpoint may be necessary for predicting Angus X animals with serum IgG >25g/L. However, further research is required to confirm what value cut-off should be used.

Comparing the use of a STS cut-point of 62g/L to predict if a calf has IgG > 25g/L

Show the code 
df_holstein <- data %>% filter(
  breed == "Holstein") %>%
  group_by(sts_gte_6.2,igg_gt_25) %>%
  summarise(n = n())

pivot_wider(df_holstein, id_cols = c(sts_gte_6.2), names_from = igg_gt_25, values_from = n)
sts_gte_6.2 IgG >= 25 IgG < 25
STS >= 6.2 112 5
STS < 6.2 106 49
Show the code 
output <- epi.tests(df_holstein, method = "exact", digits = 2, 
   conf.level = 0.95)

tab1 <- output$detail %>% 
  filter(statistic == "se" | statistic == "sp") %>% 
  select(statistic,est) %>%
  rename(holstein_cut_at_62 = est)
tab1
statistic holstein_cut_at_62
se 0.5137615
sp 0.9074074
Show the code 
df_angus <- data %>% filter(
  breed == "Angus x Holstein") %>%
  group_by(sts_gte_6.2,igg_gt_25) %>%
  summarise(n = n())

pivot_wider(df_angus, id_cols = c(sts_gte_6.2), names_from = igg_gt_25, values_from = n)
sts_gte_6.2 IgG >= 25 IgG < 25
STS >= 6.2 22 NA
STS < 6.2 66 31
Show the code 
output <- epi.tests(df_angus, method = "exact", digits = 2, 
   conf.level = 0.95)

tab2 <- output$detail %>% 
  filter(statistic == "se" | statistic == "sp") %>% 
  select(statistic,est) %>%
  rename(angus_holstein_cut_at_62 = est)

tab2
statistic angus_holstein_cut_at_62
se 0.25
sp 1.00
Show the code 
df_angus <- data %>% filter(
  breed == "Angus x Holstein") %>%
  mutate(
    sts_gte_5.9 = case_when(
      sts >= 5.9 ~ 1,
      !is.na(sts) ~ 0
    ) %>% factor(levels = c(1,0),
                 labels = c("STS >= 5.9","STS < 5.9"))) %>%
  group_by(sts_gte_5.9,igg_gt_25) %>%
  summarise(n = n())

pivot_wider(df_angus, id_cols = c(sts_gte_5.9), names_from = igg_gt_25, values_from = n)
sts_gte_5.9 IgG >= 25 IgG < 25
STS >= 5.9 48 2
STS < 5.9 40 29
Show the code 
output <- epi.tests(df_angus, method = "exact", digits = 2, 
   conf.level = 0.95)

tab3 <- output$detail %>% 
  filter(statistic == "se" | statistic == "sp") %>% 
  select(statistic,est) %>%
  rename(angus_holstein_cut_at_59 = est)

Final table summarising cut points for predicting IgG > 25g/L

Show the code 
tab1 %>% left_join(tab2) %>% left_join(tab3)
statistic holstein_cut_at_62 angus_holstein_cut_at_62 angus_holstein_cut_at_59
se 0.5137615 0.25 0.5454545
sp 0.9074074 1.00 0.9354839