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