IgA cecal content and pup weight gain analysis

Import data and ANOVA analysis

IgA level in cecal content and weight gain at weaning and between D35 and D49 of the young rabbits are presented below. Effects of age, treatment and their interactions using ANOVA are presented below.

suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(dplyr))
library(broom)
library(car)
library(ggplot2)

IGA<-read.csv("data/identifiant_ech_pds_IGA.csv", sep=";")

Preparing data

data_long <- IGA |> 
  dplyr::select("id","treatment","age","contrast","pds_j2", "pds_j35", "pds_abat", "IgA_prot") |> 
  dplyr::rename(D35_weight = pds_j35,
                D2_weight  = pds_j2,
                slaughter_weight   =pds_abat ) |> 
  dplyr::mutate(weaning_weight_gain = D35_weight-D2_weight,
                D35toD49_weight_gain = slaughter_weight-D35_weight) |> 
 # dplyr::select("id","treatment","age","contrast", D2_weight,IgA_prot:D35toD49_weight_gain) |>  
  dplyr::mutate_at(vars(treatment,age,contrast), as.factor) |> 
  pivot_longer(cols = D2_weight:D35toD49_weight_gain, 
               names_to = "Parameter",
               values_to = "Measure")

options(contrasts = c("contr.sum", "contr.poly"))

ANOVA function

fit_model1 <- function(data){
lm((Measure) ~ treatment*age , data = data) |> 
       Anova(type="III") |> 
         tidy() |> 
  dplyr::select (c(term,p.value)) |> 
  pivot_wider(names_from = term,
              values_from = p.value)
 }

fit_model2 <- function(data){
lm((Measure) ~ treatment , data = data) |> 
       Anova(type="III") |> 
         tidy() |> 
  dplyr::select (c(term,p.value)) |> 
  pivot_wider(names_from = term,
              values_from = p.value)
 }

IgA analysis

data_long |> 
  dplyr::filter(Parameter =="IgA_prot") |> 
  fit_model1()

# A tibble: 1 × 5
  `(Intercept)` treatment      age `treatment:age` Residuals
          <dbl>     <dbl>    <dbl>           <dbl>     <dbl>
1  0.0000000688     0.532 0.000338           0.447        NA

tmp<- data_long |> 
  dplyr::filter(Parameter =="IgA_prot")

  lm((Measure) ~ treatment*age, data=tmp) |>  
  lsmeans::lsmeans(~ treatment*age) |> 
  multcomp::cld(Letters=letters) |> 
  mutate(treatment_age = paste0(treatment,"_", age)) %>% 
  relocate(treatment_age, .after = age) |> 
  dplyr::select(c(treatment_age, .group)) 

  treatment_age .group
5      FFab_d49     a 
4        FF_d49     a 
6        NF_d49     a 
1        FF_d35     ab
3        NF_d35     ab
2      FFab_d35      b

p <- ggplot(IGA, aes(x = treatment, y = log(IgA_prot)))

p<-p+ geom_jitter(
  aes(shape = age, color = age), size = 3.5,
  position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.8)
  ) +
  stat_summary(
    aes(shape=age), fun.data="mean_sdl", fun.args = list(mult=1), 
    size = 0.3, position = position_dodge(0.8)
    )#+
  #scale_color_manual(values =  c("#00AFBB", "#E7B800"))

p+ scale_y_continuous("log(Caecal IgA) (AU)") +
  theme_classic()+
  theme(axis.text=element_text(size=12),
          axis.title=element_text(size=12)) #,face="bold"

Weight gain

No effect of treatment on weaning weight gain or between day 35 and day 49

data_long |> 
  dplyr::filter(Parameter =="weaning_weight_gain") |> 
  fit_model2()

# A tibble: 1 × 3
  `(Intercept)` treatment Residuals
          <dbl>     <dbl>     <dbl>
1      8.25e-48     0.100        NA

data_long |> 
  dplyr::filter(age == "d49", Parameter =="D35toD49_weight_gain") |> 
  fit_model2()

# A tibble: 1 × 3
  `(Intercept)` treatment Residuals
          <dbl>     <dbl>     <dbl>
1      1.28e-19     0.138        NA

Preparing data for weight gain barplot

#mean calculation
tmp1<-data_long |> 
  dplyr::filter(Parameter =="weaning_weight_gain") |>
  dplyr::group_by(treatment) |> 
  dplyr::summarise(moy_trt=mean(Measure/33),sd_trt=sd(Measure/33),sem = sd(Measure/33) / sqrt(length(Measure/33))) |> 
  mutate(traits = "weaning_weight_gain")

tmp2<-data_long |> 
  dplyr::filter(Parameter =="D35toD49_weight_gain") |>
  dplyr::group_by(treatment) |> 
  dplyr::summarise(moy_trt=mean(Measure/15),sd_trt=sd(Measure/15),sem = sd(Measure/15) / sqrt(length(Measure/15))) |> 
  mutate(traits = "D35toD49_weight_gain")
tmp<-rbind(tmp1,tmp2)
  tmp

# A tibble: 6 × 5
  treatment moy_trt sd_trt   sem traits              
  <fct>       <dbl>  <dbl> <dbl> <chr>               
1 FF           25.5   3.92 0.899 weaning_weight_gain 
2 FFab         25.0   2.89 0.663 weaning_weight_gain 
3 NF           23.1   3.40 0.801 weaning_weight_gain 
4 FF           27.2  29.9  6.87  D35toD49_weight_gain
5 FFab         26.4  29.6  6.79  D35toD49_weight_gain
6 NF           24.0  25.8  6.08  D35toD49_weight_gain

barplot

p <- ggplot(tmp, aes(x = traits, y = moy_trt, fill=treatment))

p<-p+ geom_bar(stat="identity", position=position_dodge()) +
  geom_errorbar(aes(ymin=moy_trt, ymax=moy_trt+sd_trt), width=.2,
                 position=position_dodge(.9))+
  labs(x="", y = "Daily weight gain")
  
  
p+ 
  theme_classic()+
  theme(axis.text=element_text(size=12),
          axis.title=element_text(size=12)) #,face="bold"