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 <- "EPOC"

# Set file path
  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:  20.15   max:  687.15 
## median:  306.425 
## mean:  311.551 
## estimated sd:  171.0271 
## estimated skewness:  0.147712 
## estimated kurtosis:  2.066642
  # Fit various distributions
    nd <- fitdist(sub_trait$value, "norm")
## $start.arg
## $start.arg$mean
## [1] 311.551
## 
## $start.arg$sd
## [1] 170.2334
## 
## 
## $fix.arg
## NULL
    gd <- fitdist(sub_trait$value, "gamma")
## $start.arg
## $start.arg$shape
## [1] 3.349413
## 
## $start.arg$rate
## [1] 0.01075077
## 
## 
## $fix.arg
## NULL
    wd <- fitdist(sub_trait$value, "weibull")
## $start.arg
## $start.arg$shape
## [1] 1.556164
## 
## $start.arg$scale
## [1] 360.2923
## 
## 
## $fix.arg
## NULL
    lnd <- fitdist(sub_trait$value, "lnorm")
## $start.arg
## $start.arg$meanlog
## [1] 5.519345
## 
## $start.arg$sdlog
## [1] 0.7675485
## 
## 
## $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.06491519 0.09755319 0.08208289 0.1382808
## Cramer-von Mises statistic   0.09472074 0.25095048 0.12686354 0.5370490
## Anderson-Darling statistic   0.70245530 1.55418026 0.92079324 3.2428572
## 
## Goodness-of-fit criteria
##                                  Normal    gamma  Weibull LogNormal
## Akaike's Information Criterion 1420.120 1422.775 1413.881  1445.526
## Bayesian Information Criterion 1425.484 1428.139 1419.245  1450.890

–> 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"), 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()
EPOC
Species mean n sd ci lower_ci upper_ci min max
Jasus edwardsii 286.2 54 158 43.1 243.0 329.3 23.1 623.3
Sagmariasus verreauxi 336.9 54 181 49.4 287.6 386.3 20.1 687.1

Grouped by species and acclimation temperature

# 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()
EPOC
Species AccTemp mean n sd ci lower_ci upper_ci min max
Jasus edwardsii 14 375.6 18 130.3 64.5 311.1 440.1 184.1 594.2
Jasus edwardsii 17.5 326.2 18 155.0 76.8 249.4 402.9 36.2 623.3
Jasus edwardsii 21.5 156.7 18 94.6 46.8 109.9 203.6 23.1 324.6
Sagmariasus verreauxi 14 398.4 18 215.5 106.7 291.7 505.1 72.2 687.1
Sagmariasus verreauxi 17.5 374.2 18 109.8 54.4 319.8 428.6 188.3 581.9
Sagmariasus verreauxi 21.5 238.3 18 167.1 82.7 155.5 321.0 20.1 567.5

By species and experimental temperature

# Summary statistics
  substrait_stats <- 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,
                  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()
