Randomized clinical trial evaluating the practical and biological limits of colostral IgG dosing in dairy calves

Author

Dini Hapukotuwa and Sam Rowe

Published

December 22, 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

  1. Determine the effect of oral IgG dose (200, 250, 300, 350g IgG) on serum IgG

    1. A secondary objective was to determine if the effect of oral IgG on serum IgG is modified by other variables.

  2. Determine the effect of oral IgG dose (200, 250, 300, 350g IgG) on apparent efficiency of absorption

  3. To identify the practically feasible dose of colostral IgG when considering IgG supply (i.e., IgG harvested and stored), demand (number of calves to be fed), and practical constraints on dose (colostral IgG concentration and volume limits per calf)

  4. Identify if volume impacts IgG independently of total dose administered

Causal diagram (DAG)

Based on this causal diagram, colostrum IgG concentration and volume should be included as covariates as they could provide a backdoor path (confounding bias). Other variables do not need to be included.

Import data and start packages

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

weights <- read_excel(here("Data in","weights.xlsx"))
weather_data_2 <- read_excel(here("Data in","weather_data_2.xlsx"), 
    col_types = c("date", "numeric", "numeric"))

data <- read_excel(here("Data in","Colostrum Study 2 Data V1 20241126.xlsx"), 
    col_types = c("text", "numeric", "date", 
        "date", "date", "numeric", 
        "numeric", "numeric", "text", "text", 
        "text", "numeric", "text", "numeric", 
        "text", "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "text")) %>% clean_names()

data <- data %>% mutate(
  time_born = as.POSIXct(
    paste(enrol_date,format(time_born, "%H:%M:%S")),
    format = "%Y-%m-%d %H:%M:%S"),
  
  time_fed = as.POSIXct(
    paste(enrol_date,format(time_fed, "%H:%M:%S")),
    format = "%Y-%m-%d %H:%M:%S"),
  
  time_fed = case_when(
    time_fed < time_born ~ time_fed + days(1),
    T ~ time_fed
  ),
  
  time_to_feed = as.numeric(time_fed - time_born, units = "mins"),
  
  tx_group = factor(tx_administered,
                           labels = c("200g",
                                      "250g",
                                      "300g",
                                      "350g")),
  
  tx_assigned = factor(tx_assigned,
                           labels = c("200g",
                                      "250g",
                                      "300g",
                                      "350g")),
  
  tx_per_protocol = case_when(
    tx_group == tx_assigned ~ tx_group,
    T ~ NA
  ),
  
  breed = case_when(
    breed == "H" ~ "Holstein",
    breed == "A" ~ "Angus x Holstein"
  )
) %>%
  rename(
    igg = serum_ig_g,
    sts = serum_total_protein_g_l
  ) %>%
  
  mutate(
    igg_gt_10 = case_when(
      igg < 10 ~ 0,
      igg >= 10 ~ 1
    ),
    
    igg_gt_25 = case_when(
      igg >= 25 ~ 1,
      igg < 25 ~ 0
    ) %>% factor(levels = c(1,0),
                 labels = c("IgG >= 25","IgG < 25")),
    
    sts_gte_6.2 = case_when(
      sts >= 6.2 ~ 1,
      !is.na(sts) ~ 0
    ) %>% factor(levels = c(1,0),
                 labels = c("STS >= 6.2","STS < 6.2")),
    
    igg_category = factor(case_when(
      igg < 10 ~ "Poor (< 10g/L)",
      igg >= 10 & igg < 18 ~ "Fair (10 to <18g/L)",
      igg >= 18 & igg < 25 ~ "Good (18 to <25g/L)",
      igg >= 25 ~ "Excellent (>= 25g/L)"
    ), levels = c("Excellent (>= 25g/L)", "Good (18 to <25g/L)", "Fair (10 to <18g/L)","Poor (< 10g/L)")),
    
    sts_category = factor(case_when(
      sts < 5.0999 ~ "Poor (< 5.1/L)",
      sts >= 5.1 & sts < 5.7999 ~ "Fair (5.1 to <5.7g/dL)",
      sts >= 5.8 & sts < 6.1999 ~ "Good (5.8 to 6.1g/dL)",
      sts >= 6.2 ~ "Excellent (>= 6.2g/dL)"
    ), levels = c("Excellent (>= 6.2g/dL)", "Good (5.8 to 6.1g/dL)", "Fair (5.1 to <5.7g/dL)","Poor (< 5.1g/dL)")),
    
    sts_grams_l = sts*10,
    

  
  
    time_to_feed_cat = ntile(time_to_feed, 4) %>% 
      factor(labels = c("Q1 (Lowest)", "Q2", "Q3", "Q4 (Highest)"))
  )
  

events <- read_csv(here("Data in/","dcomp data.CSV"), 
    col_types = cols(HERDID = col_skip(), 
        BDAT = col_date(format = "%d/%m/%y"), 
        Date = col_date(format = "%d/%m/%y"), 
        FDAT = col_skip(), LACT = col_skip())) %>% clean_names

all_events <- events
events <- events %>% filter(id %in% data$calf_id)

#events %>% tabyl(event)


data <- data %>% 
  filter(!is.na(enrol_date))

day <- data %>% group_by(enrol_date) %>%
  summarise(
    born = n()
  ) %>% rename(date = enrol_date)


trial_dates <- tibble(
  start = min(day$date),
  end = max(day$date) - day(1)
)

batch <- read_excel(here("Data in","Colostrum Harvest Recording Sheet.xlsx"), 
    col_types = c("date", "date", "numeric", 
        "numeric", "numeric", "numeric", 
        "text")) %>% clean_names

batch <- batch %>% 
  filter(between(date,trial_dates$start,trial_dates$end)) %>%
  mutate(
    igg_kg = pooled_colostrum_volume * rid_reading_g_l / 1000,
    dose_at_4l = rid_reading_g_l * 4,
    dose_at_3.8l = rid_reading_g_l * 3.8
  )

batch_day <- batch %>% 
  group_by(date) %>% summarise(
    igg_kg = sum(igg_kg),
    concentration_limited_dose = min(dose_at_4l)
  )

day <- day %>%
  left_join(batch_day,by = "date") %>%
  filter(!is.na(igg_kg)) %>%
  mutate(
    igg_per_calf = igg_kg / born * 1000
  )

day$date <- as.Date(day$date)


dose <- read_excel(here("Data in","Dose Inventory Data V1 20241127.xlsx")) %>% clean_names()

dose <- dose %>%
  rename(tx_id = dose_id, vol = volume_l)

data <- data %>% left_join(dose %>% select(tx_id,batch_id,vol))

data <- data %>% mutate(
  batch_id = case_when(
    is.na(batch_id) ~ "Missing",
    T ~ batch_id),
    
  vol_cat = case_when(
    vol < 2.5 ~ "<2.5",
    vol >= 2.5 & vol < 3.0 ~ "2.5-2.999",
    vol >= 3.0 & vol < 3.5 ~ "3.0-3.499",
    vol >= 3.5 & vol <= 4.0 ~ "3.5-4.0",
    TRUE ~ NA_character_
  ),
  
  concentration = tx_administered/vol
)

#data %>% tabyl(batch_id)

#Apparent efficiency of absorption
#AEA = serum IgG (g/L) × birth weight (kg) × serum volume (%)/IgG fed (g) × 100

# Holsteins were 38.9kg and angus calves were 42.1kg, but the difference was not statistically significant, so we will assume the same weight for all calves (39.6kg)

average_weight <- weights$weight %>%mean

data <- data %>% mutate(
  assumed_bw = case_when(
    breed == "Holstein" ~ 39.6,
    breed == "Angus x Holstein" ~ 39.6
  ),
  total_igg_absorbed = igg * assumed_bw * 0.07,
  aea = total_igg_absorbed / tx_administered
)

Describe dataset

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

Data Frame Summary

data

Dimensions: 424 x 39
Duplicates: 0
No Variable Stats / Values Freqs (% of Valid) Graph Missing
1 sample_id [character]
1. 1
2. 10
3. 100
4. 101
5. 102
6. 103
7. 104
8. 105
9. 106
10. 107
[ 412 others ]
1 ( 0.2% )
1 ( 0.2% )
1 ( 0.2% )
1 ( 0.2% )
1 ( 0.2% )
1 ( 0.2% )
1 ( 0.2% )
1 ( 0.2% )
1 ( 0.2% )
1 ( 0.2% )
412 ( 97.6% )
2 (0.5%)
2 calf_id [numeric]
Mean (sd) : 422948.6 (425010.2)
min ≤ med ≤ max:
108720 ≤ 108931.5 ≤ 996765
IQR (CV) : 887832.5 (1)
424 distinct values 0 (0.0%)
3 enrol_date [POSIXct, POSIXt]
min : 2024-10-24
med : 2024-11-03
max : 2024-11-10
range : 17d
18 distinct values 0 (0.0%)
4 time_born [POSIXct, POSIXt]
min : 2024-10-24 11:00:00
med : 2024-11-03 14:07:30
max : 2024-11-10 23:45:00
range : 17d 12H 45M 0S
404 distinct values 0 (0.0%)
5 time_fed [POSIXct, POSIXt]
min : 2024-10-24 11:40:00
med : 2024-11-03 13:05:00
max : 2024-11-11 00:30:00
range : 17d 12H 50M 0S
403 distinct values 7 (1.7%)
6 sts [numeric]
Mean (sd) : 6 (0.6)
min ≤ med ≤ max:
4 ≤ 6 ≤ 8
IQR (CV) : 0.6 (0.1)
33 distinct values 3 (0.7%)
7 tx_assigned [factor]
1. 200g
2. 250g
3. 300g
4. 350g
45 ( 10.6% )
132 ( 31.1% )
133 ( 31.4% )
114 ( 26.9% )
0 (0.0%)
8 tx_administered [numeric]
Mean (sd) : 286.7 (49)
min ≤ med ≤ max:
200 ≤ 300 ≤ 350
IQR (CV) : 100 (0.2)
200 : 46 ( 11.2% )
250 : 129 ( 31.3% )
300 : 126 ( 30.6% )
350 : 111 ( 26.9% )
12 (2.8%)
9 tx_id [character]
1. EXCLUDED
2. P108
3. B1
4. B10
5. B100
6. B101
7. B102
8. B104
9. B105
10. B106
[ 402 others ]
4 ( 1.0% )
2 ( 0.5% )
1 ( 0.2% )
1 ( 0.2% )
1 ( 0.2% )
1 ( 0.2% )
1 ( 0.2% )
1 ( 0.2% )
1 ( 0.2% )
1 ( 0.2% )
402 ( 96.6% )
8 (1.9%)
10 breed [character]
1. Angus x Holstein
2. Holstein
130 ( 30.7% )
294 ( 69.3% )
0 (0.0%)
11 sex [character]
1. F
2. M
339 ( 80.0% )
85 ( 20.0% )
0 (0.0%)
12 excluded [numeric]
Min : 0
Mean : 0.1
Max : 1
0 : 400 ( 94.3% )
1 : 24 ( 5.7% )
0 (0.0%)
13 excluded_reason [character]
1. Dose ID not recorded
2. Full breech, difficult ca
3. Premature. Not given dose
4. Too large
5. Too small
8 ( 33.3% )
1 ( 4.2% )
2 ( 8.3% )
2 ( 8.3% )
11 ( 45.8% )
400 (94.3%)
14 lost_to_follow_up [numeric]
Min : 0
Mean : 0
Max : 1
0 : 422 ( 99.5% )
1 : 2 ( 0.5% )
0 (0.0%)
15 lost_to_follow_up_reason [character]
1. Died 3/11/24
2. Euthanased 10/11/24
1 ( 50.0% )
1 ( 50.0% )
422 (99.5%)
16 dilution_factor [numeric]
Mean (sd) : 1.9 (0.4)
min ≤ med ≤ max:
1 ≤ 2 ≤ 3
IQR (CV) : 0 (0.2)
1 : 49 ( 11.8% )
2 : 346 ( 83.6% )
3 : 19 ( 4.6% )
10 (2.4%)
17 sample_1_read_1 [numeric]
Mean (sd) : 6 (0.6)
min ≤ med ≤ max:
3.7 ≤ 6 ≤ 7.5
IQR (CV) : 0.7 (0.1)
195 distinct values 10 (2.4%)
18 sample_1_read_2 [numeric]
Mean (sd) : 6 (0.6)
min ≤ med ≤ max:
3.7 ≤ 6 ≤ 7.5
IQR (CV) : 0.7 (0.1)
190 distinct values 10 (2.4%)
19 sample_2_read_1 [numeric]
Mean (sd) : 6 (0.6)
min ≤ med ≤ max:
3.9 ≤ 6 ≤ 7.5
IQR (CV) : 0.7 (0.1)
197 distinct values 11 (2.6%)
20 sample_2_read_2 [numeric]
Mean (sd) : 6 (0.6)
min ≤ med ≤ max:
3.8 ≤ 6 ≤ 7.5
IQR (CV) : 0.7 (0.1)
195 distinct values 11 (2.6%)
21 igg [numeric]
Mean (sd) : 31.6 (7.8)
min ≤ med ≤ max:
7.6 ≤ 32.1 ≤ 53.9
IQR (CV) : 10.3 (0.2)
386 distinct values 10 (2.4%)
22 comments [character] 1. <-- 350G DOSE COMMENCED
1 ( 100.0% )
423 (99.8%)
23 time_to_feed [numeric]
Mean (sd) : 61.6 (48.5)
min ≤ med ≤ max:
0 ≤ 55 ≤ 850
IQR (CV) : 35 (0.8)
38 distinct values 7 (1.7%)
24 tx_group [factor]
1. 200g
2. 250g
3. 300g
4. 350g
46 ( 11.2% )
129 ( 31.3% )
126 ( 30.6% )
111 ( 26.9% )
12 (2.8%)
25 tx_per_protocol [factor]
1. 200g
2. 250g
3. 300g
4. 350g
45 ( 11.1% )
127 ( 31.3% )
125 ( 30.8% )
109 ( 26.8% )
18 (4.2%)
26 igg_gt_10 [numeric]
Min : 0
Mean : 1
Max : 1
0 : 2 ( 0.5% )
1 : 412 ( 99.5% )
10 (2.4%)
27 igg_gt_25 [factor]
1. IgG >= 25
2. IgG < 25
327 ( 79.0% )
87 ( 21.0% )
10 (2.4%)
28 sts_gte_6.2 [factor]
1. STS >= 6.2
2. STS < 6.2
148 ( 35.2% )
273 ( 64.8% )
3 (0.7%)
29 igg_category [factor]
1. Excellent (>= 25g/L)
2. Good (18 to <25g/L)
3. Fair (10 to <18g/L)
4. Poor (< 10g/L)
327 ( 79.0% )
70 ( 16.9% )
15 ( 3.6% )
2 ( 0.5% )
10 (2.4%)
30 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)
148 ( 37.6% )
140 ( 35.5% )
106 ( 26.9% )
0 ( 0.0% )
30 (7.1%)
31 sts_grams_l [numeric]
Mean (sd) : 59.7 (6.1)
min ≤ med ≤ max:
40 ≤ 60 ≤ 80
IQR (CV) : 6 (0.1)
33 distinct values 3 (0.7%)
32 time_to_feed_cat [factor]
1. Q1 (Lowest)
2. Q2
3. Q3
4. Q4 (Highest)
105 ( 25.2% )
104 ( 24.9% )
104 ( 24.9% )
104 ( 24.9% )
7 (1.7%)
33 batch_id [character]
1. D1
2. G1
3. F1
4. C1
5. H
6. Missing
7. A1
8. B
9. B3
10. B5
[ 40 others ]
16 ( 3.8% )
15 ( 3.5% )
13 ( 3.1% )
12 ( 2.8% )
12 ( 2.8% )
12 ( 2.8% )
11 ( 2.6% )
11 ( 2.6% )
11 ( 2.6% )
11 ( 2.6% )
300 ( 70.8% )
0 (0.0%)
34 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 12 (2.8%)
35 vol_cat [character]
1. <2.5
2. 2.5-2.999
3. 3.0-3.499
4. 3.5-4.0
53 ( 12.9% )
123 ( 29.9% )
99 ( 24.0% )
137 ( 33.3% )
12 (2.8%)
36 concentration [numeric]
Mean (sd) : 93.9 (11.6)
min ≤ med ≤ max:
71.4 ≤ 92.6 ≤ 120.7
IQR (CV) : 13.8 (0.1)
40 distinct values 12 (2.8%)
37 assumed_bw [numeric] 1 distinct value
39.60 : 424 ( 100.0% )
0 (0.0%)
38 total_igg_absorbed [numeric]
Mean (sd) : 87.7 (21.6)
min ≤ med ≤ max:
21.1 ≤ 89 ≤ 149.5
IQR (CV) : 28.6 (0.2)
386 distinct values 10 (2.4%)
39 aea [numeric]
Mean (sd) : 0.3 (0.1)
min ≤ med ≤ max:
0.1 ≤ 0.3 ≤ 0.5
IQR (CV) : 0.1 (0.2)
393 distinct values 21 (5.0%)

Generated by summarytools 1.1.4 (R version 4.5.1)
2025-12-22

Check for duplicated IDs

Show the code 
data %>% get_dupes(calf_id)
calf_id dupe_count sample_id enrol_date time_born time_fed sts tx_assigned tx_administered tx_id breed sex excluded excluded_reason lost_to_follow_up lost_to_follow_up_reason dilution_factor sample_1_read_1 sample_1_read_2 sample_2_read_1 sample_2_read_2 igg comments time_to_feed tx_group tx_per_protocol igg_gt_10 igg_gt_25 sts_gte_6.2 igg_category sts_category sts_grams_l time_to_feed_cat batch_id vol vol_cat concentration assumed_bw total_igg_absorbed aea

Descriptive statistics

Comparison of treatment groups at enrollment

Show the code 
table1(~ breed + sex + time_to_feed + vol + factor(excluded) + excluded_reason +  factor(lost_to_follow_up) + lost_to_follow_up_reason | factor(tx_assigned), data=data)
200g
(N=45)		 250g
(N=132)		 300g
(N=133)		 350g
(N=114)		 Overall
(N=424)
breed
Angus x Holstein 21 (46.7%) 39 (29.5%) 44 (33.1%) 26 (22.8%) 130 (30.7%)
Holstein 24 (53.3%) 93 (70.5%) 89 (66.9%) 88 (77.2%) 294 (69.3%)
sex
F 35 (77.8%) 105 (79.5%) 110 (82.7%) 89 (78.1%) 339 (80.0%)
M 10 (22.2%) 27 (20.5%) 23 (17.3%) 25 (21.9%) 85 (20.0%)
time_to_feed
Mean (SD) 56.0 (28.2) 66.9 (76.8) 57.5 (27.9) 62.5 (24.9) 61.6 (48.5)
Median [Min, Max] 50.0 [0, 155] 50.0 [0, 850] 50.0 [0, 150] 60.0 [5.00, 130] 55.0 [0, 850]
Missing 0 (0%) 1 (0.8%) 4 (3.0%) 2 (1.8%) 7 (1.7%)
vol
Mean (SD) 2.23 (0.239) 2.78 (0.397) 3.28 (0.387) 3.55 (0.299) 3.08 (0.551)
Median [Min, Max] 2.20 [1.80, 2.60] 2.70 [1.80, 3.90] 3.30 [2.50, 3.90] 3.60 [2.90, 3.90] 3.10 [1.80, 3.90]
Missing 0 (0%) 2 (1.5%) 5 (3.8%) 5 (4.4%) 12 (2.8%)
factor(excluded)
0 44 (97.8%) 126 (95.5%) 125 (94.0%) 105 (92.1%) 400 (94.3%)
1 1 (2.2%) 6 (4.5%) 8 (6.0%) 9 (7.9%) 24 (5.7%)
excluded_reason
Too small 1 (2.2%) 3 (2.3%) 2 (1.5%) 5 (4.4%) 11 (2.6%)
Dose ID not recorded 0 (0%) 2 (1.5%) 3 (2.3%) 3 (2.6%) 8 (1.9%)
Full breech, difficult calving 5/7, too small 0 (0%) 1 (0.8%) 0 (0%) 0 (0%) 1 (0.2%)
Premature. Not given dose. 0 (0%) 0 (0%) 2 (1.5%) 0 (0%) 2 (0.5%)
Too large 0 (0%) 0 (0%) 1 (0.8%) 1 (0.9%) 2 (0.5%)
Missing 44 (97.8%) 126 (95.5%) 125 (94.0%) 105 (92.1%) 400 (94.3%)
factor(lost_to_follow_up)
0 45 (100%) 132 (100%) 131 (98.5%) 114 (100%) 422 (99.5%)
1 0 (0%) 0 (0%) 2 (1.5%) 0 (0%) 2 (0.5%)
lost_to_follow_up_reason
Died 3/11/24 0 (0%) 0 (0%) 1 (0.8%) 0 (0%) 1 (0.2%)
Euthanased 10/11/24 0 (0%) 0 (0%) 1 (0.8%) 0 (0%) 1 (0.2%)
Missing 45 (100%) 132 (100%) 131 (98.5%) 114 (100%) 422 (99.5%)

  • The number of excluded animals is low and similar between each group. This indicates that there is less bias between groups.

  • There is a higher proportion of Angus in the 200g group compared to the other groups. This may create bias if Angus was different in their response to oral IgG. When results were adjusted for breed below, it was found that there was no significant difference between both, and therefore there is no bias in the numbers above.

Comparison of treatment groups after exclusions, losses to follow-up and calves changing treatment group (n = 6)

Show the code 
data <- data %>% mutate(
  originally_assiged_to_other_dose = case_when(
    tx_group == tx_assigned ~ 0,
    T ~ 1
  ) %>% as_factor
)

data <- data %>% filter(excluded==0,
                        lost_to_follow_up == 0,
                        !is.na(tx_group))
Show the code 
table1(~ breed + sex + time_to_feed + vol + originally_assiged_to_other_dose | factor(tx_group), data=data)
200g
(N=45)		 250g
(N=125)		 300g
(N=122)		 350g
(N=107)		 Overall
(N=399)
breed
Angus x Holstein 21 (46.7%) 37 (29.6%) 42 (34.4%) 26 (24.3%) 126 (31.6%)
Holstein 24 (53.3%) 88 (70.4%) 80 (65.6%) 81 (75.7%) 273 (68.4%)
sex
F 35 (77.8%) 98 (78.4%) 101 (82.8%) 84 (78.5%) 318 (79.7%)
M 10 (22.2%) 27 (21.6%) 21 (17.2%) 23 (21.5%) 81 (20.3%)
time_to_feed
Mean (SD) 56.7 (28.3) 61.9 (34.4) 57.0 (28.3) 61.4 (24.6) 59.7 (29.5)
Median [Min, Max] 50.0 [0, 155] 50.0 [0, 195] 50.0 [0, 150] 60.0 [5.00, 130] 50.0 [0, 195]
vol
Mean (SD) 2.22 (0.247) 2.78 (0.368) 3.30 (0.378) 3.56 (0.298) 3.08 (0.552)
Median [Min, Max] 2.20 [1.80, 2.60] 2.70 [2.10, 3.50] 3.30 [2.50, 3.90] 3.60 [2.90, 3.90] 3.10 [1.80, 3.90]
originally_assiged_to_other_dose
0 44 (97.8%) 123 (98.4%) 121 (99.2%) 105 (98.1%) 393 (98.5%)
1 1 (2.2%) 2 (1.6%) 1 (0.8%) 2 (1.9%) 6 (1.5%)

Objective 1 - Effect of oral IgG on serum IgG in Holstein calves

Plots

Box and whisker plot showing serum IgG for calves randomized to 4 doses of colostrum.

Show the code 
figure <- ggplot(data %>% filter(breed == "Holstein"), 
                 aes(x = tx_group, y = igg)) +  # Removed fill from aes()
  
  # Add shaded regions
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -5, ymax = 10), 
            fill = "pink", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 10, ymax = 18), 
            fill = "#FFF2CC", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 18, ymax = 25), 
            fill = "#C3D69B", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 25, ymax = Inf), 
            fill = "#A2D9CE", alpha = 0.5, inherit.aes = FALSE) +
  
  # 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 = "Oral IgG dose",
    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

