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
# R document settings
# Echo = True means show all code chunks
knitr::opts_chunk$set(echo = TRUE)
# Select trait
trait <- "MMR"
# Set working directory of executing R. script
setwd(paste("C:/Users/oelle/Documents/Meine Dokumente/Research/Projects/Project_RockLobster/Temperature acclimation/Aerobic scope/Statistics/LinearMixedEffectModelling/", trait, "/", sep = ""))
# Load tidy raw data in long format
mo2_data_long <- read.csv("../../MO2_data_Tidy_RawData_long.csv")
# Shorten column names
names(mo2_data_long)[c(5,6)] <- c("AccTemp", "ExpTemp")
# Define fixed 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 == "MMR")
# Plot data versus a range of distrubutions
descdist(sub_trait$value, boot = 1000)
## summary statistics
## ------
## min: 62.81 max: 160.02
## median: 108.025
## mean: 108.1056
## estimated sd: 22.05236
## estimated skewness: -0.01523508
## estimated kurtosis: 2.469454
# Fit various distributions
nd <- fitdist(sub_trait$value, "norm")
## $start.arg
## $start.arg$mean
## [1] 108.1056
##
## $start.arg$sd
## [1] 21.95003
##
##
## $fix.arg
## NULL
gd <- fitdist(sub_trait$value, "gamma")
## $start.arg
## $start.arg$shape
## [1] 24.25641
##
## $start.arg$rate
## [1] 0.2243769
##
##
## $fix.arg
## NULL
wd <- fitdist(sub_trait$value, "weibull")
## $start.arg
## $start.arg$shape
## [1] 5.624781
##
## $start.arg$scale
## [1] 117.0993
##
##
## $fix.arg
## NULL
lnd <- fitdist(sub_trait$value, "lnorm")
## $start.arg
## $start.arg$meanlog
## [1] 4.661329
##
## $start.arg$sdlog
## [1] 0.2123517
##
##
## $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.04961111 0.06203287 0.04928462 0.07189584
## Cramer-von Mises statistic 0.02731938 0.07516722 0.02481451 0.12058465
## Anderson-Darling statistic 0.20138464 0.48593740 0.20839501 0.76875941
##
## Goodness-of-fit criteria
## Normal gamma Weibull LogNormal
## Akaike's Information Criterion 977.6647 979.6881 978.6070 982.6433
## Bayesian Information Criterion 983.0290 985.0524 983.9713 988.0075
# Gplot
ggplot(sub_trait, aes(x = ExpTemp, y = value, fill = as.factor(AccTemp), alpha = .5)) +
geom_boxplot(outlier.size = 0) +
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 | 101.7 | 54 | 20.2 | 5.5 | 96.1 | 107.2 | 62.8 | 148.7 |
Sagmariasus verreauxi | 114.6 | 54 | 22.1 | 6.0 | 108.5 | 120.6 | 66.9 | 160.0 |
# Summary statistics
substrait_stats <- ddply(sub_trait, c("Species", "ExpTemp"), summarise,
mean = mean(value, na.rm = TRUE),
n = length(value),
sd = sd(value, na.rm = TRUE),
ci = qt(0.975, n) * sd / sqrt(n),
lower_ci = mean-ci,
upper_ci = mean+ci,
min = min(value),
max = max(value))
# Round values
substrait_stats <- substrait_stats %>% mutate_if(is.numeric, round, digit = 3)
# Display results
substrait_stats %>% kbl(digits = 1, caption = trait) %>% kable_classic()
Species | ExpTemp | mean | n | sd | ci | lower_ci | upper_ci | min | max |
---|---|---|---|---|---|---|---|---|---|
Jasus edwardsii | 14 | 90.6 | 18 | 16.4 | 8.1 | 82.5 | 98.8 | 62.8 | 118.9 |
Jasus edwardsii | 17.5 | 100.3 | 18 | 19.3 | 9.6 | 90.7 | 109.9 | 64.3 | 144.3 |
Jasus edwardsii | 21.5 | 114.0 | 18 | 18.3 | 9.1 | 104.9 | 123.1 | 80.0 | 148.7 |
Sagmariasus verreauxi | 14 | 96.3 | 18 | 15.2 | 7.5 | 88.8 | 103.9 | 66.9 | 126.4 |
Sagmariasus verreauxi | 17.5 | 116.0 | 18 | 16.6 | 8.2 | 107.8 | 124.3 | 70.8 | 137.7 |
Sagmariasus verreauxi | 21.5 | 131.3 | 18 | 19.1 | 9.5 | 121.9 | 140.8 | 74.0 | 160.0 |
# Start with the most complex model
complex_model <- lmer(value ~ ExpTemp *
AccTemp *
Species +
Sex +
Weight_g +
(1|AntTag), data = sub_trait, REML = FALSE)
# Model summary
anova(complex_model)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ExpTemp 12645.9 6323.0 2 70.370 72.1287 < 2.2e-16 ***
## AccTemp 677.7 338.9 2 37.762 3.8657 0.029693 *
## Species 1036.8 1036.8 1 37.929 11.8275 0.001433 **
## Sex 11.6 11.6 1 39.103 0.1321 0.718198
## Weight_g 447.8 447.8 1 44.202 5.1081 0.028791 *
## ExpTemp:AccTemp 1001.4 250.4 4 70.436 2.8559 0.029724 *
## ExpTemp:Species 526.4 263.2 2 70.844 3.0024 0.056028 .
## AccTemp:Species 231.6 115.8 2 38.031 1.3209 0.278871
## ExpTemp:AccTemp:Species 186.0 46.5 4 70.900 0.5305 0.713710
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(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 -427.03 898.06
## (1 | AntTag) 0 21 -444.75 931.50 35.433 1 2.639e-09 ***
## ---
## 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 11.58 11.58 1 39.103 0.1321
## ExpTemp:AccTemp:Species 2 188.75 47.19 4 70.866 0.5389
## AccTemp:Species 3 237.97 118.98 2 38.207 1.3193
## ExpTemp:Species 4 525.24 262.62 2 70.869 2.9138
## Species 0 995.55 995.55 1 38.184 10.2588
## Weight_g 0 458.26 458.26 1 44.916 4.7222
## ExpTemp:AccTemp 0 1013.40 253.35 4 70.386 2.6107
## Pr(>F)
## Sex 0.718198
## ExpTemp:AccTemp:Species 0.707647
## AccTemp:Species 0.279222
## ExpTemp:Species 0.060799 .
## Species 0.002743 **
## Weight_g 0.035089 *
## ExpTemp:AccTemp 0.042626 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## value ~ ExpTemp + AccTemp + Species + Weight_g + (1 | AntTag) + ExpTemp:AccTemp
# Extract the model that step found:
final_model <- get_model(step_result)
final_model
## Linear mixed model fit by maximum likelihood ['lmerModLmerTest']
## Formula: value ~ ExpTemp + AccTemp + Species + Weight_g + (1 | AntTag) +
## ExpTemp:AccTemp
## Data: sub_trait
## AIC BIC logLik deviance df.resid
## 890.5019 925.3696 -432.2510 864.5019 95
## Random effects:
## Groups Name Std.Dev.
## AntTag (Intercept) 12.207
## Residual 9.851
## Number of obs: 108, groups: AntTag, 39
## Fixed Effects:
## (Intercept) ExpTemp17.5
## 120.35372 14.62665
## ExpTemp21.5 AccTemp17.5
## 21.56607 2.21054
## AccTemp21.5 SpeciesSagmariasus verreauxi
## -18.30465 14.16481
## Weight_g ExpTemp17.5:AccTemp17.5
## -0.02478 -8.30831
## ExpTemp21.5:AccTemp17.5 ExpTemp17.5:AccTemp21.5
## 2.53309 4.22518
## ExpTemp21.5:AccTemp21.5
## 13.86326
# Remodel using Maximum likelihood estimation
final_model2 <- lmer(value ~ ExpTemp *
AccTemp *
Species +
Weight_g +
(1|AntTag), data = sub_trait, REML = TRUE,
contrasts = list(
ExpTemp = "contr.sum",
AccTemp = "contr.sum",
Species = "contr.sum"
))
# Model summary
summary(final_model)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: value ~ ExpTemp + AccTemp + Species + Weight_g + (1 | AntTag) +
## ExpTemp:AccTemp
## Data: sub_trait
##
## AIC BIC logLik deviance df.resid
## 890.5 925.4 -432.3 864.5 95
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3632 -0.5462 0.0977 0.5889 1.8819
##
## Random effects:
## Groups Name Variance Std.Dev.
## AntTag (Intercept) 149.02 12.207
## Residual 97.04 9.851
## Number of obs: 108, groups: AntTag, 39
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 120.35372 13.48976 47.90078 8.922 9.41e-12
## ExpTemp17.5 14.62665 4.10981 71.31597 3.559 0.000667
## ExpTemp21.5 21.56607 4.19779 73.96192 5.137 2.19e-06
## AccTemp17.5 2.21054 6.25427 69.52310 0.353 0.724824
## AccTemp21.5 -18.30465 6.33446 67.45980 -2.890 0.005181
## SpeciesSagmariasus verreauxi 14.16481 4.42245 38.18352 3.203 0.002743
## Weight_g -0.02478 0.01140 44.91586 -2.173 0.035089
## ExpTemp17.5:AccTemp17.5 -8.30831 5.81591 71.46023 -1.429 0.157490
## ExpTemp21.5:AccTemp17.5 2.53309 5.87350 72.60468 0.431 0.667546
## ExpTemp17.5:AccTemp21.5 4.22518 5.75244 70.08435 0.735 0.465094
## ExpTemp21.5:AccTemp21.5 13.86326 5.81451 71.46809 2.384 0.019771
##
## (Intercept) ***
## ExpTemp17.5 ***
## ExpTemp21.5 ***
## AccTemp17.5
## AccTemp21.5 **
## SpeciesSagmariasus verreauxi **
## Weight_g *
## ExpTemp17.5:AccTemp17.5
## ExpTemp21.5:AccTemp17.5
## ExpTemp17.5:AccTemp21.5
## ExpTemp21.5:AccTemp21.5 *
## ---
## 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 SpcsSv Wght_g ET17.5:AT1
## ExpTemp17.5 -0.127
## ExpTemp21.5 -0.150 0.511
## AccTemp17.5 -0.245 0.326 0.335
## AccTemp21.5 -0.309 0.320 0.331 0.489
## SpcsSgmrssv -0.025 0.011 0.015 -0.019 -0.009
## Weight_g -0.931 -0.027 -0.008 0.020 0.089 -0.148
## ET17.5:AT17 0.064 -0.707 -0.361 -0.460 -0.224 -0.007 0.046
## ET21.5:AT17 0.103 -0.365 -0.715 -0.467 -0.236 -0.007 0.010 0.516
## ET17.5:AT21 0.077 -0.715 -0.365 -0.232 -0.449 -0.010 0.034 0.506
## ET21.5:AT21 0.095 -0.369 -0.722 -0.242 -0.457 -0.013 0.021 0.261
## ET21.5:AT1 ET17.5:AT2
## ExpTemp17.5
## ExpTemp21.5
## AccTemp17.5
## AccTemp21.5
## SpcsSgmrssv
## Weight_g
## ET17.5:AT17
## ET21.5:AT17
## ET17.5:AT21 0.261
## ET21.5:AT21 0.516 0.506
# Compare model with and without interaction with species
anova(final_model, final_model2)
## refitting model(s) with ML (instead of REML)
## Data: sub_trait
## Models:
## final_model: value ~ ExpTemp + AccTemp + Species + Weight_g + (1 | AntTag) + ExpTemp:AccTemp
## final_model2: value ~ ExpTemp * AccTemp * Species + Weight_g + (1 | AntTag)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## final_model 13 890.50 925.37 -432.25 864.50
## final_model2 21 896.19 952.52 -427.10 854.19 10.308 8 0.2441
## --> model with species interaction quite similar and will therefore retain interaction. Contrast analysis also showed species dependent effects of acclimation (see below)
# Replace final model
final_model <- final_model2
# Model summary
#summary(final_model)
# 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 12628.9 6314.4 2 58.485 59.7084 7.446e-15 ***
## AccTemp 655.7 327.8 2 31.272 3.0999 0.05913 .
## Species 1008.9 1008.9 1 31.212 9.5402 0.00420 **
## Weight_g 437.1 437.1 1 35.791 4.1333 0.04952 *
## ExpTemp:AccTemp 1004.3 251.1 4 58.539 2.3742 0.06241 .
## ExpTemp:Species 527.2 263.6 2 58.851 2.4927 0.09138 .
## AccTemp:Species 222.9 111.4 2 31.315 1.0538 0.36064
## ExpTemp:AccTemp:Species 188.6 47.2 4 58.897 0.4459 0.77491
## ---
## 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) = 59.71, p < .001; Eta2 (partial) = 0.67, 90% CI [0.55, 0.75])
## - The main effect of AccTemp is statistically not significant and large (F(2) = 3.10, p = 0.059; Eta2 (partial) = 0.17, 90% CI [0.00, 0.34])
## - The main effect of Species is statistically significant and large (F(1) = 9.54, p = 0.004; Eta2 (partial) = 0.23, 90% CI [0.05, 0.43])
## - The main effect of Weight_g is statistically significant and medium (F(1) = 4.13, p = 0.050; Eta2 (partial) = 0.10, 90% CI [1.25e-04, 0.28])
## - The interaction between ExpTemp and AccTemp is statistically not significant and medium (F(4) = 2.37, p = 0.062; Eta2 (partial) = 0.14, 90% CI [0.00, 0.27])
## - The interaction between ExpTemp and Species is statistically not significant and medium (F(2) = 2.49, p = 0.091; Eta2 (partial) = 0.08, 90% CI [0.00, 0.19])
## - The interaction between AccTemp and Species is statistically not significant and medium (F(2) = 1.05, p = 0.361; Eta2 (partial) = 0.06, 90% CI [0.00, 0.21])
## - The interaction between ExpTemp, AccTemp and Species is statistically not significant and small (F(4) = 0.45, p = 0.775; Eta2 (partial) = 0.03, 90% CI [0.00, 0.07])
##
## Effect sizes were labelled following Field's (2013) recommendations.
# Table output from sjPlot
tab_model(final_model, digits.p = 3)
value | |||
---|---|---|---|
Predictors | Estimates | CI | p |
(Intercept) | 135.95 | 108.73 – 163.17 | <0.001 |
ExpTemp1 | -13.51 | -16.32 – -10.70 | <0.001 |
ExpTemp2 | -0.09 | -2.85 – 2.68 | 0.952 |
AccTemp1 | 4.03 | -2.36 – 10.42 | 0.216 |
AccTemp2 | 4.31 | -2.13 – 10.76 | 0.189 |
Species1 | -7.27 | -11.88 – -2.66 | 0.002 |
Weight_g | -0.03 | -0.05 – -0.00 | 0.042 |
ExpTemp1 * AccTemp1 | 1.47 | -2.53 – 5.47 | 0.471 |
ExpTemp2 * AccTemp1 | 2.77 | -1.16 – 6.69 | 0.168 |
ExpTemp1 * AccTemp2 | 3.11 | -0.87 – 7.10 | 0.126 |
ExpTemp2 * AccTemp2 | -3.61 | -7.52 – 0.30 | 0.071 |
ExpTemp1 * Species1 | 2.98 | 0.17 – 5.79 | 0.038 |
ExpTemp2 * Species1 | -0.47 | -3.25 – 2.30 | 0.739 |
AccTemp1 * Species1 | 4.25 | -2.17 – 10.66 | 0.195 |
AccTemp2 * Species1 | -0.16 | -6.61 – 6.28 | 0.961 |
ExpTemp1 * AccTemp1 * Species1 |
-1.00 | -5.01 – 3.02 | 0.626 |
ExpTemp2 * AccTemp1 * Species1 |
0.57 | -3.37 – 4.50 | 0.777 |
ExpTemp1 * AccTemp2 * Species1 |
1.65 | -2.35 – 5.64 | 0.419 |
ExpTemp2 * AccTemp2 * Species1 |
-2.44 | -6.38 – 1.49 | 0.223 |
Random Effects | |||
σ2 | 105.75 | ||
τ00 AntTag | 169.32 | ||
ICC | 0.62 | ||
N AntTag | 39 | ||
Observations | 108 | ||
Marginal R2 / Conditional R2 | 0.470 / 0.796 |
# 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, Species and Weight_g (formula: value ~ ExpTemp * AccTemp * Species + Weight_g). The model included AntTag as random effect (formula: ~1 | AntTag). The model's total explanatory power is substantial (conditional R2 = 0.80) and the part related to the fixed effects alone (marginal R2) is of 0.47. The model's intercept, corresponding to ExpTemp = 14, AccTemp = 14, Species = Jasus edwardsii and Weight_g = 0, is at 135.95 (95% CI [108.73, 163.17], t(87) = 9.79, p < .001). Within this model:
##
## - The effect of ExpTemp1 is statistically significant and negative (beta = -13.51, 95% CI [-16.32, -10.70], t(87) = -9.43, p < .001; Std. beta = -0.61, 95% CI [-0.74, -0.49])
## - The effect of ExpTemp2 is statistically non-significant and negative (beta = -0.09, 95% CI [-2.85, 2.68], t(87) = -0.06, p = 0.952; Std. beta = -3.89e-03, 95% CI [-0.13, 0.12])
## - The effect of AccTemp1 is statistically non-significant and positive (beta = 4.03, 95% CI [-2.36, 10.42], t(87) = 1.24, p = 0.216; Std. beta = 0.18, 95% CI [-0.11, 0.47])
## - The effect of AccTemp2 is statistically non-significant and positive (beta = 4.31, 95% CI [-2.13, 10.76], t(87) = 1.31, p = 0.189; Std. beta = 0.20, 95% CI [-0.10, 0.49])
## - The effect of Species1 is statistically significant and negative (beta = -7.27, 95% CI [-11.88, -2.66], t(87) = -3.09, p = 0.002; Std. beta = -0.33, 95% CI [-0.54, -0.12])
## - The effect of Weight_g is statistically significant and negative (beta = -0.03, 95% CI [-0.05, -9.03e-04], t(87) = -2.03, p = 0.042; Std. beta = -0.22, 95% CI [-0.43, -7.90e-03])
## - The interaction effect of AccTemp1 on ExpTemp1 is statistically non-significant and positive (beta = 1.47, 95% CI [-2.53, 5.47], t(87) = 0.72, p = 0.471; Std. beta = 0.07, 95% CI [-0.11, 0.25])
## - The interaction effect of AccTemp1 on ExpTemp2 is statistically non-significant and positive (beta = 2.77, 95% CI [-1.16, 6.69], t(87) = 1.38, p = 0.168; Std. beta = 0.13, 95% CI [-0.05, 0.30])
## - The interaction effect of AccTemp2 on ExpTemp1 is statistically non-significant and positive (beta = 3.11, 95% CI [-0.87, 7.10], t(87) = 1.53, p = 0.126; Std. beta = 0.14, 95% CI [-0.04, 0.32])
## - The interaction effect of AccTemp2 on ExpTemp2 is statistically non-significant and negative (beta = -3.61, 95% CI [-7.52, 0.30], t(87) = -1.81, p = 0.071; Std. beta = -0.16, 95% CI [-0.34, 0.01])
## - The interaction effect of Species1 on ExpTemp1 is statistically significant and positive (beta = 2.98, 95% CI [0.17, 5.79], t(87) = 2.08, p = 0.038; Std. beta = 0.14, 95% CI [7.62e-03, 0.26])
## - The interaction effect of Species1 on ExpTemp2 is statistically non-significant and negative (beta = -0.47, 95% CI [-3.25, 2.30], t(87) = -0.33, p = 0.739; Std. beta = -0.02, 95% CI [-0.15, 0.10])
## - The interaction effect of Species1 on AccTemp1 is statistically non-significant and positive (beta = 4.25, 95% CI [-2.17, 10.66], t(87) = 1.30, p = 0.195; Std. beta = 0.19, 95% CI [-0.10, 0.48])
## - The interaction effect of Species1 on AccTemp2 is statistically non-significant and negative (beta = -0.16, 95% CI [-6.61, 6.28], t(87) = -0.05, p = 0.961; Std. beta = -7.36e-03, 95% CI [-0.30, 0.28])
## - The interaction effect of Species1 on ExpTemp1 * AccTemp1 is statistically non-significant and negative (beta = -1.00, 95% CI [-5.01, 3.02], t(87) = -0.49, p = 0.626; Std. beta = -0.05, 95% CI [-0.23, 0.14])
## - The interaction effect of Species1 on ExpTemp2 * AccTemp1 is statistically non-significant and positive (beta = 0.57, 95% CI [-3.37, 4.50], t(87) = 0.28, p = 0.777; Std. beta = 0.03, 95% CI [-0.15, 0.20])
## - The interaction effect of Species1 on ExpTemp1 * AccTemp2 is statistically non-significant and positive (beta = 1.65, 95% CI [-2.35, 5.64], t(87) = 0.81, p = 0.419; Std. beta = 0.07, 95% CI [-0.11, 0.26])
## - The interaction effect of Species1 on ExpTemp2 * AccTemp2 is statistically non-significant and negative (beta = -2.44, 95% CI [-6.38, 1.49], t(87) = -1.22, p = 0.223; Std. beta = -0.11, 95% CI [-0.29, 0.07])
##
## 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 | 12628.9 | 6314.4 | 58.5 | 59.71 (2) | 7.446e-15 | 0.67 [0.55-0.75] |
AccTemp | 655.7 | 327.8 | 31.3 | 3.1 (2) | 5.913e-02 | 0.17 [0-0.34] |
Species | 1008.9 | 1008.9 | 31.2 | 9.54 (1) | 4.200e-03 | 0.23 [0.05-0.43] |
Weight_g | 437.1 | 437.1 | 35.8 | 4.13 (1) | 4.952e-02 | 0.1 [0-0.28] |
ExpTemp:AccTemp | 1004.3 | 251.1 | 58.5 | 2.37 (4) | 6.241e-02 | 0.14 [0-0.27] |
ExpTemp:Species | 527.2 | 263.6 | 58.9 | 2.49 (2) | 9.138e-02 | 0.08 [0-0.19] |
AccTemp:Species | 222.9 | 111.4 | 31.3 | 1.05 (2) | 3.606e-01 | 0.06 [0-0.21] |
ExpTemp:AccTemp:Species | 188.6 | 47.2 | 58.9 | 0.45 (4) | 7.749e-01 | 0.03 [0-0.07] |
# 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_model2, 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 | 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 | -0.1587081 | -0.0166972 | 1 | No | MMR |
14 - 21.5 | 14 | Jasus edwardsii | 26.4164781 | 2.7533828 | 0.024 | Yes | MMR |
17.5 - 21.5 | 14 | Jasus edwardsii | 26.5751862 | 2.7690258 | 0.023 | Yes | MMR |
14 - 17.5 | 17.5 | Jasus edwardsii | 13.5136817 | 1.4264036 | 0.478 | No | MMR |
14 - 21.5 | 17.5 | Jasus edwardsii | 21.3262011 | 2.2253668 | 0.091 | No | MMR |
17.5 - 21.5 | 17.5 | Jasus edwardsii | 7.8125194 | 0.8119432 | 1 | No | MMR |
14 - 17.5 | 21.5 | Jasus edwardsii | -0.9730142 | -0.1029332 | 1 | No | MMR |
14 - 21.5 | 21.5 | Jasus edwardsii | 14.3868949 | 1.5113248 | 0.41 | No | MMR |
17.5 - 21.5 | 21.5 | Jasus edwardsii | 15.3599092 | 1.6011369 | 0.346 | No | MMR |
14 - 17.5 | 14 | Sagmariasus verreauxi | -3.6852236 | -0.3968126 | 1 | No | MMR |
14 - 21.5 | 14 | Sagmariasus verreauxi | 10.4543134 | 1.1089070 | 0.817 | No | MMR |
17.5 - 21.5 | 14 | Sagmariasus verreauxi | 14.1395370 | 1.5003347 | 0.418 | No | MMR |
14 - 17.5 | 17.5 | Sagmariasus verreauxi | -1.3256659 | -0.1420389 | 1 | No | MMR |
14 - 21.5 | 17.5 | Sagmariasus verreauxi | 7.2819029 | 0.7702734 | 1 | No | MMR |
17.5 - 21.5 | 17.5 | Sagmariasus verreauxi | 8.6075688 | 0.9104056 | 1 | No | MMR |
14 - 17.5 | 21.5 | Sagmariasus verreauxi | -9.0530080 | -0.9696316 | 1 | No | MMR |
14 - 21.5 | 21.5 | Sagmariasus verreauxi | -5.5780971 | -0.5900463 | 1 | No | MMR |
17.5 - 21.5 | 21.5 | Sagmariasus verreauxi | 3.4749110 | 0.3671210 | 1 | No | MMR |
# 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", "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.8350000 | -2.1617550 | 0.105 | No | MMR |
14 - 21.5 | 14 | Jasus edwardsii | -17.3390272 | -2.7825252 | 0.021 | Yes | MMR |
17.5 - 21.5 | 14 | Jasus edwardsii | -4.5040272 | -0.7227954 | 1 | No | MMR |
14 - 17.5 | 17.5 | Jasus edwardsii | 0.8373898 | 0.1407746 | 1 | No | MMR |
14 - 21.5 | 17.5 | Jasus edwardsii | -18.1533333 | -3.0575036 | 0.01 | Yes | MMR |
17.5 - 21.5 | 17.5 | Jasus edwardsii | -18.9907231 | -3.1925530 | 0.007 | Yes | MMR |
14 - 17.5 | 21.5 | Jasus edwardsii | -17.9252770 | -3.0189098 | 0.011 | Yes | MMR |
14 - 21.5 | 21.5 | Jasus edwardsii | -29.3686103 | -4.9461543 | 0 | Yes | MMR |
17.5 - 21.5 | 21.5 | Jasus edwardsii | -11.4433333 | -1.9273614 | 0.177 | No | MMR |
14 - 17.5 | 14 | Sagmariasus verreauxi | -16.6007584 | -2.6631218 | 0.03 | Yes | MMR |
14 - 21.5 | 14 | Sagmariasus verreauxi | -25.4524250 | -4.0831212 | 0 | Yes | MMR |
17.5 - 21.5 | 14 | Sagmariasus verreauxi | -8.8516667 | -1.4908558 | 0.425 | No | MMR |
14 - 17.5 | 17.5 | Sagmariasus verreauxi | -14.2412006 | -2.2708359 | 0.08 | No | MMR |
14 - 21.5 | 17.5 | Sagmariasus verreauxi | -30.8202095 | -4.9477974 | 0 | Yes | MMR |
17.5 - 21.5 | 17.5 | Sagmariasus verreauxi | -16.5790088 | -2.7756200 | 0.022 | Yes | MMR |
14 - 17.5 | 21.5 | Sagmariasus verreauxi | -19.7731688 | -3.3273515 | 0.005 | Yes | MMR |
14 - 21.5 | 21.5 | Sagmariasus verreauxi | -41.4848355 | -6.9809057 | 0 | Yes | MMR |
17.5 - 21.5 | 21.5 | Sagmariasus verreauxi | -21.7116667 | -3.6568215 | 0.002 | Yes | MMR |
# 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 | -8.580352 | -1.556536 | 0.125 | No | MMR |
Jasus edwardsii - Sagmariasus verreauxi | 17.5 | -15.477765 | -2.825278 | 0.007 | Yes | MMR |
Jasus edwardsii - Sagmariasus verreauxi | 21.5 | -19.545851 | -3.539815 | 0.001 | Yes | MMR |
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)