1. Load libraries

library(tidyverse)
library(readxl)
library(dplyr)
library(car)
library(plyr)
library(kableExtra)
library(hrbrthemes)

knitr::opts_chunk$set(echo = TRUE)

2. Load data

file_path <- "C:/Users/oelle/Documents/Meine Dokumente/Research/Projects/Project_RockLobster/Temperature acclimation/Aerobic scope/Statistics/"

# Read combined data file
  o2_data <- read.csv(paste(file_path, "MO2_data_Tidy_RawData.csv", sep = ""))

3. Tidy data

# Reduce data set to volume and weight
  mass <- o2_data[,c(1:5, 7, 10:11)]
  
# Remove duplicates
  mass <- mass[!duplicated(mass$AntTag),]

4. Plot body mass distribution

# Gplot
   ggplot(mass, aes(x = AcclimationTemperature_C, y = Weight_g, fill = as.factor(AcclimationTemperature_C), alpha = .5)) +
        geom_boxplot(outlier.size = 0) +
        ylab("Body mass (g)") +
        facet_wrap(~ Species, scales = "free")+
        geom_point(pch = 21, position = position_jitterdodge(seed = 100))+
        geom_text(aes(label=AntTag, color = as.factor(AcclimationTemperature_C)), alpha = .8, position = position_jitterdodge(seed = 100))+
        theme_classic()

# 5. Summary statistics

5A Body mass

Grouped by species

# Summary statistics
  weight_stats <- ddply(mass, c("Species"), summarise,
                  mean = mean(Weight_g, na.rm = TRUE),
                  n = length(Weight_g),
                  sd = sd(Weight_g, na.rm = TRUE),
                  ci = qt(0.975, n) * sd / sqrt(n),
                  lower_ci = mean-ci, 
                  upper_ci = mean+ci,
                  min = min(Weight_g),
                  max = max(Weight_g))
   
# Round values
  weight_stats <- weight_stats %>% mutate_if(is.numeric, round, digit = 3)
   
# Display results
  weight_stats %>% kbl(digits = 1, caption = "Descriptvive stats for body mass (g)") %>%  kable_classic()
Descriptvive stats for body mass (g)
Species mean n sd ci lower_ci upper_ci min max
Jasus edwardsii 1074.6 19 205.3 98.6 976.0 1173.2 732 1525
Sagmariasus verreauxi 1138.4 20 180.6 84.3 1054.1 1222.6 716 1452

Grouped by species and acclimation

# Summary statistics
  weight_stats <- ddply(mass, c("Species", "AcclimationTemperature_C"), summarise,
                  mean = mean(Weight_g, na.rm = TRUE),
                  n = length(Weight_g),
                  sd = sd(Weight_g, na.rm = TRUE),
                  ci = qt(0.975, n) * sd / sqrt(n),
                  lower_ci = mean-ci, 
                  upper_ci = mean+ci,
                  min = min(Weight_g),
                  max = max(Weight_g))
   
# Round values
  weight_stats <- weight_stats %>% mutate_if(is.numeric, round, digit = 3)
   
# Display results
  weight_stats %>% kbl(digits = 1, caption = "Descriptvive stats for body mass (g)") %>%  kable_classic()
Descriptvive stats for body mass (g)
Species AcclimationTemperature_C mean n sd ci lower_ci upper_ci min max
Jasus edwardsii 14.0 1133.6 7 227.5 203.4 930.2 1336.9 933 1525
Jasus edwardsii 17.5 1066.1 6 242.4 242.2 824.0 1308.3 732 1291
Jasus edwardsii 21.5 1014.2 6 146.6 146.4 867.8 1160.7 818 1213
Sagmariasus verreauxi 14.0 1135.9 7 151.7 135.6 1000.2 1271.5 973 1452
Sagmariasus verreauxi 17.5 1150.8 7 178.1 159.2 991.6 1310.0 849 1366
Sagmariasus verreauxi 21.5 1126.8 6 240.6 240.4 886.5 1367.2 716 1403

