Exploratory analyses

Code

rm(list=ls())
gc()

          used (Mb) gc trigger (Mb) max used (Mb)
Ncells  587057 31.4    1346851   72   660908 35.3
Vcells 1066839  8.2    8388608   64  1800310 13.8

Code

library(dplyr)
library(tidyverse)
library(tidyr)
library(ggplot2)
library(knitr)
library(pbapply)
library(data.table)

1 Source code

Code

try(source(here::here("R/settings.R")))

Error in eval(ei, envir) : object 'realm' not found

2 Load data

Code

realm <- "Mar"

# Bioshifts data
sp_bioshifts <- read.csv(here::here("Output",paste(realm,"sps_bioshifts.csv")))
# SDM outputs data
sp_SDM_out <- read.csv(here::here("Output",paste(realm,"sp_SDM_out.csv")))
# range position
sp_range_pos <- read.csv(here::here("Output",paste(realm,"sp_range_pos.csv")))
sp_range_pos <- sp_range_pos[,-which(names(sp_range_pos) %in% c("START","END", "BG"))]
# Exposure
sp_exposure <- read.csv(here::here("Output",paste(realm,"sp_exposure.csv")))
sp_exposure_SA <- sp_exposure[,-which(names(sp_exposure) == "exp_BG")]
sp_exposure_BG <- sp_exposure[,-which(names(sp_exposure) == "exp_SA")]
sp_exposure_BG <- sp_exposure_BG %>% spread(bio, exp_BG)
names(sp_exposure_BG)[6:9] <- paste(names(sp_exposure_BG)[6:9],
                                    "BG", sep = "_")
sp_exposure_SA <- sp_exposure_SA %>% spread(bio, exp_SA)
names(sp_exposure_SA)[6:9] <- paste(names(sp_exposure_SA)[6:9],
                                    "SA", sep = "_")
sp_exposure <- cbind(sp_exposure_BG, sp_exposure_SA[,6:9])
sp_exposure <- sp_exposure[,-which(names(sp_exposure) %in% c("START","END"))]

# SDM shifts
sp_SDM_shift <- read.csv(here::here("Output",paste(realm,"sp_SDM_shift.csv")))
sp_SDM_shift <- sp_SDM_shift %>% filter(!quant %in% c("0%","100%","mean"))
sp_SDM_shift$time_period <- paste(round(sp_SDM_shift$START),
                                  round(sp_SDM_shift$END),
                                  sep = "-")

3 Merge SDM results with bioshifts

3.1 Leading edge

Code

LE <- sp_bioshifts %>% filter(Param == "LE")

# add SDM output data
sp_SDM_out_CV <- sp_SDM_out %>% filter(Mod == "CV")
LE <- left_join(LE, sp_SDM_out_CV, by = join_by(New_name == Species))
LE$time_period <- paste(round(LE$START), round(LE$END),sep="-")

# add range position
LE <- left_join(LE, sp_range_pos, 
                by = join_by(New_name == Species,
                             ID == ID, 
                             time_period == time_period,
                             Type == type))

# add Exposure
LE <- left_join(LE, sp_exposure, 
                by = join_by(New_name == Species,
                             ID == ID, 
                             time_period == Time_periods))

# add SDM shift
tmp <- sp_SDM_shift %>% 
    filter(quant %in% c("5%","25%","50%","75%","95%","99%")) %>%
    select(c("quant","shift_BG","shift_SA","ID","Species", "Type", "time_period"))
tmp$quant <- gsub("%","",tmp$quant)
tmp_BG <- tmp[,-3]
tmp_BG <- tmp_BG %>% spread(quant, value = shift_BG)
names(tmp_BG)[5:9] <- paste0("BG_", names(tmp_BG)[5:9])
tmp_SA <- tmp[,-2]
tmp_SA <- tmp_SA %>% spread(quant, value = shift_SA)
names(tmp_SA)[5:9] <- paste0("SA_", names(tmp_SA)[5:9])
tmp <- cbind(tmp_BG,tmp_SA[,5:9])

LE <- left_join(LE, tmp, 
                by = join_by(New_name == Species,
                             ID == ID, 
                             time_period == time_period,
                             Type == Type))

3.2 Trailing edge

Code

TE <- sp_bioshifts %>% filter(Param == "TE")

# add SDM output data
sp_SDM_out_CV <- sp_SDM_out %>% filter(Mod == "CV")
TE <- left_join(TE, sp_SDM_out_CV, by = join_by(New_name == Species))
TE$time_period <- paste(round(TE$START), round(TE$END),sep="-")

