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

# 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)

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:  9   max:  169 
## median:  79 
## mean:  82.50962 
## estimated sd:  38.24423 
## estimated skewness:  0.1752918 
## estimated kurtosis:  2.322914
  # Fit various distributions
    nd <- fitdist(sub_trait$value, "norm")
## $start.arg
## $start.arg$mean
## [1] 82.50962
## 
## $start.arg$sd
## [1] 38.05992
## 
## 
## $fix.arg
## NULL
    gd <- fitdist(sub_trait$value, "gamma")
## $start.arg
## $start.arg$shape
## [1] 4.699735
## 
## $start.arg$rate
## [1] 0.05695984
## 
## 
## $fix.arg
## NULL
    wd <- fitdist(sub_trait$value, "weibull")
## $start.arg
## $start.arg$shape
## [1] 2.07303
## 
## $start.arg$scale
## [1] 94.73581
## 
## 
## $fix.arg
## NULL
    lnd <- fitdist(sub_trait$value, "lnorm")
## $start.arg
## $start.arg$meanlog
## [1] 4.275167
## 
## $start.arg$sdlog
## [1] 0.5760732
## 
## 
## $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.07070856 0.09984929 0.08395626 0.1031488
## Cramer-von Mises statistic   0.07139295 0.16226300 0.06405192 0.3380288
## Anderson-Darling statistic   0.46821302 1.00919041 0.41801271 2.0888024
## 
## Goodness-of-fit criteria
##                                  Normal    gamma  Weibull LogNormal
## Akaike's Information Criterion 1056.085 1058.960 1051.292  1073.658
## Bayesian Information Criterion 1061.374 1064.249 1056.581  1078.947

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

–> Data fit best to a 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_total_Acclimation.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)
    
  # 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_total_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_total_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_total_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_total
Species mean n sd ci lower_ci upper_ci min max
Jasus edwardsii 68.0 53 34.3 9.5 58.5 77.4 9 148
Sagmariasus verreauxi 97.6 51 36.5 10.2 87.4 107.9 27 169

Grouped by species, experimental temperature, and acclimation

# 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()
escape_total
Species ExpTemp AccTemp mean n sd ci lower_ci upper_ci min max
Jasus edwardsii 14 14 69.3 6 30.4 30.4 39.0 99.7 24 109
Jasus edwardsii 14 17.5 79.5 6 30.8 30.8 48.7 110.3 55 131
Jasus edwardsii 14 21.5 37.8 6 22.9 22.8 15.0 60.7 19 70
Jasus edwardsii 17.5 14 85.5 6 31.5 31.5 54.0 117.0 46 124
Jasus edwardsii 17.5 17.5 64.7 6 38.7 38.7 26.0 103.4 9 110
Jasus edwardsii 17.5 21.5 69.8 5 47.3 54.3 15.5 124.1 30 148
Jasus edwardsii 21.5 14 62.3 6 42.6 42.6 19.8 104.9 15 117
Jasus edwardsii 21.5 17.5 87.5 6 27.7 27.7 59.8 115.2 59 134
Jasus edwardsii 21.5 21.5 55.5 6 26.1 26.0 29.5 81.5 30 102
Sagmariasus verreauxi 14 14 114.0 6 16.9 16.9 97.1 130.9 89 135
Sagmariasus verreauxi 14 17.5 86.7 6 40.3 40.3 46.4 127.0 27 145
Sagmariasus verreauxi 14 21.5 84.8 5 37.8 43.4 41.4 128.2 43 127
Sagmariasus verreauxi 17.5 14 125.8 5 38.0 43.7 82.1 169.5 73 168
Sagmariasus verreauxi 17.5 17.5 119.8 6 39.7 39.6 80.2 159.5 80 169
Sagmariasus verreauxi 17.5 21.5 86.8 6 35.4 35.4 51.4 122.2 57 150
Sagmariasus verreauxi 21.5 14 82.6 5 43.4 49.9 32.7 132.5 31 130
Sagmariasus verreauxi 21.5 17.5 95.7 6 37.8 37.8 57.9 133.5 29 143
Sagmariasus verreauxi 21.5 21.5 82.5 6 25.5 25.4 57.1 107.9 42 108

Decreased number of total escapes by trend in Jasus edwardsii between 17.5 and 21.5 acclimation temperature at 14 and 21.5 experimental temperature

# At 14°C Experimental temperature
  print(paste("decrease of", round((79.5-37.8)/79.5*100, 0), "% from 17.5 to 21.5°C acclimation temperature at 17.5°C experimental temperature"))
