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