Clean up the environment

rm(list = ls())

Load libraries for data reading, plotting and statistical analysis

library(readxl)       # For reading Excel files
library(ggpubr)       # For creating publication-ready plots
library(dplyr)        # For data manipulation
library(rstatix)      # For statistical tests and adding significance markers
library(knitr)        # For displaying tibbles in a nicely formatted table
library(kableExtra)   # For enhancing table formatting

Load data

# Load pig information, weight and intake data
PigInfo <- read_excel("Data/Data.xlsx", sheet = "PigInfo")
Weight <- read_excel("Data/Data.xlsx", sheet = "Weight")
Intake <- read_excel("Data/Data.xlsx", sheet = "Formula_intake")

Weight

Plot weight over time

Fig1 <- Weight %>% inner_join(PigInfo, by = "PigID") %>% 
                    filter(Timepoint_PD > 3) %>%
                    mutate(Diet = factor(Diet, levels = c("ALAC", "WPI", "SF"))) %>%
                    na.omit() %>%
                    ggline(y = "Weight_kg", x = "Timepoint_PD", color = "Diet", palette = c("turquoise3","purple", "orange"),
                           xlab = "Postnatal day", ylab = "Body weight (kg)", position = position_dodge(0.3), size = 0.7,
                           add = "mean_se") 
Fig1

Statistical analysis: overall difference via Tukey HSD

Weight.during.feeding <- Weight %>% inner_join(PigInfo, by = "PigID") %>% 
                           filter(Timepoint_PD > 4) %>%
                           mutate(Diet = factor(Diet, levels = c("ALAC", "WPI", "SF"))) %>%
                           mutate(Timepoint_PD = factor(Timepoint_PD))%>%
                           na.omit()

# Perform Tukey's HSD test
Weight.during.feeding %>% 
  tukey_hsd(Weight_kg ~ Diet + Timepoint_PD + SowID) %>% # Perform Tukey's HSD test
  filter(term == "Diet") %>%
  kable() %>% # Print the tibble in a nicely formatted table
  kable_styling(bootstrap_options = c("striped", "hover")) # Enhance table formatting with kableExtra
term group1 group2 null.value estimate conf.low conf.high p.adj p.adj.signif
Diet ALAC WPI 0 0.0445663 -0.0690802 0.1582127 0.626 ns
Diet ALAC SF 0 0.8462586 0.7080992 0.9844179 0.000 ****
Diet WPI SF 0 0.8016923 0.6656135 0.9377711 0.000 ****

Statistical analysis: difference at enrollment, adjusting for litter effect

Weight.at.baseline <- Weight %>% inner_join(PigInfo, by = "PigID") %>% 
                        filter(Timepoint_PD == 4) %>%
                        na.omit()

Weight.at.baseline %>% 
  tukey_hsd(Weight_kg ~ Diet + SowID) %>% # Perform Tukey's HSD test
  filter(term == "Diet") %>%
  kable() %>% # Print the tibble in a nicely formatted table
  kable_styling(bootstrap_options = c("striped", "hover")) # Enhance table formatting with kableExtra
term group1 group2 null.value estimate conf.low conf.high p.adj p.adj.signif
Diet ALAC SF 0 0.1976667 -0.0467093 0.4420426 0.129 ns
Diet ALAC WPI 0 0.1510000 -0.0516258 0.3536258 0.171 ns
Diet SF WPI 0 -0.0466667 -0.2832827 0.1899493 0.875 ns

Statistical analysis: cross-sectional difference during feeding period

Weight.during.feeding %>% 
  group_by(Timepoint_PD) %>% 
  tukey_hsd(Weight_kg ~ Diet + SowID) %>% # Perform Tukey's HSD test
  filter(term == "Diet") %>% 
  filter(p.adj.signif != "ns") %>% #only show pairs that are not significantly different
  kable() %>% # Print the tibble in a nicely formatted table
  kable_styling(bootstrap_options = c("striped", "hover")) # Enhance table formatting with kableExtra
Timepoint_PD term group1 group2 null.value estimate conf.low conf.high p.adj p.adj.signif
6 Diet ALAC SF 0 0.4157576 0.1199856 0.7115296 4.92e-03 **
7 Diet ALAC SF 0 0.6793939 0.4080429 0.9507449 5.30e-06 ****
7 Diet WPI SF 0 0.5425000 0.2751692 0.8098308 1.01e-04 ***
8 Diet ALAC SF 0 0.5919697 0.2736257 0.9103137 2.93e-04 ***
8 Diet WPI SF 0 0.5116667 0.1980391 0.8252942 1.23e-03 **
9 Diet ALAC SF 0 0.7075758 0.3403773 1.0747742 1.92e-04 ***
9 Diet WPI SF 0 0.5166667 0.1549085 0.8784248 4.29e-03 **
10 Diet ALAC SF 0 0.7507576 0.3456542 1.1558609 3.05e-04 ***
10 Diet WPI SF 0 0.6708333 0.2717318 1.0699348 9.00e-04 ***
11 Diet ALAC SF 0 0.8015152 0.3179058 1.2851245 1.04e-03 **
11 Diet WPI SF 0 0.8333333 0.3568890 1.3097777 5.87e-04 ***
12 Diet ALAC SF 0 0.8992424 0.3925593 1.4059256 4.99e-04 ***
12 Diet WPI SF 0 1.0291667 0.5299903 1.5283430 8.21e-05 ****
13 Diet ALAC SF 0 1.1287879 0.5896587 1.6679171 6.72e-05 ****
13 Diet WPI SF 0 1.2166667 0.6855250 1.7478083 1.97e-05 ****
14 Diet ALAC SF 0 1.3492424 0.7342385 1.9642463 3.59e-05 ****
14 Diet WPI SF 0 1.4583333 0.8524411 2.0642255 9.70e-06 ****
15 Diet ALAC SF 0 1.9466667 1.3020551 2.5912782 9.00e-07 ****
15 Diet WPI SF 0 1.8100000 1.1770040 2.4429960 2.00e-06 ****

