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 <- "escape_speed_CLSec"
# File path
file_path <- "C:/Users/oelle/Documents/Meine Dokumente/Research/Projects/Project_RockLobster/Temperature acclimation/Exhaustion/Analysis/LinearMixedEffectModels/EscapesTotal/"
# Load tidy raw data in long format
Escape_MO2_data_long <- read_xlsx("C:/Users/oelle/Documents/Meine Dokumente/Research/Projects/Project_RockLobster/Temperature acclimation/Exhaustion/Analysis/Escape_MO2_MergedData_long.xlsx", sheet = 1)
# Shorten column names
names(Escape_MO2_data_long)[c(4,5)] <- c("AccTemp", "ExpTemp")
# Define fixed effects as factors
Escape_MO2_data_long$AntTag <- as.factor(Escape_MO2_data_long$AntTag)
Escape_MO2_data_long$Species <- as.factor(Escape_MO2_data_long$Species)
Escape_MO2_data_long$AccTemp <- as.factor(Escape_MO2_data_long$AccTemp)
Escape_MO2_data_long$ExpTemp <- as.factor(Escape_MO2_data_long$ExpTemp)
Escape_MO2_data_long$Sex <- as.factor(Escape_MO2_data_long$Sex)
# Subset data by trait
sub_trait <- subset(Escape_MO2_data_long, Trait == trait)
# Remove any columns with NAs
sub_trait <- sub_trait %>% drop_na("value")
# Plot data versus a range of distrubutions
descdist(sub_trait$value, boot = 1000)
## summary statistics
## ------
## min: 2.37285 max: 7.847453
## median: 5.1999
## mean: 5.247999
## estimated sd: 0.9990123
## estimated skewness: 0.2011247
## estimated kurtosis: 3.410981
# Fit various distributions
nd <- fitdist(sub_trait$value, "norm")
## $start.arg
## $start.arg$mean
## [1] 5.247999
##
## $start.arg$sd
## [1] 0.9941032
##
##
## $fix.arg
## NULL
gd <- fitdist(sub_trait$value, "gamma")
## $start.arg
## $start.arg$shape
## [1] 27.86921
##
## $start.arg$rate
## [1] 5.310444
##
##
## $fix.arg
## NULL
wd <- fitdist(sub_trait$value, "weibull")
## $start.arg
## $start.arg$shape
## [1] 6.07508
##
## $start.arg$scale
## [1] 5.659492
##
##
## $fix.arg
## NULL
lnd <- fitdist(sub_trait$value, "lnorm")
## $start.arg
## $start.arg$meanlog
## [1] 1.639179
##
## $start.arg$sdlog
## [1] 0.1965576
##
##
## $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.04624564 0.04440644 0.06594791 0.05185765
## Cramer-von Mises statistic 0.03048918 0.02165731 0.11649361 0.04036483
## Anderson-Darling statistic 0.23512578 0.18069049 0.84665220 0.30091269
##
## Goodness-of-fit criteria
## Normal gamma Weibull LogNormal
## Akaike's Information Criterion 292.2569 293.1328 298.6669 295.9888
## Bayesian Information Criterion 297.5069 298.3827 303.9168 301.2388
–> Anderson-Darling statistics gives good measure for fit for both middle and tail of distribution, the lower the better
# Experimental temperature versus acclimation
ggplot(sub_trait, aes(x = ExpTemp, y = value, fill = AccTemp, alpha = .5)) +
geom_boxplot(outlier.size = 0) +
ylab(trait) +
facet_wrap(~ Species, scales = "free")+
geom_point(pch = 21, position = position_jitterdodge(seed = 100))+
geom_text(aes(label=AntTag, color = AccTemp), alpha = .8, position = position_jitterdodge(seed = 100))
ggsave("Escapes_speed_AccTemp.pdf", width = 14, height = 8)
# Line Plot comparing both species
ggline(sub_trait, x = "ExpTemp", y = "value",
add = c("mean_ci", "jitter"), # Add mean_se and jitter points
# add.params = list(size = 0.7), # Add point size
# label = "AntTag", # Add point labels
# #label.select = list(top.up = 2), # show only labels for the top 2 points
# font.label = list(color = "Species"), # Color labels by .y., here gene names
# repel = TRUE, # Use repel to avoid labels overplotting
color = "Species",
palette = "jco")
# Calculate means across groups
value_means <- ddply(sub_trait, c("Species", "ExpTemp", "AccTemp"), summarise,
mean = mean(value, na.rm = TRUE),
n = length(value),
sd = sd(value, na.rm = TRUE),
ci = qt(0.975, n) * sd / sqrt(n),
lower_ci = mean-ci,
upper_ci = mean+ci)
print(value_means)
## Species ExpTemp AccTemp mean n sd ci
## 1 Jasus edwardsii 14 14 5.032607 6 0.5782466 0.5776380
## 2 Jasus edwardsii 14 17.5 4.839765 6 0.7271401 0.7263748
## 3 Jasus edwardsii 14 21.5 4.495885 6 0.3680634 0.3676761
## 4 Jasus edwardsii 17.5 14 5.096463 6 0.9752342 0.9742078
## 5 Jasus edwardsii 17.5 17.5 5.104246 6 1.4954577 1.4938839
## 6 Jasus edwardsii 17.5 21.5 4.615829 5 0.7168912 0.8241375
## 7 Jasus edwardsii 21.5 14 5.395649 6 1.7527376 1.7508929
## 8 Jasus edwardsii 21.5 17.5 6.196764 5 1.0016904 1.1515424
## 9 Jasus edwardsii 21.5 21.5 5.384381 6 0.6674706 0.6667681
## 10 Sagmariasus verreauxi 14 14 5.491680 6 0.8945832 0.8936417
## 11 Sagmariasus verreauxi 14 17.5 5.863798 6 1.2729140 1.2715744
## 12 Sagmariasus verreauxi 14 21.5 4.596102 5 1.1564672 1.3294737
## 13 Sagmariasus verreauxi 17.5 14 4.934275 5 0.5698294 0.6550754
## 14 Sagmariasus verreauxi 17.5 17.5 5.472298 5 1.1384884 1.3088053
## 15 Sagmariasus verreauxi 17.5 21.5 5.026087 6 0.8987795 0.8978336
## 16 Sagmariasus verreauxi 21.5 14 5.619175 5 0.9409606 1.0817275
## 17 Sagmariasus verreauxi 21.5 17.5 5.694452 6 0.8109461 0.8100927
## 18 Sagmariasus verreauxi 21.5 21.5 5.595601 6 0.5152544 0.5147121
## lower_ci upper_ci
## 1 4.454969 5.610245
## 2 4.113390 5.566139
## 3 4.128209 4.863561
## 4 4.122255 6.070670
## 5 3.610362 6.598130
## 6 3.791692 5.439967
## 7 3.644756 7.146542
## 8 5.045221 7.348306
## 9 4.717613 6.051149
## 10 4.598038 6.385321
## 11 4.592224 7.135373
## 12 3.266628 5.925576
## 13 4.279200 5.589351
## 14 4.163493 6.781103
## 15 4.128253 5.923920
## 16 4.537448 6.700903
## 17 4.884359 6.504545
## 18 5.080889 6.110314
# Experimental temperature versus acclimation - Means +/- CI
pd <- position_dodge(0.3) # move them .05 to the left and right
ggplot(value_means, aes(x=ExpTemp, y=mean, group = AccTemp, colour=AccTemp)) +
facet_wrap(~ Species) +
geom_line(position=pd) +
geom_point(position=pd) +
geom_errorbar(aes(ymin=lower_ci, ymax=upper_ci), width=.1, position=pd)
# Acclimation temperature only
ggplot(sub_trait, aes(x = Species, y = value, fill = AccTemp, alpha = .5)) +
geom_boxplot(outlier.size = 0) +
ylab(trait) +
geom_point(pch = 21, position = position_jitterdodge(seed = 100))+
geom_text(aes(label=AntTag, color = Species), alpha = .8, position = position_jitterdodge(seed = 100))
ggsave("Escapes_speed_AccTemp.pdf", width = 14, height = 8)
# Experimental temperature only
ggplot(sub_trait, aes(x = ExpTemp, y = value, fill = Species, alpha = .5)) +
geom_boxplot(outlier.size = 0) +
ylab(trait) +
geom_point(pch = 21, position = position_jitterdodge(seed = 100))+
geom_text(aes(label=AntTag, color = Species), alpha = .8, position = position_jitterdodge(seed = 100))
ggsave("Escapes_speed_ExpTemp.pdf", width = 14, height = 8)
# Species
ggplot(sub_trait, aes(x = Species, y = value, alpha = .5)) +
geom_boxplot(outlier.size = 0) +
ylab(trait) +
geom_point(pch = 21, position = position_jitterdodge(seed = 100))+
geom_text(aes(label=AntTag, color = Species), alpha = .8, position = position_jitterdodge(seed = 100))
ggsave("Escapes_speed_Species.pdf", width = 14, height = 8)
# 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 | 5.1 | 52 | 1.0 | 0.3 | 4.8 | 5.4 | 2.4 | 7.8 |
Sagmariasus verreauxi | 5.4 | 50 | 0.9 | 0.3 | 5.1 | 5.7 | 3.2 | 7.8 |
# 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 | 5.2 | 18 | 1.1 | 0.6 | 4.6 | 5.7 | 3.5 | 7.8 |
Jasus edwardsii | 17.5 | 5.3 | 17 | 1.2 | 0.6 | 4.7 | 6.0 | 2.4 | 7.7 |
Jasus edwardsii | 21.5 | 4.8 | 17 | 0.7 | 0.4 | 4.5 | 5.2 | 3.9 | 6.1 |
Sagmariasus verreauxi | 14 | 5.4 | 16 | 0.8 | 0.4 | 4.9 | 5.8 | 4.3 | 7.0 |
Sagmariasus verreauxi | 17.5 | 5.7 | 17 | 1.0 | 0.5 | 5.2 | 6.2 | 4.0 | 7.8 |
Sagmariasus verreauxi | 21.5 | 5.1 | 17 | 0.9 | 0.5 | 4.6 | 5.6 | 3.2 | 6.2 |
# Pooled data
ggscatter(sub_trait, x = "value", y = "EPOC", fill = "Species",color = "Species",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Escape speed", ylab = "EPOC")
## `geom_smooth()` using formula 'y ~ x'
ggsave("Escapes_speed_Corelation_EPOC.pdf", width = 8, height = 6)
## `geom_smooth()` using formula 'y ~ x'
### By species
# Subset
ERL <- subset(sub_trait, Species == "Sagmariasus verreauxi")
SRL <- subset(sub_trait, Species == "Jasus edwardsii")
# Normality tests
shapiro.test(ERL$value)
##
## Shapiro-Wilk normality test
##
## data: ERL$value
## W = 0.9942, p-value = 0.9972
shapiro.test(SRL$value)
##
## Shapiro-Wilk normality test
##
## data: SRL$value
## W = 0.97421, p-value = 0.3162
# Correlation plots
ggscatter(ERL, x = "value", y = "EPOC",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Escape speed", ylab = "EPOC")+
ggtitle("Sagmariasus verreauxi")
## `geom_smooth()` using formula 'y ~ x'
ggscatter(SRL, x = "value", y = "EPOC",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Escape speed", ylab = "EPOC")+
ggtitle("Jasus edwardsii")
## `geom_smooth()` using formula 'y ~ x'
# Correlation tests
cor.test(ERL$value, ERL$EPOC, method = "pearson")
##
## Pearson's product-moment correlation
##
## data: ERL$value and ERL$EPOC
## t = 2.0465, df = 48, p-value = 0.0462
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.005365486 0.520587842
## sample estimates:
## cor
## 0.2832901
cor.test(SRL$value, SRL$EPOC, method = "pearson")
##
## Pearson's product-moment correlation
##
## data: SRL$value and SRL$EPOC
## t = 2.5634, df = 50, p-value = 0.01342
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.07488617 0.56149668
## sample estimates:
## cor
## 0.3408213
# Pooled data
ggscatter(sub_trait, x = "value", y = "MMR", fill = "Species",color = "Species",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Escape speed", ylab = "MMR")
## `geom_smooth()` using formula 'y ~ x'
ggsave("Escapes_speed_Corelation_MMR.pdf", width = 8, height = 6)
## `geom_smooth()` using formula 'y ~ x'
### By species
# Subset
ERL <- subset(sub_trait, Species == "Sagmariasus verreauxi")
SRL <- subset(sub_trait, Species == "Jasus edwardsii")
# Normality tests
shapiro.test(ERL$value)
##
## Shapiro-Wilk normality test
##
## data: ERL$value
## W = 0.9942, p-value = 0.9972
shapiro.test(SRL$value)
##
## Shapiro-Wilk normality test
##
## data: SRL$value
## W = 0.97421, p-value = 0.3162
# Correlation plots
ggscatter(ERL, x = "value", y = "MMR",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Escape speed", ylab = "MMR")+
ggtitle("Sagmariasus verreauxi")
## `geom_smooth()` using formula 'y ~ x'
ggscatter(SRL, x = "value", y = "MMR",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Escape speed", ylab = "MMR")+
ggtitle("Jasus edwardsii")
## `geom_smooth()` using formula 'y ~ x'
# Correlation tests
cor.test(ERL$value, ERL$MMR, method = "pearson")
##
## Pearson's product-moment correlation
##
## data: ERL$value and ERL$MMR
## t = 1.4862, df = 48, p-value = 0.1438
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.07286133 0.46116470
## sample estimates:
## cor
## 0.2097401
cor.test(SRL$value, SRL$MMR, method = "pearson")
##
## Pearson's product-moment correlation
##
## data: SRL$value and SRL$MMR
## t = 3.193, df = 50, p-value = 0.002436
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1561909 0.6153395
## sample estimates:
## cor
## 0.4115491
names(sub_trait)
## [1] "Species" "AntTag" "AnimalID"
## [4] "AccTemp" "ExpTemp" "Sex"
## [7] "Weight_g" "Time_since_moult_d" "Acclimation_time"
## [10] "MMR" "SMR" "AerobicScope"
## [13] "recovery_time" "EPOC" "Recovery_Rate"
## [16] "Trait" "value"
# Start with the most complex model
complex_model <- lmer(value ~ ExpTemp *
AccTemp *
Species +
Sex +
Weight_g +
(1|AntTag), data = sub_trait, REML = FALSE)
# Model summary
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
## 282.9 340.7 -119.5 238.9 80
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.24810 -0.65638 -0.01561 0.38840 2.00547
##
## Random effects:
## Groups Name Variance Std.Dev.
## AntTag (Intercept) 0.2701 0.5197
## Residual 0.4205 0.6484
## Number of obs: 102, groups: AntTag, 39
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) 6.6409538 0.7302508
## ExpTemp17.5 0.0638552 0.3743724
## ExpTemp21.5 0.1757457 0.3858516
## AccTemp17.5 -0.4219128 0.4812347
## AccTemp21.5 -0.7958619 0.4840122
## SpeciesSagmariasus verreauxi 0.3198782 0.4749665
## SexMale -0.5882381 0.2378635
## Weight_g -0.0009436 0.0005936
## ExpTemp17.5:AccTemp17.5 0.2281950 0.5297265
## ExpTemp21.5:AccTemp17.5 0.9210081 0.5539125
## ExpTemp17.5:AccTemp21.5 0.0653483 0.5461824
## ExpTemp21.5:AccTemp21.5 0.7177046 0.5377035
## ExpTemp17.5:SpeciesSagmariasus verreauxi -0.7078266 0.5576536
## ExpTemp21.5:SpeciesSagmariasus verreauxi -0.0232111 0.5560048
## AccTemp17.5:SpeciesSagmariasus verreauxi 0.7204110 0.6782367
## AccTemp21.5:SpeciesSagmariasus verreauxi -0.0334204 0.6926546
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi 0.3420914 0.7884879
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi -1.0746224 0.7845790
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi 0.9286007 0.7920824
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi 0.0492534 0.7798092
## df t value
## (Intercept) 40.4346580 9.094
## ExpTemp17.5 56.9599755 0.171
## ExpTemp21.5 64.4586306 0.455
## AccTemp17.5 78.0294532 -0.877
## AccTemp21.5 77.4559012 -1.644
## SpeciesSagmariasus verreauxi 82.0967265 0.673
## SexMale 35.0921021 -2.473
## Weight_g 32.9412131 -1.590
## ExpTemp17.5:AccTemp17.5 57.1018556 0.431
## ExpTemp21.5:AccTemp17.5 61.7928656 1.663
## ExpTemp17.5:AccTemp21.5 58.1987799 0.120
## ExpTemp21.5:AccTemp21.5 60.7496177 1.335
## ExpTemp17.5:SpeciesSagmariasus verreauxi 64.3762657 -1.269
## ExpTemp21.5:SpeciesSagmariasus verreauxi 63.0626043 -0.042
## AccTemp17.5:SpeciesSagmariasus verreauxi 80.3470632 1.062
## AccTemp21.5:SpeciesSagmariasus verreauxi 80.0995022 -0.048
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi 63.7887733 0.434
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi 62.4663915 -1.370
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi 61.7684417 1.172
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi 60.6782422 0.063
## Pr(>|t|)
## (Intercept) 2.51e-11 ***
## ExpTemp17.5 0.8652
## ExpTemp21.5 0.6503
## AccTemp17.5 0.3833
## AccTemp21.5 0.1042
## SpeciesSagmariasus verreauxi 0.5025
## SexMale 0.0184 *
## Weight_g 0.1215
## ExpTemp17.5:AccTemp17.5 0.6683
## ExpTemp21.5:AccTemp17.5 0.1014
## ExpTemp17.5:AccTemp21.5 0.9052
## ExpTemp21.5:AccTemp21.5 0.1869
## ExpTemp17.5:SpeciesSagmariasus verreauxi 0.2089
## ExpTemp21.5:SpeciesSagmariasus verreauxi 0.9668
## AccTemp17.5:SpeciesSagmariasus verreauxi 0.2913
## AccTemp21.5:SpeciesSagmariasus verreauxi 0.9616
## ExpTemp17.5:AccTemp17.5:SpeciesSagmariasus verreauxi 0.6659
## ExpTemp21.5:AccTemp17.5:SpeciesSagmariasus verreauxi 0.1757
## ExpTemp17.5:AccTemp21.5:SpeciesSagmariasus verreauxi 0.2456
## ExpTemp21.5:AccTemp21.5:SpeciesSagmariasus verreauxi 0.9498
## ---
## 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 6.1260 3.06299 2 61.364 7.2848 0.00145 **
## AccTemp 2.1503 1.07515 2 31.972 2.5570 0.09329 .
## Species 0.9881 0.98809 1 31.904 2.3500 0.13514
## Sex 2.5715 2.57146 1 35.092 6.1158 0.01839 *
## Weight_g 1.0624 1.06238 1 32.941 2.5267 0.12149
## ExpTemp:AccTemp 1.7055 0.42637 4 61.170 1.0140 0.40727
## ExpTemp:Species 0.6085 0.30427 2 61.469 0.7237 0.48906
## AccTemp:Species 0.3554 0.17771 2 32.128 0.4227 0.65889
## ExpTemp:AccTemp:Species 1.9973 0.49932 4 61.365 1.1876 0.32526
## ---
## 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) = 7.28, p = 0.001; Eta2 (partial) = 0.19, 90% CI [0.05, 0.32])
## - The main effect of AccTemp is statistically not significant and medium (F(2) = 2.56, p = 0.093; Eta2 (partial) = 0.14, 90% CI [0.00, 0.31])
## - The main effect of Species is statistically not significant and medium (F(1) = 2.35, p = 0.135; Eta2 (partial) = 0.07, 90% CI [0.00, 0.24])
## - The main effect of Sex is statistically significant and large (F(1) = 6.12, p = 0.018; Eta2 (partial) = 0.15, 90% CI [0.01, 0.33])
## - The main effect of Weight_g is statistically not significant and medium (F(1) = 2.53, p = 0.121; Eta2 (partial) = 0.07, 90% CI [0.00, 0.24])
## - The interaction between ExpTemp and AccTemp is statistically not significant and medium (F(4) = 1.01, p = 0.407; Eta2 (partial) = 0.06, 90% CI [0.00, 0.13])
## - The interaction between ExpTemp and Species is statistically not significant and small (F(2) = 0.72, p = 0.489; Eta2 (partial) = 0.02, 90% CI [0.00, 0.10])
## - The interaction between AccTemp and Species is statistically not significant and small (F(2) = 0.42, p = 0.659; Eta2 (partial) = 0.03, 90% CI [0.00, 0.13])
## - The interaction between ExpTemp, AccTemp and Species is statistically not significant and medium (F(4) = 1.19, p = 0.325; Eta2 (partial) = 0.07, 90% CI [0.00, 0.15])
##
## Effect sizes were labelled following Field's (2013) recommendations.
# Table output from sjPlot
tab_model(complex_model)
value | |||
---|---|---|---|
Predictors | Estimates | CI | p |
(Intercept) | 6.64 | 5.21 – 8.07 | <0.001 |
ExpTemp [17.5] | 0.06 | -0.67 – 0.80 | 0.865 |
ExpTemp [21.5] | 0.18 | -0.58 – 0.93 | 0.649 |
AccTemp [17.5] | -0.42 | -1.37 – 0.52 | 0.381 |
AccTemp [21.5] | -0.80 | -1.74 – 0.15 | 0.100 |
Species [Sagmariasus verreauxi] |
0.32 | -0.61 – 1.25 | 0.501 |
Sex [Male] | -0.59 | -1.05 – -0.12 | 0.013 |
Weight_g | -0.00 | -0.00 – 0.00 | 0.112 |
ExpTemp [17.5] * AccTemp [17.5] |
0.23 | -0.81 – 1.27 | 0.667 |
ExpTemp [21.5] * AccTemp [17.5] |
0.92 | -0.16 – 2.01 | 0.096 |
ExpTemp [17.5] * AccTemp [21.5] |
0.07 | -1.01 – 1.14 | 0.905 |
ExpTemp [21.5] * AccTemp [21.5] |
0.72 | -0.34 – 1.77 | 0.182 |
ExpTemp [17.5] * Species [Sagmariasus verreauxi] |
-0.71 | -1.80 – 0.39 | 0.204 |
ExpTemp [21.5] * Species [Sagmariasus verreauxi] |
-0.02 | -1.11 – 1.07 | 0.967 |
AccTemp [17.5] * Species [Sagmariasus verreauxi] |
0.72 | -0.61 – 2.05 | 0.288 |
AccTemp [21.5] * Species [Sagmariasus verreauxi] |
-0.03 | -1.39 – 1.32 | 0.962 |
(ExpTemp [17.5] * AccTemp [17.5]) * Species [Sagmariasus verreauxi] |
0.34 | -1.20 – 1.89 | 0.664 |
(ExpTemp [21.5] * AccTemp [17.5]) * Species [Sagmariasus verreauxi] |
-1.07 | -2.61 – 0.46 | 0.171 |
(ExpTemp [17.5] * AccTemp [21.5]) * Species [Sagmariasus verreauxi] |
0.93 | -0.62 – 2.48 | 0.241 |
(ExpTemp [21.5] * AccTemp [21.5]) * Species [Sagmariasus verreauxi] |
0.05 | -1.48 – 1.58 | 0.950 |
Random Effects | |||
σ2 | 0.42 | ||
τ00 AntTag | 0.27 | ||
ICC | 0.39 | ||
N AntTag | 39 | ||
Observations | 102 | ||
Marginal R2 / Conditional R2 | 0.318 / 0.585 |
(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 -119.47 282.94
## (1 | AntTag) 0 21 -124.03 290.07 9.1361 1 0.002506 **
## ---
## 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)
## ExpTemp:AccTemp:Species 1 1.9973 0.4993 4 61.365 1.1876 0.325264
## AccTemp:Species 2 0.3838 0.1919 2 33.174 0.4240 0.657893
## ExpTemp:AccTemp 3 1.5324 0.3831 4 61.488 0.8353 0.508013
## ExpTemp:Species 4 0.7336 0.3668 2 63.546 0.7429 0.479815
## Weight_g 5 1.1427 1.1427 1 33.799 2.2915 0.139383
## Species 6 0.9036 0.9036 1 33.845 1.8156 0.186794
## AccTemp 7 1.9503 0.9751 2 34.091 1.9729 0.154608
## ExpTemp 0 5.8724 2.9362 2 64.236 5.9796 0.004154
## Sex 0 3.9604 3.9604 1 37.104 8.0655 0.007282
##
## ExpTemp:AccTemp:Species
## AccTemp:Species
## ExpTemp:AccTemp
## ExpTemp:Species
## Weight_g
## Species
## AccTemp
## ExpTemp **
## Sex **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## value ~ ExpTemp + Sex + (1 | AntTag)
# Extract the model that step found:
final_model <- get_model(step_result)
# Remodel using Maximum likelihood estimation
final_model <- lmer(value ~ ExpTemp +
Sex +
(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 + Sex + (1 | AntTag)
## Data: sub_trait
##
## REML criterion at convergence: 263.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9734 -0.6440 -0.1214 0.4985 2.2912
##
## Random effects:
## Groups Name Variance Std.Dev.
## AntTag (Intercept) 0.3670 0.6058
## Residual 0.5061 0.7114
## Number of obs: 102, groups: AntTag, 39
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.5595 0.2311 53.8508 24.060 < 2e-16 ***
## ExpTemp17.5 0.0266 0.1762 62.6949 0.151 0.88048
## ExpTemp21.5 0.5290 0.1744 62.2701 3.033 0.00353 **
## SexMale -0.7084 0.2563 35.1966 -2.764 0.00903 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ET17.5 ET21.5
## ExpTemp17.5 -0.364
## ExpTemp21.5 -0.376 0.491
## SexMale -0.737 -0.004 0.009
# 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 5.8473 2.9236 2 62.315 5.7770 0.004994 **
## Sex 3.8657 3.8657 1 35.197 7.6385 0.009029 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
report(anova(final_model))
## The ANOVA suggests that:
##
## - The main effect of ExpTemp is statistically significant and large (F(2) = 5.78, p = 0.005; Eta2 (partial) = 0.16, 90% CI [0.03, 0.29])
## - The main effect of Sex is statistically significant and large (F(1) = 7.64, p = 0.009; Eta2 (partial) = 0.18, 90% CI [0.03, 0.36])
##
## Effect sizes were labelled following Field's (2013) recommendations.
# Table output from sjPlot
tab_model(final_model)
value | |||
---|---|---|---|
Predictors | Estimates | CI | p |
(Intercept) | 5.56 | 5.11 – 6.01 | <0.001 |
ExpTemp [17.5] | 0.03 | -0.32 – 0.37 | 0.880 |
ExpTemp [21.5] | 0.53 | 0.19 – 0.87 | 0.002 |
Sex [Male] | -0.71 | -1.21 – -0.21 | 0.006 |
Random Effects | |||
σ2 | 0.51 | ||
τ00 AntTag | 0.37 | ||
ICC | 0.42 | ||
N AntTag | 39 | ||
Observations | 102 | ||
Marginal R2 / Conditional R2 | 0.164 / 0.515 |
# 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 Sex (formula: value ~ ExpTemp + Sex). The model included AntTag as random effect (formula: ~1 | AntTag). The model's total explanatory power is substantial (conditional R2 = 0.52) and the part related to the fixed effects alone (marginal R2) is of 0.16. The model's intercept, corresponding to ExpTemp = 14 and Sex = Female, is at 5.56 (95% CI [5.11, 6.01], t(96) = 24.06, p < .001). Within this model:
##
## - The effect of ExpTemp [17.5] is statistically non-significant and positive (beta = 0.03, 95% CI [-0.32, 0.37], t(96) = 0.15, p = 0.880; Std. beta = 0.03, 95% CI [-0.32, 0.37])
## - The effect of ExpTemp [21.5] is statistically significant and positive (beta = 0.53, 95% CI [0.19, 0.87], t(96) = 3.03, p = 0.002; Std. beta = 0.53, 95% CI [0.19, 0.87])
## - The effect of Sex [Male] is statistically significant and negative (beta = -0.71, 95% CI [-1.21, -0.21], t(96) = -2.76, p = 0.006; Std. beta = -0.71, 95% CI [-1.21, -0.21])
##
## 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, paste("Model_summary_", trait, ".csv", sep = ""))
# Print table
kable(lme_table,
format = "html",
row.names = TRUE,
caption = "Model summary for main effects",
col.names = c("Sum of sqares", "Mean squares", "Den df", "F(df)", "p value", "Partical squared eta [CI 90%]")) %>%
kable_classic(full_width = F, html_font = "Cambria")
Sum of sqares | Mean squares | Den df | F(df) | p value | Partical squared eta [CI 90%] | |
---|---|---|---|---|---|---|
ExpTemp | 5.8 | 2.9 | 62.3 | 5.78 (2) | 0.004994 | 0.16 [0.03-0.29] |
Sex | 3.9 | 3.9 | 35.2 | 7.64 (1) | 0.009029 | 0.18 [0.03-0.36] |
# 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$ExpTemp)
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)