1. Load libraries

library(tidyverse)
library(readxl)
library(dplyr)
library(plyr)
library(openxlsx)
library(psycho)
library(pander)
library(ggpubr)
library(knitr)
library(kableExtra) # For Markdown tables
library(report) # textual summaries of analysis

library(fitdistrplus) # to fit probability distributions
library(emmeans)
library(glmmTMB) # for GLMM
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 <- "SMR"

 # 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 and random 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 response variable

–> to check normality assumptions for the LMM or choose the best fitting distribution for the GLMM

  # Plot data versus a range of distrubutions
    # Skewness-kurtosis plot for continuous data
    descdist(sub_trait$value, boot = 1000)

## summary statistics
## ------
## min:  5.58   max:  57.79 
## median:  23.07 
## mean:  25.6712 
## estimated sd:  11.90557 
## estimated skewness:  0.7880276 
## estimated kurtosis:  3.038012

  # Fit various distributions
    nd <- fitdist(sub_trait$value, "norm")

## $start.arg
## $start.arg$mean
## [1] 25.6712
## 
## $start.arg$sd
## [1] 11.85032
## 
## 
## $fix.arg
## NULL

    gd <- fitdist(sub_trait$value, "gamma")

## $start.arg
## $start.arg$shape
## [1] 4.692801
## 
## $start.arg$rate
## [1] 0.1828041
## 
## 
## $fix.arg
## NULL

    wd <- fitdist(sub_trait$value, "weibull")

## $start.arg
## $start.arg$shape
## [1] 2.519369
## 
## $start.arg$scale
## [1] 28.9234
## 
## 
## $fix.arg
## NULL

    lnd <- fitdist(sub_trait$value, "lnorm")

## $start.arg
## $start.arg$meanlog
## [1] 3.13761
## 
## $start.arg$sdlog
## [1] 0.4740995
## 
## 
## $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
    g <- gofstat(list(nd, gd, wd, lnd), fitnames = c("Normal", "gamma", "Weibull", "LogNormal"))
    g

## Goodness-of-fit statistics
##                                 Normal      gamma    Weibull  LogNormal
## Kolmogorov-Smirnov statistic 0.1149477 0.06536997 0.08557554 0.06042317
## Cramer-von Mises statistic   0.3310982 0.09988048 0.17101880 0.06072103
## Anderson-Darling statistic   2.0129375 0.58101517 1.07488632 0.41214586
## 
## Goodness-of-fit criteria
##                                  Normal    gamma  Weibull LogNormal
## Akaike's Information Criterion 844.5194 826.2741 832.4556  827.0055
## Bayesian Information Criterion 849.8837 831.6383 837.8199  832.3697

–> Anderson-Darling statistics gives good measure for fit for both middle and tail of distribution

–> SMR data are not normal distributed and fit best to a log normal or gamma distribution **

–> use of GLMM as modelling approach for non-normal data

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

Grouped 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()

SMR
Species	mean	n	sd	ci	lower_ci	upper_ci	min	max
Jasus edwardsii	26.0	54	12.0	3.3	22.7	29.3	5.6	52.1
Sagmariasus verreauxi	25.3	54	11.9	3.3	22.1	28.6	9.2	57.8

Grouped 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()

SMR
Species	ExpTemp	mean	n	sd	ci	lower_ci	upper_ci	min	max
Jasus edwardsii	14	15.9	18	4.9	2.4	13.5	18.3	5.6	23.2
Jasus edwardsii	17.5	24.9	18	8.7	4.3	20.6	29.2	12.4	40.6
Jasus edwardsii	21.5	37.2	18	10.3	5.1	32.1	42.3	20.7	52.1
Sagmariasus verreauxi	14	15.0	18	3.6	1.8	13.2	16.7	9.2	21.1
Sagmariasus verreauxi	17.5	23.7	18	7.5	3.7	20.0	27.4	14.4	37.7
Sagmariasus verreauxi	21.5	37.3	18	10.4	5.1	32.1	42.4	24.8	57.8

Grouped by species, experimental temperature, and acclimation 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()

