0.setup(package, directory)

library(readxl)
library(ggplot2)
library(rstatix)
library(tidyverse)
library(lawstat) ##levene test
library(ggpubr)


df<-read_xlsx("Homepen.xlsx")
df <- df %>% 
  mutate(df, unique_pen=as.factor(paste(Replicate,Pen)),pen_pos=as.numeric(substr(df$Pen,3,3))) %>% 
  convert_as_factor(Replicate,Pen,Stress, unique_pen)
head(df)
## # A tibble: 6 × 26
##   Replicate Pen   Pig   Stress   FT2  INC2   NC2   BN2  SLP2   SP2   FE2   LD2
##   <fct>     <fct> <chr> <fct>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1         ML5   3     SV       670     0    17     0     0     0  1531  9547
## 2 1         ML6   2     SV       296     5    27     0     0     0   789 12105
## 3 1         ML6   3     SR       246     5    32     0     0     0   233 11089
## 4 1         ML6   5     SV       393     5    30     0     0     0  2386  9848
## 5 1         ML6   6     SR        39     7   103     0     0     0   242  9118
## 6 1         MR6   7     SR         0    15   101     0     0     0  2788  8889
## # … with 14 more variables: EX2 <dbl>, DR2 <dbl>, FT5 <dbl>, INC5 <dbl>,
## #   NC5 <dbl>, BN5 <dbl>, SLP5 <dbl>, SP5 <dbl>, FE5 <dbl>, LD5 <dbl>,
## #   EX5 <dbl>, DR5 <dbl>, unique_pen <fct>, pen_pos <dbl>
df2 <- read_xlsx("Homepen2.xlsx")
df2 %>% 
  convert_as_factor(ID,Stress,Age) 
## # A tibble: 104 × 13
##    ID    Stress Age      FT   INC    NC    BN   SLP    SP    FE    LD    EX
##    <fct> <fct>  <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1 1ML53 SV     2       670     0    17     0     0     0  1531  9547  1405
##  2 1ML62 SV     2       296     5    27     0     0     0   789 12105   329
##  3 1ML63 SR     2       246     5    32     0     0     0   233 11089   860
##  4 1ML65 SV     2       393     5    30     0     0     0  2386  9848   353
##  5 1ML66 SR     2        39     7   103     0     0     0   242  9118  2908
##  6 1MR67 SR     2         0    15   101     0     0     0  2788  8889   980
##  7 1MR72 SR     2      1053     0     7     0     3   608     0 11754   274
##  8 1MR76 SV     2         0     0    33    12     0     0  1210  7937   535
##  9 2ML28 SV     2       670    20   127     0     0     0  1403  6413  3212
## 10 2MR22 SV     2        55     7   126     0     0     0   133  6760  4966
## # … with 94 more rows, and 1 more variable: DR <dbl>
str(df2)
## tibble [104 × 13] (S3: tbl_df/tbl/data.frame)
##  $ ID    : chr [1:104] "1ML53" "1ML62" "1ML63" "1ML65" ...
##  $ Stress: chr [1:104] "SV" "SV" "SR" "SV" ...
##  $ Age   : num [1:104] 2 2 2 2 2 2 2 2 2 2 ...
##  $ FT    : num [1:104] 670 296 246 393 39 ...
##  $ INC   : num [1:104] 0 5 5 5 7 15 0 0 20 7 ...
##  $ NC    : num [1:104] 17 27 32 30 103 101 7 33 127 126 ...
##  $ BN    : num [1:104] 0 0 0 0 0 0 0 12 0 0 ...
##  $ SLP   : num [1:104] 0 0 0 0 0 0 3 0 0 0 ...
##  $ SP    : num [1:104] 0 0 0 0 0 0 608 0 0 0 ...
##  $ FE    : num [1:104] 1531 789 233 2386 242 ...
##  $ LD    : num [1:104] 9547 12105 11089 9848 9118 ...
##  $ EX    : num [1:104] 1405 329 860 353 2908 ...
##  $ DR    : num [1:104] 15 6 24 42 15 22 13 17 16 25 ...
head(df2)
## # A tibble: 6 × 13
##   ID    Stress   Age    FT   INC    NC    BN   SLP    SP    FE    LD    EX    DR
##   <chr> <chr>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1ML53 SV         2   670     0    17     0     0     0  1531  9547  1405    15
## 2 1ML62 SV         2   296     5    27     0     0     0   789 12105   329     6
## 3 1ML63 SR         2   246     5    32     0     0     0   233 11089   860    24
## 4 1ML65 SV         2   393     5    30     0     0     0  2386  9848   353    42
## 5 1ML66 SR         2    39     7   103     0     0     0   242  9118  2908    15
## 6 1MR67 SR         2     0    15   101     0     0     0  2788  8889   980    22

