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)
# 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")
# 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)
–> 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
# 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")
# 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()
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 |
# 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()
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 |
# 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()
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 |
# 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
# 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 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.
# 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 |
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)