## [1] "decrease of 52 % from 17.5 to 21.5°C acclimation temperature at 17.5°C experimental temperature"
# At 21.5°C Experimental temperature
print(paste("decrease of", round((87.5-55.5)/87.5*100, 0), "% from 17.5 to 21.5°C acclimation temperature at 21.5°C experimental temperature"))
## [1] "decrease of 37 % from 17.5 to 21.5°C acclimation temperature at 21.5°C experimental temperature"

6. Correlations

6.B 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 = "Total escapes", ylab = "EPOC")
## `geom_smooth()` using formula 'y ~ x'

ggsave("Total escapes_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.98224, p-value = 0.6372
    shapiro.test(SRL$value)  
## 
##  Shapiro-Wilk normality test
## 
## data:  SRL$value
## W = 0.96634, p-value = 0.1397
  # Correlation plots
    ggscatter(ERL, x = "value", y = "EPOC", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Total escapes", 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 = "Total escapes", 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.6862, df = 49, p-value = 0.09812
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.04425861  0.47886099
## sample estimates:
##       cor 
## 0.2341814
    cor.test(SRL$value, SRL$EPOC, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  SRL$value and SRL$EPOC
## t = 3.0299, df = 51, p-value = 0.003835
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1344801 0.5977643
## sample estimates:
##       cor 
## 0.3905768

6.B 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 = "Total escapes", 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.98224, p-value = 0.6372
    shapiro.test(SRL$value)  
## 
##  Shapiro-Wilk normality test
## 
## data:  SRL$value
## W = 0.96634, p-value = 0.1397
  # Correlation plots
    ggscatter(ERL, x = "value", y = "MMR", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Total escapes", 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 = "Total escapes", 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.2604, df = 49, p-value = 0.2135
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1034275  0.4317081
## sample estimates:
##       cor 
## 0.1772071
    cor.test(SRL$value, SRL$MMR, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  SRL$value and SRL$MMR
## t = -0.57954, df = 51, p-value = 0.5648
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3436663  0.1936413
## sample estimates:
##         cor 
## -0.08088605

6.C Weight correlations

# Pooled data
ggscatter(sub_trait, x = "value", y = "Weight_g", fill = "Species",color = "Species",
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Total escapes", ylab = "Weight_g")
## `geom_smooth()` using formula 'y ~ x'

ggsave("Escapes_speed_Corelation_Weight_g.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.98224, p-value = 0.6372
    shapiro.test(SRL$value)  
## 
##  Shapiro-Wilk normality test
## 
## data:  SRL$value
## W = 0.96634, p-value = 0.1397
  # Correlation plots
    ggscatter(ERL, x = "value", y = "Weight_g", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Total escapes", ylab = "Weight_g")+
      ggtitle("Sagmariasus verreauxi")
## `geom_smooth()` using formula 'y ~ x'

    ggscatter(SRL, x = "value", y = "Weight_g", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Total escapes", ylab = "Weight_g")+
      ggtitle("Jasus edwardsii")
## `geom_smooth()` using formula 'y ~ x'

  # Correlation tests
    cor.test(ERL$value, ERL$Weight_g, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  ERL$value and ERL$Weight_g
## t = 0.56767, df = 49, p-value = 0.5729
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1991903  0.3486478
## sample estimates:
##        cor 
## 0.08083047
    cor.test(SRL$value, SRL$Weight_g, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  SRL$value and SRL$Weight_g
## t = 4.1276, df = 51, p-value = 0.0001358
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2661080 0.6788791
## sample estimates:
##       cor 
## 0.5004087

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
    anova(complex_model)
## Type III Analysis of Variance Table with Satterthwaite's method
##                         Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## ExpTemp                 3632.5  1816.2     2 65.831  4.7614 0.011714 * 
## AccTemp                 1553.9   776.9     2 35.904  2.0368 0.145247   
## Species                 3938.0  3938.0     1 35.939 10.3237 0.002771 **
## Sex                     1201.9  1201.9     1 38.184  3.1510 0.083853 . 
## Weight_g                1338.4  1338.4     1 40.472  3.5086 0.068285 . 
## ExpTemp:AccTemp         3751.4   937.8     4 65.684  2.4586 0.053996 . 
## ExpTemp:Species         2405.6  1202.8     2 66.285  3.1532 0.049188 * 
## AccTemp:Species          370.4   185.2     2 36.117  0.4855 0.619372   
## ExpTemp:AccTemp:Species 3640.5   910.1     4 66.052  2.3860 0.059942 . 
## ---
## 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 -484.93 1013.9                         
## (1 | AntTag)          0   21 -497.50 1037.0 25.155  1  5.289e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                         Eliminated Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)
## Sex                              1 1201.9  1201.9     1 38.184  3.1510 0.083853
## ExpTemp:AccTemp:Species          2 3660.8   915.2     4 66.193  2.4133 0.057586
## AccTemp:Species                  3  405.9   203.0     2 36.232  0.4777 0.624045
## ExpTemp:AccTemp                  4 3850.0   962.5     4 65.359  2.2604 0.072038
## AccTemp                          5 1694.7   847.3     2 36.143  1.7383 0.190196
## ExpTemp:Species                  6 2442.9  1221.5     2 66.565  2.5203 0.088088
## ExpTemp                          7 3081.4  1540.7     2 66.395  2.9190 0.060951
## Species                          0 5076.7  5076.7     1 36.470  8.9186 0.005021
## Weight_g                         0 2356.4  2356.4     1 39.741  4.1396 0.048601
##                           
## Sex                     . 
## ExpTemp:AccTemp:Species . 
## AccTemp:Species           
## ExpTemp:AccTemp         . 
## AccTemp                   
## ExpTemp:Species         . 
## ExpTemp                 . 
## Species                 **
## Weight_g                * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## value ~ Species + Weight_g + (1 | AntTag)
  # Extract the model that step found:
    final_model <- get_model(step_result)

–> effect of species and weak effect of weight

–> no acclimation effects

7.C Re-run model fit using REML estimation

Best model: Species + Weight_g + (1 | AntTag)

    # Remodel using Maximum likelyhood estimation
      final_model <- lmer(value ~ Species +
                            Weight_g +
                            (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 ~ Species + Weight_g + (1 | AntTag)
##    Data: sub_trait
## 
## REML criterion at convergence: 999.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.93924 -0.57972  0.05684  0.56562  2.17566 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  AntTag   (Intercept) 647.3    25.44   
##  Residual             568.2    23.84   
## Number of obs: 104, groups:  AntTag, 39
## 
## Fixed effects:
##                              Estimate Std. Error       df t value Pr(>|t|)   
## (Intercept)                  16.34591   26.83250 37.26126   0.609   0.5461   
## SpeciesSagmariasus verreauxi 27.59436    9.61042 33.89699   2.871   0.0070 **
## Weight_g                      0.04677    0.02425 37.45211   1.929   0.0614 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) SpcsSv
## SpcsSgmrssv -0.016       
## Weight_g    -0.968 -0.164
    # 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)   
## Species  4684.5  4684.5     1 33.897  8.2443 0.006998 **
## Weight_g 2113.4  2113.4     1 37.452  3.7194 0.061390 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      report(anova(final_model))
## The ANOVA suggests that:
## 
##   - The main effect of Species is statistically significant and large (F(1) = 8.24, p = 0.007; Eta2 (partial) = 0.20, 90% CI [0.04, 0.38])
##   - The main effect of Weight_g is statistically not significant and medium (F(1) = 3.72, p = 0.061; Eta2 (partial) = 0.09, 90% CI [0.00, 0.26])
## 
## Effect sizes were labelled following Field's (2013) recommendations.
    # Table output from sjPlot
      tab_model(final_model)
  value
Predictors Estimates CI p
(Intercept) 16.35 -36.24 – 68.94 0.542
Species [Sagmariasus
verreauxi]
27.59 8.76 – 46.43 0.004
Weight_g 0.05 -0.00 – 0.09 0.054
Random Effects
σ2 568.21
τ00 AntTag 647.28
ICC 0.53
N AntTag 39
Observations 104
Marginal R2 / Conditional R2 0.208 / 0.630
    # Textual report of the model
      report(final_model)
## We fitted a linear mixed model (estimated using REML and nloptwrap optimizer) to predict value with Species and Weight_g (formula: value ~ Species + Weight_g). The model included AntTag as random effect (formula: ~1 | AntTag). The model's total explanatory power is substantial (conditional R2 = 0.63) and the part related to the fixed effects alone (marginal R2) is of 0.21. The model's intercept, corresponding to Species = Jasus edwardsii and Weight_g = 0, is at 16.35 (95% CI [-36.24, 68.94], t(99) = 0.61, p = 0.542). Within this model:
## 
##   - The effect of Species [Sagmariasus verreauxi] is statistically significant and positive (beta = 27.59, 95% CI [8.76, 46.43], t(99) = 2.87, p = 0.004; Std. beta = 0.72, 95% CI [0.23, 1.21])
##   - The effect of Weight_g is statistically non-significant and positive (beta = 0.05, 95% CI [-7.61e-04, 0.09], t(99) = 1.93, p = 0.054; Std. beta = 0.24, 95% CI [-3.94e-03, 0.49])
## 
## 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%]
Species 4684.5 4684.5 33.9 8.24 (1) 0.006998 0.2 [0.04-0.38]
Weight_g 2113.4 2113.4 37.5 3.72 (1) 0.061390 0.09 [0-0.26]

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$Species)

    plotResiduals(simulationOutput, sub_trait$Weight_g)

–> Model diagnostics do not show apparent violations of model assumptions

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)