1. SV vs. SR all factors

SV<-subset(df, Stress=="SV")
SR<-subset(df, Stress=="SR")


# 1.1 Agonistic Behaviors(Day2)
# 1.1.1 FT
#boxplot
p <- ggplot(df, aes(x=Stress, y=FT2)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=2)
p

#blocked anova test/equal-variance/block o
levene.test(df$FT2,df$Stress, location="mean")
## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  df$FT2
## Test Statistic = 0.60036, p-value = 0.4421
summary(aov(FT2~Stress+unique_pen, data=df))
##             Df  Sum Sq Mean Sq F value  Pr(>F)   
## Stress       1   31853   31853   0.239 0.62842   
## unique_pen  20 6819120  340956   2.559 0.00973 **
## Residuals   30 3996740  133225                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 1.1.2 INC
#boxplot
p <- ggplot(df, aes(x=Stress, y=INC2)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=2)
p

#blocked anova test/equal-variance/block o
levene.test(df$INC2,df$Stress, location="mean")
## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  df$INC2
## Test Statistic = 0.76593, p-value = 0.3857
summary(aov(INC2~Stress+unique_pen, data=df))
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## Stress       1    185   184.7   2.199 0.14856    
## unique_pen  20   5953   297.6   3.543 0.00089 ***
## Residuals   30   2520    84.0                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 1.2 Non-Agonistic Behaviors(Day2)
# 1.2.1 NC
#boxplot
p <- ggplot(df, aes(x=Stress, y=NC2)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=2)
p

#blocked anova test/equal-variance/block o & stress o *****
levene.test(df$NC2,df$Stress, location="mean")
## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  df$NC2
## Test Statistic = 0.74617, p-value = 0.3918
summary(aov(NC2~Stress+unique_pen, data=df))
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Stress       1   9694    9694   5.396 0.027151 *  
## unique_pen  20 135404    6770   3.768 0.000534 ***
## Residuals   30  53897    1797                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 1.2.2 BN - zero inflation check
#boxplot
p <- ggplot(df, aes(x=Stress, y=BN2)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=2)
p

#blocked anova test/equal-variance/block o
levene.test(df$BN2,df$Stress, location="mean")
## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  df$BN2
## Test Statistic = 0.010931, p-value = 0.9172
summary(aov(BN2~Stress+unique_pen, data=df))
##             Df Sum Sq Mean Sq F value Pr(>F)  
## Stress       1    0.0   0.000   0.000  1.000  
## unique_pen  20  242.3  12.113   2.344  0.017 *
## Residuals   30  155.1   5.168                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 1.3 Daily Behaviors(Day2)
# 1.3.1 SLP -zero inflation
#boxplot
p <- ggplot(df, aes(x=Stress, y=SLP2)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=2)
p

#blocked anova test/unequal-variance/
levene.test(df$SLP2,df$Stress, location="mean")
## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  df$SLP2
## Test Statistic = 5.3523, p-value = 0.02485
summary(aov(SLP2~Stress+unique_pen, data=df)) ##unequal
##             Df Sum Sq Mean Sq F value Pr(>F)
## Stress       1    246   245.6   1.307  0.262
## unique_pen  20   3221   161.1   0.857  0.634
## Residuals   30   5636   187.9
oneway.test(SLP2~Stress, data=df)
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  SLP2 and Stress
## F = 1.3862, num df = 1.000, denom df = 25.094, p-value = 0.2501
# 1.3.2 SP 
#boxplot
p <- ggplot(df, aes(x=Stress, y=SP2)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=2)
p

#blocked anova test/equal-variance/stress, block x
levene.test(df$SP2,df$Stress, location="mean")
## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  df$SP2
## Test Statistic = 1.5379, p-value = 0.2207
summary(aov(SP2~Stress+unique_pen, data=df)) ##unequal
##             Df Sum Sq Mean Sq F value Pr(>F)  
## Stress       1   7416    7416   0.795 0.3798  
## unique_pen  20 349005   17450   1.870 0.0588 .
## Residuals   30 279970    9332                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 1.3.3 FE 
#boxplot
p <- ggplot(df, aes(x=Stress, y=FE2)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=2)
p