EPOC
Species ExpTemp AccTemp mean n sd ci lower_ci upper_ci min max
Jasus edwardsii 14 14 272.4 6 113.8 113.7 158.7 386.0 184.1 479.0
Jasus edwardsii 14 17.5 290.7 6 83.9 83.8 206.8 374.5 151.8 378.2
Jasus edwardsii 14 21.5 84.8 6 55.5 55.5 29.4 140.3 23.1 164.5
Jasus edwardsii 17.5 14 436.7 6 147.0 146.9 289.8 583.6 237.3 594.2
Jasus edwardsii 17.5 17.5 216.0 6 137.7 137.5 78.5 353.6 36.2 414.4
Jasus edwardsii 17.5 21.5 204.6 6 86.3 86.2 118.4 290.8 90.3 324.6
Jasus edwardsii 21.5 14 417.7 6 61.1 61.0 356.6 478.7 311.4 484.3
Jasus edwardsii 21.5 17.5 471.8 6 119.0 118.8 353.0 590.7 289.6 623.3
Jasus edwardsii 21.5 21.5 180.8 6 101.2 101.1 79.7 281.9 69.1 312.2
Sagmariasus verreauxi 14 14 263.0 6 190.6 190.4 72.6 453.4 74.3 483.1
Sagmariasus verreauxi 14 17.5 314.4 6 56.8 56.7 257.7 371.1 238.2 393.1
Sagmariasus verreauxi 14 21.5 174.2 6 143.4 143.3 31.0 317.5 20.1 428.8
Sagmariasus verreauxi 17.5 14 457.1 6 212.1 211.9 245.2 669.0 72.2 633.3
Sagmariasus verreauxi 17.5 17.5 373.5 6 98.7 98.6 274.9 472.2 268.8 536.7
Sagmariasus verreauxi 17.5 21.5 198.9 6 106.1 106.0 92.9 304.9 35.4 331.2
Sagmariasus verreauxi 21.5 14 475.1 6 208.2 207.9 267.1 683.0 121.1 687.1
Sagmariasus verreauxi 21.5 17.5 434.6 6 139.2 139.0 295.6 573.6 188.3 581.9
Sagmariasus verreauxi 21.5 21.5 341.6 6 208.5 208.3 133.3 549.9 55.3 567.5

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 
##   1372.2   1431.2   -664.1   1328.2       86 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.27183 -0.51859  0.03203  0.56554  2.06072 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  AntTag   (Intercept) 6146     78.40   
##  Residual             8727     93.42   
## Number of obs: 108, groups:  AntTag, 39
## 
## Fixed effects:
##                                                        Estimate Std. Error
## (Intercept)                                           404.04837  109.05957
## ExpTemp17.5                                           164.34167   53.93567
## ExpTemp21.5                                           126.33407   55.67596
## AccTemp17.5                                             5.76447   70.46407
## AccTemp21.5                                          -205.08643   71.13506
## SpeciesSagmariasus verreauxi                          -10.37149   69.45235
## SexMale                                               -24.77882   34.55668
## Weight_g                                               -0.09626    0.08857
## ExpTemp17.5:AccTemp17.5                              -236.18767   76.32045
## ExpTemp21.5:AccTemp17.5                                54.80760   77.51690
## ExpTemp17.5:AccTemp21.5                               -44.07131   76.27798
## ExpTemp21.5:AccTemp21.5                               -29.86205   77.52935
## ExpTemp17.5:SpeciesSagmariasus verreauxi               13.94273   77.50517
## ExpTemp21.5:SpeciesSagmariasus verreauxi               69.91533   78.76059
## AccTemp17.5:SpeciesSagmariasus verreauxi               40.86027   99.18984
## AccTemp21.5:SpeciesSagmariasus verreauxi              112.56804   99.93895
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi  115.67153  109.96846
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi -127.14657  110.51267
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi -111.45928  108.79288
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi   -0.95355  109.72861
##                                                              df t value
## (Intercept)                                            49.14455   3.705
## ExpTemp17.5                                            69.35486   3.047
## ExpTemp21.5                                            76.64227   2.269
## AccTemp17.5                                            84.48444   0.082
## AccTemp21.5                                            83.66140  -2.883
## SpeciesSagmariasus verreauxi                           88.90528  -0.149
## SexMale                                                39.56152  -0.717
## Weight_g                                               40.97006  -1.087
## ExpTemp17.5:AccTemp17.5                                69.50981  -3.095
## ExpTemp21.5:AccTemp17.5                                73.07288   0.707
## ExpTemp17.5:AccTemp21.5                                69.35987  -0.578
## ExpTemp21.5:AccTemp21.5                                73.08895  -0.385
## ExpTemp17.5:SpeciesSagmariasus verreauxi               73.17693   0.180
## ExpTemp21.5:SpeciesSagmariasus verreauxi               76.84449   0.888
## AccTemp17.5:SpeciesSagmariasus verreauxi               86.89592   0.412
## AccTemp21.5:SpeciesSagmariasus verreauxi               84.49161   1.126
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi   74.09935   1.052
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi   75.01382  -1.151
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi   71.43833  -1.025
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi   73.35665  -0.009
##                                                      Pr(>|t|)    
## (Intercept)                                          0.000536 ***
## ExpTemp17.5                                          0.003268 ** 
## ExpTemp21.5                                          0.026076 *  
## AccTemp17.5                                          0.934994    
## AccTemp21.5                                          0.005005 ** 
## SpeciesSagmariasus verreauxi                         0.881629    
## SexMale                                              0.477556    
## Weight_g                                             0.283487    
## ExpTemp17.5:AccTemp17.5                              0.002838 ** 
## ExpTemp21.5:AccTemp17.5                              0.481790    
## ExpTemp17.5:AccTemp21.5                              0.565289    
## ExpTemp21.5:AccTemp21.5                              0.701229    
## ExpTemp17.5:SpeciesSagmariasus verreauxi             0.857733    
## ExpTemp21.5:SpeciesSagmariasus verreauxi             0.377475    
## AccTemp17.5:SpeciesSagmariasus verreauxi             0.681398    
## AccTemp21.5:SpeciesSagmariasus verreauxi             0.263199    
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi 0.296282    
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi 0.253586    
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi 0.309052    
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi 0.993090    
## ---
## 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                 383511  191756     2 71.826 21.9722 3.586e-08 ***
## AccTemp                 241527  120764     2 38.043 13.8377 3.048e-05 ***
## Species                  26472   26472     1 38.196  3.0333   0.08962 .  
## Sex                       4487    4487     1 39.562  0.5142   0.47756    
## Weight_g                 10308   10308     1 40.970  1.1811   0.28349    
## ExpTemp:AccTemp         121937   30484     4 71.841  3.4930   0.01157 *  
## ExpTemp:Species           3244    1622     2 72.144  0.1859   0.83078    
## AccTemp:Species           8264    4132     2 38.364  0.4735   0.62642    
## ExpTemp:AccTemp:Species 101204   25301     4 72.120  2.8991   0.02772 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  # Textual report of the model
    report(complex_model)
