Exploratory analyses

Code

rm(list=ls())
gc()

          used (Mb) gc trigger (Mb) max used (Mb)
Ncells  576905 30.9    1318102 70.4   660394 35.3
Vcells 1055737  8.1    8388608 64.0  1769918 13.6

Code

library(dplyr)
library(tidyverse)
library(tidyr)
library(ggplot2)
library(knitr)
library(pbapply)
library(data.table)
library(GGally)
library(glmmTMB)
library(here)
library(terra)
library(tidyterra)

computer = "personal"
source(here("R/settings.R"))
source(here("R/my_functions.R"))
source(here("R/velocity_functions.R"))

1 Load data

Code

bioshifts <- read.csv(here(Bioshifts_dir,Bioshifts_DB_v3))
bioshifts <- bioshifts_sdms_selection(bioshifts)

bioshifts$new_ID <- paste(bioshifts$ID,
                          bioshifts$sp_name_std,
                          round(bioshifts$Start,0),
                          round(bioshifts$End,0),
                          bioshifts$Type,
                          bioshifts$Param)

# any(duplicated(bioshifts$new_ID))
# test <- bioshifts$new_ID[which(duplicated(bioshifts$new_ID))]
# tocheck <- bioshifts %>% filter(new_ID %in% test)
# write.csv(tocheck,here("tocheck.csv"),row.names = FALSE)
# head(sort(table(tocheck$new_ID)))

Code

sp_bioshifts <- read.csv(here("Data/Output/ExposureData.csv"))

sp_bioshifts$bv.lat.mat <- NA
sp_bioshifts$wv.lat.mat <- NA

North <- (sp_bioshifts$lat.01 > 0 & sp_bioshifts$lat.99 > 0)
South <- (!sp_bioshifts$lat.01 > 0 & !sp_bioshifts$lat.99 > 0)

# LAT + O
pos <- which(sp_bioshifts$Param=="O" & sp_bioshifts$Type == "LAT")
sp_bioshifts$bv.lat.mat[pos] <- sp_bioshifts$bv.lat.median[pos]
sp_bioshifts$wv.lat.mat[pos] <- sp_bioshifts$vel_edge5[pos]
# LAT + LE + North
pos <- which(sp_bioshifts$Param=="LE" & sp_bioshifts$Type == "LAT" & North)
sp_bioshifts$bv.lat.mat[pos] <- sp_bioshifts$bv.lat.95[pos]
sp_bioshifts$wv.lat.mat[pos] <- sp_bioshifts$vel_edge95[pos]
# LAT + LE + South
pos <- which(sp_bioshifts$Param=="LE" & sp_bioshifts$Type == "LAT" & South)
sp_bioshifts$bv.lat.mat[pos] <- sp_bioshifts$bv.lat.05[pos]
sp_bioshifts$wv.lat.mat[pos] <- sp_bioshifts$vel_edge05[pos]
# LAT + TE + North
pos <- which(sp_bioshifts$Param=="TE" & sp_bioshifts$Type == "LAT" & North)
sp_bioshifts$bv.lat.mat[pos] <- sp_bioshifts$bv.lat.05[pos]
sp_bioshifts$wv.lat.mat[pos] <- sp_bioshifts$vel_edge05[pos]
# LAT + TE + South
pos <- which(sp_bioshifts$Param=="TE" & sp_bioshifts$Type == "LAT" & South)
sp_bioshifts$bv.lat.mat[pos] <- sp_bioshifts$bv.lat.95[pos]
sp_bioshifts$wv.lat.mat[pos] <- sp_bioshifts$vel_edge95[pos]

sp_bioshifts$new_ID <- paste(sp_bioshifts$ID,
                             sp_bioshifts$sp_name_std,
                             round(sp_bioshifts$Start,0),
                             round(sp_bioshifts$End,0),
                             sp_bioshifts$Type,
                             sp_bioshifts$Param)

keep_col <- c("new_ID", names(sp_bioshifts)[!names(sp_bioshifts) %in% names(bioshifts)])
sp_bioshifts <- sp_bioshifts %>% dplyr::select(keep_col)
sp_bioshifts <- unique(sp_bioshifts)

bioshifts <- merge(bioshifts, sp_bioshifts, by = "new_ID", all.x = TRUE)

class_select <- sort(table(na.omit(bioshifts)$class),decreasing = TRUE)
class_select <- names(class_select)[which(class_select>20)]
bioshifts <- bioshifts %>% filter(class %in% class_select)

