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)
# Select trait
trait <- "Recovery_Rate"
# Set working directory of executing R. script
file_path <- 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("C:/Users/oelle/Documents/Meine Dokumente/Research/Projects/Project_RockLobster/Temperature acclimation/Aerobic scope/Statistics/MO2_data_Tidy_RawData_long.csv")
# Shorten column names
names(mo2_data_long)[c(5,6)] <- c("AccTemp", "ExpTemp")
# Define fixed effects as factors
mo2_data_long$AntTag <- as.factor(mo2_data_long$AntTag)
mo2_data_long$Species <- as.factor(mo2_data_long$Species)
mo2_data_long$AccTemp <- as.factor(mo2_data_long$AccTemp)
mo2_data_long$ExpTemp <- as.factor(mo2_data_long$ExpTemp)
mo2_data_long$Sex <- as.factor(mo2_data_long$Sex)
# Subset data by trait
sub_trait <- subset(mo2_data_long, Trait == trait)
# Plot data versus a range of distrubutions
descdist(sub_trait$value, boot = 1000)
## summary statistics
## ------
## min: 6.65 max: 53.93
## median: 30.37
## mean: 29.92694
## estimated sd: 10.0487
## estimated skewness: -0.01235046
## estimated kurtosis: 2.536347
# Fit various distributions
nd <- fitdist(sub_trait$value, "norm")
## $start.arg
## $start.arg$mean
## [1] 29.92694
##
## $start.arg$sd
## [1] 10.00207
##
##
## $fix.arg
## NULL
gd <- fitdist(sub_trait$value, "gamma")
## $start.arg
## $start.arg$shape
## [1] 8.952508
##
## $start.arg$rate
## [1] 0.2991454
##
##
## $fix.arg
## NULL
wd <- fitdist(sub_trait$value, "weibull")
## $start.arg
## $start.arg$shape
## [1] 3.081664
##
## $start.arg$scale
## [1] 33.70825
##
##
## $fix.arg
## NULL
lnd <- fitdist(sub_trait$value, "lnorm")
## $start.arg
## $start.arg$meanlog
## [1] 3.332128
##
## $start.arg$sdlog
## [1] 0.3875931
##
##
## $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.04456903 0.09565532 0.04708087 0.1187847
## Cramer-von Mises statistic 0.02658116 0.16156035 0.02516870 0.3020099
## Anderson-Darling statistic 0.19316881 0.94949365 0.17188950 1.7792257
##
## Goodness-of-fit criteria
## Normal gamma Weibull LogNormal
## Akaike's Information Criterion 807.8939 814.9342 806.1422 825.5058
## Bayesian Information Criterion 813.2581 820.2985 811.5065 830.8701
–> 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 | 27.8 | 54 | 9.2 | 2.5 | 25.3 | 30.3 | 11.6 | 46.8 |
Sagmariasus verreauxi | 32.1 | 54 | 10.5 | 2.9 | 29.2 | 34.9 | 6.7 | 53.9 |
# Summary statistics
substrait_stats <- ddply(sub_trait, c("Species", "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 | AccTemp | mean | n | sd | ci | lower_ci | upper_ci | min | max |
---|---|---|---|---|---|---|---|---|---|
Jasus edwardsii | 14 | 32.4 | 18 | 8.8 | 4.3 | 28.1 | 36.8 | 15.2 | 44.7 |
Jasus edwardsii | 17.5 | 28.9 | 18 | 8.4 | 4.1 | 24.8 | 33.1 | 11.6 | 46.8 |
Jasus edwardsii | 21.5 | 22.0 | 18 | 7.5 | 3.7 | 18.3 | 25.7 | 12.8 | 33.1 |
Sagmariasus verreauxi | 14 | 32.1 | 18 | 13.1 | 6.5 | 25.6 | 38.6 | 6.7 | 51.0 |
Sagmariasus verreauxi | 17.5 | 34.7 | 18 | 9.3 | 4.6 | 30.1 | 39.3 | 20.2 | 53.9 |
Sagmariasus verreauxi | 21.5 | 29.4 | 18 | 8.5 | 4.2 | 25.2 | 33.6 | 12.3 | 42.6 |
# 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 | 29.4 | 6 | 6.9 | 6.9 | 22.4 | 36.3 | 18.3 | 38.0 |
Jasus edwardsii | 14 | 17.5 | 26.8 | 6 | 5.7 | 5.7 | 21.1 | 32.5 | 21.7 | 34.8 |
Jasus edwardsii | 14 | 21.5 | 16.2 | 6 | 4.1 | 4.1 | 12.1 | 20.3 | 12.8 | 22.0 |
Jasus edwardsii | 17.5 | 14 | 34.1 | 6 | 10.2 | 10.2 | 23.9 | 44.3 | 15.2 | 41.1 |
Jasus edwardsii | 17.5 | 17.5 | 23.5 | 6 | 7.5 | 7.5 | 16.0 | 31.0 | 11.6 | 34.1 |
Jasus edwardsii | 17.5 | 21.5 | 23.3 | 6 | 7.4 | 7.4 | 16.0 | 30.7 | 13.3 | 33.1 |
Jasus edwardsii | 21.5 | 14 | 33.9 | 6 | 9.6 | 9.5 | 24.3 | 43.4 | 16.7 | 44.7 |
Jasus edwardsii | 21.5 | 17.5 | 36.4 | 6 | 6.3 | 6.3 | 30.1 | 42.8 | 30.3 | 46.8 |
Jasus edwardsii | 21.5 | 21.5 | 26.5 | 6 | 7.1 | 7.1 | 19.4 | 33.7 | 13.4 | 32.7 |
Sagmariasus verreauxi | 14 | 14 | 29.6 | 6 | 8.3 | 8.2 | 21.3 | 37.8 | 17.6 | 38.1 |
Sagmariasus verreauxi | 14 | 17.5 | 29.0 | 6 | 5.5 | 5.5 | 23.5 | 34.5 | 20.2 | 37.4 |
Sagmariasus verreauxi | 14 | 21.5 | 22.3 | 6 | 6.2 | 6.2 | 16.1 | 28.6 | 12.3 | 30.6 |
Sagmariasus verreauxi | 17.5 | 14 | 36.0 | 6 | 13.5 | 13.5 | 22.5 | 49.4 | 10.7 | 48.5 |
Sagmariasus verreauxi | 17.5 | 17.5 | 35.1 | 6 | 7.6 | 7.6 | 27.5 | 42.7 | 23.5 | 45.1 |
Sagmariasus verreauxi | 17.5 | 21.5 | 28.9 | 6 | 4.9 | 4.9 | 23.9 | 33.8 | 18.9 | 32.1 |
Sagmariasus verreauxi | 21.5 | 14 | 30.7 | 6 | 17.6 | 17.5 | 13.2 | 48.3 | 6.7 | 51.0 |
Sagmariasus verreauxi | 21.5 | 17.5 | 39.9 | 6 | 11.5 | 11.5 | 28.4 | 51.4 | 25.0 | 53.9 |
Sagmariasus verreauxi | 21.5 | 21.5 | 37.1 | 6 | 7.1 | 7.1 | 29.9 | 44.2 | 23.6 | 42.6 |
# Relative decrease of recovery rate of warm acclimated lobsters compared to cold acclimated lobsters at 14°C experimental temperature
paste("Recovery rate was", round((1-16.2/29.4)*100, 2), "% lower in warm acclimated SRL")
## [1] "Recovery rate was 44.9 % lower in warm acclimated SRL"
# 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
summary(complex_model)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: value ~ ExpTemp * AccTemp * Species + Sex + Weight_g + (1 | AntTag)
## Data: sub_trait
##
## AIC BIC logLik deviance df.resid
## 765.8 824.8 -360.9 721.8 86
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.44065 -0.44013 -0.03318 0.60452 2.19281
##
## Random effects:
## Groups Name Variance Std.Dev.
## AntTag (Intercept) 33.57 5.794
## Residual 27.65 5.258
## Number of obs: 108, groups: AntTag, 39
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) 39.747868 7.288236
## ExpTemp17.5 4.758333 3.035913
## ExpTemp21.5 2.406650 3.162835
## AccTemp17.5 -3.686538 4.502355
## AccTemp21.5 -14.845524 4.549526
## SpeciesSagmariasus verreauxi -0.419718 4.410568
## SexMale 2.604575 2.356429
## Weight_g -0.010295 0.006005
## ExpTemp17.5:AccTemp17.5 -7.769204 4.297012
## ExpTemp21.5:AccTemp17.5 7.213350 4.384095
## ExpTemp17.5:AccTemp21.5 2.422384 4.293545
## ExpTemp21.5:AccTemp21.5 7.972400 4.384892
## ExpTemp17.5:SpeciesSagmariasus verreauxi 1.235820 4.384237
## ExpTemp21.5:SpeciesSagmariasus verreauxi -1.629164 4.475575
## AccTemp17.5:SpeciesSagmariasus verreauxi 2.838752 6.317436
## AccTemp21.5:SpeciesSagmariasus verreauxi 7.913093 6.385986
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi 8.012806 6.230063
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi 3.598234 6.265651
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi -2.077494 6.140789
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi 5.772490 6.208664
## df t value
## (Intercept) 49.641660 5.454
## ExpTemp17.5 69.667080 1.567
## ExpTemp21.5 75.765772 0.761
## AccTemp17.5 73.070064 -0.819
## AccTemp21.5 72.503374 -3.263
## SpeciesSagmariasus verreauxi 77.602399 -0.095
## SexMale 40.200074 1.105
## Weight_g 43.895066 -1.715
## ExpTemp17.5:AccTemp17.5 69.884010 -1.808
## ExpTemp21.5:AccTemp17.5 72.811401 1.645
## ExpTemp17.5:AccTemp21.5 69.674090 0.564
## ExpTemp21.5:AccTemp21.5 72.825970 1.818
## ExpTemp17.5:SpeciesSagmariasus verreauxi 72.978041 0.282
## ExpTemp21.5:SpeciesSagmariasus verreauxi 76.011573 -0.364
## AccTemp17.5:SpeciesSagmariasus verreauxi 75.362036 0.449
## AccTemp21.5:SpeciesSagmariasus verreauxi 72.992835 1.239
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi 74.306257 1.286
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi 74.484012 0.574
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi 71.579811 -0.338
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi 73.185982 0.930
## Pr(>|t|)
## (Intercept) 1.55e-06 ***
## ExpTemp17.5 0.12156
## ExpTemp21.5 0.44907
## AccTemp17.5 0.41556
## AccTemp21.5 0.00168 **
## SpeciesSagmariasus verreauxi 0.92443
## SexMale 0.27560
## Weight_g 0.09348 .
## ExpTemp17.5:AccTemp17.5 0.07490 .
## ExpTemp21.5:AccTemp17.5 0.10421
## ExpTemp17.5:AccTemp21.5 0.57444
## ExpTemp21.5:AccTemp21.5 0.07315 .
## ExpTemp17.5:SpeciesSagmariasus verreauxi 0.77883
## ExpTemp21.5:SpeciesSagmariasus verreauxi 0.71686
## AccTemp17.5:SpeciesSagmariasus verreauxi 0.65447
## AccTemp21.5:SpeciesSagmariasus verreauxi 0.21927
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi 0.20238
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi 0.56751
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi 0.73612
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi 0.35556
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 20 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
# Anova
anova(complex_model)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ExpTemp 1169.75 584.87 2 71.778 21.1526 5.999e-08 ***
## AccTemp 223.95 111.97 2 38.807 4.0497 0.0252738 *
## Species 132.96 132.96 1 38.969 4.8088 0.0343499 *
## Sex 33.78 33.78 1 40.200 1.2217 0.2755984
## Weight_g 81.29 81.29 1 43.895 2.9398 0.0934769 .
## ExpTemp:AccTemp 653.58 163.40 4 71.822 5.9094 0.0003615 ***
## ExpTemp:Species 45.09 22.54 2 72.187 0.8153 0.4465454
## AccTemp:Species 90.35 45.17 2 39.103 1.6337 0.2082451
## ExpTemp:AccTemp:Species 147.17 36.79 4 72.211 1.3307 0.2668184
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
report(anova(complex_model))
## Type 3 ANOVAs only give sensible and informative results when covariates are
## mean-centered and factors are coded with orthogonal contrasts (such as those
## produced by 'contr.sum', 'contr.poly', or 'contr.helmert', but *not* by the
## default 'contr.treatment').
## The ANOVA suggests that:
##
## - The main effect of ExpTemp is statistically significant and large (F(2) = 21.15, p < .001; Eta2 (partial) = 0.37, 90% CI [0.22, 0.49])
## - The main effect of AccTemp is statistically significant and large (F(2) = 4.05, p = 0.025; Eta2 (partial) = 0.17, 90% CI [0.01, 0.33])
## - The main effect of Species is statistically significant and medium (F(1) = 4.81, p = 0.034; Eta2 (partial) = 0.11, 90% CI [4.61e-03, 0.28])
## - The main effect of Sex is statistically not significant and small (F(1) = 1.22, p = 0.276; Eta2 (partial) = 0.03, 90% CI [0.00, 0.16])
## - The main effect of Weight_g is statistically not significant and medium (F(1) = 2.94, p = 0.093; Eta2 (partial) = 0.06, 90% CI [0.00, 0.21])
## - The interaction between ExpTemp and AccTemp is statistically significant and large (F(4) = 5.91, p < .001; Eta2 (partial) = 0.25, 90% CI [0.09, 0.36])
## - The interaction between ExpTemp and Species is statistically not significant and small (F(2) = 0.82, p = 0.447; Eta2 (partial) = 0.02, 90% CI [0.00, 0.09])
## - The interaction between AccTemp and Species is statistically not significant and medium (F(2) = 1.63, p = 0.208; Eta2 (partial) = 0.08, 90% CI [0.00, 0.21])
## - The interaction between ExpTemp, AccTemp and Species is statistically not significant and medium (F(4) = 1.33, p = 0.267; Eta2 (partial) = 0.07, 90% CI [0.00, 0.14])
##
## Effect sizes were labelled following Field's (2013) recommendations.
# Textual report of the model
report(complex_model)
## We fitted a linear mixed model (estimated using ML and nloptwrap optimizer) to predict value with ExpTemp, AccTemp, Species, Sex and Weight_g (formula: value ~ ExpTemp * AccTemp * Species + Sex + Weight_g). The model included AntTag as random effect (formula: ~1 | AntTag). The model's total explanatory power is substantial (conditional R2 = 0.72) and the part related to the fixed effects alone (marginal R2) is of 0.39. The model's intercept, corresponding to ExpTemp = 14, AccTemp = 14, Species = Jasus edwardsii, Sex = Female and Weight_g = 0, is at 39.75 (95% CI [25.46, 54.03], t(86) = 5.45, p < .001). Within this model:
##
## - The effect of ExpTemp [17.5] is statistically non-significant and positive (beta = 4.76, 95% CI [-1.19, 10.71], t(86) = 1.57, p = 0.117; Std. beta = 0.47, 95% CI [-0.12, 1.07])
## - The effect of ExpTemp [21.5] is statistically non-significant and positive (beta = 2.41, 95% CI [-3.79, 8.61], t(86) = 0.76, p = 0.447; Std. beta = 0.24, 95% CI [-0.38, 0.86])
## - The effect of AccTemp [17.5] is statistically non-significant and negative (beta = -3.69, 95% CI [-12.51, 5.14], t(86) = -0.82, p = 0.413; Std. beta = -0.37, 95% CI [-1.25, 0.51])
## - The effect of AccTemp [21.5] is statistically significant and negative (beta = -14.85, 95% CI [-23.76, -5.93], t(86) = -3.26, p = 0.001; Std. beta = -1.48, 95% CI [-2.36, -0.59])
## - The effect of Species [Sagmariasus verreauxi] is statistically non-significant and negative (beta = -0.42, 95% CI [-9.06, 8.22], t(86) = -0.10, p = 0.924; Std. beta = -0.04, 95% CI [-0.90, 0.82])
## - The effect of Sex [Male] is statistically non-significant and positive (beta = 2.60, 95% CI [-2.01, 7.22], t(86) = 1.11, p = 0.269; Std. beta = 0.26, 95% CI [-0.20, 0.72])
## - The effect of Weight_g is statistically non-significant and negative (beta = -0.01, 95% CI [-0.02, 1.47e-03], t(86) = -1.71, p = 0.086; Std. beta = -0.20, 95% CI [-0.42, 0.03])
## - The interaction effect of AccTemp [17.5] on ExpTemp [17.5] is statistically non-significant and negative (beta = -7.77, 95% CI [-16.19, 0.65], t(86) = -1.81, p = 0.071; Std. beta = -0.77, 95% CI [-1.61, 0.06])
## - The interaction effect of AccTemp [17.5] on ExpTemp [21.5] is statistically non-significant and positive (beta = 7.21, 95% CI [-1.38, 15.81], t(86) = 1.65, p = 0.100; Std. beta = 0.72, 95% CI [-0.14, 1.57])
## - The interaction effect of AccTemp [21.5] on ExpTemp [17.5] is statistically non-significant and positive (beta = 2.42, 95% CI [-5.99, 10.84], t(86) = 0.56, p = 0.573; Std. beta = 0.24, 95% CI [-0.60, 1.08])
## - The interaction effect of AccTemp [21.5] on ExpTemp [21.5] is statistically non-significant and positive (beta = 7.97, 95% CI [-0.62, 16.57], t(86) = 1.82, p = 0.069; Std. beta = 0.79, 95% CI [-0.06, 1.65])
## - The interaction effect of Species [Sagmariasus verreauxi] on ExpTemp [17.5] is statistically non-significant and positive (beta = 1.24, 95% CI [-7.36, 9.83], t(86) = 0.28, p = 0.778; Std. beta = 0.12, 95% CI [-0.73, 0.98])
## - The interaction effect of Species [Sagmariasus verreauxi] on ExpTemp [21.5] is statistically non-significant and negative (beta = -1.63, 95% CI [-10.40, 7.14], t(86) = -0.36, p = 0.716; Std. beta = -0.16, 95% CI [-1.04, 0.71])
## - The interaction effect of Species [Sagmariasus verreauxi] on AccTemp [17.5] is statistically non-significant and positive (beta = 2.84, 95% CI [-9.54, 15.22], t(86) = 0.45, p = 0.653; Std. beta = 0.28, 95% CI [-0.95, 1.51])
## - The interaction effect of Species [Sagmariasus verreauxi] on AccTemp [21.5] is statistically non-significant and positive (beta = 7.91, 95% CI [-4.60, 20.43], t(86) = 1.24, p = 0.215; Std. beta = 0.79, 95% CI [-0.46, 2.03])
## - The interaction effect of Species [Sagmariasus verreauxi] on (ExpTemp [17.5] * AccTemp [17.5]) is statistically non-significant and positive (beta = 8.01, 95% CI [-4.20, 20.22], t(86) = 1.29, p = 0.198; Std. beta = 0.80, 95% CI [-0.42, 2.01])
## - The interaction effect of Species [Sagmariasus verreauxi] on (ExpTemp [21.5] * AccTemp [17.5]) is statistically non-significant and positive (beta = 3.60, 95% CI [-8.68, 15.88], t(86) = 0.57, p = 0.566; Std. beta = 0.36, 95% CI [-0.86, 1.58])
## - The interaction effect of Species [Sagmariasus verreauxi] on (ExpTemp [17.5] * AccTemp [21.5]) is statistically non-significant and negative (beta = -2.08, 95% CI [-14.11, 9.96], t(86) = -0.34, p = 0.735; Std. beta = -0.21, 95% CI [-1.40, 0.99])
## - The interaction effect of Species [Sagmariasus verreauxi] on (ExpTemp [21.5] * AccTemp [21.5]) is statistically non-significant and positive (beta = 5.77, 95% CI [-6.40, 17.94], t(86) = 0.93, p = 0.353; Std. beta = 0.57, 95% CI [-0.64, 1.79])
##
## 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.
(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 -360.88 765.75
## (1 | AntTag) 0 21 -375.49 792.97 29.218 1 6.467e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## Sex 1 33.78 33.780 1 40.200 1.2217
## ExpTemp:AccTemp:Species 2 150.68 37.670 4 71.922 1.3585
## ExpTemp:Species 3 45.77 22.885 2 71.812 0.7662
## AccTemp:Species 4 96.25 48.123 2 38.642 1.5787
## Weight_g 5 65.42 65.422 1 43.868 2.1386
## Species 6 96.62 96.619 1 38.773 3.1003
## ExpTemp:AccTemp 0 657.85 164.463 4 71.494 5.2948
## Pr(>F)
## Sex 0.2755984
## ExpTemp:AccTemp:Species 0.2568937
## ExpTemp:Species 0.4685149
## AccTemp:Species 0.2192546
## Weight_g 0.1507646
## Species 0.0861660 .
## ExpTemp:AccTemp 0.0008594 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## value ~ ExpTemp + AccTemp + (1 | AntTag) + ExpTemp:AccTemp
# Extract the model that step found:
final_model <- get_model(step_result)
# Remodel using Maximum likelyhood estimation
final_model <- lmer(value ~ ExpTemp *
AccTemp +
(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 ~ ExpTemp * AccTemp + (1 | AntTag)
## Data: sub_trait
##
## REML criterion at convergence: 707.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.41722 -0.35111 0.02927 0.52148 1.86133
##
## Random effects:
## Groups Name Variance Std.Dev.
## AntTag (Intercept) 46.00 6.782
## Residual 33.96 5.828
## Number of obs: 108, groups: AntTag, 39
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 29.497 2.511 69.252 11.745 < 2e-16 ***
## ExpTemp17.5 5.297 2.428 66.302 2.182 0.03268 *
## ExpTemp21.5 1.483 2.478 68.891 0.599 0.55142
## AccTemp17.5 -1.752 3.572 67.578 -0.490 0.62541
## AccTemp21.5 -10.237 3.601 65.598 -2.842 0.00596 **
## ExpTemp17.5:AccTemp17.5 -3.733 3.434 66.300 -1.087 0.28094
## ExpTemp21.5:AccTemp17.5 8.953 3.469 67.616 2.580 0.01204 *
## ExpTemp17.5:AccTemp21.5 1.539 3.399 65.078 0.453 0.65220
## ExpTemp21.5:AccTemp21.5 11.044 3.435 66.434 3.215 0.00202 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ExT17.5 ExT21.5 AT17.5 AT21.5 ET17.5:AT1 ET21.5:AT1
## ExpTemp17.5 -0.481
## ExpTemp21.5 -0.493 0.510
## AccTemp17.5 -0.703 0.338 0.347
## AccTemp21.5 -0.697 0.335 0.344 0.490
## ET17.5:AT17 0.340 -0.707 -0.361 -0.478 -0.237
## ET21.5:AT17 0.352 -0.365 -0.714 -0.484 -0.246 0.515
## ET17.5:AT21 0.343 -0.714 -0.365 -0.241 -0.471 0.505 0.260
## ET21.5:AT21 0.356 -0.368 -0.721 -0.250 -0.477 0.260 0.515
## ET17.5:AT2
## ExpTemp17.5
## ExpTemp21.5
## AccTemp17.5
## AccTemp21.5
## ET17.5:AT17
## ET21.5:AT17
## ET17.5:AT21
## ET21.5:AT21 0.505
# Anova to test main effects of final model representing the combined significance for all coefficients
anova(final_model)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ExpTemp 1151.50 575.75 2 65.494 16.9520 1.167e-06 ***
## AccTemp 177.01 88.51 2 36.082 2.6059 0.087677 .
## ExpTemp:AccTemp 657.58 164.40 4 65.463 4.8404 0.001757 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
report(anova(final_model))
## Type 3 ANOVAs only give sensible and informative results when covariates are
## mean-centered and factors are coded with orthogonal contrasts (such as those
## produced by 'contr.sum', 'contr.poly', or 'contr.helmert', but *not* by the
## default 'contr.treatment').
## The ANOVA suggests that:
##
## - The main effect of ExpTemp is statistically significant and large (F(2) = 16.95, p < .001; Eta2 (partial) = 0.34, 90% CI [0.18, 0.47])
## - The main effect of AccTemp is statistically not significant and medium (F(2) = 2.61, p = 0.088; Eta2 (partial) = 0.13, 90% CI [0.00, 0.29])
## - The interaction between ExpTemp and AccTemp is statistically significant and large (F(4) = 4.84, p = 0.002; Eta2 (partial) = 0.23, 90% CI [0.06, 0.34])
##
## Effect sizes were labelled following Field's (2013) recommendations.
# Table output from sjPlot
tab_model(final_model)
value | |||
---|---|---|---|
Predictors | Estimates | CI | p |
(Intercept) | 29.50 | 24.57 – 34.42 | <0.001 |
ExpTemp [17.5] | 5.30 | 0.54 – 10.06 | 0.029 |
ExpTemp [21.5] | 1.48 | -3.37 – 6.34 | 0.549 |
AccTemp [17.5] | -1.75 | -8.75 – 5.25 | 0.624 |
AccTemp [21.5] | -10.24 | -17.30 – -3.18 | 0.004 |
ExpTemp [17.5] * AccTemp [17.5] |
-3.73 | -10.46 – 3.00 | 0.277 |
ExpTemp [21.5] * AccTemp [17.5] |
8.95 | 2.15 – 15.75 | 0.010 |
ExpTemp [17.5] * AccTemp [21.5] |
1.54 | -5.12 – 8.20 | 0.651 |
ExpTemp [21.5] * AccTemp [21.5] |
11.04 | 4.31 – 17.78 | 0.001 |
Random Effects | |||
σ2 | 33.96 | ||
τ00 AntTag | 46.00 | ||
ICC | 0.58 | ||
N AntTag | 39 | ||
Observations | 108 | ||
Marginal R2 / Conditional R2 | 0.244 / 0.679 |
# Textual report of the model
report(final_model)
## We fitted a linear mixed model (estimated using REML and nloptwrap optimizer) to predict value with ExpTemp and AccTemp (formula: value ~ ExpTemp * AccTemp). The model included AntTag as random effect (formula: ~1 | AntTag). The model's total explanatory power is substantial (conditional R2 = 0.68) and the part related to the fixed effects alone (marginal R2) is of 0.24. The model's intercept, corresponding to ExpTemp = 14 and AccTemp = 14, is at 29.50 (95% CI [24.57, 34.42], t(97) = 11.75, p < .001). Within this model:
##
## - The effect of ExpTemp [17.5] is statistically significant and positive (beta = 5.30, 95% CI [0.54, 10.06], t(97) = 2.18, p = 0.029; Std. beta = 0.53, 95% CI [0.05, 1.00])
## - The effect of ExpTemp [21.5] is statistically non-significant and positive (beta = 1.48, 95% CI [-3.37, 6.34], t(97) = 0.60, p = 0.549; Std. beta = 0.15, 95% CI [-0.34, 0.63])
## - The effect of AccTemp [17.5] is statistically non-significant and negative (beta = -1.75, 95% CI [-8.75, 5.25], t(97) = -0.49, p = 0.624; Std. beta = -0.17, 95% CI [-0.87, 0.52])
## - The effect of AccTemp [21.5] is statistically significant and negative (beta = -10.24, 95% CI [-17.30, -3.18], t(97) = -2.84, p = 0.004; Std. beta = -1.02, 95% CI [-1.72, -0.32])
## - The interaction effect of AccTemp [17.5] on ExpTemp [17.5] is statistically non-significant and negative (beta = -3.73, 95% CI [-10.46, 3.00], t(97) = -1.09, p = 0.277; Std. beta = -0.37, 95% CI [-1.04, 0.30])
## - The interaction effect of AccTemp [17.5] on ExpTemp [21.5] is statistically significant and positive (beta = 8.95, 95% CI [2.15, 15.75], t(97) = 2.58, p = 0.010; Std. beta = 0.89, 95% CI [0.21, 1.57])
## - The interaction effect of AccTemp [21.5] on ExpTemp [17.5] is statistically non-significant and positive (beta = 1.54, 95% CI [-5.12, 8.20], t(97) = 0.45, p = 0.651; Std. beta = 0.15, 95% CI [-0.51, 0.82])
## - The interaction effect of AccTemp [21.5] on ExpTemp [21.5] is statistically significant and positive (beta = 11.04, 95% CI [4.31, 17.78], t(97) = 3.21, p = 0.001; Std. beta = 1.10, 95% CI [0.43, 1.77])
##
## 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.
library(effectsize)
# Calculate anova test stats for model
aov_lme <-anova(final_model)
# Caculate partial effects size
eta = eta_squared(aov_lme)
# Create output stats table
lme_table <- round(aov_lme[, c(1,2,4)], 1)
# Attach F statistics
lme_table$F_value <- paste(round(aov_lme[,5],2), " (", aov_lme[,3], ")", sep = "")
# Attach p value
lme_table$p_value <- format(aov_lme[, 6], digits = 4)
# Create string with eta and 90% eta range and join with table
lme_table$Partial_Eta <- paste(round(eta$Eta2_partial, 2), " [", round(eta$CI_low, 2), "-", round(eta$CI_high, 2), "]", sep = "")
# Save table as csv
write.csv(lme_table, "Model_summary.csv")
# Print table
kable(lme_table,
format = "html",
row.names = TRUE,
caption = "Model summary for main effects",
col.names = c("Sum of sqares", "Mean squares", "Den df", "F(df)", "p value", "Partical squared eta (CI 90%)")) %>%
kable_classic(full_width = F, html_font = "Cambria")
Sum of sqares | Mean squares | Den df | F(df) | p value | Partical squared eta (CI 90%) | |
---|---|---|---|---|---|---|
ExpTemp | 1151.5 | 575.7 | 65.5 | 16.95 (2) | 1.167e-06 | 0.34 [0.18-0.47] |
AccTemp | 177.0 | 88.5 | 36.1 | 2.61 (2) | 8.768e-02 | 0.13 [0-0.29] |
ExpTemp:AccTemp | 657.6 | 164.4 | 65.5 | 4.84 (4) | 1.757e-03 | 0.23 [0.06-0.34] |
# 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$AccTemp)
plotResiduals(simulationOutput, sub_trait$ExpTemp)
# Perform multiple comparisions using the emmeans package
# Instructions here https://cran.r-project.org/web/packages/emmeans/vignettes/interactions.html
contrasts <- emmeans(final_model, pairwise ~ AccTemp | ExpTemp, type = "response", adjust = "bonferroni")
# Display interaction plots
emmip(final_model, AccTemp ~ ExpTemp)
# Convert results to dataframe (https://aosmith.rbind.io/2019/03/25/getting-started-with-emmeans/)
contrasts_df <- contrasts$contrasts %>%
summary(infer = TRUE) %>%
as.data.frame()
# Plot comparisons
### blue bars are confidence intervals for the EMMs, and the red arrows are for the comparisons among them. If an arrow ### from one mean overlaps an arrow from another group, the difference is not “significant,” based on the adjust setting ### (which defaults to "tukey") and the value of alpha (which defaults to 0.05)
### https://cran.r-project.org/web/packages/emmeans/vignettes/comparisons.html
plot(contrasts, comparisons = TRUE)
# Save contrasts in table
contrasts_expTemp <- contrasts_df[, c(1:3,8,9)] %>% mutate(p.value = round(p.value, 3))
# Rename column headers
names(contrasts_expTemp)[1:5] <- c("AccTemp", "ExpTemp", "Difference", "t-ratio", "p")
# Add significance descriptor
contrasts_expTemp$contrasts_sign <- ifelse(contrasts_expTemp$p < 0.05, "Yes", "No")
# Add trait columns
contrasts_expTemp$Trait <- trait
# Save contrasts results
write.csv(contrasts_expTemp,paste("Contrasts_AcclimationTemperature_", trait, ".csv", sep = ""))
# Display as HTML table colouring significant p values red
contrasts_expTemp %>% mutate(p = cell_spec(p, "html", color = ifelse(p < 0.05, "red", "grey"))) %>%
kable(escape = FALSE) %>%
kable_styling(bootstrap_options = c('hover'))
AccTemp | ExpTemp | Difference | t-ratio | p | contrasts_sign | Trait |
---|---|---|---|---|---|---|
14 - 17.5 | 14 | 1.751647 | 0.4901131 | 1 | No | Recovery_Rate |
14 - 21.5 | 14 | 10.236745 | 2.8410039 | 0.018 | Yes | Recovery_Rate |
17.5 - 21.5 | 14 | 8.485099 | 2.3428044 | 0.067 | No | Recovery_Rate |
14 - 17.5 | 17.5 | 5.484422 | 1.5294321 | 0.393 | No | Recovery_Rate |
14 - 21.5 | 17.5 | 8.697475 | 2.4104881 | 0.056 | No | Recovery_Rate |
17.5 - 21.5 | 17.5 | 3.213053 | 0.8854707 | 1 | No | Recovery_Rate |
14 - 17.5 | 21.5 | -7.201185 | -2.0109854 | 0.145 | No | Recovery_Rate |
14 - 21.5 | 21.5 | -0.807299 | -0.2240497 | 1 | No | Recovery_Rate |
17.5 - 21.5 | 21.5 | 6.393886 | 1.7620621 | 0.249 | No | Recovery_Rate |
# Perform multiple comparisions using the emmeans package
# Instructions here https://cran.r-project.org/web/packages/emmeans/vignettes/interactions.html
contrasts <- emmeans(final_model, pairwise ~ ExpTemp | AccTemp, type = "response", adjust = "bonferroni")
# Display interaction plots
emmip(final_model, ExpTemp ~ AccTemp)
# Convert results to dataframe (https://aosmith.rbind.io/2019/03/25/getting-started-with-emmeans/)
contrasts_df <- contrasts$contrasts %>%
summary(infer = TRUE) %>%
as.data.frame()
# Plot comparisons
### blue bars are confidence intervals for the EMMs, and the red arrows are for the comparisons among them. If an arrow ### from one mean overlaps an arrow from another group, the difference is not “significant,” based on the adjust setting ### (which defaults to "tukey") and the value of alpha (which defaults to 0.05)
### https://cran.r-project.org/web/packages/emmeans/vignettes/comparisons.html
plot(contrasts, comparisons = TRUE)
# Save contrasts in table
contrasts_expTemp <- contrasts_df[, c(1:3,8,9)] %>% mutate(p.value = round(p.value, 3))
# Rename column headers
names(contrasts_expTemp)[1:5] <- c("ExpTemp", "AccTemp", "Difference", "t-ratio", "p")
# Add significance descriptor
contrasts_expTemp$contrasts_sign <- ifelse(contrasts_expTemp$p < 0.05, "Yes", "No")
# Add trait columns
contrasts_expTemp$Trait <- trait
# Save contrasts results
write.csv(contrasts_expTemp,paste("Contrasts_ExperimentalTemperature_", trait, ".csv", sep = ""))
# Display as HTML table colouring significant p values red
contrasts_expTemp %>% mutate(p = cell_spec(p, "html", color = ifelse(p < 0.05, "red", "grey"))) %>%
kable(escape = FALSE) %>%
kable_styling(bootstrap_options = c('hover'))
ExpTemp | AccTemp | Difference | t-ratio | p | contrasts_sign | Trait |
---|---|---|---|---|---|---|
14 - 17.5 | 14 | -5.297396 | -2.1788600 | 0.099 | No | Recovery_Rate |
14 - 21.5 | 14 | -1.483456 | -0.5970262 | 1 | No | Recovery_Rate |
17.5 - 21.5 | 14 | 3.813940 | 1.5687032 | 0.365 | No | Recovery_Rate |
14 - 17.5 | 17.5 | -1.564621 | -0.6435458 | 1 | No | Recovery_Rate |
14 - 21.5 | 17.5 | -10.436288 | -4.2925598 | 0 | Yes | Recovery_Rate |
17.5 - 21.5 | 17.5 | -8.871667 | -3.7288555 | 0.001 | Yes | Recovery_Rate |
14 - 17.5 | 21.5 | -6.836667 | -2.8735234 | 0.017 | Yes | Recovery_Rate |
14 - 21.5 | 21.5 | -12.527500 | -5.2654409 | 0 | Yes | Recovery_Rate |
17.5 - 21.5 | 21.5 | -5.690833 | -2.3919175 | 0.059 | No | Recovery_Rate |
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)