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

Authors and affiliations

Dini Hapukotuwa¹, John House¹, Katie Denholm², and Sam Rowe¹

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

²School 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

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"))

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

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-10-28

Check for duplicated IDs

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

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)

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))

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.

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

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.

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

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)

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

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

check_model(lm)

check_outliers(lm)

OK: No outliers detected.
- Based on the following method and threshold: cook (0.8).
- For variable: (Whole model)

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

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

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

check_model(lm)

check_outliers(lm)

OK: No outliers detected.
- Based on the following method and threshold: cook (0.7).
- For variable: (Whole model)

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

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.

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

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

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

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

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_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_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_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_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_interaction_model <- lm(igg ~ vol*concentration, data = data)
# summary(lm_interaction_model)

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")

tx	estimate_univariable	estimate_maximal	estimate_adjust_vol_conc	estimate_adjust_vol	estimate_adjust_conc
200g	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)
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)
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)

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

Mixed model vs fixed effect model

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

icc(lmm)

ICC_adjusted	ICC_unadjusted	optional
0.1590343	0.1338655	FALSE

lm_fixed_effect <- 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_fixed_effect = estimate) %>%
  select(tx,estimate_fixed_effect)

lmm <- lmer(igg ~ tx_group + (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	4.8 95% CI: (2.3, 7.2)	3.9 95% CI: (1.5, 6.3)
300g	6.9 95% CI: (4.4, 9.3)	6.2 95% CI: (3.8, 8.6)
350g	11.3 95% CI: (8.8, 13.7)	9.9 95% CI: (7.2, 12.6)

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

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

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
)

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

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
)

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

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

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

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.

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

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

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

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

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

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

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

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

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 

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

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

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.

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

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

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

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.

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

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

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

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

#pairs(emm)

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)"
  )

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

#pairs(emm)

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

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

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.

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

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

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

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

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

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.

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

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

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

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

Objective 2 - Effect of oral IgG dose on apparent efficiency of absorption

Weight data used for assumptions

Holsteins were slightly heavier, but this difference was not statistically significant. Decision is to use average weight across all calves.

weights %>% group_by(breed) %>%
  summarise(
    mean_weight = mean(weight),
    sd = sd(weight)
  )

breed	mean_weight	sd
A	42.08333	6.694072
H	38.85714	5.689406

lm(weight ~ breed, data = weights) %>% summary

Call:
lm(formula = weight ~ breed, data = weights)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.0833  -4.6637   0.1429   3.8929  13.1429 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   42.083      1.708  24.641   <2e-16 ***
breedH        -3.226      1.937  -1.666    0.102    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.916 on 52 degrees of freedom
Multiple R-squared:  0.05067,   Adjusted R-squared:  0.03241 
F-statistic: 2.775 on 1 and 52 DF,  p-value: 0.1017

weights %>%
    summarise(
    mean_weight = mean(weight),
    sd = sd(weight)
  )

mean_weight	sd
39.57407	6.014454

Plot

figure <- ggplot(data, 
                 aes(x = tx_group, y = aea)) +  # Removed fill from aes()
  
  # 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 = "Apparent efficiency of absorption",
    fill = "Treatment"
  ) +
  
  # Custom theme
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.border = element_rect(
      color = "black", fill = NA, linewidth = 1
    ),
    text = element_text(family = "Times New Roman")
  )

figure

Statistical models

Linear regression

lm <- lm(aea ~ tx_group + breed, data = data)
summary(lm)

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

Residuals:
      Min        1Q    Median        3Q       Max 
-0.239142 -0.045364  0.001324  0.046861  0.213060 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)    0.341049   0.011431  29.834  < 2e-16 ***
tx_group250g  -0.017667   0.012498  -1.414    0.158    
tx_group300g  -0.052887   0.012519  -4.225 2.99e-05 ***
tx_group350g  -0.060628   0.012836  -4.723 3.25e-06 ***
breedHolstein  0.008027   0.007919   1.014    0.311    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.07126 on 387 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.09045,   Adjusted R-squared:  0.08105 
F-statistic: 9.621 on 4 and 387 DF,  p-value: 1.997e-07