# add range position
TE <- left_join(TE, sp_range_pos, 
                by = join_by(New_name == Species,
                             ID == ID, 
                             time_period == time_period,
                             Type == type))

# add Exposure
TE <- left_join(TE, sp_exposure, 
                by = join_by(New_name == Species,
                             ID == ID, 
                             time_period == Time_periods))

# add SDM shift
tmp <- sp_SDM_shift %>% 
    filter(quant %in% c("5%","25%","50%","75%","95%","99%")) %>%
    select(c("quant","shift_BG","shift_SA","ID","Species", "Type", "time_period"))
tmp$quant <- gsub("%","",tmp$quant)
tmp_BG <- tmp[,-3]
tmp_BG <- tmp_BG %>% spread(quant, value = shift_BG)
names(tmp_BG)[5:9] <- paste0("BG_", names(tmp_BG)[5:9])
tmp_SA <- tmp[,-2]
tmp_SA <- tmp_SA %>% spread(quant, value = shift_SA)
names(tmp_SA)[5:9] <- paste0("SA_", names(tmp_SA)[5:9])
tmp <- cbind(tmp_BG,tmp_SA[,5:9])

TE <- left_join(TE, tmp, 
                by = join_by(New_name == Species,
                             ID == ID, 
                             time_period == time_period,
                             Type == Type))

3.3 Optimum

Code

O <- sp_bioshifts %>% filter(Param == "O")

# add SDM output data
sp_SDM_out_CV <- sp_SDM_out %>% filter(Mod == "CV")
O <- left_join(O, sp_SDM_out_CV, by = join_by(New_name == Species))
O$time_period <- paste(round(O$START), round(O$END),sep="-")

# add range position
O <- left_join(O, sp_range_pos, 
               by = join_by(New_name == Species,
                            ID == ID, 
                            time_period == time_period,
                            Type == type))

# add Exposure
O <- left_join(O, sp_exposure, 
               by = join_by(New_name == Species,
                            ID == ID, 
                            time_period == Time_periods))

# add SDM shift
tmp <- sp_SDM_shift %>% 
    filter(quant %in% c("5%","25%","50%","75%","95%","99%")) %>%
    select(c("quant","shift_BG","shift_SA","ID","Species", "Type", "time_period"))
tmp$quant <- gsub("%","",tmp$quant)
tmp_BG <- tmp[,-3]
tmp_BG <- tmp_BG %>% spread(quant, value = shift_BG)
names(tmp_BG)[5:9] <- paste0("BG_", names(tmp_BG)[5:9])
tmp_SA <- tmp[,-2]
tmp_SA <- tmp_SA %>% spread(quant, value = shift_SA)
names(tmp_SA)[5:9] <- paste0("SA_", names(tmp_SA)[5:9])
tmp <- cbind(tmp_BG,tmp_SA[,5:9])

O <- left_join(O, tmp, 
               by = join_by(New_name == Species,
                            ID == ID, 
                            time_period == time_period,
                            Type == Type))

There are 173 shifts at the leading edge, 66 shifts at the trailing edge, and 286 at the optimum (center) of the species range.

3.4 Group

Code

SDM_shifts <- rbind(LE,TE,O)

4 Explore relationships

4.1 N shifts per parameter

Code

tmp <- data.frame(table(SDM_shifts$Param))

ggplot(tmp, aes(x = Var1, y = Freq))+
    geom_col()+
    theme_classic()+
    geom_text(aes(label = Freq), vjust = -0.5)+
    xlab("N shifts")

5 How closely related are observed and potential shifts?

5.1 Shifts at the study area

Code

tmp_s_SA <-  SDM_shifts %>% 
    mutate(SA_5 = (SA_5*111)/DUR, 
           SA_25 = (SA_25*111)/DUR,
           SA_50 = (SA_50*111)/DUR, 
           SA_75 = (SA_75*111)/DUR, 
           SA_95 = (SA_95*111)/DUR) %>%
    select(SHIFT, Param, Param, New_name, Type, ID, time_period,
           SA_5, SA_25, SA_50, SA_75, SA_95) 

tmp_s_SA <- tmp_s_SA %>% 
    gather(type, SDM_shift, -c(SHIFT,Param,Param,New_name,Type,ID,time_period))  
tmp_s_SA$type <- factor(tmp_s_SA$type, 
                        levels = c("SA_5", "SA_25", "SA_50", "SA_75", "SA_95"))

