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 <- "AerobicScope"
# 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("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: 16.18 max: 121.76
## median: 84.69
## mean: 82.43444
## estimated sd: 20.76008
## estimated skewness: -0.7854203
## estimated kurtosis: 3.608739
# Fit various distributions
nd <- fitdist(sub_trait$value, "norm")
## $start.arg
## $start.arg$mean
## [1] 82.43444
##
## $start.arg$sd
## [1] 20.66374
##
##
## $fix.arg
## NULL
gd <- fitdist(sub_trait$value, "gamma")
## $start.arg
## $start.arg$shape
## [1] 15.91474
##
## $start.arg$rate
## [1] 0.1930593
##
##
## $fix.arg
## NULL
wd <- fitdist(sub_trait$value, "weibull")
## $start.arg
## $start.arg$shape
## [1] 3.702822
##
## $start.arg$scale
## [1] 92.17977
##
##
## $fix.arg
## NULL
lnd <- fitdist(sub_trait$value, "lnorm")
## $start.arg
## $start.arg$meanlog
## [1] 4.369264
##
## $start.arg$sdlog
## [1] 0.3225733
##
##
## $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.08510534 0.1354993 0.06386793 0.1597644
## Cramer-von Mises statistic 0.17227356 0.5082014 0.08967312 0.7525458
## Anderson-Darling statistic 1.11153080 3.0383902 0.66887090 4.4307737
##
## Goodness-of-fit criteria
## Normal gamma Weibull LogNormal
## Akaike's Information Criterion 964.6209 990.1547 959.4582 1009.864
## Bayesian Information Criterion 969.9852 995.5190 964.8224 1015.228
–> 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 | 75.6 | 54 | 19.5 | 5.3 | 70.3 | 81.0 | 23.7 | 120.5 |
Sagmariasus verreauxi | 89.2 | 54 | 19.9 | 5.4 | 83.8 | 94.6 | 16.2 | 121.8 |
# Summary statistics
substrait_stats <- ddply(sub_trait, c("Species", "ExpTemp"), summarise,
mean = mean(value, na.rm = TRUE),
n = length(value),
sd = sd(value, na.rm = TRUE),
ci = qt(0.975, n) * sd / sqrt(n),
lower_ci = mean-ci,
upper_ci = mean+ci,
min = min(value),
max = max(value))
# Round values
substrait_stats <- substrait_stats %>% mutate_if(is.numeric, round, digit = 3)
# Display results
substrait_stats %>% kbl(digits = 1, caption = trait) %>% kable_classic()
Species | ExpTemp | mean | n | sd | ci | lower_ci | upper_ci | min | max |
---|---|---|---|---|---|---|---|---|---|
Jasus edwardsii | 14 | 74.7 | 18 | 16.0 | 7.9 | 66.8 | 82.7 | 44.3 | 101.3 |
Jasus edwardsii | 17.5 | 75.4 | 18 | 24.5 | 12.1 | 63.2 | 87.5 | 23.7 | 120.5 |
Jasus edwardsii | 21.5 | 76.8 | 18 | 18.2 | 9.0 | 67.8 | 85.8 | 32.7 | 102.2 |
Sagmariasus verreauxi | 14 | 81.3 | 18 | 14.9 | 7.4 | 73.9 | 88.7 | 51.2 | 105.6 |
Sagmariasus verreauxi | 17.5 | 92.3 | 18 | 18.8 | 9.3 | 83.0 | 101.6 | 34.6 | 113.0 |
Sagmariasus verreauxi | 21.5 | 94.0 | 18 | 23.5 | 11.6 | 82.4 | 105.7 | 16.2 | 121.8 |
# REML has to be put to FALSE as models can only be compared using ML
# 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 613.03 306.51 2 71.161 2.3141 0.106251
## AccTemp 230.44 115.22 2 38.252 0.8699 0.427127
## Species 1311.38 1311.38 1 38.414 9.9006 0.003187 **
## Sex 127.91 127.91 1 39.627 0.9657 0.331718
## Weight_g 545.48 545.48 1 43.566 4.1183 0.048560 *
## ExpTemp:AccTemp 2222.29 555.57 4 71.210 4.1944 0.004169 **
## ExpTemp:Species 517.32 258.66 2 71.583 1.9528 0.149360
## AccTemp:Species 381.61 190.80 2 38.541 1.4405 0.249270
## ExpTemp:AccTemp:Species 213.04 53.26 4 71.614 0.4021 0.806518
## ---
## 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 -446.22 936.44
## (1 | AntTag) 0 21 -461.38 964.76 30.327 1 3.65e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## ExpTemp:AccTemp:Species 1 213.04 53.26 4 71.614 0.4021
## Sex 2 139.21 139.21 1 39.864 1.0277
## AccTemp:Species 3 364.98 182.49 2 38.701 1.3489
## ExpTemp:Species 4 520.12 260.06 2 71.446 1.9211
## Weight_g 5 409.26 409.26 1 44.395 2.8763
## Species 0 933.09 933.09 1 38.348 6.5437
## ExpTemp:AccTemp 0 2250.83 562.71 4 70.629 3.9462
## Pr(>F)
## ExpTemp:AccTemp:Species 0.806518
## Sex 0.316814
## AccTemp:Species 0.271448
## ExpTemp:Species 0.153938
## Weight_g 0.096893 .
## Species 0.014596 *
## ExpTemp:AccTemp 0.006007 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## value ~ ExpTemp + AccTemp + Species + (1 | AntTag) + ExpTemp:AccTemp
# Extract the model that step found:
final_model <- get_model(step_result)
# Model summary
anova(final_model)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ExpTemp 589.05 294.53 2 70.664 2.0655 0.134342
## AccTemp 186.95 93.48 2 38.297 0.6555 0.524886
## Species 933.09 933.09 1 38.348 6.5437 0.014596 *
## ExpTemp:AccTemp 2250.83 562.71 4 70.629 3.9462 0.006007 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 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: 854.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9881 -0.5056 0.1133 0.5103 1.7286
##
## Random effects:
## Groups Name Variance Std.Dev.
## AntTag (Intercept) 231.7 15.22
## Residual 155.9 12.49
## Number of obs: 108, groups: AntTag, 39
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 74.879 6.176 57.525 12.125 < 2e-16 ***
## ExpTemp17.5 3.012 5.206 65.451 0.579 0.56490
## ExpTemp21.5 -5.868 5.318 67.967 -1.103 0.27374
## AccTemp17.5 4.792 7.852 63.268 0.610 0.54383
## AccTemp21.5 -12.350 7.921 61.391 -1.559 0.12413
## SpeciesSagmariasus verreauxi 13.231 5.471 34.398 2.418 0.02104 *
## ExpTemp17.5:AccTemp17.5 -4.976 7.362 65.448 -0.676 0.50154
## ExpTemp21.5:AccTemp17.5 8.471 7.442 66.726 1.138 0.25907
## ExpTemp17.5:AccTemp21.5 9.941 7.286 64.266 1.364 0.17724
## ExpTemp21.5:AccTemp21.5 25.207 7.367 65.586 3.422 0.00108 **
## ---
## 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 ET17.5:AT1 ET21.5:AT1
## ExpTemp17.5 -0.422
## ExpTemp21.5 -0.437 0.511
## AccTemp17.5 -0.621 0.329 0.338
## AccTemp21.5 -0.625 0.327 0.336 0.489
## SpcsSgmrssv -0.449 0.007 0.014 -0.016 0.005
## ET17.5:AT17 0.296 -0.707 -0.361 -0.466 -0.231 0.000
## ET21.5:AT17 0.310 -0.365 -0.715 -0.472 -0.240 -0.005 0.516
## ET17.5:AT21 0.302 -0.715 -0.365 -0.235 -0.459 -0.005 0.505 0.261
## ET21.5:AT21 0.315 -0.369 -0.722 -0.244 -0.465 -0.010 0.261 0.516
## ET17.5:AT2
## ExpTemp17.5
## ExpTemp21.5
## AccTemp17.5
## AccTemp21.5
## SpcsSgmrssv
## ET17.5:AT17
## ET21.5:AT17
## ET17.5:AT21
## ET21.5:AT21 0.506
# 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 585.13 292.57 2 64.671 1.8764 0.16138
## AccTemp 182.79 91.40 2 34.352 0.5862 0.56192
## Species 911.77 911.77 1 34.398 5.8477 0.02104 *
## ExpTemp:AccTemp 2256.78 564.20 4 64.639 3.6185 0.01008 *
## ---
## 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 not significant and small (F(2) = 1.88, p = 0.161; Eta2 (partial) = 0.05, 90% CI [0.00, 0.15])
## - The main effect of AccTemp is statistically not significant and small (F(2) = 0.59, p = 0.562; Eta2 (partial) = 0.03, 90% CI [0.00, 0.14])
## - The main effect of Species is statistically significant and large (F(1) = 5.85, p = 0.021; Eta2 (partial) = 0.15, 90% CI [0.01, 0.33])
## - The interaction between ExpTemp and AccTemp is statistically significant and large (F(4) = 3.62, p = 0.010; Eta2 (partial) = 0.18, 90% CI [0.03, 0.29])
##
## Effect sizes were labelled following Field's (2013) recommendations.
# Table output from sjPlot
tab_model(final_model)
value | |||
---|---|---|---|
Predictors | Estimates | CI | p |
(Intercept) | 74.88 | 62.78 – 86.98 | <0.001 |
ExpTemp [17.5] | 3.01 | -7.19 – 13.22 | 0.563 |
ExpTemp [21.5] | -5.87 | -16.29 – 4.56 | 0.270 |
AccTemp [17.5] | 4.79 | -10.60 – 20.18 | 0.542 |
AccTemp [21.5] | -12.35 | -27.88 – 3.18 | 0.119 |
Species [Sagmariasus verreauxi] |
13.23 | 2.51 – 23.95 | 0.016 |
ExpTemp [17.5] * AccTemp [17.5] |
-4.98 | -19.41 – 9.45 | 0.499 |
ExpTemp [21.5] * AccTemp [17.5] |
8.47 | -6.11 – 23.06 | 0.255 |
ExpTemp [17.5] * AccTemp [21.5] |
9.94 | -4.34 – 24.22 | 0.172 |
ExpTemp [21.5] * AccTemp [21.5] |
25.21 | 10.77 – 39.65 | 0.001 |
Random Effects | |||
σ2 | 155.92 | ||
τ00 AntTag | 231.67 | ||
ICC | 0.60 | ||
N AntTag | 39 | ||
Observations | 108 | ||
Marginal R2 / Conditional R2 | 0.172 / 0.667 |
# 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.67) and the part related to the fixed effects alone (marginal R2) is of 0.17. The model's intercept, corresponding to ExpTemp = 14, AccTemp = 14 and Species = Jasus edwardsii, is at 74.88 (95% CI [62.78, 86.98], t(96) = 12.13, p < .001). Within this model:
##
## - The effect of ExpTemp [17.5] is statistically non-significant and positive (beta = 3.01, 95% CI [-7.19, 13.22], t(96) = 0.58, p = 0.563; Std. beta = 0.15, 95% CI [-0.35, 0.64])
## - The effect of ExpTemp [21.5] is statistically non-significant and negative (beta = -5.87, 95% CI [-16.29, 4.56], t(96) = -1.10, p = 0.270; Std. beta = -0.28, 95% CI [-0.78, 0.22])
## - The effect of AccTemp [17.5] is statistically non-significant and positive (beta = 4.79, 95% CI [-10.60, 20.18], t(96) = 0.61, p = 0.542; Std. beta = 0.23, 95% CI [-0.51, 0.97])
## - The effect of AccTemp [21.5] is statistically non-significant and negative (beta = -12.35, 95% CI [-27.88, 3.18], t(96) = -1.56, p = 0.119; Std. beta = -0.59, 95% CI [-1.34, 0.15])
## - The effect of Species [Sagmariasus verreauxi] is statistically significant and positive (beta = 13.23, 95% CI [2.51, 23.95], t(96) = 2.42, p = 0.016; Std. beta = 0.64, 95% CI [0.12, 1.15])
## - The interaction effect of AccTemp [17.5] on ExpTemp [17.5] is statistically non-significant and negative (beta = -4.98, 95% CI [-19.41, 9.45], t(96) = -0.68, p = 0.499; Std. beta = -0.24, 95% CI [-0.93, 0.46])
## - The interaction effect of AccTemp [17.5] on ExpTemp [21.5] is statistically non-significant and positive (beta = 8.47, 95% CI [-6.11, 23.06], t(96) = 1.14, p = 0.255; Std. beta = 0.41, 95% CI [-0.29, 1.11])
## - The interaction effect of AccTemp [21.5] on ExpTemp [17.5] is statistically non-significant and positive (beta = 9.94, 95% CI [-4.34, 24.22], t(96) = 1.36, p = 0.172; Std. beta = 0.48, 95% CI [-0.21, 1.17])
## - The interaction effect of AccTemp [21.5] on ExpTemp [21.5] is statistically significant and positive (beta = 25.21, 95% CI [10.77, 39.65], t(96) = 3.42, p < .001; Std. beta = 1.21, 95% CI [0.52, 1.91])
##
## 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 | 585.1 | 292.6 | 64.7 | 1.88 (2) | 0.16138 | 0.05 [0-0.15] |
AccTemp | 182.8 | 91.4 | 34.4 | 0.59 (2) | 0.56192 | 0.03 [0-0.14] |
Species | 911.8 | 911.8 | 34.4 | 5.85 (1) | 0.02104 | 0.15 [0.01-0.33] |
ExpTemp:AccTemp | 2256.8 | 564.2 | 64.6 | 3.62 (4) | 0.01008 | 0.18 [0.03-0.29] |
# 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)
## –> QQ plot not perfect but acceptable ## –> Model diagnostics do not show apparent violations of model assumptions
# 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("AccTemp", "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_AcclimationTemperature_", trait, ".csv", sep = ""))
# Display as HTML table colouring significant p values red
contrasts_expTemp %>% mutate(p = cell_spec(p, "html", color = ifelse(p < 0.05, "red", "grey"))) %>%
kable(escape = FALSE) %>%
kable_styling(bootstrap_options = c('hover'))
AccTemp | ExpTemp | Species | Difference | t-ratio | p | contrasts_sign | Trait |
---|---|---|---|---|---|---|---|
14 - 17.5 | 14 | Jasus edwardsii | -4.7923955 | -0.6099558 | 1 | No | AerobicScope |
14 - 21.5 | 14 | Jasus edwardsii | 12.3498986 | 1.5583451 | 0.373 | No | AerobicScope |
17.5 - 21.5 | 14 | Jasus edwardsii | 17.1422940 | 2.1505725 | 0.107 | No | AerobicScope |
14 - 17.5 | 17.5 | Jasus edwardsii | 0.1832713 | 0.0232475 | 1 | No | AerobicScope |
14 - 21.5 | 17.5 | Jasus edwardsii | 2.4092909 | 0.3035954 | 1 | No | AerobicScope |
17.5 - 21.5 | 17.5 | Jasus edwardsii | 2.2260196 | 0.2787547 | 1 | No | AerobicScope |
14 - 17.5 | 21.5 | Jasus edwardsii | -13.2633251 | -1.6846128 | 0.291 | No | AerobicScope |
14 - 21.5 | 21.5 | Jasus edwardsii | -12.8573054 | -1.6223711 | 0.33 | No | AerobicScope |
17.5 - 21.5 | 21.5 | Jasus edwardsii | 0.4060196 | 0.0508441 | 1 | No | AerobicScope |
14 - 17.5 | 14 | Sagmariasus verreauxi | -4.7923955 | -0.6099558 | 1 | No | AerobicScope |
14 - 21.5 | 14 | Sagmariasus verreauxi | 12.3498986 | 1.5583451 | 0.373 | No | AerobicScope |
17.5 - 21.5 | 14 | Sagmariasus verreauxi | 17.1422940 | 2.1505725 | 0.107 | No | AerobicScope |
14 - 17.5 | 17.5 | Sagmariasus verreauxi | 0.1832713 | 0.0232475 | 1 | No | AerobicScope |
14 - 21.5 | 17.5 | Sagmariasus verreauxi | 2.4092909 | 0.3035954 | 1 | No | AerobicScope |
17.5 - 21.5 | 17.5 | Sagmariasus verreauxi | 2.2260196 | 0.2787547 | 1 | No | AerobicScope |
14 - 17.5 | 21.5 | Sagmariasus verreauxi | -13.2633251 | -1.6846128 | 0.291 | No | AerobicScope |
14 - 21.5 | 21.5 | Sagmariasus verreauxi | -12.8573054 | -1.6223711 | 0.33 | No | AerobicScope |
17.5 - 21.5 | 21.5 | Sagmariasus verreauxi | 0.4060196 | 0.0508441 | 1 | No | AerobicScope |
# 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("ExpTemp", "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'))
ExpTemp | AccTemp | Species | Difference | t-ratio | p | contrasts_sign | Trait |
---|---|---|---|---|---|---|---|
14 - 17.5 | 14 | Jasus edwardsii | -3.011892 | -0.5777632 | 1 | No | AerobicScope |
14 - 21.5 | 14 | Jasus edwardsii | 5.868037 | 1.1005499 | 0.825 | No | AerobicScope |
17.5 - 21.5 | 14 | Jasus edwardsii | 8.879930 | 1.7034130 | 0.28 | No | AerobicScope |
14 - 17.5 | 17.5 | Jasus edwardsii | 1.963774 | 0.3767108 | 1 | No | AerobicScope |
14 - 21.5 | 17.5 | Jasus edwardsii | -2.602892 | -0.4993128 | 1 | No | AerobicScope |
17.5 - 21.5 | 17.5 | Jasus edwardsii | -4.566667 | -0.8958322 | 1 | No | AerobicScope |
14 - 17.5 | 21.5 | Jasus edwardsii | -12.952500 | -2.5408612 | 0.041 | Yes | AerobicScope |
14 - 21.5 | 21.5 | Jasus edwardsii | -19.339167 | -3.7937184 | 0.001 | Yes | AerobicScope |
17.5 - 21.5 | 21.5 | Jasus edwardsii | -6.386667 | -1.2528572 | 0.645 | No | AerobicScope |
14 - 17.5 | 14 | Sagmariasus verreauxi | -3.011892 | -0.5777632 | 1 | No | AerobicScope |
14 - 21.5 | 14 | Sagmariasus verreauxi | 5.868037 | 1.1005499 | 0.825 | No | AerobicScope |
17.5 - 21.5 | 14 | Sagmariasus verreauxi | 8.879930 | 1.7034130 | 0.28 | No | AerobicScope |
14 - 17.5 | 17.5 | Sagmariasus verreauxi | 1.963774 | 0.3767108 | 1 | No | AerobicScope |
14 - 21.5 | 17.5 | Sagmariasus verreauxi | -2.602892 | -0.4993128 | 1 | No | AerobicScope |
17.5 - 21.5 | 17.5 | Sagmariasus verreauxi | -4.566667 | -0.8958322 | 1 | No | AerobicScope |
14 - 17.5 | 21.5 | Sagmariasus verreauxi | -12.952500 | -2.5408612 | 0.041 | Yes | AerobicScope |
14 - 21.5 | 21.5 | Sagmariasus verreauxi | -19.339167 | -3.7937184 | 0.001 | Yes | AerobicScope |
17.5 - 21.5 | 21.5 | Sagmariasus verreauxi | -6.386667 | -1.2528572 | 0.645 | No | AerobicScope |
# 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")
# Display interaction plots
emmip(final_model, 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: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 | -13.23081 | -2.417294 | 0.021 | Yes | AerobicScope |
Jasus edwardsii - Sagmariasus verreauxi | 17.5 | -13.23081 | -2.417294 | 0.021 | Yes | AerobicScope |
Jasus edwardsii - Sagmariasus verreauxi | 21.5 | -13.23081 | -2.417294 | 0.021 | Yes | AerobicScope |
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)