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 <- "EPOC"
# Set file path
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: 20.15 max: 687.15
## median: 306.425
## mean: 311.551
## estimated sd: 171.0271
## estimated skewness: 0.147712
## estimated kurtosis: 2.066642
# Fit various distributions
nd <- fitdist(sub_trait$value, "norm")
## $start.arg
## $start.arg$mean
## [1] 311.551
##
## $start.arg$sd
## [1] 170.2334
##
##
## $fix.arg
## NULL
gd <- fitdist(sub_trait$value, "gamma")
## $start.arg
## $start.arg$shape
## [1] 3.349413
##
## $start.arg$rate
## [1] 0.01075077
##
##
## $fix.arg
## NULL
wd <- fitdist(sub_trait$value, "weibull")
## $start.arg
## $start.arg$shape
## [1] 1.556164
##
## $start.arg$scale
## [1] 360.2923
##
##
## $fix.arg
## NULL
lnd <- fitdist(sub_trait$value, "lnorm")
## $start.arg
## $start.arg$meanlog
## [1] 5.519345
##
## $start.arg$sdlog
## [1] 0.7675485
##
##
## $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.06491519 0.09755319 0.08208289 0.1382808
## Cramer-von Mises statistic 0.09472074 0.25095048 0.12686354 0.5370490
## Anderson-Darling statistic 0.70245530 1.55418026 0.92079324 3.2428572
##
## Goodness-of-fit criteria
## Normal gamma Weibull LogNormal
## Akaike's Information Criterion 1420.120 1422.775 1413.881 1445.526
## Bayesian Information Criterion 1425.484 1428.139 1419.245 1450.890
–> 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 | 286.2 | 54 | 158 | 43.1 | 243.0 | 329.3 | 23.1 | 623.3 |
Sagmariasus verreauxi | 336.9 | 54 | 181 | 49.4 | 287.6 | 386.3 | 20.1 | 687.1 |
# 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 | 375.6 | 18 | 130.3 | 64.5 | 311.1 | 440.1 | 184.1 | 594.2 |
Jasus edwardsii | 17.5 | 326.2 | 18 | 155.0 | 76.8 | 249.4 | 402.9 | 36.2 | 623.3 |
Jasus edwardsii | 21.5 | 156.7 | 18 | 94.6 | 46.8 | 109.9 | 203.6 | 23.1 | 324.6 |
Sagmariasus verreauxi | 14 | 398.4 | 18 | 215.5 | 106.7 | 291.7 | 505.1 | 72.2 | 687.1 |
Sagmariasus verreauxi | 17.5 | 374.2 | 18 | 109.8 | 54.4 | 319.8 | 428.6 | 188.3 | 581.9 |
Sagmariasus verreauxi | 21.5 | 238.3 | 18 | 167.1 | 82.7 | 155.5 | 321.0 | 20.1 | 567.5 |
# 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 | 272.4 | 6 | 113.8 | 113.7 | 158.7 | 386.0 | 184.1 | 479.0 |
Jasus edwardsii | 14 | 17.5 | 290.7 | 6 | 83.9 | 83.8 | 206.8 | 374.5 | 151.8 | 378.2 |
Jasus edwardsii | 14 | 21.5 | 84.8 | 6 | 55.5 | 55.5 | 29.4 | 140.3 | 23.1 | 164.5 |
Jasus edwardsii | 17.5 | 14 | 436.7 | 6 | 147.0 | 146.9 | 289.8 | 583.6 | 237.3 | 594.2 |
Jasus edwardsii | 17.5 | 17.5 | 216.0 | 6 | 137.7 | 137.5 | 78.5 | 353.6 | 36.2 | 414.4 |
Jasus edwardsii | 17.5 | 21.5 | 204.6 | 6 | 86.3 | 86.2 | 118.4 | 290.8 | 90.3 | 324.6 |
Jasus edwardsii | 21.5 | 14 | 417.7 | 6 | 61.1 | 61.0 | 356.6 | 478.7 | 311.4 | 484.3 |
Jasus edwardsii | 21.5 | 17.5 | 471.8 | 6 | 119.0 | 118.8 | 353.0 | 590.7 | 289.6 | 623.3 |
Jasus edwardsii | 21.5 | 21.5 | 180.8 | 6 | 101.2 | 101.1 | 79.7 | 281.9 | 69.1 | 312.2 |
Sagmariasus verreauxi | 14 | 14 | 263.0 | 6 | 190.6 | 190.4 | 72.6 | 453.4 | 74.3 | 483.1 |
Sagmariasus verreauxi | 14 | 17.5 | 314.4 | 6 | 56.8 | 56.7 | 257.7 | 371.1 | 238.2 | 393.1 |
Sagmariasus verreauxi | 14 | 21.5 | 174.2 | 6 | 143.4 | 143.3 | 31.0 | 317.5 | 20.1 | 428.8 |
Sagmariasus verreauxi | 17.5 | 14 | 457.1 | 6 | 212.1 | 211.9 | 245.2 | 669.0 | 72.2 | 633.3 |
Sagmariasus verreauxi | 17.5 | 17.5 | 373.5 | 6 | 98.7 | 98.6 | 274.9 | 472.2 | 268.8 | 536.7 |
Sagmariasus verreauxi | 17.5 | 21.5 | 198.9 | 6 | 106.1 | 106.0 | 92.9 | 304.9 | 35.4 | 331.2 |
Sagmariasus verreauxi | 21.5 | 14 | 475.1 | 6 | 208.2 | 207.9 | 267.1 | 683.0 | 121.1 | 687.1 |
Sagmariasus verreauxi | 21.5 | 17.5 | 434.6 | 6 | 139.2 | 139.0 | 295.6 | 573.6 | 188.3 | 581.9 |
Sagmariasus verreauxi | 21.5 | 21.5 | 341.6 | 6 | 208.5 | 208.3 | 133.3 | 549.9 | 55.3 | 567.5 |
# 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
## 1372.2 1431.2 -664.1 1328.2 86
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.27183 -0.51859 0.03203 0.56554 2.06072
##
## Random effects:
## Groups Name Variance Std.Dev.
## AntTag (Intercept) 6146 78.40
## Residual 8727 93.42
## Number of obs: 108, groups: AntTag, 39
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) 404.04837 109.05957
## ExpTemp17.5 164.34167 53.93567
## ExpTemp21.5 126.33407 55.67596
## AccTemp17.5 5.76447 70.46407
## AccTemp21.5 -205.08643 71.13506
## SpeciesSagmariasus verreauxi -10.37149 69.45235
## SexMale -24.77882 34.55668
## Weight_g -0.09626 0.08857
## ExpTemp17.5:AccTemp17.5 -236.18767 76.32045
## ExpTemp21.5:AccTemp17.5 54.80760 77.51690
## ExpTemp17.5:AccTemp21.5 -44.07131 76.27798
## ExpTemp21.5:AccTemp21.5 -29.86205 77.52935
## ExpTemp17.5:SpeciesSagmariasus verreauxi 13.94273 77.50517
## ExpTemp21.5:SpeciesSagmariasus verreauxi 69.91533 78.76059
## AccTemp17.5:SpeciesSagmariasus verreauxi 40.86027 99.18984
## AccTemp21.5:SpeciesSagmariasus verreauxi 112.56804 99.93895
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi 115.67153 109.96846
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi -127.14657 110.51267
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi -111.45928 108.79288
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi -0.95355 109.72861
## df t value
## (Intercept) 49.14455 3.705
## ExpTemp17.5 69.35486 3.047
## ExpTemp21.5 76.64227 2.269
## AccTemp17.5 84.48444 0.082
## AccTemp21.5 83.66140 -2.883
## SpeciesSagmariasus verreauxi 88.90528 -0.149
## SexMale 39.56152 -0.717
## Weight_g 40.97006 -1.087
## ExpTemp17.5:AccTemp17.5 69.50981 -3.095
## ExpTemp21.5:AccTemp17.5 73.07288 0.707
## ExpTemp17.5:AccTemp21.5 69.35987 -0.578
## ExpTemp21.5:AccTemp21.5 73.08895 -0.385
## ExpTemp17.5:SpeciesSagmariasus verreauxi 73.17693 0.180
## ExpTemp21.5:SpeciesSagmariasus verreauxi 76.84449 0.888
## AccTemp17.5:SpeciesSagmariasus verreauxi 86.89592 0.412
## AccTemp21.5:SpeciesSagmariasus verreauxi 84.49161 1.126
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi 74.09935 1.052
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi 75.01382 -1.151
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi 71.43833 -1.025
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi 73.35665 -0.009
## Pr(>|t|)
## (Intercept) 0.000536 ***
## ExpTemp17.5 0.003268 **
## ExpTemp21.5 0.026076 *
## AccTemp17.5 0.934994
## AccTemp21.5 0.005005 **
## SpeciesSagmariasus verreauxi 0.881629
## SexMale 0.477556
## Weight_g 0.283487
## ExpTemp17.5:AccTemp17.5 0.002838 **
## ExpTemp21.5:AccTemp17.5 0.481790
## ExpTemp17.5:AccTemp21.5 0.565289
## ExpTemp21.5:AccTemp21.5 0.701229
## ExpTemp17.5:SpeciesSagmariasus verreauxi 0.857733
## ExpTemp21.5:SpeciesSagmariasus verreauxi 0.377475
## AccTemp17.5:SpeciesSagmariasus verreauxi 0.681398
## AccTemp21.5:SpeciesSagmariasus verreauxi 0.263199
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi 0.296282
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi 0.253586
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi 0.309052
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi 0.993090
## ---
## 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 to test main effects of final model representing the combined significance for all coefficients
anova(complex_model)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ExpTemp 383511 191756 2 71.826 21.9722 3.586e-08 ***
## AccTemp 241527 120764 2 38.043 13.8377 3.048e-05 ***
## Species 26472 26472 1 38.196 3.0333 0.08962 .
## Sex 4487 4487 1 39.562 0.5142 0.47756
## Weight_g 10308 10308 1 40.970 1.1811 0.28349
## ExpTemp:AccTemp 121937 30484 4 71.841 3.4930 0.01157 *
## ExpTemp:Species 3244 1622 2 72.144 0.1859 0.83078
## AccTemp:Species 8264 4132 2 38.364 0.4735 0.62642
## ExpTemp:AccTemp:Species 101204 25301 4 72.120 2.8991 0.02772 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 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.70) and the part related to the fixed effects alone (marginal R2) is of 0.48. The model's intercept, corresponding to ExpTemp = 14, AccTemp = 14, Species = Jasus edwardsii, Sex = Female and Weight_g = 0, is at 404.05 (95% CI [190.30, 617.80], t(86) = 3.70, p < .001). Within this model:
##
## - The effect of ExpTemp [17.5] is statistically significant and positive (beta = 164.34, 95% CI [58.63, 270.05], t(86) = 3.05, p = 0.002; Std. beta = 0.96, 95% CI [0.34, 1.58])
## - The effect of ExpTemp [21.5] is statistically significant and positive (beta = 126.33, 95% CI [17.21, 235.46], t(86) = 2.27, p = 0.023; Std. beta = 0.74, 95% CI [0.10, 1.38])
## - The effect of AccTemp [17.5] is statistically non-significant and positive (beta = 5.76, 95% CI [-132.34, 143.87], t(86) = 0.08, p = 0.935; Std. beta = 0.03, 95% CI [-0.77, 0.84])
## - The effect of AccTemp [21.5] is statistically significant and negative (beta = -205.09, 95% CI [-344.51, -65.66], t(86) = -2.88, p = 0.004; Std. beta = -1.20, 95% CI [-2.01, -0.38])
## - The effect of Species [Sagmariasus verreauxi] is statistically non-significant and negative (beta = -10.37, 95% CI [-146.50, 125.75], t(86) = -0.15, p = 0.881; Std. beta = -0.06, 95% CI [-0.86, 0.74])
## - The effect of Sex [Male] is statistically non-significant and negative (beta = -24.78, 95% CI [-92.51, 42.95], t(86) = -0.72, p = 0.473; Std. beta = -0.14, 95% CI [-0.54, 0.25])
## - The effect of Weight_g is statistically non-significant and negative (beta = -0.10, 95% CI [-0.27, 0.08], t(86) = -1.09, p = 0.277; Std. beta = -0.11, 95% CI [-0.30, 0.09])
## - The interaction effect of AccTemp [17.5] on ExpTemp [17.5] is statistically significant and negative (beta = -236.19, 95% CI [-385.77, -86.60], t(86) = -3.09, p = 0.002; Std. beta = -1.38, 95% CI [-2.26, -0.51])
## - The interaction effect of AccTemp [17.5] on ExpTemp [21.5] is statistically non-significant and positive (beta = 54.81, 95% CI [-97.12, 206.74], t(86) = 0.71, p = 0.480; Std. beta = 0.32, 95% CI [-0.57, 1.21])
## - The interaction effect of AccTemp [21.5] on ExpTemp [17.5] is statistically non-significant and negative (beta = -44.07, 95% CI [-193.57, 105.43], t(86) = -0.58, p = 0.563; Std. beta = -0.26, 95% CI [-1.13, 0.62])
## - The interaction effect of AccTemp [21.5] on ExpTemp [21.5] is statistically non-significant and negative (beta = -29.86, 95% CI [-181.82, 122.09], t(86) = -0.39, p = 0.700; Std. beta = -0.17, 95% CI [-1.06, 0.71])
## - The interaction effect of Species [Sagmariasus verreauxi] on ExpTemp [17.5] is statistically non-significant and positive (beta = 13.94, 95% CI [-137.96, 165.85], t(86) = 0.18, p = 0.857; Std. beta = 0.08, 95% CI [-0.81, 0.97])
## - The interaction effect of Species [Sagmariasus verreauxi] on ExpTemp [21.5] is statistically non-significant and positive (beta = 69.92, 95% CI [-84.45, 224.28], t(86) = 0.89, p = 0.375; Std. beta = 0.41, 95% CI [-0.49, 1.31])
## - The interaction effect of Species [Sagmariasus verreauxi] on AccTemp [17.5] is statistically non-significant and positive (beta = 40.86, 95% CI [-153.55, 235.27], t(86) = 0.41, p = 0.680; Std. beta = 0.24, 95% CI [-0.90, 1.38])
## - The interaction effect of Species [Sagmariasus verreauxi] on AccTemp [21.5] is statistically non-significant and positive (beta = 112.57, 95% CI [-83.31, 308.44], t(86) = 1.13, p = 0.260; Std. beta = 0.66, 95% CI [-0.49, 1.80])
## - The interaction effect of Species [Sagmariasus verreauxi] on (ExpTemp [17.5] * AccTemp [17.5]) is statistically non-significant and positive (beta = 115.67, 95% CI [-99.86, 331.21], t(86) = 1.05, p = 0.293; Std. beta = 0.68, 95% CI [-0.58, 1.94])
## - The interaction effect of Species [Sagmariasus verreauxi] on (ExpTemp [21.5] * AccTemp [17.5]) is statistically non-significant and negative (beta = -127.15, 95% CI [-343.75, 89.45], t(86) = -1.15, p = 0.250; Std. beta = -0.74, 95% CI [-2.01, 0.52])
## - The interaction effect of Species [Sagmariasus verreauxi] on (ExpTemp [17.5] * AccTemp [21.5]) is statistically non-significant and negative (beta = -111.46, 95% CI [-324.69, 101.77], t(86) = -1.02, p = 0.306; Std. beta = -0.65, 95% CI [-1.90, 0.60])
## - The interaction effect of Species [Sagmariasus verreauxi] on (ExpTemp [21.5] * AccTemp [21.5]) is statistically non-significant and negative (beta = -0.95, 95% CI [-216.02, 214.11], t(86) = -8.69e-03, p = 0.993; Std. beta = -5.58e-03, 95% CI [-1.26, 1.25])
##
## 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 -664.11 1372.2
## (1 | AntTag) 0 21 -672.10 1386.2 15.997 1 6.343e-05 ***
## ---
## 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 4487 4487.2 1 39.562 0.5142 0.47756
## Weight_g 2 16675 16675.0 1 40.813 1.9205 0.17333
## ExpTemp:AccTemp:Species 0 105765 26441.1 4 72.421 2.9797 0.02458
##
## Sex
## Weight_g
## ExpTemp:AccTemp:Species *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## value ~ ExpTemp + AccTemp + Species + (1 | AntTag) + ExpTemp:AccTemp + ExpTemp:Species + AccTemp:Species + ExpTemp:AccTemp:Species
# Extract the model that step found:
final_model <- get_model(step_result)
# Remodel using Maximum likelyhood estimation
final_model <- lmer(value ~ ExpTemp *
AccTemp *
Species +
(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 * Species + (1 | AntTag)
## Data: sub_trait
##
## REML criterion at convergence: 1157.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.09029 -0.48842 0.02245 0.54816 1.96764
##
## Random effects:
## Groups Name Variance Std.Dev.
## AntTag (Intercept) 7591 87.13
## Residual 10687 103.38
## Number of obs: 108, groups: AntTag, 39
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) 277.726 54.524
## ExpTemp17.5 164.342 59.685
## ExpTemp21.5 128.970 61.582
## AccTemp17.5 12.946 77.584
## AccTemp21.5 -192.913 77.584
## SpeciesSagmariasus verreauxi -6.515 76.912
## ExpTemp17.5:AccTemp17.5 -239.000 84.408
## ExpTemp21.5:AccTemp17.5 52.172 85.760
## ExpTemp17.5:AccTemp21.5 -44.577 84.408
## ExpTemp21.5:AccTemp21.5 -33.003 85.760
## ExpTemp17.5:SpeciesSagmariasus verreauxi 11.787 85.760
## ExpTemp21.5:SpeciesSagmariasus verreauxi 65.124 87.090
## AccTemp17.5:SpeciesSagmariasus verreauxi 25.837 109.107
## AccTemp21.5:SpeciesSagmariasus verreauxi 95.933 109.582
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi 126.674 121.282
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi -121.366 122.227
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi -106.857 120.330
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi 6.284 121.282
## df t value Pr(>|t|)
## (Intercept) 74.155 5.094 2.59e-06
## ExpTemp17.5 58.324 2.753 0.00785
## ExpTemp21.5 64.328 2.094 0.04018
## AccTemp17.5 71.542 0.167 0.86795
## AccTemp21.5 71.542 -2.486 0.01524
## SpeciesSagmariasus verreauxi 74.709 -0.085 0.93272
## ExpTemp17.5:AccTemp17.5 58.324 -2.831 0.00635
## ExpTemp21.5:AccTemp17.5 61.388 0.608 0.54520
## ExpTemp17.5:AccTemp21.5 58.324 -0.528 0.59943
## ExpTemp21.5:AccTemp21.5 61.388 -0.385 0.70169
## ExpTemp17.5:SpeciesSagmariasus verreauxi 61.388 0.137 0.89113
## ExpTemp21.5:SpeciesSagmariasus verreauxi 64.328 0.748 0.45732
## AccTemp17.5:SpeciesSagmariasus verreauxi 73.400 0.237 0.81347
## AccTemp21.5:SpeciesSagmariasus verreauxi 71.816 0.875 0.38425
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi 61.388 1.044 0.30037
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi 62.875 -0.993 0.32454
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi 59.870 -0.888 0.37808
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi 61.388 0.052 0.95885
##
## (Intercept) ***
## ExpTemp17.5 **
## ExpTemp21.5 *
## AccTemp17.5
## AccTemp21.5 *
## SpeciesSagmariasus verreauxi
## ExpTemp17.5:AccTemp17.5 **
## ExpTemp21.5:AccTemp17.5
## ExpTemp17.5:AccTemp21.5
## ExpTemp21.5:AccTemp21.5
## ExpTemp17.5:SpeciesSagmariasus verreauxi
## ExpTemp21.5:SpeciesSagmariasus verreauxi
## AccTemp17.5:SpeciesSagmariasus verreauxi
## AccTemp21.5:SpeciesSagmariasus verreauxi
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
# 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 386129 193064 2 60.352 18.0654 7.102e-07 ***
## AccTemp 233698 116849 2 32.806 10.9338 0.0002298 ***
## Species 21595 21595 1 32.794 2.0207 0.1646089
## ExpTemp:AccTemp 118592 29648 4 60.317 2.7742 0.0349077 *
## ExpTemp:Species 3263 1632 2 60.352 0.1527 0.8587350
## AccTemp:Species 5820 2910 2 32.806 0.2723 0.7633360
## ExpTemp:AccTemp:Species 105752 26438 4 60.317 2.4739 0.0537655 .
## ---
## 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) = 18.07, p < .001; Eta2 (partial) = 0.37, 90% CI [0.21, 0.50])
## - The main effect of AccTemp is statistically significant and large (F(2) = 10.93, p < .001; Eta2 (partial) = 0.40, 90% CI [0.17, 0.56])
## - The main effect of Species is statistically not significant and small (F(1) = 2.02, p = 0.165; Eta2 (partial) = 0.06, 90% CI [0.00, 0.22])
## - The interaction between ExpTemp and AccTemp is statistically significant and large (F(4) = 2.77, p = 0.035; Eta2 (partial) = 0.16, 90% CI [6.69e-03, 0.26])
## - The interaction between ExpTemp and Species is statistically not significant and very small (F(2) = 0.15, p = 0.859; Eta2 (partial) = 5.03e-03, 90% CI [0.00, 0.04])
## - The interaction between AccTemp and Species is statistically not significant and small (F(2) = 0.27, p = 0.763; Eta2 (partial) = 0.02, 90% CI [0.00, 0.10])
## - The interaction between ExpTemp, AccTemp and Species is statistically not significant and large (F(4) = 2.47, p = 0.054; Eta2 (partial) = 0.14, 90% CI [0.00, 0.25])
##
## Effect sizes were labelled following Field's (2013) recommendations.
# Table output from sjPlot
tab_model(final_model)
value | |||
---|---|---|---|
Predictors | Estimates | CI | p |
(Intercept) | 277.73 | 170.86 – 384.59 | <0.001 |
ExpTemp [17.5] | 164.34 | 47.36 – 281.32 | 0.006 |
ExpTemp [21.5] | 128.97 | 8.27 – 249.67 | 0.036 |
AccTemp [17.5] | 12.95 | -139.12 – 165.01 | 0.867 |
AccTemp [21.5] | -192.91 | -344.97 – -40.85 | 0.013 |
Species [Sagmariasus verreauxi] |
-6.51 | -157.26 – 144.23 | 0.932 |
ExpTemp [17.5] * AccTemp [17.5] |
-239.00 | -404.44 – -73.56 | 0.005 |
ExpTemp [21.5] * AccTemp [17.5] |
52.17 | -115.91 – 220.26 | 0.543 |
ExpTemp [17.5] * AccTemp [21.5] |
-44.58 | -210.01 – 120.86 | 0.597 |
ExpTemp [21.5] * AccTemp [21.5] |
-33.00 | -201.09 – 135.08 | 0.700 |
ExpTemp [17.5] * Species [Sagmariasus verreauxi] |
11.79 | -156.30 – 179.87 | 0.891 |
ExpTemp [21.5] * Species [Sagmariasus verreauxi] |
65.12 | -105.57 – 235.82 | 0.455 |
AccTemp [17.5] * Species [Sagmariasus verreauxi] |
25.84 | -188.01 – 239.68 | 0.813 |
AccTemp [21.5] * Species [Sagmariasus verreauxi] |
95.93 | -118.84 – 310.71 | 0.381 |
(ExpTemp [17.5] * AccTemp [17.5]) * Species [Sagmariasus verreauxi] |
126.67 | -111.03 – 364.38 | 0.296 |
(ExpTemp [21.5] * AccTemp [17.5]) * Species [Sagmariasus verreauxi] |
-121.37 | -360.93 – 118.19 | 0.321 |
(ExpTemp [17.5] * AccTemp [21.5]) * Species [Sagmariasus verreauxi] |
-106.86 | -342.70 – 128.99 | 0.375 |
(ExpTemp [21.5] * AccTemp [21.5]) * Species [Sagmariasus verreauxi] |
6.28 | -231.42 – 243.99 | 0.959 |
Random Effects | |||
σ2 | 10686.96 | ||
τ00 AntTag | 7591.36 | ||
ICC | 0.42 | ||
N AntTag | 39 | ||
Observations | 108 | ||
Marginal R2 / Conditional R2 | 0.422 / 0.662 |
# 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, AccTemp and Species (formula: value ~ ExpTemp * AccTemp * Species). The model included AntTag as random effect (formula: ~1 | AntTag). The model's total explanatory power is substantial (conditional R2 = 0.66) and the part related to the fixed effects alone (marginal R2) is of 0.42. The model's intercept, corresponding to ExpTemp = 14, AccTemp = 14 and Species = Jasus edwardsii, is at 277.73 (95% CI [170.86, 384.59], t(88) = 5.09, p < .001). Within this model:
##
## - The effect of ExpTemp [17.5] is statistically significant and positive (beta = 164.34, 95% CI [47.36, 281.32], t(88) = 2.75, p = 0.006; Std. beta = 0.96, 95% CI [0.28, 1.64])
## - The effect of ExpTemp [21.5] is statistically significant and positive (beta = 128.97, 95% CI [8.27, 249.67], t(88) = 2.09, p = 0.036; Std. beta = 0.75, 95% CI [0.05, 1.46])
## - The effect of AccTemp [17.5] is statistically non-significant and positive (beta = 12.95, 95% CI [-139.12, 165.01], t(88) = 0.17, p = 0.867; Std. beta = 0.08, 95% CI [-0.81, 0.96])
## - The effect of AccTemp [21.5] is statistically significant and negative (beta = -192.91, 95% CI [-344.97, -40.85], t(88) = -2.49, p = 0.013; Std. beta = -1.13, 95% CI [-2.02, -0.24])
## - The effect of Species [Sagmariasus verreauxi] is statistically non-significant and negative (beta = -6.51, 95% CI [-157.26, 144.23], t(88) = -0.08, p = 0.932; Std. beta = -0.04, 95% CI [-0.92, 0.84])
## - The interaction effect of AccTemp [17.5] on ExpTemp [17.5] is statistically significant and negative (beta = -239.00, 95% CI [-404.44, -73.56], t(88) = -2.83, p = 0.005; Std. beta = -1.40, 95% CI [-2.36, -0.43])
## - The interaction effect of AccTemp [17.5] on ExpTemp [21.5] is statistically non-significant and positive (beta = 52.17, 95% CI [-115.91, 220.26], t(88) = 0.61, p = 0.543; Std. beta = 0.31, 95% CI [-0.68, 1.29])
## - The interaction effect of AccTemp [21.5] on ExpTemp [17.5] is statistically non-significant and negative (beta = -44.58, 95% CI [-210.01, 120.86], t(88) = -0.53, p = 0.597; Std. beta = -0.26, 95% CI [-1.23, 0.71])
## - The interaction effect of AccTemp [21.5] on ExpTemp [21.5] is statistically non-significant and negative (beta = -33.00, 95% CI [-201.09, 135.08], t(88) = -0.38, p = 0.700; Std. beta = -0.19, 95% CI [-1.18, 0.79])
## - The interaction effect of Species [Sagmariasus verreauxi] on ExpTemp [17.5] is statistically non-significant and positive (beta = 11.79, 95% CI [-156.30, 179.87], t(88) = 0.14, p = 0.891; Std. beta = 0.07, 95% CI [-0.91, 1.05])
## - The interaction effect of Species [Sagmariasus verreauxi] on ExpTemp [21.5] is statistically non-significant and positive (beta = 65.12, 95% CI [-105.57, 235.82], t(88) = 0.75, p = 0.455; Std. beta = 0.38, 95% CI [-0.62, 1.38])
## - The interaction effect of Species [Sagmariasus verreauxi] on AccTemp [17.5] is statistically non-significant and positive (beta = 25.84, 95% CI [-188.01, 239.68], t(88) = 0.24, p = 0.813; Std. beta = 0.15, 95% CI [-1.10, 1.40])
## - The interaction effect of Species [Sagmariasus verreauxi] on AccTemp [21.5] is statistically non-significant and positive (beta = 95.93, 95% CI [-118.84, 310.71], t(88) = 0.88, p = 0.381; Std. beta = 0.56, 95% CI [-0.69, 1.82])
## - The interaction effect of Species [Sagmariasus verreauxi] on (ExpTemp [17.5] * AccTemp [17.5]) is statistically non-significant and positive (beta = 126.67, 95% CI [-111.03, 364.38], t(88) = 1.04, p = 0.296; Std. beta = 0.74, 95% CI [-0.65, 2.13])
## - The interaction effect of Species [Sagmariasus verreauxi] on (ExpTemp [21.5] * AccTemp [17.5]) is statistically non-significant and negative (beta = -121.37, 95% CI [-360.93, 118.19], t(88) = -0.99, p = 0.321; Std. beta = -0.71, 95% CI [-2.11, 0.69])
## - The interaction effect of Species [Sagmariasus verreauxi] on (ExpTemp [17.5] * AccTemp [21.5]) is statistically non-significant and negative (beta = -106.86, 95% CI [-342.70, 128.99], t(88) = -0.89, p = 0.375; Std. beta = -0.62, 95% CI [-2.00, 0.75])
## - The interaction effect of Species [Sagmariasus verreauxi] on (ExpTemp [21.5] * AccTemp [21.5]) is statistically non-significant and positive (beta = 6.28, 95% CI [-231.42, 243.99], t(88) = 0.05, p = 0.959; Std. beta = 0.04, 95% CI [-1.35, 1.43])
##
## 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 | 386128.9 | 193064.5 | 60.4 | 18.07 (2) | 7.102e-07 | 0.37 [0.21-0.5] |
AccTemp | 233698.1 | 116849.1 | 32.8 | 10.93 (2) | 2.298e-04 | 0.4 [0.17-0.56] |
Species | 21595.1 | 21595.1 | 32.8 | 2.02 (1) | 1.646e-01 | 0.06 [0-0.22] |
ExpTemp:AccTemp | 118591.6 | 29647.9 | 60.3 | 2.77 (4) | 3.491e-02 | 0.16 [0.01-0.26] |
ExpTemp:Species | 3263.4 | 1631.7 | 60.4 | 0.15 (2) | 8.587e-01 | 0.01 [0-0.04] |
AccTemp:Species | 5820.0 | 2910.0 | 32.8 | 0.27 (2) | 7.633e-01 | 0.02 [0-0.1] |
ExpTemp:AccTemp:Species | 105752.4 | 26438.1 | 60.3 | 2.47 (4) | 5.377e-02 | 0.14 [0-0.25] |
# 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)
plotResiduals(simulationOutput, sub_trait$Species)
# 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 | Species, type = "response", adjust = "bonferroni")
# Display interaction plots
emmip(final_model, AccTemp ~ ExpTemp | Species)
# 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:4,9,10)] %>% mutate(p.value = round(p.value, 3))
# Rename column headers
names(contrasts_expTemp)[1:6] <- c("Contrast", "ExpTemp", "Species", "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'))
Contrast | ExpTemp | Species | Difference | t-ratio | p | contrasts_sign | Trait |
---|---|---|---|---|---|---|---|
14 - 17.5 | 14 | Jasus edwardsii | -12.94558 | -0.1666996 | 1 | No | EPOC |
14 - 21.5 | 14 | Jasus edwardsii | 192.91275 | 2.4841280 | 0.046 | Yes | EPOC |
17.5 - 21.5 | 14 | Jasus edwardsii | 205.85833 | 2.6373083 | 0.031 | Yes | EPOC |
14 - 17.5 | 17.5 | Jasus edwardsii | 226.05442 | 2.9108916 | 0.014 | Yes | EPOC |
14 - 21.5 | 17.5 | Jasus edwardsii | 237.48942 | 3.0581396 | 0.009 | Yes | EPOC |
17.5 - 21.5 | 17.5 | Jasus edwardsii | 11.43500 | 0.1464970 | 1 | No | EPOC |
14 - 17.5 | 21.5 | Jasus edwardsii | -65.11770 | -0.8408618 | 1 | No | EPOC |
14 - 21.5 | 21.5 | Jasus edwardsii | 225.91563 | 2.9172380 | 0.014 | Yes | EPOC |
17.5 - 21.5 | 21.5 | Jasus edwardsii | 291.03333 | 3.7285090 | 0.001 | Yes | EPOC |
14 - 17.5 | 14 | Sagmariasus verreauxi | -38.78229 | -0.5048325 | 1 | No | EPOC |
14 - 21.5 | 14 | Sagmariasus verreauxi | 96.97970 | 1.2522943 | 0.644 | No | EPOC |
17.5 - 21.5 | 14 | Sagmariasus verreauxi | 135.76200 | 1.7530883 | 0.252 | No | EPOC |
14 - 17.5 | 17.5 | Sagmariasus verreauxi | 73.54351 | 0.9519209 | 1 | No | EPOC |
14 - 21.5 | 17.5 | Sagmariasus verreauxi | 248.41362 | 3.1988100 | 0.006 | Yes | EPOC |
17.5 - 21.5 | 17.5 | Sagmariasus verreauxi | 174.87011 | 2.2517939 | 0.082 | No | EPOC |
14 - 17.5 | 21.5 | Sagmariasus verreauxi | 30.41184 | 0.3936400 | 1 | No | EPOC |
14 - 21.5 | 21.5 | Sagmariasus verreauxi | 123.69862 | 1.5928610 | 0.347 | No | EPOC |
17.5 - 21.5 | 21.5 | Sagmariasus verreauxi | 93.28678 | 1.2012493 | 0.701 | No | EPOC |
# 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 | Species, type = "response", adjust = "bonferroni")
# Display interaction plots
emmip(final_model, AccTemp ~ ExpTemp | Species)
# 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:4,9,10)] %>% mutate(p.value = round(p.value, 3))
# Rename column headers
names(contrasts_expTemp)[1:6] <- c("Contrast", "AccTemp", "Species", "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'))
Contrast | AccTemp | Species | Difference | t-ratio | p | contrasts_sign | Trait |
---|---|---|---|---|---|---|---|
14 - 17.5 | 14 | Jasus edwardsii | -164.34167 | -2.7534752 | 0.024 | Yes | EPOC |
14 - 21.5 | 14 | Jasus edwardsii | -128.96955 | -2.0874881 | 0.123 | No | EPOC |
17.5 - 21.5 | 14 | Jasus edwardsii | 35.37212 | 0.5725296 | 1 | No | EPOC |
14 - 17.5 | 17.5 | Jasus edwardsii | 74.65833 | 1.2508688 | 0.648 | No | EPOC |
14 - 21.5 | 17.5 | Jasus edwardsii | -181.14167 | -3.0349521 | 0.011 | Yes | EPOC |
17.5 - 21.5 | 17.5 | Jasus edwardsii | -255.80000 | -4.2858209 | 0 | Yes | EPOC |
14 - 17.5 | 21.5 | Jasus edwardsii | -119.76500 | -2.0066119 | 0.149 | No | EPOC |
14 - 21.5 | 21.5 | Jasus edwardsii | -95.96667 | -1.6078809 | 0.34 | No | EPOC |
17.5 - 21.5 | 21.5 | Jasus edwardsii | 23.79833 | 0.3987310 | 1 | No | EPOC |
14 - 17.5 | 14 | Sagmariasus verreauxi | -176.12892 | -2.8508050 | 0.018 | Yes | EPOC |
14 - 21.5 | 14 | Sagmariasus verreauxi | -194.09392 | -3.1415847 | 0.008 | Yes | EPOC |
17.5 - 21.5 | 14 | Sagmariasus verreauxi | -17.96500 | -0.3009960 | 1 | No | EPOC |
14 - 17.5 | 17.5 | Sagmariasus verreauxi | -63.80312 | -1.0327109 | 0.917 | No | EPOC |
14 - 21.5 | 17.5 | Sagmariasus verreauxi | -124.89979 | -2.0216154 | 0.142 | No | EPOC |
17.5 - 21.5 | 17.5 | Sagmariasus verreauxi | -61.09667 | -1.0236488 | 0.931 | No | EPOC |
14 - 17.5 | 21.5 | Sagmariasus verreauxi | -24.69500 | -0.4137543 | 1 | No | EPOC |
14 - 21.5 | 21.5 | Sagmariasus verreauxi | -167.37500 | -2.8042974 | 0.021 | Yes | EPOC |
17.5 - 21.5 | 21.5 | Sagmariasus verreauxi | -142.68000 | -2.3905431 | 0.06 | No | EPOC |
# 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 ~ Species | ExpTemp , type = "response", adjust = "bonferroni")
## NOTE: Results may be misleading due to involvement in interactions
# Display interaction plots
emmip(final_model, ExpTemp ~ Species)
## NOTE: Results may be misleading due to involvement in interactions
# 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("Contrast", "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_Species_", 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'))
Contrast | ExpTemp | Difference | t-ratio | p | contrasts_sign | Trait |
---|---|---|---|---|---|---|
Jasus edwardsii - Sagmariasus verreauxi | 14 | -34.07520 | -0.761412 | 0.449 | No | EPOC |
Jasus edwardsii - Sagmariasus verreauxi | 17.5 | -52.46810 | -1.170224 | 0.246 | No | EPOC |
Jasus edwardsii - Sagmariasus verreauxi | 21.5 | -60.83881 | -1.358181 | 0.179 | No | EPOC |
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)