Clean up the environment

rm(list = ls())

Load necessary libraries for data reading, plotting, and statistical analysis

library(readxl)       # For reading Excel files
library(dplyr)        # For data manipulation
library(rstatix)      # For statistical tests
library(tidyr)        # For reshaping data

Load data

# Load Organ weight data
Organ <- read_excel("Data/Data.xlsx", sheet = "Necropsy_organ_wt")
PigInfo <- read_excel("Data/Data.xlsx", sheet = "PigInfo")

Data organization

Organ_weight <- inner_join(Organ, PigInfo, by = "PigID") %>%
                  mutate(Brain_perc = Brain_g/Nec_weight_kg/10) %>%
                  mutate(Liver_perc = Liver_g/Nec_weight_kg/10) %>%
                  mutate(Left_kidney_perc = Left_kidney_g/Nec_weight_kg/10) %>%
                  mutate(Right_kidney_perc = Right_kidney_g/Nec_weight_kg/10) %>%
                  mutate(Diet = factor(Diet, levels = c("ALAC", "WPI", "SF"))) %>%
                  select(Diet, SowID, Brain_g, Liver_g, Left_kidney_g, Right_kidney_g, Brain_perc, Liver_perc, Left_kidney_perc, Right_kidney_perc) %>%
                  pivot_longer(cols = c("Brain_g", "Liver_g", "Left_kidney_g", "Right_kidney_g", "Brain_perc", "Liver_perc", "Left_kidney_perc", "Right_kidney_perc"), names_to = "variable", values_to = "value") %>%
                  mutate(variable = factor(variable, levels = c("Brain_g", "Brain_perc", "Liver_g", "Liver_perc", "Left_kidney_g", "Left_kidney_perc", "Right_kidney_g", "Right_kidney_perc")))
# Compute sample sizes for each Diet
sample_sizes <- Organ_weight %>%
  group_by(Diet, variable) %>%
  summarise(n = n() / length(unique(variable)), .groups = 'drop')
Summary_stats <- Organ_weight %>%
                  group_by(Diet, variable) %>%
                  summarise(
                    Mean_SE = sprintf("%.2f ± %.2f", mean(value, na.rm = TRUE), sd(value, na.rm = TRUE) / sqrt(n())),
                    .groups = 'drop'
                  )
Summary_stats_wide <- Summary_stats %>% pivot_wider(names_from = Diet, values_from = Mean_SE)
# Modify column names to include sample sizes
names(Summary_stats_wide)[-1] <- paste(names(Summary_stats_wide)[-1], 
                               "(n =", 
                               sample_sizes$n[match(names(Summary_stats_wide)[-1], sample_sizes$Diet)],
                               ")", sep = "")

Summary_stats_wide
## # A tibble: 8 × 4
##   variable          `ALAC(n =11)` `WPI(n =12)` `SF(n =6)`   
##   <fct>             <chr>         <chr>        <chr>        
## 1 Brain_g           40.18 ± 0.80  39.84 ± 1.07 43.63 ± 1.77 
## 2 Brain_perc        1.46 ± 0.09   1.47 ± 0.06  0.92 ± 0.06  
## 3 Liver_g           83.58 ± 5.31  75.08 ± 5.09 132.77 ± 5.90
## 4 Liver_perc        2.96 ± 0.11   2.71 ± 0.08  2.76 ± 0.08  
## 5 Left_kidney_g     12.23 ± 1.05  13.78 ± 1.89 16.37 ± 0.88 
## 6 Left_kidney_perc  0.43 ± 0.02   0.49 ± 0.05  0.34 ± 0.02  
## 7 Right_kidney_g    11.88 ± 0.99  13.12 ± 1.57 15.92 ± 0.91 
## 8 Right_kidney_perc 0.42 ± 0.02   0.47 ± 0.04  0.33 ± 0.02
# Conduct ANOVA and Tukey HSD for each value variable
TukeyHSD_results <- Organ_weight %>%
                        group_by(variable) %>%
                        tukey_hsd(value ~ Diet + SowID) %>%
                        filter(term == "Diet") %>%
                        add_significance(p.col = "p.adj",
                          cutpoints = c(0, 0.001, 0.01, 0.05, 1),
                          symbols = c("***", "**", "*", "ns")) %>%
                        select(variable, group1, group2, p.adj.signif) %>%
                        pivot_wider(names_from = c("group1", "group2"), names_sep = " vs ", values_from = "p.adj.signif")
left_join(Summary_stats_wide, TukeyHSD_results, by = "variable")
## # A tibble: 8 × 7
##   variable      `ALAC(n =11)` `WPI(n =12)` `SF(n =6)` `ALAC vs WPI` `ALAC vs SF`
##   <fct>         <chr>         <chr>        <chr>      <chr>         <chr>       
## 1 Brain_g       40.18 ± 0.80  39.84 ± 1.07 43.63 ± 1… ns            ns          
## 2 Brain_perc    1.46 ± 0.09   1.47 ± 0.06  0.92 ± 0.… ns            ***         
## 3 Liver_g       83.58 ± 5.31  75.08 ± 5.09 132.77 ± … ns            ***         
## 4 Liver_perc    2.96 ± 0.11   2.71 ± 0.08  2.76 ± 0.… *             ns          
## 5 Left_kidney_g 12.23 ± 1.05  13.78 ± 1.89 16.37 ± 0… ns            ns          
## 6 Left_kidney_… 0.43 ± 0.02   0.49 ± 0.05  0.34 ± 0.… ns            ns          
## 7 Right_kidney… 11.88 ± 0.99  13.12 ± 1.57 15.92 ± 0… ns            *           
## 8 Right_kidney… 0.42 ± 0.02   0.47 ± 0.04  0.33 ± 0.… ns            ns          
## # ℹ 1 more variable: `WPI vs SF` <chr>
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] tidyr_1.3.0   rstatix_0.7.2 dplyr_1.1.3   readxl_1.4.3 
## 
## loaded via a namespace (and not attached):
##  [1] rstudioapi_0.15.0 knitr_1.45        magrittr_2.0.3    tidyselect_1.2.0 
##  [5] R6_2.5.1          rlang_1.1.2       fastmap_1.1.1     carData_3.0-5    
##  [9] fansi_1.0.6       car_3.1-2         tools_4.2.2       broom_1.0.5      
## [13] xfun_0.40         utf8_1.2.4        cli_3.6.2         withr_3.0.1      
## [17] jquerylib_0.1.4   htmltools_0.5.6.1 abind_1.4-5       yaml_2.3.7       
## [21] digest_0.6.33     tibble_3.2.1      lifecycle_1.0.4   purrr_1.0.2      
## [25] sass_0.4.7        vctrs_0.6.5       glue_1.6.2        cachem_1.0.8     
## [29] evaluate_1.0.1    rmarkdown_2.28    compiler_4.2.2    bslib_0.5.1      
## [33] pillar_1.9.0      cellranger_1.1.0  backports_1.4.1   generics_0.1.3   
## [37] jsonlite_1.8.8    pkgconfig_2.0.3