Exploratory analyses

rm(list=ls())
gc()

          used (Mb) gc trigger (Mb) max used (Mb)
Ncells  580680 31.1    1307771 69.9   660394 35.3
Vcells 1073040  8.2    8388608 64.0  1769918 13.6

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)
library(rcompanion)

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

1 Load data

bioshifts <- read.csv(here(Bioshifts_dir,Bioshifts_DB_v3))
# remove shifts coded as non-significant
bioshifts <- bioshifts %>% filter(!Significance == "N")

# subset usefull columns from bioshifts
bioshifts <- bioshifts_sdms_selection(bioshifts)
# remove shifts == 0
# bioshifts <- bioshifts %>% filter(abs(ShiftRate) > 0)

# create ID for merging
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)))

# N species
length(unique(bioshifts$sp_name_std))

[1] 4512

bioshifts %>% 
  group_by(Eco) %>%
  summarise(N=length(unique(sp_name_std)))

# A tibble: 2 × 2
  Eco       N
  <chr> <int>
1 Mar     707
2 Ter    3805

# N species with GBIF data
n_occ <- read.csv(here("Data/N_OCC.csv"))
n_occ$scientificName <- gsub(" ","_",n_occ$scientificName)
n_occ <- n_occ %>% filter(N_OCC >= 30)

sps <- unique(bioshifts$sp_name_std[which(bioshifts$sp_name_std %in% n_occ$scientificName)])
length(sps)

[1] 4255

bioshifts %>%
  filter(sp_name_std %in% sps) %>%
  group_by(Eco) %>%
  summarise(N=length(unique(sp_name_std)))

# A tibble: 2 × 2
  Eco       N
  <chr> <int>
1 Mar     611
2 Ter    3644

Add exposure metrics

sp_bioshifts <- read.csv(here("Data/Output/ExposureData.csv"))
# names(sp_bioshifts)
sp_bioshifts %>% group_by(Eco) %>% summarise(N=length(unique(sp_name_std)))

# A tibble: 2 × 2
  Eco             N
  <chr>       <int>
1 Marine        766
2 Terrestrial  3871

# Variables to get
# bioclimatic velocity (slope from lm for latitude)
sp_bioshifts$bv.mat <- NA 
# bioclimatic velocity (std error from lm for latitude)
sp_bioshifts$bv.mat.error <- NA 
# bioclimatic velocity (slope from lm for latitude). 
# Use the centroid shift when study estimated shifts of one edge only. We assume the suty area is located at the edge position.
sp_bioshifts$bv.mat2 <- NA 
# bioclimatic angle (angle of the vector from the magnitudes lat and long)
sp_bioshifts$bv.mat.angle <- NA 

# edge velocities
# climate change velocity at the edge
sp_bioshifts$ev.mat <- NA 
# climate change velocity at the edge. Use the centroid shift when study estimated shifts of one edge only
sp_bioshifts$ev.mat2 <- NA 

# Change shift direction for studies located at the south hemisphere
# Find studies located at the south hemisphere
North <- (sp_bioshifts$SA_y_min > 0 & sp_bioshifts$SA_y_max > 0)
South <- (!sp_bioshifts$SA_y_min > 0 & !sp_bioshifts$SA_y_max > 0)
sp_bioshifts$Hemisphere <- ifelse(North,"North hemisphere","South hemisphere")

# Populate bioclimatic velocities
# O
pos <- which(sp_bioshifts$Param=="O")
sp_bioshifts$bv.mat[pos] <- sp_bioshifts$centroid_shift_lat_lm[pos] 
sp_bioshifts$bv.mat.error[pos] <- sp_bioshifts$centroid_shift_lat_lm_error[pos] 
sp_bioshifts$bv.mat.angle[pos] <- sp_bioshifts$centroid_shift_xy_angle[pos] 
sp_bioshifts$ev.mat[pos] <- sp_bioshifts$vel_edge5[pos] 
# LE + North
pos <- which(sp_bioshifts$Param=="LE" & North)
sp_bioshifts$bv.mat[pos] <- sp_bioshifts$edge_shift_lat_lm_0.75[pos]
sp_bioshifts$bv.mat.error[pos] <- sp_bioshifts$edge_shift_lat_lm_error_0.75[pos] 
sp_bioshifts$bv.mat.angle[pos] <- sp_bioshifts$edge_shift_xy_angle_0.75[pos] 
sp_bioshifts$ev.mat[pos] <- sp_bioshifts$vel_edge75[pos] 
# TE + North
pos <- which(sp_bioshifts$Param=="TE" & North)
sp_bioshifts$bv.mat[pos] <- sp_bioshifts$edge_shift_lat_lm_0.25[pos]
sp_bioshifts$bv.mat.error[pos] <- sp_bioshifts$edge_shift_lat_lm_error_0.25[pos] 
sp_bioshifts$bv.mat.angle[pos] <- sp_bioshifts$edge_shift_xy_angle_0.25[pos] 
sp_bioshifts$ev.mat[pos] <- sp_bioshifts$vel_edge25[pos] 
# LE + South
pos <- which(sp_bioshifts$Param=="LE" & South)
sp_bioshifts$bv.mat[pos] <- sp_bioshifts$edge_shift_lat_lm_0.25[pos]
sp_bioshifts$bv.mat.error[pos] <- sp_bioshifts$edge_shift_lat_lm_error_0.25[pos] 
sp_bioshifts$bv.mat.angle[pos] <- sp_bioshifts$edge_shift_xy_angle_0.25[pos] 
sp_bioshifts$ev.mat[pos] <- sp_bioshifts$vel_edge25[pos] 
# TE + South
pos <- which(sp_bioshifts$Param=="TE" & South)
sp_bioshifts$bv.mat[pos] <- sp_bioshifts$edge_shift_lat_lm_0.75[pos]
sp_bioshifts$bv.mat.error[pos] <- sp_bioshifts$edge_shift_lat_lm_error_0.75[pos] 
sp_bioshifts$bv.mat.angle[pos] <- sp_bioshifts$edge_shift_xy_angle_0.75[pos] 
sp_bioshifts$ev.mat[pos] <- sp_bioshifts$vel_edge75[pos] 

# Change the sign of velocities in the South hemisphere (positive values mean edges shift towards the poles)
pos <- which(South)
sp_bioshifts$bv.mat[pos] <- sp_bioshifts$bv.mat[pos] * -1

# Use centroid shift for studies with only one edge
# if study estimated shifts of one edge only, we assume the SA is located at the edge of species range and thus we can use the SA centroid shift directly.
one_edge_study <- bioshifts %>%
  filter(!Param=="O") %>%
  group_by(ID) %>%
  summarise(N = length(unique(Param))) %>%
  filter(N == 1) %>%
  select(ID)

pos <- which(sp_bioshifts$ID %in% one_edge_study$ID)
sp_bioshifts$bv.mat2 <- sp_bioshifts$bv.mat
sp_bioshifts$bv.mat2[pos] <- sp_bioshifts$centroid_shift_lat_lm[pos] 
sp_bioshifts$ev.mat2 <- sp_bioshifts$bv.mat
sp_bioshifts$ev.mat2[pos] <- sp_bioshifts$vel_edge5[pos] 

# shift angle in the south hemisphere
# pos <- which(South)
# sp_bioshifts$bv.mat.angle[pos] <- (sp_bioshifts$bv.mat.angle[pos] + 180) %% 360


# sp_bioshifts$Param[pos] 
# View(sp_bioshifts[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)

sp_bioshifts <- sp_bioshifts %>%
  select(new_ID, Hemisphere, bv.mat, bv.mat.error, bv.mat2, bv.mat.angle,
         ev.mat, ev.mat2,
         vel_mat_lat_median, vel_mat_lat_median_25km, 
         sensitivity, validation, calibration,
         impeded_prop, diffuse_prop, intensified_prop, channelized_prop,
         ecto.endo, Mobility, Locomotion_mode)

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)

# add methods from bioshifts to the data
bioshifts <- merge(bioshifts, sp_bioshifts, by = "new_ID", all.x = TRUE)

# shifts are in meters
# convert to km
bioshifts$bv.mat <- bioshifts$bv.mat/1000
bioshifts$bv.mat2 <- bioshifts$bv.mat2/1000
bioshifts$bv.mat.error <- bioshifts$bv.mat.error/1000

# calculate signal to noise ration
bioshifts$signoise <- abs(bioshifts$bv.mat) / bioshifts$bv.mat.error
# bioshifts$signoise <- sqrt(bioshifts$bv.mat) / sqrt(bioshifts$bv.mat.error)

# calculate mismatch between observed and predicted
bioshifts$mismatch <- bioshifts$ShiftRate - bioshifts$bv.mat
bioshifts$mismatch_percent <- (bioshifts$mismatch/bioshifts$ShiftRate)*100

# get mismatch binary
bioshifts$mismatch_bin <- ifelse(bioshifts$mismatch < 0, 
                                 "observed shift > predicted shift", 
                                 "observed shift < predicted shift") 

# get mismatch direction
bioshifts$match_direction <- "mismatch" 
bioshifts$match_direction[which(bioshifts$ShiftRate > 0 & bioshifts$bv.mat  > 0 | 
                                  bioshifts$ShiftRate < 0 & bioshifts$bv.mat  < 0)] <- "match"

