1. Load libraries
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)
2. Load data
# Select trait
trait <- "escape_total"
# 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)
3. Process data
# 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)
4. Find best probability distribution of data to check assumptions for linear model
If not normal distributed find best distribution fit using fitdistrplus package
# Plot data versus a range of distrubutions
descdist(sub_trait$value, boot = 1000)

## summary statistics
## ------
## min: 9 max: 169
## median: 79
## mean: 82.50962
## estimated sd: 38.24423
## estimated skewness: 0.1752918
## estimated kurtosis: 2.322914
# Fit various distributions
nd <- fitdist(sub_trait$value, "norm")
## $start.arg
## $start.arg$mean
## [1] 82.50962
##
## $start.arg$sd
## [1] 38.05992
##
##
## $fix.arg
## NULL
gd <- fitdist(sub_trait$value, "gamma")
## $start.arg
## $start.arg$shape
## [1] 4.699735
##
## $start.arg$rate
## [1] 0.05695984
##
##
## $fix.arg
## NULL
wd <- fitdist(sub_trait$value, "weibull")
## $start.arg
## $start.arg$shape
## [1] 2.07303
##
## $start.arg$scale
## [1] 94.73581
##
##
## $fix.arg
## NULL
lnd <- fitdist(sub_trait$value, "lnorm")
## $start.arg
## $start.arg$meanlog
## [1] 4.275167
##
## $start.arg$sdlog
## [1] 0.5760732
##
##
## $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.07070856 0.09984929 0.08395626 0.1031488
## Cramer-von Mises statistic 0.07139295 0.16226300 0.06405192 0.3380288
## Anderson-Darling statistic 0.46821302 1.00919041 0.41801271 2.0888024
##
## Goodness-of-fit criteria
## Normal gamma Weibull LogNormal
## Akaike's Information Criterion 1056.085 1058.960 1051.292 1073.658
## Bayesian Information Criterion 1061.374 1064.249 1056.581 1078.947
–> Anderson-Darling statistics gives good measure for fit for both middle and tail of distribution, the lower the better
–> Data fit best to a normal distribution
5. Plot data
# 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_total_Acclimation.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)
# 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_total_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_total_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_total_Species.pdf", width = 14, height = 8)
6. Summary statistics
By species
# 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()
escape_total
Species
|
mean
|
n
|
sd
|
ci
|
lower_ci
|
upper_ci
|
min
|
max
|
Jasus edwardsii
|
68.0
|
53
|
34.3
|
9.5
|
58.5
|
77.4
|
9
|
148
|
Sagmariasus verreauxi
|
97.6
|
51
|
36.5
|
10.2
|
87.4
|
107.9
|
27
|
169
|
Grouped by species, experimental temperature, and acclimation
# Summary statistics
substrait_stats <- ddply(sub_trait, c("Species", "ExpTemp", "AccTemp"), summarise,
mean = mean(value, na.rm = TRUE),
n = length(value),
sd = sd(value, na.rm = TRUE),
ci = qt(0.975, n) * sd / sqrt(n),
lower_ci = mean-ci,
upper_ci = mean+ci,
min = min(value),
max = max(value))
# Round values
substrait_stats <- substrait_stats %>% mutate_if(is.numeric, round, digit = 3)
# Display results
substrait_stats %>% kbl(digits = 1, caption = trait) %>% kable_classic()
escape_total
Species
|
ExpTemp
|
AccTemp
|
mean
|
n
|
sd
|
ci
|
lower_ci
|
upper_ci
|
min
|
max
|
Jasus edwardsii
|
14
|
14
|
69.3
|
6
|
30.4
|
30.4
|
39.0
|
99.7
|
24
|
109
|
Jasus edwardsii
|
14
|
17.5
|
79.5
|
6
|
30.8
|
30.8
|
48.7
|
110.3
|
55
|
131
|
Jasus edwardsii
|
14
|
21.5
|
37.8
|
6
|
22.9
|
22.8
|
15.0
|
60.7
|
19
|
70
|
Jasus edwardsii
|
17.5
|
14
|
85.5
|
6
|
31.5
|
31.5
|
54.0
|
117.0
|
46
|
124
|
Jasus edwardsii
|
17.5
|
17.5
|
64.7
|
6
|
38.7
|
38.7
|
26.0
|
103.4
|
9
|
110
|
Jasus edwardsii
|
17.5
|
21.5
|
69.8
|
5
|
47.3
|
54.3
|
15.5
|
124.1
|
30
|
148
|
Jasus edwardsii
|
21.5
|
14
|
62.3
|
6
|
42.6
|
42.6
|
19.8
|
104.9
|
15
|
117
|
Jasus edwardsii
|
21.5
|
17.5
|
87.5
|
6
|
27.7
|
27.7
|
59.8
|
115.2
|
59
|
134
|
Jasus edwardsii
|
21.5
|
21.5
|
55.5
|
6
|
26.1
|
26.0
|
29.5
|
81.5
|
30
|
102
|
Sagmariasus verreauxi
|
14
|
14
|
114.0
|
6
|
16.9
|
16.9
|
97.1
|
130.9
|
89
|
135
|
Sagmariasus verreauxi
|
14
|
17.5
|
86.7
|
6
|
40.3
|
40.3
|
46.4
|
127.0
|
27
|
145
|
Sagmariasus verreauxi
|
14
|
21.5
|
84.8
|
5
|
37.8
|
43.4
|
41.4
|
128.2
|
43
|
127
|
Sagmariasus verreauxi
|
17.5
|
14
|
125.8
|
5
|
38.0
|
43.7
|
82.1
|
169.5
|
73
|
168
|
Sagmariasus verreauxi
|
17.5
|
17.5
|
119.8
|
6
|
39.7
|
39.6
|
80.2
|
159.5
|
80
|
169
|
Sagmariasus verreauxi
|
17.5
|
21.5
|
86.8
|
6
|
35.4
|
35.4
|
51.4
|
122.2
|
57
|
150
|
Sagmariasus verreauxi
|
21.5
|
14
|
82.6
|
5
|
43.4
|
49.9
|
32.7
|
132.5
|
31
|
130
|
Sagmariasus verreauxi
|
21.5
|
17.5
|
95.7
|
6
|
37.8
|
37.8
|
57.9
|
133.5
|
29
|
143
|
Sagmariasus verreauxi
|
21.5
|
21.5
|
82.5
|
6
|
25.5
|
25.4
|
57.1
|
107.9
|
42
|
108
|
Decreased number of total escapes by trend in Jasus edwardsii between 17.5 and 21.5 acclimation temperature at 14 and 21.5 experimental temperature
# At 14°C Experimental temperature
print(paste("decrease of", round((79.5-37.8)/79.5*100, 0), "% from 17.5 to 21.5°C acclimation temperature at 17.5°C experimental temperature"))
## [1] "decrease of 52 % from 17.5 to 21.5°C acclimation temperature at 17.5°C experimental temperature"
# At 21.5°C Experimental temperature
print(paste("decrease of", round((87.5-55.5)/87.5*100, 0), "% from 17.5 to 21.5°C acclimation temperature at 21.5°C experimental temperature"))
## [1] "decrease of 37 % from 17.5 to 21.5°C acclimation temperature at 21.5°C experimental temperature"
6. Correlations
6.B EPOC correlations
# 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 = "Total escapes", ylab = "EPOC")
## `geom_smooth()` using formula 'y ~ x'

