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 <- "recovery_time"
# File path to folder
file_path <- paste("C:/Users/oelle/Documents/Meine Dokumente/Research/Projects/Project_RockLobster/Temperature acclimation/Aerobic scope/Statistics/LinearMixedEffectModelling/", trait, "/", sep = "")
# Load tidy raw data in long format
mo2_data_long <- read.csv("C:/Users/oelle/Documents/Meine Dokumente/Research/Projects/Project_RockLobster/Temperature acclimation/Aerobic scope/Statistics/MO2_data_Tidy_RawData_long.csv")
3. Process data
# Shorten column names
names(mo2_data_long)[c(5,6)] <- c("AccTemp", "ExpTemp")
# Define fixed effects as factors
mo2_data_long$AntTag <- as.factor(mo2_data_long$AntTag)
mo2_data_long$Species <- as.factor(mo2_data_long$Species)
mo2_data_long$AccTemp <- as.factor(mo2_data_long$AccTemp)
mo2_data_long$ExpTemp <- as.factor(mo2_data_long$ExpTemp)
mo2_data_long$Sex <- as.factor(mo2_data_long$Sex)
# Subset data by trait
sub_trait <- subset(mo2_data_long, Trait == trait)
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: 1.64 max: 24.27
## median: 10.27
## mean: 10.23565
## estimated sd: 4.662695
## estimated skewness: 0.1954576
## estimated kurtosis: 2.926359
# Fit various distributions
nd <- fitdist(sub_trait$value, "norm")
## $start.arg
## $start.arg$mean
## [1] 10.23565
##
## $start.arg$sd
## [1] 4.641058
##
##
## $fix.arg
## NULL
gd <- fitdist(sub_trait$value, "gamma")
## $start.arg
## $start.arg$shape
## [1] 4.864035
##
## $start.arg$rate
## [1] 0.4752054
##
##
## $fix.arg
## NULL
wd <- fitdist(sub_trait$value, "weibull")
## $start.arg
## $start.arg$shape
## [1] 2.036858
##
## $start.arg$scale
## [1] 11.79931
##
##
## $fix.arg
## NULL
lnd <- fitdist(sub_trait$value, "lnorm")
## $start.arg
## $start.arg$meanlog
## [1] 2.187216
##
## $start.arg$sdlog
## [1] 0.5864089
##
##
## $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.06203552 0.1200286 0.07758804 0.1523719
## Cramer-von Mises statistic 0.06039320 0.3115573 0.11945781 0.5722539
## Anderson-Darling statistic 0.41020508 1.8311508 0.74995768 3.4632159
##
## Goodness-of-fit criteria
## Normal gamma Weibull LogNormal
## Akaike's Information Criterion 642.0383 649.2826 639.1851 667.6421
## Bayesian Information Criterion 647.4025 654.6468 644.5494 673.0063
–> Anderson-Darling statistics gives good measure for fit for both middle and tail of distribution
–> Data fit best to a normal distribution
5. Plot data
# Gplot
ggplot(sub_trait, aes(x = ExpTemp, y = value, fill = as.factor(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 = as.factor(AccTemp)), alpha = .8, position = position_jitterdodge(seed = 100))

# 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")

6. Summary statistics
By species
# Summary statistics
substrait_stats <- ddply(sub_trait, c("Species", "ExpTemp"), 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()
recovery_time
Species
|
ExpTemp
|
mean
|
n
|
sd
|
ci
|
lower_ci
|
upper_ci
|
min
|
max
|
Jasus edwardsii
|
14
|
8.6
|
18
|
3.8
|
1.9
|
6.7
|
10.5
|
1.6
|
14.0
|
Jasus edwardsii
|
17.5
|
10.8
|
18
|
5.2
|
2.6
|
8.2
|
13.3
|
1.9
|
20.8
|
Jasus edwardsii
|
21.5
|
11.2
|
18
|
4.6
|
2.3
|
9.0
|
13.5
|
2.1
|
18.7
|
Sagmariasus verreauxi
|
14
|
8.7
|
18
|
4.1
|
2.0
|
6.7
|
10.7
|
1.6
|
15.2
|
Sagmariasus verreauxi
|
17.5
|
9.8
|
18
|
3.8
|
1.9
|
7.9
|
11.7
|
1.9
|
17.7
|
Sagmariasus verreauxi
|
21.5
|
12.3
|
18
|
5.6
|
2.8
|
9.5
|
15.1
|
2.3
|
24.3
|
# Summary statistics for acclimation temperature
substrait_stats_acc <- 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_acc <- substrait_stats_acc %>% mutate_if(is.numeric, round, digit = 3)
# Display results
substrait_stats_acc %>% kbl(digits = 1, caption = trait) %>% kable_classic()
recovery_time
Species
|
AccTemp
|
mean
|
n
|
sd
|
ci
|
lower_ci
|
upper_ci
|
min
|
max
|
Jasus edwardsii
|
14
|
12.0
|
18
|
3.9
|
1.9
|
10.1
|
14.0
|
5.8
|
18.7
|
Jasus edwardsii
|
17.5
|
11.2
|
18
|
4.5
|
2.2
|
9.0
|
13.5
|
1.9
|
20.8
|
Jasus edwardsii
|
21.5
|
7.3
|
18
|
4.2
|
2.1
|
5.2
|
9.4
|
1.6
|
16.6
|
Sagmariasus verreauxi
|
14
|
12.4
|
18
|
5.4
|
2.7
|
9.7
|
15.1
|
2.6
|
24.3
|
Sagmariasus verreauxi
|
17.5
|
11.1
|
18
|
3.1
|
1.5
|
9.6
|
12.6
|
7.0
|
17.7
|
Sagmariasus verreauxi
|
21.5
|
7.4
|
18
|
4.1
|
2.0
|
5.4
|
9.4
|
1.6
|
14.0
|
–> Warm acclimated J. edwardsii recover 4.7 hours, and S. verreauxi 5 h faster than cold acclimated lobsters
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
# 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
## 621.4 680.4 -288.7 577.4 86
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2847 -0.6662 -0.1295 0.6145 2.9688
##
## Random effects:
## Groups Name Variance Std.Dev.
## AntTag (Intercept) 2.929 1.711
## Residual 9.924 3.150
## Number of obs: 108, groups: AntTag, 39
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) 8.745418 2.977257
## ExpTemp17.5 3.963333 1.818755
## ExpTemp21.5 3.740013 1.852708
## AccTemp17.5 1.449575 2.078070
## AccTemp21.5 -3.905287 2.094626
## SpeciesSagmariasus verreauxi -1.180706 2.061866
## SexMale -1.844171 0.917818
## Weight_g 0.001738 0.002355
## ExpTemp17.5:AccTemp17.5 -4.990775 2.573028
## ExpTemp21.5:AccTemp17.5 -1.525013 2.596227
## ExpTemp17.5:AccTemp21.5 -0.355791 2.572138
## ExpTemp21.5:AccTemp21.5 -1.494137 2.596571
## ExpTemp17.5:SpeciesSagmariasus verreauxi -0.490364 2.595515
## ExpTemp21.5:SpeciesSagmariasus verreauxi 4.906290 2.620143
## AccTemp17.5:SpeciesSagmariasus verreauxi 1.384899 2.936110
## AccTemp21.5:SpeciesSagmariasus verreauxi 2.468335 2.947474
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi 1.406433 3.678142
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi -7.170423 3.688560
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi -3.435465 3.655011
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi -5.292119 3.673618
## df t value
## (Intercept) 52.070994 2.937
## ExpTemp17.5 70.764971 2.179
## ExpTemp21.5 78.893798 2.019
## AccTemp17.5 99.590701 0.698
## AccTemp21.5 98.819832 -1.864
## SpeciesSagmariasus verreauxi 101.941888 -0.573
## SexMale 40.501138 -2.009
## Weight_g 39.817102 0.738
## ExpTemp17.5:AccTemp17.5 70.862814 -1.940
## ExpTemp21.5:AccTemp17.5 74.863750 -0.587
## ExpTemp17.5:AccTemp21.5 70.768132 -0.138
## ExpTemp21.5:AccTemp21.5 74.880424 -0.575
## ExpTemp17.5:SpeciesSagmariasus verreauxi 74.908363 -0.189
## ExpTemp21.5:SpeciesSagmariasus verreauxi 79.027780 1.873
## AccTemp17.5:SpeciesSagmariasus verreauxi 100.976679 0.472
## AccTemp21.5:SpeciesSagmariasus verreauxi 99.605891 0.837
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi 75.456961 0.382
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi 76.972520 -1.944
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi 72.927934 -0.940
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi 75.049184 -1.441
## Pr(>|t|)
## (Intercept) 0.00492 **
## ExpTemp17.5 0.03265 *
## ExpTemp21.5 0.04692 *
## AccTemp17.5 0.48708
## AccTemp21.5 0.06523 .
## SpeciesSagmariasus verreauxi 0.56815
## SexMale 0.05120 .
## Weight_g 0.46480
## ExpTemp17.5:AccTemp17.5 0.05640 .
## ExpTemp21.5:AccTemp17.5 0.55871
## ExpTemp17.5:AccTemp21.5 0.89038
## ExpTemp21.5:AccTemp21.5 0.56673
## ExpTemp17.5:SpeciesSagmariasus verreauxi 0.85066
## ExpTemp21.5:SpeciesSagmariasus verreauxi 0.06483 .
## AccTemp17.5:SpeciesSagmariasus verreauxi 0.63817
## AccTemp21.5:SpeciesSagmariasus verreauxi 0.40435
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi 0.70326
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi 0.05555 .
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi 0.35035
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi 0.15386
## ---
## 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 to test main effects of final model representing the combined significance for all coefficients
anova(complex_model)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ExpTemp 170.908 85.454 2 73.471 8.6112 0.0004367 ***
## AccTemp 237.132 118.566 2 38.705 11.9479 9.098e-05 ***
## Species 0.017 0.017 1 38.853 0.0017 0.9672298
## Sex 40.064 40.064 1 40.501 4.0373 0.0512016 .
## Weight_g 5.406 5.406 1 39.817 0.5448 0.4647960
## ExpTemp:AccTemp 101.493 25.373 4 73.466 2.5569 0.0457301 *
## ExpTemp:Species 16.674 8.337 2 73.696 0.8401 0.4357557
## AccTemp:Species 0.792 0.396 2 39.012 0.0399 0.9609073
## ExpTemp:AccTemp:Species 79.778 19.944 4 73.639 2.0098 0.1019061
## ---
## 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) = 8.61, p < .001; Eta2 (partial) = 0.19, 90% CI [0.06, 0.31])
## - The main effect of AccTemp is statistically significant and large (F(2) = 11.95, p < .001; Eta2 (partial) = 0.38, 90% CI [0.17, 0.53])
## - The main effect of Species is statistically not significant and very small (F(1) = 1.71e-03, p = 0.967; Eta2 (partial) = 4.40e-05, 90% CI [0.00, 0.00])
## - The main effect of Sex is statistically not significant and medium (F(1) = 4.04, p = 0.051; Eta2 (partial) = 0.09, 90% CI [0.00, 0.25])
## - The main effect of Weight_g is statistically not significant and small (F(1) = 0.54, p = 0.465; Eta2 (partial) = 0.01, 90% CI [0.00, 0.12])
## - The interaction between ExpTemp and AccTemp is statistically significant and medium (F(4) = 2.56, p = 0.046; Eta2 (partial) = 0.12, 90% CI [1.33e-03, 0.22])
## - The interaction between ExpTemp and Species is statistically not significant and small (F(2) = 0.84, p = 0.436; Eta2 (partial) = 0.02, 90% CI [0.00, 0.09])
## - The interaction between AccTemp and Species is statistically not significant and very small (F(2) = 0.04, p = 0.961; Eta2 (partial) = 2.04e-03, 90% CI [0.00, 0.00])
## - The interaction between ExpTemp, AccTemp and Species is statistically not significant and medium (F(4) = 2.01, p = 0.102; Eta2 (partial) = 0.10, 90% CI [0.00, 0.18])
##
## 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, reduce.random=FALSE)
step_result
## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 22 -288.72 621.44
## (1 | AntTag) 0 21 -291.24 624.47 5.0367 1 0.02482 *
## ---
## 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 5.406 5.406 1 39.817 0.5448
## ExpTemp:AccTemp:Species 2 78.769 19.692 4 73.715 1.9982
## AccTemp:Species 3 0.376 0.188 2 38.700 0.0173
## ExpTemp:Species 4 17.544 8.772 2 73.389 0.8092
## Species 5 0.082 0.082 1 38.859 0.0074
## Sex 6 35.906 35.906 1 39.947 3.2439
## ExpTemp:AccTemp 7 101.772 25.443 4 73.191 2.3017
## ExpTemp 0 170.143 85.072 2 72.877 6.8585
## AccTemp 0 272.327 136.163 2 38.326 10.9775
## Pr(>F)
## Weight_g 0.4647960
## ExpTemp:AccTemp:Species 0.1036202
## AccTemp:Species 0.9828439
## ExpTemp:Species 0.4491553
## Species 0.9320584
## Sex 0.0792432 .
## ExpTemp:AccTemp 0.0666062 .
## ExpTemp 0.0018658 **
## AccTemp 0.0001702 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## value ~ ExpTemp + AccTemp + (1 | AntTag)
# Extract the model that step found:
final_model <- get_model(step_result)
–> Experimental and acclimation temperature + species + acclimation time have effect on data
7.C Re-run model fit using REML estimation
Best model: ExpTemp + AccTemp + (1 | AntTag)
# Remodel using Maximum likelihood estimation
final_model <- lmer(value ~ ExpTemp +
AccTemp +
(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 + AccTemp + (1 | AntTag)
## Data: sub_trait
##
## REML criterion at convergence: 592.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.51225 -0.60882 -0.03267 0.48846 2.43002
##
## Random effects:
## Groups Name Variance Std.Dev.
## AntTag (Intercept) 3.681 1.919
## Residual 12.766 3.573
## Number of obs: 108, groups: AntTag, 39
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 10.6813 0.9336 65.0307 11.441 < 2e-16 ***
## ExpTemp17.5 1.5766 0.8470 70.8194 1.861 0.066827 .
## ExpTemp21.5 3.0949 0.8494 72.0733 3.644 0.000503 ***
## AccTemp17.5 -1.0660 1.1311 36.6461 -0.942 0.352137
## AccTemp21.5 -4.8821 1.1376 35.3434 -4.292 0.000131 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ET17.5 ET21.5 AT17.5
## ExpTemp17.5 -0.452
## ExpTemp21.5 -0.456 0.504
## AccTemp17.5 -0.600 0.000 0.003
## AccTemp21.5 -0.595 -0.003 0.000 0.492
# 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 169.49 84.746 2 70.787 6.6386 0.0022782 **
## AccTemp 258.09 129.045 2 35.379 10.1088 0.0003369 ***
## ---
## 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) = 6.64, p = 0.002; Eta2 (partial) = 0.16, 90% CI [0.04, 0.28])
## - The main effect of AccTemp is statistically significant and large (F(2) = 10.11, p < .001; Eta2 (partial) = 0.36, 90% CI [0.14, 0.52])
##
## Effect sizes were labelled following Field's (2013) recommendations.
# Table output from sjPlot
tab_model(final_model)
|
value
|
Predictors
|
Estimates
|
CI
|
p
|
(Intercept)
|
10.68
|
8.85 – 12.51
|
<0.001
|
ExpTemp [17.5]
|
1.58
|
-0.08 – 3.24
|
0.063
|
ExpTemp [21.5]
|
3.09
|
1.43 – 4.76
|
<0.001
|
AccTemp [17.5]
|
-1.07
|
-3.28 – 1.15
|
0.346
|
AccTemp [21.5]
|
-4.88
|
-7.11 – -2.65
|
<0.001
|
Random Effects
|
σ2
|
12.77
|
τ00 AntTag
|
3.68
|
ICC
|
0.22
|
N AntTag
|
39
|
Observations
|
108
|
Marginal R2 / Conditional R2
|
0.269 / 0.432
|
# 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 AccTemp (formula: value ~ ExpTemp + AccTemp). The model included AntTag as random effect (formula: ~1 | AntTag). The model's total explanatory power is substantial (conditional R2 = 0.43) and the part related to the fixed effects alone (marginal R2) is of 0.27. The model's intercept, corresponding to ExpTemp = 14 and AccTemp = 14, is at 10.68 (95% CI [8.85, 12.51], t(101) = 11.44, p < .001). Within this model:
##
## - The effect of ExpTemp [17.5] is statistically non-significant and positive (beta = 1.58, 95% CI [-0.08, 3.24], t(101) = 1.86, p = 0.063; Std. beta = 0.34, 95% CI [-0.02, 0.69])
## - The effect of ExpTemp [21.5] is statistically significant and positive (beta = 3.09, 95% CI [1.43, 4.76], t(101) = 3.64, p < .001; Std. beta = 0.66, 95% CI [0.31, 1.02])
## - The effect of AccTemp [17.5] is statistically non-significant and negative (beta = -1.07, 95% CI [-3.28, 1.15], t(101) = -0.94, p = 0.346; Std. beta = -0.23, 95% CI [-0.70, 0.25])
## - The effect of AccTemp [21.5] is statistically significant and negative (beta = -4.88, 95% CI [-7.11, -2.65], t(101) = -4.29, p < .001; Std. beta = -1.05, 95% CI [-1.53, -0.57])
##
## 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, "Model_summary.csv")
# 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
|
169.5
|
84.7
|
70.8
|
6.64 (2)
|
0.0022782
|
0.16 [0.04-0.28]
|
AccTemp
|
258.1
|
129.0
|
35.4
|
10.11 (2)
|
0.0003369
|
0.36 [0.14-0.52]
|
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$AccTemp)

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 acclimation temperature 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 ~ AccTemp | ExpTemp, type = "response", adjust = "bonferroni")
# Display interaction plots
emmip(final_model, AccTemp ~ 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("AccTemp", "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'))
AccTemp
|
ExpTemp
|
Difference
|
t-ratio
|
p
|
contrasts_sign
|
Trait
|
14 - 17.5
|
14
|
1.066036
|
0.9409613
|
1
|
No
|
recovery_time
|
14 - 21.5
|
14
|
4.882122
|
4.2874294
|
0
|
Yes
|
recovery_time
|
17.5 - 21.5
|
14
|
3.816086
|
3.3347182
|
0.006
|
Yes
|
recovery_time
|
14 - 17.5
|
17.5
|
1.066036
|
0.9409613
|
1
|
No
|
recovery_time
|
14 - 21.5
|
17.5
|
4.882122
|
4.2874294
|
0
|
Yes
|
recovery_time
|
17.5 - 21.5
|
17.5
|
3.816086
|
3.3347182
|
0.006
|
Yes
|
recovery_time
|
14 - 17.5
|
21.5
|
1.066036
|
0.9409613
|
1
|
No
|
recovery_time
|
14 - 21.5
|
21.5
|
4.882122
|
4.2874294
|
0
|
Yes
|
recovery_time
|
17.5 - 21.5
|
21.5
|
3.816086
|
3.3347182
|
0.006
|
Yes
|
recovery_time
|
9.B Multiple comparisons between experimental temperatures for each acclimation 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 ~ ExpTemp | AccTemp, type = "response", adjust = "bonferroni")
# Display interaction plots
emmip(final_model, ExpTemp ~ AccTemp)

# 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", "AccTemp", "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_ExperimentalTemperature_", 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
|
AccTemp
|
Difference
|
t-ratio
|
p
|
contrasts_sign
|
Trait
|
14 - 17.5
|
14
|
-1.576612
|
-1.859695
|
0.201
|
No
|
recovery_time
|
14 - 21.5
|
14
|
-3.094877
|
-3.638514
|
0.002
|
Yes
|
recovery_time
|
17.5 - 21.5
|
14
|
-1.518265
|
-1.796872
|
0.23
|
No
|
recovery_time
|
14 - 17.5
|
17.5
|
-1.576612
|
-1.859695
|
0.201
|
No
|
recovery_time
|
14 - 21.5
|
17.5
|
-3.094877
|
-3.638514
|
0.002
|
Yes
|
recovery_time
|
17.5 - 21.5
|
17.5
|
-1.518265
|
-1.796872
|
0.23
|
No
|
recovery_time
|
14 - 17.5
|
21.5
|
-1.576612
|
-1.859695
|
0.201
|
No
|
recovery_time
|
14 - 21.5
|
21.5
|
-3.094877
|
-3.638514
|
0.002
|
Yes
|
recovery_time
|
17.5 - 21.5
|
21.5
|
-1.518265
|
-1.796872
|
0.23
|
No
|
recovery_time
|