Show the code 
ggsave(
  filename = "figure 1.tiff",    # File name
  plot = figure,              # ggplot object to save
  device = "tiff",             # Output format
  width = 4,                   # Width in inches
  height = 4,                  # Height in inches
  dpi = 200                    # Resolution in dots per inch
)

  • This box and whisker plot is comparing the IgG levels across treatment groups (250, 300 and 350g), relative to baseline (200g) for Holstein animals.

  • There appears to be a dose dependent effect of oral IgG on serum IgG in Holstein calves.

Show the code 
prop_data <- data %>%
  filter(breed == "Holstein") %>%
  filter(!is.na(tx_group), !is.na(igg_category)) %>%
  group_by(tx_group, igg_category) %>%
  summarize(count = n(), .groups = "drop") %>%
  group_by(tx_group) %>%
  mutate(proportion = count / sum(count))

figure <- ggplot(prop_data, aes(x = tx_group, y = proportion, fill = igg_category)) +
  geom_bar(stat = "identity", position = "stack") +
  
  # Add text labels for proportions
  geom_text(
    aes(label = scales::percent(proportion, accuracy = 1)), # Format proportion as percentage
    position = position_stack(vjust = 0.5), # Place text in the middle of the stacked bars
    size = 3, # Adjust text size
    family = "Times New Roman" # Ensure consistent font
  ) +
  
  labs(
    title = "",
    x = "Oral IgG dose",
    y = "Proportion",
    fill = "IgG Levels"
  ) +
  scale_y_continuous(labels = scales::percent) + # Show percentages
  theme_minimal() +
  scale_fill_manual(values = c('#A2D9CE', '#C3D69B', '#FFF2CC', "pink")) +
  theme(
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1),
    text = element_text(family = "Times New Roman")
  )

figure

Show the code 
ggsave(
  filename = "figure 2.tiff",    # File name
  plot = figure,              # ggplot object to save
  device = "tiff",             # Output format
  width = 6,                   # Width in inches
  height = 5,                  # Height in inches
  dpi = 200                    # Resolution in dots per inch
)

  • This graph shows the proportion of Holstein calves from each dose group that fall into the new proposed categories by Godden and colleagues (2019) for passive immunity in calves.

  • Lombard and collaborators (2020) proposed that >40% calves should be able to achieve an excellent level of passive immunity when given 300g if feeding once.

  • It is evident from this graph that this is also possible at a dose of 250g as 75% Holstein calves in the study were able to achieve a serum IgG >=25g/L IgG.

Statistical models

Main effects linear regression

Model includes concentration and volume as covariates (multi-colinearity. model rejected)

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

Call:
lm(formula = igg ~ tx_group + concentration + vol, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-23.2936  -4.7233  -0.0212   4.6516  19.5518 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept)    45.6277    31.1431   1.465   0.1437  
tx_group250g    7.5217     3.7355   2.014   0.0447 *
tx_group300g   12.2245     6.9846   1.750   0.0809 .
tx_group350g   18.6577     9.8676   1.891   0.0594 .
concentration  -0.1084     0.1916  -0.566   0.5719  
vol            -4.8773     6.1978  -0.787   0.4318  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.064 on 386 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.1989,    Adjusted R-squared:  0.1885 
F-statistic: 19.16 on 5 and 386 DF,  p-value: < 2.2e-16
Show the code 
broom::tidy(lm, conf.int=T) %>% 
  select(term,estimate,conf.low,conf.high)