Intake

# Plot formula intake normalized by weight

Fig2 <- inner_join(Weight, Intake, by = join_by(Timepoint_PD, PigID)) %>%
            mutate(Weighted_intake = Formula_intake_g/Weight_kg) %>% # Normalize intake by weight
            select(PigID, Timepoint_PD, Weighted_intake) %>%
            inner_join(PigInfo, by = "PigID") %>% 
            mutate(Diet = factor(Diet, levels = c("ALAC", "WPI"))) %>%
            na.omit() %>%
            filter(Timepoint_PD > 6) %>%
            ggline(y = "Weighted_intake", x = "Timepoint_PD", color = "Diet", palette = c("turquoise3","purple"), size = 0.7,
                   xlab = "Postnatal day", ylab = "Formula intake (g per kg body weight)", position = position_dodge(0.2),
                   add = "mean_se") 
Fig2

Statistical analysis: cross-sectional difference, adjusting for litter effect

inner_join(Weight, Intake, by = join_by(Timepoint_PD, PigID)) %>%
              mutate(Weighted_intake = Formula_intake_g/Weight_kg) %>%
              select(PigID, Timepoint_PD, Weighted_intake) %>%
              inner_join(PigInfo, by = "PigID") %>% 
              na.omit() %>%
              filter(Timepoint_PD > 6) %>%
              mutate(Timepoint_PD = factor(Timepoint_PD)) %>% 
              group_by(Timepoint_PD) %>% 
              anova_test(Weighted_intake ~ Diet + SowID, type = 1) %>%
              as_tibble() %>%
              filter(Effect == "Diet") %>%
              filter(`p<.05` == "*") %>% # Only displace significance
              kable() %>% # Print the tibble in a nicely formatted table
              kable_styling(bootstrap_options = c("striped", "hover")) # Enhance table formatting with kableExtra
Timepoint_PD Effect DFn DFd F p p<.05 ges
12 Diet 1 19 9.753 0.006000
0.339
13 Diet 1 19 22.368 0.000146
0.541

Final plot

Fig.all <- ggarrange(Fig1, Fig2, common.legend = TRUE)

Fig.all

# Save the figure in pdf:
ggsave(plot=Fig.all, height=4, width=10, dpi=300, filename="SI Figure 1/SI.Fig1.pdf", useDingbats=FALSE)
sessionInfo()
## R version 4.2.2 (2022-10-31)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur ... 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] kableExtra_1.3.4 knitr_1.45       rstatix_0.7.2    dplyr_1.1.3     
## [5] ggpubr_0.6.0     ggplot2_3.5.1    readxl_1.4.3    
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.0  xfun_0.40         bslib_0.5.1       purrr_1.0.2      
##  [5] carData_3.0-5     colorspace_2.1-0  vctrs_0.6.5       generics_0.1.3   
##  [9] viridisLite_0.4.2 htmltools_0.5.6.1 yaml_2.3.7        utf8_1.2.4       
## [13] rlang_1.1.2       jquerylib_0.1.4   pillar_1.9.0      glue_1.6.2       
## [17] withr_3.0.1       lifecycle_1.0.4   stringr_1.5.1     munsell_0.5.1    
## [21] ggsignif_0.6.4    gtable_0.3.5      cellranger_1.1.0  ragg_1.2.6       
## [25] rvest_1.0.3       evaluate_1.0.1    labeling_0.4.3    fastmap_1.1.1    
## [29] fansi_1.0.6       highr_0.10        broom_1.0.5       scales_1.3.0     
## [33] backports_1.4.1   cachem_1.0.8      webshot_0.5.5     jsonlite_1.8.8   
## [37] abind_1.4-5       farver_2.1.1      systemfonts_1.0.5 textshaping_0.3.7
## [41] gridExtra_2.3     digest_0.6.33     stringi_1.7.12    cowplot_1.1.1    
## [45] grid_4.2.2        cli_3.6.2         tools_4.2.2       magrittr_2.0.3   
## [49] sass_0.4.7        tibble_3.2.1      tidyr_1.3.0       car_3.1-2        
## [53] pkgconfig_2.0.3   xml2_1.3.5        rmarkdown_2.28    svglite_2.1.2    
## [57] httr_1.4.7        rstudioapi_0.15.0 R6_2.5.1          compiler_4.2.2