# get mismatch classification
bioshifts$match_direction_class <- NA
bioshifts$match_direction_class[which(bioshifts$ShiftRate > 0 & bioshifts$bv.mat  > 0)] <- "pospos"
bioshifts$match_direction_class[which(bioshifts$ShiftRate > 0 & bioshifts$bv.mat  < 0)] <- "posneg"
bioshifts$match_direction_class[which(bioshifts$ShiftRate < 0 & bioshifts$bv.mat  > 0)] <- "negpos"
bioshifts$match_direction_class[which(bioshifts$ShiftRate < 0 & bioshifts$bv.mat  < 0)] <- "negneg"

# classify grain
bioshifts$Grain_size <- factor(bioshifts$Grain_size, levels = c("small","moderate","large","very_large"))

# classify category
bioshifts$Category <- factor(bioshifts$Category, levels = c("TimeSeries","CensusPeriods","Survey-Resurvey"))

N species predict to no shift

# bioshifts %>% 
#   filter(!is.na(bv.mat)) %>% 
#   group_by(Eco, Param) %>%
#   summarise(Shift = length(which(edge_sd > 0)),
#             NoShift = length(which(edge_sd == 0)),
#             Proportion = round(Shift/NoShift,2),
#             Total = length(Eco))

Conversions and other metrics calculation

# remove species with SA outside edge (these species don't shift their ranges)
# bioshifts <- bioshifts %>%
#   filter(edge_sd > 0 | is.na(edge_sd))

# select classes with more than X observations per edge
N_obs_class = 5
class_param_select <- bioshifts %>%
  filter(!is.na(bv.mat)) %>%
  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)

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

# exclude observations without predicted shift
bioshifts <- bioshifts %>% filter(!is.na(bv.mat))

# N species
length(unique(bioshifts$sp_name_std))

[1] 922

bioshifts %>% 
  group_by(Param) %>% 
  summarise(N=length(unique(sp_name_std)))

# A tibble: 3 × 2
  Param     N
  <chr> <int>
1 LE      352
2 O       670
3 TE       93

bioshifts %>% 
  group_by(Eco) %>% 
  summarise(N=length(unique(sp_name_std)))

# A tibble: 2 × 2
  Eco       N
  <chr> <int>
1 Mar     530
2 Ter     392

# N range shifts
length(bioshifts$sp_name_std)

[1] 1543

bioshifts %>% 
  group_by(Param) %>% 
  summarise(N=length(sp_name_std))

# A tibble: 3 × 2
  Param     N
  <chr> <int>
1 LE      393
2 O      1051
3 TE       99

bioshifts %>% 
  group_by(Eco) %>%
  summarise(N=length(sp_name_std))

# A tibble: 2 × 2
  Eco       N
  <chr> <int>
1 Mar     951
2 Ter     592

2 SDM evaluation results

sdms_CV <- bioshifts %>%
  select("sensitivity",
         "calibration","validation",
         "Eco","sp_name_std","class") %>%
  group_by(sp_name_std,Eco,class) %>%
  summarise(sensitivity = mean(sensitivity,na.rm=TRUE),
            calibration = mean(calibration,na.rm=TRUE),
            validation = mean(validation,na.rm=TRUE))

`summarise()` has grouped output by 'sp_name_std', 'Eco'. You can override
using the `.groups` argument.

ggplot(sdms_CV, aes(x = validation))+
  geom_histogram(fill="white",color="black")+
  facet_wrap(.~Eco, scales = "free")+
  labs(x="TSS validation",y="")+
  geom_vline(xintercept = .5,color='red', linewidth = 1, alpha = .5)+
  geom_vline(xintercept = .7,color='blue', linewidth = 1, alpha = .5)+
  theme_classic()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: Removed 61 rows containing non-finite outside the scale range
(`stat_bin()`).

ggplot(sdms_CV, aes(x = calibration, y = class))+
  geom_boxplot(outlier.alpha = 0)+
  facet_wrap(.~Eco, scales = "free")+
  labs(x="TSS calibration",y="")+
  geom_vline(xintercept = .5,color='red', linewidth = 1, alpha = .5)+
  geom_vline(xintercept = .7,color='blue', linewidth = 1, alpha = .5)+
  theme_classic()

Warning: Removed 1 row containing non-finite outside the scale range
(`stat_boxplot()`).

ggplot(sdms_CV, aes(x = validation, y = class))+
  geom_boxplot(outlier.alpha = 0)+
  facet_wrap(.~Eco, scales = "free")+
  labs(x="TSS validation",y="")+
  geom_vline(xintercept = .5,color='red', linewidth = 1, alpha = .5)+
  geom_vline(xintercept = .7,color='blue', linewidth = 1, alpha = .5)+
  theme_classic()

Warning: Removed 61 rows containing non-finite outside the scale range
(`stat_boxplot()`).

ggplot(sdms_CV, aes(x = sensitivity, y = class))+
  geom_boxplot(outlier.alpha = 0)+
  facet_wrap(.~Eco, scales = "free")+
  geom_vline(xintercept = 50,color='red', linewidth = 1, alpha = .5)+
  geom_vline(xintercept = 70,color='blue', linewidth = 1, alpha = .5)+
  theme_classic()

Warning: Removed 1 row containing non-finite outside the scale range
(`stat_boxplot()`).

# Average TSS calibration
range(sdms_CV$calibration,na.rm=TRUE)

[1] 0.0900000 0.9836667

mean(sdms_CV$calibration,na.rm=TRUE)

[1] 0.6529396

sd(sdms_CV$calibration,na.rm=TRUE)

[1] 0.1305807

# Average TSS validation
mean(sdms_CV$validation,na.rm=TRUE)

[1] 0.6152985

sd(sdms_CV$validation,na.rm=TRUE)

[1] 0.133801

# Average sensitivity
mean(sdms_CV$sensitivity,na.rm=TRUE)

[1] 86.4879

sd(sdms_CV$sensitivity,na.rm=TRUE)

[1] 7.442999

3 Explore relationships

3.1 N shifts per parameter

ggplot(bioshifts %>% 
         group_by(Param) %>%
         tally, aes(x = n, y = Param))+
  geom_col()+
  theme_classic()+
  geom_text(aes(label = n), vjust = -0.5)+
  labs(x="N shifts",y="N obs")

ggplot(bioshifts %>% 
         group_by(Param, Eco) %>%
         tally, aes(x = n, y = Param))+
  geom_col()+
  theme_classic()+
  geom_text(aes(label = n), vjust = -0.5)+
  labs(x="N shifts",y="N obs")+
  facet_wrap(Eco~., scales = "free")

bioshifts %>% 
  group_by(Param, class) %>%
  tally

# A tibble: 39 × 3
# Groups:   Param [3]
   Param class               n
   <chr> <chr>           <int>
 1 LE    Actinopteri        74
 2 LE    Asteroidea         14
 3 LE    Aves               14
 4 LE    Bivalvia           12
 5 LE    Echinoidea          6
 6 LE    Elasmobranchii      7
 7 LE    Florideophyceae    11
 8 LE    Gastropoda         18
 9 LE    Insecta           132
10 LE    Liliopsida          5
# ℹ 29 more rows

4 How related are observed shifts and climate exposure metrics?

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

mydata$Metric <- factor(mydata$Metric, levels = rev(c("Biotic vel.",  "Bioclimatic vel.")))

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

