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

# R document settings
# Echo = True means show all code chunks
knitr::opts_chunk$set(echo = TRUE)

2. Load data

  # Select trait
    trait <- "MMR"

 # Set working directory of executing R. script
    setwd(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("../../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 == "MMR")

4. Find best probability distribution of data to check assumptions for linear model

Find best distribution fit using fitdistrplus package

  # Plot data versus a range of distrubutions
    descdist(sub_trait$value, boot = 1000)

## summary statistics
## ------
## min:  62.81   max:  160.02 
## median:  108.025 
## mean:  108.1056 
## estimated sd:  22.05236 
## estimated skewness:  -0.01523508 
## estimated kurtosis:  2.469454
  # Fit various distributions
    nd <- fitdist(sub_trait$value, "norm")
## $start.arg
## $start.arg$mean
## [1] 108.1056
## 
## $start.arg$sd
## [1] 21.95003
## 
## 
## $fix.arg
## NULL
    gd <- fitdist(sub_trait$value, "gamma")
## $start.arg
## $start.arg$shape
## [1] 24.25641
## 
## $start.arg$rate
## [1] 0.2243769
## 
## 
## $fix.arg
## NULL
    wd <- fitdist(sub_trait$value, "weibull")
## $start.arg
## $start.arg$shape
## [1] 5.624781
## 
## $start.arg$scale
## [1] 117.0993
## 
## 
## $fix.arg
## NULL
    lnd <- fitdist(sub_trait$value, "lnorm")
## $start.arg
## $start.arg$meanlog
## [1] 4.661329
## 
## $start.arg$sdlog
## [1] 0.2123517
## 
## 
## $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.04961111 0.06203287 0.04928462 0.07189584
## Cramer-von Mises statistic   0.02731938 0.07516722 0.02481451 0.12058465
## Anderson-Darling statistic   0.20138464 0.48593740 0.20839501 0.76875941
## 
## Goodness-of-fit criteria
##                                  Normal    gamma  Weibull LogNormal
## Akaike's Information Criterion 977.6647 979.6881 978.6070  982.6433
## Bayesian Information Criterion 983.0290 985.0524 983.9713  988.0075

–> 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) +
        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()
MMR
Species mean n sd ci lower_ci upper_ci min max
Jasus edwardsii 101.7 54 20.2 5.5 96.1 107.2 62.8 148.7
Sagmariasus verreauxi 114.6 54 22.1 6.0 108.5 120.6 66.9 160.0

By species and experimental temperature

# 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()
MMR
Species ExpTemp mean n sd ci lower_ci upper_ci min max
Jasus edwardsii 14 90.6 18 16.4 8.1 82.5 98.8 62.8 118.9
Jasus edwardsii 17.5 100.3 18 19.3 9.6 90.7 109.9 64.3 144.3
Jasus edwardsii 21.5 114.0 18 18.3 9.1 104.9 123.1 80.0 148.7
Sagmariasus verreauxi 14 96.3 18 15.2 7.5 88.8 103.9 66.9 126.4
Sagmariasus verreauxi 17.5 116.0 18 16.6 8.2 107.8 124.3 70.8 137.7
Sagmariasus verreauxi 21.5 131.3 18 19.1 9.5 121.9 140.8 74.0 160.0

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
    anova(complex_model)
## Type III Analysis of Variance Table with Satterthwaite's method
##                          Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## ExpTemp                 12645.9  6323.0     2 70.370 72.1287 < 2.2e-16 ***
## AccTemp                   677.7   338.9     2 37.762  3.8657  0.029693 *  
## Species                  1036.8  1036.8     1 37.929 11.8275  0.001433 ** 
## Sex                        11.6    11.6     1 39.103  0.1321  0.718198    
## Weight_g                  447.8   447.8     1 44.202  5.1081  0.028791 *  
## ExpTemp:AccTemp          1001.4   250.4     4 70.436  2.8559  0.029724 *  
## ExpTemp:Species           526.4   263.2     2 70.844  3.0024  0.056028 .  
## AccTemp:Species           231.6   115.8     2 38.031  1.3209  0.278871    
## ExpTemp:AccTemp:Species   186.0    46.5     4 70.900  0.5305  0.713710    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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 -427.03 898.06                         
## (1 | AntTag)          0   21 -444.75 931.50 35.433  1  2.639e-09 ***
## ---
## 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
## Sex                              1   11.58   11.58     1 39.103  0.1321
## ExpTemp:AccTemp:Species          2  188.75   47.19     4 70.866  0.5389
## AccTemp:Species                  3  237.97  118.98     2 38.207  1.3193
## ExpTemp:Species                  4  525.24  262.62     2 70.869  2.9138
## Species                          0  995.55  995.55     1 38.184 10.2588
## Weight_g                         0  458.26  458.26     1 44.916  4.7222
## ExpTemp:AccTemp                  0 1013.40  253.35     4 70.386  2.6107
##                           Pr(>F)   
## Sex                     0.718198   
## ExpTemp:AccTemp:Species 0.707647   
## AccTemp:Species         0.279222   
## ExpTemp:Species         0.060799 . 
## Species                 0.002743 **
## Weight_g                0.035089 * 
## ExpTemp:AccTemp         0.042626 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## value ~ ExpTemp + AccTemp + Species + Weight_g + (1 | AntTag) + ExpTemp:AccTemp
  # Extract the model that step found:
    final_model <- get_model(step_result)
    final_model
## Linear mixed model fit by maximum likelihood  ['lmerModLmerTest']
## Formula: value ~ ExpTemp + AccTemp + Species + Weight_g + (1 | AntTag) +  
##     ExpTemp:AccTemp
##    Data: sub_trait
##       AIC       BIC    logLik  deviance  df.resid 
##  890.5019  925.3696 -432.2510  864.5019        95 
## Random effects:
##  Groups   Name        Std.Dev.
##  AntTag   (Intercept) 12.207  
##  Residual              9.851  
## Number of obs: 108, groups:  AntTag, 39
## Fixed Effects:
##                  (Intercept)                   ExpTemp17.5  
##                    120.35372                      14.62665  
##                  ExpTemp21.5                   AccTemp17.5  
##                     21.56607                       2.21054  
##                  AccTemp21.5  SpeciesSagmariasus verreauxi  
##                    -18.30465                      14.16481  
##                     Weight_g       ExpTemp17.5:AccTemp17.5  
##                     -0.02478                      -8.30831  
##      ExpTemp21.5:AccTemp17.5       ExpTemp17.5:AccTemp21.5  
##                      2.53309                       4.22518  
##      ExpTemp21.5:AccTemp21.5  
##                     13.86326

7.C Re-run model fit using REML estimation

Best model: ExpTemp * AccTemp + Species + Weight_g + (1 | AntTag), however retained interaction term with species as close to significant

  # Remodel using Maximum likelihood estimation
   final_model2 <- lmer(value ~ ExpTemp * 
                          AccTemp * 
                          Species +
                          Weight_g +
                          (1|AntTag), data = sub_trait, REML = TRUE,
                        contrasts = list(
                                        ExpTemp = "contr.sum",
                                        AccTemp = "contr.sum",
                                        Species = "contr.sum"
                                      ))

  # Model summary  
    summary(final_model)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: value ~ ExpTemp + AccTemp + Species + Weight_g + (1 | AntTag) +  
##     ExpTemp:AccTemp
##    Data: sub_trait
## 
##      AIC      BIC   logLik deviance df.resid 
##    890.5    925.4   -432.3    864.5       95 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3632 -0.5462  0.0977  0.5889  1.8819 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  AntTag   (Intercept) 149.02   12.207  
##  Residual              97.04    9.851  
## Number of obs: 108, groups:  AntTag, 39
## 
## Fixed effects:
##                               Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)                  120.35372   13.48976  47.90078   8.922 9.41e-12
## ExpTemp17.5                   14.62665    4.10981  71.31597   3.559 0.000667
## ExpTemp21.5                   21.56607    4.19779  73.96192   5.137 2.19e-06
## AccTemp17.5                    2.21054    6.25427  69.52310   0.353 0.724824
## AccTemp21.5                  -18.30465    6.33446  67.45980  -2.890 0.005181
## SpeciesSagmariasus verreauxi  14.16481    4.42245  38.18352   3.203 0.002743
## Weight_g                      -0.02478    0.01140  44.91586  -2.173 0.035089
## ExpTemp17.5:AccTemp17.5       -8.30831    5.81591  71.46023  -1.429 0.157490
## ExpTemp21.5:AccTemp17.5        2.53309    5.87350  72.60468   0.431 0.667546
## ExpTemp17.5:AccTemp21.5        4.22518    5.75244  70.08435   0.735 0.465094
## ExpTemp21.5:AccTemp21.5       13.86326    5.81451  71.46809   2.384 0.019771
##                                 
## (Intercept)                  ***
## ExpTemp17.5                  ***
## ExpTemp21.5                  ***
## AccTemp17.5                     
## AccTemp21.5                  ** 
## SpeciesSagmariasus verreauxi ** 
## Weight_g                     *  
## ExpTemp17.5:AccTemp17.5         
## ExpTemp21.5:AccTemp17.5         
## ExpTemp17.5:AccTemp21.5         
## ExpTemp21.5:AccTemp21.5      *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) ExT17.5 ExT21.5 AT17.5 AT21.5 SpcsSv Wght_g ET17.5:AT1
## ExpTemp17.5 -0.127                                                       
## ExpTemp21.5 -0.150  0.511                                                
## AccTemp17.5 -0.245  0.326   0.335                                        
## AccTemp21.5 -0.309  0.320   0.331   0.489                                
## SpcsSgmrssv -0.025  0.011   0.015  -0.019 -0.009                         
## Weight_g    -0.931 -0.027  -0.008   0.020  0.089 -0.148                  
## ET17.5:AT17  0.064 -0.707  -0.361  -0.460 -0.224 -0.007  0.046           
## ET21.5:AT17  0.103 -0.365  -0.715  -0.467 -0.236 -0.007  0.010  0.516    
## ET17.5:AT21  0.077 -0.715  -0.365  -0.232 -0.449 -0.010  0.034  0.506    
## ET21.5:AT21  0.095 -0.369  -0.722  -0.242 -0.457 -0.013  0.021  0.261    
##             ET21.5:AT1 ET17.5:AT2
## ExpTemp17.5                      
## ExpTemp21.5                      
## AccTemp17.5                      
## AccTemp21.5                      
## SpcsSgmrssv                      
## Weight_g                         
## ET17.5:AT17                      
## ET21.5:AT17                      
## ET17.5:AT21  0.261               
## ET21.5:AT21  0.516      0.506
  # Compare model with and without interaction with species
    anova(final_model, final_model2)
## refitting model(s) with ML (instead of REML)
## Data: sub_trait
## Models:
## final_model: value ~ ExpTemp + AccTemp + Species + Weight_g + (1 | AntTag) + ExpTemp:AccTemp
## final_model2: value ~ ExpTemp * AccTemp * Species + Weight_g + (1 | AntTag)
##              npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## final_model    13 890.50 925.37 -432.25   864.50                     
## final_model2   21 896.19 952.52 -427.10   854.19 10.308  8     0.2441
    ## --> model with species interaction quite similar and will therefore retain interaction. Contrast analysis also showed species dependent effects of acclimation (see below)
    # Replace final model
      final_model <- final_model2
      
  # Model summary
    #summary(final_model)

    # 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                 12628.9  6314.4     2 58.485 59.7084 7.446e-15 ***
## AccTemp                   655.7   327.8     2 31.272  3.0999   0.05913 .  
## Species                  1008.9  1008.9     1 31.212  9.5402   0.00420 ** 
## Weight_g                  437.1   437.1     1 35.791  4.1333   0.04952 *  
## ExpTemp:AccTemp          1004.3   251.1     4 58.539  2.3742   0.06241 .  
## ExpTemp:Species           527.2   263.6     2 58.851  2.4927   0.09138 .  
## AccTemp:Species           222.9   111.4     2 31.315  1.0538   0.36064    
## ExpTemp:AccTemp:Species   188.6    47.2     4 58.897  0.4459   0.77491    
## ---
## 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) = 59.71, p < .001; Eta2 (partial) = 0.67, 90% CI [0.55, 0.75])
##   - The main effect of AccTemp is statistically not significant and large (F(2) = 3.10, p = 0.059; Eta2 (partial) = 0.17, 90% CI [0.00, 0.34])
##   - The main effect of Species is statistically significant and large (F(1) = 9.54, p = 0.004; Eta2 (partial) = 0.23, 90% CI [0.05, 0.43])
##   - The main effect of Weight_g is statistically significant and medium (F(1) = 4.13, p = 0.050; Eta2 (partial) = 0.10, 90% CI [1.25e-04, 0.28])
##   - The interaction between ExpTemp and AccTemp is statistically not significant and medium (F(4) = 2.37, p = 0.062; Eta2 (partial) = 0.14, 90% CI [0.00, 0.27])
##   - The interaction between ExpTemp and Species is statistically not significant and medium (F(2) = 2.49, p = 0.091; Eta2 (partial) = 0.08, 90% CI [0.00, 0.19])
##   - The interaction between AccTemp and Species is statistically not significant and medium (F(2) = 1.05, p = 0.361; Eta2 (partial) = 0.06, 90% CI [0.00, 0.21])
##   - The interaction between ExpTemp, AccTemp and Species is statistically not significant and small (F(4) = 0.45, p = 0.775; Eta2 (partial) = 0.03, 90% CI [0.00, 0.07])
## 
## Effect sizes were labelled following Field's (2013) recommendations.
  # Table output from sjPlot
    tab_model(final_model, digits.p = 3)
  value
Predictors Estimates CI p
(Intercept) 135.95 108.73 – 163.17 <0.001
ExpTemp1 -13.51 -16.32 – -10.70 <0.001
ExpTemp2 -0.09 -2.85 – 2.68 0.952
AccTemp1 4.03 -2.36 – 10.42 0.216
AccTemp2 4.31 -2.13 – 10.76 0.189
Species1 -7.27 -11.88 – -2.66 0.002
Weight_g -0.03 -0.05 – -0.00 0.042
ExpTemp1 * AccTemp1 1.47 -2.53 – 5.47 0.471
ExpTemp2 * AccTemp1 2.77 -1.16 – 6.69 0.168
ExpTemp1 * AccTemp2 3.11 -0.87 – 7.10 0.126
ExpTemp2 * AccTemp2 -3.61 -7.52 – 0.30 0.071
ExpTemp1 * Species1 2.98 0.17 – 5.79 0.038
ExpTemp2 * Species1 -0.47 -3.25 – 2.30 0.739
AccTemp1 * Species1 4.25 -2.17 – 10.66 0.195
AccTemp2 * Species1 -0.16 -6.61 – 6.28 0.961
ExpTemp1 * AccTemp1 *
Species1
-1.00 -5.01 – 3.02 0.626
ExpTemp2 * AccTemp1 *
Species1
0.57 -3.37 – 4.50 0.777
ExpTemp1 * AccTemp2 *
Species1
1.65 -2.35 – 5.64 0.419
ExpTemp2 * AccTemp2 *
Species1
-2.44 -6.38 – 1.49 0.223
Random Effects
σ2 105.75
τ00 AntTag 169.32
ICC 0.62
N AntTag 39
Observations 108
Marginal R2 / Conditional R2 0.470 / 0.796
  # 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, Species and Weight_g (formula: value ~ ExpTemp * AccTemp * Species + Weight_g). The model included AntTag as random effect (formula: ~1 | AntTag). The model's total explanatory power is substantial (conditional R2 = 0.80) and the part related to the fixed effects alone (marginal R2) is of 0.47. The model's intercept, corresponding to ExpTemp = 14, AccTemp = 14, Species = Jasus edwardsii and Weight_g = 0, is at 135.95 (95% CI [108.73, 163.17], t(87) = 9.79, p < .001). Within this model:
## 
##   - The effect of ExpTemp1 is statistically significant and negative (beta = -13.51, 95% CI [-16.32, -10.70], t(87) = -9.43, p < .001; Std. beta = -0.61, 95% CI [-0.74, -0.49])
##   - The effect of ExpTemp2 is statistically non-significant and negative (beta = -0.09, 95% CI [-2.85, 2.68], t(87) = -0.06, p = 0.952; Std. beta = -3.89e-03, 95% CI [-0.13, 0.12])
##   - The effect of AccTemp1 is statistically non-significant and positive (beta = 4.03, 95% CI [-2.36, 10.42], t(87) = 1.24, p = 0.216; Std. beta = 0.18, 95% CI [-0.11, 0.47])
##   - The effect of AccTemp2 is statistically non-significant and positive (beta = 4.31, 95% CI [-2.13, 10.76], t(87) = 1.31, p = 0.189; Std. beta = 0.20, 95% CI [-0.10, 0.49])
##   - The effect of Species1 is statistically significant and negative (beta = -7.27, 95% CI [-11.88, -2.66], t(87) = -3.09, p = 0.002; Std. beta = -0.33, 95% CI [-0.54, -0.12])
##   - The effect of Weight_g is statistically significant and negative (beta = -0.03, 95% CI [-0.05, -9.03e-04], t(87) = -2.03, p = 0.042; Std. beta = -0.22, 95% CI [-0.43, -7.90e-03])
##   - The interaction effect of AccTemp1 on ExpTemp1 is statistically non-significant and positive (beta = 1.47, 95% CI [-2.53, 5.47], t(87) = 0.72, p = 0.471; Std. beta = 0.07, 95% CI [-0.11, 0.25])
##   - The interaction effect of AccTemp1 on ExpTemp2 is statistically non-significant and positive (beta = 2.77, 95% CI [-1.16, 6.69], t(87) = 1.38, p = 0.168; Std. beta = 0.13, 95% CI [-0.05, 0.30])
##   - The interaction effect of AccTemp2 on ExpTemp1 is statistically non-significant and positive (beta = 3.11, 95% CI [-0.87, 7.10], t(87) = 1.53, p = 0.126; Std. beta = 0.14, 95% CI [-0.04, 0.32])
##   - The interaction effect of AccTemp2 on ExpTemp2 is statistically non-significant and negative (beta = -3.61, 95% CI [-7.52, 0.30], t(87) = -1.81, p = 0.071; Std. beta = -0.16, 95% CI [-0.34, 0.01])
##   - The interaction effect of Species1 on ExpTemp1 is statistically significant and positive (beta = 2.98, 95% CI [0.17, 5.79], t(87) = 2.08, p = 0.038; Std. beta = 0.14, 95% CI [7.62e-03, 0.26])
##   - The interaction effect of Species1 on ExpTemp2 is statistically non-significant and negative (beta = -0.47, 95% CI [-3.25, 2.30], t(87) = -0.33, p = 0.739; Std. beta = -0.02, 95% CI [-0.15, 0.10])
##   - The interaction effect of Species1 on AccTemp1 is statistically non-significant and positive (beta = 4.25, 95% CI [-2.17, 10.66], t(87) = 1.30, p = 0.195; Std. beta = 0.19, 95% CI [-0.10, 0.48])
##   - The interaction effect of Species1 on AccTemp2 is statistically non-significant and negative (beta = -0.16, 95% CI [-6.61, 6.28], t(87) = -0.05, p = 0.961; Std. beta = -7.36e-03, 95% CI [-0.30, 0.28])
##   - The interaction effect of Species1 on ExpTemp1 * AccTemp1 is statistically non-significant and negative (beta = -1.00, 95% CI [-5.01, 3.02], t(87) = -0.49, p = 0.626; Std. beta = -0.05, 95% CI [-0.23, 0.14])
##   - The interaction effect of Species1 on ExpTemp2 * AccTemp1 is statistically non-significant and positive (beta = 0.57, 95% CI [-3.37, 4.50], t(87) = 0.28, p = 0.777; Std. beta = 0.03, 95% CI [-0.15, 0.20])
##   - The interaction effect of Species1 on ExpTemp1 * AccTemp2 is statistically non-significant and positive (beta = 1.65, 95% CI [-2.35, 5.64], t(87) = 0.81, p = 0.419; Std. beta = 0.07, 95% CI [-0.11, 0.26])
##   - The interaction effect of Species1 on ExpTemp2 * AccTemp2 is statistically non-significant and negative (beta = -2.44, 95% CI [-6.38, 1.49], t(87) = -1.22, p = 0.223; Std. beta = -0.11, 95% CI [-0.29, 0.07])
## 
## 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 12628.9 6314.4 58.5 59.71 (2) 7.446e-15 0.67 [0.55-0.75]
AccTemp 655.7 327.8 31.3 3.1 (2) 5.913e-02 0.17 [0-0.34]
Species 1008.9 1008.9 31.2 9.54 (1) 4.200e-03 0.23 [0.05-0.43]
Weight_g 437.1 437.1 35.8 4.13 (1) 4.952e-02 0.1 [0-0.28]
ExpTemp:AccTemp 1004.3 251.1 58.5 2.37 (4) 6.241e-02 0.14 [0-0.27]
ExpTemp:Species 527.2 263.6 58.9 2.49 (2) 9.138e-02 0.08 [0-0.19]
AccTemp:Species 222.9 111.4 31.3 1.05 (2) 3.606e-01 0.06 [0-0.21]
ExpTemp:AccTemp:Species 188.6 47.2 58.9 0.45 (4) 7.749e-01 0.03 [0-0.07]

–> Experimental temperature species and weight have all significant effects on MMR

–> trend for acclimation temperature

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_model2, 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 strong violations of model assumptions

–> some outliers

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 -0.1587081 -0.0166972 1 No MMR
14 - 21.5 14 Jasus edwardsii 26.4164781 2.7533828 0.024 Yes MMR
17.5 - 21.5 14 Jasus edwardsii 26.5751862 2.7690258 0.023 Yes MMR
14 - 17.5 17.5 Jasus edwardsii 13.5136817 1.4264036 0.478 No MMR
14 - 21.5 17.5 Jasus edwardsii 21.3262011 2.2253668 0.091 No MMR
17.5 - 21.5 17.5 Jasus edwardsii 7.8125194 0.8119432 1 No MMR
14 - 17.5 21.5 Jasus edwardsii -0.9730142 -0.1029332 1 No MMR
14 - 21.5 21.5 Jasus edwardsii 14.3868949 1.5113248 0.41 No MMR
17.5 - 21.5 21.5 Jasus edwardsii 15.3599092 1.6011369 0.346 No MMR
14 - 17.5 14 Sagmariasus verreauxi -3.6852236 -0.3968126 1 No MMR
14 - 21.5 14 Sagmariasus verreauxi 10.4543134 1.1089070 0.817 No MMR
17.5 - 21.5 14 Sagmariasus verreauxi 14.1395370 1.5003347 0.418 No MMR
14 - 17.5 17.5 Sagmariasus verreauxi -1.3256659 -0.1420389 1 No MMR
14 - 21.5 17.5 Sagmariasus verreauxi 7.2819029 0.7702734 1 No MMR
17.5 - 21.5 17.5 Sagmariasus verreauxi 8.6075688 0.9104056 1 No MMR
14 - 17.5 21.5 Sagmariasus verreauxi -9.0530080 -0.9696316 1 No MMR
14 - 21.5 21.5 Sagmariasus verreauxi -5.5780971 -0.5900463 1 No MMR
17.5 - 21.5 21.5 Sagmariasus verreauxi 3.4749110 0.3671210 1 No MMR

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", "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.8350000 -2.1617550 0.105 No MMR
14 - 21.5 14 Jasus edwardsii -17.3390272 -2.7825252 0.021 Yes MMR
17.5 - 21.5 14 Jasus edwardsii -4.5040272 -0.7227954 1 No MMR
14 - 17.5 17.5 Jasus edwardsii 0.8373898 0.1407746 1 No MMR
14 - 21.5 17.5 Jasus edwardsii -18.1533333 -3.0575036 0.01 Yes MMR
17.5 - 21.5 17.5 Jasus edwardsii -18.9907231 -3.1925530 0.007 Yes MMR
14 - 17.5 21.5 Jasus edwardsii -17.9252770 -3.0189098 0.011 Yes MMR
14 - 21.5 21.5 Jasus edwardsii -29.3686103 -4.9461543 0 Yes MMR
17.5 - 21.5 21.5 Jasus edwardsii -11.4433333 -1.9273614 0.177 No MMR
14 - 17.5 14 Sagmariasus verreauxi -16.6007584 -2.6631218 0.03 Yes MMR
14 - 21.5 14 Sagmariasus verreauxi -25.4524250 -4.0831212 0 Yes MMR
17.5 - 21.5 14 Sagmariasus verreauxi -8.8516667 -1.4908558 0.425 No MMR
14 - 17.5 17.5 Sagmariasus verreauxi -14.2412006 -2.2708359 0.08 No MMR
14 - 21.5 17.5 Sagmariasus verreauxi -30.8202095 -4.9477974 0 Yes MMR
17.5 - 21.5 17.5 Sagmariasus verreauxi -16.5790088 -2.7756200 0.022 Yes MMR
14 - 17.5 21.5 Sagmariasus verreauxi -19.7731688 -3.3273515 0.005 Yes MMR
14 - 21.5 21.5 Sagmariasus verreauxi -41.4848355 -6.9809057 0 Yes MMR
17.5 - 21.5 21.5 Sagmariasus verreauxi -21.7116667 -3.6568215 0.002 Yes MMR

9.C 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")
## 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 -8.580352 -1.556536 0.125 No MMR
Jasus edwardsii - Sagmariasus verreauxi 17.5 -15.477765 -2.825278 0.007 Yes MMR
Jasus edwardsii - Sagmariasus verreauxi 21.5 -19.545851 -3.539815 0.001 Yes MMR

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)