# usable data
# clean_obs <- apply(sp_bioshifts[,24:67], 1, function(x) !any(is.na(x)))
# sp_bioshifts <- sp_bioshifts[which(clean_obs),]

2 Concept figure pieces

This figures will be improved! ## Fig 1a The idea of this figure is to show how the weighted climate velocities are calculated. First we have the climate velocity metrics for the study area, then we overlap and average weight climate velocity by a species suitability map to the same area. Areas of high suitability (where species is more likely to occur will have higher weight). Here we can see that weighted climate velocities can be very different for two species with different distribution ranges in a given study area.

Code

sps <- c("Leuconotopicus_borealis","Gymnorhinus_cyanocephalus")
sps_names <- gsub("_"," ",sps)
# IDtogo <- "A7_P1"
IDtogo <- "A10_P1"

vel_SA <- terra::rast(here("Data/Concept_figure",paste0(IDtogo,"_mat_gVelLat.tif")))

map_sp1 <- terra::rast(here("Data/Concept_figure",paste("sdms",sps[1],".tif")))[[1]]
map_sp2 <- terra::rast(here("Data/Concept_figure",paste("sdms",sps[2],".tif")))[[1]]

map_sp1 <- project(map_sp1, vel_SA)

|---------|---------|---------|---------|
=========================================
                                          

Code

map_sp2 <- project(map_sp2, vel_SA)

|---------|---------|---------|---------|
=========================================
                                          

Code

map_sp1 <- map_sp1/minmax(map_sp1)[2]
map_sp2 <- map_sp2/minmax(map_sp2)[2]


w_vel_sp1 <- vel_SA*map_sp1
w_vel_sp2 <- vel_SA*map_sp2

mydata <- c(vel_SA,w_vel_sp1,w_vel_sp2)
names(mydata) <- c("Velocity SA",paste("Weighted vel.\n",sps_names))
mydata <- terra::values(mydata, nrows=100, na.rm=TRUE, mat=FALSE, dataframe=TRUE)
mydata <- mydata %>% tidyr::pivot_longer(cols = colnames(mydata), names_to = "Velocity", values_to = "value")
mydata$Velocity <- factor(mydata$Velocity, levels = c("Velocity SA",paste("Weighted vel.\n",sps_names)))

p1 <- ggplot(mydata, aes(x = value, y = Velocity)) +
  geom_boxplot(outlier.colour = NA)+
  scale_x_continuous(limits = quantile(mydata$value, c(0.1, 0.9)))+
  labs(x="Climate velocity (km/yr)",y="")+
  theme_classic()


mvel <- velocity_map(vel_SA, 
                     main = paste("Climate velocity\nBioshifts Study area",IDtogo), 
                     ggplot = TRUE)
msp1 <- velocity_map(map_sp1, 
                     main = paste("Suitability map\nSpecies:",sps_names[1]), 
                     ggplot = TRUE)
msp2 <- velocity_map(map_sp2, 
                     main = paste("Suitability map\nSpecies:",sps_names[2]), 
                     ggplot = TRUE)
mwvel_sp1 <- velocity_map(w_vel_sp1, 
                          main = paste("Weighted climate velocity\nBioshifts Study area:",IDtogo,"\nSpecies:",sps_names[1]), 
                          ggplot = TRUE)
mwvel_sp2 <- velocity_map(w_vel_sp2, 
                          main = paste("Weighted climate velocity\nBioshifts Study area:",IDtogo,"\nSpecies:",sps_names[2]), 
                          ggplot = TRUE)

cowplot::plot_grid(
  plotlist = list(
    mvel, msp1, msp2, mwvel_sp1, mwvel_sp2, p1), 
  align = "v", 
  ncol=3, nrow=2)

2.1 Fig 1b

The idea here is to show how the bioclimatic velocities are calculated. We start with time-series of SDM maps for a given species. This time-series refers to projections to multiple years from the beginning to the end a observed range shift. At each grid cell, the bioclimate velocity is calculated using the same formula for climate velocity (Trend/Gradient). We can see here, for a given grid-cell highlighted in blue, bioclimate velocities vary from species to species.

Code

my_point <- terra::vect(data.frame(x=-10000000,y=5000000), geom=c("x", "y"), crs=crs(map_sp1))

