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)

–> 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
σ2 0.04
τ00 AntTag 6.96
ICC 0.99
N AntTag 39
Observations 108
Marginal R2 / Conditional R2 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 p.value Trait
14 14 18.49663 1.338843 96 15.23433 21.75892 13.81539 0 SMR
17.5 14 16.23560 1.236852 96 13.22182 19.24938 13.12656 0 SMR
21.5 14 13.48194 1.122152 96 10.74764 16.21623 12.01436 0 SMR
14 17.5 29.68958 1.909179 96 25.03757 34.34159 15.55097 0 SMR
17.5 17.5 24.30480 1.614343 96 20.37120 28.23839 15.05554 0 SMR
21.5 17.5 19.68852 1.381729 96 16.32172 23.05531 14.24919 0 SMR
14 21.5 45.46888 2.794752 96 38.65903 52.27872 16.26938 0 SMR
17.5 21.5 37.49196 2.344896 96 31.77827 43.20566 15.98876 0 SMR
21.5 21.5 29.97639 1.918562 96 25.30152 34.65126 15.62441 0 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 p.value Trait
14 14 18.49663 1.338843 96 15.23433 21.75892 13.81539 0 SMR
17.5 14 29.68958 1.909179 96 25.03757 34.34159 15.55097 0 SMR
21.5 14 45.46888 2.794752 96 38.65903 52.27872 16.26938 0 SMR
14 17.5 16.23560 1.236852 96 13.22182 19.24938 13.12656 0 SMR
17.5 17.5 24.30480 1.614343 96 20.37120 28.23839 15.05554 0 SMR
21.5 17.5 37.49196 2.344896 96 31.77827 43.20566 15.98876 0 SMR
14 21.5 13.48194 1.122152 96 10.74764 16.21623 12.01436 0 SMR
17.5 21.5 19.68852 1.381729 96 16.32172 23.05531 14.24919 0 SMR
21.5 21.5 29.97639 1.918562 96 25.30152 34.65126 15.62441 0 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)