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

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)