ts_sp1 <- terra::rast(here("Data/Concept_figure",paste("sdms",sps[1],".tif")))
ts_sp1 <- terra::extract(ts_sp1,my_point)
ts_sp2 <- terra::rast(here("Data/Concept_figure",paste("sdms",sps[2],".tif")))
ts_sp2 <- terra::extract(ts_sp2,my_point)
ts_sp1 <- data.frame(year = names(ts_sp1)[-1], value = as.numeric(ts_sp1[1,-1]))
ts_sp2 <- data.frame(year = names(ts_sp2)[-1], value = as.numeric(ts_sp2[1,-1]))

ts1 <- ggplot(ts_sp1, aes(x = as.numeric(year), y = value/max(value))) +
  geom_point()+
  labs(x="Year",y="SDM suitability",title = sps_names[1])+
  theme_classic()+
  geom_smooth(color = "black", method = "lm")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

ts2 <- ggplot(ts_sp2, aes(x = as.numeric(year), y = value/max(value))) +
  geom_point()+
  labs(x="Year",y="SDM suitability",title = sps_names[2])+
  theme_classic()+
  geom_smooth(color = "black", method = "lm")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

bvel_sp1 <- terra::rast(here("Data/Concept_figure",paste("sdms_vel",sps[1],".tif")))
bvel_sp2 <- terra::rast(here("Data/Concept_figure",paste("sdms_vel",sps[2],".tif")))

mydata <- c(bvel_sp1,bvel_sp2)
names(mydata) <- paste("Bioclimatic vel.\n",sps_names)
mydata <- terra::values(mydata, nrows=100, na.rm=TRUE, mat=FALSE, dataframe=TRUE)
mydata <- mydata %>% tidyr::pivot_longer(cols = colnames(mydata), names_to = "bVel", values_to = "value")
mydata$bVel <- factor(mydata$bVel, levels = paste("Bioclimatic vel.\n",sps_names))

p1 <- ggplot(mydata, aes(x = value, y = bVel)) +
  geom_boxplot(outlier.colour = NA)+
  scale_x_continuous(limits = quantile(mydata$value, c(0.1, 0.9)))+
  labs(x="Bioclimatic velocity (km/yr)",y="")+
  theme_classic()

mbvel_sp1 <- velocity_map(map_sp1, 
                          main = paste("Suitability map\nSpecies:",sps_names[1]), 
                          ggplot = TRUE) +
  geom_spatvector(data = my_point, color = "blue")
mbvel_sp2 <- velocity_map(map_sp2, 
                          main = paste("Suitability map\nSpecies:",sps_names[2]), 
                          ggplot = TRUE) +
  geom_spatvector(data = my_point, color = "blue")

cowplot::plot_grid(
  plotlist = list(
    mbvel_sp1, mbvel_sp2, NA, ts1, ts2, p1), 
  align = "v", 
  ncol=3, nrow=2)

3 Explore relationships

3.1 N shifts per parameter

Code

# tmp <- data.frame(table(sp_bioshifts$Param))
# 
# ggplot(tmp, aes(x = Var1, y = Freq))+
#     geom_col()+
#     theme_classic()+
#     geom_text(aes(label = Freq), vjust = -0.5)+
#     xlab("N shifts")
# 
# tmp <- data.frame(table(sp_bioshifts$class))
# 
# ggplot(tmp, aes(x = Var1, y = Freq))+
#     geom_col()+
#     theme_classic()+
#     geom_text(aes(label = Freq), vjust = -0.5)+
#     xlab("N shifts")+
#   coord_flip()

4 How related are observed shifts and climate exposure metrics?

Code

mydata <- bioshifts %>%
  mutate(`Bioclimatic vel.` = bv.lat.mat,
         `Weighted vel.` = wv.lat.mat,
         `Biotic vel.` = ShiftRate) %>%
  select(Eco, Param, `Bioclimatic vel.`, `Weighted vel.`, `Biotic vel.`) %>%
  pivot_longer(cols = -c(Eco, Param), names_to = "Metric", values_to = "Velocity") %>%
  na.omit()

plot_cat <- expand.grid(Param=c("LE","O","TE"),
                        realms=c("Ter","Mar"))