term estimate conf.low conf.high
(Intercept) 45.6277205 -15.6036391 106.8590801
tx_group250g 7.5216737 0.1771622 14.8661851
tx_group300g 12.2245125 -1.5081703 25.9571952
tx_group350g 18.6577267 -0.7431776 38.0586310
concentration -0.1083789 -0.4850355 0.2682777
vol -4.8773485 -17.0629477 7.3082507

Model diagnostics

Show the code 
check_model(lm)

Show the code 
check_outliers(lm)
OK: No outliers detected.
- Based on the following method and threshold: cook (0.8).
- For variable: (Whole model)
Show the code 
check_collinearity(lm)
Term VIF VIF_CI_low VIF_CI_high SE_factor Tolerance Tolerance_CI_low Tolerance_CI_high
tx_group 82.54447 68.10437 100.09191 2.086642 0.0121147 0.0099908 0.0146833
concentration 37.91441 31.33392 45.92243 6.157468 0.0263752 0.0217759 0.0319143
vol 92.55015 76.34800 112.23627 9.620299 0.0108050 0.0089098 0.0130979

Decision - model is unstable due to multi-colinearity (volume and concentration). Will include volume only as more data support the link between volume and AEA.

Model includes volume as covariate

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

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

Residuals:
     Min       1Q   Median       3Q      Max 
-23.1171  -4.6460  -0.0524   4.7329  19.6052 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    28.068      2.549  11.009  < 2e-16 ***
tx_group250g    5.553      1.359   4.086 5.34e-05 ***
tx_group300g    8.388      1.674   5.011 8.24e-07 ***
tx_group350g   13.178      1.884   6.993 1.19e-11 ***
vol            -1.422      1.047  -1.358    0.175    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.058 on 387 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.1982,    Adjusted R-squared:  0.1899 
F-statistic: 23.92 on 4 and 387 DF,  p-value: < 2.2e-16
Show the code 
broom::tidy(lm, conf.int=T) %>% 
  select(term,estimate,conf.low,conf.high)
term estimate conf.low conf.high
(Intercept) 28.068299 23.055741 33.0808572
tx_group250g 5.553479 2.881255 8.2257038
tx_group300g 8.388434 5.097434 11.6794340
tx_group350g 13.178252 9.473223 16.8832809
vol -1.421567 -3.480291 0.6371571

Model diagnostics

Show the code 
check_model(lm)

Show the code 
check_outliers(lm)
OK: No outliers detected.
- Based on the following method and threshold: cook (0.7).
- For variable: (Whole model)
Show the code 
check_collinearity(lm)
Term VIF VIF_CI_low VIF_CI_high SE_factor Tolerance Tolerance_CI_low Tolerance_CI_high
tx_group 2.646368 2.2781 3.120749 1.176093 0.3778764 0.3204359 0.4389624
vol 2.646368 2.2781 3.120749 1.626766 0.3778764 0.3204359 0.4389624

Estimated marginal means

Show the code 
pairwise_results <- emmeans(lm, pairwise ~ tx_group,infer = c(TRUE, TRUE))
summary(pairwise_results)
$emmeans
 tx_group emmean    SE  df lower.CL upper.CL t.ratio p.value
 200g       23.7 1.380 387     21.0     26.4  17.108  <.0001
 250g       29.2 0.713 387     27.8     30.6  41.018  <.0001
 300g       32.1 0.686 387     30.7     33.4  46.740  <.0001
 350g       36.9 0.853 387     35.2     38.5  43.233  <.0001

Confidence level used: 0.95 

$contrasts
 contrast    estimate    SE  df lower.CL upper.CL t.ratio p.value
 200g - 250g    -5.55 1.360 387    -9.06  -2.0466  -4.086  0.0003
 200g - 300g    -8.39 1.670 387   -12.71  -4.0695  -5.011  <.0001
 200g - 350g   -13.18 1.880 387   -18.04  -8.3160  -6.993  <.0001
 250g - 300g    -2.83 1.060 387    -5.57  -0.0959  -2.671  0.0393
 250g - 350g    -7.62 1.250 387   -10.85  -4.4029  -6.106  <.0001
 300g - 350g    -4.79 0.984 387    -7.33  -2.2512  -4.868  <.0001

Confidence level used: 0.95 
Conf-level adjustment: tukey method for comparing a family of 4 estimates 
P value adjustment: tukey method for comparing a family of 4 estimates

Conclusion:

  • There is a significant difference between the means of the 250g vs 350g and 300g and 350g groups. This shows that there is a linear effect of oral IgG dose and serum IgG. This is consistent with Godden’s principle of a dose-dependant effect of oral IgG and serum IgG.
Show the code 
figure <- ggplot(data, 
                 aes(x = tx_group, y = igg)) +  # Removed fill from aes()
  
  # Add shaded regions
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -5, ymax = 10), 
            fill = "pink", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 10, ymax = 18), 
            fill = "#FFF2CC", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 18, ymax = 25), 
            fill = "#C3D69B", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 25, ymax = Inf), 
            fill = "#A2D9CE", alpha = 0.5, inherit.aes = FALSE) +
  
  # 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 = "Oral IgG dose",
    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

Show the code 
ggsave(
  filename = "figure 1a.tiff",    # File name
  plot = figure,              # ggplot object to save
  device = "tiff",             # Output format
  width = 5,                   # Width in inches
  height = 4,                  # Height in inches
  dpi = 200                    # Resolution in dots per inch
)

Table Objective 1

Final model includes volume as a covariate

Show the code 
lm_all_breed <- lm(igg ~ tx_group + vol, 
                  data = data %>%
                    mutate(tx_group = fct_relevel(tx_group,c("250g","200g","300g","350g")))
                    )


lm_all_breed_tab <- lm_all_breed %>%
    broom::tidy(conf.int=T) %>% 
  mutate(
    difference = case_when(
      term == "(Intercept)" ~ "Ref",
      T ~ paste0(round(estimate,1),
                        " (",
                        round(conf.low,1),
                        " to ",
                        round(conf.high,1),
                        ")"
  )),
  term = case_when(
    term == "(Intercept)" ~ "250g",
    term != "(Intercept)" ~ str_remove(term, "tx_group"))
  )%>%
  select(term,difference)

emmean_all_breed_tab <- emmeans(lm_all_breed, ~tx_group) %>% 
  as_tibble %>%
    mutate(
      term = tx_group,
      emmean = paste0(round(emmean,1),
                        " (",
                        round(lower.CL,1),
                        " to ",
                        round(upper.CL,1),
                        ")"
  )) %>%
  select(term,emmean)


lm_holstein <- lm(igg ~tx_group + vol, 
                  data = data %>% 
                    filter(breed == "Holstein") %>%
                    mutate(tx_group = fct_relevel(tx_group,c("250g","200g","300g","350g")))
                    )


lm_holstein_tab <- lm_holstein %>%
    broom::tidy(conf.int=T) %>% 
  mutate(
    difference = case_when(
      term == "(Intercept)" ~ "Ref",
      T ~ paste0(round(estimate,1),
                        " (",
                        round(conf.low,1),
                        " to ",
                        round(conf.high,1),
                        ")"
  )),
  term = case_when(
    term == "(Intercept)" ~ "250g",
    term != "(Intercept)" ~ str_remove(term, "tx_group"))
  )%>%
  select(term,difference)

emmean_holstein_tab <- emmeans(lm_holstein, ~tx_group) %>% 
  as_tibble %>%
    mutate(
      term = tx_group,
      emmean = paste0(round(emmean,1),
                        " (",
                        round(lower.CL,1),
                        " to ",
                        round(upper.CL,1),
                        ")"
  )) %>%
  select(term,emmean)


lm_angus <- lm(igg ~ tx_group + vol,data = data %>% 
                    filter(breed == "Angus x Holstein") %>%
                    mutate(tx_group = fct_relevel(tx_group,c("250g","200g","300g","350g"))))

lm_angus_tab <- lm_angus %>%
    broom::tidy(conf.int=T) %>% 
  mutate(
    difference = case_when(
      term == "(Intercept)" ~ "Ref",
      T ~ paste0(round(estimate,1),
                        " (",
                        round(conf.low,1),
                        " to ",
                        round(conf.high,1),
                        ")"
  )),
  term = case_when(
    term == "(Intercept)" ~ "250g",
    term != "(Intercept)" ~ str_remove(term, "tx_group"))
  )%>%
  select(term,difference)

emmean_angus_tab <- emmeans(lm_angus, ~tx_group) %>% 
  as_tibble %>%
    mutate(
      term = tx_group,
      emmean = paste0(round(emmean,1),
                        " (",
                        round(lower.CL,1),
                        " to ",
                        round(upper.CL,1),
                        ")"
  )) %>%
  select(term,emmean)

emmean_all_breed_tab %>% 
  left_join(lm_all_breed_tab, by = "term") %>%
  left_join(emmean_holstein_tab, by ="term",suffix = c("_overall","_holstein")) %>%
  left_join(lm_holstein_tab, by = "term", suffix = c("_overall","_holstein")) %>% 
  left_join(emmean_angus_tab %>% rename(emmean_angus = emmean), by ="term") %>%
  left_join(lm_angus_tab %>% rename(difference_angus = difference), by = "term") %>%
  rename(tx = term) %>%
  arrange(tx)
tx emmean_overall difference_overall emmean_holstein difference_holstein emmean_angus difference_angus
200g 23.7 (21 to 26.4) -5.6 (-8.2 to -2.9) 22.9 (19.2 to 26.6) -6 (-9.7 to -2.4) 25.4 (21.4 to 29.4) -4.6 (-8.6 to -0.6)
250g 29.2 (27.8 to 30.6) Ref 29 (27.2 to 30.7) Ref 30 (27.5 to 32.5) Ref
300g 32.1 (30.7 to 33.4) 2.8 (0.7 to 4.9) 32.9 (31.2 to 34.5) 3.9 (1.4 to 6.4) 30 (27.6 to 32.5) 0 (-3.8 to 3.9)
350g 36.9 (35.2 to 38.5) 7.6 (5.2 to 10.1) 37.4 (35.5 to 39.4) 8.5 (5.5 to 11.4) 35 (31.7 to 38.3) 5 (0.4 to 9.6)

Sensitivity analysis

Intention to treat vs per protocol vs as-treated

Show the code 
lm_as_treated <- lm(igg ~ tx_group,data = data) %>% 
  broom::tidy(conf.int = T) %>%
  mutate(
    tx = c("200g","250g","300g","350g"),
    estimate = paste0(round(estimate,1)," 95% CI: (",round(conf.low,1),", ",round(conf.high,1),")"),
    estimate = case_when(
      tx == "200g" ~ "Ref",
      T ~ as.character(estimate)
    )
  ) %>%
  select(tx,estimate)

lm_per_protocol <- lm(igg ~ tx_assigned, data = data) %>%
  broom::tidy(conf.int = T) %>%
  mutate(
    tx = c("200g","250g","300g","350g"),
    estimate = paste0(round(estimate,1)," 95% CI: (",round(conf.low,1),", ",round(conf.high,1),")"),
    estimate = case_when(
      tx == "200g" ~ "Ref",
      T ~ as.character(estimate)
    )
  ) %>%
  select(tx,estimate)


lm_intention_to_treat <- lm(igg ~ tx_per_protocol, data = data) %>%
  broom::tidy(conf.int = T) %>%
  mutate(
    tx = c("200g","250g","300g","350g"),
    estimate = paste0(round(estimate,1)," 95% CI: (",round(conf.low,1),", ",round(conf.high,1),")"),
    estimate = case_when(
      tx == "200g" ~ "Ref",
      T ~ as.character(estimate)
    )
  ) %>%
  select(tx,estimate)


lm_intention_to_treat %>% 
  left_join(lm_per_protocol,by = "tx",suffix = c("_intention_to_treat","_per_protocol")) %>%
  left_join(lm_as_treated %>% rename(estimate_as_treated = estimate), by = "tx")
tx estimate_intention_to_treat estimate_per_protocol estimate_as_treated
200g Ref Ref Ref
250g 4.9 95% CI: (2.5, 7.4) 5 95% CI: (2.6, 7.5) 4.8 95% CI: (2.3, 7.2)
300g 6.9 95% CI: (4.4, 9.4) 6.8 95% CI: (4.3, 9.2) 6.9 95% CI: (4.4, 9.3)
350g 11.3 95% CI: (8.8, 13.8) 11.3 95% CI: (8.8, 13.8) 11.3 95% CI: (8.8, 13.7)

Comparison of different covariate structures

estimate_univariable = No covariates

estimate_maximal = Model controls for sex, breed, time to feeding, volume, concentration

estimate_adjust_vol_conc = Model controls for volume and concentration

estimate_adjust_vol = Model controls for volume only

estimate_adjust_conc = Model controls for concentration only

estimate_adjust_sex_breed_time_volume = Model controls for sex, breed, time to feeding and volume

Show the code 
lm_univariable <- lm(igg ~ tx_group,data = data) %>% 
  broom::tidy(conf.int = T) %>%
  slice_head(n=4) %>%
  mutate(
    tx = c("200g","250g","300g","350g"),
    estimate = paste0(round(estimate,1)," 95% CI: (",round(conf.low,1),", ",round(conf.high,1),")"),
    estimate = case_when(
      tx == "200g" ~ "Ref",
      T ~ as.character(estimate)
    )
  ) %>%
  rename(estimate_univariable = estimate) %>%
  select(tx,estimate_univariable)

lm(igg ~ tx_group,data = data) %>% summary

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

Residuals:
     Min       1Q   Median       3Q      Max 