SMR
Species	ExpTemp	AccTemp	mean	n	sd	ci	lower_ci	upper_ci	min	max
Jasus edwardsii	14	14	18.1	6	3.2	3.2	15.0	21.3	13.6	23.2
Jasus edwardsii	14	17.5	15.6	6	3.5	3.5	12.2	19.1	10.9	20.0
Jasus edwardsii	14	21.5	13.9	6	7.0	7.0	6.9	20.9	5.6	21.7
Jasus edwardsii	17.5	14	28.1	6	8.3	8.3	19.8	36.4	15.4	36.6
Jasus edwardsii	17.5	17.5	25.3	6	10.8	10.8	14.5	36.1	15.5	40.6
Jasus edwardsii	17.5	21.5	21.4	6	6.6	6.6	14.8	27.9	12.4	27.8
Jasus edwardsii	21.5	14	44.3	6	6.2	6.2	38.1	50.5	34.3	51.4
Jasus edwardsii	21.5	17.5	38.4	6	10.9	10.9	27.5	49.3	25.4	52.1
Jasus edwardsii	21.5	21.5	29.0	6	7.5	7.5	21.4	36.5	20.7	41.3
Sagmariasus verreauxi	14	14	17.1	6	3.1	3.1	14.0	20.1	12.0	20.9
Sagmariasus verreauxi	14	17.5	15.4	6	4.3	4.3	11.0	19.7	9.2	21.1
Sagmariasus verreauxi	14	21.5	12.5	6	1.5	1.5	11.0	14.0	10.3	14.8
Sagmariasus verreauxi	17.5	14	30.7	6	4.0	4.0	26.8	34.7	25.1	36.2
Sagmariasus verreauxi	17.5	17.5	23.2	6	7.3	7.3	15.9	30.5	18.6	37.7
Sagmariasus verreauxi	17.5	21.5	17.2	6	3.2	3.2	14.0	20.4	14.4	23.0
Sagmariasus verreauxi	21.5	14	45.9	6	12.5	12.5	33.4	58.3	30.3	57.8
Sagmariasus verreauxi	21.5	17.5	36.0	6	6.3	6.3	29.7	42.2	26.8	46.0
Sagmariasus verreauxi	21.5	21.5	30.0	6	4.2	4.2	25.8	34.2	24.8	34.8

7. Generalised mixed effect model

# Instructions to use glmm package followed by https://rdrr.io/cran/glmm/f/inst/doc/glmm.pdf

# Fit GLMM to progressively more complex fixed effect terms
# Animal ID as random factor

  # Fixed Effects: Experimental and Acclimation temperature 
    glmm1 <- glmmTMB(value ~ ExpTemp + AccTemp + (1| AntTag),
                    data=sub_trait,
                    family=Gamma(link = "identity"))

  # Fixed Effects: Experimental and Acclimation temperature + interaction
    glmm1b <- glmmTMB(value ~ ExpTemp * AccTemp + (1| AntTag),
                    data=sub_trait,
                    family=Gamma(link = "identity"))

  # Fixed Effects: Experimental and Acclimation temperature and species 
    glmm2 <- glmmTMB(value ~ ExpTemp + AccTemp + Species + (1|AntTag),
                    data=sub_trait,
                    family=Gamma(link = "identity"))

  # Fixed Effects: Experimental and Acclimation temperature (+ interaction) and species 
    glmm2b <- glmmTMB(value ~ ExpTemp * AccTemp + Species + (1|AntTag),
                    data=sub_trait,
                    family=Gamma(link = "identity"))

  # Fixed Effects: Experimental and Acclimation temperature and species + interaction between all terms
    glmm2c <- glmmTMB(value ~ ExpTemp * AccTemp * Species + (1|AntTag),
                data=sub_trait,
                family=Gamma(link = "identity"))
    
  # Fixed Effects: Experimental and Acclimation temperature and species + interaction between all terms + Sex
    glmm3 <- glmmTMB(value ~ ExpTemp + AccTemp + Sex + (1|AntTag),
                data=sub_trait,
                family=Gamma(link = "identity"))
    
  # Fixed Effects: Experimental and Acclimation temperature and species + interaction between all terms + Sex
    glmm3b <- glmmTMB(value ~ ExpTemp * AccTemp + Sex + (1|AntTag),
                data=sub_trait,
                family=Gamma(link = "identity"))
    
  # Fixed Effects: Experimental and Acclimation temperature and species + interaction between all terms + Sex + Weight_g
    glmm4 <- glmmTMB(value ~ ExpTemp + AccTemp + Weight_g + (1|AntTag),
                data=sub_trait,
                family=Gamma(link = "identity"))
    
  # Fixed Effects: Experimental and Acclimation temperature and species + interaction between all terms + Sex + Weight_g
    glmm4b <- glmmTMB(value ~ ExpTemp * AccTemp + Weight_g + (1|AntTag),
                data=sub_trait,
                family=Gamma(link = "identity"))