plotlist <- lapply(1:nrow(plot_cat), function(i){
  ggplot(mydata%>%
           filter(Eco == plot_cat$realms[i],
                  Param == plot_cat$Param[i]), 
         aes(x = Velocity, y = Metric))+
    ggtitle(paste(plot_cat$realm[i],plot_cat$Param[i],sep = " - "))+
    geom_boxplot(outlier.colour = NA)+
    theme_classic()+
    scale_x_continuous(limits = quantile(mydata$Velocity, c(0.1, 0.9), 
                                         na.rm = TRUE))+
    theme(aspect.ratio=1,
          plot.title = element_text(hjust = 0.5))
})

ggpubr::ggarrange(
  plotlist = plotlist, 
  align = "v", 
  ncol=3, nrow=2,common.legend = TRUE)

4.1 Bioclimatic velocity

Code

plot_cat <- expand.grid(Param=c("LE","O","TE"),
                        realms=c("Ter","Mar"))

plotlist <- lapply(1:nrow(plot_cat), function(i){
  ggplot(bioshifts %>%
           filter(Eco == plot_cat$realms[i],
                  Param == plot_cat$Param[i]) %>% 
           select(bv.lat.mat,ShiftRate,Eco,Duration,Param,class) %>%
           na.omit, 
         aes(x = bv.lat.mat, y = ShiftRate, size = Duration, color = class))+
    geom_point(alpha = .5)+
    labs(x="Bioclimatic velocity (km/yr SDM suitability)", 
         y="Observed range shift velocity (km/yr)",
         title = paste(plot_cat$realm[i],plot_cat$Param[i],sep = " - "))+
    theme_classic() +
    theme(aspect.ratio=1,
          plot.title = element_text(hjust = 0.5)) +
    tune::coord_obs_pred()+
    # geom_smooth(method = "lm")+
    geom_abline(intercept = 0, slope = 1, linetype = "dashed")
})

ggpubr::ggarrange(
  plotlist = plotlist, 
  align = "v", 
  ncol=3, nrow=2,common.legend = TRUE)

Compare across taxonomic classes

Code

plot_cat <- expand.grid(Param=c("LE","O","TE"),
                        realms=c("Ter"),
                        class=unique(bioshifts$class))

plotlist <- lapply(1:nrow(plot_cat), function(i){
  datatogo <- bioshifts %>%
           filter(Eco == plot_cat$realms[i],
                  Param == plot_cat$Param[i],
                  class == plot_cat$class[i]) %>% 
           select(bv.lat.mat,ShiftRate,Eco,Duration,Param,class) %>%
           na.omit
  if(nrow(datatogo)>10){
    ggplot(datatogo, 
         aes(x = bv.lat.mat, y = ShiftRate, size = Duration, color = class))+
    geom_point(alpha = .5)+
    labs(x="Bioclimatic velocity (km/yr SDM suitability)", 
         y="Observed range shift velocity (km/yr)",
         title = paste(plot_cat$realm[i],plot_cat$Param[i],plot_cat$class[i],
                       sep = " - "))+
    theme_classic() +
    theme(aspect.ratio=1,
          plot.title = element_text(hjust = 0.5)) +
    tune::coord_obs_pred()+
    # geom_smooth(method = "lm")+
    geom_abline(intercept = 0, slope = 1, linetype = "dashed")
  }
})
keep <- sapply(plotlist, is.null)
plotlist <- plotlist[!keep]

ggpubr::ggarrange(
  plotlist = plotlist, 
  align = "v", common.legend = TRUE)

Code

plot_cat <- expand.grid(Param=c("LE","O","TE"),
                        realms=c("Mar"),
                        class=unique(bioshifts$class))

plotlist <- lapply(1:nrow(plot_cat), function(i){
  datatogo <- bioshifts %>%
           filter(Eco == plot_cat$realms[i],
                  Param == plot_cat$Param[i],
                  class == plot_cat$class[i]) %>% 
           select(bv.lat.mat,ShiftRate,Eco,Duration,Param,class) %>%
           na.omit
  if(nrow(datatogo)>10){
    ggplot(datatogo, 
         aes(x = bv.lat.mat, y = ShiftRate, size = Duration, color = class))+
    geom_point(alpha = .5)+
    labs(x="Bioclimatic velocity (km/yr SDM suitability)", 
         y="Observed range shift velocity (km/yr)",
         title = paste(plot_cat$realm[i],plot_cat$Param[i],plot_cat$class[i],
                       sep = " - "))+
    theme_classic() +
    theme(aspect.ratio=1,
          plot.title = element_text(hjust = 0.5)) +
    tune::coord_obs_pred()+
    # geom_smooth(method = "lm")+
    geom_abline(intercept = 0, slope = 1, linetype = "dashed")
  }
})
keep <- sapply(plotlist, is.null)
plotlist <- plotlist[!keep]