## We fitted a linear mixed model (estimated using ML and nloptwrap optimizer) to predict value with ExpTemp, AccTemp, Species, Sex and Weight_g (formula: value ~ ExpTemp * AccTemp * Species + Sex + Weight_g). The model included AntTag as random effect (formula: ~1 | AntTag). The model's total explanatory power is substantial (conditional R2 = 0.70) and the part related to the fixed effects alone (marginal R2) is of 0.48. The model's intercept, corresponding to ExpTemp = 14, AccTemp = 14, Species = Jasus edwardsii, Sex = Female and Weight_g = 0, is at 404.05 (95% CI [190.30, 617.80], t(86) = 3.70, p < .001). Within this model:
## 
##   - The effect of ExpTemp [17.5] is statistically significant and positive (beta = 164.34, 95% CI [58.63, 270.05], t(86) = 3.05, p = 0.002; Std. beta = 0.96, 95% CI [0.34, 1.58])
##   - The effect of ExpTemp [21.5] is statistically significant and positive (beta = 126.33, 95% CI [17.21, 235.46], t(86) = 2.27, p = 0.023; Std. beta = 0.74, 95% CI [0.10, 1.38])
##   - The effect of AccTemp [17.5] is statistically non-significant and positive (beta = 5.76, 95% CI [-132.34, 143.87], t(86) = 0.08, p = 0.935; Std. beta = 0.03, 95% CI [-0.77, 0.84])
##   - The effect of AccTemp [21.5] is statistically significant and negative (beta = -205.09, 95% CI [-344.51, -65.66], t(86) = -2.88, p = 0.004; Std. beta = -1.20, 95% CI [-2.01, -0.38])
##   - The effect of Species [Sagmariasus verreauxi] is statistically non-significant and negative (beta = -10.37, 95% CI [-146.50, 125.75], t(86) = -0.15, p = 0.881; Std. beta = -0.06, 95% CI [-0.86, 0.74])
##   - The effect of Sex [Male] is statistically non-significant and negative (beta = -24.78, 95% CI [-92.51, 42.95], t(86) = -0.72, p = 0.473; Std. beta = -0.14, 95% CI [-0.54, 0.25])
##   - The effect of Weight_g is statistically non-significant and negative (beta = -0.10, 95% CI [-0.27, 0.08], t(86) = -1.09, p = 0.277; Std. beta = -0.11, 95% CI [-0.30, 0.09])
##   - The interaction effect of AccTemp [17.5] on ExpTemp [17.5] is statistically significant and negative (beta = -236.19, 95% CI [-385.77, -86.60], t(86) = -3.09, p = 0.002; Std. beta = -1.38, 95% CI [-2.26, -0.51])
##   - The interaction effect of AccTemp [17.5] on ExpTemp [21.5] is statistically non-significant and positive (beta = 54.81, 95% CI [-97.12, 206.74], t(86) = 0.71, p = 0.480; Std. beta = 0.32, 95% CI [-0.57, 1.21])
##   - The interaction effect of AccTemp [21.5] on ExpTemp [17.5] is statistically non-significant and negative (beta = -44.07, 95% CI [-193.57, 105.43], t(86) = -0.58, p = 0.563; Std. beta = -0.26, 95% CI [-1.13, 0.62])
##   - The interaction effect of AccTemp [21.5] on ExpTemp [21.5] is statistically non-significant and negative (beta = -29.86, 95% CI [-181.82, 122.09], t(86) = -0.39, p = 0.700; Std. beta = -0.17, 95% CI [-1.06, 0.71])
##   - The interaction effect of Species [Sagmariasus verreauxi] on ExpTemp [17.5] is statistically non-significant and positive (beta = 13.94, 95% CI [-137.96, 165.85], t(86) = 0.18, p = 0.857; Std. beta = 0.08, 95% CI [-0.81, 0.97])
##   - The interaction effect of Species [Sagmariasus verreauxi] on ExpTemp [21.5] is statistically non-significant and positive (beta = 69.92, 95% CI [-84.45, 224.28], t(86) = 0.89, p = 0.375; Std. beta = 0.41, 95% CI [-0.49, 1.31])
##   - The interaction effect of Species [Sagmariasus verreauxi] on AccTemp [17.5] is statistically non-significant and positive (beta = 40.86, 95% CI [-153.55, 235.27], t(86) = 0.41, p = 0.680; Std. beta = 0.24, 95% CI [-0.90, 1.38])
##   - The interaction effect of Species [Sagmariasus verreauxi] on AccTemp [21.5] is statistically non-significant and positive (beta = 112.57, 95% CI [-83.31, 308.44], t(86) = 1.13, p = 0.260; Std. beta = 0.66, 95% CI [-0.49, 1.80])
##   - The interaction effect of Species [Sagmariasus verreauxi] on (ExpTemp [17.5] * AccTemp [17.5]) is statistically non-significant and positive (beta = 115.67, 95% CI [-99.86, 331.21], t(86) = 1.05, p = 0.293; Std. beta = 0.68, 95% CI [-0.58, 1.94])
##   - The interaction effect of Species [Sagmariasus verreauxi] on (ExpTemp [21.5] * AccTemp [17.5]) is statistically non-significant and negative (beta = -127.15, 95% CI [-343.75, 89.45], t(86) = -1.15, p = 0.250; Std. beta = -0.74, 95% CI [-2.01, 0.52])
##   - The interaction effect of Species [Sagmariasus verreauxi] on (ExpTemp [17.5] * AccTemp [21.5]) is statistically non-significant and negative (beta = -111.46, 95% CI [-324.69, 101.77], t(86) = -1.02, p = 0.306; Std. beta = -0.65, 95% CI [-1.90, 0.60])
##   - The interaction effect of Species [Sagmariasus verreauxi] on (ExpTemp [21.5] * AccTemp [21.5]) is statistically non-significant and negative (beta = -0.95, 95% CI [-216.02, 214.11], t(86) = -8.69e-03, p = 0.993; Std. beta = -5.58e-03, 95% CI [-1.26, 1.25])
## 
## 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.

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 -664.11 1372.2                         
## (1 | AntTag)          0   21 -672.10 1386.2 15.997  1  6.343e-05 ***
## ---
## 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  Pr(>F)
## Sex                              1   4487  4487.2     1 39.562  0.5142 0.47756
## Weight_g                         2  16675 16675.0     1 40.813  1.9205 0.17333
## ExpTemp:AccTemp:Species          0 105765 26441.1     4 72.421  2.9797 0.02458
##                          
## Sex                      
## Weight_g                 
## ExpTemp:AccTemp:Species *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## value ~ ExpTemp + AccTemp + Species + (1 | AntTag) + ExpTemp:AccTemp + ExpTemp:Species + AccTemp:Species + ExpTemp:AccTemp:Species
  # 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 * Species + (1 | AntTag)

    # Remodel using Maximum likelyhood estimation
      final_model <- lmer(value ~ ExpTemp *
                            AccTemp *
                            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 * AccTemp * Species + (1 | AntTag)
##    Data: sub_trait
## 
## REML criterion at convergence: 1157.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.09029 -0.48842  0.02245  0.54816  1.96764 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  AntTag   (Intercept)  7591     87.13  
##  Residual             10687    103.38  
## Number of obs: 108, groups:  AntTag, 39
## 
## Fixed effects:
##                                                      Estimate Std. Error
## (Intercept)                                           277.726     54.524
## ExpTemp17.5                                           164.342     59.685
## ExpTemp21.5                                           128.970     61.582
## AccTemp17.5                                            12.946     77.584
## AccTemp21.5                                          -192.913     77.584
## SpeciesSagmariasus verreauxi                           -6.515     76.912
## ExpTemp17.5:AccTemp17.5                              -239.000     84.408
## ExpTemp21.5:AccTemp17.5                                52.172     85.760
## ExpTemp17.5:AccTemp21.5                               -44.577     84.408
## ExpTemp21.5:AccTemp21.5                               -33.003     85.760
## ExpTemp17.5:SpeciesSagmariasus verreauxi               11.787     85.760
## ExpTemp21.5:SpeciesSagmariasus verreauxi               65.124     87.090
## AccTemp17.5:SpeciesSagmariasus verreauxi               25.837    109.107
## AccTemp21.5:SpeciesSagmariasus verreauxi               95.933    109.582
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi  126.674    121.282
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi -121.366    122.227
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi -106.857    120.330
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi    6.284    121.282
##                                                            df t value Pr(>|t|)
## (Intercept)                                            74.155   5.094 2.59e-06
## ExpTemp17.5                                            58.324   2.753  0.00785
## ExpTemp21.5                                            64.328   2.094  0.04018
## AccTemp17.5                                            71.542   0.167  0.86795
## AccTemp21.5                                            71.542  -2.486  0.01524
## SpeciesSagmariasus verreauxi                           74.709  -0.085  0.93272
## ExpTemp17.5:AccTemp17.5                                58.324  -2.831  0.00635
## ExpTemp21.5:AccTemp17.5                                61.388   0.608  0.54520
## ExpTemp17.5:AccTemp21.5                                58.324  -0.528  0.59943
## ExpTemp21.5:AccTemp21.5                                61.388  -0.385  0.70169
## ExpTemp17.5:SpeciesSagmariasus verreauxi               61.388   0.137  0.89113
## ExpTemp21.5:SpeciesSagmariasus verreauxi               64.328   0.748  0.45732
## AccTemp17.5:SpeciesSagmariasus verreauxi               73.400   0.237  0.81347
## AccTemp21.5:SpeciesSagmariasus verreauxi               71.816   0.875  0.38425
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi   61.388   1.044  0.30037
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi   62.875  -0.993  0.32454
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi   59.870  -0.888  0.37808
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi   61.388   0.052  0.95885
##                                                         
## (Intercept)                                          ***
## ExpTemp17.5                                          ** 
## ExpTemp21.5                                          *  
## AccTemp17.5                                             
## AccTemp21.5                                          *  
## SpeciesSagmariasus verreauxi                            
## ExpTemp17.5:AccTemp17.5                              ** 
## ExpTemp21.5:AccTemp17.5                                 
## ExpTemp17.5:AccTemp21.5                                 
## ExpTemp21.5:AccTemp21.5                                 
## ExpTemp17.5:SpeciesSagmariasus verreauxi                
## ExpTemp21.5:SpeciesSagmariasus verreauxi                
## AccTemp17.5:SpeciesSagmariasus verreauxi                
## AccTemp21.5:SpeciesSagmariasus verreauxi                
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi    
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi    
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi    
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 18 > 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(final_model)
## Type III Analysis of Variance Table with Satterthwaite's method
##                         Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## ExpTemp                 386129  193064     2 60.352 18.0654 7.102e-07 ***
## AccTemp                 233698  116849     2 32.806 10.9338 0.0002298 ***
## Species                  21595   21595     1 32.794  2.0207 0.1646089    
## ExpTemp:AccTemp         118592   29648     4 60.317  2.7742 0.0349077 *  
## ExpTemp:Species           3263    1632     2 60.352  0.1527 0.8587350    
## AccTemp:Species           5820    2910     2 32.806  0.2723 0.7633360    
## ExpTemp:AccTemp:Species 105752   26438     4 60.317  2.4739 0.0537655 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      report(anova(final_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) = 18.07, p < .001; Eta2 (partial) = 0.37, 90% CI [0.21, 0.50])
##   - The main effect of AccTemp is statistically significant and large (F(2) = 10.93, p < .001; Eta2 (partial) = 0.40, 90% CI [0.17, 0.56])
##   - The main effect of Species is statistically not significant and small (F(1) = 2.02, p = 0.165; Eta2 (partial) = 0.06, 90% CI [0.00, 0.22])
##   - The interaction between ExpTemp and AccTemp is statistically significant and large (F(4) = 2.77, p = 0.035; Eta2 (partial) = 0.16, 90% CI [6.69e-03, 0.26])
##   - The interaction between ExpTemp and Species is statistically not significant and very small (F(2) = 0.15, p = 0.859; Eta2 (partial) = 5.03e-03, 90% CI [0.00, 0.04])
##   - The interaction between AccTemp and Species is statistically not significant and small (F(2) = 0.27, p = 0.763; Eta2 (partial) = 0.02, 90% CI [0.00, 0.10])
##   - The interaction between ExpTemp, AccTemp and Species is statistically not significant and large (F(4) = 2.47, p = 0.054; Eta2 (partial) = 0.14, 90% CI [0.00, 0.25])
## 
## Effect sizes were labelled following Field's (2013) recommendations.
    # Table output from sjPlot
      tab_model(final_model)
  value
Predictors Estimates CI p
(Intercept) 277.73 170.86 – 384.59 <0.001
ExpTemp [17.5] 164.34 47.36 – 281.32 0.006
ExpTemp [21.5] 128.97 8.27 – 249.67 0.036
AccTemp [17.5] 12.95 -139.12 – 165.01 0.867
AccTemp [21.5] -192.91 -344.97 – -40.85 0.013
Species [Sagmariasus
verreauxi]
-6.51 -157.26 – 144.23 0.932
ExpTemp [17.5] * AccTemp
[17.5]
-239.00 -404.44 – -73.56 0.005
ExpTemp [21.5] * AccTemp
[17.5]
52.17 -115.91 – 220.26 0.543
ExpTemp [17.5] * AccTemp
[21.5]
-44.58 -210.01 – 120.86 0.597
ExpTemp [21.5] * AccTemp
[21.5]
-33.00 -201.09 – 135.08 0.700
ExpTemp [17.5] * Species
[Sagmariasus verreauxi]
11.79 -156.30 – 179.87 0.891
ExpTemp [21.5] * Species
[Sagmariasus verreauxi]
65.12 -105.57 – 235.82 0.455
AccTemp [17.5] * Species
[Sagmariasus verreauxi]
25.84 -188.01 – 239.68 0.813
AccTemp [21.5] * Species
[Sagmariasus verreauxi]
95.93 -118.84 – 310.71 0.381
(ExpTemp [17.5] * AccTemp
[17.5]) * Species
[Sagmariasus verreauxi]
126.67 -111.03 – 364.38 0.296
(ExpTemp [21.5] * AccTemp
[17.5]) * Species
[Sagmariasus verreauxi]
-121.37 -360.93 – 118.19 0.321
(ExpTemp [17.5] * AccTemp
[21.5]) * Species
[Sagmariasus verreauxi]
-106.86 -342.70 – 128.99 0.375
(ExpTemp [21.5] * AccTemp
[21.5]) * Species
[Sagmariasus verreauxi]
6.28 -231.42 – 243.99 0.959
Random Effects
σ2 10686.96
τ00 AntTag 7591.36
ICC 0.42
N AntTag 39
Observations 108
Marginal R2 / Conditional R2 0.422 / 0.662
    # 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, AccTemp and Species (formula: value ~ ExpTemp * AccTemp * Species). The model included AntTag as random effect (formula: ~1 | AntTag). The model's total explanatory power is substantial (conditional R2 = 0.66) and the part related to the fixed effects alone (marginal R2) is of 0.42. The model's intercept, corresponding to ExpTemp = 14, AccTemp = 14 and Species = Jasus edwardsii, is at 277.73 (95% CI [170.86, 384.59], t(88) = 5.09, p < .001). Within this model:
## 
##   - The effect of ExpTemp [17.5] is statistically significant and positive (beta = 164.34, 95% CI [47.36, 281.32], t(88) = 2.75, p = 0.006; Std. beta = 0.96, 95% CI [0.28, 1.64])
##   - The effect of ExpTemp [21.5] is statistically significant and positive (beta = 128.97, 95% CI [8.27, 249.67], t(88) = 2.09, p = 0.036; Std. beta = 0.75, 95% CI [0.05, 1.46])
##   - The effect of AccTemp [17.5] is statistically non-significant and positive (beta = 12.95, 95% CI [-139.12, 165.01], t(88) = 0.17, p = 0.867; Std. beta = 0.08, 95% CI [-0.81, 0.96])
##   - The effect of AccTemp [21.5] is statistically significant and negative (beta = -192.91, 95% CI [-344.97, -40.85], t(88) = -2.49, p = 0.013; Std. beta = -1.13, 95% CI [-2.02, -0.24])
##   - The effect of Species [Sagmariasus verreauxi] is statistically non-significant and negative (beta = -6.51, 95% CI [-157.26, 144.23], t(88) = -0.08, p = 0.932; Std. beta = -0.04, 95% CI [-0.92, 0.84])
##   - The interaction effect of AccTemp [17.5] on ExpTemp [17.5] is statistically significant and negative (beta = -239.00, 95% CI [-404.44, -73.56], t(88) = -2.83, p = 0.005; Std. beta = -1.40, 95% CI [-2.36, -0.43])
##   - The interaction effect of AccTemp [17.5] on ExpTemp [21.5] is statistically non-significant and positive (beta = 52.17, 95% CI [-115.91, 220.26], t(88) = 0.61, p = 0.543; Std. beta = 0.31, 95% CI [-0.68, 1.29])
##   - The interaction effect of AccTemp [21.5] on ExpTemp [17.5] is statistically non-significant and negative (beta = -44.58, 95% CI [-210.01, 120.86], t(88) = -0.53, p = 0.597; Std. beta = -0.26, 95% CI [-1.23, 0.71])
##   - The interaction effect of AccTemp [21.5] on ExpTemp [21.5] is statistically non-significant and negative (beta = -33.00, 95% CI [-201.09, 135.08], t(88) = -0.38, p = 0.700; Std. beta = -0.19, 95% CI [-1.18, 0.79])
##   - The interaction effect of Species [Sagmariasus verreauxi] on ExpTemp [17.5] is statistically non-significant and positive (beta = 11.79, 95% CI [-156.30, 179.87], t(88) = 0.14, p = 0.891; Std. beta = 0.07, 95% CI [-0.91, 1.05])
##   - The interaction effect of Species [Sagmariasus verreauxi] on ExpTemp [21.5] is statistically non-significant and positive (beta = 65.12, 95% CI [-105.57, 235.82], t(88) = 0.75, p = 0.455; Std. beta = 0.38, 95% CI [-0.62, 1.38])
##   - The interaction effect of Species [Sagmariasus verreauxi] on AccTemp [17.5] is statistically non-significant and positive (beta = 25.84, 95% CI [-188.01, 239.68], t(88) = 0.24, p = 0.813; Std. beta = 0.15, 95% CI [-1.10, 1.40])
##   - The interaction effect of Species [Sagmariasus verreauxi] on AccTemp [21.5] is statistically non-significant and positive (beta = 95.93, 95% CI [-118.84, 310.71], t(88) = 0.88, p = 0.381; Std. beta = 0.56, 95% CI [-0.69, 1.82])
##   - The interaction effect of Species [Sagmariasus verreauxi] on (ExpTemp [17.5] * AccTemp [17.5]) is statistically non-significant and positive (beta = 126.67, 95% CI [-111.03, 364.38], t(88) = 1.04, p = 0.296; Std. beta = 0.74, 95% CI [-0.65, 2.13])
##   - The interaction effect of Species [Sagmariasus verreauxi] on (ExpTemp [21.5] * AccTemp [17.5]) is statistically non-significant and negative (beta = -121.37, 95% CI [-360.93, 118.19], t(88) = -0.99, p = 0.321; Std. beta = -0.71, 95% CI [-2.11, 0.69])
##   - The interaction effect of Species [Sagmariasus verreauxi] on (ExpTemp [17.5] * AccTemp [21.5]) is statistically non-significant and negative (beta = -106.86, 95% CI [-342.70, 128.99], t(88) = -0.89, p = 0.375; Std. beta = -0.62, 95% CI [-2.00, 0.75])
##   - The interaction effect of Species [Sagmariasus verreauxi] on (ExpTemp [21.5] * AccTemp [21.5]) is statistically non-significant and positive (beta = 6.28, 95% CI [-231.42, 243.99], t(88) = 0.05, p = 0.959; Std. beta = 0.04, 95% CI [-1.35, 1.43])
## 
## 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 386128.9 193064.5 60.4 18.07 (2) 7.102e-07 0.37 [0.21-0.5]
AccTemp 233698.1 116849.1 32.8 10.93 (2) 2.298e-04 0.4 [0.17-0.56]
Species 21595.1 21595.1 32.8 2.02 (1) 1.646e-01 0.06 [0-0.22]
ExpTemp:AccTemp 118591.6 29647.9 60.3 2.77 (4) 3.491e-02 0.16 [0.01-0.26]
ExpTemp:Species 3263.4 1631.7 60.4 0.15 (2) 8.587e-01 0.01 [0-0.04]
AccTemp:Species 5820.0 2910.0 32.8 0.27 (2) 7.633e-01 0.02 [0-0.1]
ExpTemp:AccTemp:Species 105752.4 26438.1 60.3 2.47 (4) 5.377e-02 0.14 [0-0.25]

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)

    plotResiduals(simulationOutput, sub_trait$Species)