# Compare model performances using AIC
  AIC(glmm1,glmm1b,glmm2, glmm2b,glmm2c, glmm3, glmm3b,glmm4,glmm4b)

##        df      AIC
## glmm1   7 702.0261
## glmm1b 11 699.2096
## glmm2   8 703.8759
## glmm2b 12 700.9846
## glmm2c 20 712.3681
## glmm3   8 700.2934
## glmm3b 12 697.3871
## glmm4   8 703.3554
## glmm4b 12 700.4961

  AICctab(glmm1,glmm1b,glmm2, glmm2b,glmm2c, glmm3, glmm3b,glmm4,glmm4b, nobs=nrow(sub_trait)) # ranked

##        dAICc df
## glmm3b  0.0  12
## glmm3   1.1  8 
## glmm1b  1.3  11
## glmm1   2.5  7 
## glmm4b  3.1  12
## glmm2b  3.6  12
## glmm4   4.1  8 
## glmm2   4.7  8 
## glmm2c 21.4  20

–> Model with experimental, acclimation temperature and sex as fixed factors without interaction performed best

–> Species as fixed factor did not improve model

–> Best model: ExpTemp * AccTemp + Sex + (1|AntTag)

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 = glmm3b, 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$Sex)

–> Model with gamma distribution and “idendity” link function provides a good fit

–> no overdispersion

9. Model summary

  # Model summary
    summary(glmm3b)

##  Family: Gamma  ( identity )
## Formula:          value ~ ExpTemp * AccTemp + Sex + (1 | AntTag)
## Data: sub_trait
## 
##      AIC      BIC   logLik deviance df.resid 
##    697.4    729.6   -336.7    673.4       96 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  AntTag (Intercept) 6.96     2.638   
## Number of obs: 108, groups:  AntTag, 39
## 
## Dispersion estimate for Gamma family (sigma^2): 0.0429 
## 
## Conditional model:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               19.806      1.617  12.245  < 2e-16 ***
## ExpTemp17.5               11.193      2.065   5.420 5.95e-08 ***
## ExpTemp21.5               26.972      2.913   9.260  < 2e-16 ***
## AccTemp17.5               -2.261      1.761  -1.284  0.19912    
## AccTemp21.5               -5.015      1.690  -2.967  0.00301 ** 
## SexMale                   -2.619      1.325  -1.978  0.04797 *  
## ExpTemp17.5:AccTemp17.5   -3.124      2.671  -1.169  0.24228    
## ExpTemp21.5:AccTemp17.5   -5.716      3.774  -1.515  0.12989    
## ExpTemp17.5:AccTemp21.5   -4.986      2.477  -2.013  0.04408 *  
## ExpTemp21.5:AccTemp21.5  -10.478      3.469  -3.020  0.00253 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  # Plot estimates and confidence intervals
    plot_model(glmm3b, 
                 type = "est",
                 sort.est = TRUE,
                 show.values = TRUE)

  # Table output from sjPlot
    tab_model(glmm3b)

	value
Predictors	Estimates	CI	p
(Intercept)	19.81	16.64 – 22.98	<0.001
ExpTemp [17.5]	11.19	7.15 – 15.24	<0.001
ExpTemp [21.5]	26.97	21.26 – 32.68	<0.001
AccTemp [17.5]	-2.26	-5.71 – 1.19	0.199
AccTemp [21.5]	-5.01	-8.33 – -1.70	0.003
Sex [Male]	-2.62	-5.22 – -0.02	0.048
ExpTemp [17.5] * AccTemp [17.5]	-3.12	-8.36 – 2.11	0.242
ExpTemp [21.5] * AccTemp [17.5]	-5.72	-13.11 – 1.68	0.130
ExpTemp [17.5] * AccTemp [21.5]	-4.99	-9.84 – -0.13	0.044
ExpTemp [21.5] * AccTemp [21.5]	-10.48	-17.28 – -3.68	0.003
Random Effects
σ²	0.04
τ₀₀ _AntTag	6.96
ICC	0.99
N _AntTag	39
Observations	108
Marginal R² / Conditional R²	0.935 / 1.000

  # Textual report of the model
    report(glmm3b)

