1. Load libraries
library(tidyverse)
library(readxl)
library(dplyr)
library(plyr)
library(openxlsx)
library(pbkrtest)
library(psycho)
library(pander)
library(ggpubr)
library(knitr)
library(kableExtra) # For Markdown tables
library(report) # textual summaries of analysis
library(lme4) # to fit the linear mixed effect model
library(lmerTest)
library(fitdistrplus) # to fit probability distributions
library(emmeans) # to compare contrasts
library(bbmle) # for AICtab
library(sjPlot) # for plot_model function
library(DHARMa) # for GLMM model diagnostics
knitr::opts_chunk$set(echo = TRUE)
2. Load data
# Select trait
trait <- "escape_speed_cmSec"
# File path
file_path <- "C:/Users/oelle/Documents/Meine Dokumente/Research/Projects/Project_RockLobster/Temperature acclimation/Exhaustion/Analysis/LinearMixedEffectModels/EscapesTotal/"
# Load tidy raw data in long format
Escape_MO2_data_long <- read_xlsx("C:/Users/oelle/Documents/Meine Dokumente/Research/Projects/Project_RockLobster/Temperature acclimation/Exhaustion/Analysis/Escape_MO2_MergedData_long.xlsx", sheet = 1)
3. Process data
# Shorten column names
names(Escape_MO2_data_long)[c(4,5)] <- c("AccTemp", "ExpTemp")
# Define fixed effects as factors
Escape_MO2_data_long$AntTag <- as.factor(Escape_MO2_data_long$AntTag)
Escape_MO2_data_long$Species <- as.factor(Escape_MO2_data_long$Species)
Escape_MO2_data_long$AccTemp <- as.factor(Escape_MO2_data_long$AccTemp)
Escape_MO2_data_long$ExpTemp <- as.factor(Escape_MO2_data_long$ExpTemp)
Escape_MO2_data_long$Sex <- as.factor(Escape_MO2_data_long$Sex)
# Subset data by trait
sub_trait <- subset(Escape_MO2_data_long, Trait == trait)
# Remove any columns with NAs
sub_trait <- sub_trait %>% drop_na("value")
4. Find best probability distribution of data to check assumptions for linear model
If not normal distributed find best distribution fit using fitdistrplus package
# Plot data versus a range of distrubutions
descdist(sub_trait$value, boot = 1000)

## summary statistics
## ------
## min: 32.98261 max: 117.4075
## median: 71.01666
## mean: 72.0916
## estimated sd: 14.46943
## estimated skewness: 0.3166321
## estimated kurtosis: 3.181273
# Fit various distributions
nd <- fitdist(sub_trait$value, "norm")
## $start.arg
## $start.arg$mean
## [1] 72.0916
##
## $start.arg$sd
## [1] 14.39833
##
##
## $fix.arg
## NULL
gd <- fitdist(sub_trait$value, "gamma")
## $start.arg
## $start.arg$shape
## [1] 25.06947
##
## $start.arg$rate
## [1] 0.3477446
##
##
## $fix.arg
## NULL
wd <- fitdist(sub_trait$value, "weibull")
## $start.arg
## $start.arg$shape
## [1] 5.84722
##
## $start.arg$scale
## [1] 77.8969
##
##
## $fix.arg
## NULL
lnd <- fitdist(sub_trait$value, "lnorm")
## $start.arg
## $start.arg$meanlog
## [1] 4.257562
##
## $start.arg$sdlog
## [1] 0.2042172
##
##
## $fix.arg
## NULL
# Plot distributions for the data
plot.legend <- c("Normal", "gamma", "Weibull", "LogNormal")
denscomp(list(nd, gd, wd, lnd), legendtext = plot.legend)

# QQ plot
qqcomp(list(nd, gd, wd, lnd), legendtext = plot.legend)

# Perform goodness of fit test to select best distribution
gofstat(list(nd, gd, wd, lnd), fitnames = c("Normal", "gamma", "Weibull", "LogNormal"))
## Goodness-of-fit statistics
## Normal gamma Weibull LogNormal
## Kolmogorov-Smirnov statistic 0.06611752 0.03928513 0.08176727 0.04166618
## Cramer-von Mises statistic 0.05791803 0.02367328 0.12194106 0.02667363
## Anderson-Darling statistic 0.38959069 0.20361406 0.84976333 0.23322871
##
## Goodness-of-fit criteria
## Normal gamma Weibull LogNormal
## Akaike's Information Criterion 837.5544 836.2694 844.1822 837.9376
## Bayesian Information Criterion 842.8043 841.5194 849.4321 843.1876
–> Anderson-Darling statistics gives good measure for fit for both middle and tail of distribution, the lower the better
–> Data fit best to various distribution. Visual analysis of histogram and QQ plot show that data still correspond well to normal distribution
5. Plot data
# Experimental temperature versus acclimation
ggplot(sub_trait, aes(x = ExpTemp, y = value, fill = AccTemp, alpha = .5)) +
geom_boxplot(outlier.size = 0) +
ylab(trait) +
facet_wrap(~ Species, scales = "free")+
geom_point(pch = 21, position = position_jitterdodge(seed = 100))+
geom_text(aes(label=AntTag, color = AccTemp), alpha = .8, position = position_jitterdodge(seed = 100))