broom::tidy(lm, conf.int=T) %>% 
  select(term,estimate,conf.low,conf.high)

term	estimate	conf.low	conf.high
(Intercept)	0.3410488	0.3185733	0.3635244
tx_group250g	-0.0176667	-0.0422395	0.0069061
tx_group300g	-0.0528872	-0.0775008	-0.0282737
tx_group350g	-0.0606284	-0.0858649	-0.0353919
breedHolstein	0.0080271	-0.0075425	0.0235967

Estimated marginal means

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      0.345 0.01060 387    0.324    0.366  32.474  <.0001
 250g      0.327 0.00665 387    0.314    0.340  49.248  <.0001
 300g      0.292 0.00667 387    0.279    0.305  43.783  <.0001
 350g      0.284 0.00728 387    0.270    0.299  39.077  <.0001

Results are averaged over the levels of: breed 
Confidence level used: 0.95 

$contrasts
 contrast    estimate      SE  df lower.CL upper.CL t.ratio p.value
 200g - 250g  0.01767 0.01250 387  -0.0146   0.0499   1.414  0.4916
 200g - 300g  0.05289 0.01250 387   0.0206   0.0852   4.225  0.0002
 200g - 350g  0.06063 0.01280 387   0.0275   0.0937   4.723  <.0001
 250g - 300g  0.03522 0.00917 387   0.0116   0.0589   3.841  0.0008
 250g - 350g  0.04296 0.00948 387   0.0185   0.0674   4.533  <.0001
 300g - 350g  0.00774 0.00957 387  -0.0170   0.0324   0.809  0.8504

Results are averaged over the levels of: breed 
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 

Model diagnostics

check_model(lm)

check_outliers(lm)

OK: No outliers detected.
- Based on the following method and threshold: cook (0.5).
- For variable: (Whole model)

Table out

lm_all_breed <- lm(aea ~ tx_group + breed, 
                  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,3),
                        " (",
                        round(conf.low,3),
                        " to ",
                        round(conf.high,3),
                        ")"
  )),
  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,3),
                        " (",
                        round(lower.CL,3),
                        " to ",
                        round(upper.CL,3),
                        ")"
  )) %>%
  select(term,emmean)


lm_holstein <- lm(aea ~ tx_group, 
                  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,3),
                        " (",
                        round(conf.low,3),
                        " to ",
                        round(conf.high,3),
                        ")"
  )),
  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,3),
                        " (",
                        round(lower.CL,3),
                        " to ",
                        round(upper.CL,3),
                        ")"
  )) %>%
  select(term,emmean)


lm_angus <- lm(aea ~ tx_group,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,3),
                        " (",
                        round(conf.low,3),
                        " to ",
                        round(conf.high,3),
                        ")"
  )),
  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,3),
                        " (",
                        round(lower.CL,3),
                        " to ",
                        round(upper.CL,3),
                        ")"
  )) %>%
  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	0.345 (0.324 to 0.366)	0.018 (-0.007 to 0.042)	0.344 (0.316 to 0.373)	0.016 (-0.016 to 0.048)	0.346 (0.316 to 0.377)	0.015 (-0.023 to 0.054)
250g	0.327 (0.314 to 0.34)	Ref	0.328 (0.313 to 0.343)	Ref	0.331 (0.307 to 0.355)	Ref
300g	0.292 (0.279 to 0.305)	-0.035 (-0.053 to -0.017)	0.301 (0.285 to 0.316)	-0.028 (-0.05 to -0.006)	0.279 (0.257 to 0.302)	-0.052 (-0.085 to -0.019)
350g	0.284 (0.27 to 0.299)	-0.043 (-0.062 to -0.024)	0.289 (0.273 to 0.304)	-0.04 (-0.061 to -0.018)	0.279 (0.25 to 0.308)	-0.052 (-0.089 to -0.014)