plotlist <- lapply(1:nrow(plot_cat), function(i){
  toplot <- mydata%>%
    filter(Eco == plot_cat$realms[i],
           Param == plot_cat$Param[i])
  ggplot(toplot, 
         aes(x = Velocity, fill = Metric, color = Metric))+
    ggtitle(paste(plot_cat$realm[i],plot_cat$Param[i],sep = " - "))+
    geom_histogram(alpha = .5)+
    theme_classic()+
    # scale_x_continuous(limits = quantile(toplot$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)

ggplot(bioshifts, 
       aes(x = log(signoise), fill = Param, color = Param))+
  geom_histogram(alpha = .5)+
  theme_classic()+
  geom_hline(yintercept = 0)+
  # scale_x_continuous(limits = quantile(toplot$Velocity, c(0.1, 0.9), na.rm = TRUE))+
  theme(aspect.ratio=1,
        plot.title = element_text(hjust = 0.5))

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: Removed 28 rows containing non-finite outside the scale range
(`stat_bin()`).

quantile(bioshifts$signoise, na.rm = TRUE)

          0%          25%          50%          75%         100% 
1.527593e-13 4.591534e-01 9.892635e-01 1.703393e+00 7.977586e+00 

5 Bioclimatic velocity

5.1 Angle

5.1.1 Frequency

5.1.1.1 Hemisphere

# Define the breaks
breaks <- seq(0, 360, length.out = 25)

# Calculate bin width
bin_width <- diff(breaks)[1]

# Bin the data and calculate midpoints and counts
hist_data <- bioshifts %>%
  group_by(Hemisphere) %>%
  mutate(bin = cut(bv.mat.angle, breaks = breaks, include.lowest = TRUE, right = FALSE, labels = FALSE)) %>%
  group_by(bin,Hemisphere) %>%
  summarise(counts = n()) %>%
  mutate(midpoints = (breaks[bin] + breaks[bin + 1]) / 2)

`summarise()` has grouped output by 'bin'. You can override using the `.groups`
argument.

# Select and arrange the relevant columns
hist_data <- hist_data %>%
  select(midpoints, counts, Hemisphere) %>%
  arrange(midpoints)

Adding missing grouping variables: `bin`

ggplot(hist_data, aes(x = midpoints, y = counts))+
  ggtitle("Frequence of bioclimatic direction")+
  geom_bar(stat = "identity", aes(fill=Hemisphere)) +
  scale_x_continuous(breaks = seq(from=0, to=359, by=45), limits=c(0, 360)) +
  coord_polar() +
  theme_minimal()+
  labs(x="",y="")+
  theme(legend.position = "none",
        panel.grid.major.x = element_blank()
  )+
  facet_wrap(.~Hemisphere, scales = "free")

5.1.1.2 Hemisphere vs Eco

# Bin the data and calculate midpoints and counts
hist_data <- bioshifts %>%
  group_by(Eco, Hemisphere) %>%
  mutate(bin = cut(bv.mat.angle, breaks = breaks, include.lowest = TRUE, right = FALSE, labels = FALSE)) %>%
  group_by(Eco,bin, Hemisphere) %>%
  summarise(counts = n()) %>%
  mutate(midpoints = (breaks[bin] + breaks[bin + 1]) / 2)

`summarise()` has grouped output by 'Eco', 'bin'. You can override using the
`.groups` argument.

# Select and arrange the relevant columns
hist_data <- hist_data %>%
  select(Eco,midpoints, counts, Hemisphere) %>%
  arrange(midpoints)

Adding missing grouping variables: `bin`

ggplot(hist_data, aes(x = midpoints, y = counts,
                      fill = Eco))+
  ggtitle("Frequence of bioclimatic direction")+
  geom_bar(stat = "identity") +
  scale_x_continuous(breaks = seq(from=0, to=359, by=45), limits=c(0, 360)) +
  coord_polar() +
  scale_fill_manual(values = c("dodgerblue","orange4"))+
  theme_minimal()+
  labs(x="",y="")+
  theme(legend.position = "none",
        panel.grid.major.x = element_blank()
  )+
  facet_wrap(Hemisphere~Eco, scales = "free")

5.1.1.3 Hemisphere vs Eco vs Param

# Bin the data and calculate midpoints and counts
hist_data <- bioshifts %>%
  group_by(Eco,Param,Hemisphere) %>%
  mutate(bin = cut(bv.mat.angle, breaks = breaks, include.lowest = TRUE, right = FALSE, labels = FALSE)) %>%
  group_by(Eco,Param,bin,Hemisphere) %>%
  summarise(counts = n()) %>%
  mutate(midpoints = (breaks[bin] + breaks[bin + 1]) / 2)

`summarise()` has grouped output by 'Eco', 'Param', 'bin'. You can override
using the `.groups` argument.

# Select and arrange the relevant columns
hist_data <- hist_data %>%
  select(Eco,midpoints, counts,Hemisphere) %>%
  arrange(midpoints)

Adding missing grouping variables: `Param`, `bin`

ggplot(hist_data, aes(x = midpoints, y = counts,
                      fill = Eco))+
  ggtitle("Frequence of bioclimatic direction")+
  geom_bar(stat = "identity") +
  scale_x_continuous(breaks = seq(from=0, to=359, by=45), limits=c(0, 360)) +
  coord_polar() +
  scale_fill_manual(values = c("dodgerblue","orange4"))+
  theme_minimal()+
  labs(x="",y="")+
  theme(legend.position = "none",
        panel.grid.major.x = element_blank()
  )+
  facet_wrap(Hemisphere+Eco~Param, scales = "free", drop = FALSE, ncol = 3)

5.1.2 Magnitude

5.1.2.1 Hemisphere

ggplot(bioshifts, aes(x = bv.mat.angle, 
                      y = abs(bv.mat2), color = Hemisphere))+
  ggtitle("Absolute magnitude of bioclimatic direction (km/year)")+
  geom_point(alpha=.1) +
  scale_x_continuous(breaks = seq(from=0, to=359, by=45), limits=c(0, 360)) +
  coord_polar() +
  theme_minimal()+
  labs(x="",y="")+
  theme(legend.position = "none",
        panel.grid.major.x = element_blank()
  )+
  facet_wrap(.~Hemisphere, scales = "free")

5.1.2.2 Hemisphere vs Eco

ggplot(bioshifts, aes(x = bv.mat.angle, 
                      y = abs(bv.mat2),
                      color = Eco))+
  ggtitle("Absolute magnitude of bioclimatic direction (km/year)")+
  geom_point(alpha=.1) +
  scale_color_manual(values = c("dodgerblue","orange4"))+
  scale_x_continuous(breaks = seq(from=0, to=359, by=45), limits=c(0, 360)) +
  coord_polar() +
  theme_minimal()+
  labs(x="",y="")+
  theme(legend.position = "none",
        panel.grid.major.x = element_blank()
  )+
  facet_wrap(Hemisphere~Eco, scales = "free")

5.1.2.3 Hemisphere vs Eco vs Param

ggplot(bioshifts, aes(x = bv.mat.angle, 
                      y = abs(bv.mat2),
                      color = Eco))+
  ggtitle("Absolute magnitude of bioclimatic direction (km/year)")+
  geom_point(alpha=.1) +
  scale_color_manual(values = c("dodgerblue","orange4"))+
  scale_x_continuous(breaks = seq(from=0, to=359, by=45), limits=c(0, 360)) +
  coord_polar() +
  theme_minimal()+
  labs(x="",y="")+
  theme(legend.position = "none",
        panel.grid.major.x = element_blank()
  )+
  facet_wrap(Hemisphere+Eco~Param, scales = "free", drop = FALSE, ncol = 3)

5.2 Direction match

ggplot(bioshifts %>%
         group_by(match_direction_class) %>%
         count %>%
         na.omit %>% 
         data.frame %>%
         mutate(n = (n/sum(n))*100,
                bioclimatic = c("-","-","+","+"),
                biotic = c("-","+","-","+"),
                match_direction_class = match_direction_class), 
       aes(x = bioclimatic, y = biotic)) + 
  geom_point(aes(size = n, fill = n), alpha = 0.75, shape = 22) + 
  scale_size_continuous(limits = c(1, 100), range = c(1,50), breaks = c(1,10,50,75)) +
  labs( x= "Bioclimatic direction", y = "Biotic direction", 
        size = "Relative frequency (%)", fill = "")  + 
  theme(legend.key=element_blank(), 
        axis.text = element_text(size = 18),
        panel.background = element_blank(), 
        panel.border = element_rect(colour = "black", fill = NA, size = 1), 
        legend.position = "none") 

Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
ℹ Please use the `linewidth` argument instead.

ggplot(bioshifts %>%
         group_by(Eco, match_direction_class) %>%
         count %>%
         na.omit %>% 
         data.frame() %>%
         mutate(prop = (n/sum(n))*100,
                biotic = substr(match_direction_class, 1, 3),
                bioclimatic = substr(match_direction_class, 4, 6))%>%
         mutate(biotic = ifelse(grepl("pos",.$biotic),"Poleward","Equatoward"),
                bioclimatic = ifelse(grepl("pos",.$bioclimatic),"Poleward","Equatoward")), 
       aes(x = bioclimatic, y = biotic)) + 
  geom_point(aes(size = prop, fill = n), alpha = 0.75, shape = 22) + 
  geom_text(size = 5, aes(label = n))+
  scale_size_continuous(range=c(1,30), breaks = c(10, 15, 20, 30)) +
  labs(x= "Bioclimatic direction", y = "Biotic direction", 
       size = "Relative frequency (%)", fill = "") +
  theme(panel.background = element_blank(), 
        panel.border = element_rect(colour = "black", fill = NA, size = 1)) +
  guides(fill = FALSE)+
  facet_wrap(.~Eco)

Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
of ggplot2 3.3.4.

5.3 Scatter plots

5.3.1 Eco

new_data <- bioshifts %>%
  group_by(ID, class, Param, Eco) %>%
  select(ID, class, Param, bv.mat2,ShiftRate, Eco) %>%
  summarise(bv_sa = mean(bv.mat2),
            shift_sa  = mean(ShiftRate),
            bv_sa_sd = sd(bv.mat2),
            shift_sa_sd  = sd(ShiftRate))


p1 <- ggplot(bioshifts, 
             aes(x = bv.mat2, y = ShiftRate, color = class))+
  geom_point(alpha = .5, size =2)+
  labs(x="Bioclimatic velocity (km/yr SDM suitability)", 
       y="Observed range shift velocity (km/yr)")+
  theme_classic() +
  theme(aspect.ratio=1,
        plot.title = element_text(hjust = 0.5)) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed")+
  facet_wrap(Eco~Param, scales = "free")

p1

p1 + geom_point(data = new_data,
                aes(x=bv_sa,y=shift_sa), color = "black", alpha = .5, size =5)

p1 + geom_point(data = new_data,
                aes(x=bv_sa,y=shift_sa), color = "black", alpha = .15, size =5)+
  geom_pointrange(data = new_data,
                  aes(x=bv_sa,y=shift_sa,
                      xmin = bv_sa-bv_sa_sd,
                      xmax = bv_sa+bv_sa_sd), color = "black", alpha = .15, size =1) +
  geom_pointrange(data = new_data,
                  aes(x=bv_sa,y=shift_sa,
                      ymin = shift_sa-shift_sa_sd,
                      ymax = shift_sa+shift_sa_sd), color = "black", alpha = .15, size =1)

5.3.2 Taxonomic class

ggplot(bioshifts %>%
         filter(Eco == "Ter"), 
       aes(x = bv.mat, y = ShiftRate, color = Duration))+
  ggtitle("Terrestrials")+
  geom_point(alpha = .5, size = 3)+
  scale_color_gradient(low="blue",high = "red")+
  labs(x="Bioclimatic velocity (km/yr SDM suitability)", 
       y="Observed range shift velocity (km/yr)")+
  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")+
  facet_grid(Param~class, scales = "free")

5.3.3 Mismatch

plot_data <- bioshifts %>%
  filter(Eco == "Ter")

ggplot(plot_data, aes(x = mismatch_percent, y = class))+
  geom_boxplot(alpha = .5)+
  labs(x="Mismatch between observed shift and bioclimatic velocity (%)\n<< Faster bioclimatic velocity           Faster observed shift velocity >>",
       title = "Terrestrial")+
  scale_x_continuous(limits = quantile(plot_data$mismatch_percent, c(0.1, 0.9), na.rm = TRUE))+
  geom_vline(xintercept = c(0,100), linetype = "dashed")+
  theme_classic() +
  facet_wrap(.~Param, scales = "free_x", ncol = 4)

Warning: Removed 120 rows containing non-finite outside the scale range
(`stat_boxplot()`).

plot_data <- bioshifts %>%
  filter(Eco == "Mar")

ggplot(plot_data, aes(x = mismatch_percent, y = class))+
  geom_boxplot(alpha = .5)+
  labs(x="Mismatch between observed shift and bioclimatic velocity (%)\n<< Faster bioclimatic velocity           Faster observed shift velocity >>",
       title = "Marine")+
  scale_x_continuous(limits = quantile(plot_data$mismatch_percent, c(0.1, 0.9), na.rm = TRUE))+
  geom_vline(xintercept = c(0,100), linetype = "dashed")+
  theme_classic() +
  facet_wrap(.~Param, scales = "free_x", ncol = 4)

Warning: Removed 199 rows containing non-finite outside the scale range
(`stat_boxplot()`).

Warning: Removed 3 rows containing missing values or values outside the scale range
(`geom_vline()`).

ggplot(bioshifts, aes(x = ShiftRate, y = mismatch))+
  geom_point()+
  scale_x_continuous(limits = quantile(bioshifts$mismatch_percent, c(0.1, 0.9), na.rm = TRUE))+
  theme_classic()

The mismatch is highly connected with the magnitude of shift rate. This is because most of the bioclimatic velocities are very small.

ggplot(bioshifts, aes(x = mismatch_percent, y = bv.mat))+
  geom_point()+
  scale_x_continuous(limits = quantile(bioshifts$mismatch_percent, c(0.1, 0.9), na.rm = TRUE))+
  theme_classic()

Warning: Removed 187 rows containing missing values or values outside the scale range
(`geom_point()`).

Likewise, we can see that mismatch percentage tends to be closer to 100% when bioclimatic velocities are very low or equal to zero. In the section How related are observed shifts and climate exposure metrics? we can see there is a strong spike of bioclimatic velocities close to zero.

5.3.4 Duration

Can duration explain the decoupling between observed shifts and bioclimatic velocities?

ggplot(bioshifts, 
       aes(x = Duration, y = (abs(mismatch_percent))))+
  geom_point(alpha = .5)+
  labs(x="Shift duration", 
       y="Mismatch between observed shift and bioclimatic velocity (%)")+
  scale_y_continuous(transform = "log",
                     labels = function(x) round(exp(log(x)),0)
  )+
  theme_classic() +
  geom_smooth(method = "lm")+
  facet_wrap(Eco~Param, scales = "free")

5.3.5 N periods

Can N periods explain the decoupling between observed shifts and bioclimatic velocities?

ggplot(bioshifts %>%
         select(mismatch_percent,N_periodes,Param,Eco) %>%
         na.omit, 
       aes(x = N_periodes, y = (abs(mismatch_percent))))+
  geom_point(alpha = .5)+
  labs(x="N periodes", 
       y="Mismatch between observed shift and bioclimatic velocity (%)")+
  scale_y_continuous(transform = "log",
                     labels = function(x) round(exp(log(x)),0)
  )+
  scale_x_continuous(transform = "log",
                     labels = function(x) round(exp(log(x)),0)
  )+
  theme_classic() +
  geom_smooth(method = "lm")+
  facet_wrap(Eco~Param, scales = "free")

5.3.6 Category

Can Category explain the decoupling between observed shifts and bioclimatic velocities?

ggplot(bioshifts, 
       aes(y = Category, x = abs(mismatch_percent)))+
  geom_boxplot(alpha = .5)+
  labs(y="Samplying category", 
       x="Mismatch between observed shift and bioclimatic velocity (%)")+
  theme_classic() +
  scale_x_continuous(transform = "log",
                     labels = function(x) round(exp(log(x)),0)
  )+
  geom_vline(xintercept = 100, linetype = "dashed")+
  facet_wrap(Eco~Param, scales = "free_x")

5.3.7 SDM performance

Can SDM performance explain the decoupling between observed shifts and bioclimatic velocities?

ggplot(bioshifts, 
       aes(x = calibration, y = abs(mismatch_percent)))+
  geom_point(alpha = .5)+
  labs(x="TSS", 
       y="Mismatch between observed shift and bioclimatic velocity (%)")+
  theme_classic() +
  scale_y_continuous(transform = "log",
                     labels = function(x) round(exp(log(x)),0)
  )+
  geom_smooth(method = "lm")+
  geom_hline(yintercept = 0, linetype = "dashed")+
  facet_grid(Eco~Param, scales = "free_y")

5.3.8 Signal to noise

Can signal to noise explain the decoupling between observed shifts and bioclimatic velocities?

ggplot(bioshifts, 
       aes(x = signoise, y = abs(mismatch_percent)))+
  geom_point(alpha = .5)+
  labs(x="Sig:Noise", 
       y="Mismatch between observed shift and bioclimatic velocity (%)")+
  theme_classic() +
  geom_smooth(method = "lm")+
  scale_y_continuous(transform = "log",
                     labels = function(x) round(exp(log(x)),0)
  )+
  facet_wrap(Eco~Param, scales = "free")

ggplot(bioshifts %>%
         select(mismatch,signoise,Param,Eco) %>%
         mutate(signoise = ifelse(signoise>1,">sig", ">noise")) %>%
         na.omit, 
       aes(x = signoise, y = abs(mismatch)))+
  geom_boxplot()+
  labs(x="Sig:Noise", 
       y="Mismatch between observed shift and bioclimatic velocity (%)")+
  theme_classic() +
  # geom_smooth(method = "lm")+
  # geom_hline(yintercept = 0, linetype = "dashed")+
  # geom_vline(xintercept = 1, linetype = "dashed")+
  facet_wrap(Eco~Param, scales = "free_y")

There is a better match in observations and predictions when Std error of models are low, and the mismatch increases with increasing std error!

Plot only the data with low std error (std. error < 100)

hist(bioshifts$signoise)

quantile(bioshifts$signoise, na.rm =TRUE)

          0%          25%          50%          75%         100% 
1.527593e-13 4.591534e-01 9.892635e-01 1.703393e+00 7.977586e+00 

ggplot(bioshifts %>%
         select(Param,Eco,bv.mat,bv.mat.error,ShiftRate,Eco,Duration,Param,class,signoise) %>%
         filter(signoise >= 1) %>%
         na.omit, 
       aes(x = bv.mat, y = ShiftRate, color = class))+
  geom_point(alpha = .5, size =3)+
  labs(x="Bioclimatic velocity (km/yr SDM suitability)", 
       y="Observed range shift velocity (km/yr)")+
  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")+
  facet_wrap(Eco~Param, scales = "free")

5.3.9 Noise

Can signal to noise explain the decoupling between observed shifts and bioclimatic velocities?

ggplot(bioshifts %>%
         select(mismatch,bv.mat.error,Param,Eco) %>%
         na.omit, 
       aes(x = (bv.mat.error), y = abs(mismatch)))+
  geom_point(alpha = .5)+
  labs(x="Noise (log+1)", 
       y="Mismatch between observed shift and bioclimatic velocity (%)")+
  theme_classic() +
  geom_smooth(method = "lm")+
  facet_wrap(Eco~Param, scales = "free")

5.3.10 Traits

Can traits explain the mismatch between observed shifts and bioclimatic velocities?

plot_data <- bioshifts %>%
  filter(Eco == "Ter") %>%
  select(Param,mismatch,bv.mat,ShiftRate,Eco,Locomotion_mode,Param,class) %>%
  na.omit

ggplot(plot_data, aes(x = abs(mismatch), y = Locomotion_mode))+
  geom_boxplot(alpha = .5)+
  scale_x_continuous(limits = quantile(plot_data$mismatch, c(0.1, 0.9), na.rm = TRUE))+
  labs(x="Mismatch between observed shift and bioclimatic velocity (%)", title="Terrestrial")+
  theme_classic() +
  geom_vline(xintercept = 0, linetype = "dashed")+
  facet_wrap(.~Param, scales = "free_x", ncol = 4)

Warning: Removed 48 rows containing non-finite outside the scale range
(`stat_boxplot()`).

plot_data <- bioshifts %>%
  filter(Eco == "Mar") %>%
  select(Param,mismatch,bv.mat,ShiftRate,Eco,Locomotion_mode,Param,class) %>%
  na.omit

ggplot(plot_data, aes(x = abs(mismatch), y = Locomotion_mode))+
  geom_boxplot(alpha = .5)+
  scale_x_continuous(limits = quantile(plot_data$mismatch, c(0.1, 0.9), na.rm = TRUE))+
  labs(x="Mismatch between observed shift and bioclimatic velocity (%)", title="Marine")+
  theme_classic() +
  geom_vline(xintercept = 0, linetype = "dashed")+
  facet_wrap(.~Param, scales = "free_x", ncol = 4)

Warning: Removed 96 rows containing non-finite outside the scale range
(`stat_boxplot()`).

5.3.11 Landscape connectivity

ggplot(bioshifts, 
       aes(x = impeded_prop*100, y = (abs(mismatch_percent))))+
  geom_point(alpha = .5)+
  labs(x="Landscape connectivity\nImpeded (%)", 
       y="Mismatch between observed shift and bioclimatic velocity (%)")+
  scale_y_continuous(transform = "log",
                     labels = function(x) round(exp(log(x)),0)
  )+
  theme_classic() +
  geom_smooth(method = "lm")+
  facet_wrap(Eco~Param, scales = "free")

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 1221 rows containing non-finite outside the scale range
(`stat_smooth()`).

Warning: Removed 1168 rows containing missing values or values outside the scale range
(`geom_point()`).

ggplot(bioshifts, 
       aes(x = diffuse_prop*100, y = (abs(mismatch_percent))))+
  geom_point(alpha = .5)+
  labs(x="Landscape connectivity\nDiffused (%)", 
       y="Mismatch between observed shift and bioclimatic velocity (%)")+
  scale_y_continuous(transform = "log",
                     labels = function(x) round(exp(log(x)),0)
  )+
  theme_classic() +
  geom_smooth(method = "lm")+
  facet_wrap(Eco~Param, scales = "free")

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 1221 rows containing non-finite outside the scale range
(`stat_smooth()`).
Removed 1168 rows containing missing values or values outside the scale range
(`geom_point()`).

ggplot(bioshifts, 
       aes(x = intensified_prop*100, y = (abs(mismatch_percent))))+
  geom_point(alpha = .5)+
  labs(x="Landscape connectivity\nIntensified (%)", 
       y="Mismatch between observed shift and bioclimatic velocity (%)")+
  scale_y_continuous(transform = "log",
                     labels = function(x) round(exp(log(x)),0)
  )+
  theme_classic() +
  geom_smooth(method = "lm")+
  facet_wrap(Eco~Param, scales = "free")

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 1221 rows containing non-finite outside the scale range
(`stat_smooth()`).
Removed 1168 rows containing missing values or values outside the scale range
(`geom_point()`).

ggplot(bioshifts, 
       aes(x = channelized_prop*100, y = (abs(mismatch_percent))))+
  geom_point(alpha = .5)+
  labs(x="Landscape connectivity\nChannelized (%)", 
       y="Mismatch between observed shift and bioclimatic velocity (%)")+
  scale_y_continuous(transform = "log",
                     labels = function(x) round(exp(log(x)),0)
  )+
  theme_classic() +
  geom_smooth(method = "lm")+
  facet_wrap(Eco~Param, scales = "free")

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 1221 rows containing non-finite outside the scale range
(`stat_smooth()`).
Removed 1168 rows containing missing values or values outside the scale range
(`geom_point()`).

5.4 Climate velocity

ggplot(bioshifts %>%
         select(Param,Eco,vel_mat_lat_median,vel_mat_lat_median_25km,ShiftRate,Eco,Duration,Param,class,sensitivity) %>%
         na.omit, 
       aes(x = vel_mat_lat_median, y = ShiftRate, color = class))+
  geom_point(alpha = .5, size =3)+
  ggtitle("1km resolution for terrestrial species")+
  labs(x="Climatic velocity (km/yr Temperature)", 
       y="Observed range shift velocity (km/yr)")+
  theme_classic() +
  theme(aspect.ratio=1,
        plot.title = element_text(hjust = 0.5)) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed")+
  facet_wrap(Eco~Param, scales = "free")

ggplot(bioshifts %>%
         select(Param,Eco,vel_mat_lat_median,vel_mat_lat_median_25km,ShiftRate,Eco,Duration,Param,class,sensitivity) %>%
         na.omit, 
       aes(x = vel_mat_lat_median_25km, y = ShiftRate, color = class))+
  geom_point(alpha = .5, size =3)+
  ggtitle("25km resolution for terrestrial species")+
  labs(x="Climatic velocity (km/yr Temperature)", 
       y="Observed range shift velocity (km/yr)")+
  theme_classic() +
  theme(aspect.ratio=1,
        plot.title = element_text(hjust = 0.5)) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed")+
  facet_wrap(Eco~Param, scales = "free")

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_median,ShiftRate,Eco,Duration,Param,class) %>%
           na.omit, 
         aes(x = vel_mat_lat_median, 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)

5.5 Edge climate velocity

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(ev.mat,ShiftRate,Eco,Duration,Param,class) %>%
           na.omit, 
         aes(x = ev.mat, y = ShiftRate, size = Duration, color = class))+
    geom_point(alpha = .5)+
    labs(x="Edge 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)

Compare across taxonomic classes

ggplot(bioshifts %>%
         filter(Eco == "Mar") %>% 
         select(ev.mat,ShiftRate,Eco,Duration,Param,class) %>%
         na.omit, 
       aes(x = ev.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")

ggplot(bioshifts %>%
         filter(Eco == "Ter") %>% 
         select(ev.mat,ShiftRate,Eco,Duration,Param,class) %>%
         na.omit, 
       aes(x = ev.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")

5.6 Edge climate velocity vs Bioclimatic velocity

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(ev.mat,bv.mat,Eco,Duration,Param,class) %>%
           na.omit, 
         aes(x = ev.mat, y = bv.mat, size = Duration, color = class))+
    geom_point(alpha = .5)+
    labs(x="Edge 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)

6 Models

Run separate models for marine and terrestrials

# hist(log(abs(bioshifts$ShiftRate-bioshifts$bv.mat)))
# hist(log(abs(bioshifts$ShiftRate/bioshifts$bv.mat)))

ggplot(bioshifts, aes(mismatch/(ShiftRate)))+
  geom_histogram(bins=100)+
  xlim(c(-10,10))+
  theme_classic()+
  geom_vline(xintercept = c(-1,1))

Warning: Removed 185 rows containing non-finite outside the scale range
(`stat_bin()`).

ggplot(bioshifts, aes((ShiftRate-ev.mat)/abs(ShiftRate)))+
  geom_histogram(bins=100)+
  xlim(c(-10,10))+
  theme_classic()+
  geom_vline(xintercept = c(-1,1))

Warning: Removed 554 rows containing non-finite outside the scale range
(`stat_bin()`).

7 0) Bioclimatic velocities as predictors of observed range shifts

model_formula <- c("ShiftRate ~  bv.mat2 * Param + N_periodes + log(LatExtentk) + Grain_size + Obs_type + (1 | class) + (1 | Eco)")
# model_formula <- c("ShiftRate ~  bv.mat2 * Param + (1 | class) + (1 | Eco)")

my_model <- glmmTMB(as.formula(model_formula),
                    # weight = calibration,
                    data = bioshifts)

summary(my_model)

 Family: gaussian  ( identity )
Formula:          ShiftRate ~ bv.mat2 * Param + N_periodes + log(LatExtentk) +  
    Grain_size + Obs_type + (1 | class) + (1 | Eco)
Data: bioshifts

     AIC      BIC   logLik deviance df.resid 
  9346.5   9426.7  -4658.3   9316.5     1528 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 class    (Intercept)  1.195   1.093   
 Eco      (Intercept)  1.033   1.016   
 Residual             24.106   4.910   
Number of obs: 1543, groups:  class, 22; Eco, 2

Dispersion estimate for gaussian family (sigma^2): 24.1 

Conditional model:
                       Estimate Std. Error z value Pr(>|z|)    
(Intercept)          -11.730439   1.902875  -6.165 7.07e-10 ***
bv.mat2               -0.776543   0.283130  -2.743  0.00609 ** 
ParamO                -0.713945   0.478691  -1.491  0.13584    
ParamTE               -1.578648   0.600914  -2.627  0.00861 ** 
N_periodes            -0.001845   0.001290  -1.430  0.15283    
log(LatExtentk)        1.823607   0.234469   7.778 7.39e-15 ***
Grain_sizemoderate     1.728861   0.565386   3.058  0.00223 ** 
Grain_sizelarge        0.659908   0.540732   1.220  0.22231    
Grain_sizevery_large   0.002985   0.376570   0.008  0.99368    
Obs_typeoccurrence     0.928906   0.563527   1.648  0.09927 .  
bv.mat2:ParamO         1.017678   0.311721   3.265  0.00110 ** 
bv.mat2:ParamTE        1.304351   0.698537   1.867  0.06187 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Actual    = bioshifts$ShiftRate
Predicted = predict(my_model, type="response")
Residuals = residuals(my_model)

efronRSquared(residual = Residuals, 
              predicted = Predicted, 
              statistic = "EfronRSquared")

EfronRSquared 
        0.115 

plot(effects::Effect(c("bv.mat2","Param"), my_model))

## terrestrial
model_formula <- c("ShiftRate ~  bv.mat2 * Param + N_periodes + log(LatExtentk) + Grain_size + Obs_type + (1 | class)")

subdata <- bioshifts %>%
  filter(Eco=="Ter")

my_model <- glmmTMB(as.formula(model_formula),
                    # weight = calibration,
                    data = subdata)

summary(my_model)

 Family: gaussian  ( identity )
Formula:          ShiftRate ~ bv.mat2 * Param + N_periodes + log(LatExtentk) +  
    Grain_size + Obs_type + (1 | class)
Data: subdata

     AIC      BIC   logLik deviance df.resid 
  2737.6   2799.0  -1354.8   2709.6      578 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 class    (Intercept) 0.04946  0.2224  
 Residual             5.66064  2.3792  
Number of obs: 592, groups:  class, 7

Dispersion estimate for gaussian family (sigma^2): 5.66 

Conditional model:
                       Estimate Std. Error z value Pr(>|z|)   
(Intercept)          -0.6448469  1.8964223  -0.340  0.73383   
bv.mat2              -0.7295613  0.2664620  -2.738  0.00618 **
ParamO               -0.6892371  0.4087140  -1.686  0.09173 . 
ParamTE              -0.8348913  0.5543296  -1.506  0.13203   
N_periodes            0.0002355  0.0008161   0.289  0.77290   
log(LatExtentk)       0.2319926  0.2436277   0.952  0.34097   
Grain_sizemoderate   -0.6835747  0.4579309  -1.493  0.13550   
Grain_sizelarge      -1.0604148  0.4666653  -2.272  0.02307 * 
Grain_sizevery_large -0.0680719  0.5152468  -0.132  0.89489   
Obs_typeoccurrence    0.7649229  0.5724574   1.336  0.18148   
bv.mat2:ParamO        0.6426952  0.3399177   1.891  0.05866 . 
bv.mat2:ParamTE      -0.1880607  1.5089677  -0.125  0.90082   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Actual    = log1p(abs(subdata$mismatch))
Predicted = predict(my_model, type="response")
Residuals = residuals(my_model)

efronRSquared(residual = Residuals, 
              predicted = Predicted, 
              statistic = "EfronRSquared")

EfronRSquared 
       0.0942 

plot(effects::Effect(c("bv.mat2","Param"), my_model))

## marine

subdata <- bioshifts %>%
  filter(Eco=="Mar")

my_model <- glmmTMB(as.formula(model_formula),
                    # weight = calibration,
                    # family = Gamma(link = "log"),
                    data = subdata)

summary(my_model)

 Family: gaussian  ( identity )
Formula:          ShiftRate ~ bv.mat2 * Param + N_periodes + log(LatExtentk) +  
    Grain_size + Obs_type + (1 | class)
Data: subdata

     AIC      BIC   logLik deviance df.resid 
  6118.5   6186.5  -3045.2   6090.5      937 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 class    (Intercept)  0.8723  0.9339  
 Residual             35.0304  5.9186  
Number of obs: 951, groups:  class, 15

Dispersion estimate for gaussian family (sigma^2):   35 

Conditional model:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)          -13.0129     2.4846  -5.237 1.63e-07 ***
bv.mat2               -0.6093     0.4068  -1.498   0.1342    
ParamO                -1.0902     0.7882  -1.383   0.1666    
ParamTE               -2.3625     0.8806  -2.683   0.0073 ** 
N_periodes             0.0254     0.0173   1.468   0.1422    
log(LatExtentk)        2.1951     0.3770   5.823 5.78e-09 ***
Grain_sizemoderate     2.3700     1.3684   1.732   0.0833 .  
Grain_sizelarge        1.7596     0.7694   2.287   0.0222 *  
Grain_sizevery_large  -0.3687     0.5784  -0.637   0.5238    
Obs_typeoccurrence     0.6579     0.7919   0.831   0.4061    
bv.mat2:ParamO         0.8985     0.4368   2.057   0.0397 *  
bv.mat2:ParamTE        1.0753     0.8970   1.199   0.2306    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Actual    = log1p(abs(subdata$mismatch))
Predicted = predict(my_model, type="response")
Residuals = residuals(my_model)

efronRSquared(residual = Residuals, 
              predicted = Predicted, 
              statistic = "EfronRSquared")

EfronRSquared 
        0.112 

plot(effects::Effect(c("bv.mat2","Param"), my_model))

8 1) SDM performance

## overall
model_formula <- c("log1p(abs(mismatch)) ~  validation * Param + N_periodes + log(LatExtentk) + Grain_size + Obs_type + (1 | class) + (1 | Eco)")

my_model <- glmmTMB(as.formula(model_formula),
                    # weight = calibration,
                    # family = Gamma(link = "log"),
                    data = bioshifts)

summary(my_model)

 Family: gaussian  ( identity )
Formula:          
log1p(abs(mismatch)) ~ validation * Param + N_periodes + log(LatExtentk) +  
    Grain_size + Obs_type + (1 | class) + (1 | Eco)
Data: bioshifts

     AIC      BIC   logLik deviance df.resid 
  2932.1   3011.3  -1451.1   2902.1     1435 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 class    (Intercept) 0.02043  0.1429  
 Eco      (Intercept) 0.13391  0.3659  
 Residual             0.42455  0.6516  
Number of obs: 1450, groups:  class, 21; Eco, 2

Dispersion estimate for gaussian family (sigma^2): 0.425 

Conditional model:
                       Estimate Std. Error z value Pr(>|z|)    
(Intercept)          -2.6405640  0.3884408  -6.798 1.06e-11 ***
validation            0.6219625  0.2796196   2.224 0.026127 *  
ParamO                0.2756460  0.2073732   1.329 0.183773    
ParamTE              -0.2331722  0.3797459  -0.614 0.539201    
N_periodes           -0.0006639  0.0001774  -3.743 0.000182 ***
log(LatExtentk)       0.4357221  0.0311318  13.996  < 2e-16 ***
Grain_sizemoderate    0.4395040  0.0773363   5.683 1.32e-08 ***
Grain_sizelarge       0.1632741  0.0763347   2.139 0.032442 *  
Grain_sizevery_large  0.2450071  0.0514198   4.765 1.89e-06 ***
Obs_typeoccurrence    0.2369733  0.0769316   3.080 0.002068 ** 
validation:ParamO    -0.6241761  0.3180401  -1.963 0.049696 *  
validation:ParamTE   -0.0194135  0.5662261  -0.034 0.972649    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Actual    = log1p(abs(bioshifts$mismatch))
Predicted = predict(my_model, type="response")
Residuals = residuals(my_model)

efronRSquared(residual = Residuals, 
              predicted = Predicted, 
              statistic = "EfronRSquared")

EfronRSquared 
        0.325 

plot(effects::Effect(c("validation","Param"), my_model))

## terrestrial
model_formula <- c("log1p(abs(mismatch)) ~  validation * Param + N_periodes + log(LatExtentk) + Grain_size + Obs_type + (1 | class)")

subdata <- bioshifts %>%
  filter(Eco=="Ter")

my_model <- glmmTMB(as.formula(model_formula),
                    # weight = calibration,
                    # family = Gamma(link = "log"),
                    data = subdata)

summary(my_model)

 Family: gaussian  ( identity )
Formula:          
log1p(abs(mismatch)) ~ validation * Param + N_periodes + log(LatExtentk) +  
    Grain_size + Obs_type + (1 | class)
Data: subdata

     AIC      BIC   logLik deviance df.resid 
   824.6    884.4   -398.3    796.6      513 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 class    (Intercept) 0.02558  0.1599  
 Residual             0.25999  0.5099  
Number of obs: 527, groups:  class, 6

Dispersion estimate for gaussian family (sigma^2): 0.26 

Conditional model:
                       Estimate Std. Error z value Pr(>|z|)    
(Intercept)          -1.8276077  0.5530877  -3.304 0.000952 ***
validation            0.9181573  0.4187716   2.193 0.028343 *  
ParamO                0.5442508  0.3025893   1.799 0.072075 .  
ParamTE              -0.0077040  0.9938380  -0.008 0.993815    
N_periodes           -0.0004950  0.0001789  -2.768 0.005646 ** 
log(LatExtentk)       0.3234443  0.0621751   5.202 1.97e-07 ***
Grain_sizemoderate    0.0841308  0.1034659   0.813 0.416146    
Grain_sizelarge      -0.0945054  0.1073790  -0.880 0.378799    
Grain_sizevery_large -0.2357395  0.1125153  -2.095 0.036155 *  
Obs_typeoccurrence   -0.0739081  0.1395651  -0.530 0.596417    
validation:ParamO    -1.1478448  0.4846973  -2.368 0.017876 *  
validation:ParamTE   -0.2671198  1.5673273  -0.170 0.864672    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Actual    = log1p(abs(subdata$mismatch))
Predicted = predict(my_model, type="response")
Residuals = residuals(my_model)

efronRSquared(residual = Residuals, 
              predicted = Predicted,
              statistic = "EfronRSquared")

EfronRSquared 
        0.253 

plot(effects::Effect(c("validation","Param"), my_model))

## marine
model_formula <- c("log1p(abs(mismatch)) ~  validation * Param + (1 | class)")

subdata <- bioshifts %>%
  filter(Eco=="Mar")

my_model <- glmmTMB(as.formula(model_formula),
                    # weight = calibration,
                    # family = Gamma(link = "log"),
                    data = subdata)

summary(my_model)

 Family: gaussian  ( identity )
Formula:          log1p(abs(mismatch)) ~ validation * Param + (1 | class)
Data: subdata

     AIC      BIC   logLik deviance df.resid 
  2208.7   2247.3  -1096.4   2192.7      915 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 class    (Intercept) 0.2544   0.5044  
 Residual             0.6046   0.7776  
Number of obs: 923, groups:  class, 15

Dispersion estimate for gaussian family (sigma^2): 0.605 

Conditional model:
                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)          1.1156     0.2802   3.981 6.85e-05 ***
validation           0.3379     0.4019   0.841    0.400    
ParamO              -0.1300     0.2952  -0.440    0.660    
ParamTE             -0.4293     0.5015  -0.856    0.392    
validation:ParamO   -0.2205     0.4543  -0.485    0.627    
validation:ParamTE   0.1625     0.7340   0.221    0.825    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Actual    = log1p(abs(subdata$mismatch))
Predicted = predict(my_model, type="response")
Residuals = residuals(my_model)

efronRSquared(residual = Residuals, 
              predicted = Predicted, 
              statistic = "EfronRSquared")

EfronRSquared 
        0.122 

plot(effects::Effect(c("validation","Param"), my_model))

At the global model (combining terrestrial and marine shifts), we found a counterintuitive pattern that higher TSS leads to greater mismatches between observations and predictions. The same patterns seems to emerge when looking at terrestrial species, but marine species show no significant effect of TSS.

9 2) Quality of observed shifts