–> 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 | Species, type = "response", adjust = "bonferroni")

# Display interaction plots
  emmip(final_model, AccTemp ~ 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:4,9,10)] %>% mutate(p.value = round(p.value, 3))
    # Rename column headers
      names(contrasts_expTemp)[1:6] <- c("Contrast", "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_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'))
Contrast ExpTemp Species Difference t-ratio p contrasts_sign Trait
14 - 17.5 14 Jasus edwardsii -12.94558 -0.1666996 1 No EPOC
14 - 21.5 14 Jasus edwardsii 192.91275 2.4841280 0.046 Yes EPOC
17.5 - 21.5 14 Jasus edwardsii 205.85833 2.6373083 0.031 Yes EPOC
14 - 17.5 17.5 Jasus edwardsii 226.05442 2.9108916 0.014 Yes EPOC
14 - 21.5 17.5 Jasus edwardsii 237.48942 3.0581396 0.009 Yes EPOC
17.5 - 21.5 17.5 Jasus edwardsii 11.43500 0.1464970 1 No EPOC
14 - 17.5 21.5 Jasus edwardsii -65.11770 -0.8408618 1 No EPOC
14 - 21.5 21.5 Jasus edwardsii 225.91563 2.9172380 0.014 Yes EPOC
17.5 - 21.5 21.5 Jasus edwardsii 291.03333 3.7285090 0.001 Yes EPOC
14 - 17.5 14 Sagmariasus verreauxi -38.78229 -0.5048325 1 No EPOC
14 - 21.5 14 Sagmariasus verreauxi 96.97970 1.2522943 0.644 No EPOC
17.5 - 21.5 14 Sagmariasus verreauxi 135.76200 1.7530883 0.252 No EPOC
14 - 17.5 17.5 Sagmariasus verreauxi 73.54351 0.9519209 1 No EPOC
14 - 21.5 17.5 Sagmariasus verreauxi 248.41362 3.1988100 0.006 Yes EPOC
17.5 - 21.5 17.5 Sagmariasus verreauxi 174.87011 2.2517939 0.082 No EPOC
14 - 17.5 21.5 Sagmariasus verreauxi 30.41184 0.3936400 1 No EPOC
14 - 21.5 21.5 Sagmariasus verreauxi 123.69862 1.5928610 0.347 No EPOC
17.5 - 21.5 21.5 Sagmariasus verreauxi 93.28678 1.2012493 0.701 No EPOC

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 | Species, type = "response", adjust = "bonferroni")