#blocked anova test/equal-variance/stress, block x
levene.test(df$FE2,df$Stress, location="mean")
## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  df$FE2
## Test Statistic = 2.9629, p-value = 0.09138
summary(aov(FE2~Stress+unique_pen, data=df)) ##unequal
##             Df   Sum Sq Mean Sq F value Pr(>F)
## Stress       1   119904  119904   0.102  0.751
## unique_pen  20 34534613 1726731   1.476  0.163
## Residuals   30 35100900 1170030
# 1.3.4 LD 
#boxplot
p <- ggplot(df, aes(x=Stress, y=LD2)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=2)
p

#blocked anova test/equal-variance/ block o
levene.test(df$LD2,df$Stress, location="mean")
## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  df$LD2
## Test Statistic = 0.0099597, p-value = 0.9209
summary(aov(LD2~Stress+unique_pen, data=df))
##             Df   Sum Sq Mean Sq F value Pr(>F)  
## Stress       1  5235577 5235577   2.405 0.1314  
## unique_pen  20 97074164 4853708   2.230 0.0228 *
## Residuals   30 65309753 2176992                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 1.3.5 EX 
#boxplot
p <- ggplot(df, aes(x=Stress, y=EX2)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=2)
p

#blocked anova test/equal-variance/block o
levene.test(df$EX2,df$Stress, location="mean")
## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  df$EX2
## Test Statistic = 2.4928, p-value = 0.1207
summary(aov(EX2~Stress+unique_pen, data=df)) 
##             Df   Sum Sq Mean Sq F value Pr(>F)  
## Stress       1  1151176 1151176   1.077  0.308  
## unique_pen  20 42521908 2126095   1.989  0.043 *
## Residuals   30 32069717 1068991                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 1.3.6 DR 
#boxplot
p <- ggplot(df, aes(x=Stress, y=DR2)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=2)
p

#blocked anova test/equal-variance/block,Stress x
levene.test(df$DR2,df$Stress, location="mean")
## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  df$DR2
## Test Statistic = 0.32203, p-value = 0.5729
summary(aov(DR2~Stress+unique_pen, data=df)) 
##             Df Sum Sq Mean Sq F value Pr(>F)
## Stress       1   32.3   32.33   0.440  0.512
## unique_pen  20 1600.6   80.03   1.088  0.408
## Residuals   30 2205.8   73.53
##################################################

# 1.4 Agonistic Behaviors(Day5)
# 1.4.1 FT
#boxplot
p <- ggplot(df, aes(x=Stress, y=FT5)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=2)
p

#blocked anova test/equal-variance/stress,block o
levene.test(df$FT5,df$Stress, location="mean")
## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  df$FT5
## Test Statistic = 0.15955, p-value = 0.6913
summary(aov(FT5~Stress+unique_pen, data=df)) 
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Stress       1   3233    3233   4.805   0.0363 *  
## unique_pen  20  61483    3074   4.569 9.69e-05 ***
## Residuals   30  20184     673                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 1.4.2 INC
#boxplot
p <- ggplot(df, aes(x=Stress, y=INC5)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=2)
p

#blocked anova test/unequal-variance/
levene.test(df$INC5,df$Stress, location="mean")
## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  df$INC5
## Test Statistic = 5.4916, p-value = 0.02313
summary(aov(INC5~Stress+unique_pen, data=df)) ## unequal
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Stress       1    291   290.9   4.554   0.0411 *  
## unique_pen  20   8312   415.6   6.505 2.85e-06 ***
## Residuals   30   1917    63.9                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
oneway.test(INC5~Stress, data=df)
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  INC5 and Stress
## F = 1.4221, num df = 1.000, denom df = 36.694, p-value = 0.2407
# 1.5 Non-Agonistic Behaviors(Day5)
# 1.5.1 NC
#boxplot
p <- ggplot(df, aes(x=Stress, y=NC5)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=2)
p