model_formula <- c("log1p(abs(mismatch)) ~  Param + N_periodes + log(LatExtentk) + Category + Grain_size + Obs_type + (1 | class) + (1 | Eco)")
# model_formula <- c("ShiftRate ~  bv.mat2 * Param + (1 | class) + (1 | Eco)")

my_model <- glmmTMB(as.formula(model_formula),
                    # weight = calibration,
                    data = bioshifts)

summary(my_model)

 Family: gaussian  ( identity )
Formula:          
log1p(abs(mismatch)) ~ Param + N_periodes + log(LatExtentk) +  
    Category + Grain_size + Obs_type + (1 | class) + (1 | Eco)
Data: bioshifts

     AIC      BIC   logLik deviance df.resid 
  3099.0   3173.8  -1535.5   3071.0     1529 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 class    (Intercept) 0.04722  0.2173  
 Eco      (Intercept) 0.14613  0.3823  
 Residual             0.41712  0.6458  
Number of obs: 1543, groups:  class, 22; Eco, 2

Dispersion estimate for gaussian family (sigma^2): 0.417 

Conditional model:
                          Estimate Std. Error z value Pr(>|z|)    
(Intercept)             -2.1799366  0.3656605  -5.962 2.50e-09 ***
ParamO                  -0.1503592  0.0660932  -2.275 0.022909 *  
ParamTE                 -0.2367238  0.0794443  -2.980 0.002885 ** 
N_periodes              -0.0008480  0.0001853  -4.577 4.73e-06 ***
log(LatExtentk)          0.4358987  0.0315068  13.835  < 2e-16 ***
CategoryCensusPeriods   -0.1113329  0.0622261  -1.789 0.073588 .  
CategorySurvey-Resurvey -0.3008071  0.0912967  -3.295 0.000985 ***
Grain_sizemoderate       0.4308645  0.0748375   5.757 8.55e-09 ***
Grain_sizelarge          0.1900252  0.0775894   2.449 0.014321 *  
Grain_sizevery_large     0.2207535  0.0587414   3.758 0.000171 ***
Obs_typeoccurrence       0.2767142  0.0802920   3.446 0.000568 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Actual    = bioshifts$ShiftRate
Predicted = predict(my_model, type="response")
Residuals = residuals(my_model)