# Display interaction plots
  emmip(final_model, AccTemp ~ 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:4,9,10)] %>% mutate(p.value = round(p.value, 3))
    # Rename column headers
      names(contrasts_expTemp)[1:6] <- c("Contrast", "AccTemp", "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_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'))
Contrast AccTemp Species Difference t-ratio p contrasts_sign Trait
14 - 17.5 14 Jasus edwardsii -164.34167 -2.7534752 0.024 Yes EPOC
14 - 21.5 14 Jasus edwardsii -128.96955 -2.0874881 0.123 No EPOC
17.5 - 21.5 14 Jasus edwardsii 35.37212 0.5725296 1 No EPOC
14 - 17.5 17.5 Jasus edwardsii 74.65833 1.2508688 0.648 No EPOC
14 - 21.5 17.5 Jasus edwardsii -181.14167 -3.0349521 0.011 Yes EPOC
17.5 - 21.5 17.5 Jasus edwardsii -255.80000 -4.2858209 0 Yes EPOC
14 - 17.5 21.5 Jasus edwardsii -119.76500 -2.0066119 0.149 No EPOC
14 - 21.5 21.5 Jasus edwardsii -95.96667 -1.6078809 0.34 No EPOC
17.5 - 21.5 21.5 Jasus edwardsii 23.79833 0.3987310 1 No EPOC
14 - 17.5 14 Sagmariasus verreauxi -176.12892 -2.8508050 0.018 Yes EPOC
14 - 21.5 14 Sagmariasus verreauxi -194.09392 -3.1415847 0.008 Yes EPOC
17.5 - 21.5 14 Sagmariasus verreauxi -17.96500 -0.3009960 1 No EPOC
14 - 17.5 17.5 Sagmariasus verreauxi -63.80312 -1.0327109 0.917 No EPOC
14 - 21.5 17.5 Sagmariasus verreauxi -124.89979 -2.0216154 0.142 No EPOC
17.5 - 21.5 17.5 Sagmariasus verreauxi -61.09667 -1.0236488 0.931 No EPOC
14 - 17.5 21.5 Sagmariasus verreauxi -24.69500 -0.4137543 1 No EPOC
14 - 21.5 21.5 Sagmariasus verreauxi -167.37500 -2.8042974 0.021 Yes EPOC
17.5 - 21.5 21.5 Sagmariasus verreauxi -142.68000 -2.3905431 0.06 No EPOC

9.C 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 ~  Species | ExpTemp , type = "response", adjust = "bonferroni")
## NOTE: Results may be misleading due to involvement in interactions
# Display interaction plots
  emmip(final_model, ExpTemp ~ Species)