#blocked anova test/equal-variance/block o
levene.test(df$NC5,df$Stress, location="mean")
## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  df$NC5
## Test Statistic = 2.2443, p-value = 0.1404
summary(aov(NC5~Stress+unique_pen, data=df)) 
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Stress       1    438     438   0.760     0.39    
## unique_pen  20  63970    3198   5.545 1.49e-05 ***
## Residuals   30  17306     577                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 1.5.2 BN - zero inflation check
#boxplot
p <- ggplot(df, aes(x=Stress, y=BN5)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=2)
p

#blocked anova test/equal-variance/block o
levene.test(df$BN5,df$Stress, location="mean")
## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  df$BN5
## Test Statistic = 0.73966, p-value = 0.3939
summary(aov(BN5~Stress+unique_pen, data=df)) 
##             Df Sum Sq Mean Sq F value Pr(>F)  
## Stress       1   0.17  0.1731   0.252  0.620  
## unique_pen  20  33.25  1.6624   2.416  0.014 *
## Residuals   30  20.64  0.6879                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 1.6 Daily Behaviors(Day2)
# 1.6.1 SLP -zero inflation
#boxplot
p <- ggplot(df, aes(x=Stress, y=SLP5)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=2)
p

#blocked anova test/unequal-variance/
levene.test(df$SLP5,df$Stress, location="mean")
## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  df$SLP5
## Test Statistic = 10.203, p-value = 0.002428
summary(aov(SLP5~Stress+unique_pen, data=df)) ##unequal
##             Df Sum Sq Mean Sq F value Pr(>F)
## Stress       1   96.9   96.94   2.788  0.105
## unique_pen  20  670.8   33.54   0.964  0.524
## Residuals   30 1043.3   34.78
oneway.test(SLP5~Stress, data=df)
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  SLP5 and Stress
## F = 2.8279, num df = 1.000, denom df = 28.682, p-value = 0.1035
# 1.6.2 SP - zero inflated
#boxplot
p <- ggplot(df, aes(x=Stress, y=SP5)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=2)
p

#blocked anova test/equal-variance
levene.test(df$SP5,df$Stress, location="mean")
## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  df$SP5
## Test Statistic = 2.1816, p-value = 0.1459
summary(aov(SP5~Stress+unique_pen, data=df)) 
##             Df Sum Sq Mean Sq F value Pr(>F)  
## Stress       1    746   746.3   0.880 0.3557  
## unique_pen  20  29630  1481.5   1.747 0.0813 .
## Residuals   30  25443   848.1                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 1.6.3 FE 
#boxplot
p <- ggplot(df, aes(x=Stress, y=FE5)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=2)
p

#blocked anova test/equal-variance/
levene.test(df$FE5,df$Stress, location="mean")
## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  df$FE5
## Test Statistic = 3.7686, p-value = 0.05787
summary(aov(FE5~Stress+unique_pen, data=df)) 
##             Df   Sum Sq Mean Sq F value Pr(>F)
## Stress       1   497646  497646   0.443  0.511
## unique_pen  20 25740003 1287000   1.145  0.361
## Residuals   30 33728657 1124289
# 1.6.4 LD 
#boxplot
p <- ggplot(df, aes(x=Stress, y=LD5)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=2)
p

#blocked anova test/equal-variance/block o
levene.test(df$LD5,df$Stress, location="mean")
## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  df$LD5
## Test Statistic = 0.62986, p-value = 0.4312
summary(aov(LD5~Stress+unique_pen, data=df)) 
##             Df    Sum Sq Mean Sq F value   Pr(>F)    
## Stress       1    772261  772261   0.329 0.570467    
## unique_pen  20 170133142 8506657   3.625 0.000738 ***
## Residuals   30  70398039 2346601                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 1.6.5 EX 
#boxplot
p <- ggplot(df, aes(x=Stress, y=EX5)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=2)
p

#blocked anova test/equal-variance/block o
levene.test(df$EX5,df$Stress, location="mean")
## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  df$EX5
## Test Statistic = 0.83665, p-value = 0.3647
summary(aov(EX5~Stress+unique_pen, data=df)) 
##             Df   Sum Sq Mean Sq F value  Pr(>F)   
## Stress       1    18545   18545   0.065 0.80042   
## unique_pen  20 18978879  948944   3.329 0.00147 **
## Residuals   30  8551725  285058                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 1.6.6 DR 
#boxplot
p <- ggplot(df, aes(x=Stress, y=DR5)) + 
  geom_boxplot(outlier.colour="red", outlier.shape=8,
               outlier.size=2)