ggplot(tmp_s_SA, aes(x = SHIFT, y = SDM_shift))+
    geom_point()+
    xlab("Observed shift (km/yr)")+
    ylab("Potential shift (km/yr)")+
    theme_classic() +
    geom_smooth(method = "gam")+
    facet_wrap(.~type, nrow = 1)

5.1.1 Per parameter

Code

ggplot(tmp_s_SA, aes(x = SHIFT, y = SDM_shift))+
    geom_point()+
    xlab("Observed shift (km/yr)")+
    ylab("Potential shift (km/yr)")+
    theme_classic() +
    geom_smooth(method = "gam")+
    facet_wrap(Param~type, nrow = 3)

5.2 Shifts at the background area

Code

tmp_s_BG <-  SDM_shifts %>% 
    mutate(BG_5 = BG_5/DUR, 
           BG_25 = BG_25/DUR,
           BG_50 = BG_50/DUR, 
           BG_75 = BG_75/DUR, 
           BG_95 = BG_95/DUR) %>%
    select(SHIFT, Param, Param, New_name, Type, ID, time_period,
           BG_5, BG_25, BG_50, BG_75, BG_95) 

tmp_s_BG <- tmp_s_BG %>% 
    gather(type, SDM_shift, -c(SHIFT,Param,Param,New_name,Type,ID,time_period))  
tmp_s_BG$type <- factor(tmp_s_BG$type, 
                        levels = c("BG_5", "BG_25", "BG_50", "BG_75", "BG_95"))

ggplot(tmp_s_BG, aes(x = SHIFT, y = SDM_shift))+
    geom_point()+
    xlab("Observed shift (km/yr)")+
    ylab("Potential shift (km/yr)")+
    theme_classic() +
    geom_smooth(method = "gam")+
    facet_wrap(.~type, nrow = 1)

5.2.1 Per parameter

Code

ggplot(tmp_s_BG, aes(x = SHIFT, y = SDM_shift))+
    geom_point()+
    xlab("Observed shift (km/yr)")+
    ylab("Potential shift (km/yr)")+
    theme_classic() +
    geom_smooth(method = "gam")+
    facet_wrap(Param~type, nrow = 3)

Observed shifts are unrelated to shifts at the background area. This suggests that observed shifts represent local responses rather than a response that could be extrapolated to the complete species range.

6 Is range shift sensitive to climate exposure?

6.1 Exposure at the study area

Code

tmp_e_SA <-  SDM_shifts %>% 
    select(SHIFT, Param, New_name, Type, ID, time_period,
           SST_max_SA, SST_mean_SA, SST_min_SA, SST_sd_SA) 

tmp_e_SA <- tmp_e_SA %>% 
    gather(type, Exposure, -c(SHIFT,Param,New_name,Type,ID,time_period))  
tmp_e_SA$type <- factor(tmp_e_SA$type, 
                        levels = c("SST_max_SA", "SST_mean_SA", "SST_min_SA", "SST_sd_SA"))

ggplot(tmp_e_SA, aes(x = Exposure, y = SHIFT))+
    geom_point()+
    xlab("Exposure")+
    ylab("Observed shift (km/yr)")+
    theme_classic() +
    geom_smooth(method = "gam")+
    facet_wrap(.~type, nrow = 1)

6.1.1 Per parameter

Code

ggplot(tmp_e_SA, aes(x = Exposure, y = SHIFT))+
    geom_point()+
    xlab("Exposure")+
    ylab("Observed shift (km/yr)")+
    theme_classic() +
    geom_smooth(method = "gam")+
    facet_wrap(Param~type, nrow = 3)

6.2 Exposure at the background area

Code

tmp_e_BG <-  SDM_shifts %>% 
    select(SHIFT, Param, New_name, Type, ID, time_period,
           SST_max_BG, SST_mean_BG, SST_min_BG, SST_sd_BG) 

tmp_e_BG <- tmp_e_BG %>% 
    gather(type, Exposure, -c(SHIFT,Param,New_name,Type,ID,time_period))  
tmp_e_BG$type <- factor(tmp_e_BG$type, 
                        levels = c("SST_max_BG", "SST_mean_BG", "SST_min_BG", "SST_sd_BG"))

ggplot(tmp_e_BG, aes(x = Exposure, y = SHIFT))+
    geom_point()+
    xlab("Exposure")+
    ylab("Observed shift (km/yr)")+
    theme_classic() +
    geom_smooth(method = "gam")+
    facet_wrap(.~type, nrow = 1)

6.2.1 Per parameter

Code