-22.9797  -4.6840  -0.1127   4.6102  19.4214 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    24.916      1.053  23.655  < 2e-16 ***
tx_group250g    4.767      1.231   3.873 0.000126 ***
tx_group300g    6.855      1.237   5.544 5.48e-08 ***
tx_group350g   11.273      1.259   8.954  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.066 on 388 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.1944,    Adjusted R-squared:  0.1882 
F-statistic: 31.21 on 3 and 388 DF,  p-value: < 2.2e-16
Show the code 
lm_maximal <- lm(igg ~ tx_group + sex + breed + time_to_feed + vol + concentration,data = data) %>% 
  broom::tidy(conf.int = T) %>%
  slice_head(n=4) %>%
  mutate(
    tx = c("200g","250g","300g","350g"),
    estimate = paste0(round(estimate,1)," 95% CI: (",round(conf.low,1),", ",round(conf.high,1),")"),
    estimate = case_when(
      tx == "200g" ~ "Ref",
      T ~ as.character(estimate)
    )
  ) %>%
  rename(estimate_maximal = estimate) %>%
  select(tx,estimate_maximal)

lm(igg ~ tx_group + sex + breed + time_to_feed + vol + concentration,data = data) %>% summary

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

Residuals:
     Min       1Q   Median       3Q      Max 
-21.8347  -4.7677   0.0907   4.7099  18.7628 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept)   44.92278   31.00293   1.449   0.1482  
tx_group250g   7.28103    3.72438   1.955   0.0513 .
tx_group300g  11.51533    6.95483   1.656   0.0986 .
tx_group350g  17.87119    9.82640   1.819   0.0697 .
sexM          -1.29987    1.04648  -1.242   0.2149  
breedHolstein  0.38650    0.90943   0.425   0.6711  
time_to_feed  -0.02575    0.01209  -2.130   0.0338 *
vol           -4.35849    6.16722  -0.707   0.4802  
concentration -0.09635    0.19059  -0.506   0.6135  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.023 on 383 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.2142,    Adjusted R-squared:  0.1978 
F-statistic: 13.05 on 8 and 383 DF,  p-value: < 2.2e-16
Show the code 
lm_adjusting_for_volume_concentration <- lm(igg ~ tx_group + vol + concentration,data = data) %>% 
  broom::tidy(conf.int = T) %>%
  slice_head(n=4) %>%
  mutate(
    tx = c("200g","250g","300g","350g"),
    estimate = paste0(round(estimate,1)," 95% CI: (",round(conf.low,1),", ",round(conf.high,1),")"),
    estimate = case_when(
      tx == "200g" ~ "Ref",
      T ~ as.character(estimate)
    )
  ) %>%
  rename(estimate_adjust_vol_conc = estimate) %>%
  select(tx,estimate_adjust_vol_conc)

lm(igg ~ tx_group + vol + concentration,data = data) %>% summary

Call:
lm(formula = igg ~ tx_group + vol + concentration, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-23.2936  -4.7233  -0.0212   4.6516  19.5518 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept)    45.6277    31.1431   1.465   0.1437  
tx_group250g    7.5217     3.7355   2.014   0.0447 *
tx_group300g   12.2245     6.9846   1.750   0.0809 .
tx_group350g   18.6577     9.8676   1.891   0.0594 .
vol            -4.8773     6.1978  -0.787   0.4318  
concentration  -0.1084     0.1916  -0.566   0.5719  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.064 on 386 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.1989,    Adjusted R-squared:  0.1885 
F-statistic: 19.16 on 5 and 386 DF,  p-value: < 2.2e-16
Show the code 
lm_adjusting_for_volume_only <- lm(igg ~ tx_group + vol,data = data) %>% 
  broom::tidy(conf.int = T) %>%
  slice_head(n=4) %>%
  mutate(
    tx = c("200g","250g","300g","350g"),
    estimate = paste0(round(estimate,1)," 95% CI: (",round(conf.low,1),", ",round(conf.high,1),")"),
    estimate = case_when(
      tx == "200g" ~ "Ref",
      T ~ as.character(estimate)
    )
  ) %>%
  rename(estimate_adjust_vol = estimate) %>%
  select(tx,estimate_adjust_vol)

lm(igg ~ tx_group + vol,data = data) %>% summary

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

Residuals:
     Min       1Q   Median       3Q      Max 
-23.1171  -4.6460  -0.0524   4.7329  19.6052 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    28.068      2.549  11.009  < 2e-16 ***
tx_group250g    5.553      1.359   4.086 5.34e-05 ***
tx_group300g    8.388      1.674   5.011 8.24e-07 ***
tx_group350g   13.178      1.884   6.993 1.19e-11 ***
vol            -1.422      1.047  -1.358    0.175    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.058 on 387 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.1982,    Adjusted R-squared:  0.1899 
F-statistic: 23.92 on 4 and 387 DF,  p-value: < 2.2e-16
Show the code 
lm_adjusting_for_conc_only <- lm(igg ~ tx_group + concentration,data = data) %>% 
  broom::tidy(conf.int = T) %>%
  slice_head(n=4) %>%
  mutate(
    tx = c("200g","250g","300g","350g"),
    estimate = paste0(round(estimate,1)," 95% CI: (",round(conf.low,1),", ",round(conf.high,1),")"),
    estimate = case_when(
      tx == "200g" ~ "Ref",
      T ~ as.character(estimate)
    )
  ) %>%
  rename(estimate_adjust_conc = estimate) %>%
  select(tx,estimate_adjust_conc)

lm(igg ~ tx_group + concentration,data = data) %>% summary

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

Residuals:
     Min       1Q   Median       3Q      Max 
-23.0381  -4.5991  -0.0886   4.7512  19.6476 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   21.24441    3.13803   6.770 4.79e-11 ***
tx_group250g   4.74616    1.23024   3.858 0.000134 ***
tx_group300g   6.81478    1.23606   5.513 6.45e-08 ***
tx_group350g  10.95843    1.28328   8.539 3.12e-16 ***
concentration  0.04021    0.03238   1.242 0.215048    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.061 on 387 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.1976,    Adjusted R-squared:  0.1893 
F-statistic: 23.82 on 4 and 387 DF,  p-value: < 2.2e-16
Show the code 
estimate_adjust_sex_breed_time_volume <- lm(igg ~ tx_group + sex + breed + time_to_feed + vol,data = data) %>% 
  broom::tidy(conf.int = T) %>%
  slice_head(n=4) %>%
  mutate(
    tx = c("200g","250g","300g","350g"),
    estimate = paste0(round(estimate,1)," 95% CI: (",round(conf.low,1),", ",round(conf.high,1),")"),
    estimate = case_when(
      tx == "200g" ~ "Ref",
      T ~ as.character(estimate)
    )
  ) %>%
  rename(estimate_adjust_sex_breed_time_volume = estimate) %>%
  select(tx,estimate_adjust_sex_breed_time_volume)
# lm_interaction_model <- lm(igg ~ vol*concentration, data = data)
# summary(lm_interaction_model)

lm(igg ~ tx_group + sex + breed + time_to_feed + vol,data = data) %>% summary

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

Residuals:
     Min       1Q   Median       3Q      Max 
-21.6754  -4.7492   0.1258   4.7825  18.7664 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   29.30768    2.68010  10.935  < 2e-16 ***
tx_group250g   5.52887    1.36231   4.058 5.99e-05 ***
tx_group300g   8.10243    1.67101   4.849 1.81e-06 ***
tx_group350g  12.99606    1.88781   6.884 2.38e-11 ***
sexM          -1.29738    1.04545  -1.241   0.2154    
breedHolstein  0.40039    0.90813   0.441   0.6595    
time_to_feed  -0.02585    0.01207  -2.141   0.0329 *  
vol           -1.28553    1.04239  -1.233   0.2182    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.017 on 384 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.2137,    Adjusted R-squared:  0.1994 
F-statistic: 14.91 on 7 and 384 DF,  p-value: < 2.2e-16
Show the code 
lm_univariable %>% 
  left_join(lm_maximal,by = "tx") %>%
  left_join(lm_adjusting_for_volume_concentration,by = "tx") %>%
  left_join(lm_adjusting_for_volume_only,by = "tx") %>%
  left_join(lm_adjusting_for_conc_only,by = "tx") %>%
  left_join(estimate_adjust_sex_breed_time_volume,by = "tx")
tx estimate_univariable estimate_maximal estimate_adjust_vol_conc estimate_adjust_vol estimate_adjust_conc estimate_adjust_sex_breed_time_volume
200g Ref Ref Ref Ref Ref Ref
250g 4.8 95% CI: (2.3, 7.2) 7.3 95% CI: (0, 14.6) 7.5 95% CI: (0.2, 14.9) 5.6 95% CI: (2.9, 8.2) 4.7 95% CI: (2.3, 7.2) 5.5 95% CI: (2.9, 8.2)
300g 6.9 95% CI: (4.4, 9.3) 11.5 95% CI: (-2.2, 25.2) 12.2 95% CI: (-1.5, 26) 8.4 95% CI: (5.1, 11.7) 6.8 95% CI: (4.4, 9.2) 8.1 95% CI: (4.8, 11.4)
350g 11.3 95% CI: (8.8, 13.7) 17.9 95% CI: (-1.4, 37.2) 18.7 95% CI: (-0.7, 38.1) 13.2 95% CI: (9.5, 16.9) 11 95% CI: (8.4, 13.5) 13 95% CI: (9.3, 16.7)

These findings indicate that the linear relationship between dose and serum IgG is robust to covariate selection.

Mixed model vs fixed effect model

Show the code 
library(lme4)
lmm <- lmer(igg ~ tx_group + (1|batch_id), data = data)

broom.mixed::tidy(lmm, conf.int=T)
effect group term estimate std.error statistic conf.low conf.high
fixed NA (Intercept) 25.854202 1.145124 22.577644 23.609800 28.098604
fixed NA tx_group250g 3.934679 1.220385 3.224130 1.542769 6.326589
fixed NA tx_group300g 6.163747 1.218187 5.059772 3.776145 8.551349
fixed NA tx_group350g 9.895257 1.365982 7.244059 7.217981 12.572533
ran_pars batch_id sd__(Intercept) 2.835955 NA NA NA NA
ran_pars Residual sd__Observation 6.521435 NA NA NA NA
Show the code 
icc(lmm)
ICC_adjusted ICC_unadjusted optional
0.1590343 0.1338655 FALSE
Show the code 
lm_fixed_effect <- lm(igg ~ tx_group + vol,data = data) %>% 
  broom::tidy(conf.int = T) %>%
  slice_head(n=4) %>%
  mutate(
    tx = c("200g","250g","300g","350g"),
    estimate = paste0(round(estimate,1)," 95% CI: (",round(conf.low,1),", ",round(conf.high,1),")"),
    estimate = case_when(
      tx == "200g" ~ "Ref",
      T ~ as.character(estimate)
    )
  ) %>%
  rename(estimate_fixed_effect = estimate) %>%
  select(tx,estimate_fixed_effect)

lmm <- lmer(igg ~ tx_group + vol + (1|batch_id), data = data) %>% 
  broom.mixed::tidy(conf.int = T) %>%
  slice_head(n=4) %>%
  mutate(
    tx = c("200g","250g","300g","350g"),
    estimate = paste0(round(estimate,1)," 95% CI: (",round(conf.low,1),", ",round(conf.high,1),")"),
    estimate = case_when(
      tx == "200g" ~ "Ref",
      T ~ as.character(estimate)
    )
  ) %>%
  rename(estimate_mixed_model = estimate) %>%
  select(tx,estimate_mixed_model)

lm_fixed_effect %>% 
  left_join(lmm,by = "tx")
tx estimate_fixed_effect estimate_mixed_model
200g Ref Ref
250g 5.6 95% CI: (2.9, 8.2) 5 95% CI: (2.1, 7.9)
300g 8.4 95% CI: (5.1, 11.7) 8.1 95% CI: (4.2, 12.1)
350g 13.2 95% CI: (9.5, 16.9) 12.5 95% CI: (7.5, 17.4)

These findings indicate that the linear relationship between dose and serum IgG is robust to potential clustering of calves within batches of colostrum.

Objective 1b - Investigating factors that may impact IgG absorption and therefore modify the effect of IgG dose on serum IgG.

Breed (Holstein vs Holstein x Angus)

Breed as an effect modifier

Box and whisker plot showing serum IgG for calves randomized to 4 doses of colostrum, stratified by breed

Show the code 
figure <- ggplot(data %>% filter(!is.na(tx_group)), 
                 aes(x = tx_group, y = igg)) +  # Removed fill from aes()
  
  # Add shaded regions
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -5, ymax = 10), 
            fill = "pink", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 10, ymax = 18), 
            fill = "#FFF2CC", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 18, ymax = 25), 
            fill = "#C3D69B", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 25, ymax = Inf), 
            fill = "#A2D9CE", alpha = 0.5, inherit.aes = FALSE) +
  
  # 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 = "Oral IgG dose",
    y = "Serum IgG (g/L)",
    fill = "Treatment"
  ) +
  
  # Custom theme
  theme_minimal() +
  facet_wrap(~ breed) +
  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

Show the code 
ggsave(
  filename = "figure holstein and angus.tiff",    # File name
  plot = figure,              # ggplot object to save
  device = "tiff",             # Output format
  width = 5,                   # Width in inches
  height = 3,                  # Height in inches
  dpi = 200                    # Resolution in dots per inch
)
Show the code 
prop_data <- data %>%
  filter(!is.na(tx_group), !is.na(igg_category)) %>%
  group_by(tx_group, igg_category, breed) %>%
  summarize(count = n(), .groups = "drop") %>%
  group_by(tx_group,breed) %>%
  mutate(proportion = count / sum(count))

figure <- ggplot(prop_data, aes(x = tx_group, y = proportion, fill = igg_category)) +
  geom_bar(stat = "identity", position = "stack") +
  
  # Add text labels for proportions
  geom_text(
    aes(label = scales::percent(proportion, accuracy = 1)), # Format proportion as percentage
    position = position_stack(vjust = 0.5), # Place text in the middle of the stacked bars
    size = 3, # Adjust text size
    family = "Times New Roman" # Ensure consistent font
  ) +
  
  labs(
    title = "",
    x = "Oral IgG dose",
    y = "Proportion",
    fill = "IgG Levels"
  ) +
  scale_y_continuous(labels = scales::percent) + # Show percentages
  theme_minimal() +
  facet_wrap(~ breed) +
  scale_fill_manual(values = c('#A2D9CE', '#C3D69B', '#FFF2CC', "pink")) +
  theme(
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1),
    text = element_text(family = "Times New Roman")
  )