p

#blocked anova test/equal-variance/
levene.test(df$DR5,df$Stress, location="mean")
## 
##  Classical Levene's test based on the absolute deviations from the mean
##  ( none not applied because the location is not set to median )
## 
## data:  df$DR5
## Test Statistic = 0.039832, p-value = 0.8426
summary(aov(DR5~Stress+unique_pen, data=df)) 
##             Df Sum Sq Mean Sq F value Pr(>F)
## Stress       1    0.1    0.08   0.002  0.969
## unique_pen  20 1415.2   70.76   1.425  0.186
## Residuals   30 1490.0   49.67

2. SV vs. SR all factors with ages

# 2.1 
# 2.1.1 FT *** interaction x, group x, age o
summary(aov(FT~Stress*Age+Error(ID/Age),data=df2)) 
## 
## Error: ID
##           Df  Sum Sq Mean Sq F value Pr(>F)
## Stress     1   27690   27690   0.244  0.624
## Residuals 50 5678396  113568               
## 
## Error: ID:Age
##            Df  Sum Sq Mean Sq F value   Pr(>F)    
## Age         1 2024145 2024145  19.392 5.61e-05 ***
## Stress:Age  1    7395    7395   0.071    0.791    
## Residuals  50 5219130  104383                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
interaction.plot(df2$Age,df2$Stress,df2$FT,
                 fun=mean,type="b",pch=c(2,4),col=c(3,4),xlab="Time",ylab="Mean of FT", trace.label = "Stress", legend = F) 
legend("topright",legend=c("SR","SV"),pch=c(2,4),col=c(3,4), bg="gray90")

ggplot(df2, aes(x=Stress, y=FT, color=as.factor(Age)))+
  geom_boxplot()+
  ggtitle("FT distribution by Stress and Age")+
  theme(plot.title=element_text(face="bold", hjust=0.5))+
  scale_color_discrete(name="Age")

# 2.1.2 INC *** interaction o, stress x, age x
summary(aov(INC~Stress*Age+Error(ID/Age),data=df2)) 
## 
## Error: ID
##           Df Sum Sq Mean Sq F value Pr(>F)
## Stress     1      6    6.01   0.024  0.879
## Residuals 50  12754  255.09               
## 
## Error: ID:Age
##            Df Sum Sq Mean Sq F value Pr(>F)  
## Age         1     83    83.2   0.699 0.4071  
## Stress:Age  1    470   469.6   3.948 0.0524 .
## Residuals  50   5948   119.0                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
interaction.plot(df2$Age,df2$Stress,df2$INC,
                 fun=mean,type="b",pch=c(2,4),col=c(3,4),xlab="Time",ylab="Mean of INC", trace.label = "Stress", legend = F) 
legend("topleft",legend=c("SR","SV"),pch=c(2,4),col=c(3,4), bg="gray90")

ggplot(df2, aes(x=Stress, y=INC, color=as.factor(Age)))+
  geom_boxplot()+
  ggtitle("INC distribution by Stress and Age")+
  theme(plot.title=element_text(face="bold", hjust=0.5))+
  scale_color_discrete(name="Age")

# 2.2
# 2.2.1 NC *interaction x, stress x, age o
summary(aov(NC~Stress*Age+Error(ID/Age),data=df2)) 
## 
## Error: ID
##           Df Sum Sq Mean Sq F value Pr(>F)
## Stress     1   7128    7128   2.052  0.158
## Residuals 50 173722    3474               
## 
## Error: ID:Age
##            Df Sum Sq Mean Sq F value Pr(>F)  
## Age         1  12916   12916   6.668 0.0128 *
## Stress:Age  1   3005    3005   1.551 0.2188  
## Residuals  50  96855    1937                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
interaction.plot(df2$Age,df2$Stress,df2$NC,
                 fun=mean,type="b",pch=c(2,4),col=c(3,4),xlab="Time",ylab="Mean of NC", trace.label = "Stress", legend = F) 
legend("topright",legend=c("SR","SV"),pch=c(2,4),col=c(3,4), bg="gray90")

ggplot(df2, aes(x=Stress, y=NC, color=as.factor(Age)))+
  geom_boxplot()+
  ggtitle("NC distribution by Stress and Age")+
  theme(plot.title=element_text(face="bold", hjust=0.5))+
  scale_color_discrete(name="Age")