efronRSquared(residual = Residuals, 
              predicted = Predicted, 
              statistic = "EfronRSquared")

EfronRSquared 
        0.334 

plot(effects::Effect(c("Category"), my_model))

plot(effects::Effect(c("Grain_size"), my_model))

## terrestrial
model_formula <- c("log1p(abs(mismatch)) ~  Param + N_periodes + log(LatExtentk) + Category + Grain_size + Obs_type + (1 | class)")

subdata <- bioshifts %>%
  filter(Eco=="Ter")

my_model <- glmmTMB(as.formula(model_formula),
                    # weight = calibration,
                    data = subdata)

summary(my_model)

 Family: gaussian  ( identity )
Formula:          
log1p(abs(mismatch)) ~ Param + N_periodes + log(LatExtentk) +  
    Category + Grain_size + Obs_type + (1 | class)
Data: subdata

     AIC      BIC   logLik deviance df.resid 
   923.4    980.4   -448.7    897.4      579 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 class    (Intercept) 0.02233  0.1494  
 Residual             0.26150  0.5114  
Number of obs: 592, groups:  class, 7

Dispersion estimate for gaussian family (sigma^2): 0.261 

Conditional model:
                          Estimate Std. Error z value Pr(>|z|)    
(Intercept)             -0.6383857  0.5077025  -1.257  0.20861    
ParamO                  -0.1986633  0.0915506  -2.170  0.03001 *  
ParamTE                 -0.1014641  0.1251060  -0.811  0.41735    
N_periodes              -0.0005998  0.0001887  -3.179  0.00148 ** 
log(LatExtentk)          0.2632098  0.0608897   4.323 1.54e-05 ***
CategoryCensusPeriods   -0.2555933  0.1378212  -1.855  0.06366 .  
CategorySurvey-Resurvey -0.2407791  0.1605248  -1.500  0.13363    
Grain_sizemoderate       0.0961725  0.0980901   0.980  0.32686    
Grain_sizelarge         -0.1171856  0.1038748  -1.128  0.25926    
Grain_sizevery_large    -0.2387291  0.1398906  -1.707  0.08791 .  
Obs_typeoccurrence      -0.0844251  0.1259999  -0.670  0.50283    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Actual    = log1p(abs(subdata$mismatch))
Predicted = predict(my_model, type="response")
Residuals = residuals(my_model)