## NOTE: Results may be misleading due to involvement in interactions

# 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("Contrast", "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_Species_", 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'))
Contrast ExpTemp Difference t-ratio p contrasts_sign Trait
Jasus edwardsii - Sagmariasus verreauxi 14 -34.07520 -0.761412 0.449 No EPOC
Jasus edwardsii - Sagmariasus verreauxi 17.5 -52.46810 -1.170224 0.246 No EPOC
Jasus edwardsii - Sagmariasus verreauxi 21.5 -60.83881 -1.358181 0.179 No EPOC

10. R Session Information

report_system()
## Analyses were conducted using the R Statistical language (version 4.1.2; R Core Team, 2021) on Windows 10 x64 (build 19043)
report_packages()
##   - ggpubr (version 0.4.0; Alboukadel Kassambara, 2020)
##   - effectsize (version 0.4.5; Ben-Shachar M et al., 2020)
##   - bbmle (version 1.0.24; Ben Bolker and R Development Core Team, 2021)
##   - Matrix (version 1.3.4; Douglas Bates and Martin Maechler, 2021)
##   - lme4 (version 1.1.27.1; Douglas Bates et al., 2015)
##   - DHARMa (version 0.4.4; Florian Hartig, 2021)
##   - pander (version 0.6.4; Gergely Daróczi and Roman Tsegelskyi, 2021)
##   - ggplot2 (version 3.3.5; Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016.)
##   - plyr (version 1.8.6; Hadley Wickham, 2011)
##   - stringr (version 1.4.0; Hadley Wickham, 2019)
##   - forcats (version 0.5.1; Hadley Wickham, 2021)
##   - tidyr (version 1.1.3; Hadley Wickham, 2021)
##   - readxl (version 1.3.1; Hadley Wickham and Jennifer Bryan, 2019)
##   - readr (version 1.4.0; Hadley Wickham and Jim Hester, 2020)
##   - dplyr (version 1.0.7; Hadley Wickham et al., 2021)
##   - kableExtra (version 1.3.4; Hao Zhu, 2021)
##   - tibble (version 3.1.2; Kirill Müller and Hadley Wickham, 2021)
##   - lmerTest (version 3.1.3; Kuznetsova A et al., 2017)
##   - purrr (version 0.3.4; Lionel Henry and Hadley Wickham, 2020)
##   - sjPlot (version 2.8.10; Lüdecke D, 2021)
##   - psycho (version 0.6.1; Makowski, 2018)
##   - report (version 0.4.0; Makowski et al., 2020)
##   - fitdistrplus (version 1.1.5; Marie Laure Delignette-Muller, Christophe Dutang, 2015)
##   - openxlsx (version 4.2.4; Philipp Schauberger and Alexander Walker, 2021)
##   - R (version 4.1.2; R Core Team, 2021)
##   - emmeans (version 1.6.2.1; Russell Lenth, 2021)
##   - survival (version 3.2.13; Therneau T, 2021)
##   - pbkrtest (version 0.5.1; Ulrich Halekoh, Søren Højsgaard, 2014)
##   - MASS (version 7.3.54; Venables et al., 2002)
##   - tidyverse (version 1.3.1; Wickham et al., 2019)
##   - knitr (version 1.33; Yihui Xie, 2021)