ggpubr::ggarrange(
  plotlist = plotlist, 
  align = "v", common.legend = TRUE)

Check if duration can explain the decoupling between observed shifts and bioclimatic velocities

Code

ggplot(bioshifts %>%
         filter(Eco == "Ter") %>%
         select(bv.lat.mat,ShiftRate,Duration,Param) %>%
         na.omit, 
       aes(x = Duration, y = bv.lat.mat-ShiftRate))+
  geom_point(alpha = .5)+
  labs(x="Shift duration", 
       y="Difference Bioclimatic velocity - Observed range shift velocity (km/yr)",
       title = "Terrestrial")+
  theme_classic() +
  geom_smooth(method = "lm")+
  geom_hline(yintercept = 0, linetype = "dashed")+
  facet_wrap(Param~., ncol = 3, scales = "free_y")

Code

ggplot(bioshifts %>%
         filter(Eco == "Mar") %>%
         select(bv.lat.mat,ShiftRate,Duration,Param) %>%
         na.omit, 
       aes(x = Duration, y = bv.lat.mat-ShiftRate))+
  geom_point(alpha = .5)+
  labs(x="Shift duration", 
       y="Difference Bioclimatic velocity - Observed range shift velocity (km/yr)",
       title = "Marine")+
  theme_classic() +
  geom_smooth(method = "lm")+
  geom_hline(yintercept = 0, linetype = "dashed")+
  facet_wrap(Param~., ncol = 3, scales = "free_y")

Check if SDM performance can explain the decoupling between observed shifts and bioclimatic velocities

Code

ggplot(bioshifts %>%
         filter(Eco == "Ter") %>%
         select(bv.lat.mat,ShiftRate,sdm_sensitivity_mean,Param) %>%
         na.omit, 
       aes(x = sdm_sensitivity_mean, y = bv.lat.mat-ShiftRate))+
  geom_point(alpha = .5)+
  labs(x="SDM sensitivity", 
       y="Difference Bioclimatic velocity - Observed range shift velocity (km/yr)",
       title = "Terrestrial")+
  theme_classic() +
  geom_smooth(method = "lm")+
  geom_hline(yintercept = 0, linetype = "dashed")+
  facet_wrap(Param~., ncol = 3, scales = "free_y")

Code

ggplot(bioshifts %>%
         filter(Eco == "Mar") %>%
         select(bv.lat.mat,ShiftRate,sdm_sensitivity_mean,Param) %>%
         na.omit, 
       aes(x = sdm_sensitivity_mean, y = bv.lat.mat-ShiftRate))+
  geom_point(alpha = .5)+
  labs(x="SDM sensitivity", 
       y="Difference Bioclimatic velocity - Observed range shift velocity (km/yr)",
       title = "Marine")+
  theme_classic() +
  geom_smooth(method = "lm")+
  geom_hline(yintercept = 0, linetype = "dashed")+
  facet_wrap(Param~., ncol = 3, scales = "free_y")

4.2 Climate velocity

Code

plot_cat <- expand.grid(Param=c("LE","O","TE"),
                        realms=c("Ter","Mar"))

plotlist <- lapply(1:nrow(plot_cat), function(i){
  ggplot(bioshifts %>%
           filter(Eco == plot_cat$realms[i],
                  Param == plot_cat$Param[i]) %>% 
           select(vel_mat_lat,ShiftRate,Eco,Duration,Param,class) %>%
           na.omit, 
         aes(x = vel_mat_lat, y = ShiftRate, size = Duration, color = class))+
    geom_point(alpha = .5)+
    labs(x="Climatic velocity (km/yr Temperature)", 
         y="Observed range shift velocity (km/yr)",
         title = paste(plot_cat$realm[i],plot_cat$Param[i],sep = " - "))+
    theme_classic() +
    theme(aspect.ratio=1,
          plot.title = element_text(hjust = 0.5)) +
    tune::coord_obs_pred()+
    # geom_smooth(method = "lm")+
    geom_abline(intercept = 0, slope = 1, linetype = "dashed")
})

ggpubr::ggarrange(
  plotlist = plotlist, 
  align = "v", 
  ncol=3, nrow=2,common.legend = TRUE)

4.3 Weighted climate velocity

Code