Predicting AEA slope and point of diminishing returns

Need to identify if the relationship between IgG dose and AEA is linear or quadratic.

Scatter plot

ggplot(data, aes(x = tx_administered, y = aea)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", formula = y ~ x, color = "blue", se = FALSE, linetype = "dashed") +  # Linear
  geom_smooth(method = "lm", formula = y ~ poly(x, 2, raw = TRUE), color = "red", se = FALSE) +  # Quadratic
  labs(
    title = "Linear vs Quadratic Fit",
    x = "Oral IgG Dose (tx_administered)",
    y = "Apparent Efficiency of Absorption (AEA)"
  ) +
  theme_minimal()

Linear model

lm <- lm(aea ~ tx_administered + breed, data = data)
summary(lm)

Call:
lm(formula = aea ~ tx_administered + breed, data = data)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.234220 -0.044388 -0.002447  0.045777  0.217056 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)      4.282e-01  2.153e-02  19.889  < 2e-16 ***
tx_administered -4.388e-04  7.385e-05  -5.942 6.27e-09 ***
breedHolstein    8.953e-03  7.888e-03   1.135    0.257    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.07134 on 389 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.08366,   Adjusted R-squared:  0.07895 
F-statistic: 17.76 on 2 and 389 DF,  p-value: 4.167e-08

Quadratic model

lm_quad <- lm(aea ~ tx_administered + I(tx_administered^2) + breed, data = data)
summary(lm_quad)

Call:
lm(formula = aea ~ tx_administered + I(tx_administered^2) + breed, 
    data = data)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.232312 -0.044640 -0.002283  0.045461  0.218769 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           5.283e-01  1.222e-01   4.324 1.95e-05 ***
tx_administered      -1.171e-03  8.824e-04  -1.327    0.185    
I(tx_administered^2)  1.295e-06  1.556e-06   0.832    0.406    
breedHolstein         9.149e-03  7.895e-03   1.159    0.247    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.07137 on 388 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.0853,    Adjusted R-squared:  0.07822 
F-statistic: 12.06 on 3 and 388 DF,  p-value: 1.454e-07

Compare AIC

AIC(lm,lm_quad)

	df	AIC
lm	4	-952.5559
lm_quad	5	-951.2554

Compare using LRT

anova(lm,lm_quad)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
389	1.979783	NA	NA	NA	NA
388	1.976253	1	0.0035299	0.693039	0.4056447

No evidence that quadratic model is superior. Decision is to use linear model.

lm <- lm(aea ~ tx_administered + breed, data = data)
summary(lm)

Call:
lm(formula = aea ~ tx_administered + breed, data = data)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.234220 -0.044388 -0.002447  0.045777  0.217056 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)      4.282e-01  2.153e-02  19.889  < 2e-16 ***
tx_administered -4.388e-04  7.385e-05  -5.942 6.27e-09 ***
breedHolstein    8.953e-03  7.888e-03   1.135    0.257    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.07134 on 389 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.08366,   Adjusted R-squared:  0.07895 
F-statistic: 17.76 on 2 and 389 DF,  p-value: 4.167e-08

predictions <- emmeans(
  lm,
  specs = ~ tx_administered,
  at = list(tx_administered = c(seq(from = 200, to = 500, by = 1)))
) %>% 
  as_tibble %>%
  rename(aea = emmean,
         dose = tx_administered) %>%
  select(dose,aea) %>%
  mutate(
    predicted_serum_igg = (aea * dose) / (39 * 0.07)
  )

#predictions

# Models
mod_aea <- lm(aea * 100 ~ dose, data = predictions)
mod_igg <- lm(predicted_serum_igg ~ aea, data = predictions)

# Labels
coef_aea <- coef(mod_aea)
coef_igg <- coef(mod_igg)

label_aea <- sprintf("AEA = %.3f*Dose + %.1f", coef_aea[2], coef_aea[1])
#label_igg <- paste0("Serum IgG = (AEA*Dose) / (39 * 0.07)")