## We fitted a general linear mixed model (Gamma family with a identity link) (estimated using ML and nlminb optimizer) to predict value with ExpTemp, AccTemp and Sex (formula: value ~ ExpTemp * AccTemp + Sex). The model included AntTag as random effect (formula: ~1 | AntTag). The model's total explanatory power is substantial (conditional R2 = 1.00) and the part related to the fixed effects alone (marginal R2) is of 0.94. The model's intercept, corresponding to ExpTemp = 14, AccTemp = 14 and Sex = Female, is at 19.81 (95% CI [16.64, 22.98], p < .001). Within this model:
## 
##   - The effect of ExpTemp [17.5] is statistically significant and positive (beta = 11.19, 95% CI [7.15, 15.24], p < .001; Std. beta = 11.19, 95% CI [7.15, 15.24])
##   - The effect of ExpTemp [21.5] is statistically significant and positive (beta = 26.97, 95% CI [21.26, 32.68], p < .001; Std. beta = 26.97, 95% CI [21.26, 32.68])
##   - The effect of AccTemp [17.5] is statistically non-significant and negative (beta = -2.26, 95% CI [-5.71, 1.19], p = 0.199; Std. beta = -2.26, 95% CI [-5.71, 1.19])
##   - The effect of AccTemp [21.5] is statistically significant and negative (beta = -5.01, 95% CI [-8.33, -1.70], p = 0.003; Std. beta = -5.01, 95% CI [-8.33, -1.70])
##   - The effect of Sex [Male] is statistically significant and negative (beta = -2.62, 95% CI [-5.22, -0.02], p = 0.048; Std. beta = -2.62, 95% CI [-5.22, -0.02])
##   - The interaction effect of AccTemp [17.5] on ExpTemp [17.5] is statistically non-significant and negative (beta = -3.12, 95% CI [-8.36, 2.11], p = 0.242; Std. beta = -3.12, 95% CI [-8.36, 2.11])
##   - The interaction effect of AccTemp [17.5] on ExpTemp [21.5] is statistically non-significant and negative (beta = -5.72, 95% CI [-13.11, 1.68], p = 0.130; Std. beta = -5.72, 95% CI [-13.11, 1.68])
##   - The interaction effect of AccTemp [21.5] on ExpTemp [17.5] is statistically significant and negative (beta = -4.99, 95% CI [-9.84, -0.13], p = 0.044; Std. beta = -4.99, 95% CI [-9.84, -0.13])
##   - The interaction effect of AccTemp [21.5] on ExpTemp [21.5] is statistically significant and negative (beta = -10.48, 95% CI [-17.28, -3.68], p = 0.003; Std. beta = -10.48, 95% CI [-17.28, -3.68])
## 
## Standardized parameters were obtained by fitting the model on a standardized version of the dataset.

–> SMR increases significantly with experimental temperature irrespective of species and acclimation temperature

–> SMR is significantly lower at 21.5°C and 17.5°C compared to 14°C irrespective of species and experimental temperature

10. Contrasts between levels

# Display interaction plots
  emmip(glmm3b, AccTemp ~ ExpTemp | Sex)

## Between acclimation temperature

# Calculate contrasts
  contrasts1 <- emmeans(glmm3b, ~ AccTemp | ExpTemp, component = "cond", type = "response", adjust = "bonferroni")

  contrasts_df <- contrasts1 %>% 
                  summary(infer = TRUE) %>% 
                  as.data.frame() 

# Plot contrasts
  plot(contrasts1, comparisons = TRUE)

# Save contrasts in table
  # Add trait columns
    contrasts_df$Trait <- trait

  # Save contrasts results
    write.csv(contrasts_df, paste("Contrasts_AcclimationTemperature_", trait, ".csv", sep = ""))
      