plot_cat <- expand.grid(Param=c("LE","O","TE"),
                        realms=c("Ter","Mar"))

plotlist <- lapply(1:nrow(plot_cat), function(i){
  ggplot(bioshifts %>%
           filter(Eco == plot_cat$realms[i],
                  Param == plot_cat$Param[i]) %>% 
           select(wv.lat.mat,ShiftRate,Eco,Duration,Param,class) %>%
           na.omit, 
         aes(x = wv.lat.mat, y = ShiftRate, size = Duration, color = class))+
    geom_point(alpha = .5)+
    labs(x="Weighted climatic velocity (km/yr)", 
         y="Observed range shift velocity (km/yr)",
         title = paste(plot_cat$realm[i],plot_cat$Param[i],sep = " - "))+
    theme_classic() +
    theme(aspect.ratio=1,
          plot.title = element_text(hjust = 0.5)) +
    tune::coord_obs_pred()+
    # geom_smooth(method = "lm")+
    geom_abline(intercept = 0, slope = 1, linetype = "dashed")
})

ggpubr::ggarrange(
  plotlist = plotlist, 
  align = "v", 
  ncol=3, nrow=2,common.legend = TRUE)

Compare across taxonomic classes

Code

ggplot(bioshifts %>%
         filter(Eco == "Mar") %>% 
         select(wv.lat.mat,ShiftRate,Eco,Duration,Param,class) %>%
         na.omit, 
       aes(x = wv.lat.mat, y = ShiftRate, size = Duration, color = class))+
  geom_point(alpha = .5)+
  labs(x="Weighted climate velocity (km/yr)", 
       y="Observed range shift velocity (km/yr)",
       title = "Marine")+
  theme_classic() +
  # geom_smooth(method = "lm")+
  geom_abline(intercept = 0, slope = 1, linetype = "dashed")+
  facet_wrap(class~Param, ncol = 3, scales = "free")

Code

ggplot(bioshifts %>%
         filter(Eco == "Ter") %>% 
         select(wv.lat.mat,ShiftRate,Eco,Duration,Param,class) %>%
         na.omit, 
       aes(x = wv.lat.mat, y = ShiftRate, size = Duration, color = class))+
  geom_point(alpha = .5)+
  labs(x="Weighted climate velocity (km/yr)", 
       y="Observed range shift velocity (km/yr)",
       title = "Terrestrial")+
  theme_classic() +
  # geom_smooth(method = "lm")+
  geom_abline(intercept = 0, slope = 1, linetype = "dashed")+
  facet_wrap(class~Param, ncol = 3, scales = "free")

4.4 Weighted climate velocity vs Bioclimatic velocity

Code

plot_cat <- expand.grid(Param=c("LE","O","TE"),
                        realms=c("Ter","Mar"))

plotlist <- lapply(1:nrow(plot_cat), function(i){
  ggplot(bioshifts %>%
           filter(Eco == plot_cat$realms[i],
                  Param == plot_cat$Param[i]) %>% 
           select(wv.lat.mat,bv.lat.mat,Eco,Duration,Param,class) %>%
           na.omit, 
         aes(x = wv.lat.mat, y = bv.lat.mat, size = Duration, color = class))+
    geom_point(alpha = .5)+
    labs(x="Weighted climatic velocity (km/yr)", 
         y="Bioclimatic velocity (km/yr SDM suitability)", 
         title = paste(plot_cat$realm[i],plot_cat$Param[i],sep = " - "))+
    theme_classic() +
    theme(aspect.ratio=1,
          plot.title = element_text(hjust = 0.5)) +
    tune::coord_obs_pred()+
    # geom_smooth(method = "lm")+
    geom_abline(intercept = 0, slope = 1, linetype = "dashed")
})

ggpubr::ggarrange(
  plotlist = plotlist, 
  align = "v", 
  ncol=3, nrow=2,common.legend = TRUE)

5 Difference observed shift and bioclimatic velocity

Code

ggplot(bioshifts %>%
         filter(Eco == "Mar") %>% 
         select(bv.lat.mat,ShiftRate,Eco,Duration,Param,class) %>%
         mutate(r.obs.bv = ShiftRate-bv.lat.mat) %>%
         filter(r.obs.bv > quantile(r.obs.bv,na.rm = TRUE,probs = .05) & r.obs.bv < quantile(r.obs.bv,na.rm = TRUE,probs = .95)) %>% 
         na.omit, 
       aes(x = r.obs.bv, size = Duration, fill = class))+
  geom_density(alpha = .5)+
  labs(x="Difference observed shift - bioclimatic velocity\n<< Faster bioclimatic velocity           Faster observed shift velocity >>",
       title = "Marine")+
  theme_classic() +
  geom_vline(xintercept = 0, linetype = "dashed")+
  facet_wrap(.~Param, nrow = 1, scales = "free")