ggsave("Escapes_speed_AccTemp.pdf", width = 14, height = 8)
# Line Plot comparing both species
ggline(sub_trait, x = "ExpTemp", y = "value",
add = c("mean_ci", "jitter"), # Add mean_se and jitter points
# add.params = list(size = 0.7), # Add point size
# label = "AntTag", # Add point labels
# #label.select = list(top.up = 2), # show only labels for the top 2 points
# font.label = list(color = "Species"), # Color labels by .y., here gene names
# repel = TRUE, # Use repel to avoid labels overplotting
color = "Species",
palette = "jco")

# Calculate means across groups
value_means <- ddply(sub_trait, c("Species", "ExpTemp", "AccTemp"), summarise,
mean = mean(value, na.rm = TRUE),
n = length(value),
sd = sd(value, na.rm = TRUE),
ci = qt(0.975, n) * sd / sqrt(n),
lower_ci = mean-ci,
upper_ci = mean+ci)
print(value_means)
## Species ExpTemp AccTemp mean n sd ci
## 1 Jasus edwardsii 14 14 66.82045 6 8.842155 8.832849
## 2 Jasus edwardsii 14 17.5 62.59641 6 6.853407 6.846194
## 3 Jasus edwardsii 14 21.5 57.72742 6 5.053764 5.048445
## 4 Jasus edwardsii 17.5 14 67.24931 6 10.756348 10.745028
## 5 Jasus edwardsii 17.5 17.5 65.96031 6 19.479264 19.458764
## 6 Jasus edwardsii 17.5 21.5 58.87299 5 6.349076 7.298893
## 7 Jasus edwardsii 21.5 14 70.66670 6 19.902947 19.882001
## 8 Jasus edwardsii 21.5 17.5 79.01077 5 8.388997 9.643984
## 9 Jasus edwardsii 21.5 21.5 69.19166 6 9.502744 9.492743
## 10 Sagmariasus verreauxi 14 14 80.07882 6 14.428551 14.413366
## 11 Sagmariasus verreauxi 14 17.5 85.53435 6 19.141255 19.121111
## 12 Sagmariasus verreauxi 14 21.5 65.60422 5 17.534303 20.157420
## 13 Sagmariasus verreauxi 17.5 14 72.04296 5 6.030945 6.933170
## 14 Sagmariasus verreauxi 17.5 17.5 79.38831 5 17.778802 20.438496
## 15 Sagmariasus verreauxi 17.5 21.5 72.42871 6 13.475079 13.460898
## 16 Sagmariasus verreauxi 21.5 14 81.45283 5 10.353960 11.902904
## 17 Sagmariasus verreauxi 21.5 17.5 82.72420 6 10.575225 10.564095
## 18 Sagmariasus verreauxi 21.5 21.5 80.93551 6 11.096866 11.085187
## lower_ci upper_ci
## 1 57.98760 75.65330
## 2 55.75022 69.44261
## 3 52.67897 62.77586
## 4 56.50428 77.99433
## 5 46.50155 85.41908
## 6 51.57410 66.17188
## 7 50.78470 90.54870
## 8 69.36678 88.65475
## 9 59.69892 78.68440
## 10 65.66545 94.49218
## 11 66.41324 104.65546
## 12 45.44680 85.76164
## 13 65.10979 78.97613
## 14 58.94982 99.82681
## 15 58.96781 85.88960
## 16 69.54992 93.35573
## 17 72.16011 93.28830
## 18 69.85032 92.02070
# Experimental temperature versus acclimation - Means +/- CI
pd <- position_dodge(0.3) # move them .05 to the left and right
ggplot(value_means, aes(x=ExpTemp, y=mean, group = AccTemp, colour=AccTemp)) +
facet_wrap(~ Species) +
geom_line(position=pd) +
geom_point(position=pd) +
geom_errorbar(aes(ymin=lower_ci, ymax=upper_ci), width=.1, position=pd)

# Acclimation temperature only
ggplot(sub_trait, aes(x = Species, y = value, fill = AccTemp, alpha = .5)) +
geom_boxplot(outlier.size = 0) +
ylab(trait) +
geom_point(pch = 21, position = position_jitterdodge(seed = 100))+
geom_text(aes(label=AntTag, color = Species), alpha = .8, position = position_jitterdodge(seed = 100))