plot_data <- predictions %>%
  mutate(
    observed = if_else(dose <= 350, "Observed", "Extrapolated"),
    `Apparent efficiency of absorption` = aea * 100,  # convert to percentage
    `Predicted Serum IgG` = predicted_serum_igg         # no scaling now
  ) %>%
  pivot_longer(cols = c(`Apparent efficiency of absorption`, `Predicted Serum IgG`), 
               names_to = "Variable", values_to = "Value")


plot <- ggplot(plot_data, aes(x = dose, y = Value, color = Variable, linetype = observed)) +
  geom_line(size = 1) +
  scale_color_manual(
    values = c("Apparent efficiency of absorption" = "blue", "Predicted Serum IgG" = "red")
  ) +
  scale_linetype_manual(
    values = c("Observed" = "solid", "Extrapolated" = "dashed")
  ) +
  scale_y_continuous(
    limits = c(20, 40),
    name = "Apparent efficiency of absorption (%)",
    sec.axis = sec_axis(~ ., name = "Predicted Serum IgG (g/L)")
  ) +
  labs(
    x = "Oral IgG dose administered (g)",
    color = NULL
  ) +
  guides(
    linetype = "none"
  ) +
    annotate(
  "text",
  x = 370, y = 27.5,
  label = "AEA==42.6-0.043%*%Dose",
  parse = TRUE,
  color = "blue",
  hjust = 0,
  size = 4,
  family = "Times New Roman"
) +
  annotate(
  "text",
  x = 205, y = 37,
  label = "Predicted~Serum~IgG==frac(AEA%*%Dose, 39.6%*%0.07)",
  parse = TRUE,
  color = "red",
  hjust = 0,
  size = 4,
  family = "Times New Roman"
) +
  theme_minimal(base_size = 14) +
  theme(
    text = element_text(family = "Times New Roman"),
    legend.position = c(0.35, 0.15),
    legend.direction = "vertical",
    legend.background = element_rect(fill = "white", color = "black"),
    panel.border = element_rect(color = "black", fill = NA),
    axis.title = element_text(face = "bold")
  )


plot

ggsave(
  filename = "figure 5.tiff",    # File name
  plot = plot,              # 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
)

plot_data_blue <- plot_data %>%
  filter(Variable == "Apparent efficiency of absorption")

plot_blue <- ggplot(plot_data_blue, aes(x = dose, y = Value, linetype = observed)) +
  geom_line(size = 1, color = "blue") +  # color fixed to blue
  scale_linetype_manual(
    values = c("Observed" = "solid", "Extrapolated" = "dashed")
  ) +
  scale_y_continuous(
    limits = c(20, 40),
    name = "Apparent efficiency of absorption (%)",
    sec.axis = sec_axis(~ ., name = "Predicted Serum IgG (g/L)")  # keep the axis for flick match
  ) +
  labs(
    x = "Oral IgG dose administered (g)",
    color = NULL
  ) +
  guides(
    linetype = "none"
  ) +
  annotate(
    "text",
    x = 370, y = 27.5,
    label = "AEA==42.6-0.043%*%Dose",
    parse = TRUE,
    color = "blue",
    hjust = 0,
    size = 4,
    family = "Times New Roman"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    text = element_text(family = "Times New Roman"),
    legend.position = c(0.35, 0.15),
    legend.direction = "vertical",
    legend.background = element_rect(fill = "white", color = "black"),
    panel.border = element_rect(color = "black", fill = NA),
    axis.title = element_text(face = "bold")
  )

ggsave(
  filename = "figure_5a_blue_only.tiff",
  plot = plot_blue,
  device = "tiff",
  width = 6,
  height = 5,
  dpi = 200
)