Code

ggplot(bioshifts %>%
         filter(Eco == "Ter") %>% 
         select(bv.lat.mat,ShiftRate,Eco,Duration,Param,class) %>%
         mutate(r.obs.bv = ShiftRate-bv.lat.mat) %>%
         filter(r.obs.bv > quantile(r.obs.bv,na.rm = TRUE,probs = .05) & r.obs.bv < quantile(r.obs.bv,na.rm = TRUE,probs = .95)) %>% 
         na.omit, 
       aes(x = r.obs.bv, size = Duration, fill = class))+
  geom_density(alpha = .5)+
  labs(x="Difference observed shift - bioclimatic velocity\n<< Faster bioclimatic velocity           Faster observed shift velocity >>",
       title = "Terrestrial")+
  theme_classic() +
  # geom_smooth(method = "lm")+
  geom_vline(xintercept = 0, linetype = "dashed")+
  facet_wrap(.~Param, nrow = 1, scales = "free")

6 Models

Run separate models for marine and terrestrials

Code

models_to_go <- c("ShiftRate ~ bv.lat.mat * Param + 
                  N_periodes + LatExtentk + Grain_size + Obs_type + Uncertainty_distribution + 
                  (bv.lat.mat * Param | class) + 
                  (1 | ID)",
                  
                  "ShiftRate ~ wv.lat.mat * Param + 
                  N_periodes + LatExtentk + Grain_size + Obs_type + Uncertainty_distribution + 
                  (wv.lat.mat * Param | class) + 
                  (1 | ID)")
models_to_go <- expand.grid(models_to_go,Eco = c("Ter","Mar"))


vel_models <- lapply(1:nrow(models_to_go), function(i){
  glmmTMB(as.formula(as.character(models_to_go[i,1])), 
          weights = sdm_sensitivity_mean,
          data = bioshifts %>% 
            select(ShiftRate,sdm_sensitivity_mean,bv.lat.mat,wv.lat.mat,Param,N_periodes,LatExtentk,Grain_size,Obs_type,Uncertainty_distribution,class,Eco,ID) %>% 
            filter(Eco==models_to_go[i,2]) %>% 
            na.omit())
})

vel_models_summary <- lapply(vel_models, function(x){
  tmp <- summary(x)
  tmp <- tmp$coefficients
  tmp <- round(tmp$cond,3)
  tmp <- data.frame(variable = rownames(tmp), tmp)
  tmp
})


models_to_go[1,]

                                                                                                                                                                                                               Var1
1 ShiftRate ~ bv.lat.mat * Param + \n                  N_periodes + LatExtentk + Grain_size + Obs_type + Uncertainty_distribution + \n                  (bv.lat.mat * Param | class) + \n                  (1 | ID)
  Eco
1 Ter

Code

nice_table(vel_models_summary[[1]])

Code

models_to_go[2,]

                                                                                                                                                                                                               Var1
2 ShiftRate ~ wv.lat.mat * Param + \n                  N_periodes + LatExtentk + Grain_size + Obs_type + Uncertainty_distribution + \n                  (wv.lat.mat * Param | class) + \n                  (1 | ID)
  Eco
2 Ter

Code

nice_table(vel_models_summary[[2]])

Code

models_to_go[3,]

                                                                                                                                                                                                               Var1
3 ShiftRate ~ bv.lat.mat * Param + \n                  N_periodes + LatExtentk + Grain_size + Obs_type + Uncertainty_distribution + \n                  (bv.lat.mat * Param | class) + \n                  (1 | ID)
  Eco
3 Mar

Code

nice_table(vel_models_summary[[3]])

Code

models_to_go[4,]

                                                                                                                                                                                                               Var1
4 ShiftRate ~ wv.lat.mat * Param + \n                  N_periodes + LatExtentk + Grain_size + Obs_type + Uncertainty_distribution + \n                  (wv.lat.mat * Param | class) + \n                  (1 | ID)
  Eco
4 Mar

Code

nice_table(vel_models_summary[[4]])