ggsave("Escapes_speed_AccTemp.pdf", width = 14, height = 8)
# Experimental temperature only
ggplot(sub_trait, aes(x = ExpTemp, y = value, fill = Species, alpha = .5)) +
geom_boxplot(outlier.size = 0) +
ylab(trait) +
geom_point(pch = 21, position = position_jitterdodge(seed = 100))+
geom_text(aes(label=AntTag, color = Species), alpha = .8, position = position_jitterdodge(seed = 100))

ggsave("Escapes_speed_ExpTemp.pdf", width = 14, height = 8)
# Species
ggplot(sub_trait, aes(x = Species, y = value, alpha = .5)) +
geom_boxplot(outlier.size = 0) +
ylab(trait) +
geom_point(pch = 21, position = position_jitterdodge(seed = 100))+
geom_text(aes(label=AntTag, color = Species), alpha = .8, position = position_jitterdodge(seed = 100))

ggsave("Escapes_speed_Species.pdf", width = 14, height = 8)
6. Summary statistics
By species
# Summary statistics
substrait_stats <- ddply(sub_trait, c("Species"), summarise,
mean = mean(value, na.rm = TRUE),
n = length(value),
sd = sd(value, na.rm = TRUE),
ci = qt(0.975, n) * sd / sqrt(n),
lower_ci = mean-ci,
upper_ci = mean+ci,
min = min(value),
max = max(value))
# Round values
substrait_stats <- substrait_stats %>% mutate_if(is.numeric, round, digit = 3)
# Display results
substrait_stats %>% kbl(digits = 1, caption = trait) %>% kable_classic()
escape_speed_cmSec
Species
|
mean
|
n
|
sd
|
ci
|
lower_ci
|
upper_ci
|
min
|
max
|
Jasus edwardsii
|
66.4
|
52
|
12.4
|
3.5
|
62.9
|
69.8
|
33.0
|
98.8
|
Sagmariasus verreauxi
|
78.1
|
50
|
14.1
|
4.0
|
74.0
|
82.1
|
47.8
|
117.4
|
By acclimation temperature and species
# Summary statistics
substrait_stats <- ddply(sub_trait, c("Species", "AccTemp"), summarise,
mean = mean(value, na.rm = TRUE),
n = length(value),
sd = sd(value, na.rm = TRUE),
ci = qt(0.975, n) * sd / sqrt(n),
lower_ci = mean-ci,
upper_ci = mean+ci,
min = min(value),
max = max(value))
# Round values
substrait_stats <- substrait_stats %>% mutate_if(is.numeric, round, digit = 3)
# Display results
substrait_stats %>% kbl(digits = 1, caption = trait) %>% kable_classic()
escape_speed_cmSec
Species
|
AccTemp
|
mean
|
n
|
sd
|
ci
|
lower_ci
|
upper_ci
|
min
|
max
|
Jasus edwardsii
|
14
|
68.2
|
18
|
13.3
|
6.6
|
61.7
|
74.8
|
49.2
|
98.8
|
Jasus edwardsii
|
17.5
|
68.6
|
17
|
14.2
|
7.3
|
61.4
|
75.9
|
33.0
|
90.9
|
Jasus edwardsii
|
21.5
|
62.1
|
17
|
8.7
|
4.4
|
57.7
|
66.6
|
53.4
|
84.7
|
Sagmariasus verreauxi
|
14
|
78.0
|
16
|
11.2
|
5.9
|
72.1
|
83.9
|
63.9
|
102.1
|
Sagmariasus verreauxi
|
17.5
|
82.7
|
17
|
15.3
|
7.8
|
74.9
|
90.6
|
58.2
|
117.4
|
Sagmariasus verreauxi
|
21.5
|
73.4
|
17
|
14.6
|
7.5
|
66.0
|
80.9
|
47.8
|
93.7
|
6. Correlations
6A EPOC correlations
# Pooled data
ggscatter(sub_trait, x = "value", y = "EPOC", fill = "Species",color = "Species",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Escape speed", ylab = "EPOC")
## `geom_smooth()` using formula 'y ~ x'

ggsave("Escapes_speed_Corelation_EPOC.pdf", width = 8, height = 6)
## `geom_smooth()` using formula 'y ~ x'
### By species
# Subset
ERL <- subset(sub_trait, Species == "Sagmariasus verreauxi")
SRL <- subset(sub_trait, Species == "Jasus edwardsii")
# Normality tests
shapiro.test(ERL$value)
##
## Shapiro-Wilk normality test
##
## data: ERL$value
## W = 0.99153, p-value = 0.9752
shapiro.test(SRL$value)
##
## Shapiro-Wilk normality test
##
## data: SRL$value
## W = 0.97349, p-value = 0.2954
# Correlation plots
ggscatter(ERL, x = "value", y = "EPOC",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Escape speed", ylab = "EPOC")+
ggtitle("Sagmariasus verreauxi")
## `geom_smooth()` using formula 'y ~ x'

ggscatter(SRL, x = "value", y = "EPOC",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Escape speed", ylab = "EPOC")+
ggtitle("Jasus edwardsii")
## `geom_smooth()` using formula 'y ~ x'

# Correlation tests
cor.test(ERL$value, ERL$EPOC, method = "pearson")
##
## Pearson's product-moment correlation
##
## data: ERL$value and ERL$EPOC
## t = 1.9561, df = 48, p-value = 0.05628
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.007171638 0.511388740
## sample estimates:
## cor
## 0.2717185
cor.test(SRL$value, SRL$EPOC, method = "pearson")
##
## Pearson's product-moment correlation
##
## data: SRL$value and SRL$EPOC
## t = 2.5465, df = 50, p-value = 0.01401
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.0726445 0.5599514
## sample estimates:
## cor
## 0.3388277
6B MMR correlations
# Pooled data
ggscatter(sub_trait, x = "value", y = "MMR", fill = "Species",color = "Species",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Escape speed", ylab = "MMR")
## `geom_smooth()` using formula 'y ~ x'

ggsave("Escapes_speed_Corelation_MMR.pdf", width = 8, height = 6)
## `geom_smooth()` using formula 'y ~ x'
### By species
# Subset
ERL <- subset(sub_trait, Species == "Sagmariasus verreauxi")
SRL <- subset(sub_trait, Species == "Jasus edwardsii")
# Normality tests
shapiro.test(ERL$value)
##
## Shapiro-Wilk normality test
##
## data: ERL$value
## W = 0.99153, p-value = 0.9752
shapiro.test(SRL$value)
##
## Shapiro-Wilk normality test
##
## data: SRL$value
## W = 0.97349, p-value = 0.2954
# Correlation plots
ggscatter(ERL, x = "value", y = "MMR",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Escape speed", ylab = "MMR")+
ggtitle("Sagmariasus verreauxi")
## `geom_smooth()` using formula 'y ~ x'

ggscatter(SRL, x = "value", y = "MMR",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Escape speed", ylab = "MMR")+
ggtitle("Jasus edwardsii")
## `geom_smooth()` using formula 'y ~ x'

# Correlation tests
cor.test(ERL$value, ERL$MMR, method = "pearson")
##
## Pearson's product-moment correlation
##
## data: ERL$value and ERL$MMR
## t = 1.1227, df = 48, p-value = 0.2672
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1239046 0.4196241
## sample estimates:
## cor
## 0.1599599
cor.test(SRL$value, SRL$MMR, method = "pearson")
##
## Pearson's product-moment correlation
##
## data: SRL$value and SRL$MMR
## t = 3.1535, df = 50, p-value = 0.002728
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1512057 0.6121570
## sample estimates:
## cor
## 0.4072992
7. Linear mixed effect modelling
7.A Find best model
Fit the most complex model with all meaningful fixed factors
Animal ID (AntTag) as random factor
names(sub_trait)
## [1] "Species" "AntTag" "AnimalID"
## [4] "AccTemp" "ExpTemp" "Sex"
## [7] "Weight_g" "Time_since_moult_d" "Acclimation_time"
## [10] "MMR" "SMR" "AerobicScope"
## [13] "recovery_time" "EPOC" "Recovery_Rate"
## [16] "Trait" "value"
# Start with the most complex model
complex_model <- lmer(value ~ ExpTemp *
AccTemp *
Species +
Sex +
Weight_g +
(1|AntTag), data = sub_trait, REML = FALSE)
# Model summary
summary(complex_model)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: value ~ ExpTemp * AccTemp * Species + Sex + Weight_g + (1 | AntTag)
## Data: sub_trait
##
## AIC BIC logLik deviance df.resid
## 820.0 877.8 -388.0 776.0 80
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.15573 -0.69706 0.06521 0.39413 2.22231
##
## Random effects:
## Groups Name Variance Std.Dev.
## AntTag (Intercept) 58.24 7.631
## Residual 78.65 8.869
## Number of obs: 102, groups: AntTag, 39
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) 65.844355 10.404383
## ExpTemp17.5 0.428856 5.120298
## ExpTemp21.5 1.716893 5.289890
## AccTemp17.5 -5.323763 6.770103
## AccTemp21.5 -9.991944 6.810285
## SpeciesSagmariasus verreauxi 12.075511 6.673776
## SexMale -6.774605 3.398479
## Weight_g 0.006302 0.008491
## ExpTemp17.5:AccTemp17.5 2.750931 7.245444
## ExpTemp21.5:AccTemp17.5 11.962417 7.587658
## ExpTemp17.5:AccTemp21.5 0.940718 7.473152
## ExpTemp21.5:AccTemp21.5 9.714262 7.363260
## ExpTemp17.5:SpeciesSagmariasus verreauxi -10.507042 7.645173
## ExpTemp21.5:SpeciesSagmariasus verreauxi -0.397893 7.619401
## AccTemp17.5:SpeciesSagmariasus verreauxi 8.328517 9.535023
## AccTemp21.5:SpeciesSagmariasus verreauxi -3.261935 9.738428
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi 7.184614 10.807793
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi -13.214850 10.749709
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi 14.315722 10.850214
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi 2.651796 10.678402
## df t value
## (Intercept) 40.706536 6.329
## ExpTemp17.5 57.454238 0.084
## ExpTemp21.5 64.655868 0.325
## AccTemp17.5 75.374323 -0.786
## AccTemp21.5 74.833607 -1.467
## SpeciesSagmariasus verreauxi 79.627783 1.809
## SexMale 35.476485 -1.993
## Weight_g 33.705049 0.742
## ExpTemp17.5:AccTemp17.5 57.608367 0.380
## ExpTemp21.5:AccTemp17.5 62.061289 1.577
## ExpTemp17.5:AccTemp21.5 58.601006 0.126
## ExpTemp21.5:AccTemp21.5 61.107761 1.319
## ExpTemp17.5:SpeciesSagmariasus verreauxi 64.685005 -1.374
## ExpTemp21.5:SpeciesSagmariasus verreauxi 63.341635 -0.052
## AccTemp17.5:SpeciesSagmariasus verreauxi 77.745570 0.873
## AccTemp21.5:SpeciesSagmariasus verreauxi 77.496155 -0.335
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi 64.079607 0.665
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi 62.739691 -1.229
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi 62.105356 1.319
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi 61.046947 0.248
## Pr(>|t|)
## (Intercept) 1.52e-07 ***
## ExpTemp17.5 0.9335
## ExpTemp21.5 0.7466
## AccTemp17.5 0.4341
## AccTemp21.5 0.1465
## SpeciesSagmariasus verreauxi 0.0742 .
## SexMale 0.0539 .
## Weight_g 0.4632
## ExpTemp17.5:AccTemp17.5 0.7056
## ExpTemp21.5:AccTemp17.5 0.1200
## ExpTemp17.5:AccTemp21.5 0.9003
## ExpTemp21.5:AccTemp21.5 0.1920
## ExpTemp17.5:SpeciesSagmariasus verreauxi 0.1741
## ExpTemp21.5:SpeciesSagmariasus verreauxi 0.9585
## AccTemp17.5:SpeciesSagmariasus verreauxi 0.3851
## AccTemp21.5:SpeciesSagmariasus verreauxi 0.7386
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi 0.5086
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi 0.2235
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi 0.1919
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi 0.8047
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 20 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
# Anova
anova(complex_model)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ExpTemp 1061.32 530.66 2 61.666 6.7469 0.0022398 **
## AccTemp 361.04 180.52 2 32.543 2.2951 0.1168285
## Species 1050.39 1050.39 1 32.487 13.3549 0.0009007 ***
## Sex 312.54 312.54 1 35.476 3.9737 0.0539487 .
## Weight_g 43.32 43.32 1 33.705 0.5507 0.4631566
## ExpTemp:AccTemp 390.38 97.60 4 61.489 1.2408 0.3030680
## ExpTemp:Species 73.78 36.89 2 61.787 0.4690 0.6278123
## AccTemp:Species 57.40 28.70 2 32.715 0.3649 0.6970608
## ExpTemp:AccTemp:Species 399.97 99.99 4 61.681 1.2713 0.2909532
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
report(anova(complex_model))
## Type 3 ANOVAs only give sensible and informative results when covariates are
## mean-centered and factors are coded with orthogonal contrasts (such as those
## produced by 'contr.sum', 'contr.poly', or 'contr.helmert', but *not* by the
## default 'contr.treatment').
## The ANOVA suggests that:
##
## - The main effect of ExpTemp is statistically significant and large (F(2) = 6.75, p = 0.002; Eta2 (partial) = 0.18, 90% CI [0.05, 0.31])
## - The main effect of AccTemp is statistically not significant and medium (F(2) = 2.30, p = 0.117; Eta2 (partial) = 0.12, 90% CI [0.00, 0.29])
## - The main effect of Species is statistically significant and large (F(1) = 13.35, p < .001; Eta2 (partial) = 0.29, 90% CI [0.09, 0.48])
## - The main effect of Sex is statistically not significant and medium (F(1) = 3.97, p = 0.054; Eta2 (partial) = 0.10, 90% CI [0.00, 0.28])
## - The main effect of Weight_g is statistically not significant and small (F(1) = 0.55, p = 0.463; Eta2 (partial) = 0.02, 90% CI [0.00, 0.14])
## - The interaction between ExpTemp and AccTemp is statistically not significant and medium (F(4) = 1.24, p = 0.303; Eta2 (partial) = 0.07, 90% CI [0.00, 0.15])
## - The interaction between ExpTemp and Species is statistically not significant and small (F(2) = 0.47, p = 0.628; Eta2 (partial) = 0.01, 90% CI [0.00, 0.08])
## - The interaction between AccTemp and Species is statistically not significant and small (F(2) = 0.36, p = 0.697; Eta2 (partial) = 0.02, 90% CI [0.00, 0.12])
## - The interaction between ExpTemp, AccTemp and Species is statistically not significant and medium (F(4) = 1.27, p = 0.291; Eta2 (partial) = 0.08, 90% CI [0.00, 0.16])
##
## Effect sizes were labelled following Field's (2013) recommendations.
7.B Find the best model by backward elimination of non-significant effects of linear mixed effects model using the step function of the LmerTest package
(the model with the lowest AIC and significant p value is the best)
# Backward elimination
step_result <- step(complex_model)
step_result
## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 22 -388.01 820.03
## (1 | AntTag) 0 21 -393.54 829.09 11.06 1 0.000882 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## Weight_g 1 43.32 43.32 1 33.705 0.5507
## ExpTemp:AccTemp:Species 2 393.81 98.45 4 61.875 1.2530
## AccTemp:Species 3 76.34 38.17 2 34.010 0.4480
## ExpTemp:Species 4 80.51 40.25 2 62.583 0.4673
## ExpTemp:AccTemp 5 360.01 90.00 4 62.216 1.0397
## AccTemp 6 449.37 224.68 2 34.036 2.3853
## Sex 7 305.63 305.63 1 37.525 3.2737
## ExpTemp 0 1014.60 507.30 2 63.723 5.4097
## Species 0 1176.09 1176.09 1 34.893 12.5416
## Pr(>F)
## Weight_g 0.463157
## ExpTemp:AccTemp:Species 0.298134
## AccTemp:Species 0.642605
## ExpTemp:Species 0.628875
## ExpTemp:AccTemp 0.394044
## AccTemp 0.107280
## Sex 0.078414 .
## ExpTemp 0.006761 **
## Species 0.001152 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## value ~ ExpTemp + Species + (1 | AntTag)
# Extract the model that step found:
final_model <- get_model(step_result)
–> Species had strong effect, and experimental temperature
–> no acclimation effects
7.C Re-run model fit using REML estimation
Best model: ExpTemp + Species + (1 | AntTag)
# Remodel using Maximum likelihood estimation
final_model <- lmer(value ~ ExpTemp +
Species +
(1|AntTag), data = sub_trait, REML = TRUE)
# Model summary
summary(final_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: value ~ ExpTemp + Species + (1 | AntTag)
## Data: sub_trait
##
## REML criterion at convergence: 780.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.78444 -0.66668 -0.05271 0.48793 2.69398
##
## Random effects:
## Groups Name Variance Std.Dev.
## AntTag (Intercept) 76.93 8.771
## Residual 96.60 9.828
## Number of obs: 102, groups: AntTag, 39
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 64.33248 2.80295 53.38791 22.952 < 2e-16 ***
## ExpTemp17.5 0.07923 2.43647 62.25737 0.033 0.97416
## ExpTemp21.5 6.83380 2.41174 61.85325 2.834 0.00621 **
## SpeciesSagmariasus verreauxi 11.88535 3.45276 33.15498 3.442 0.00158 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ET17.5 ET21.5
## ExpTemp17.5 -0.420
## ExpTemp21.5 -0.423 0.491
## SpcsSgmrssv -0.620 0.004 0.003
# Anova to test main effects of final model representing the combined significance for all coefficients
anova(final_model)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ExpTemp 1010.4 505.17 2 61.891 5.2297 0.00797 **
## Species 1144.6 1144.61 1 33.155 11.8493 0.00158 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
report(anova(final_model))
## The ANOVA suggests that:
##
## - The main effect of ExpTemp is statistically significant and large (F(2) = 5.23, p = 0.008; Eta2 (partial) = 0.14, 90% CI [0.02, 0.27])
## - The main effect of Species is statistically significant and large (F(1) = 11.85, p = 0.002; Eta2 (partial) = 0.26, 90% CI [0.07, 0.45])
##
## Effect sizes were labelled following Field's (2013) recommendations.
# Table output from sjPlot
tab_model(final_model)
|
value
|
Predictors
|
Estimates
|
CI
|
p
|
(Intercept)
|
64.33
|
58.84 – 69.83
|
<0.001
|
ExpTemp [17.5]
|
0.08
|
-4.70 – 4.85
|
0.974
|
ExpTemp [21.5]
|
6.83
|
2.11 – 11.56
|
0.005
|
Species [Sagmariasus verreauxi]
|
11.89
|
5.12 – 18.65
|
0.001
|
Random Effects
|
σ2
|
96.60
|
τ00 AntTag
|
76.93
|
ICC
|
0.44
|
N AntTag
|
39
|
Observations
|
102
|
Marginal R2 / Conditional R2
|
0.212 / 0.561
|
# Textual report of the model
report(final_model)
## We fitted a linear mixed model (estimated using REML and nloptwrap optimizer) to predict value with ExpTemp and Species (formula: value ~ ExpTemp + Species). The model included AntTag as random effect (formula: ~1 | AntTag). The model's total explanatory power is substantial (conditional R2 = 0.56) and the part related to the fixed effects alone (marginal R2) is of 0.21. The model's intercept, corresponding to ExpTemp = 14 and Species = Jasus edwardsii, is at 64.33 (95% CI [58.84, 69.83], t(96) = 22.95, p < .001). Within this model:
##
## - The effect of ExpTemp [17.5] is statistically non-significant and positive (beta = 0.08, 95% CI [-4.70, 4.85], t(96) = 0.03, p = 0.974; Std. beta = 5.48e-03, 95% CI [-0.32, 0.34])
## - The effect of ExpTemp [21.5] is statistically significant and positive (beta = 6.83, 95% CI [2.11, 11.56], t(96) = 2.83, p = 0.005; Std. beta = 0.47, 95% CI [0.15, 0.80])
## - The effect of Species [Sagmariasus verreauxi] is statistically significant and positive (beta = 11.89, 95% CI [5.12, 18.65], t(96) = 3.44, p < .001; Std. beta = 0.82, 95% CI [0.35, 1.29])
##
## Standardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed using the Wald approximation.
Summary table of main effects
library(effectsize)
# Calculate anova test stats for model
aov_lme <-anova(final_model)
# Caculate partial effects size
eta = eta_squared(aov_lme)
# Create output stats table
lme_table <- round(aov_lme[, c(1,2,4)], 1)
# Attach F statistics
lme_table$F_value <- paste(round(aov_lme[,5],2), " (", aov_lme[,3], ")", sep = "")
# Attach p value
lme_table$p_value <- format(aov_lme[, 6], digits = 4)
# Create string with eta and 90% eta range and join with table
lme_table$Partial_Eta <- paste(round(eta$Eta2_partial, 2), " [", round(eta$CI_low, 2), "-", round(eta$CI_high, 2), "]", sep = "")
# Save table as csv
write.csv(lme_table, paste("Model_summary_", trait, ".csv", sep = ""))
# Print table
kable(lme_table,
format = "html",
row.names = TRUE,
caption = "Model summary for main effects",
col.names = c("Sum of sqares", "Mean squares", "Den df", "F(df)", "p value", "Partical squared eta [CI 90%]")) %>%
kable_classic(full_width = F, html_font = "Cambria")
Model summary for main effects
|
Sum of sqares
|
Mean squares
|
Den df
|
F(df)
|
p value
|
Partical squared eta [CI 90%]
|
ExpTemp
|
1010.3
|
505.2
|
61.9
|
5.23 (2)
|
0.00797
|
0.14 [0.02-0.27]
|
Species
|
1144.6
|
1144.6
|
33.2
|
11.85 (1)
|
0.00158
|
0.26 [0.07-0.45]
|
8. Model diagnostics
# Use of DHARMa package see https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html
# Set see for reproducibility of residual simulation
set.seed(1000)
# Simulate residuals
# (scaled residual value of 0.5 means that half of the simulated data are higher than the observed value)
simulationOutput <- simulateResiduals(fittedModel = final_model, plot = F)
# Plot scaled residuals
# Plot 1 = QQ plot with test for KS, overdispersion and outliers
# Plot 2 = Residuals plot
plot(simulationOutput)

# Plot against predictors to check for further deviations
plotResiduals(simulationOutput, sub_trait$ExpTemp)

–> Model diagnostics do not show apparent violations of model assumptions
9. Multiple comparisons across treatment levels
9.A Multiple comparisons between species for each experimental temperature
# Perform multiple comparisions using the emmeans package
# Instructions here https://cran.r-project.org/web/packages/emmeans/vignettes/interactions.html
contrasts <- emmeans(final_model, pairwise ~ Species | ExpTemp, type = "response", adjust = "bonferroni")
# Display interaction plots
emmip(final_model, Species ~ ExpTemp)

# Convert results to dataframe (https://aosmith.rbind.io/2019/03/25/getting-started-with-emmeans/)
contrasts_df <- contrasts$contrasts %>%
summary(infer = TRUE) %>%
as.data.frame()
# Plot comparisons
### blue bars are confidence intervals for the EMMs, and the red arrows are for the comparisons among them. If an arrow ### from one mean overlaps an arrow from another group, the difference is not “significant,” based on the adjust setting ### (which defaults to "tukey") and the value of alpha (which defaults to 0.05)
### https://cran.r-project.org/web/packages/emmeans/vignettes/comparisons.html
plot(contrasts, comparisons = TRUE)

# Save contrasts in table
contrasts_expTemp <- contrasts_df[, c(1:3,8,9)] %>% mutate(p.value = round(p.value, 3))
# Rename column headers
names(contrasts_expTemp)[1:5] <- c("Species", "ExpTemp", "Difference", "t-ratio", "p")
# Add significance descriptor
contrasts_expTemp$contrasts_sign <- ifelse(contrasts_expTemp$p < 0.05, "Yes", "No")
# Add trait columns
contrasts_expTemp$Trait <- trait
# Save contrasts results
write.csv(contrasts_expTemp,paste("Contrasts_AcclimationTemperature_", trait, ".csv", sep = ""))
# Display as HTML table colouring significant p values red
contrasts_expTemp %>% mutate(p = cell_spec(p, "html", color = ifelse(p < 0.05, "red", "grey"))) %>%
kable(escape = FALSE) %>%
kable_styling(bootstrap_options = c('hover'))
Species
|
ExpTemp
|
Difference
|
t-ratio
|
p
|
contrasts_sign
|
Trait
|
Jasus edwardsii - Sagmariasus verreauxi
|
14
|
-11.88535
|
-3.439063
|
0.001
|
Yes
|
escape_speed_cmSec
|
Jasus edwardsii - Sagmariasus verreauxi
|
17.5
|
-11.88535
|
-3.439063
|
0.001
|
Yes
|
escape_speed_cmSec
|
Jasus edwardsii - Sagmariasus verreauxi
|
21.5
|
-11.88535
|
-3.439063
|
0.001
|
Yes
|
escape_speed_cmSec
|
9.B Multiple comparisons between experimental temperature for each species
# Perform multiple comparisions using the emmeans package
# Instructions here https://cran.r-project.org/web/packages/emmeans/vignettes/interactions.html
contrasts <- emmeans(final_model, pairwise ~ ExpTemp | Species, type = "response", adjust = "bonferroni")
# Display interaction plots
emmip(final_model, ExpTemp ~ Species)

# Convert results to dataframe (https://aosmith.rbind.io/2019/03/25/getting-started-with-emmeans/)
contrasts_df <- contrasts$contrasts %>%
summary(infer = TRUE) %>%
as.data.frame()
# Plot comparisons
### blue bars are confidence intervals for the EMMs, and the red arrows are for the comparisons among them. If an arrow ### from one mean overlaps an arrow from another group, the difference is not “significant,” based on the adjust setting ### (which defaults to "tukey") and the value of alpha (which defaults to 0.05)
### https://cran.r-project.org/web/packages/emmeans/vignettes/comparisons.html
plot(contrasts, comparisons = TRUE)

# Save contrasts in table
contrasts_expTemp <- contrasts_df[, c(1:3,8,9)] %>% mutate(p.value = round(p.value, 3))
# Rename column headers
names(contrasts_expTemp)[1:5] <- c("ExpTemp", "Species", "Difference", "t-ratio", "p")
# Add significance descriptor
contrasts_expTemp$contrasts_sign <- ifelse(contrasts_expTemp$p < 0.05, "Yes", "No")
# Add trait columns
contrasts_expTemp$Trait <- trait
# Save contrasts results
write.csv(contrasts_expTemp,paste("Contrasts_AcclimationTemperature_", trait, ".csv", sep = ""))
# Display as HTML table colouring significant p values red
contrasts_expTemp %>% mutate(p = cell_spec(p, "html", color = ifelse(p < 0.05, "red", "grey"))) %>%
kable(escape = FALSE) %>%
kable_styling(bootstrap_options = c('hover'))
ExpTemp
|
Species
|
Difference
|
t-ratio
|
p
|
contrasts_sign
|
Trait
|
14 - 17.5
|
Jasus edwardsii
|
-0.0792253
|
-0.0324524
|
1
|
No
|
escape_speed_cmSec
|
14 - 21.5
|
Jasus edwardsii
|
-6.8337956
|
-2.8286253
|
0.019
|
Yes
|
escape_speed_cmSec
|
17.5 - 21.5
|
Jasus edwardsii
|
-6.7545703
|
-2.7564300
|
0.023
|
Yes
|
escape_speed_cmSec
|
14 - 17.5
|
Sagmariasus verreauxi
|
-0.0792253
|
-0.0324524
|
1
|
No
|
escape_speed_cmSec
|
14 - 21.5
|
Sagmariasus verreauxi
|
-6.8337956
|
-2.8286253
|
0.019
|
Yes
|
escape_speed_cmSec
|
17.5 - 21.5
|
Sagmariasus verreauxi
|
-6.7545703
|
-2.7564300
|
0.023
|
Yes
|
escape_speed_cmSec
|