figure

Show the code 
ggsave(
  filename = "figure angus holstein groups.tiff",    # File name
  plot = figure,              # ggplot object to save
  device = "tiff",             # Output format
  width = 8,                   # Width in inches
  height = 5,                  # Height in inches
  dpi = 200                    # Resolution in dots per inch
)
Show the code 
lm <- lm(igg ~ tx_group*breed + vol, data = data)

summary(lm)

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

Residuals:
    Min      1Q  Median      3Q     Max 
-23.820  -4.827   0.034   4.775  19.631 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                27.90179    2.82934   9.862  < 2e-16 ***
tx_group250g                5.46452    2.01331   2.714  0.00694 ** 
tx_group300g                6.65560    2.24264   2.968  0.00319 ** 
tx_group350g               11.95730    2.52244   4.740 3.01e-06 ***
breedHolstein              -0.23146    2.11271  -0.110  0.91282    
vol                        -1.29083    1.05318  -1.226  0.22108    
tx_group250g:breedHolstein  0.08221    2.54388   0.032  0.97424    
tx_group300g:breedHolstein  2.41561    2.52290   0.957  0.33893    
tx_group350g:breedHolstein  1.42700    2.67584   0.533  0.59414    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.067 on 383 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.2045,    Adjusted R-squared:  0.1879 
F-statistic: 12.31 on 8 and 383 DF,  p-value: 1.026e-15

Type 2 P-values for assessing significance of an interaction term

Show the code 
car::Anova(lm)
Sum Sq Df F value Pr(>F)
tx_group 2536.59849 3 16.9311836 0.0000000
breed 65.55740 1 1.3127395 0.2526146
vol 75.02029 1 1.5022271 0.2210821
tx_group:breed 86.36277 3 0.5764507 0.6307912
Residuals 19126.78291 383 NA NA

Estimated marginal means

Show the code 
emm <- emmeans(lm, ~ tx_group | breed)
summary(emm)
tx_group breed emmean SE df lower.CL upper.CL
200g Angus x Holstein 23.92888 1.7702482 383 20.44825 27.40950
250g Angus x Holstein 29.39340 1.2507369 383 26.93423 31.85257
300g Angus x Holstein 30.58448 1.1710728 383 28.28194 32.88701
350g Angus x Holstein 35.88617 1.5301584 383 32.87761 38.89474
200g Holstein 23.69742 1.7204628 383 20.31468 27.08015
250g Holstein 29.24415 0.8125525 383 27.64653 30.84177
300g Holstein 32.76863 0.8139977 383 31.16817 34.36909
350g Holstein 37.08172 0.9332789 383 35.24673 38.91671
Show the code 
emm_df <- as.data.frame(emm)

# Plot the results
ggplot(emm_df, aes(x = tx_group, y = emmean, ymin = lower.CL, ymax = upper.CL)) +
  geom_point() + # Points for the means
  geom_errorbar(width = 0.2) + # Error bars for confidence intervals
  facet_wrap(~ breed) + 
  labs(
    title = "Estimated Marginal Means of IgG by Treatment group and breed",
    x = "Treatment Group",
    y = "Estimated Mean igg (with 95% CI)"
  )

Breed as a main effect

Univariable model for effect of volume on IgG.

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

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

Residuals:
     Min       1Q   Median       3Q      Max 
-23.2917  -5.7233   0.6448   4.9886  21.8532 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    30.2033     0.7154  42.219   <2e-16 ***
breedHolstein   1.8790     0.8572   2.192    0.029 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.804 on 390 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.01217,   Adjusted R-squared:  0.009637 
F-statistic: 4.805 on 1 and 390 DF,  p-value: 0.02897
Show the code 
broom::tidy(lm_not_control_for_dose, conf.int=T)
term estimate std.error statistic p.value conf.low conf.high
(Intercept) 30.20329 0.7153875 42.219482 0.000000 28.7967931 31.609790
breedHolstein 1.87905 0.8572413 2.191973 0.028972 0.1936574 3.564442

Model for effect of volume on IgG after controlling for IgG dose

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

summary(lm_control_for_dose)

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

Residuals:
     Min       1Q   Median       3Q      Max 
-23.2853  -4.7318   0.1331   4.8803  19.1561 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    24.4182     1.1329  21.554  < 2e-16 ***
tx_group250g    4.5976     1.2386   3.712 0.000236 ***
tx_group300g    6.7252     1.2407   5.421 1.05e-07 ***
tx_group350g   11.0509     1.2721   8.687  < 2e-16 ***
breedHolstein   0.9325     0.7848   1.188 0.235480    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.062 on 387 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.1973,    Adjusted R-squared:  0.189 
F-statistic: 23.78 on 4 and 387 DF,  p-value: < 2.2e-16
Show the code 
broom::tidy(lm_control_for_dose, conf.int=T)
term estimate std.error statistic p.value conf.low conf.high
(Intercept) 24.4182497 1.1328875 21.553993 0.0000000 22.1908652 26.645634
tx_group250g 4.5975982 1.2385995 3.711933 0.0002359 2.1623719 7.032825
tx_group300g 6.7252044 1.2406529 5.420698 0.0000001 4.2859408 9.164468
tx_group350g 11.0508788 1.2720549 8.687423 0.0000000 8.5498753 13.551882
breedHolstein 0.9324941 0.7847895 1.188209 0.2354796 -0.6104905 2.475479

Time to feeding

Show the code 
figure <- ggplot(data %>% filter(!is.na(tx_group)), 
                 aes(x = tx_group, y = igg)) +  # Removed fill from aes()
  
  # Add shaded regions
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -5, ymax = 10), 
            fill = "pink", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 10, ymax = 18), 
            fill = "#FFF2CC", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 18, ymax = 25), 
            fill = "#C3D69B", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 25, ymax = Inf), 
            fill = "#A2D9CE", alpha = 0.5, inherit.aes = FALSE) +
  
  # 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 = "Stratifed by time to feeding quartile",
    x = "Oral IgG dose",
    y = "Serum IgG (g/L)",
    fill = "Treatment"
  ) +
  
  # Custom theme
  theme_minimal() +
  facet_wrap(~ time_to_feed_cat) +
  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

Time to feeding as an effect modifier

Linear model evaluating effect modification. Statistical interaction between tx and time to feeding

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

summary(lm)

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

Residuals:
     Min       1Q   Median       3Q      Max 
-21.5379  -4.7311   0.1043   4.7294  19.1700 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)               28.879079   3.332775   8.665  < 2e-16 ***
tx_group250g               5.857875   2.759071   2.123   0.0344 *  
tx_group300g               8.733863   3.024292   2.888   0.0041 ** 
tx_group350g              15.103089   3.311455   4.561 6.87e-06 ***
time_to_feed              -0.016053   0.037600  -0.427   0.6697    
vol                       -1.376495   1.048032  -1.313   0.1898    
tx_group250g:time_to_feed -0.003930   0.041877  -0.094   0.9253    
tx_group300g:time_to_feed -0.006984   0.044137  -0.158   0.8744    
tx_group350g:time_to_feed -0.031003   0.046819  -0.662   0.5082    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.047 on 383 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.2089,    Adjusted R-squared:  0.1924 
F-statistic: 12.65 on 8 and 383 DF,  p-value: 3.757e-16
Show the code 
car::Anova(lm)
Sum Sq Df F value Pr(>F)
tx_group 2693.30237 3 18.077639 0.0000000
time_to_feed 220.42335 1 4.438492 0.0357871
vol 85.66860 1 1.725041 0.1898317
tx_group:time_to_feed 37.82364 3 0.253875 0.8585660
Residuals 19020.45609 383 NA NA

Estimated marginal means

Show the code 
data %>% 
  group_by(time_to_feed_cat) %>%
  summarise(
    mean_time_to_feed = mean(time_to_feed,na.rm=T),
    median_time_to_feed = median(time_to_feed)
  )
time_to_feed_cat mean_time_to_feed median_time_to_feed
Q1 (Lowest) 30.5000 30
Q2 45.5534 45
Q3 63.3299 60
Q4 (Highest) 100.3838 95
Show the code 
emm <- emmeans(lm, ~ tx_group | time_to_feed, 
               at = list(time_to_feed = c(30, 45, 60, 95)))

emm
time_to_feed = 30:
 tx_group emmean    SE  df lower.CL upper.CL
 200g       24.2 1.710 383     20.8     27.5
 250g       29.9 0.937 383     28.1     31.7
 300g       32.7 0.927 383     30.9     34.5
 350g       38.3 1.220 383     35.9     40.7

time_to_feed = 45:
 tx_group emmean    SE  df lower.CL upper.CL
 200g       23.9 1.450 383     21.1     26.8
 250g       29.6 0.786 383     28.1     31.1
 300g       32.3 0.739 383     30.9     33.8
 350g       37.6 0.967 383     35.7     39.5

time_to_feed = 60:
 tx_group emmean    SE  df lower.CL upper.CL
 200g       23.7 1.390 383     20.9     26.4
 250g       29.3 0.714 383     27.9     30.7
 300g       32.0 0.689 383     30.6     33.3
 350g       36.9 0.853 383     35.2     38.6

time_to_feed = 95:
 tx_group emmean    SE  df lower.CL upper.CL
 200g       23.1 2.000 383     19.2     27.0
 250g       28.6 0.924 383     26.8     30.4
 300g       31.2 1.110 383     29.0     33.4
 350g       35.3 1.270 383     32.8     37.8

Confidence level used: 0.95
Show the code 
emm_df <- as.data.frame(emm)

# Plot the results
ggplot(emm_df, aes(x = tx_group, y = emmean, ymin = lower.CL, ymax = upper.CL)) +
  geom_point() + # Points for the means
  geom_errorbar(width = 0.2) + # Error bars for confidence intervals
  facet_wrap(~ time_to_feed) + 
  labs(
    title = "Estimated Marginal Means of IgG by Treatment group and time to feeding",
    x = "Treatment Group",
    y = "Estimated Mean igg (with 95% CI)"
  )

Time to feeding as a main effect

Show the code 
figure_scatter <- ggplot(data %>% filter(!is.na(tx_group)),
                         aes(x = time_to_feed, y = igg)) +
  # Background zones
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -5, ymax = 10),
            fill = "pink", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 10, ymax = 18),
            fill = "#FFF2CC", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 18, ymax = 25),
            fill = "#C3D69B", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 25, ymax = Inf),
            fill = "#A2D9CE", alpha = 0.5, inherit.aes = FALSE) +
  
  # Scatter points
  geom_point(size = 1, alpha = 0.8, colour = "black") +
  
  # Optionally add a smoothed line (remove if not desired)
  geom_smooth(method = "lm", se = T, linewidth = 0.5, linetype = "dashed",colour = "black") +
  
  # Labels
  labs(
    title = "Relationship between time to feeding and Serum IgG",
    x = "Time to feeding (min)",
    y = "Serum IgG (g/L)"
  ) +
  
  # Facet by treatment group if relevant
  facet_wrap(~ tx_group) +
  
  coord_cartesian(ylim = c(0, 60)) +
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.border = element_rect(
      color = "black", fill = NA, linewidth = 1
    ),
    text = element_text(family = "Times New Roman")
  )

figure_scatter

Show the code 
ggsave(
  filename = "figure scatter time.tiff",    # File name
  plot = figure_scatter,              # ggplot object to save
  device = "tiff",             # Output format
  width = 6,                   # Width in inches
  height = 5,                  # Height in inches
  dpi = 200                    # Resolution in dots per inch
)

Univariable model for effect of volume on IgG.

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

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

Residuals:
     Min       1Q   Median       3Q      Max 
-23.3501  -5.6577   0.4578   5.1479  21.9856 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  32.84105    0.89147  36.839   <2e-16 ***
time_to_feed -0.02228    0.01339  -1.663   0.0971 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.824 on 390 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.007044,  Adjusted R-squared:  0.004498 
F-statistic: 2.767 on 1 and 390 DF,  p-value: 0.09706
Show the code 
broom::tidy(lm_not_control_for_dose, conf.int=T)
term estimate std.error statistic p.value conf.low conf.high
(Intercept) 32.8410482 0.8914709 36.839172 0.0000000 31.0883582 34.5937382
time_to_feed -0.0222772 0.0133934 -1.663298 0.0970562 -0.0486096 0.0040551

Model for effect of volume on IgG after controlling for IgG dose

Show the code 
lm_control_for_dose <- lm(igg ~ tx_group + time_to_feed, data = data)

summary(lm_control_for_dose)

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

Residuals:
     Min       1Q   Median       3Q      Max 
-21.4920  -4.6019   0.0523   4.8640  18.8429 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  26.38345    1.25273  21.061  < 2e-16 ***
tx_group250g  4.90516    1.22702   3.998 7.66e-05 ***
tx_group300g  6.85231    1.23084   5.567 4.85e-08 ***
tx_group350g 11.39734    1.25450   9.085  < 2e-16 ***
time_to_feed -0.02587    0.01208  -2.141   0.0329 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.033 on 387 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.2038,    Adjusted R-squared:  0.1956 
F-statistic: 24.77 on 4 and 387 DF,  p-value: < 2.2e-16
Show the code 
broom::tidy(lm_control_for_dose, conf.int=T)
term estimate std.error statistic p.value conf.low conf.high
(Intercept) 26.3834503 1.2527301 21.060762 0.0000000 23.9204417 28.8464589
tx_group250g 4.9051560 1.2270208 3.997615 0.0000766 2.4926948 7.3176172
tx_group300g 6.8523114 1.2308374 5.567195 0.0000000 4.4323462 9.2722766
tx_group350g 11.3973387 1.2544971 9.085185 0.0000000 8.9308559 13.8638215
time_to_feed -0.0258732 0.0120848 -2.140965 0.0329011 -0.0496333 -0.0021131

