# rm(list=ls())
# gc()
library(dplyr)
library(tidyverse)
library(tidyr)
library(ggplot2)
library(knitr)
library(pbapply)
library(data.table)
library(GGally)
library(here)
library(terra)
library(tidyterra)
library(rcompanion)
library(lme4)
library(glmmTMB)
library(MuMIn)
library(emmeans)
library(parallel)
library(sjPlot)
library(jtools)
library(deming)
library(cowplot)
library(DescTools)
library(gtools)
library(ggpp)
computer = "personal"
source(here("R/settings.R"))
source(here("R/my_functions.R"))
source(here("R/velocity_functions.R"))
source(here("R/emm_extra_functions.R"))
options(na.action = "na.fail")Exploratory analyses
Data amount: all
1 Setup
2 Load data
# data amount
if(params$data_amount=="all"){
bioshifts <- read.csv(here("Data/Output/Bioshifts_merge_Exposure_all.csv"))
} else {
bioshifts <- read.csv(here("Data/Output/Bioshifts_merge_Exposure.csv"))
}# select classes with more than X observations per edge
N_obs_class = 5
class_param_select <- bioshifts %>%
group_by(class, Param) %>%
tally() %>%
filter(n >= N_obs_class) %>%
mutate(class_param = paste(class,Param))
bioshifts <- bioshifts %>%
mutate(class_param = paste(class,Param)) %>%
filter(class_param %in% class_param_select$class_param) %>%
select(!class_param)# N species
length(unique(bioshifts$sp_name_std))[1] 3623bioshifts %>%
group_by(Param) %>%
summarise(N=length(unique(sp_name_std)))# A tibble: 3 × 2
Param N
<chr> <int>
1 LE 1772
2 O 2647
3 TE 568bioshifts %>%
group_by(Eco) %>%
summarise(N=length(unique(sp_name_std)))# A tibble: 2 × 2
Eco N
<chr> <int>
1 Mar 519
2 Ter 3104bioshifts %>%
group_by(Eco,Param) %>%
summarise(N=length(unique(sp_name_std)))`summarise()` has grouped output by 'Eco'. You can override using the `.groups`
argument.# A tibble: 6 × 3
# Groups: Eco [2]
Eco Param N
<chr> <chr> <int>
1 Mar LE 206
2 Mar O 388
3 Mar TE 75
4 Ter LE 1566
5 Ter O 2259
6 Ter TE 493# N range shifts
length(bioshifts$sp_name_std)[1] 9579bioshifts %>%
group_by(Param) %>%
summarise(N=length(sp_name_std))# A tibble: 3 × 2
Param N
<chr> <int>
1 LE 2468
2 O 6351
3 TE 760bioshifts %>%
group_by(Eco) %>%
summarise(N=length(sp_name_std))# A tibble: 2 × 2
Eco N
<chr> <int>
1 Mar 937
2 Ter 8642bioshifts %>%
group_by(Eco,Param) %>%
summarise(N=length(sp_name_std))`summarise()` has grouped output by 'Eco'. You can override using the `.groups`
argument.# A tibble: 6 × 3
# Groups: Eco [2]
Eco Param N
<chr> <chr> <int>
1 Mar LE 211
2 Mar O 650
3 Mar TE 76
4 Ter LE 2257
5 Ter O 5701
6 Ter TE 6843 Make dataset for modeling
mod_data <- bioshifts %>%
mutate(ShiftRate_raw = ShiftRate,
ShiftMod_raw = bvel_lat,
ShiftMod_raw95 = bvel95_lat,
ShiftRate = scale(ShiftRate),
ShiftMod = scale(bvel_lat),
ShiftMod95 = scale(bvel95_lat),
Lat_extent = scale(log(LatExtentk)),
Study_duration = scale(Duration),
TSS = scale(calibration),
Interannual_var = scale(log1p(bvel_lat_error)),
N_periods = scale(N_periodes),
Connectivity = scale((connectivity)),
match_direction = relevel(factor(match_direction), ref = "mismatch"),
mismatch = abs(mismatch),
Mobility = as.character(str_to_sentence(Mobility)),
Locomotion_mode = as.character(str_to_sentence(Locomotion_mode)),
Life_form = as.character(str_to_sentence(lifeform)),
Match_dir = case_when(match_direction_class=="pospos" ~ "Poleward",
match_direction_class=="negneg" ~ "Equatorward",
match_direction == "mismatch" ~ "Mismatch"),
mismatch_binary = case_when(match_direction=="mismatch" ~ 0,
match_direction=="match" ~ 1),
Sampling_category = Category,
Random_ID_Param_Class = paste(ID,Param,class),
Random_ID_Class = paste(ID,class),
Random_Param_Class = paste(Param,class)) %>%
select(ShiftRate_raw, ShiftMod_raw, ShiftMod_raw95, ShiftRate, ShiftMod, ShiftMod95,
mismatch, match_direction, mismatch_binary, Match_dir, mismatch_trend,
Lat_extent, Sampling_category, TSS, Interannual_var, N_periods, Study_duration,
Life_form, Mobility, Locomotion_mode, Connectivity,
ID, Param, class, Eco, sp_name_std,
Random_ID_Param_Class, Random_Param_Class, Random_ID_Class)Extra fixes for variables
# fixing Locomotion_mode
unique(mod_data$Locomotion_mode) [1] "Flying" "Running/flying" "Running"
[4] "Slithering" "Swimming" "Planktonic"
[7] "Crawling" "Sessile" "Swimming/walking"
[10] "Walking/hopping" "Walking" "Burrowing"
[13] "Swimming/burrowing"pos <- grep("lying",mod_data$Locomotion_mode)
mod_data$Locomotion_mode[pos] <- "Flying"
pos <- grep("wimming",mod_data$Locomotion_mode)
mod_data$Locomotion_mode[pos] <- "Swimming"
pos <- grep("alk",mod_data$Locomotion_mode)
mod_data$Locomotion_mode[pos] <- "Run/Walk"
pos <- grep("unning",mod_data$Locomotion_mode)
mod_data$Locomotion_mode[pos] <- "Run/Walk"
pos <- grep("lithering",mod_data$Locomotion_mode)
mod_data$Locomotion_mode[pos] <- "Slithering/Crawling"
pos <- grep("rawling",mod_data$Locomotion_mode)
mod_data$Locomotion_mode[pos] <- "Slithering/Crawling"
pos <- grep("Burrowing",mod_data$Locomotion_mode)
mod_data$Locomotion_mode[pos] <- "Sessile"
unique(mod_data$Locomotion_mode)[1] "Flying" "Run/Walk" "Slithering/Crawling"
[4] "Swimming" "Planktonic" "Sessile" mod_data$Locomotion_mode <- factor(
mod_data$Locomotion_mode,
levels = c("Flying", "Run/Walk", "Swimming", "Slithering/Crawling", "Planktonic", "Sessile"))
table(mod_data$Locomotion_mode[which(mod_data$Eco=="Ter")])
Flying Run/Walk Swimming Slithering/Crawling
5477 91 12 27
Planktonic Sessile
0 3035 table(mod_data$Locomotion_mode[which(mod_data$Eco=="Mar")])
Flying Run/Walk Swimming Slithering/Crawling
0 0 794 50
Planktonic Sessile
28 65 # fixing Sampling_category
table(mod_data$Sampling_category)
Census-Resurvey CensusPeriods Survey-Resurvey Synthesis TimeSeries
34 6525 1010 4 2006 # remove synthresis because there is a only a single study
mod_data <- mod_data %>% filter(!Sampling_category == "Synthesis")
# combine census-resurvey and survey-resurvey
pos <- grep("Census-Resurvey",mod_data$Sampling_category)
mod_data$Sampling_category[pos] <- "Survey-Resurvey"
table(mod_data$Sampling_category)
CensusPeriods Survey-Resurvey TimeSeries
6525 1044 2006 # Fix direction levels
mod_data$Match_dir <- factor(mod_data$Match_dir, levels = c("Poleward", "Equatorward", "Mismatch"))
table(mod_data$Mobility[which(mod_data$Eco=="Ter")])
Mobile Sessile
5606 3032 table(mod_data$Mobility[which(mod_data$Eco=="Mar")])
Drift Facultatively mobile Mobile
28 44 842
Sessile
23 table(mod_data$Life_form[which(mod_data$Eco=="Ter")])
Cryptogam Ectotherm Endotherm Phanerogam
94 1688 3918 2938 table(mod_data$Life_form[which(mod_data$Eco=="Mar")])
Cryptogam Ectotherm Endotherm
45 884 8 # N classes
mod_data %>%
group_by(Eco) %>%
summarise(N = length(unique(class)))# A tibble: 2 × 2
Eco N
<chr> <int>
1 Mar 15
2 Ter 10# N studies
mod_data %>%
group_by(Eco) %>%
summarise(N = length(unique(ID)))# A tibble: 2 × 2
Eco N
<chr> <int>
1 Mar 44
2 Ter 94# make numeric as.numeric
pos <- which(sapply(mod_data, is.numeric))
tmp <- lapply(mod_data[pos], as.numeric)
tmp <- do.call(cbind,tmp)
mod_data[pos] <- tmp4 Check for correlations
4.1 Scatters
# is Study_duration correlated with range shift mismatch?
ggplot(mod_data, aes(Study_duration,mismatch))+
geom_point()+
theme_classic()+
geom_smooth()+
facet_wrap(.~Eco)`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
# is NPeriods correlated with range shift mismatch?
ggplot(mod_data, aes((N_periods),mismatch))+
geom_point()+
theme_classic()+
geom_smooth()+
facet_wrap(.~Eco)`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
# is env variability correlated with range shift mismatch?
ggplot(mod_data, aes(Interannual_var,mismatch))+
geom_point()+
theme_classic()+
geom_smooth()+
facet_wrap(.~Eco)`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
# is NPeriods correlated with Study_duration?
ggplot(mod_data, aes((N_periods),Study_duration))+
geom_point()+
theme_classic()+
geom_smooth()+
facet_wrap(.~Eco)`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
# is NPeriods correlated with Sampling_category?
ggplot(mod_data, aes(Sampling_category,(N_periods)))+
geom_point()+
theme_classic()+
geom_smooth()+
facet_wrap(.~Eco)`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
# is Study_duration correlated with Sampling_category?
ggplot(mod_data, aes(Sampling_category,Study_duration))+
geom_point()+
theme_classic()+
geom_smooth()+
facet_wrap(.~Eco)`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
# is Study_duration correlated with env variability?
ggplot(mod_data, aes(Interannual_var,Study_duration))+
geom_point()+
theme_classic()+
geom_smooth()+
facet_wrap(.~Eco)`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
Removed N_periods as it is highly correlated with Sampling_category
5 VIF
5.1 Terrestrial
subdata_ter <- mod_data %>%
filter(Eco=="Ter") %>%
select("mismatch", "TSS", "Study_duration", "Interannual_var", "Lat_extent", "Connectivity", "Sampling_category", "Life_form", "Locomotion_mode") %>%
na.omit
mm_formula_ter <- "mismatch ~ TSS + Study_duration + Interannual_var + Lat_extent + Connectivity + Sampling_category + Life_form + Locomotion_mode"
m1 <- glm(mm_formula_ter,
data = subdata_ter,
family = Gamma(link = "log"))
summary(m1)
Call:
glm(formula = mm_formula_ter, family = Gamma(link = "log"), data = subdata_ter)
Coefficients: (1 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.14314 0.13441 1.065 0.286932
TSS -0.04776 0.01584 -3.016 0.002571 **
Study_duration -0.24158 0.02010 -12.017 < 2e-16 ***
Interannual_var 0.22560 0.01960 11.508 < 2e-16 ***
Lat_extent 0.43081 0.02698 15.968 < 2e-16 ***
Connectivity -0.11750 0.01720 -6.830 9.04e-12 ***
Sampling_categorySurvey-Resurvey -0.20582 0.05738 -3.587 0.000336 ***
Sampling_categoryTimeSeries -0.08206 0.04638 -1.769 0.076873 .
Life_formEctotherm 0.80750 0.14125 5.717 1.12e-08 ***
Life_formEndotherm 0.48740 0.13970 3.489 0.000488 ***
Life_formPhanerogam -0.22967 0.13378 -1.717 0.086064 .
Locomotion_modeRun/Walk -0.49722 0.14027 -3.545 0.000395 ***
Locomotion_modeSwimming -0.85041 0.37103 -2.292 0.021928 *
Locomotion_modeSlithering/Crawling -0.07291 0.25403 -0.287 0.774096
Locomotion_modeSessile NA NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 1.618215)
Null deviance: 18518 on 8603 degrees of freedom
Residual deviance: 12358 on 8590 degrees of freedom
AIC: 24224
Number of Fisher Scoring iterations: 8# locomotion mode is problematic >> remove
mm_formula_ter <- "mismatch ~ TSS + Study_duration + Interannual_var + Lat_extent + Connectivity + Sampling_category + Life_form"
m1 <- glm(mm_formula_ter,
data = subdata_ter,
family = Gamma(link = "log"))
summary(m1)
Call:
glm(formula = mm_formula_ter, family = Gamma(link = "log"), data = subdata_ter)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.14743 0.13385 1.102 0.270709
TSS -0.04402 0.01575 -2.794 0.005220 **
Study_duration -0.25226 0.01990 -12.675 < 2e-16 ***
Interannual_var 0.22450 0.01944 11.547 < 2e-16 ***
Lat_extent 0.44050 0.02675 16.468 < 2e-16 ***
Connectivity -0.12090 0.01706 -7.085 1.50e-12 ***
Sampling_categorySurvey-Resurvey -0.18436 0.05606 -3.288 0.001011 **
Sampling_categoryTimeSeries -0.07412 0.04598 -1.612 0.107022
Life_formEctotherm 0.77679 0.14003 5.547 2.99e-08 ***
Life_formEndotherm 0.46906 0.13908 3.373 0.000748 ***
Life_formPhanerogam -0.22572 0.13323 -1.694 0.090250 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 1.604798)
Null deviance: 18518 on 8603 degrees of freedom
Residual deviance: 12382 on 8593 degrees of freedom
AIC: 24238
Number of Fisher Scoring iterations: 8VIF(m1) GVIF Df GVIF^(1/(2*Df))
TSS 1.232935 1 1.110376
Study_duration 2.236088 1 1.495356
Interannual_var 2.053056 1 1.432849
Lat_extent 3.912820 1 1.978085
Connectivity 1.560710 1 1.249284
Sampling_category 2.396487 2 1.244210
Life_form 4.416607 3 1.280899# high VIF value for life form
# try keeping locomotion mode and removing life form
mm_formula_ter <- "mismatch ~ TSS + Study_duration + Interannual_var + Lat_extent + Connectivity + Sampling_category + Locomotion_mode"
m1 <- glm(mm_formula_ter,
data = subdata_ter,
family = Gamma(link = "log"))
summary(m1)
Call:
glm(formula = mm_formula_ter, family = Gamma(link = "log"), data = subdata_ter)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.73919 0.02480 29.805 < 2e-16 ***
TSS -0.02997 0.01580 -1.897 0.05780 .
Study_duration -0.19629 0.01793 -10.947 < 2e-16 ***
Interannual_var 0.24505 0.01917 12.781 < 2e-16 ***
Lat_extent 0.37246 0.02484 14.994 < 2e-16 ***
Connectivity -0.10212 0.01705 -5.988 2.2e-09 ***
Sampling_categorySurvey-Resurvey -0.02824 0.04958 -0.570 0.56895
Sampling_categoryTimeSeries -0.10138 0.04638 -2.186 0.02885 *
Locomotion_modeRun/Walk -0.45106 0.13895 -3.246 0.00117 **
Locomotion_modeSwimming -0.70875 0.37192 -1.906 0.05673 .
Locomotion_modeSlithering/Crawling 0.10591 0.25237 0.420 0.67474
Locomotion_modeSessile -0.86913 0.04373 -19.873 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 1.630934)
Null deviance: 18518 on 8603 degrees of freedom
Residual deviance: 12429 on 8592 degrees of freedom
AIC: 24280
Number of Fisher Scoring iterations: 8VIF(m1) # Ok GVIF Df GVIF^(1/(2*Df))
TSS 1.219583 1 1.104347
Study_duration 1.785977 1 1.336405
Interannual_var 1.964511 1 1.401610
Lat_extent 3.320521 1 1.822230
Connectivity 1.533969 1 1.238535
Sampling_category 1.764746 2 1.152579
Locomotion_mode 2.405843 4 1.1159855.2 Marine
subdata_mar <- mod_data %>%
filter(Eco=="Mar") %>%
select("mismatch", "TSS", "Study_duration", "N_periods", "Interannual_var", "Lat_extent", "Sampling_category", "Life_form", "Locomotion_mode") %>%
na.omit
mm_formula_mar <- "mismatch ~ TSS + Study_duration + Interannual_var + Lat_extent + Sampling_category + Life_form + Locomotion_mode"
m1 <- glm(mm_formula_mar,
data = subdata_mar,
family = Gamma(link = "log"))
summary(m1)
Call:
glm(formula = mm_formula_mar, family = Gamma(link = "log"), data = subdata_mar)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.13607 0.34325 -0.396 0.69190
TSS 0.03874 0.03126 1.239 0.21549
Study_duration -0.43185 0.08895 -4.855 1.41e-06 ***
Interannual_var 0.13187 0.04804 2.745 0.00616 **
Lat_extent 0.65795 0.05584 11.782 < 2e-16 ***
Sampling_categorySurvey-Resurvey -0.16426 0.13620 -1.206 0.22812
Sampling_categoryTimeSeries 0.25666 0.08830 2.907 0.00374 **
Life_formEctotherm 1.33154 0.33919 3.926 9.29e-05 ***
Life_formEndotherm 0.84553 0.52782 1.602 0.10952
Locomotion_modeSlithering/Crawling -0.26091 0.17229 -1.514 0.13029
Locomotion_modePlanktonic 0.94827 0.39571 2.396 0.01676 *
Locomotion_modeSessile 0.32037 0.18420 1.739 0.08232 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 1.262171)
Null deviance: 1615.7 on 936 degrees of freedom
Residual deviance: 1278.0 on 925 degrees of freedom
AIC: 4075.7
Number of Fisher Scoring iterations: 11VIF(m1) # Ok! GVIF Df GVIF^(1/(2*Df))
TSS 1.071991 1 1.035370
Study_duration 2.024081 1 1.422702
Interannual_var 1.529302 1 1.236649
Lat_extent 1.309589 1 1.144373
Sampling_category 2.079162 2 1.200804
Life_form 3.965445 2 1.411149
Locomotion_mode 4.799972 3 1.2987936 Models
Set random and fixed effects
# Random effects
## Nested
RE_nest <- "(1 | Param/class/ID)"
## custom
RE_cust <- "(1 | Random_ID_Param_Class)"
########################
# fixed effects
cont_vars <-
c(
"TSS",
"Study_duration",
"Interannual_var",
"Lat_extent",
"Connectivity"
)
cont_vars_names <-
c(
"SDM performance",
"Study duration",
"Climate variability",
"Latitudinal extent",
"Landscape connectivity"
)
cat_vars <- c(
"Life_form",
"Sampling_category",
"Locomotion_mode"
)
cat_names <- c(
"Life form",
"Sampling_category",
"Locomotion mode"
)6.1 Documented velocities in function of modeled velocities
6.1.1 Simple model - species-level (Shift ~ ShiftMod)
shiftmodtype <- ""
# shiftmodtype <- "95"
if(shiftmodtype=="95"){
sm_fixed_eff <- "ShiftRate ~ ShiftMod95"
sm_formula <- paste(sm_fixed_eff, RE_nest, sep = " + ")
} else {
sm_fixed_eff <- "ShiftRate ~ ShiftMod"
sm_formula <- paste(sm_fixed_eff, RE_nest, sep = " + ")
}6.1.1.1 Terrestrial
subdata_ter <- mod_data %>%
filter(Eco=="Ter")
# length(unique(mod_data$Random_ID_Param_Class))
sm_ter <- glmmTMB(as.formula(sm_formula),
data = subdata_ter)
tab_model(sm_ter)| Shift Rate | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | -0.01 | -0.24 – 0.22 | 0.939 |
| ShiftMod | -0.02 | -0.04 – -0.00 | 0.025 |
| Random Effects | |||
| σ2 | 0.74 | ||
| τ00 ID:class:Param | 0.48 | ||
| τ00 class:Param | 0.02 | ||
| τ00 Param | 0.02 | ||
| ICC | 0.41 | ||
| N ID | 94 | ||
| N class | 10 | ||
| N Param | 3 | ||
| Observations | 8638 | ||
| Marginal R2 / Conditional R2 | 0.000 / 0.414 | ||
plot(effects::allEffects(sm_ter))
No association between modeled shift velocity and observed shift velocity.
6.1.1.2 Marine
subdata_mar <- mod_data %>%
filter(Eco=="Mar")
# length(unique(subdata_mar$Random_ID_Param_Class))
sm_mar <- glmmTMB(as.formula(sm_formula),
data = subdata_mar)
tab_model(sm_mar)| Shift Rate | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 0.31 | -0.03 – 0.65 | 0.071 |
| ShiftMod | -0.01 | -0.12 – 0.09 | 0.781 |
| Random Effects | |||
| σ2 | 1.63 | ||
| τ00 ID:class:Param | 0.79 | ||
| τ00 class:Param | 0.03 | ||
| τ00 Param | 0.02 | ||
| ICC | 0.34 | ||
| N ID | 44 | ||
| N class | 15 | ||
| N Param | 3 | ||
| Observations | 937 | ||
| Marginal R2 / Conditional R2 | 0.000 / 0.339 | ||
plot(effects::allEffects(sm_mar))
Same here…
6.1.1.3 Save model results
if(shiftmodtype=="95"){
sm_mar_raw <- glmmTMB(
as.formula(paste("ShiftRate_raw~ShiftMod_raw95", RE_nest, sep = " + ")),
data = subdata_mar)
} else {
sm_mar_raw <- glmmTMB(
as.formula(paste("ShiftRate_raw~ShiftMod_raw", RE_nest, sep = " + ")),
data = subdata_mar)
}
mar_sm_res <- summary(sm_mar_raw)
mar_sm_res <- data.frame(mar_sm_res$coefficients$cond)
mar_sm_res$Param <- rownames(mar_sm_res)
names(mar_sm_res) <- c("Estimate", "Std_Error", "z-value", "P-value", "Param")
if(shiftmodtype=="95"){
sm_ter_raw <- glmmTMB(
as.formula(paste("ShiftRate_raw~ShiftMod_raw95", RE_nest, sep = " + ")),
data = subdata_ter)
} else {
sm_ter_raw <- glmmTMB(
as.formula(paste("ShiftRate_raw~ShiftMod_raw", RE_nest, sep = " + ")),
data = subdata_ter)
}
ter_sm_res <- summary(sm_ter_raw)
ter_sm_res <- data.frame(ter_sm_res$coefficients$cond)
ter_sm_res$Param <- rownames(ter_sm_res)
names(ter_sm_res) <- c("Estimate", "Std_Error", "z-value", "P-value", "Param")
res_glmm_simple_model <- rbind(
data.frame(Eco = "Ter",
ter_sm_res),
data.frame(Eco = "Mar",
mar_sm_res)
)
# data amount
if(params$data_amount=="all"){
if(shiftmodtype=="95"){
write.csv(res_glmm_simple_model,
here("Data/Output/res_glmm_sm95_all.csv"),
row.names = FALSE)
} else {
write.csv(res_glmm_simple_model,
here("Data/Output/res_glmm_sm_all.csv"),
row.names = FALSE)
}
} else {
if(shiftmodtype=="95"){
write.csv(res_glmm_simple_model,
here("Data/Output/res_glmm_sm95.csv"),
row.names = FALSE)
} else {
write.csv(res_glmm_simple_model,
here("Data/Output/res_glmm_sm.csv"),
row.names = FALSE)
}
}6.1.1.4 Combine
tab_model(sm_ter,sm_mar,
dv.labels = c("Terrestrial","Marine"),
show.icc = FALSE,
show.r2 = FALSE)| Terrestrial | Marine | |||||
|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | p | Estimates | CI | p |
| (Intercept) | -0.01 | -0.24 – 0.22 | 0.939 | 0.31 | -0.03 – 0.65 | 0.071 |
| ShiftMod | -0.02 | -0.04 – -0.00 | 0.025 | -0.01 | -0.12 – 0.09 | 0.781 |
| Random Effects | ||||||
| σ2 | 0.74 | 1.63 | ||||
| τ00 | 0.48 ID:class:Param | 0.79 ID:class:Param | ||||
| 0.02 class:Param | 0.03 class:Param | |||||
| 0.02 Param | 0.02 Param | |||||
| N | 94 ID | 44 ID | ||||
| 10 class | 15 class | |||||
| 3 Param | 3 Param | |||||
| Observations | 8638 | 937 | ||||
6.2 Simple model - Study-level (Shift ~ ShiftMod)
6.2.1 Terrestrial
new_data <- as_tibble(bioshifts) %>%
filter(Eco=="Ter") %>%
select(ID, sp_name_std, class, Param, bvel_lat, bvel95_lat, ShiftRate, Eco)%>%
group_by(ID, class, Param, Eco) %>%
mutate(Nsps = length(unique(sp_name_std)),
ShiftMod_raw = bvel_lat,
ShiftMod_raw95 = bvel95_lat,
ShiftRate_raw = ShiftRate) %>%
filter(Nsps>3) %>%
summarise(Nsps = mean(Nsps),
ShiftMod = mean(ShiftMod_raw),
ShiftMod95 = mean(ShiftMod_raw95),
ShiftRate = mean(ShiftRate_raw),
ShiftMod95_sd = sd(ShiftMod_raw95),
ShiftMod_sd = sd(ShiftMod_raw),
ShiftRate_sd = sd(ShiftRate_raw)) %>%
plyr::rename(replace = c(class = "Class"), warn_missing = FALSE) `summarise()` has grouped output by 'ID', 'class', 'Param'. You can override
using the `.groups` argument.# Fit a deming model
if(shiftmodtype=="95"){
sm_SA_deming_ter <- deming(
formula = as.formula("ShiftRate ~ ShiftMod95"),
xstd = ShiftMod95_sd,
ystd = ShiftRate_sd,
weights = Nsps,
data = new_data)
} else {
sm_SA_deming_ter <- deming(
formula = as.formula(sm_fixed_eff),
xstd = ShiftMod_sd,
ystd = ShiftRate_sd,
weights = Nsps,
data = new_data)
}
sm_SA_deming_ter
Call:
deming(formula = as.formula(sm_fixed_eff), data = new_data, weights = Nsps, xstd = ShiftMod_sd, ystd = ShiftRate_sd)
n= 128
Coef se(coef) lower 0.95 upper 0.95
Intercept 0.1150676 0.07690119 -0.03565591 0.2657912
Slope 4.9399763 1.35867386 2.27702445 7.6029281
Scale= 1.04201 # Fit a GLMM (more appropriate because contains the random effects)
if(shiftmodtype=="95"){
sm_SA_glmm_ter <- glmmTMB(
as.formula("ShiftRate ~ ShiftMod95"),
weights = Nsps,
new_data)
} else {
sm_SA_glmm_ter <- glmmTMB(
as.formula(sm_fixed_eff),
weights = Nsps,
new_data)
}
tab_model(sm_SA_glmm_ter)| Shift Rate | |||
| Predictors | Estimates | CI | p |
| (Intercept) | 0.61 | 0.58 – 0.64 | <0.001 |
| ShiftMod | 0.94 | 0.86 – 1.02 | <0.001 |
| Observations | 128 | ||
| R2 / R2 adjusted | 0.062 / 0.047 | ||
tmp <- summary(sm_SA_glmm_ter)
{
plot(ShiftRate~ShiftMod, new_data)
abline(h = 0, v = 0, lty = "dashed")
abline(sm_SA_deming_ter, col = "blue")
abline(b = tmp$coefficients$cond[2,1], a = tmp$coefficients$cond[1,1], col = "red")
}
No association between modeled shift velocity and observed shift velocity.
6.2.2 Marine
new_data <- as_tibble(bioshifts) %>%
filter(Eco=="Mar") %>%
select(ID, sp_name_std, class, Param, bvel_lat, bvel95_lat, ShiftRate, Eco)%>%
group_by(ID, class, Param, Eco) %>%
mutate(Nsps = length(unique(sp_name_std)),
ShiftMod_raw = bvel_lat,
ShiftMod_raw95 = bvel95_lat,
ShiftRate_raw = ShiftRate) %>%
filter(Nsps>3) %>%
summarise(Nsps = mean(Nsps),
ShiftMod = mean(ShiftMod_raw),
ShiftMod95 = mean(ShiftMod_raw95),
ShiftRate = mean(ShiftRate_raw),
ShiftMod95_sd = sd(ShiftMod_raw95),
ShiftMod_sd = sd(ShiftMod_raw),
ShiftRate_sd = sd(ShiftRate_raw)) %>%
plyr::rename(replace = c(class = "Class"), warn_missing = FALSE) `summarise()` has grouped output by 'ID', 'class', 'Param'. You can override
using the `.groups` argument.# new_data$ShiftMod_sd[which(new_data$ShiftMod_sd==0)] <- 0.0001
# new_data$ShiftRate_sd[which(new_data$ShiftMod_sd==0)] <- 0.0001
# new_data[is.na(new_data)] <- 0.0001
# Fit a deming model
# if(shiftmodtype=="95"){
# sm_SA_deming_mar <- deming(
# formula = as.formula("ShiftRate ~ ShiftMod95"),
# xstd = ShiftMod95_sd,
# ystd = ShiftRate_sd,
# weights = Nsps,
# data = new_data)
# } else {
# sm_SA_deming_mar <- deming(
# formula = as.formula(sm_fixed_eff),
# xstd = ShiftMod_sd,
# ystd = ShiftRate_sd,
# weights = Nsps,
# data = new_data)
# }
#
# sm_SA_deming_mar
# Fit a GLMM (more appropriate because contains the random effects)
if(shiftmodtype=="95"){
sm_SA_glmm_mar <- glmmTMB(
as.formula("ShiftRate ~ ShiftMod95"),
weights = Nsps,
new_data)
} else {
sm_SA_glmm_mar <- glmmTMB(
as.formula(sm_fixed_eff),
weights = Nsps,
new_data)
}
tab_model(sm_SA_glmm_mar)| Shift Rate | |||
| Predictors | Estimates | CI | p |
| (Intercept) | 1.75 | 1.56 – 1.94 | <0.001 |
| ShiftMod | 1.59 | 1.22 – 1.97 | <0.001 |
| Observations | 59 | ||
| R2 / R2 adjusted | 0.077 / 0.044 | ||
tmp <- summary(sm_SA_glmm_mar)
{
plot(ShiftRate~ShiftMod, new_data)
abline(h = 0, v = 0, lty = "dashed")
# abline(sm_SA_deming_mar, col = "blue")
abline(b = tmp$coefficients$cond[2,1], a = tmp$coefficients$cond[1,1], col = "red")
}
Same here…
6.2.3 Save model results
ter_sm_SA_glmm_res <- summary(sm_SA_glmm_ter)
ter_sm_SA_glmm_res <- data.frame(ter_sm_SA_glmm_res$coefficients$cond)
ter_sm_SA_glmm_res$Param <- rownames(ter_sm_SA_glmm_res)
names(ter_sm_SA_glmm_res) <- c("Estimate", "Std_Error", "z-value", "P-value", "Param")
mar_sm_SA_glmm_res <- summary(sm_SA_glmm_mar)
mar_sm_SA_glmm_res <- data.frame(mar_sm_SA_glmm_res$coefficients$cond)
mar_sm_SA_glmm_res$Param <- rownames(mar_sm_SA_glmm_res)
names(mar_sm_SA_glmm_res) <- c("Estimate", "Std_Error", "z-value", "P-value", "Param")
sm_SA_glmm_res <- rbind(data.frame(Eco = "Ter",ter_sm_SA_glmm_res),
data.frame(Eco = "Mar",mar_sm_SA_glmm_res))
# data amount
if(params$data_amount=="all"){
if(shiftmodtype=="95"){
write.csv(sm_SA_glmm_res,
here("Data/Output/res_sm_SA_glmm95_all.csv"),
row.names = FALSE)
} else {
write.csv(sm_SA_glmm_res,
here("Data/Output/res_sm_SA_glmm_all.csv"),
row.names = FALSE)
}
} else {
if(shiftmodtype=="95"){
write.csv(sm_SA_glmm_res,
here("Data/Output/res_sm_SA_glmm95.csv"),
row.names = FALSE)
} else {
write.csv(sm_SA_glmm_res,
here("Data/Output/res_sm_SA_glmm.csv"),
row.names = FALSE)
}
}6.2.4 Combine
tab_model(sm_SA_glmm_ter,sm_SA_glmm_mar,
dv.labels = c("Terrestrial","Marine"),
show.icc = FALSE,
show.r2 = FALSE)| Terrestrial | Marine | |||||
| Predictors | Estimates | CI | p | Estimates | CI | p |
| (Intercept) | 0.61 | 0.58 – 0.64 | <0.001 | 1.75 | 1.56 – 1.94 | <0.001 |
| ShiftMod | 0.94 | 0.86 – 1.02 | <0.001 | 1.59 | 1.22 – 1.97 | <0.001 |
| Observations | 128 | 59 | ||||
6.3 Simple model - direction (Shift ~ ShiftMod * dir)
sm_fixed_eff_dir <- "ShiftRate ~ ShiftMod * Match_dir"
sm_formula_dir <- paste(sm_fixed_eff_dir, RE_nest, sep = " + ")6.3.1 Terrestrial
var_to_go <- strsplit(sm_formula_dir,"~|[+]|[/]|[|]|[*]")[[1]]
var_to_go <- gsub(" |[(]|[)]","",var_to_go)
var_to_go <- var_to_go[-which(var_to_go=="1")]
subdata_ter <- mod_data %>%
filter(Eco=="Ter") %>%
select(var_to_go) %>%
na.omitWarning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
ℹ Please use `all_of()` or `any_of()` instead.
# Was:
data %>% select(var_to_go)
# Now:
data %>% select(all_of(var_to_go))
See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.sm_dir_ter <- glmmTMB(as.formula(sm_formula_dir),
data = subdata_ter)
tab_model(sm_dir_ter)| Shift Rate | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 0.37 | 0.14 – 0.60 | 0.002 |
| ShiftMod | 0.06 | 0.03 – 0.09 | <0.001 |
| Match dir [Equatorward] | -0.71 | -0.77 – -0.65 | <0.001 |
| Match dir [Mismatch] | -0.52 | -0.56 – -0.48 | <0.001 |
| ShiftMod × Match dir [Equatorward] |
0.16 | 0.09 – 0.23 | <0.001 |
| ShiftMod × Match dir [Mismatch] |
-0.35 | -0.39 – -0.31 | <0.001 |
| Random Effects | |||
| σ2 | 0.63 | ||
| τ00 ID:class:Param | 0.50 | ||
| τ00 class:Param | 0.00 | ||
| τ00 Param | 0.03 | ||
| N ID | 94 | ||
| N class | 10 | ||
| N Param | 3 | ||
| Observations | 8215 | ||
| Marginal R2 / Conditional R2 | 0.176 / NA | ||
plot(effects::allEffects(sm_dir_ter), main = "Terrestrial")
Modeled shift velocity is positively related to observed shift velocity only when they shift in the same direction.
6.3.2 Marine
var_to_go <- strsplit(sm_formula_dir,"~|[+]|[/]|[|]|[*]")[[1]]
var_to_go <- gsub(" |[(]|[)]","",var_to_go)
var_to_go <- var_to_go[-which(var_to_go=="1")]
subdata_mar <- mod_data %>%
filter(Eco=="Mar") %>%
select(var_to_go) %>%
na.omit
sm_dir_mar <- glmmTMB(as.formula(sm_formula_dir),
data = subdata_mar)
tab_model(sm_dir_mar)| Shift Rate | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 0.55 | 0.30 – 0.80 | <0.001 |
| ShiftMod | 0.53 | 0.31 – 0.75 | <0.001 |
| Match dir [Equatorward] | -0.75 | -1.01 – -0.50 | <0.001 |
| Match dir [Mismatch] | -0.23 | -0.45 – -0.01 | 0.044 |
| ShiftMod × Match dir [Equatorward] |
-0.14 | -0.52 – 0.24 | 0.464 |
| ShiftMod × Match dir [Mismatch] |
-1.14 | -1.40 – -0.88 | <0.001 |
| Random Effects | |||
| σ2 | 1.35 | ||
| τ00 ID:class:Param | 0.98 | ||
| τ00 class:Param | 0.00 | ||
| τ00 Param | 0.00 | ||
| N ID | 43 | ||
| N class | 15 | ||
| N Param | 3 | ||
| Observations | 852 | ||
| Marginal R2 / Conditional R2 | 0.235 / NA | ||
plot(effects::allEffects(sm_dir_mar), main = "Marine")
Modeled shift velocity is positively related to observed shift velocity only when they shift in the same direction.
6.3.3 Combine
tab_model(sm_dir_ter,sm_dir_mar,
dv.labels = c("Terrestrial","Marine"),
show.icc = FALSE,
show.r2 = FALSE)| Terrestrial | Marine | |||||
|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | p | Estimates | CI | p |
| (Intercept) | 0.37 | 0.14 – 0.60 | 0.002 | 0.55 | 0.30 – 0.80 | <0.001 |
| ShiftMod | 0.06 | 0.03 – 0.09 | <0.001 | 0.53 | 0.31 – 0.75 | <0.001 |
| Match dir [Equatorward] | -0.71 | -0.77 – -0.65 | <0.001 | -0.75 | -1.01 – -0.50 | <0.001 |
| Match dir [Mismatch] | -0.52 | -0.56 – -0.48 | <0.001 | -0.23 | -0.45 – -0.01 | 0.044 |
| ShiftMod × Match dir [Equatorward] |
0.16 | 0.09 – 0.23 | <0.001 | -0.14 | -0.52 – 0.24 | 0.464 |
| ShiftMod × Match dir [Mismatch] |
-0.35 | -0.39 – -0.31 | <0.001 | -1.14 | -1.40 – -0.88 | <0.001 |
| Random Effects | ||||||
| σ2 | 0.63 | 1.35 | ||||
| τ00 | 0.50 ID:class:Param | 0.98 ID:class:Param | ||||
| 0.00 class:Param | 0.00 class:Param | |||||
| 0.03 Param | 0.00 Param | |||||
| N | 94 ID | 43 ID | ||||
| 10 class | 15 class | |||||
| 3 Param | 3 Param | |||||
| Observations | 8215 | 852 | ||||
6.3.4 Chi-squared test for direction frequency
# Terrestrial
dir_table <- mod_data %>%
filter(Eco == "Ter",
ShiftRate!=0) %>%
group_by(match_direction) %>%
count %>%
na.omit
dir_table# A tibble: 2 × 2
# Groups: match_direction [2]
match_direction n
<fct> <int>
1 mismatch 3814
2 match 4459chisq.test(dir_table$n)
Chi-squared test for given probabilities
data: dir_table$n
X-squared = 50.287, df = 1, p-value = 1.328e-12# Marine
dir_table <- mod_data %>%
filter(Eco == "Mar",
ShiftRate!=0) %>%
group_by(match_direction) %>%
count %>%
na.omit
dir_table# A tibble: 2 × 2
# Groups: match_direction [2]
match_direction n
<fct> <int>
1 mismatch 230
2 match 696chisq.test(dir_table$n)
Chi-squared test for given probabilities
data: dir_table$n
X-squared = 234.51, df = 1, p-value < 2.2e-166.4 Binomial model
6.4.0.1 Terrestrial
subdata_ter <- mod_data %>%
filter(Eco=="Ter") %>%
select(mismatch_binary,Param,class,ID) %>%
na.omit
# fit model
mm01_ter <- glmmTMB(mismatch_binary ~ 1 + (1 | Param/class/ID),
family = binomial,
data = subdata_ter)
tab_model(mm01_ter, transform = NULL, string.est = "Estimate")| mismatch binary | |||
|---|---|---|---|
| Predictors | Estimate | CI | p |
| (Intercept) | 0.11 | 0.01 – 0.21 | 0.032 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID:class:Param | 0.20 | ||
| τ00 class:Param | 0.00 | ||
| τ00 Param | 0.00 | ||
| N ID | 94 | ||
| N class | 10 | ||
| N Param | 3 | ||
| Observations | 8273 | ||
| Marginal R2 / Conditional R2 | 0.000 / NA | ||
6.4.0.2 Marine
subdata_mar <- mod_data %>%
filter(Eco=="Mar") %>%
select(mismatch_binary,Param,class,ID) %>%
na.omit
# fit model
mm01_mar <- glmmTMB(mismatch_binary ~ 1 + (1 | Param/class/ID),
family = binomial,
data = subdata_mar)
tab_model(mm01_mar, transform = NULL, string.est = "Estimate")| mismatch binary | |||
|---|---|---|---|
| Predictors | Estimate | CI | p |
| (Intercept) | 1.21 | NaN – NaN | NaN |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID:class:Param | 0.70 | ||
| τ00 class:Param | 0.00 | ||
| τ00 Param | 0.00 | ||
| N ID | 44 | ||
| N class | 15 | ||
| N Param | 3 | ||
| Observations | 926 | ||
| Marginal R2 / Conditional R2 | 0.000 / NA | ||
6.5 Mismatch between documented and modeled shifts
6.5.1 Terrestrial
y <- strsplit(mm_formula_ter," ")[[1]][1]
var_to_go <- strsplit(paste(mm_formula_ter, "+", RE_nest),"~|[+]|[/]|[|]")[[1]]
var_to_go <- gsub(" |[(]|[)]","",var_to_go)
var_to_go <- var_to_go[-which(var_to_go=="1")]
subdata_ter <- mod_data %>%
filter(Eco=="Ter") %>%
select(var_to_go) %>%
na.omit
nrow(subdata_ter)[1] 8604# fit model
mm_ter <- glmmTMB(as.formula(paste(mm_formula_ter, "+", RE_nest)),
family = Gamma(link = "log"),
data = subdata_ter)
# simulationOutput <- DHARMa::simulateResiduals(fittedModel = mm_ter, plot = F)
# DHARMa::testResiduals(simulationOutput)
tab_model(mm_ter, transform = NULL, string.est = "Estimate")| mismatch | |||
|---|---|---|---|
| Predictors | Estimate | CI | p |
| (Intercept) | 0.83 | 0.59 – 1.07 | <0.001 |
| TSS | -0.02 | -0.04 – 0.01 | 0.243 |
| Study duration | -0.19 | -0.29 – -0.09 | <0.001 |
| Interannual var | 0.16 | 0.12 – 0.20 | <0.001 |
| Lat extent | 0.38 | 0.24 – 0.51 | <0.001 |
| Connectivity | -0.11 | -0.15 – -0.06 | <0.001 |
| Sampling category [Survey-Resurvey] |
-0.00 | -0.44 – 0.43 | 0.987 |
| Sampling category [TimeSeries] |
0.16 | -0.19 – 0.50 | 0.369 |
| Locomotion mode [Run/Walk] |
-0.36 | -0.79 – 0.07 | 0.101 |
| Locomotion mode [Swimming] |
-0.38 | -1.13 – 0.38 | 0.326 |
| Locomotion mode [Slithering/Crawling] |
-0.08 | -0.68 – 0.52 | 0.796 |
| Locomotion mode [Sessile] | -0.97 | -1.33 – -0.61 | <0.001 |
| Random Effects | |||
| σ2 | 1.08 | ||
| τ00 ID:class:Param | 0.28 | ||
| τ00 class:Param | 0.05 | ||
| τ00 Param | 0.00 | ||
| N ID | 94 | ||
| N class | 10 | ||
| N Param | 3 | ||
| Observations | 8604 | ||
| Marginal R2 / Conditional R2 | 0.440 / NA | ||
# get emmeans for categorical variables
files_cat_cont_names <- gsub("_"," ",cat_vars)
mm_ter_emm <- lapply(cat_vars, function(x){
if(any(grepl(x,var_to_go))){
emmeans(object = mm_ter,
specs = x,
type = "response")
}
})
names(mm_ter_emm) <- files_cat_cont_names
test <- sapply(mm_ter_emm,is.null)
if(any(test)){
mm_ter_emm <- mm_ter_emm[-which(test)]
}
mm_ter_emm$`Sampling category`
Sampling_category response SE df asymp.LCL asymp.UCL
CensusPeriods 1.62 0.218 Inf 1.245 2.11
Survey-Resurvey 1.61 0.404 Inf 0.989 2.64
TimeSeries 1.90 0.416 Inf 1.235 2.92
Results are averaged over the levels of: Locomotion_mode
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$`Locomotion mode`
Locomotion_mode response SE df asymp.LCL asymp.UCL
Flying 2.439 0.368 Inf 1.814 3.28
Run/Walk 1.698 0.379 Inf 1.096 2.63
Swimming 1.670 0.640 Inf 0.788 3.54
Slithering/Crawling 2.254 0.693 Inf 1.234 4.12
Sessile 0.927 0.138 Inf 0.693 1.24
Results are averaged over the levels of: Sampling_category
Confidence level used: 0.95
Intervals are back-transformed from the log scale mm_ter_cont <- lapply(mm_ter_emm,
contrast,
method = "eff",
type = "response")
names(mm_ter_cont) <- names(mm_ter_cont)
mm_ter_cont$`Sampling category`
contrast ratio SE df null z.ratio p.value
CensusPeriods effect 0.950 0.0944 Inf 1 -0.518 0.7214
(Survey-Resurvey) effect 0.946 0.1463 Inf 1 -0.357 0.7214
TimeSeries effect 1.112 0.1472 Inf 1 0.805 0.7214
Results are averaged over the levels of: Locomotion_mode
P value adjustment: fdr method for 3 tests
Tests are performed on the log scale
$`Locomotion mode`
contrast ratio SE df null z.ratio p.value
Flying effect 1.430 0.2186 Inf 1 2.337 0.0486
(Run/Walk) effect 0.995 0.1414 Inf 1 -0.035 0.9720
Swimming effect 0.979 0.2770 Inf 1 -0.075 0.9720
(Slithering/Crawling) effect 1.321 0.3132 Inf 1 1.174 0.4007
Sessile effect 0.544 0.0898 Inf 1 -3.692 0.0011
Results are averaged over the levels of: Sampling_category
P value adjustment: fdr method for 5 tests
Tests are performed on the log scale mm_ter_cont_pairs <- lapply(mm_ter_emm,
contrast,
method = "pairwise",
type = "response")
names(mm_ter_cont_pairs) <- names(mm_ter_cont)
mm_ter_cont_pairs$`Sampling category`
contrast ratio SE df null z.ratio p.value
CensusPeriods / (Survey-Resurvey) 1.004 0.224 Inf 1 0.016 0.9999
CensusPeriods / TimeSeries 0.854 0.150 Inf 1 -0.899 0.6411
(Survey-Resurvey) / TimeSeries 0.851 0.230 Inf 1 -0.599 0.8209
Results are averaged over the levels of: Locomotion_mode
P value adjustment: tukey method for comparing a family of 3 estimates
Tests are performed on the log scale
$`Locomotion mode`
contrast ratio SE df null z.ratio p.value
Flying / (Run/Walk) 1.437 0.317 Inf 1 1.642 0.4704
Flying / Swimming 1.460 0.563 Inf 1 0.982 0.8636
Flying / (Slithering/Crawling) 1.082 0.330 Inf 1 0.259 0.9990
Flying / Sessile 2.630 0.482 Inf 1 5.275 <.0001
(Run/Walk) / Swimming 1.016 0.337 Inf 1 0.049 1.0000
(Run/Walk) / (Slithering/Crawling) 0.753 0.231 Inf 1 -0.923 0.8881
(Run/Walk) / Sessile 1.831 0.430 Inf 1 2.577 0.0747
Swimming / (Slithering/Crawling) 0.741 0.329 Inf 1 -0.674 0.9620
Swimming / Sessile 1.801 0.701 Inf 1 1.513 0.5540
(Slithering/Crawling) / Sessile 2.430 0.782 Inf 1 2.760 0.0457
Results are averaged over the levels of: Sampling_category
P value adjustment: tukey method for comparing a family of 5 estimates
Tests are performed on the log scale # marginal R2 model avg
Actual = subdata_ter[,y]
Predicted = predict(mm_ter,
type="response",
re.form=NA)
Residuals = Actual - Predicted
mm_ter_R2_avg_marginal <- efronRSquared(residual = Residuals,
predicted = Predicted,
statistic = "EfronRSquared")
mm_ter_R2_avg_marginalEfronRSquared
0.196 # conditional R2 model avg
Actual = subdata_ter[,y]
Predicted = predict(mm_ter, type="response")
Residuals = Actual - Predicted
mm_ter_R2_avg_conditional <- efronRSquared(residual = Residuals,
predicted = Predicted,
statistic = "EfronRSquared")
mm_ter_R2_avg_conditionalEfronRSquared
0.301 6.5.2 Marine
y <- strsplit(mm_formula_mar," ")[[1]][1]
var_to_go <- strsplit(paste(mm_formula_mar, "+", RE_nest),"~|[+]|[/]|[|]")[[1]]
var_to_go <- gsub(" |[(]|[)]","",var_to_go)
var_to_go <- var_to_go[-which(var_to_go=="1")]
subdata_mar <- mod_data %>%
filter(Eco=="Mar") %>%
select(var_to_go) %>%
na.omit
# fit model
mm_mar <- glmmTMB(as.formula(paste(mm_formula_mar, "+", RE_nest)),
family = Gamma(link = "log"),
data = subdata_mar)
# simulationOutput <- DHARMa::simulateResiduals(fittedModel = mm_mar, plot = F)
# DHARMa::testResiduals(simulationOutput)
tab_model(mm_mar, transform = NULL, string.est = "Estimate")| mismatch | |||
|---|---|---|---|
| Predictors | Estimate | CI | p |
| (Intercept) | -0.30 | -1.27 – 0.66 | 0.538 |
| TSS | 0.06 | -0.01 – 0.12 | 0.088 |
| Study duration | -0.35 | -0.60 – -0.09 | 0.007 |
| Interannual var | 0.16 | 0.07 – 0.26 | 0.001 |
| Lat extent | 0.63 | 0.48 – 0.78 | <0.001 |
| Sampling category [Survey-Resurvey] |
-0.25 | -0.69 – 0.18 | 0.256 |
| Sampling category [TimeSeries] |
0.28 | 0.00 – 0.56 | 0.048 |
| Life form [Ectotherm] | 1.57 | 0.65 – 2.49 | 0.001 |
| Life form [Endotherm] | 1.00 | -0.21 – 2.21 | 0.107 |
| Locomotion mode [Slithering/Crawling] |
-0.49 | -0.96 – -0.02 | 0.041 |
| Locomotion mode [Planktonic] |
1.21 | 0.13 – 2.29 | 0.028 |
| Locomotion mode [Sessile] | 0.31 | -0.19 – 0.81 | 0.231 |
| Random Effects | |||
| σ2 | 1.04 | ||
| τ00 ID:class:Param | 0.16 | ||
| τ00 class:Param | 0.00 | ||
| τ00 Param | 0.06 | ||
| N ID | 44 | ||
| N class | 15 | ||
| N Param | 3 | ||
| Observations | 937 | ||
| Marginal R2 / Conditional R2 | 0.276 / NA | ||
# get emmeans for categorical variables
files_cat_cont_names <- gsub("_"," ",cat_vars)
mm_mar_emm <- lapply(cat_vars, function(x){
if(any(grepl(x,var_to_go))){
emmeans(object = mm_mar,
specs = x,
type = "response")
}
})
names(mm_mar_emm) <- files_cat_cont_names
test <- sapply(mm_mar_emm,is.null)
if(any(test)){
mm_mar_emm <- mm_mar_emm[-which(test)]
}
mm_mar_emm$`Life form`
Life_form response SE df asymp.LCL asymp.UCL
Cryptogam 0.861 0.333 Inf 0.404 1.84
Ectotherm 4.152 0.966 Inf 2.631 6.55
Endotherm 2.331 1.072 Inf 0.947 5.74
Results are averaged over the levels of: Sampling_category, Locomotion_mode
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$`Sampling category`
Sampling_category response SE df asymp.LCL asymp.UCL
CensusPeriods 2.01 0.486 Inf 1.248 3.23
Survey-Resurvey 1.56 0.449 Inf 0.887 2.74
TimeSeries 2.66 0.651 Inf 1.649 4.30
Results are averaged over the levels of: Life_form, Locomotion_mode
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$`Locomotion mode`
Locomotion_mode response SE df asymp.LCL asymp.UCL
Swimming 1.567 0.413 Inf 0.935 2.63
Slithering/Crawling 0.961 0.334 Inf 0.486 1.90
Planktonic 5.274 2.509 Inf 2.076 13.40
Sessile 2.127 0.615 Inf 1.207 3.75
Results are averaged over the levels of: Sampling_category, Life_form
Confidence level used: 0.95
Intervals are back-transformed from the log scale mm_mar_cont <- lapply(mm_mar_emm,
contrast,
method = "eff",
type = "response")
names(mm_mar_cont) <- names(mm_mar_cont)
mm_mar_cont$`Life form`
contrast ratio SE df null z.ratio p.value
Cryptogam effect 0.425 0.145 Inf 1 -2.515 0.0179
Ectotherm effect 2.048 0.417 Inf 1 3.515 0.0013
Endotherm effect 1.150 0.355 Inf 1 0.452 0.6512
Results are averaged over the levels of: Sampling_category, Locomotion_mode
P value adjustment: fdr method for 3 tests
Tests are performed on the log scale
$`Sampling category`
contrast ratio SE df null z.ratio p.value
CensusPeriods effect 0.99 0.0988 Inf 1 -0.103 0.9180
(Survey-Resurvey) effect 0.77 0.1076 Inf 1 -1.874 0.0915
TimeSeries effect 1.31 0.1314 Inf 1 2.719 0.0196
Results are averaged over the levels of: Life_form, Locomotion_mode
P value adjustment: fdr method for 3 tests
Tests are performed on the log scale
$`Locomotion mode`
contrast ratio SE df null z.ratio p.value
Swimming effect 0.773 0.146 Inf 1 -1.367 0.2287
(Slithering/Crawling) effect 0.474 0.113 Inf 1 -3.119 0.0073
Planktonic effect 2.601 1.033 Inf 1 2.407 0.0321
Sessile effect 1.049 0.193 Inf 1 0.261 0.7941
Results are averaged over the levels of: Sampling_category, Life_form
P value adjustment: fdr method for 4 tests
Tests are performed on the log scale mm_mar_cont_pairs <- lapply(mm_mar_emm,
contrast,
method = "pairwise",
type = "response")
names(mm_mar_cont_pairs) <- names(mm_mar_cont)
mm_mar_cont_pairs$`Life form`
contrast ratio SE df null z.ratio p.value
Cryptogam / Ectotherm 0.207 0.0972 Inf 1 -3.356 0.0023
Cryptogam / Endotherm 0.369 0.2280 Inf 1 -1.613 0.2399
Ectotherm / Endotherm 1.781 0.7075 Inf 1 1.453 0.3139
Results are averaged over the levels of: Sampling_category, Locomotion_mode
P value adjustment: tukey method for comparing a family of 3 estimates
Tests are performed on the log scale
$`Sampling category`
contrast ratio SE df null z.ratio p.value
CensusPeriods / (Survey-Resurvey) 1.286 0.285 Inf 1 1.137 0.4913
CensusPeriods / TimeSeries 0.754 0.108 Inf 1 -1.976 0.1181
(Survey-Resurvey) / TimeSeries 0.586 0.130 Inf 1 -2.409 0.0423
Results are averaged over the levels of: Life_form, Locomotion_mode
P value adjustment: tukey method for comparing a family of 3 estimates
Tests are performed on the log scale
$`Locomotion mode`
contrast ratio SE df null z.ratio p.value
Swimming / (Slithering/Crawling) 1.631 0.390 Inf 1 2.045 0.1714
Swimming / Planktonic 0.297 0.164 Inf 1 -2.201 0.1230
Swimming / Sessile 0.737 0.188 Inf 1 -1.198 0.6278
(Slithering/Crawling) / Planktonic 0.182 0.109 Inf 1 -2.852 0.0225
(Slithering/Crawling) / Sessile 0.452 0.143 Inf 1 -2.512 0.0581
Planktonic / Sessile 2.479 1.251 Inf 1 1.799 0.2739
Results are averaged over the levels of: Sampling_category, Life_form
P value adjustment: tukey method for comparing a family of 4 estimates
Tests are performed on the log scale # marginal R2 model avg
Actual = subdata_mar[,y]
Predicted = predict(mm_mar,
type="response",
re.form=NA)
Residuals = Actual - Predicted
mm_mar_R2_avg_marginal <- efronRSquared(residual = Residuals,
predicted = Predicted,
statistic = "EfronRSquared")
mm_mar_R2_avg_marginalEfronRSquared
0.163 # conditional R2 model avg
Actual = subdata_mar[,y]
Predicted = predict(mm_mar, type="response")
Residuals = Actual - Predicted
mm_mar_R2_avg_conditional <- efronRSquared(residual = Residuals,
predicted = Predicted,
statistic = "EfronRSquared")
mm_mar_R2_avg_conditionalEfronRSquared
0.318 6.5.3 Combine
tab_model(mm_ter, mm_mar,
transform = NULL, string.est = "Estimate",
dv.labels = c("Terrestrial","Marine"),
show.icc = FALSE,
show.r2 = FALSE)| Terrestrial | Marine | |||||
|---|---|---|---|---|---|---|
| Predictors | Estimate | CI | p | Estimate | CI | p |
| (Intercept) | 0.83 | 0.59 – 1.07 | <0.001 | -0.30 | -1.27 – 0.66 | 0.538 |
| TSS | -0.02 | -0.04 – 0.01 | 0.243 | 0.06 | -0.01 – 0.12 | 0.088 |
| Study duration | -0.19 | -0.29 – -0.09 | <0.001 | -0.35 | -0.60 – -0.09 | 0.007 |
| Interannual var | 0.16 | 0.12 – 0.20 | <0.001 | 0.16 | 0.07 – 0.26 | 0.001 |
| Lat extent | 0.38 | 0.24 – 0.51 | <0.001 | 0.63 | 0.48 – 0.78 | <0.001 |
| Connectivity | -0.11 | -0.15 – -0.06 | <0.001 | |||
| Sampling category [Survey-Resurvey] |
-0.00 | -0.44 – 0.43 | 0.987 | -0.25 | -0.69 – 0.18 | 0.256 |
| Sampling category [TimeSeries] |
0.16 | -0.19 – 0.50 | 0.369 | 0.28 | 0.00 – 0.56 | 0.048 |
| Locomotion mode [Run/Walk] |
-0.36 | -0.79 – 0.07 | 0.101 | |||
| Locomotion mode [Swimming] |
-0.38 | -1.13 – 0.38 | 0.326 | |||
| Locomotion mode [Slithering/Crawling] |
-0.08 | -0.68 – 0.52 | 0.796 | -0.49 | -0.96 – -0.02 | 0.041 |
| Locomotion mode [Sessile] | -0.97 | -1.33 – -0.61 | <0.001 | 0.31 | -0.19 – 0.81 | 0.231 |
| Life form [Ectotherm] | 1.57 | 0.65 – 2.49 | 0.001 | |||
| Life form [Endotherm] | 1.00 | -0.21 – 2.21 | 0.107 | |||
| Locomotion mode [Planktonic] |
1.21 | 0.13 – 2.29 | 0.028 | |||
| Random Effects | ||||||
| σ2 | 1.08 | 1.04 | ||||
| τ00 | 0.28 ID:class:Param | 0.16 ID:class:Param | ||||
| 0.05 class:Param | 0.00 class:Param | |||||
| 0.00 Param | 0.06 Param | |||||
| N | 94 ID | 44 ID | ||||
| 10 class | 15 class | |||||
| 3 Param | 3 Param | |||||
| Observations | 8604 | 937 | ||||
7 Compile model results
7.1 Terrestrial
# get results from continuous variables
cont_vars_res_ter <- summary(mm_ter)
cont_vars_res_ter <- data.frame(cont_vars_res_ter$coefficients$cond)
confint_mm_ter <- confint(mm_ter)[,1:2]
confint_mm_ter <- confint_mm_ter[rownames(cont_vars_res_ter),]
cont_vars_res_ter <- cbind(cont_vars_res_ter, confint_mm_ter)
cont_vars_res_ter$var <- rownames(cont_vars_res_ter)
rownames(cont_vars_res_ter) <- NULL
cont_vars_res_ter <- cont_vars_res_ter %>% filter(var %in% c("Int",cont_vars))
names(cont_vars_res_ter) <- c("Estimate","Std.Error","z.value","P-value","Lower.CI","Upper.CI","var")
# data amount
if(params$data_amount=="all"){
write.csv(cont_vars_res_ter,
here("Data/Output/cont_vars_res_ter_all.csv"),
row.names = TRUE)
} else {
write.csv(cont_vars_res_ter,
here("Data/Output/cont_vars_res_ter.csv"),
row.names = TRUE)
}
# get results from categorical variables
# contrasts
cat_vars_cont_res_ter <- lapply(1:length(mm_ter_cont), function(x){
data.frame(confint(mm_ter_cont[[x]]),var=names(mm_ter_cont)[x])
})
rem <- which(sapply(cat_vars_cont_res_ter,is.null))
cat_vars_cont_res_ter <- rbindlist(cat_vars_cont_res_ter)
cat_vars_cont_res_ter$level <- cat_vars_cont_res_ter$contrast
cat_vars_cont_res_ter$level <- gsub(" effect","",cat_vars_cont_res_ter$level)
cat_vars_cont_res_ter <- cat_vars_cont_res_ter[,-1]
names(cat_vars_cont_res_ter) <- c("Estimate","SE","df","Lower.CI","Upper.CI","var","level")
# data amount
if(params$data_amount=="all"){
write.csv(cat_vars_cont_res_ter,
here("Data/Output/cat_vars_cont_res_ter_all.csv"),
row.names = TRUE)
} else {
write.csv(cat_vars_cont_res_ter,
here("Data/Output/cat_vars_cont_res_ter.csv"),
row.names = TRUE)
}
# pairwise contrasts
cat_vars_cont_pairs_res_ter <- lapply(1:length(mm_ter_cont_pairs), function(x){
data.frame(mm_ter_cont_pairs[[x]],var=names(mm_ter_cont_pairs)[x])
})
rem <- which(sapply(cat_vars_cont_pairs_res_ter,is.null))
cat_vars_cont_pairs_res_ter <- rbindlist(cat_vars_cont_pairs_res_ter)
cat_vars_cont_pairs_res_ter$level <- cat_vars_cont_pairs_res_ter$contrast
cat_vars_cont_pairs_res_ter$level <- gsub(" effect","",cat_vars_cont_pairs_res_ter$level)
cat_vars_cont_pairs_res_ter <- cat_vars_cont_pairs_res_ter[,-1]
names(cat_vars_cont_pairs_res_ter) <- c("Estimate","SE","df","null","z-ratio","P-value","var","level")
# data amount
if(params$data_amount=="all"){
write.csv(cat_vars_cont_pairs_res_ter,
here("Data/Output/cat_vars_cont_pairs_res_ter_all.csv"),
row.names = TRUE)
} else {
write.csv(cat_vars_cont_pairs_res_ter,
here("Data/Output/cat_vars_cont_pairs_res_ter.csv"),
row.names = TRUE)
}
# emmeans
cat_emm_vars_res_ter <- lapply(1:length(mm_ter_emm), function(x){
data.frame(confint(mm_ter_emm[[x]]),var=names(mm_ter_emm)[x])
})
cat_emm_vars_res_ter <- rbindlist(cat_emm_vars_res_ter)Column 1 ['Locomotion_mode'] of item 2 is missing in item 1. Use fill=TRUE to fill with NA (NULL for list columns), or use.names=FALSE to ignore column names. use.names='check' (default from v1.12.2) emits this message and proceeds as if use.names=FALSE for backwards compatibility. See news item 5 in v1.12.2 for options to control this message.cat_emm_vars_res_ter$level <- cat_emm_vars_res_ter[,1]
cat_emm_vars_res_ter <- cat_emm_vars_res_ter[,-1]
cat_emm_vars_res_ter$level <- gsub(" effect","",cat_emm_vars_res_ter$level)
names(cat_emm_vars_res_ter) <- c("Estimate","SE","df","Lower.CI","Upper.CI","var","level")
# data amount
if(params$data_amount=="all"){
write.csv(cat_emm_vars_res_ter,
here("Data/Output/cat_emm_vars_res_ter_all.csv"),
row.names = TRUE)
} else {
write.csv(cat_emm_vars_res_ter,
here("Data/Output/cat_emm_vars_res_ter.csv"),
row.names = TRUE)
}7.2 Marine
# get results from continuous variables
cont_vars_res_mar <- summary(mm_mar)
cont_vars_res_mar <- data.frame(cont_vars_res_mar$coefficients$cond)
confint_mm_mar <- confint(mm_mar)[,1:2]
confint_mm_mar <- confint_mm_mar[rownames(cont_vars_res_mar),]
cont_vars_res_mar <- cbind(cont_vars_res_mar, confint_mm_mar)
cont_vars_res_mar$var <- rownames(cont_vars_res_mar)
rownames(cont_vars_res_mar) <- NULL
cont_vars_res_mar <- cont_vars_res_mar %>% filter(var %in% c("Int",cont_vars))
names(cont_vars_res_mar) <- c("Estimate","Std.Error","z.value","P-value","Lower.CI","Upper.CI","var")
# data amount
if(params$data_amount=="all"){
write.csv(cont_vars_res_mar,
here("Data/Output/cont_vars_res_mar_all.csv"),
row.names = TRUE)
} else {
write.csv(cont_vars_res_mar,
here("Data/Output/cont_vars_res_mar.csv"),
row.names = TRUE)
}
# get results from categorical variables
# contrasts
cat_vars_cont_res_mar <- lapply(1:length(mm_mar_cont), function(x){
data.frame(confint(mm_mar_cont[[x]]),var=names(mm_mar_cont)[x])
})
rem <- which(sapply(cat_vars_cont_res_mar,is.null))
cat_vars_cont_res_mar <- rbindlist(cat_vars_cont_res_mar)
cat_vars_cont_res_mar$level <- cat_vars_cont_res_mar$contrast
cat_vars_cont_res_mar$level <- gsub(" effect","",cat_vars_cont_res_mar$level)
cat_vars_cont_res_mar <- cat_vars_cont_res_mar[,-1]
names(cat_vars_cont_res_mar) <- c("Estimate","SE","df","Lower.CI","Upper.CI","var","level")
# data amount
if(params$data_amount=="all"){
write.csv(cat_vars_cont_res_mar,
here("Data/Output/cat_vars_cont_res_mar_all.csv"),
row.names = TRUE)
} else {
write.csv(cat_vars_cont_res_mar,
here("Data/Output/cat_vars_cont_res_mar.csv"),
row.names = TRUE)
}
# pairwise contrasts
cat_vars_cont_pairs_res_mar <- lapply(1:length(mm_mar_cont_pairs), function(x){
data.frame(mm_mar_cont_pairs[[x]],var=names(mm_mar_cont_pairs)[x])
})
rem <- which(sapply(cat_vars_cont_pairs_res_mar,is.null))
cat_vars_cont_pairs_res_mar <- rbindlist(cat_vars_cont_pairs_res_mar)
cat_vars_cont_pairs_res_mar$level <- cat_vars_cont_pairs_res_mar$contrast
cat_vars_cont_pairs_res_mar$level <- gsub(" effect","",cat_vars_cont_pairs_res_mar$level)
cat_vars_cont_pairs_res_mar <- cat_vars_cont_pairs_res_mar[,-1]
names(cat_vars_cont_pairs_res_mar) <- c("Estimate","SE","df","null","z-ratio","P-value","var","level")
# data amount
if(params$data_amount=="all"){
write.csv(cat_vars_cont_pairs_res_mar,
here("Data/Output/cat_vars_cont_pairs_res_mar_all.csv"),
row.names = TRUE)
} else {
write.csv(cat_vars_cont_pairs_res_mar,
here("Data/Output/cat_vars_cont_pairs_res_mar.csv"),
row.names = TRUE)
}
# emmeans
cat_emm_vars_res_mar <- lapply(1:length(mm_mar_emm), function(x){
data.frame(confint(mm_mar_emm[[x]]),var=names(mm_mar_emm)[x])
})
cat_emm_vars_res_mar <- rbindlist(cat_emm_vars_res_mar)Column 1 ['Sampling_category'] of item 2 is missing in item 1. Use fill=TRUE to fill with NA (NULL for list columns), or use.names=FALSE to ignore column names. use.names='check' (default from v1.12.2) emits this message and proceeds as if use.names=FALSE for backwards compatibility. See news item 5 in v1.12.2 for options to control this message.cat_emm_vars_res_mar$level <- cat_emm_vars_res_mar[,1]
cat_emm_vars_res_mar <- cat_emm_vars_res_mar[,-1]
cat_emm_vars_res_mar$level <- gsub(" effect","",cat_emm_vars_res_mar$level)
names(cat_emm_vars_res_mar) <- c("Estimate","SE","df","Lower.CI","Upper.CI","var","level")
# data amount
if(params$data_amount=="all"){
write.csv(cat_emm_vars_res_mar,
here("Data/Output/cat_emm_vars_res_mar_all.csv"),
row.names = TRUE)
} else {
write.csv(cat_emm_vars_res_mar,
here("Data/Output/cat_emm_vars_res_mar.csv"),
row.names = TRUE)
}8 Plot coeffs
8.1 Effects contrast
mm_mag_coefs_cont_ter <- data.frame(Eco = "Terrestrial",
rbindlist(list(cont_vars_res_ter,
cat_vars_cont_res_ter),
fill = TRUE))
mm_mag_coefs_cont_mar <- data.frame(Eco = "Marine",
rbindlist(list(cont_vars_res_mar,
cat_vars_cont_res_mar),
fill = TRUE))
mag_coefs_cont <- rbind(mm_mag_coefs_cont_ter,
mm_mag_coefs_cont_mar)
mag_coefs_cont$level[which(is.na(mag_coefs_cont$level))] <- mag_coefs_cont$var[which(is.na(mag_coefs_cont$level))]
sign_test <- sapply(1:nrow(mag_coefs_cont), function(x){
sign(mag_coefs_cont$Lower.CI[x])==sign(mag_coefs_cont$Upper.CI[x])
})
mag_coefs_cont$sig <- ifelse(sign_test, 1, 0)
mag_coefs_cont$level <- factor(mag_coefs_cont$level, levels = rev(unique(mag_coefs_cont$level)))
mag_coefs_cont$var <- factor(mag_coefs_cont$var, levels = (unique(mag_coefs_cont$var)))
mag_coefs_cont_plot_data <- mag_coefs_cont %>% filter(!var == "Int")
mag_coefs_cont_plot <- ggplot(mag_coefs_cont_plot_data,
aes(x = Estimate, y = level, alpha = sig)) +
geom_vline(xintercept=0, color = "gray",linetype="dashed", lwd=1.5)+
geom_point(size=6, shape = 21, color = "black", fill = "gray")+
geom_errorbar(aes(xmin=Lower.CI,xmax=Upper.CI), width=0, lwd=1.5) +
theme_classic(base_size = 20)+
# coord_cartesian(xlim = c(-1,1))+
xlab("Coefficient\n<<< Match - Mismatch >>>")+ylab("")+
theme(legend.position = "none")+
facet_grid(var~Eco, scales = "free_y", space = "free")
mag_coefs_cont_plot
8.2 Effects pairwise
mm_mag_coefs_cont_pair_ter <- data.frame(Eco = "Terrestrial",
rbindlist(list(cont_vars_res_ter,
cat_vars_cont_pairs_res_ter),
fill = TRUE))
mm_mag_coefs_cont_pair_mar <- data.frame(Eco = "Marine",
rbindlist(list(cont_vars_res_mar,
cat_vars_cont_pairs_res_mar),
fill = TRUE))
mag_coefs_cont_pair <- rbind(mm_mag_coefs_cont_pair_ter,
mm_mag_coefs_cont_pair_mar)
mag_coefs_cont_pair$level[which(is.na(mag_coefs_cont_pair$level))] <- mag_coefs_cont_pair$var[which(is.na(mag_coefs_cont_pair$level))]
sign_test <- sapply(1:nrow(mag_coefs_cont_pair), function(x){
mag_coefs_cont_pair$P.value[x] < 0.05
})
mag_coefs_cont_pair$sig <- ifelse(sign_test, 1, 0)
mag_coefs_cont_pair$level <- factor(mag_coefs_cont_pair$level, levels = rev(unique(mag_coefs_cont_pair$level)))
mag_coefs_cont_pair$var <- factor(mag_coefs_cont_pair$var, levels = (unique(mag_coefs_cont_pair$var)))
mag_coefs_cont_pair_plot_data <- mag_coefs_cont_pair %>% filter(!var == "Int")
mag_coefs_cont_pair_plot <- ggplot(mag_coefs_cont_pair_plot_data,
aes(x = Estimate, y = level, alpha = sig)) +
geom_vline(xintercept=0, color = "gray",linetype="dashed", lwd=1.5)+
geom_point(size=6, shape = 21, color = "black", fill = "gray")+
geom_errorbar(aes(xmin=Lower.CI,xmax=Upper.CI), width=0, lwd=1.5) +
theme_classic(base_size = 20)+
coord_cartesian(xlim = c(-1,1))+
xlab("Coefficient\n<<< Match - Mismatch >>>")+ylab("")+
theme(legend.position = "none")+
facet_grid(var~Eco, scales = "free_y", space = "free")
mag_coefs_cont_pair_plot
8.3 Effects Emmeans
mm_mag_coefs_emm_ter <- data.frame(Eco = "Terrestrial",
rbindlist(list(cont_vars_res_ter,
cat_emm_vars_res_ter),
fill = TRUE))
mm_mag_coefs_emm_mar <- data.frame(Eco = "Marine",
rbindlist(list(cont_vars_res_mar,
cat_emm_vars_res_mar),
fill = TRUE))
mag_coefs_emm <- rbind(mm_mag_coefs_emm_ter,
mm_mag_coefs_emm_mar)
mag_coefs_emm$level[which(is.na(mag_coefs_emm$level))] <- mag_coefs_emm$var[which(is.na(mag_coefs_emm$level))]
sign_test <- sapply(1:nrow(mag_coefs_emm), function(x){
sign(mag_coefs_emm$Lower.CI[x])==sign(mag_coefs_emm$Upper.CI[x])
})
mag_coefs_emm$sig <- ifelse(sign_test, 1, 0)
mag_coefs_emm$level <- factor(mag_coefs_emm$level, levels = rev(unique(mag_coefs_emm$level)))
mag_coefs_emm$var <- factor(mag_coefs_emm$var, levels = (unique(mag_coefs_emm$var)))
mag_coefs_emm_plot_data <- mag_coefs_emm %>% filter(!var == "Int")
mag_coefs_emm_plot <- ggplot(mag_coefs_emm_plot_data,
aes(x = Estimate, y = level, alpha = sig)) +
geom_vline(xintercept=0, color = "gray",linetype="dashed", lwd=1.5)+
geom_point(size=6, shape = 21, color = "black", fill = "gray")+
geom_errorbar(aes(xmin=Lower.CI,xmax=Upper.CI), width=0, lwd=1.5) +
theme_classic(base_size = 20)+
# coord_cartesian(xlim = c(-1,1))+
xlab("Coefficient\n<<< Match - Mismatch >>>")+ylab("")+
theme(legend.position = "none")+
facet_grid(var~Eco, scales = "free_y", space = "free")
mag_coefs_emm_plot
9 Plot significant effects
9.1 Effects contrast
# get significant effects
mag_coefs <- mag_coefs_cont
vars_sig <- mag_coefs %>%
filter(sig==1,
!var == "Int") %>%
select(var,Eco) %>%
unique %>%
mutate(var = as.character(var))
plot_list_ter <- list()
plot_list_mar <- list()
count_ter = 0
count_mar = 0
for(i in 1:nrow(vars_sig)){
var2go <- as.character(vars_sig[i,"var"])
eco2go <- as.character(vars_sig[i,2])
if(eco2go=="Terrestrial"){
eco2go_ <- "Ter"
}else{
eco2go_ <- "Mar"
}
subdata <- mod_data %>%
filter(Eco==eco2go_) %>%
select(cont_vars, cat_vars, gsub(" ","_",var2go), Eco, mismatch, Param,class,ID) %>%
na.omit()
# is factor or numeric
if(is.numeric(subdata[,gsub(" ","_",var2go)])){
if(eco2go == "Terrestrial"){
p_i <- plot_model(mm_ter, type = "eff", terms = gsub(" ","_",var2go),
line.size = 1.5,
colors = "bw",
title = "") +
theme_classic(base_size = 20)+
ylab("Mismatch")+
xlab(gsub("_"," ",var2go))+
theme(legend.position = "none",
aspect.ratio=1)
count_ter <- count_ter+1
plot_list_ter[[count_ter]] <- p_i
} else {
p_i <- plot_model(mm_mar, type = "eff", terms = gsub(" ","_",var2go),
line.size = 1.5,
colors = "bw",
title = "") +
theme_classic(base_size = 20)+
ylab("Mismatch")+
xlab(gsub("_"," ",var2go))+
theme(legend.position = "none",
aspect.ratio=1)
count_mar <- count_mar+1
plot_list_mar[[count_mar]] <- p_i
}
} else {
subdata_coeffs <- mag_coefs %>%
filter(Eco == eco2go,
var == var2go)%>%
mutate(sig = as.character(sig))
p_i <- ggplot(subdata_coeffs,
aes_string(x="Estimate",y="level", alpha = "sig"))+
geom_vline(xintercept=0, color = "gray",linetype="dashed", lwd=1.5)+
scale_alpha_manual(values = c("1"=1,"0"=0.3))+
geom_point(size=6, shape = 21, color = "black", fill = "gray")+
geom_errorbar(aes(xmin=Lower.CI,xmax=Upper.CI), width=0, lwd=1.5) +
theme_classic(base_size = 20)+
# coord_cartesian(xlim = c(-1,1))+
xlab("Mismatch")+
ylab(gsub("_"," ",var2go))+
theme(legend.position = "none",
aspect.ratio=1)
if(eco2go == "Terrestrial"){
count_ter <- count_ter+1
plot_list_ter[[count_ter]] <- p_i
} else {
count_mar <- count_mar+1
plot_list_mar[[count_mar]] <- p_i
}
}
}Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
ℹ Please use `all_of()` or `any_of()` instead.
# Was:
data %>% select(cont_vars)
# Now:
data %>% select(all_of(cont_vars))
See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
ℹ Please use `all_of()` or `any_of()` instead.
# Was:
data %>% select(cat_vars)
# Now:
data %>% select(all_of(cat_vars))
See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
ℹ Please use tidy evaluation idioms with `aes()`.
ℹ See also `vignette("ggplot2-in-packages")` for more information.plot_list_ter[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
[[6]]
plot_list_mar[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
[[6]]
9.2 Effects pairwise
# get significant effects
mag_coefs <- mag_coefs_cont_pair
vars_sig <- mag_coefs %>%
filter(sig==1,
!var == "Int") %>%
select(var,Eco) %>%
unique %>%
mutate(var = as.character(var))
mag_coefs <- mag_coefs_emm
plot_list_ter <- list()
plot_list_mar <- list()
count_ter = 0
count_mar = 0
for(i in 1:nrow(vars_sig)){
var2go <- as.character(vars_sig[i,"var"])
eco2go <- as.character(vars_sig[i,2])
if(eco2go=="Terrestrial"){
eco2go_ <- "Ter"
}else{
eco2go_ <- "Mar"
}
subdata <- mod_data %>%
filter(Eco==eco2go_) %>%
select(cont_vars, cat_vars, gsub(" ","_",var2go), Eco, mismatch, Param,class,ID) %>%
na.omit()
# is factor or numeric
if(is.numeric(subdata[,gsub(" ","_",var2go)])){
if(eco2go == "Terrestrial"){
mm_2go <- mm_ter
} else {
mm_2go <- mm_mar
}
preds <- ggeffects::ggpredict(mm_2go,terms = gsub(" ","_",var2go))
p_i <- plot_model(mm_2go,
type = "eff",
terms = gsub(" ","_",var2go),
line.size = 1.5,
colors = "bw",
title = "") +
scale_y_continuous(
trans="log",
labels=function(x) round(x,1),
breaks = exp(seq(min(log(preds$conf.low),na.rm = TRUE),
max(log(preds$conf.high),na.rm = TRUE),
length.out = 5)))+
theme_classic(base_size = 20)+
ylab("Mismatch")+
xlab(gsub("_"," ",var2go))+
theme(legend.position = "none",
aspect.ratio=1)
if(eco2go == "Terrestrial"){
count_ter <- count_ter+1
plot_list_ter[[count_ter]] <- p_i
} else {
count_mar <- count_mar+1
plot_list_mar[[count_mar]] <- p_i
}
} else {
subdata_coeffs <- mag_coefs %>%
filter(Eco == eco2go,
var == var2go)%>%
mutate(sig = as.character(sig))
if(eco2go == "Terrestrial"){
cat_vars_cont_pairs_res <- cat_vars_cont_pairs_res_ter
} else {
cat_vars_cont_pairs_res <- cat_vars_cont_pairs_res_mar
}
cat_vars_cont_pairs_res <- cat_vars_cont_pairs_res %>%
filter(var == var2go,
`P-value` <= 0.05)
cat_vars_cont_pairs_res$from <- sapply(cat_vars_cont_pairs_res$level, function(x){
tmp <- strsplit(x," / ")[[1]][1]
})
cat_vars_cont_pairs_res$to <- sapply(cat_vars_cont_pairs_res$level, function(x){
tmp <- strsplit(x," / ")[[1]][2]
})
cat_vars_cont_pairs_res$from <- gsub("(","",cat_vars_cont_pairs_res$from,
fixed = TRUE)
cat_vars_cont_pairs_res$from <- gsub(")","",cat_vars_cont_pairs_res$from,
fixed = TRUE)
cat_vars_cont_pairs_res$to <- gsub("(","",cat_vars_cont_pairs_res$to,
fixed = TRUE)
cat_vars_cont_pairs_res$to <- gsub(")","",cat_vars_cont_pairs_res$to,
fixed = TRUE)
t.test.df <- data.frame(
p.value = as.vector(cat_vars_cont_pairs_res$`P-value`),
right.tip = cat_vars_cont_pairs_res$to,
left.tip = cat_vars_cont_pairs_res$from,
y.bar = max(subdata_coeffs$Upper.CI) + seq(1,nrow(cat_vars_cont_pairs_res))/10
)
t.test.df$label_y = as.character(stars.pval(t.test.df$p.value))
p_i <- ggplot(subdata_coeffs,
aes(x=level,y=Estimate))+
geom_errorbar(aes(ymin=Lower.CI,ymax=Upper.CI), width=0, lwd=1.5) +
geom_point(size=6, shape = 21, color = "black", fill = "gray")+
theme_classic(base_size = 20)+
scale_y_continuous(
trans="log",
labels=function(x) round(x,1),
breaks = exp(seq(min(log(subdata_coeffs$Lower.CI)),
max(log(subdata_coeffs$Upper.CI)),
length.out = 5)))+
ylab("Mismatch")+
xlab(gsub("_"," ",var2go))+
theme(legend.position = "none",
aspect.ratio=1,
axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))+
geom_text_pairwise(data = t.test.df,
aes(xmin = left.tip,
xmax = right.tip,
y = (y.bar),
label = label_y),
na.rm = TRUE,
size = 10, vjust = 0.3)+
expand_limits(y = (max(t.test.df$y.bar)+1))
if(eco2go == "Terrestrial"){
count_ter <- count_ter+1
plot_list_ter[[count_ter]] <- p_i
} else {
count_mar <- count_mar+1
plot_list_mar[[count_mar]] <- p_i
}
}
}Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.plot_list_ter[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
plot_list_mar[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
[[6]]
if(params$data_amount=="all"){
pg <- plot_grid(
plot_grid(plotlist = c(plot_list_ter[1:5],
plot_list_ter[6]),
ncol = 3, labels = LETTERS[1:4], align = "h"),
plot_grid(plotlist = c(plot_list_mar[1:3],
plot_list_mar[4:6]),
ncol = 3, labels = LETTERS[1:6+4], align = "h"),
nrow = 2)
} else {
pg <- plot_grid(
plot_grid(plotlist = c(plot_list_ter[1:3],
plot_list_ter[4]),
ncol = 3, labels = LETTERS[1:4], align = "vh"),
plot_grid(plotlist = c(plot_list_mar[1:3],
plot_list_mar[4:6]),
ncol = 3, labels = LETTERS[1:6+4], align = "vh"),
nrow = 2)
}
pg
# data amount
if(params$data_amount=="all"){
ggsave(plot = pg,
filename = here("Figures/EffectPlots_Mar_Ter_EMM_pairwise_all.svg"),
device = "svg",
width = 4.5*3, height = 6*4,
units = "in")
} else {
ggsave(plot = pg,
filename = here("Figures/EffectPlots_Mar_Ter_EMM_pairwise.svg"),
device = "svg",
width = 4.5*3, height = 6*4,
units = "in")
}10 Further comparisons
# compare flying organisms
pred_values <- predict(mm_ter)
var_to_go <- strsplit(paste(mm_formula_ter, "+", RE_nest),"~|[+]|[/]|[|]")[[1]]
var_to_go <- gsub(" |[(]|[)]","",var_to_go)
var_to_go <- var_to_go[-which(var_to_go=="1")]
subdata_ter <- mod_data %>%
filter(Eco=="Ter") %>%
select(var_to_go) %>%
na.omit
subdata_ter$pred_values <- pred_values
subdata_ter <- subdata_ter %>%
filter(Locomotion_mode=="Flying" | Locomotion_mode=="Sessile")
ggplot(subdata_ter, aes(x = pred_values, y = class))+
geom_boxplot()+
theme_bw()
11 Map study areas
mundi <- rnaturalearth::ne_coastline(scale = 110, returnclass = "sv")
mundi <- crop(mundi, ext(c(-180,180,-60,90)))
# get raster bioshifts shp files study areas
get_raster_bioshifts = "NO"
if(get_raster_bioshifts=="YES"){
my_ext = terra::ext(mundi)
my_crs = crs(mundi)
rast_bioshifts <- terra::rast(my_ext, crs = my_crs, res = 2)
values(rast_bioshifts) <- 0
fgdb <- "C:/Users/brunn/ShadowDrive/CreateGeodatabaseBioShifts/Data/Study_Areas_v1/Study_Areas.gdb"
fc_list <- terra::vector_layers(fgdb)
fc_list <- fc_list[which(fc_list %in% unique(mod_data$ID))]
for(i in 1:length(fc_list)){ cat("\r",i,"from",length(fc_list))
tmp = terra::vect(fgdb, layer = fc_list[i])
tmp = terra::cells(rast_bioshifts,tmp)
tmp_cell = tmp[,2]
tmp_vals = rast_bioshifts[tmp_cell][,1]
rast_bioshifts[tmp_cell] = tmp_vals+1
}
names(rast_bioshifts) <- "SA"
rast_bioshifts[rast_bioshifts==0] <- NA
# data amount
if(params$data_amount=="all"){
writeRaster(rast_bioshifts, here::here("Data/raster_bioshifts_SA_SDMs_all.tif"), overwrite = TRUE)
} else {
writeRaster(rast_bioshifts, here::here("Data/raster_bioshifts_SA_SDMs.tif"), overwrite = TRUE)
}
} else {
# data amount
if(params$data_amount=="all"){
rast_bioshifts <- terra::rast(here::here("Data/raster_bioshifts_SA_SDMs_all.tif"))
} else {
rast_bioshifts <- terra::rast(here::here("Data/raster_bioshifts_SA_SDMs.tif"))
}
rast_bioshifts <- crop(rast_bioshifts, ext(mundi))
}
bio_fig <- ggplot() +
geom_spatraster(data=rast_bioshifts)+
scale_fill_viridis_c(
name = "Number of studies",
alpha=.7,
na.value = "white",
# trans="log",
labels=function(x) round(x,0),
breaks=seq(
as.numeric(global(rast_bioshifts,min,na.rm=TRUE)),
as.numeric(global(rast_bioshifts,max,na.rm=TRUE)),
length.out=4)
)+
geom_spatvector(data=mundi)+
theme_void()bio_fig 
# data amount
if(params$data_amount=="all"){
filetosave <- here("Figures/Number_of_studies_all.png")
} else {
filetosave <- here("Figures/Number_of_studies.png")
}
ggsave(filename = filetosave,
device = "png",
width = 8, height = 3,
units = "in", bg = "white")12 Sensitivity
Compare models with all life form vs models with locomotion mode
cont_vars2 <-
c(
"TSS",
"Study_duration",
"Interannual_var",
"SA_Extent",
"Life_form",
"Connectivity"
)
cat_vars2 <- c(
# "Locomotion_mode",
"Category"
)
## mismatch model magnitude
mm_2_fixed_eff <- paste("mismatch ~ ",
paste(cont_vars2,collapse = " + "),
paste(cat_vars2,collapse = " + "),
sep = " + ")
mm_2_fixed_eff <- gsub(" ~ [+] " , " ~ ", mm_2_fixed_eff)mm_2_formula_ter <- paste(mm_2_fixed_eff, RE_nest, sep = " + ")
# mm_2_formula_mar <- paste(mm_2_fixed_eff, RE_nest, sep = " + ")
# mm_2_formula_ter <- paste(gsub(" [+] Locomotion_mode","",mm_2_fixed_eff), RE_nest, sep = " + ")
mm_2_formula_mar <- paste(gsub(" [+] Connectivity","",mm_2_fixed_eff), RE_nest, sep = " + ")
y <- strsplit(mm_2_formula_ter," ")[[1]][1]Terrestrial: I had to remove Locomotion_mode for the model to fit without singularity issues. We can keep Locomotion_mode removing Life_form works
12.0.1 Terrestrial
var_to_go <- strsplit(mm_2_formula_ter,"~|[+]|[/]|[|]")[[1]]
var_to_go <- gsub(" |[(]|[)]","",var_to_go)
var_to_go <- var_to_go[-which(var_to_go=="1")]
subdata_ter <- mod_data %>%
filter(Eco=="Ter") %>%
select(var_to_go) %>%
na.omit
# fit model
mm_2_ter <- glmmTMB(as.formula(mm_2_formula_ter),
family = Gamma(link = "log"),
data = subdata_ter)
tab_model(mm_2_ter, transform = NULL, string.est = "Estimate")
# plot(effects::allEffects(mm_2_ter))
# make avg model
clust <- makeCluster(detectCores()-3)
clusterExport(clust, c("mm_2_ter","subdata_ter","glmmTMB"))
mm_2_ter_dredge <- dredge(mm_2_ter, cluster = clust)
stopCluster(clust)
# get best models
# remove models that could not calculate AICc (?)
if(any(is.na(mm_2_ter_dredge$AICc))){
pos <- which(is.na(mm_2_ter_dredge$AICc))
mm_2_ter_dredge <- mm_2_ter_dredge[-pos,]
}
mm_2_ter_best_models <- subset(mm_2_ter_dredge, cumsum(weight) <= .95)
# make model avg
mm_2_ter_avg_model <- model.avg(mm_2_ter_best_models, fit = TRUE)
# tab_model(mm_2_ter_avg_model, transform = NULL, string.est = "Estimate")
# marginal R2 model avg
Actual = subdata_ter[,y]
Predicted = predict(mm_2_ter_avg_model, type="response", re.form=NA)
Residuals = Actual - Predicted
mm_2_ter_R2_avg_marginal <- efronRSquared(residual = Residuals,
predicted = Predicted,
statistic = "EfronRSquared")
mm_2_ter_R2_avg_marginal
# conditional R2 model avg
Actual = subdata_ter[,y]
Predicted = predict(mm_2_ter_avg_model, type="response")
Residuals = Actual - Predicted
mm_2_ter_R2_avg_conditional <- efronRSquared(residual = Residuals,
predicted = Predicted,
statistic = "EfronRSquared")
mm_2_ter_R2_avg_conditional12.0.2 Marine
var_to_go <- strsplit(mm_2_formula_mar,"~|[+]|[/]|[|]")[[1]]
var_to_go <- gsub(" |[(]|[)]","",var_to_go)
var_to_go <- var_to_go[-which(var_to_go=="1")]
subdata_mar <- mod_data %>%
filter(Eco=="Mar") %>%
select(var_to_go) %>%
na.omit
# fit model
mm_2_mar <- glmmTMB(as.formula(mm_2_formula_mar),
family = Gamma(link = "log"),
data = subdata_mar)
tab_model(mm_2_mar, transform = NULL, string.est = "Estimate")
# plot(effects::allEffects(mm_2_mar))
# make avg model
clust <- makeCluster(detectCores()-3)
clusterExport(clust, c("mm_2_mar","subdata_mar","glmmTMB"))
mm_2_mar_dredge <- dredge(mm_2_mar, cluster = clust)
stopCluster(clust)
# get best models
# remove models that could not calculate AICc (?)
if(any(is.na(mm_2_mar_dredge$AICc))){
pos <- which(is.na(mm_2_mar_dredge$AICc))
mm_2_mar_dredge <- mm_2_mar_dredge[-pos,]
}
mm_2_mar_best_models <- subset(mm_2_mar_dredge, cumsum(weight) <= .95)
# make model avg
mm_2_mar_avg_model <- model.avg(mm_2_mar_best_models, fit = TRUE)
# tab_model(mm_2_mar_avg_model, transform = NULL, string.est = "Estimate")
# marginal R2 model avg
Actual = subdata_mar[,y]
Predicted = predict(mm_2_mar_avg_model, type="response", re.form=NA)
Residuals = Actual - Predicted
mm_2_mar_R2_avg_marginal <- efronRSquared(residual = Residuals,
predicted = Predicted,
statistic = "EfronRSquared")
mm_2_mar_R2_avg_marginal
# conditional R2 model avg
Actual = subdata_mar[,y]
Predicted = predict(mm_2_mar_avg_model, type="response")
Residuals = Actual - Predicted
mm_2_mar_R2_avg_conditional <- efronRSquared(residual = Residuals,
predicted = Predicted,
statistic = "EfronRSquared")
mm_2_mar_R2_avg_conditional12.0.3 Compare
tab_model(mm_mar,mm_2_mar_avg_model, dv.labels = c("Marine model 1","Marine model 2"))
data.frame(R2_type = c("Marginal","Conditional"),
R2_model_1 = c(mm_mar_R2_avg_marginal,mm_mar_R2_avg_conditional),
R2_model_2 = c(mm_2_mar_R2_avg_marginal,mm_2_mar_R2_avg_conditional))tab_model(mm_ter,mm_2_ter_avg_model, dv.labels = c("Terrestrial model 1","Terrestrial model 2"))
data.frame(R2_type = c("Marginal","Conditional"),
R2_model_1 = c(mm_ter_R2_avg_marginal,mm_ter_R2_avg_conditional),
R2_model_2 = c(mm_2_ter_R2_avg_marginal,mm_2_ter_R2_avg_conditional))