efronRSquared(residual = Residuals, 
              predicted = Predicted, 
              statistic = "EfronRSquared")

EfronRSquared 
        0.229 

plot(effects::Effect(c("Category"), my_model))

## marine

subdata <- bioshifts %>%
  filter(Eco=="Mar")

my_model <- glmmTMB(as.formula(model_formula),
                    # weight = calibration,
                    # family = Gamma(link = "log"),
                    data = subdata)

summary(my_model)

 Family: gaussian  ( identity )
Formula:          
log1p(abs(mismatch)) ~ Param + N_periodes + log(LatExtentk) +  
    Category + Grain_size + Obs_type + (1 | class)
Data: subdata

     AIC      BIC   logLik deviance df.resid 
  2083.7   2146.8  -1028.8   2057.7      938 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 class    (Intercept) 0.02984  0.1728  
 Residual             0.50106  0.7079  
Number of obs: 951, groups:  class, 15

Dispersion estimate for gaussian family (sigma^2): 0.501 

Conditional model:
                          Estimate Std. Error z value Pr(>|z|)    
(Intercept)             -2.061e+00  3.198e-01  -6.443 1.17e-10 ***
ParamO                  -1.936e-01  1.004e-01  -1.929  0.05372 .  
ParamTE                 -3.369e-01  1.056e-01  -3.192  0.00141 ** 
N_periodes              -4.968e-05  2.823e-03  -0.018  0.98596    
log(LatExtentk)          4.671e-01  4.574e-02  10.214  < 2e-16 ***
CategoryCensusPeriods   -1.119e-01  9.890e-02  -1.131  0.25806    
CategorySurvey-Resurvey -1.129e-02  2.183e-01  -0.052  0.95874    
Grain_sizemoderate       5.396e-01  1.648e-01   3.273  0.00106 ** 
Grain_sizelarge          1.272e-01  1.598e-01   0.796  0.42599    
Grain_sizevery_large     3.050e-01  7.179e-02   4.248 2.15e-05 ***
Obs_typeoccurrence       3.096e-01  1.132e-01   2.735  0.00624 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Actual    = log1p(abs(subdata$mismatch))
Predicted = predict(my_model, type="response")
Residuals = residuals(my_model)