Volume

Volume as an effect modifier.

Show the code 
figure <- ggplot(data %>% filter(!is.na(tx_group)), 
                 aes(x = tx_group, y = igg)) +  # Removed fill from aes()
  
  # Add shaded regions
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -5, ymax = 10), 
            fill = "pink", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 10, ymax = 18), 
            fill = "#FFF2CC", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 18, ymax = 25), 
            fill = "#C3D69B", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 25, ymax = Inf), 
            fill = "#A2D9CE", alpha = 0.5, inherit.aes = FALSE) +
  
  # 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 = "Stratifed by volume fed (L)",
    x = "Oral IgG dose",
    y = "Serum IgG (g/L)",
    fill = "Treatment"
  ) +
  
  # Custom theme
  theme_minimal() +
  facet_wrap(~ vol_cat) +
  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

Statistical interaction between tx and volume fed

Show the code 
lm <- lm(igg ~ tx_group*vol_cat, data = data)
summary(lm)

Call:
lm(formula = igg ~ tx_group * vol_cat, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-22.0128  -4.5394  -0.0687   4.6713  19.4430 

Coefficients: (4 not defined because of singularities)
                              Estimate Std. Error t value Pr(>|t|)    
(Intercept)                     25.400      1.243  20.439  < 2e-16 ***
tx_group250g                     4.428      2.036   2.175  0.03026 *  
tx_group300g                     5.193      3.269   1.588  0.11302    
tx_group350g                     9.561      3.209   2.980  0.00307 ** 
vol_cat2.5-2.999                -1.676      2.312  -0.725  0.46899    
vol_cat3.0-3.499                 2.475      3.332   0.743  0.45804    
vol_cat3.5-4.0                   1.008      2.845   0.354  0.72324    
tx_group250g:vol_cat2.5-2.999    2.377      2.935   0.810  0.41843    
tx_group300g:vol_cat2.5-2.999    5.260      4.079   1.290  0.19798    
tx_group350g:vol_cat2.5-2.999    1.821      4.426   0.411  0.68097    
tx_group250g:vol_cat3.0-3.499   -6.227      3.295  -1.890  0.05952 .  
tx_group300g:vol_cat3.0-3.499   -2.265      2.252  -1.006  0.31529    
tx_group350g:vol_cat3.0-3.499       NA         NA      NA       NA    
tx_group250g:vol_cat3.5-4.0         NA         NA      NA       NA    
tx_group300g:vol_cat3.5-4.0         NA         NA      NA       NA    
tx_group350g:vol_cat3.5-4.0         NA         NA      NA       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.03 on 380 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.219, Adjusted R-squared:  0.1964 
F-statistic: 9.687 on 11 and 380 DF,  p-value: 1.727e-15
Show the code 
car::Anova(lm)
Sum Sq Df F value Pr(>F)
tx_group 3032.1625 3 20.452593 0.0000000
vol_cat 278.1013 3 1.875854 0.1331445
tx_group:vol_cat 313.6797 5 1.269502 0.2763275
Residuals 18778.7389 380 NA NA
Show the code 
emm <- emmeans(lm, ~ tx_group | vol_cat)

summary(emm)
tx_group vol_cat emmean SE df lower.CL upper.CL
200g <2.5 25.39973 1.2427004 380 22.95630 27.84316
250g <2.5 29.82771 1.6127409 380 26.65670 32.99873
300g <2.5 NA NA NA NA NA
350g <2.5 NA NA NA NA NA
200g 2.5-2.999 23.72383 1.9497087 380 19.89026 27.55740
250g 2.5-2.999 30.52929 0.8171947 380 28.92250 32.13608
300g 2.5-2.999 34.17682 1.4658094 380 31.29471 37.05894
350g 2.5-2.999 35.10599 2.3432583 380 30.49861 39.71336
200g 3.0-3.499 NA NA NA NA NA
250g 3.0-3.499 26.07583 1.5340226 380 23.05959 29.09206
300g 3.0-3.499 30.80342 1.0042535 380 28.82883 32.77801
350g 3.0-3.499 37.43623 1.5340226 380 34.41999 40.45246
200g 3.5-4.0 NA NA NA NA NA
250g 3.5-4.0 30.83581 2.3432583 380 26.22843 35.44318
300g 3.5-4.0 31.60076 1.0253980 380 29.58460 33.61693
350g 3.5-4.0 35.96901 0.8117285 380 34.37297 37.56505
Show the code 
#pairs(emm)
Show the code 
emm_df <- as.data.frame(emm)

# Plot the results
ggplot(emm_df, aes(x = tx_group, y = emmean, ymin = lower.CL, ymax = upper.CL)) +
  geom_point() + # Points for the means
  geom_errorbar(width = 0.2) + # Error bars for confidence intervals
  facet_wrap(~ vol_cat) + # Separate panels for each value of vol
  labs(
    title = "Estimated Marginal Means of igg by tx_group and vol",
    x = "Treatment Group (tx_group)",
    y = "Estimated Mean igg (with 95% CI)"
  )

Show the code 
emm <- emmeans(lm, ~ vol_cat | tx_group)

summary(emm)
vol_cat tx_group emmean SE df lower.CL upper.CL
<2.5 200g 25.39973 1.2427004 380 22.95630 27.84316
2.5-2.999 200g 23.72383 1.9497087 380 19.89026 27.55740
3.0-3.499 200g NA NA NA NA NA
3.5-4.0 200g NA NA NA NA NA
<2.5 250g 29.82771 1.6127409 380 26.65670 32.99873
2.5-2.999 250g 30.52929 0.8171947 380 28.92250 32.13608
3.0-3.499 250g 26.07583 1.5340226 380 23.05959 29.09206
3.5-4.0 250g 30.83581 2.3432583 380 26.22843 35.44318
<2.5 300g NA NA NA NA NA
2.5-2.999 300g 34.17682 1.4658094 380 31.29471 37.05894
3.0-3.499 300g 30.80342 1.0042535 380 28.82883 32.77801
3.5-4.0 300g 31.60076 1.0253980 380 29.58460 33.61693
<2.5 350g NA NA NA NA NA
2.5-2.999 350g 35.10599 2.3432583 380 30.49861 39.71336
3.0-3.499 350g 37.43623 1.5340226 380 34.41999 40.45246
3.5-4.0 350g 35.96901 0.8117285 380 34.37297 37.56505
Show the code 
#pairs(emm)
Show the code 
emm_df <- as.data.frame(emm)

# Plot the results
ggplot(emm_df, aes(x = vol_cat, y = emmean, ymin = lower.CL, ymax = upper.CL)) +
  geom_point() + # Points for the means
  geom_errorbar(width = 0.2) + # Error bars for confidence intervals
  facet_wrap(~ tx_group) + # Separate panels for each value of vol
  labs(
    title = "Estimated Marginal Means of igg by tx_group and vol",
    x = "Treatment Group (tx_group)",
    y = "Estimated Mean igg (with 95% CI)"
  )

Volume as a main effect

Show the code 
figure_scatter <- ggplot(data %>% filter(!is.na(tx_group)),
                         aes(x = vol, y = igg)) +
  # Background zones
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -5, ymax = 10),
            fill = "pink", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 10, ymax = 18),
            fill = "#FFF2CC", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 18, ymax = 25),
            fill = "#C3D69B", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 25, ymax = Inf),
            fill = "#A2D9CE", alpha = 0.5, inherit.aes = FALSE) +
  
  # Scatter points
  geom_point(size = 1, alpha = 0.8, colour = "black") +
  
  # Optionally add a smoothed line (remove if not desired)
  geom_smooth(method = "lm", se = T, linewidth = 0.5, linetype = "dashed",colour = "black") +
  
  # Labels
  labs(
    title = "Relationship between Volume Fed and Serum IgG",
    x = "Volume fed (L)",
    y = "Serum IgG (g/L)"
  ) +
  
  # Facet by treatment group if relevant
  facet_wrap(~ tx_group) +
  
  coord_cartesian(ylim = c(0, 60)) +
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.border = element_rect(
      color = "black", fill = NA, linewidth = 1
    ),
    text = element_text(family = "Times New Roman")
  )

figure_scatter

Show the code 
ggsave(
  filename = "figure scatter vol.tiff",    # File name
  plot = figure_scatter,              # ggplot object to save
  device = "tiff",             # Output format
  width = 6,                   # Width in inches
  height = 5,                  # Height in inches
  dpi = 200                    # Resolution in dots per inch
)

Univariable model for effect of volume on IgG.

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

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

Residuals:
     Min       1Q   Median       3Q      Max 
-24.0073  -5.0283   0.3549   5.1163  20.6597 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  18.6529     2.1394   8.719  < 2e-16 ***
vol           4.1780     0.6841   6.107 2.45e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.501 on 390 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.08729,   Adjusted R-squared:  0.08495 
F-statistic:  37.3 on 1 and 390 DF,  p-value: 2.453e-09
Show the code 
broom::tidy(lm_not_control_for_dose, conf.int=T)
term estimate std.error statistic p.value conf.low conf.high
(Intercept) 18.652909 2.1393531 8.718948 0 14.446801 22.85902
vol 4.177978 0.6841031 6.107235 0 2.832987 5.52297

Model for effect of volume on IgG after controlling for IgG dose

Show the code 
lm_control_for_dose <- lm(igg ~ tx_group + vol, data = data)

summary(lm_control_for_dose)

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

Residuals:
     Min       1Q   Median       3Q      Max 
-23.1171  -4.6460  -0.0524   4.7329  19.6052 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    28.068      2.549  11.009  < 2e-16 ***
tx_group250g    5.553      1.359   4.086 5.34e-05 ***
tx_group300g    8.388      1.674   5.011 8.24e-07 ***
tx_group350g   13.178      1.884   6.993 1.19e-11 ***
vol            -1.422      1.047  -1.358    0.175    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.058 on 387 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.1982,    Adjusted R-squared:  0.1899 
F-statistic: 23.92 on 4 and 387 DF,  p-value: < 2.2e-16
Show the code 
broom::tidy(lm_control_for_dose, conf.int=T)
term estimate std.error statistic p.value conf.low conf.high
(Intercept) 28.068299 2.549476 11.009436 0.0000000 23.055741 33.0808572
tx_group250g 5.553479 1.359141 4.086021 0.0000534 2.881255 8.2257038
tx_group300g 8.388434 1.673861 5.011427 0.0000008 5.097434 11.6794340
tx_group350g 13.178252 1.884444 6.993178 0.0000000 9.473223 16.8832809
vol -1.421567 1.047104 -1.357618 0.1753760 -3.480291 0.6371571

Colostrum concentration

Concentration as an effect modifier

Concentration as a main effect

Show the code 
figure_scatter <- ggplot(data %>% filter(!is.na(tx_group)),
                         aes(x = concentration, y = igg)) +
  # Background zones
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -5, ymax = 10),
            fill = "pink", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 10, ymax = 18),
            fill = "#FFF2CC", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 18, ymax = 25),
            fill = "#C3D69B", alpha = 0.5, inherit.aes = FALSE) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 25, ymax = Inf),
            fill = "#A2D9CE", alpha = 0.5, inherit.aes = FALSE) +
  
  # Scatter points
  geom_point(size = 1, alpha = 0.8, colour = "black") +
  
  # Optionally add a smoothed line (remove if not desired)
  geom_smooth(method = "lm", se = T, linewidth = 0.5, linetype = "dashed",colour = "black") +
  
  # Labels
  labs(
    title = "Relationship between colostrum colostrum and Serum IgG",
    x = "Concentration (g/L)",
    y = "Serum IgG (g/L)"
  ) +
  
  # Facet by treatment group if relevant
  facet_wrap(~ tx_group) +
  
  coord_cartesian(ylim = c(0, 60)) +
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.border = element_rect(
      color = "black", fill = NA, linewidth = 1
    ),
    text = element_text(family = "Times New Roman")
  )

figure_scatter

Show the code 
ggsave(
  filename = "figure scatter conc.tiff",    # File name
  plot = figure_scatter,              # ggplot object to save
  device = "tiff",             # Output format
  width = 6,                   # Width in inches
  height = 5,                  # Height in inches
  dpi = 200                    # Resolution in dots per inch
)

Univariable model for effect of concentration on IgG.

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

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

Residuals:
    Min      1Q  Median      3Q     Max 
-22.709  -5.684   0.743   4.972  21.750 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   21.19922    3.22751   6.568 1.63e-10 ***
concentration  0.10987    0.03413   3.219  0.00139 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.75 on 390 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.02588,   Adjusted R-squared:  0.02338 
F-statistic: 10.36 on 1 and 390 DF,  p-value: 0.001394
Show the code 
broom::tidy(lm_not_control_for_dose, conf.int=T)
term estimate std.error statistic p.value conf.low conf.high
(Intercept) 21.199223 3.2275128 6.568285 0.0000000 14.8537223 27.5447244
concentration 0.109867 0.0341308 3.219004 0.0013941 0.0427637 0.1769703

Controlling for IgG dose

Show the code 
lm_control_for_dose <- lm(igg ~ tx_group + concentration, data = data)

summary(lm_control_for_dose)

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

Residuals:
     Min       1Q   Median       3Q      Max 