# Display as HTML table colouring significant p values red
  contrasts_df   %>%   kable(escape = FALSE)  %>%
                kable_styling(bootstrap_options = c('hover'))

AccTemp	ExpTemp	emmean	SE	df	lower.CL	upper.CL	t.ratio	Trait
14	14	18.49663	1.338843	96	15.23433	21.75892	13.81539	SMR
17.5	14	16.23560	1.236852	96	13.22182	19.24938	13.12656	SMR
21.5	14	13.48194	1.122152	96	10.74764	16.21623	12.01436	SMR
14	17.5	29.68958	1.909179	96	25.03757	34.34159	15.55097	SMR
17.5	17.5	24.30480	1.614343	96	20.37120	28.23839	15.05554	SMR
21.5	17.5	19.68852	1.381729	96	16.32172	23.05531	14.24919	SMR
14	21.5	45.46888	2.794752	96	38.65903	52.27872	16.26938	SMR
17.5	21.5	37.49196	2.344896	96	31.77827	43.20566	15.98876	SMR
21.5	21.5	29.97639	1.918562	96	25.30152	34.65126	15.62441	SMR

## Between experimental temperature
  
# Calculate contrasts
  contrasts1 <- emmeans(glmm3b, ~ ExpTemp | AccTemp, component = "cond", type = "response", adjust = "bonferroni")
  
  contrasts_df <- contrasts1 %>% 
                  summary(infer = TRUE) %>% 
                  as.data.frame() 

# Plot contrasts
  plot(contrasts1, comparisons = TRUE)

# Save contrasts in table
  # Add trait columns
    contrasts_df$Trait <- trait

  # Save contrasts results
    write.csv(contrasts_df, paste("Contrasts_ExperimentalnTemperature_", trait, ".csv", sep = ""))
      
# Display as HTML table colouring significant p values red
  contrasts_df   %>%   kable(escape = FALSE)  %>%
                kable_styling(bootstrap_options = c('hover'))

ExpTemp	AccTemp	emmean	SE	df	lower.CL	upper.CL	t.ratio	Trait
14	14	18.49663	1.338843	96	15.23433	21.75892	13.81539	SMR
17.5	14	29.68958	1.909179	96	25.03757	34.34159	15.55097	SMR
21.5	14	45.46888	2.794752	96	38.65903	52.27872	16.26938	SMR
14	17.5	16.23560	1.236852	96	13.22182	19.24938	13.12656	SMR
17.5	17.5	24.30480	1.614343	96	20.37120	28.23839	15.05554	SMR
21.5	17.5	37.49196	2.344896	96	31.77827	43.20566	15.98876	SMR
14	21.5	13.48194	1.122152	96	10.74764	16.21623	12.01436	SMR
17.5	21.5	19.68852	1.381729	96	16.32172	23.05531	14.24919	SMR
21.5	21.5	29.97639	1.918562	96	25.30152	34.65126	15.62441	SMR

11. 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)
##   - bbmle (version 1.0.24; Ben Bolker and R Development Core Team, 2021)
##   - 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)
##   - 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)
##   - glmmTMB (version 1.1.2.3; Mollie Brooks et al., 2017)
##   - 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)
##   - 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)

Effects of acclimation temperature, experimental temperatures on SMR

Michael Oellermann

07 December, 2021

1. Load libraries

2. Load data

3. Process data

4. Find best probability distribution of response variable

–> SMR data are not normal distributed and fit best to a log normal or gamma distribution **

–> use of GLMM as modelling approach for non-normal data

5. Plot data

6. Summary statistics

Grouped by species

Grouped by species and experimental temperature

Grouped by species, experimental temperature, and acclimation temperature

7. Generalised mixed effect model

–> Model with experimental, acclimation temperature and sex as fixed factors without interaction performed best

–> Species as fixed factor did not improve model

–> Best model: ExpTemp * AccTemp + Sex + (1|AntTag)

8. Model diagnostics

–> Model with gamma distribution and “idendity” link function provides a good fit

–> no overdispersion

9. Model summary

–> SMR increases significantly with experimental temperature irrespective of species and acclimation temperature

–> SMR is significantly lower at 21.5°C and 17.5°C compared to 14°C irrespective of species and experimental temperature

10. Contrasts between levels

11. R Session Information