efronRSquared(residual = Residuals, 
              predicted = Predicted, 
              statistic = "EfronRSquared")

EfronRSquared 
        0.267 

plot(effects::Effect(c("Category"), my_model))

Grain size moderatem large and very large have higher mismatches than small. Against our expectations, time-series have higher mismatches than census or survey-resurvey data.

10 3) The magnitude of noise species range shift?

## overall
model_formula <- c("log1p(abs(mismatch)) ~  bv.mat.error * Param + (1 | class) + (1 | Eco)")

my_model <- glmmTMB(as.formula(model_formula),
                    weight = calibration,
                    # family = Gamma(link = "log"),
                    data = bioshifts)

summary(my_model)

 Family: gaussian  ( identity )
Formula:          
log1p(abs(mismatch)) ~ bv.mat.error * Param + (1 | class) + (1 |      Eco)
Data: bioshifts
Weights: calibration

     AIC      BIC   logLik deviance df.resid 
  2093.4   2141.5  -1037.7   2075.4     1533 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 class    (Intercept) 0.09832  0.3136  
 Eco      (Intercept) 0.04544  0.2132  
 Residual             0.44444  0.6667  
Number of obs: 1542, groups:  class, 22; Eco, 2

Dispersion estimate for gaussian family (sigma^2): 0.444 

Conditional model:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)           0.91328    0.17918   5.097 3.45e-07 ***
bv.mat.error          0.18851    0.03915   4.815 1.47e-06 ***
ParamO               -0.22528    0.06950  -3.241 0.001190 ** 
ParamTE              -0.00526    0.13020  -0.040 0.967779    
bv.mat.error:ParamO   0.22986    0.06718   3.421 0.000623 ***
bv.mat.error:ParamTE -0.22588    0.11956  -1.889 0.058853 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Actual    = log1p(abs(bioshifts$mismatch))
Predicted = predict(my_model, type="response")
Residuals = residuals(my_model)