# 2.2.2 BN *interaction x, stress x, age x - zero-inflated
summary(aov(BN~Stress*Age+Error(ID/Age),data=df2)) 
## 
## Error: ID
##           Df Sum Sq Mean Sq F value Pr(>F)
## Stress     1   0.09   0.087   0.016    0.9
## Residuals 50 268.79   5.376               
## 
## Error: ID:Age
##            Df Sum Sq Mean Sq F value Pr(>F)
## Age         1   7.01   7.010   1.921  0.172
## Stress:Age  1   0.09   0.087   0.024  0.878
## Residuals  50 182.40   3.648
interaction.plot(df2$Age,df2$Stress,df2$BN,
                 fun=mean,type="b",pch=c(2,4),col=c(3,4),xlab="Time",ylab="Mean of BN", trace.label = "Stress", legend = F) 
legend("topright",legend=c("SR","SV"),pch=c(2,4),col=c(3,4), bg="gray90")

ggplot(df2, aes(x=Stress, y=BN, color=as.factor(Age)))+
  geom_boxplot()+
  ggtitle("BN distribution by Stress and Age")+
  theme(plot.title=element_text(face="bold", hjust=0.5))+
  scale_color_discrete(name="Age")

#2.3
#2.3.1 SLP *interaction o, stress x, age x - zero-inflated
summary(aov(SLP~Stress*Age+Error(ID/Age),data=df2)) 
## 
## Error: ID
##           Df Sum Sq Mean Sq F value Pr(>F)
## Stress     1     17   16.96   0.162  0.689
## Residuals 50   5231  104.62               
## 
## Error: ID:Age
##            Df Sum Sq Mean Sq F value Pr(>F)  
## Age         1      5     4.7   0.044  0.835  
## Stress:Age  1    326   325.5   3.048  0.087 .
## Residuals  50   5340   106.8                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
interaction.plot(df2$Age,df2$Stress,df2$SLP,
                 fun=mean,type="b",pch=c(2,4),col=c(3,4),xlab="Time",ylab="Mean of SLP", trace.label = "Stress", legend = F) 
legend("topright",legend=c("SR","SV"),pch=c(2,4),col=c(3,4), bg="gray90")

ggplot(df2, aes(x=Stress, y=SLP, color=as.factor(Age)))+
  geom_boxplot()+
  ggtitle("SLP distribution by Stress and Age")+
  theme(plot.title=element_text(face="bold", hjust=0.5))+
  scale_color_discrete(name="Age")

#2.3.2 SP *interaction x, stress x, age o : zero-inflated
summary(aov(SP~Stress*Age+Error(ID/Age),data=df2)) 
## 
## Error: ID
##           Df Sum Sq Mean Sq F value Pr(>F)
## Stress     1   6434    6434   0.743  0.393
## Residuals 50 433130    8663               
## 
## Error: ID:Age
##            Df Sum Sq Mean Sq F value Pr(>F)  
## Age         1  17889   17889   3.565 0.0648 .
## Stress:Age  1   1729    1729   0.344 0.5599  
## Residuals  50 250917    5018                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
interaction.plot(df2$Age,df2$Stress,df2$SP,
                 fun=mean,type="b",pch=c(2,4),col=c(3,4),xlab="Time",ylab="Mean of SP", trace.label = "Stress", legend = F) 
legend("topright",legend=c("SR","SV"),pch=c(2,4),col=c(3,4), bg="gray90")

ggplot(df2, aes(x=Stress, y=SP, color=as.factor(Age)))+
  geom_boxplot()+
  ggtitle("SP distribution by Stress and Age")+
  theme(plot.title=element_text(face="bold", hjust=0.5))+
  scale_color_discrete(name="Age")

#2.3.3 FE * interaction x, stress x, age o
summary(aov(FE~Stress*Age+Error(ID/Age),data=df2)) 
## 
## Error: ID
##           Df   Sum Sq Mean Sq F value Pr(>F)
## Stress     1   553049  553049   0.386  0.537
## Residuals 50 71725761 1434515               
## 
## Error: ID:Age
##            Df   Sum Sq  Mean Sq F value  Pr(>F)   
## Age         1 10261521 10261521   8.942 0.00432 **
## Stress:Age  1    64501    64501   0.056 0.81356   
## Residuals  50 57378412  1147568                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
interaction.plot(df2$Age,df2$Stress,df2$FE,
                 fun=mean,type="b",pch=c(2,4),col=c(3,4),xlab="Time",ylab="Mean of FE", trace.label = "Stress", legend = F) 