ggsave("Total escapes_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.98224, p-value = 0.6372
shapiro.test(SRL$value)
##
## Shapiro-Wilk normality test
##
## data: SRL$value
## W = 0.96634, p-value = 0.1397
# Correlation plots
ggscatter(ERL, x = "value", y = "EPOC",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Total escapes", 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 = "Total escapes", 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 = 1.6862, df = 49, p-value = 0.09812
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.04425861 0.47886099
## sample estimates:
## cor
## 0.2341814
cor.test(SRL$value, SRL$EPOC, method = "pearson")
##
## Pearson's product-moment correlation
##
## data: SRL$value and SRL$EPOC
## t = 3.0299, df = 51, p-value = 0.003835
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1344801 0.5977643
## sample estimates:
## cor
## 0.3905768
6.B MMR correlations
# 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 = "Total escapes", 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.98224, p-value = 0.6372
shapiro.test(SRL$value)
##
## Shapiro-Wilk normality test
##
## data: SRL$value
## W = 0.96634, p-value = 0.1397
# Correlation plots
ggscatter(ERL, x = "value", y = "MMR",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Total escapes", 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 = "Total escapes", 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.2604, df = 49, p-value = 0.2135
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1034275 0.4317081
## sample estimates:
## cor
## 0.1772071
cor.test(SRL$value, SRL$MMR, method = "pearson")
##
## Pearson's product-moment correlation
##
## data: SRL$value and SRL$MMR
## t = -0.57954, df = 51, p-value = 0.5648
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3436663 0.1936413
## sample estimates:
## cor
## -0.08088605
6.C Weight correlations
# Pooled data
ggscatter(sub_trait, x = "value", y = "Weight_g", fill = "Species",color = "Species",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Total escapes", ylab = "Weight_g")
## `geom_smooth()` using formula 'y ~ x'

ggsave("Escapes_speed_Corelation_Weight_g.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.98224, p-value = 0.6372
shapiro.test(SRL$value)
##
## Shapiro-Wilk normality test
##
## data: SRL$value
## W = 0.96634, p-value = 0.1397
# Correlation plots
ggscatter(ERL, x = "value", y = "Weight_g",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Total escapes", ylab = "Weight_g")+
ggtitle("Sagmariasus verreauxi")
## `geom_smooth()` using formula 'y ~ x'

ggscatter(SRL, x = "value", y = "Weight_g",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Total escapes", ylab = "Weight_g")+
ggtitle("Jasus edwardsii")
## `geom_smooth()` using formula 'y ~ x'

# Correlation tests
cor.test(ERL$value, ERL$Weight_g, method = "pearson")
##
## Pearson's product-moment correlation
##
## data: ERL$value and ERL$Weight_g
## t = 0.56767, df = 49, p-value = 0.5729
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1991903 0.3486478
## sample estimates:
## cor
## 0.08083047
cor.test(SRL$value, SRL$Weight_g, method = "pearson")
##
## Pearson's product-moment correlation
##
## data: SRL$value and SRL$Weight_g
## t = 4.1276, df = 51, p-value = 0.0001358
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2661080 0.6788791
## sample estimates:
## cor
## 0.5004087
7. Linear mixed effect modelling
7.A Find best model
Fit the most complex model with all meaningful fixed factors
Animal ID (AntTag) as random factor
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
anova(complex_model)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ExpTemp 3632.5 1816.2 2 65.831 4.7614 0.011714 *
## AccTemp 1553.9 776.9 2 35.904 2.0368 0.145247
## Species 3938.0 3938.0 1 35.939 10.3237 0.002771 **
## Sex 1201.9 1201.9 1 38.184 3.1510 0.083853 .
## Weight_g 1338.4 1338.4 1 40.472 3.5086 0.068285 .
## ExpTemp:AccTemp 3751.4 937.8 4 65.684 2.4586 0.053996 .
## ExpTemp:Species 2405.6 1202.8 2 66.285 3.1532 0.049188 *
## AccTemp:Species 370.4 185.2 2 36.117 0.4855 0.619372
## ExpTemp:AccTemp:Species 3640.5 910.1 4 66.052 2.3860 0.059942 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
7.B Find the best model by backward elimination of non-significant effects of linear mixed effects model using the step function of the LmerTest package
(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 -484.93 1013.9
## (1 | AntTag) 0 21 -497.50 1037.0 25.155 1 5.289e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Sex 1 1201.9 1201.9 1 38.184 3.1510 0.083853
## ExpTemp:AccTemp:Species 2 3660.8 915.2 4 66.193 2.4133 0.057586
## AccTemp:Species 3 405.9 203.0 2 36.232 0.4777 0.624045
## ExpTemp:AccTemp 4 3850.0 962.5 4 65.359 2.2604 0.072038
## AccTemp 5 1694.7 847.3 2 36.143 1.7383 0.190196
## ExpTemp:Species 6 2442.9 1221.5 2 66.565 2.5203 0.088088
## ExpTemp 7 3081.4 1540.7 2 66.395 2.9190 0.060951
## Species 0 5076.7 5076.7 1 36.470 8.9186 0.005021
## Weight_g 0 2356.4 2356.4 1 39.741 4.1396 0.048601
##
## Sex .
## ExpTemp:AccTemp:Species .
## AccTemp:Species
## ExpTemp:AccTemp .
## AccTemp
## ExpTemp:Species .
## ExpTemp .
## Species **
## Weight_g *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## value ~ Species + Weight_g + (1 | AntTag)
# Extract the model that step found:
final_model <- get_model(step_result)
–> effect of species and weak effect of weight
–> no acclimation effects
7.C Re-run model fit using REML estimation
Best model: Species + Weight_g + (1 | AntTag)
# Remodel using Maximum likelyhood estimation
final_model <- lmer(value ~ Species +
Weight_g +
(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 ~ Species + Weight_g + (1 | AntTag)
## Data: sub_trait
##
## REML criterion at convergence: 999.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.93924 -0.57972 0.05684 0.56562 2.17566
##
## Random effects:
## Groups Name Variance Std.Dev.
## AntTag (Intercept) 647.3 25.44
## Residual 568.2 23.84
## Number of obs: 104, groups: AntTag, 39
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 16.34591 26.83250 37.26126 0.609 0.5461
## SpeciesSagmariasus verreauxi 27.59436 9.61042 33.89699 2.871 0.0070 **
## Weight_g 0.04677 0.02425 37.45211 1.929 0.0614 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) SpcsSv
## SpcsSgmrssv -0.016
## Weight_g -0.968 -0.164
# 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)
## Species 4684.5 4684.5 1 33.897 8.2443 0.006998 **
## Weight_g 2113.4 2113.4 1 37.452 3.7194 0.061390 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
report(anova(final_model))
## The ANOVA suggests that:
##
## - The main effect of Species is statistically significant and large (F(1) = 8.24, p = 0.007; Eta2 (partial) = 0.20, 90% CI [0.04, 0.38])
## - The main effect of Weight_g is statistically not significant and medium (F(1) = 3.72, p = 0.061; Eta2 (partial) = 0.09, 90% CI [0.00, 0.26])
##
## Effect sizes were labelled following Field's (2013) recommendations.
# Table output from sjPlot
tab_model(final_model)
|
value
|
Predictors
|
Estimates
|
CI
|
p
|
(Intercept)
|
16.35
|
-36.24 – 68.94
|
0.542
|
Species [Sagmariasus verreauxi]
|
27.59
|
8.76 – 46.43
|
0.004
|
Weight_g
|
0.05
|
-0.00 – 0.09
|
0.054
|
Random Effects
|
σ2
|
568.21
|
τ00 AntTag
|
647.28
|
ICC
|
0.53
|
N AntTag
|
39
|
Observations
|
104
|
Marginal R2 / Conditional R2
|
0.208 / 0.630
|
# Textual report of the model
report(final_model)
## We fitted a linear mixed model (estimated using REML and nloptwrap optimizer) to predict value with Species and Weight_g (formula: value ~ Species + Weight_g). The model included AntTag as random effect (formula: ~1 | AntTag). The model's total explanatory power is substantial (conditional R2 = 0.63) and the part related to the fixed effects alone (marginal R2) is of 0.21. The model's intercept, corresponding to Species = Jasus edwardsii and Weight_g = 0, is at 16.35 (95% CI [-36.24, 68.94], t(99) = 0.61, p = 0.542). Within this model:
##
## - The effect of Species [Sagmariasus verreauxi] is statistically significant and positive (beta = 27.59, 95% CI [8.76, 46.43], t(99) = 2.87, p = 0.004; Std. beta = 0.72, 95% CI [0.23, 1.21])
## - The effect of Weight_g is statistically non-significant and positive (beta = 0.05, 95% CI [-7.61e-04, 0.09], t(99) = 1.93, p = 0.054; Std. beta = 0.24, 95% CI [-3.94e-03, 0.49])
##
## 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.
Summary table of main effects
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")
Model summary for main effects
|
Sum of sqares
|
Mean squares
|
Den df
|
F(df)
|
p value
|
Partical squared eta [CI 90%]
|
Species
|
4684.5
|
4684.5
|
33.9
|
8.24 (1)
|
0.006998
|
0.2 [0.04-0.38]
|
Weight_g
|
2113.4
|
2113.4
|
37.5
|
3.72 (1)
|
0.061390
|
0.09 [0-0.26]
|
8. Model diagnostics
# 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$Species)

plotResiduals(simulationOutput, sub_trait$Weight_g)

–> Model diagnostics do not show apparent violations of model assumptions