efronRSquared(residual = Residuals, 
              predicted = Predicted, 
              statistic = "EfronRSquared")

EfronRSquared 
        0.282 

plot(effects::Effect(c("bv.mat.error","Param"), my_model))

## terrestrial
model_formula <- c("log1p(abs(mismatch)) ~  bv.mat.error * Param + (1 | class)")

subdata <- bioshifts %>%
  filter(Eco=="Ter")

my_model <- glmmTMB(as.formula(model_formula),
                    weight = calibration,
                    # family = Gamma(link = "log"),
                    data = subdata)

summary(my_model)

 Family: gaussian  ( identity )
Formula:          log1p(abs(mismatch)) ~ bv.mat.error * Param + (1 | class)
Data: subdata
Weights: calibration

     AIC      BIC   logLik deviance df.resid 
   559.3    594.3   -271.6    543.3      584 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 class    (Intercept) 0.02616  0.1617  
 Residual             0.24986  0.4999  
Number of obs: 592, groups:  class, 7

Dispersion estimate for gaussian family (sigma^2): 0.25 

Conditional model:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)           0.65759    0.09427   6.976 3.04e-12 ***
bv.mat.error          0.18311    0.03979   4.602 4.19e-06 ***
ParamO               -0.20485    0.07300  -2.806  0.00501 ** 
ParamTE              -0.15769    0.16844  -0.936  0.34917    
bv.mat.error:ParamO   0.16916    0.13613   1.243  0.21399    
bv.mat.error:ParamTE  0.74618    0.49771   1.499  0.13382    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Actual    = log1p(abs(subdata$mismatch))
Predicted = predict(my_model, type="response")
Residuals = residuals(my_model)

efronRSquared(residual = Residuals, 
              predicted = Predicted, 
              statistic = "EfronRSquared")

EfronRSquared 
         0.24 

plot(effects::Effect(c("bv.mat.error","Param"), my_model))

## marine
model_formula <- c("log1p(abs(mismatch)) ~  bv.mat.error + Param + (1 | class)")

subdata <- bioshifts %>%
  filter(Eco=="Mar")

my_model <- glmmTMB(as.formula(model_formula),
                    weight = calibration,
                    # family = Gamma(link = "log"),
                    data = subdata)

summary(my_model)

 Family: gaussian  ( identity )
Formula:          log1p(abs(mismatch)) ~ bv.mat.error + Param + (1 | class)
Data: subdata
Weights: calibration

     AIC      BIC   logLik deviance df.resid 
  1479.8   1509.0   -733.9   1467.8      944 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 class    (Intercept) 0.1101   0.3318  
 Residual             0.5675   0.7533  
Number of obs: 950, groups:  class, 15

Dispersion estimate for gaussian family (sigma^2): 0.567 

Conditional model:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)   1.09459    0.11894   9.203  < 2e-16 ***
bv.mat.error  0.27314    0.04525   6.037 1.57e-09 ***
ParamO       -0.08630    0.09050  -0.954    0.340    
ParamTE      -0.21766    0.13270  -1.640    0.101    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Actual    = log1p(abs(subdata$mismatch))
Predicted = predict(my_model, type="response")
Residuals = residuals(my_model)

efronRSquared(residual = Residuals, 
              predicted = Predicted, 
              statistic = "EfronRSquared")

EfronRSquared 
        0.152 

plot(effects::Effect(c("bv.mat.error","Param"), my_model))

According to expectatios, the higher noise the higher the mismatch between predictions and observations.

11 4) Ladscape connectivity

model_formula <- c("log1p(abs(mismatch)) ~  impeded_prop * Param + N_periodes + log(LatExtentk) + Category + Grain_size + Obs_type + (1 | class) + (1 | Eco)")
# model_formula <- c("ShiftRate ~  bv.mat2 * Param + (1 | class) + (1 | Eco)")

my_model <- glmmTMB(as.formula(model_formula),
                    # weight = calibration,
                    data = bioshifts)

summary(my_model)

 Family: gaussian  ( identity )
Formula:          
log1p(abs(mismatch)) ~ impeded_prop * Param + N_periodes + log(LatExtentk) +  
    Category + Grain_size + Obs_type + (1 | class) + (1 | Eco)
Data: bioshifts

     AIC      BIC   logLik deviance df.resid 
   680.0    746.7   -323.0    646.0      358 

Random effects:

Conditional model:
 Groups   Name        Variance  Std.Dev. 
 class    (Intercept) 4.210e-02 2.052e-01
 Eco      (Intercept) 5.674e-10 2.382e-05
 Residual             3.176e-01 5.636e-01
Number of obs: 375, groups:  class, 7; Eco, 1

Dispersion estimate for gaussian family (sigma^2): 0.318 

Conditional model:
                         Estimate Std. Error z value Pr(>|z|)   
(Intercept)             -1.133205   0.796181  -1.423  0.15465   
impeded_prop             1.369918   0.432682   3.166  0.00154 **
ParamO                   0.208051   0.284208   0.732  0.46414   
ParamTE                  3.127768   2.059985   1.518  0.12893   
N_periodes               0.003121   0.034322   0.091  0.92755   
log(LatExtentk)          0.225363   0.080035   2.816  0.00487 **
CategoryCensusPeriods   -0.105220   0.340340  -0.309  0.75720   
CategorySurvey-Resurvey -0.392925   0.547173  -0.718  0.47269   
Grain_sizemoderate       0.130111   0.175111   0.743  0.45747   
Grain_sizelarge         -0.087667   0.152665  -0.574  0.56580   
Grain_sizevery_large    -0.338662   0.306487  -1.105  0.26917   
Obs_typeoccurrence      -0.125314   0.213042  -0.588  0.55639   
impeded_prop:ParamO     -0.875220   0.466223  -1.877  0.06048 . 
impeded_prop:ParamTE    -6.507612   3.813513  -1.706  0.08792 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Actual    = bioshifts$ShiftRate
Predicted = predict(my_model, type="response")
Residuals = residuals(my_model)

efronRSquared(residual = Residuals, 
              predicted = Predicted, 
              statistic = "EfronRSquared")

EfronRSquared 
        0.296 

plot(effects::Effect(c("impeded_prop","Param"), my_model))

## terrestrial
model_formula <- c("log1p(abs(mismatch)) ~  impeded_prop * Param + N_periodes + log(LatExtentk) + Category + Grain_size + Obs_type + (1 | class)")

subdata <- bioshifts %>%
  filter(Eco=="Ter")

my_model <- glmmTMB(as.formula(model_formula),
                    # weight = calibration,
                    data = subdata)

summary(my_model)

 Family: gaussian  ( identity )
Formula:          
log1p(abs(mismatch)) ~ impeded_prop * Param + N_periodes + log(LatExtentk) +  
    Category + Grain_size + Obs_type + (1 | class)
Data: subdata

     AIC      BIC   logLik deviance df.resid 
   678.0    740.8   -323.0    646.0      359 

Random effects:

Conditional model:
 Groups   Name        Variance Std.Dev.
 class    (Intercept) 0.0421   0.2052  
 Residual             0.3176   0.5636  
Number of obs: 375, groups:  class, 7

Dispersion estimate for gaussian family (sigma^2): 0.318 

Conditional model:
                         Estimate Std. Error z value Pr(>|z|)   
(Intercept)             -1.133208   0.796181  -1.423  0.15465   
impeded_prop             1.369921   0.432682   3.166  0.00154 **
ParamO                   0.208054   0.284208   0.732  0.46414   
ParamTE                  3.127784   2.059985   1.518  0.12893   
N_periodes               0.003121   0.034322   0.091  0.92754   
log(LatExtentk)          0.225362   0.080035   2.816  0.00487 **
CategoryCensusPeriods   -0.105218   0.340340  -0.309  0.75720   
CategorySurvey-Resurvey -0.392918   0.547173  -0.718  0.47270   
Grain_sizemoderate       0.130110   0.175111   0.743  0.45747   
Grain_sizelarge         -0.087667   0.152665  -0.574  0.56580   
Grain_sizevery_large    -0.338667   0.306486  -1.105  0.26916   
Obs_typeoccurrence      -0.125312   0.213042  -0.588  0.55640   
impeded_prop:ParamO     -0.875221   0.466223  -1.877  0.06048 . 
impeded_prop:ParamTE    -6.507643   3.813513  -1.706  0.08792 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Actual    = log1p(abs(subdata$mismatch))
Predicted = predict(my_model, type="response")
Residuals = residuals(my_model)

efronRSquared(residual = Residuals, 
              predicted = Predicted, 
              statistic = "EfronRSquared")

EfronRSquared 
        0.296 

plot(effects::Effect(c("impeded_prop","Param"), my_model))

## marine

# subdata <- bioshifts %>%
#   filter(Eco=="Mar")
# 
# my_model <- glmmTMB(as.formula(model_formula),
#                     # weight = calibration,
#                     # family = Gamma(link = "log"),
#                     data = subdata)
# 
# summary(my_model)
# Actual    = log1p(abs(subdata$mismatch))
# Predicted = predict(my_model, type="response")
# Residuals = residuals(my_model)
# 
# efronRSquared(residual = Residuals, 
#               predicted = Predicted, 
#               statistic = "EfronRSquared")
# 
# plot(effects::Effect(c("impeded_prop","Param"), my_model))

# cannot fit for marine species
# Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) : 
#   contrasts can be applied only to factors with 2 or more levels

As higher the proportion of impeded cells, higher is the mismatch between observations and predictions.