ggplot(tmp_e_BG, aes(x = Exposure, y = SHIFT))+
    geom_point()+
    xlab("Exposure")+
    ylab("Observed shift (km/yr)")+
    theme_classic() +
    geom_smooth(method = "gam")+
    facet_wrap(Param~type, nrow = 3)

`geom_smooth()` using formula 'y ~ s(x, bs = "cs")'

Warning: Removed 558 rows containing non-finite values (stat_smooth).

Warning: Removed 558 rows containing missing values (geom_point).

7 Are lags in range shift related to exposure?

7.1 Lags at the study area

Code

tmp_s_SA$lag_SA <- tmp_s_SA$SHIFT/tmp_s_SA$SDM_shift
names(tmp_s_SA)[which(names(tmp_s_SA)=="type")] <- "type_shift"
tmp_e_SA2 <- tmp_e_SA[,-which(names(tmp_e_SA)=="SHIFT")]
names(tmp_e_SA2)[which(names(tmp_e_SA2)=="type")] <- "type_exposure"

tmp_lags <- left_join(tmp_s_SA,tmp_e_SA2, 
                      by = join_by(New_name == New_name,
                                   Param == Param,
                                   ID == ID, 
                                   time_period == time_period,
                                   Type == Type))

7.2 Distribution of lag values

Code

tmp_lags2 <- tmp_lags[-which(is.infinite(tmp_lags$lag_SA)),]
tmp_lags2 <- na.omit(tmp_lags2)

ggplot(tmp_lags2, 
       aes(x=Param, y=lag_SA, 
           fill = Param, color = Param))+
    geom_jitter(alpha = .3, 
                position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), 
                aes(color = Param))+
    geom_boxplot(alpha = .3)+
    theme_classic() +
    facet_wrap(.~type_shift)

Code

my_aov <- aov(lag_SA~Param, 
              tmp_lags2)
summary(my_aov)

              Df    Sum Sq   Mean Sq F value Pr(>F)
Param          2 6.998e+27 3.499e+27   1.165  0.312
Residuals   5928 1.780e+31 3.003e+27               

Code

TukeyHSD(my_aov)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = lag_SA ~ Param, data = tmp_lags2)

$Param
              diff           lwr          upr     p adj
O-LE  2.243674e+12 -1.636815e+12 6.124163e+12 0.3646132
TE-LE 2.793403e+12 -2.281384e+12 7.868190e+12 0.4006194
TE-O  5.497292e+11 -4.024667e+12 5.124126e+12 0.9571923

Remove outliers ::: {.cell}

Code

quartiles <- quantile(tmp_lags2$lag_SA, probs=c(.25, .75))
IQR <- IQR(tmp_lags2$lag_SA)
 
Lower <- quartiles[1] - 1.5*IQR
Upper <- quartiles[2] + 1.5*IQR 
 
tmp_lags2 <- subset(tmp_lags2, tmp_lags2$lag_SA > Lower & tmp_lags2$lag_SA < Upper)

ggplot(tmp_lags2, 
       aes(x=Param, y=lag_SA, 
           fill = Param, color = Param))+
    geom_jitter(alpha = .3, 
                position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8), 
                aes(color = Param))+
    geom_boxplot(alpha = .3)+
    theme_classic() +
    facet_wrap(.~type_shift)

Code

my_aov <- aov(lag_SA~Param, 
              tmp_lags2)
summary(my_aov)

              Df Sum Sq Mean Sq F value Pr(>F)    
Param          2   5190  2594.9   62.52 <2e-16 ***
Residuals   4647 192867    41.5                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

TukeyHSD(my_aov)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = lag_SA ~ Param, data = tmp_lags2)

$Param
           diff        lwr        upr   p adj
O-LE   1.136053  0.6044947  1.6676120 1.7e-06
TE-LE -1.567521 -2.2232656 -0.9117769 1.0e-07
TE-O  -2.703575 -3.2778555 -2.1292937 0.0e+00

:::

Code

ggplot(tmp_lags2, aes(x = Exposure, y = lag_SA))+
    geom_point()+
    xlab("Exposure")+
    ylab("Lag")+
    theme_classic() +
    geom_smooth(method = "gam")+
    facet_wrap(.~type_exposure, nrow = 1)

Code

ggplot(tmp_lags2, aes(x = Exposure, y = lag_SA))+
    geom_point()+
    xlab("Exposure")+
    ylab("Lag")+
    theme_classic() +
    geom_smooth(method = "gam")+
    facet_wrap(type_exposure~type_shift, nrow = 4)