-23.0381  -4.5991  -0.0886   4.7512  19.6476 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   21.24441    3.13803   6.770 4.79e-11 ***
tx_group250g   4.74616    1.23024   3.858 0.000134 ***
tx_group300g   6.81478    1.23606   5.513 6.45e-08 ***
tx_group350g  10.95843    1.28328   8.539 3.12e-16 ***
concentration  0.04021    0.03238   1.242 0.215048    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.061 on 387 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.1976,    Adjusted R-squared:  0.1893 
F-statistic: 23.82 on 4 and 387 DF,  p-value: < 2.2e-16
Show the code 
broom::tidy(lm_control_for_dose, conf.int=T)
term estimate std.error statistic p.value conf.low conf.high
(Intercept) 21.2444098 3.1380273 6.769989 0.0000000 15.0746942 27.4141254
tx_group250g 4.7461575 1.2302389 3.857915 0.0001339 2.3273691 7.1649458
tx_group300g 6.8147777 1.2360645 5.513287 0.0000001 4.3845356 9.2450199
tx_group350g 10.9584281 1.2832788 8.539398 0.0000000 8.4353573 13.4814989
concentration 0.0402088 0.0323785 1.241839 0.2150482 -0.0234509 0.1038686

Volume and Concentration in a single model

Show the code 
lm_vol_concentration <- lm(igg ~ vol*concentration, data = data)

summary(lm_vol_concentration)

Call:
lm(formula = igg ~ vol * concentration, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-23.7497  -4.7729  -0.0067   4.7928  20.1867 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)  
(Intercept)       21.57259   17.37787   1.241   0.2152  
vol               -4.43984    5.72600  -0.775   0.4386  
concentration     -0.09388    0.17920  -0.524   0.6007  
vol:concentration  0.11325    0.06023   1.880   0.0608 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.069 on 388 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.1936,    Adjusted R-squared:  0.1873 
F-statistic: 31.04 on 3 and 388 DF,  p-value: < 2.2e-16
Show the code 
broom::tidy(lm_vol_concentration, conf.int=T)
term estimate std.error statistic p.value conf.low conf.high
(Intercept) 21.5725941 17.3778657 1.2413834 0.2152141 -12.5939734 55.7391616
vol -4.4398443 5.7259976 -0.7753835 0.4385854 -15.6977104 6.8180218
concentration -0.0938767 0.1792024 -0.5238588 0.6006758 -0.4462059 0.2584525
vol:concentration 0.1132540 0.0602350 1.8802023 0.0608294 -0.0051739 0.2316818
Show the code 
check_collinearity(lm_vol_concentration)
Term VIF VIF_CI_low VIF_CI_high SE_factor Tolerance Tolerance_CI_low Tolerance_CI_high
vol 78.88442 65.02400 95.74546 8.881690 0.0126768 0.0104444 0.0153789
concentration 33.12883 27.36416 40.15397 5.755765 0.0301852 0.0249041 0.0365441
vol:concentration 68.65339 56.60318 83.31509 8.285734 0.0145659 0.0120026 0.0176669

Checking if this is equivalent to linear model where dose is on the continuous scale and is included instead of the vol*concentration interaction variable

Show the code 
lm_vol_concentration_dose <- lm(igg ~ vol + concentration + tx_administered, data = data)

summary(lm_vol_concentration_dose)

Call:
lm(formula = igg ~ vol + concentration + tx_administered, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-23.7497  -4.7729  -0.0067   4.7928  20.1867 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)  
(Intercept)     21.57259   17.37787   1.241   0.2152  
vol             -4.43984    5.72600  -0.775   0.4386  
concentration   -0.09388    0.17920  -0.524   0.6007  
tx_administered  0.11325    0.06023   1.880   0.0608 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.069 on 388 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.1936,    Adjusted R-squared:  0.1873 
F-statistic: 31.04 on 3 and 388 DF,  p-value: < 2.2e-16
Show the code 
broom::tidy(lm_vol_concentration_dose, conf.int=T)
term estimate std.error statistic p.value conf.low conf.high
(Intercept) 21.5725941 17.3778657 1.2413834 0.2152141 -12.5939734 55.7391616
vol -4.4398443 5.7259976 -0.7753835 0.4385854 -15.6977104 6.8180218
concentration -0.0938767 0.1792024 -0.5238588 0.6006758 -0.4462059 0.2584525
tx_administered 0.1132540 0.0602350 1.8802023 0.0608294 -0.0051739 0.2316818
Show the code 
check_collinearity(lm_vol_concentration_dose)
Term VIF VIF_CI_low VIF_CI_high SE_factor Tolerance Tolerance_CI_low Tolerance_CI_high
vol 78.88442 65.02400 95.74546 8.881690 0.0126768 0.0104444 0.0153789
concentration 33.12883 27.36416 40.15397 5.755765 0.0301852 0.0249041 0.0365441
tx_administered 68.65339 56.60318 83.31509 8.285734 0.0145659 0.0120026 0.0176669

Objective 2 - Practically feasible dose

Our objective is to determine the amount of IgG available to feed per calf from pooled and pasteurised batches of 40L, when there is a limiting effect of colostral IgG concentration and volume administered.

Show the code 
day
date born igg_kg concentration_limited_dose igg_per_calf
2024-10-24 14 8.798840 317.9076 628.4886
2024-10-25 22 8.254016 347.6477 375.1826
2024-10-26 15 7.864941 372.9109 524.3294
2024-10-27 28 7.254011 372.6910 259.0718
2024-10-28 18 6.837712 323.6434 379.8729
2024-10-29 17 7.940289 321.0714 467.0758
2024-10-30 20 6.121256 370.2976 306.0628
2024-10-31 24 8.923116 348.0520 371.7965
2024-11-01 26 8.435949 328.6460 324.4596
2024-11-02 15 6.238196 283.2936 415.8797
2024-11-03 21 4.299691 338.4722 204.7472
2024-11-04 24 8.903345 359.0180 370.9727
2024-11-05 19 2.424971 373.0724 127.6300
2024-11-06 32 7.826972 304.1152 244.5929
2024-11-07 38 5.596646 365.3151 147.2802
2024-11-08 35 12.133704 294.0847 346.6773

Number born per day

Show the code 
ggplot(data = day, aes(x = date, y = born)) +
  geom_col() +
  labs(
    title = "",
    x = "Date",
    y = "Number born"
  ) +
  scale_x_date(date_breaks = "1 day", date_labels = "%b %d") +  # Custom date labels
  theme_minimal() +
  theme(
    text = element_text(family = "Times New Roman"),  # Set Times New Roman font
    axis.text.x = element_text(angle = 90, hjust = 0),
    panel.grid.major.x = element_blank(),  # Remove major vertical grid lines
    panel.grid.minor.x = element_blank()   # Remove minor vertical grid lines
  )

  • From this histogram, it is clear that the number of calves born per day have increased since November 6th and that the demand for colostrum is higher. Therefore colostrum recommendations must be made to account for times of peak demand.

IgG available per calf (due to limitation of total IgG harvested on that day - mass limited dose)

Show the code 
ggplot(data = day, aes(x = date, y = igg_per_calf)) +
  geom_col() +
  labs(
    title = "Mass limited dose",
    x = "",
    y = "IgG (g) available per calf"
  ) +
  scale_x_date(date_breaks = "1 day", date_labels = "%b %d") +  # Custom date labels
  theme_minimal() +
  theme(
    text = element_text(family = "Times New Roman"),  # Set Times New Roman font
    axis.text.x = element_text(angle = 90, hjust = 0),
    panel.grid.major.x = element_blank(),  # Remove major vertical grid lines
    panel.grid.minor.x = element_blank()   # Remove minor vertical grid lines
  )

Show the code 
mass_day <- ggplot(data = day, 
       aes(y = igg_per_calf, x = 1)) +  # Add a dummy x-axis (1) for a single box
  geom_boxplot(outlier.shape = NA, width = 0.3) +  # Box plot with adjusted width and no outliers
  geom_jitter(width = 0.01, alpha = 0.6, size = 1.5) +  # Jitter with slight width, transparency, and size
  labs(
    y = "IgG (g) available per calf",
    x = NULL
  ) +
  scale_y_continuous(
    limits = c(100, 640),  # Set y-axis range
    breaks = seq(100, 640, by = 20)  # Primary y-axis gridlines every 20 units
  ) +
  theme_minimal() +
  theme(
    text = element_text(family = "Times New Roman"),  # Set Times New Roman font
    panel.grid.major.x = element_blank(),  # Remove major vertical grid lines
    panel.grid.minor.x = element_blank(),  # Remove minor vertical grid lines
    panel.grid.major.y = element_line(color = "grey80"),  # Keep primary y-axis gridlines
    panel.grid.minor.y = element_blank(),  # Remove secondary gridlines
    axis.text.x = element_blank(),  # Remove x-axis text
    axis.title.x = element_blank()  # Remove x-axis title
  )

ggsave(
  filename = "figure mass day.tiff",    # File name
  plot = mass_day,              # ggplot object to save
  device = "tiff",             # Output format
  width = 3,                   # Width in inches
  height = 5,                  # Height in inches
  dpi = 200                    # Resolution in dots per inch
)

mass_day

IgG available per calf each day (due to limitation of 4L per calf = concentration limited dose)

Show the code 
ggplot(data = day, aes(x = date, y = concentration_limited_dose)) +
  geom_col() +
  labs(
    title = "Concentration limited dose",
    x = "",
    y = "IgG (g) available per calf"
  ) +
  scale_x_date(date_breaks = "1 day", date_labels = "%b %d") +  # Custom date labels
  theme_minimal() +
  theme(
    text = element_text(family = "Times New Roman"),  # Set Times New Roman font
    axis.text.x = element_text(angle = 90, hjust = 0),
    panel.grid.major.x = element_blank(),  # Remove major vertical grid lines
    panel.grid.minor.x = element_blank()   # Remove minor vertical grid lines
  )

Concentration limited dose by batch

Show the code 
ggplot(data = batch %>% arrange(dose_at_4l) %>%
         mutate(batch = 1:nrow(batch)), 
       aes(x = batch, y = dose_at_4l)) +
  geom_col() +
  labs(
    x = "Batches sorted by concentration",
  ) +
  theme_minimal() +
  scale_y_continuous(
    name = "IgG (g) in 4L feed",
    sec.axis = sec_axis(~ . / 4, name = "IgG (g/L) in 1L")
  ) +
  theme(
    text = element_text(family = "Times New Roman"),  # Set Times New Roman font
    panel.grid.major.x = element_blank(),  # Remove major vertical grid lines
    panel.grid.minor.x = element_blank()   # Remove minor vertical grid lines
  )

Show the code 
concentration <- ggplot(data = batch, 
       aes(y = dose_at_4l, x = 1)) +  # Add a dummy x-axis (1) for a single box
  geom_boxplot(outlier.shape = NA, width = 0.3) +  # Box plot with adjusted width and no outliers
  geom_jitter(width = 0.01, alpha = 0.6, size = 1.5) +  # Jitter with slight width, transparency, and size
  labs(
    y = "IgG (g) in 4L (i.e., max dose)",
    x = NULL  # Remove x-axis title
  ) +
  scale_y_continuous(
    limits = c(260, 500),  # Set y-axis range
    breaks = seq(260, 500, by = 20),  # Primary y-axis gridlines every 20 units
    sec.axis = sec_axis(
      ~ . / 4,
      name = "IgG (g) in 1L",
      breaks = seq(260, 500, by = 20) / 4  # Align secondary axis breaks
    )
  ) +
  theme_minimal() +
  theme(
    text = element_text(family = "Times New Roman"),  # Set Times New Roman font
    panel.grid.major.x = element_blank(),  # Remove major vertical grid lines
    panel.grid.minor.x = element_blank(),  # Remove minor vertical grid lines
    axis.text.x = element_blank(),  # Remove x-axis text
    axis.title.x = element_blank()  # Remove x-axis title
  )

ggsave(
  filename = "figure concentration.tiff",    # File name
  plot = concentration,              # ggplot object to save
  device = "tiff",             # Output format
  width = 3,                   # Width in inches
  height = 5,                  # Height in inches
  dpi = 200                    # Resolution in dots per inch
)

concentration

Proportion of batches above certain thresholds

Show the code 
batch %>% summarise(
  dose_greater_equal_280 = paste0(round(100*sum(batch$dose_at_4l >= 270) / n(),1)," %"),
  dose_greater_equal_300 = paste0(round(100*sum(batch$dose_at_4l >= 300) / n(),1)," %"),
  dose_greater_equal_320 = paste0(round(100*sum(batch$dose_at_4l >= 320) / n(),1)," %"),
  dose_greater_equal_340 = paste0(round(100*sum(batch$dose_at_4l >= 340) / n(),1)," %"),
  dose_greater_equal_350 = paste0(round(100*sum(batch$dose_at_4l >= 350) / n(),1)," %")
) %>% t %>% as.table
                       A     
dose_greater_equal_280 100 % 
dose_greater_equal_300 93.2 %
dose_greater_equal_320 88.6 %
dose_greater_equal_340 72.7 %
dose_greater_equal_350 65.9 %
Show the code 
tibble(
  date_range = paste0(min(day$date),
                      " to ",
                      max(day$date)),
  number_of_days = nrow(day),
  calves_born = sum(day$born),
  igg_harvested_kg = round(sum(day$igg_kg),1),
  igg_available_per_calf_g = round(igg_harvested_kg/calves_born * 1000,0)

)
date_range number_of_days calves_born igg_harvested_kg igg_available_per_calf_g
2024-10-24 to 2024-11-08 16 368 117.9 320

Export data

Show the code 
#data %>% write_rds('igg_trial.rds')

Secondary analyses

Comparison of randomisation periods

Show the code 
weather_data_cleaned <- weather_data_2 %>% 
  clean_names %>%
    mutate(
      randomisation_period = case_when(
        day <= as.Date("2024-10-30") ~ "Days 1 - 7",
        T ~ "Days 8 - 18"
      )
    ) %>%
  filter(between(day,data$enrol_date %>% min %>% as.Date,
                 data$enrol_date %>% max %>% as.Date))