Grouped by species and sex

# Summary statistics
  weight_stats <- ddply(mass, c("Species", "Sex"), summarise,
                  mean = mean(Weight_g, na.rm = TRUE),
                  n = length(Weight_g),
                  sd = sd(Weight_g, na.rm = TRUE),
                  ci = qt(0.975, n) * sd / sqrt(n),
                  lower_ci = mean-ci, 
                  upper_ci = mean+ci,
                  min = min(Weight_g),
                  max = max(Weight_g))
   
# Round values
  weight_stats <- weight_stats %>% mutate_if(is.numeric, round, digit = 3)
   
# Display results
  weight_stats %>% kbl(digits = 1, caption = "Descriptvive stats for body mass (g)") %>%  kable_classic()
Descriptvive stats for body mass (g)
Species Sex mean n sd ci lower_ci upper_ci min max
Jasus edwardsii Female 930.9 7 144.6 129.2 801.7 1060.1 732 1180
Jasus edwardsii Male 1158.4 12 191.7 120.6 1037.8 1279.0 818 1525
Sagmariasus verreauxi Female 1126.7 7 182.0 162.7 964.0 1289.3 849 1403
Sagmariasus verreauxi Male 1144.7 13 187.0 112.1 1032.6 1256.7 716 1452

5B. Sex distribution

# Sex distribution by species
  table(mass$Species, mass$Sex) %>% kbl(digits = 1, caption = "Gender distribution by species") %>%  kable_classic()
Gender distribution by species
Female Male
Jasus edwardsii 7 12
Sagmariasus verreauxi 7 13
# Sex distribution by species and acclimation group
  table(mass$Species, mass$AcclimationTemperature_C, mass$Sex) %>%
    kbl(digits = 1, caption = "Gender distribution by species and acclimation goup") %>%  kable_classic()
Gender distribution by species and acclimation goup
Var1 Var2 Var3 Freq
Jasus edwardsii 14 Female 3
Sagmariasus verreauxi 14 Female 3
Jasus edwardsii 17.5 Female 2
Sagmariasus verreauxi 17.5 Female 2
Jasus edwardsii 21.5 Female 2
Sagmariasus verreauxi 21.5 Female 2
Jasus edwardsii 14 Male 4
Sagmariasus verreauxi 14 Male 4
Jasus edwardsii 17.5 Male 4
Sagmariasus verreauxi 17.5 Male 5
Jasus edwardsii 21.5 Male 4
Sagmariasus verreauxi 21.5 Male 4

6. Compare body mass between acclimation temperatures

# Normality
  ggplot(mass, aes(Weight_g, fill = Species))+
    geom_density( color="#e9ecef", alpha=0.6, position = 'identity')+
    scale_fill_manual(values=c("#69b3a2", "#404080")) +
    theme_ipsum() +
    labs(fill="")

# Levene test for equal variances
  with(mass, leveneTest(Weight_g, Species))
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.5154 0.4773
##       37
# Anova test for differences between acclimation groups and species
    weight_aov <- aov(Weight_g ~ AcclimationTemperature_C * Species * Sex, 
                      data = mass)
  # Test summary
    summary(weight_aov) 
##                                      Df  Sum Sq Mean Sq F value Pr(>F)  
## AcclimationTemperature_C              1   26571   26571   0.779 0.3844  
## Species                               1   39663   39663   1.162 0.2893  
## Sex                                   1  142693  142693   4.182 0.0494 *
## AcclimationTemperature_C:Species      1   19712   19712   0.578 0.4530  
## AcclimationTemperature_C:Sex          1   22622   22622   0.663 0.4217  
## Species:Sex                           1  106609  106609   3.124 0.0870 .
## AcclimationTemperature_C:Species:Sex  1    2422    2422   0.071 0.7917  
## Residuals                            31 1057870   34125                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1