legend("topleft",legend=c("SR","SV"),pch=c(2,4),col=c(3,4), bg="gray90")

ggplot(df2, aes(x=Stress, y=FE, color=as.factor(Age)))+
  geom_boxplot()+
  ggtitle("FE distribution by Stress and Age")+
  theme(plot.title=element_text(face="bold", hjust=0.5))+
  scale_color_discrete(name="Age")

#2.3.4 LD * interaction x, stress x, age o
summary(aov(LD~Stress*Age+Error(ID/Age),data=df2)) 
## 
## Error: ID
##           Df    Sum Sq Mean Sq F value Pr(>F)
## Stress     1   5014698 5014698   1.023  0.317
## Residuals 50 245107906 4902158               
## 
## Error: ID:Age
##            Df    Sum Sq  Mean Sq F value   Pr(>F)    
## Age         1  51756891 51756891  16.399 0.000178 ***
## Stress:Age  1    993140   993140   0.315 0.577335    
## Residuals  50 157807192  3156144                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
interaction.plot(df2$Age,df2$Stress,df2$LD,
                 fun=mean,type="b",pch=c(2,4),col=c(3,4),xlab="Time",ylab="Mean of LD", trace.label = "Stress", legend = F) 
legend("topleft",legend=c("SR","SV"),pch=c(2,4),col=c(3,4), bg="gray90")

ggplot(df2, aes(x=Stress, y=LD, color=as.factor(Age)))+
  geom_boxplot()+
  ggtitle("LD distribution by Stress and Age")+
  theme(plot.title=element_text(face="bold", hjust=0.5))+
  scale_color_discrete(name="Age")

#2.3.5 EX * interaction x, stress x, age o
summary(aov(EX~Stress*Age+Error(ID/Age),data=df2)) 
## 
## Error: ID
##           Df   Sum Sq Mean Sq F value Pr(>F)
## Stress     1   730971  730971   0.539  0.466
## Residuals 50 67866791 1357336               
## 
## Error: ID:Age
##            Df   Sum Sq Mean Sq F value   Pr(>F)    
## Age         1  9955378 9955378   14.53 0.000379 ***
## Stress:Age  1   438750  438750    0.64 0.427349    
## Residuals  50 34255437  685109                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
interaction.plot(df2$Age,df2$Stress,df2$EX,
                 fun=mean,type="b",pch=c(2,4),col=c(3,4),xlab="Time",ylab="Mean of EX", trace.label = "Stress", legend = F) 
legend("topright",legend=c("SR","SV"),pch=c(2,4),col=c(3,4), bg="gray90")

ggplot(df2, aes(x=Stress, y=EX, color=as.factor(Age)))+
  geom_boxplot()+
  ggtitle("EX distribution by Stress and Age")+
  theme(plot.title=element_text(face="bold", hjust=0.5))+
  scale_color_discrete(name="Age")

#2.3.6 DR *interaction x, stress x, age o
summary(aov(DR~Stress*Age+Error(ID/Age),data=df2)) 
## 
## Error: ID
##           Df Sum Sq Mean Sq F value Pr(>F)
## Stress     1     18   17.78   0.213  0.647
## Residuals 50   4180   83.61               
## 
## Error: ID:Age
##            Df Sum Sq Mean Sq F value   Pr(>F)    
## Age         1 1211.8  1211.8  23.938 1.08e-05 ***
## Stress:Age  1   14.6    14.6   0.289    0.593    
## Residuals  50 2531.1    50.6                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
interaction.plot(df2$Age,df2$Stress,df2$DR,
                 fun=mean,type="b",pch=c(2,4),col=c(3,4),xlab="Time",ylab="Mean of DR", trace.label = "Stress", legend = F) 
legend("topright",legend=c("SR","SV"),pch=c(2,4),col=c(3,4), bg="gray90")

ggplot(df2, aes(x=Stress, y=DR, color=as.factor(Age)))+
  geom_boxplot()+
  ggtitle("DR distribution by Stress and Age")+
  theme(plot.title=element_text(face="bold", hjust=0.5))+
  scale_color_discrete(name="Age")