table1::table1(~ max_temp + rainfall | randomisation_period, data = weather_data_cleaned)
Days 1 - 7
(N=7)		 Days 8 - 18
(N=11)		 Overall
(N=18)
max_temp
Mean (SD) 24.0 (4.10) 26.2 (4.64) 25.3 (4.45)
Median [Min, Max] 23.3 [18.9, 30.6] 24.5 [21.3, 36.3] 24.2 [18.9, 36.3]
rainfall
Mean (SD) 0 (0) 0 (0) 0 (0)
Median [Min, Max] 0 [0, 0] 0 [0, 0] 0 [0, 0]
Show the code 
lm(max_temp ~ randomisation_period, data = weather_data_cleaned) %>% summary

Call:
lm(formula = max_temp ~ randomisation_period, data = weather_data_cleaned)

Residuals:
   Min     1Q Median     3Q    Max 
-5.071 -3.464 -1.021  1.836 10.136 

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)    
(Intercept)                       23.971      1.680   14.27 1.61e-10 ***
randomisation_periodDays 8 - 18    2.192      2.149    1.02    0.323    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.445 on 16 degrees of freedom
Multiple R-squared:  0.06107,   Adjusted R-squared:  0.002382 
F-statistic: 1.041 on 1 and 16 DF,  p-value: 0.3229
Show the code 
lm(rainfall ~ randomisation_period, data = weather_data_cleaned) %>% summary

Call:
lm(formula = rainfall ~ randomisation_period, data = weather_data_cleaned)

Residuals:
   Min     1Q Median     3Q    Max 
     0      0      0      0      0 

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)
(Intercept)                            0          0     NaN      NaN
randomisation_periodDays 8 - 18        0          0     NaN      NaN

Residual standard error: 0 on 16 degrees of freedom
Multiple R-squared:    NaN, Adjusted R-squared:    NaN 
F-statistic:   NaN on 1 and 16 DF,  p-value: NA
Show the code 
data <- data %>%
  mutate(
    randomisation_period = case_when(
      enrol_date <= as.Date("2024-10-30") ~ "Days 1 - 7",
      T ~ "Days 8 - 18"
    )
  )


table1::table1(~ concentration + vol + sex + breed | randomisation_period, data = data)
Days 1 - 7
(N=131)		 Days 8 - 18
(N=268)		 Overall
(N=399)
concentration
Mean (SD) 93.1 (11.0) 94.1 (11.7) 93.8 (11.5)
Median [Min, Max] 90.9 [76.9, 120] 92.6 [71.4, 121] 92.6 [71.4, 121]
vol
Mean (SD) 2.71 (0.519) 3.27 (0.471) 3.08 (0.552)
Median [Min, Max] 2.60 [1.80, 3.80] 3.40 [2.10, 3.90] 3.10 [1.80, 3.90]
sex
F 98 (74.8%) 220 (82.1%) 318 (79.7%)
M 33 (25.2%) 48 (17.9%) 81 (20.3%)
breed
Angus x Holstein 57 (43.5%) 69 (25.7%) 126 (31.6%)
Holstein 74 (56.5%) 199 (74.3%) 273 (68.4%)
Show the code 
lm(concentration ~ randomisation_period, data = data) %>% summary

Call:
lm(formula = concentration ~ randomisation_period, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-22.707  -7.929  -1.543   5.864  26.938 

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)    
(Intercept)                       93.062      1.005  92.637   <2e-16 ***
randomisation_periodDays 8 - 18    1.074      1.226   0.876    0.381    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 11.5 on 397 degrees of freedom
Multiple R-squared:  0.001931,  Adjusted R-squared:  -0.0005833 
F-statistic: 0.768 on 1 and 397 DF,  p-value: 0.3814

Sample size calculations

Sample size calculations were conducted prior to the start of the trial. The goal was to demonstrate differences in proportions of calves with serum IgG > 10, 15 and 25g/L cut-off points.

The plot below is from Osaka et al. 2014.

https://doi.org/10.3168/jds.2013-7571

The B line is probably most appropriate for what we can expect the relationship between oral IgG dose and serum IgG. y = 0.0804 * oral dose. The standard deviation is not reported in the paper, but back calculation from the SE estimates in table 1 (SD = SE * sqrt(n)) give SD values of 4, 7, and 6 for groups B, C, and D, respectively. I will calculate sample sizes for SD values of 5 and 6.

Show the code 
data <- tibble(
  igg_oral = c(200,250,300,350),
  igg_serum_mean = round(igg_oral*0.0804,0),
  igg_excellent_sd5 = 1 - pnorm(25, igg_serum_mean, 5),
  igg_excellent_sd6 = 1 - pnorm(25, igg_serum_mean, 6),
  igg_ftpi_10_sd5 = 1 - pnorm(10, igg_serum_mean, 5),
  igg_ftpi_10_sd6 = 1 - pnorm(10, igg_serum_mean, 6),
  igg_ftpi_15_sd5 = 1 - pnorm(15, igg_serum_mean, 5),
  igg_ftpi_15_sd6 = 1 - pnorm(15, igg_serum_mean, 6),
  
)

data
igg_oral igg_serum_mean igg_excellent_sd5 igg_excellent_sd6 igg_ftpi_10_sd5 igg_ftpi_10_sd6 igg_ftpi_15_sd5 igg_ftpi_15_sd6
200 16 0.0359303 0.0668072 0.8849303 0.8413447 0.5792597 0.5661838
250 20 0.1586553 0.2023284 0.9772499 0.9522096 0.8413447 0.7976716
300 24 0.4207403 0.4338162 0.9974449 0.9901847 0.9640697 0.9331928
350 28 0.7257469 0.6914625 0.9998409 0.9986501 0.9953388 0.9848699

Sample size calculation for 25g/L cut-off

Assuming SD = 5

250 vs 200

Show the code 
treat <- 250; control <- 200

treat <- data %>% filter(igg_oral == treat) %>% pull(igg_excellent_sd5)
control <- data %>% filter(igg_oral == control) %>% pull(igg_excellent_sd5)

epi.sssupb(treat = treat, control = control,
           delta = 0, power = 0.8, alpha = 0.05, n = NA) %>% as_tibble
n.total n.treat n.control delta power
140 70 70 0 0.8

300 vs 250

Show the code 
treat <- 300; control <- 250

treat <- data %>% filter(igg_oral == treat) %>% pull(igg_excellent_sd5)
control <- data %>% filter(igg_oral == control) %>% pull(igg_excellent_sd5)

epi.sssupb(treat = treat, control = control,
           delta = 0, power = 0.8, alpha = 0.05, n = NA) %>% as_tibble
n.total n.treat n.control delta power
68 34 34 0 0.8

350 vs 300

Show the code 
treat <- 350; control <- 300

treat <- data %>% filter(igg_oral == treat) %>% pull(igg_excellent_sd5)
control <- data %>% filter(igg_oral == control) %>% pull(igg_excellent_sd5)

epi.sssupb(treat = treat, control = control,
           delta = 0, power = 0.8, alpha = 0.05, n = NA) %>% as_tibble
n.total n.treat n.control delta power
60 30 30 0 0.8

Assuming SD = 6

250 vs 200

Show the code 
treat <- 250; control <- 200

treat <- data %>% filter(igg_oral == treat) %>% pull(igg_excellent_sd6)
control <- data %>% filter(igg_oral == control) %>% pull(igg_excellent_sd6)

epi.sssupb(treat = treat, control = control,
           delta = 0, power = 0.8, alpha = 0.05, n = NA) %>% as_tibble
n.total n.treat n.control delta power
152 76 76 0 0.8

300 vs 250

Show the code 
treat <- 300; control <- 250

treat <- data %>% filter(igg_oral == treat) %>% pull(igg_excellent_sd6)
control <- data %>% filter(igg_oral == control) %>% pull(igg_excellent_sd6)

epi.sssupb(treat = treat, control = control,
           delta = 0, power = 0.8, alpha = 0.05, n = NA) %>% as_tibble
n.total n.treat n.control delta power
94 47 47 0 0.8

350 vs 300

Show the code 
treat <- 350; control <- 300

treat <- data %>% filter(igg_oral == treat) %>% pull(igg_excellent_sd6)
control <- data %>% filter(igg_oral == control) %>% pull(igg_excellent_sd6)

epi.sssupb(treat = treat, control = control,
           delta = 0, power = 0.8, alpha = 0.05, n = NA) %>% as_tibble
n.total n.treat n.control delta power
86 43 43 0 0.8

Sample size calculation for 10g/L

Assuming SD = 5

250 vs 200

Show the code 
treat <- 250; control <- 200

treat <- data %>% filter(igg_oral == treat) %>% pull(igg_ftpi_10_sd5)
control <- data %>% filter(igg_oral == control) %>% pull(igg_ftpi_10_sd5)

epi.sssupb(treat = treat, control = control,
           delta = 0, power = 0.8, alpha = 0.05, n = NA) %>% as_tibble
n.total n.treat n.control delta power
180 90 90 0 0.8

300 vs 250

Show the code 
treat <- 300; control <- 250

treat <- data %>% filter(igg_oral == treat) %>% pull(igg_ftpi_10_sd5)
control <- data %>% filter(igg_oral == control) %>% pull(igg_ftpi_10_sd5)

epi.sssupb(treat = treat, control = control,
           delta = 0, power = 0.8, alpha = 0.05, n = NA) %>% as_tibble
n.total n.treat n.control delta power
752 376 376 0 0.8

350 vs 300

Show the code 
treat <- 350; control <- 300

treat <- data %>% filter(igg_oral == treat) %>% pull(igg_ftpi_10_sd5)
control <- data %>% filter(igg_oral == control) %>% pull(igg_ftpi_10_sd5)

epi.sssupb(treat = treat, control = control,
           delta = 0, power = 0.8, alpha = 0.05, n = NA) %>% as_tibble
n.total n.treat n.control delta power
5832 2916 2916 0 0.8

Assuming SD = 6

250 vs 200

Show the code 
treat <- 250; control <- 200

treat <- data %>% filter(igg_oral == treat) %>% pull(igg_ftpi_10_sd6)
control <- data %>% filter(igg_oral == control) %>% pull(igg_ftpi_10_sd6)

epi.sssupb(treat = treat, control = control,
           delta = 0, power = 0.8, alpha = 0.05, n = NA) %>% as_tibble
n.total n.treat n.control delta power
182 91 91 0 0.8

300 vs 250

Show the code 
treat <- 300; control <- 250

treat <- data %>% filter(igg_oral == treat) %>% pull(igg_ftpi_10_sd6)
control <- data %>% filter(igg_oral == control) %>% pull(igg_ftpi_10_sd6)

epi.sssupb(treat = treat, control = control,
           delta = 0, power = 0.8, alpha = 0.05, n = NA) %>% as_tibble
n.total n.treat n.control delta power
474 237 237 0 0.8

350 vs 300

Show the code 
treat <- 350; control <- 300

treat <- data %>% filter(igg_oral == treat) %>% pull(igg_ftpi_10_sd6)
control <- data %>% filter(igg_oral == control) %>% pull(igg_ftpi_10_sd6)

epi.sssupb(treat = treat, control = control,
           delta = 0, power = 0.8, alpha = 0.05, n = NA) %>% as_tibble
n.total n.treat n.control delta power
1910 955 955 0 0.8

Sample size calculation for 15g/L

Assuming SD = 5

250 vs 200

Show the code 
treat <- 250; control <- 200

treat <- data %>% filter(igg_oral == treat) %>% pull(igg_ftpi_15_sd5)
control <- data %>% filter(igg_oral == control) %>% pull(igg_ftpi_15_sd5)

epi.sssupb(treat = treat, control = control,
           delta = 0, power = 0.8, alpha = 0.05, n = NA) %>% as_tibble
n.total n.treat n.control delta power
68 34 34 0 0.8

300 vs 250

Show the code 
treat <- 300; control <- 250

treat <- data %>% filter(igg_oral == treat) %>% pull(igg_ftpi_15_sd5)
control <- data %>% filter(igg_oral == control) %>% pull(igg_ftpi_15_sd5)

epi.sssupb(treat = treat, control = control,
           delta = 0, power = 0.8, alpha = 0.05, n = NA) %>% as_tibble
n.total n.treat n.control delta power
140 70 70 0 0.8

350 vs 300

Show the code 
treat <- 350; control <- 300

treat <- data %>% filter(igg_oral == treat) %>% pull(igg_excellent_sd5)
control <- data %>% filter(igg_oral == control) %>% pull(igg_excellent_sd5)

epi.sssupb(treat = treat, control = control,
           delta = 0, power = 0.8, alpha = 0.05, n = NA) %>% as_tibble
n.total n.treat n.control delta power
60 30 30 0 0.8

Assuming SD = 6

250 vs 200

Show the code 
treat <- 250; control <- 200

treat <- data %>% filter(igg_oral == treat) %>% pull(igg_excellent_sd6)
control <- data %>% filter(igg_oral == control) %>% pull(igg_excellent_sd6)

epi.sssupb(treat = treat, control = control,
           delta = 0, power = 0.8, alpha = 0.05, n = NA) %>% as_tibble
n.total n.treat n.control delta power
152 76 76 0 0.8

300 vs 250

Show the code 
treat <- 300; control <- 250

treat <- data %>% filter(igg_oral == treat) %>% pull(igg_excellent_sd6)
control <- data %>% filter(igg_oral == control) %>% pull(igg_excellent_sd6)

epi.sssupb(treat = treat, control = control,
           delta = 0, power = 0.8, alpha = 0.05, n = NA) %>% as_tibble
n.total n.treat n.control delta power
94 47 47 0 0.8

350 vs 300

Show the code 
treat <- 350; control <- 300

treat <- data %>% filter(igg_oral == treat) %>% pull(igg_excellent_sd6)
control <- data %>% filter(igg_oral == control) %>% pull(igg_excellent_sd6)

epi.sssupb(treat = treat, control = control,
           delta = 0, power = 0.8, alpha = 0.05, n = NA) %>% as_tibble
n.total n.treat n.control delta power
86 43 43 0 0.8