# scaling_factor <- max(predictions$predicted_serum_igg) / max(predictions$aea)
# 
# plot_data <- predictions %>%
#   mutate(
#     observed = if_else(dose <= 350, "Observed", "Extrapolated"),
#     `Apparent efficiency of absorption` = aea,
#     `Predicted Serum IgG` = predicted_serum_igg / scaling_factor
#   ) %>%
#   pivot_longer(cols = c(`Apparent efficiency of absorption`, `Predicted Serum IgG`), 
#                names_to = "Variable", values_to = "Value")
# 
# 
# plot <- ggplot(plot_data, aes(x = dose, y = Value, color = Variable, linetype = observed)) +
#   geom_line(size = 1) +
#   scale_color_manual(
#     values = c("Apparent efficiency of absorption" = "blue", "Predicted Serum IgG" = "red")
#   ) +
#   scale_linetype_manual(
#     values = c("Observed" = "solid", "Extrapolated" = "dashed")
#   ) +
#   scale_y_continuous(
#     name = "Apparent efficiency of absorption\n(absorbed IgG / oral IgG dose)",
#     sec.axis = sec_axis(~ . * scaling_factor, name = "Predicted Serum IgG (g/L)")
#   ) +
#   labs(
#     x = "Oral IgG dose administered (g)",
#     color = NULL
#   ) +
#   guides(
#     linetype = "none"    # <- removes the second legend
#   ) +
#   theme_minimal(base_size = 14) +
#   theme(
#     text = element_text(family = "Times New Roman"),
#     legend.position = c(0.35, 0.15),
#     legend.direction = "vertical",
#     legend.background = element_rect(fill = "white", color = "black"),
#     panel.border = element_rect(color = "black", fill = NA),
#     axis.title = element_text(face = "bold")
#   )
# 
# ggsave(
#   filename = "figure 5.tiff",    # File name
#   plot = plot,              # ggplot object to save
#   device = "tiff",             # Output format
#   width = 7,                   # Width in inches
#   height = 6,                  # Height in inches
#   dpi = 200                    # Resolution in dots per inch
# )

Other factors impacting AEA

Volume

lm <- lm(aea ~ vol + tx_group, data = data)
summary(lm)

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

Residuals:
      Min        1Q    Median        3Q       Max 
-0.240578 -0.043929 -0.000324  0.046024  0.217035 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.374341   0.025723  14.553   <2e-16 ***
vol          -0.013081   0.010565  -1.238   0.2164    
tx_group250g -0.008972   0.013713  -0.654   0.5133    
tx_group300g -0.037659   0.016889  -2.230   0.0263 *  
tx_group350g -0.041184   0.019013  -2.166   0.0309 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.07121 on 387 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.09163,   Adjusted R-squared:  0.08224 
F-statistic:  9.76 on 4 and 387 DF,  p-value: 1.572e-07

Breed

lm <- lm(aea ~ breed, data = data)
summary(lm)

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

Residuals:
      Min        1Q    Median        3Q       Max 
-0.228669 -0.049367 -0.006369  0.047162  0.234575 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.306292   0.006821  44.902   <2e-16 ***
breedHolstein 0.003602   0.008174   0.441     0.66    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.07441 on 390 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.0004977, Adjusted R-squared:  -0.002065 
F-statistic: 0.1942 on 1 and 390 DF,  p-value: 0.6597

Time to feeding

lm <- lm(aea ~ time_to_feed, data = data)
summary(lm)

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

Residuals:
      Min        1Q    Median        3Q       Max 
-0.218394 -0.048627 -0.006357  0.045703  0.229665 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.3233242  0.0084406   38.31   <2e-16 ***
time_to_feed -0.0002434  0.0001268   -1.92   0.0556 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.07408 on 390 degrees of freedom
  (7 observations deleted due to missingness)
Multiple R-squared:  0.009359,  Adjusted R-squared:  0.006819 
F-statistic: 3.685 on 1 and 390 DF,  p-value: 0.05565

Objective 3 - 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.

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

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 6^th 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)

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
  )

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)

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

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
  )

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

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 %

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

#data %>% write_rds('igg_trial.rds')

Secondary analyses

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.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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