Adaptive potential models
Bioshifts group
1 Overview
Genetic diversity (GD) is an important component of biodiversity because it provides adaptive potential (Frankham, 2005), facilitates species persistence (Frankham, 2005; Pereira et al., 2012), and provides clues about the historical demography of species (Avise, 2000). Low GD (due to inbreeding or genetic drift) likely threatens species against environmental change or fluctuation and could indicate areas were populations are going extinct. Here we test the potential of GD in modulating species responses to climate exposure by shifting their ranges. Ongoing climate change has triggered a massive redistribution of species on Earth (Parmesan & Yole, 2003; Lenoir et al. 2020), with far reaching consequences on ecosystems and human well-being (Scheffers et al. 2016). These range shifts are highly idiosyncratic, likely reflecting differences in climate exposure but also species sensitivity and capacity to adapt to environmental changes. Yet, we still lack a clear understanding of the mechanisms of species response to climate exposure across taxonomic groups and habitats.
We contrast two hypotheses by which adaptive potential from GD modulates species sensitivity to climate exposure. The first hypothesis proposes GD allows species to better track climate change. In this scenario, GD confers species’ with ability to rapidly colonize new areas that become suitable as climate change. This may occur if adaptive potential is linked to ecological generalism, competition advantages or plasticity when colonizing new areas. In contrast, low GD may promote species range shifts driven by extirpation. The second hypothesis states that GD allows species to stay in place as climate change. In this scenario GD increases plasticity or allow fast evolution of traits for tolerating new climate conditions, so that species do not have to move for tracking climate change.
We test if these mechanisms hold across different parts of a species range, which may provide deeper insights into the processes by which GD modulates species responses to climate change. At the trailing edge (or warm edge) species are typically closer to their upper thermal limit. If GD helps species to track climate change at the trailing edge, it might minimize rates of extinction and range contractions from climate change. In contrast, if GD helps species stay in place as climate change at the trailing edge, it indicates GD facilitates evolution or plasticity of thermal tolerance. At the leading edge (or cold edge), climate change often implies expansion of suitable area. If GD helps species to track climate change at the leading edge, it might facilitate colonization and range expansions. In contrast, if GD helps species stay in place as climate change at the leading edge, this indicates a decoupling between GD and colonization rates, or at least that the effect of GD on colonization rate is not enough for species to keep up with rates of climate change.
2 Load in data
2.1 Load libraries
# remotes::install_github("bbolker/broom.mixed")
rm(list=ls())
list.of.packages <- c(
"ggplot2", "data.table","pdftools","plotly","kableExtra",
"pbapply","dplyr", "tidyr", "foreach", "parallel", "doParallel",
"scales","effects","psych", "glmmTMB", "lme4", "lmerTest",
"here","rlist","ggtext","gridExtra","grid","lattice","viridis", "patchwork", "purrr", "DT")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
sapply(list.of.packages, require, character.only = TRUE)## ggplot2 data.table pdftools plotly kableExtra pbapply dplyr
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## tidyr foreach parallel doParallel scales effects psych
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## glmmTMB lme4 lmerTest here rlist ggtext gridExtra
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## grid lattice viridis patchwork purrr DT
## TRUE TRUE TRUE TRUE TRUE TRUE2.2 Load data
mydataset <- read.csv(here::here("Data/gen_data_final_fonseca.csv"))2.3 Latitude data
mydatatogo <- mydataset %>%
dplyr::filter(Type == "LAT",
shift_vel_sign == "pospos" | shift_vel_sign == "negneg", # Select only shifts in the same direction of velocity
SHIFT != 0, # remove non-significant shifts
# Nucleotide_diversity > 0 # select only GD values > 0
) %>%
mutate(
# Climate velocity
vel = as.numeric(velocity),
vel_abs = abs(vel),
vel_abs_log = log(vel_abs),
# Shift
SHIFT = SHIFT, # to deal with zero shift
SHIFT_abs = abs(SHIFT),
SHIFT_abs_log = log(SHIFT_abs),
# Genetic diversity
GD = Nucleotide_diversity,
LatGen = abs(Latitude),
# Methods
Lat = LatCentDeg,
Lon = LonCentDeg,
ID.area = scale(ID.area),
DUR = scale(DUR),
LogExtent = log(Extent),
START = scale(START),
Param = factor(Param),
Group = factor(Group),
spp = factor(spp),
ExtentF = cut(Extent,
ordered_result = TRUE,
breaks=seq(min(Extent), max(Extent), length.out=10),
include.lowest=TRUE),
NtempUnitsF = cut(NtempUnits,
ordered_result = TRUE,
breaks=seq(min(NtempUnits), max(NtempUnits), length.out=10),
include.lowest=TRUE),
Continent = case_when(
Lon < -20 ~ "NorthAmerica",
Lon > -20 & Lon < 40 & Lat > 30 ~ "Europe",
Lon > -20 & Lon < 50 & Lat < 30 ~ "African",
Lon > 50 & Lat < 0 ~ "Australian",
Lon > 50 & Lat > 0 ~ "Asian")) %>%
dplyr::select(
# Genetic diversity
GD, TajimasD,
LatGen, # latitudinal position where GD was collected
# Shift
SHIFT, SHIFT_abs, SHIFT_abs_log,
# SHIFT_cor, SHIFT_cor_abs, SHIFT_cor_abs_log, SHIFT_cor_raw, SHIFT_abs_log_scale,
# Velocity
vel, vel_abs, vel_abs_log,
trend.mean,
# Methods + Taxonomy
Article_ID,
shift_vel_sign,
Lat,Lon,Continent,
DUR, Nperiodes, LogNtempUnits, NtempUnits, LogExtent, ContinuousGrain, Quality, PrAb, ExtentF, NtempUnitsF,
Param, Group, spp, Class, Order, Family, Genus,
ECO, Grain_size, Data, Article_ID,
# Distance from GD data to Bioshifts SA
dist_cent, # distance to the centroid of the SA
dist_edge # distance to the closest edge of the SA
)
# transform continuous variables
cont_vars <- c(1:10, 13,14, 16:21, 36, 37)
mydatatogo[,cont_vars] <- lapply(mydatatogo[,cont_vars], as.numeric)
mydatatogo[,-cont_vars] <- lapply(mydatatogo[,-cont_vars], function(x) factor(x, levels = unique(x)))
# remove two range shift outliers
# summary(mydatatogo$SHIFT_abs)
mydatatogo <- mydatatogo %>% filter(SHIFT_abs < 100)2.4 Filter Classes with at least 5 species per Param
n_sps = 10
test <- mydatatogo %>%
group_by(Class,Param) %>%
dplyr::summarise(N_spp = length(unique(spp))) %>% # how many species per parameter?
dplyr::filter(N_spp >= n_sps) # select classes with > n_sps per param## `summarise()` has grouped output by 'Class'. You can override using the
## `.groups` argument.test <- mutate(test, Class_Param = paste(Class, Param))
mydatatogo <- mydatatogo %>%
mutate(Class_Param = paste(Class, Param)) %>%
filter(Class_Param %in% test$Class_Param) %>%
select(-Class_Param)
mydatatogo[,-cont_vars] <- lapply(mydatatogo[,-cont_vars], function(x) factor(x, levels = unique(x)))2.5 Extra fixes
# Fix param levels
mydatatogo$Param <- relevel(mydatatogo$Param, ref = "O")
# Calculate N obs per species
N_obs_spp <- mydatatogo %>%
group_by(spp) %>%
tally()
N_obs_spp$weight_obs_spp <- (1/N_obs_spp$n)*(1/nrow(N_obs_spp))
mydatatogo <- merge(mydatatogo,N_obs_spp,by="spp")
sum(mydatatogo$weight_obs_spp)## [1] 13 Data description
Genetic diversity (GD) data was extracted from Fonseca et al. (2023),
describing patterns of intraspecific nucleotide diversity. Nucleotide
diversity was calculated from aggregated georeferenced mitochondrial DNA
for animals (vertebrates, arachnids, and insects) and chloroplast DNA
for plants using phylogatR (Pelletier et al., 2021; https://phylogatr.osc.edu), an open-source platform
available through the Ohio Supercomputer Center (OSC).
phylogatR retrieves genetic and geographic data by
integrating information from the GBIF (https://www.gbif.org/), GenBank (https://www.ncbi.nlm.nih.gov/genbank/), and the BOLD (http://www.boldsystems.org/index.php). To minimize data
loss and duplication when downloading data from GenBank, duplicate
entries with the same accession number, but different GBIF occurrence
numbers are checked for similarities in geographic coordinates.
Sequences with high values of divergence were removed (likely resulting
from inclusion of a sequence with an incorrect species label), a
conservative approach that prevents inflation on the calculation of
genetic diversity within species. Nucleotide diversity was calculated as
the number of polymorphic sites in pairwise sequence comparison, using a
modified version of the equation proposed by Nei and Li (1979) to allow
for a comparison of sequences of different lengths (Miraldo et
al. 2016). Finally, nucleotide diversity between pairs of sequences were
averaged to obtain a global estimation of nucleotide diversity for each
species, and associated centroid of georeferenced DNA sequences.
We gather range shift estimates from the BioShifts database (Comte et al. 2020). These shift estimates were collated from multiple published literature sources that used different methodological protocols. Previous studies have shown that methodological aspects play an important role in explaining variability in species range shifts (Lenoir et al. 2020). Beyond range shift data, we also gathered from BioShifts climate velocity estimates following the method from Burrows at al. (2014), the position of a species range where range shifts were estimated (leading edge, centroid or trailing edge), and methodological variables representing the number of temporal units (continuous), grain size (small, moderate or large), data type (abundance-based vs occurrence-based), quality (balanced, raw, resurveyed), and study area extent (log km). We focus on latitudinal shifts across terrestrial and marine organisms, excluding elevational range shifts. We also excluded range shifts at the opposite direction of the climate exposure gradient, which can be related to local extinction processes or because species may be tracking climate axes other then temperature (e.g., humidity).
In order to maximize the match of species when merging range shifts
and GD datasets, we harmonize species names using a customized function
that borrows from the R packages rgbif and
bdc. Our final dataset includes 4688 shifts for 1894
species from 7 major taxonomy classes (Insecta, Aves, Magnoliopsida,
Polypodiopsida, Arachnida, Liliopsida, Actinopterygii).
# N species
length(unique(mydatatogo$spp))## [1] 1894# 1890
# N species per Class
mydatatogo %>% group_by(Class) %>% summarise(n = length(unique(spp)))## # A tibble: 7 × 2
## Class n
## <fct> <int>
## 1 Insecta 1007
## 2 Aves 343
## 3 Magnoliopsida 194
## 4 Polypodiopsida 10
## 5 Arachnida 154
## 6 Liliopsida 52
## 7 Actinopterygii 134# Class n
# 1 Insecta 1004
# 2 Aves 342
# 3 Magnoliopsida 194
# 4 Polypodiopsida 10
# 5 Arachnida 154
# 6 Liliopsida 52
# 7 Actinopterygii 134
# N species per range position
mydatatogo %>% group_by(Param) %>% summarise(n = length(unique(spp)))## # A tibble: 3 × 2
## Param n
## <fct> <int>
## 1 O 898
## 2 LE 1352
## 3 TE 199# Param n
# 1 O 898
# 2 LE 1346
# 3 TE 199
# N shifts
nrow(mydatatogo)## [1] 4688# 4675
# N shifts per Class
mydatatogo %>% group_by(Class) %>% tally## # A tibble: 7 × 2
## Class n
## <fct> <int>
## 1 Insecta 1757
## 2 Aves 1618
## 3 Magnoliopsida 664
## 4 Polypodiopsida 23
## 5 Arachnida 248
## 6 Liliopsida 143
## 7 Actinopterygii 235# Class n
# 1 Insecta 1751
# 2 Aves 1611
# 3 Magnoliopsida 664
# 4 Polypodiopsida 23
# 5 Arachnida 248
# 6 Liliopsida 143
# 7 Actinopterygii 235
# N shifts per range position
table(mydatatogo$Param)##
## O LE TE
## 2303 2090 295# O LE TE
# 2298 2082 295
# N shifts per Class and range position
mydatatogo %>% group_by(Class, Param) %>% tally## # A tibble: 17 × 3
## # Groups: Class [7]
## Class Param n
## <fct> <fct> <int>
## 1 Insecta O 194
## 2 Insecta LE 1442
## 3 Insecta TE 121
## 4 Aves O 1281
## 5 Aves LE 261
## 6 Aves TE 76
## 7 Magnoliopsida O 496
## 8 Magnoliopsida LE 100
## 9 Magnoliopsida TE 68
## 10 Polypodiopsida O 23
## 11 Arachnida LE 248
## 12 Liliopsida O 108
## 13 Liliopsida LE 23
## 14 Liliopsida TE 12
## 15 Actinopterygii O 201
## 16 Actinopterygii LE 16
## 17 Actinopterygii TE 18# Class Param n
# 1 Insecta O 195
# 2 Insecta LE 1435
# 3 Insecta TE 121
# 4 Aves O 1275
# 5 Aves LE 260
# 6 Aves TE 76
# 7 Magnoliopsida O 496
# 8 Magnoliopsida LE 100
# 9 Magnoliopsida TE 68
# 10 Polypodiopsida O 23
# 11 Arachnida LE 248
# 12 Liliopsida O 108
# 13 Liliopsida LE 23
# 14 Liliopsida TE 12
# 15 Actinopterygii O 201
# 16 Actinopterygii LE 16
# 17 Actinopterygii TE 18
# View(mydatatogo[,c(1,2,8,33:38)])
summary(mydatatogo$NtempUnits)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.000 2.000 2.000 4.885 2.000 61.000table(mydatatogo$Quality)##
## BALANCED RAW RESURVEYED
## 1874 1857 957table(mydatatogo$PrAb)##
## occurrence abundance
## 3574 1114table(mydatatogo$Order)##
## Lepidoptera Passeriformes Accipitriformes
## 910 1041 71
## Sapindales Orthoptera Asterales
## 4 17 83
## Coleoptera Charadriiformes Polypodiales
## 139 168 9
## Dipsacales Strigiformes Apiales
## 22 44 23
## Hemiptera Odonata Araneae
## 20 122 244
## Poales Anseriformes Lamiales
## 119 87 104
## Coraciiformes Galliformes Brassicales
## 4 43 19
## Clupeiformes Gruiformes Perciformes
## 16 17 48
## Hymenoptera Ericales Anguilliformes
## 401 19 3
## Gadiformes Apodiformes Pelecaniformes
## 57 4 20
## Caryophyllales Stomiiformes Pleuronectiformes
## 103 4 47
## Diptera Fagales Ophioglossales
## 143 27 1
## Psittaciformes Ranunculales Solanales
## 1 26 12
## Scorpaeniformes Tetraodontiformes Saxifragales
## 34 3 9
## Myrtales Piciformes Columbiformes
## 26 55 23
## Rosales Asparagales Cuculiformes
## 49 10 5
## Fabales Myctophiformes Blattodea
## 54 4 1
## Equisetales Malpighiales Falconiformes
## 12 32 23
## Gentianales Geraniales Beryciformes
## 25 16 1
## Boraginales Alismatales Plecoptera
## 1 12 4
## Osmeriformes Lophiiformes Malvales
## 5 5 6
## Stephanoberyciformes Opiliones Suliformes
## 2 4 1
## Caprimulgiformes Salmoniformes Osmundales
## 2 2 1
## Oxalidales Podicipediformes Batrachoidiformes
## 3 9 1
## Syngnathiformes Liliales Vitales
## 1 2 1
## Zeiformes
## 2c <- mydatatogo %>%
group_by(Class) %>%
summarise(`Range shifts` = length(Class),
`Species` = length(unique(spp)))
DT::datatable(c,
rownames = FALSE,
extensions = 'Buttons',
options = list(pageLength = 100,
dom = 'B',
buttons = c('excel'))) %>%
DT::formatStyle(columns = colnames(c), `fontSize` = '80%')c <- mydatatogo %>%
group_by(Class) %>%
mutate(Class=factor(Class, levels = sort(unique(as.character(mydatatogo$Class))))) %>%
summarise(`Avg. Range shift velocity (SD)` = paste0(round(mean(SHIFT),2), " (", round(sd(SHIFT),2), ")"),
`Avg. Genetic diversity (SD)` = paste0(formatC(mean(GD), format = "e", digits=2), " (", formatC(sd(GD), format = "e", digits=1), ")"),
`Avg. Climate change velocity (SD)` = paste0(round(mean(vel_abs),2), " (", round(sd(vel_abs),3), ")"))
DT::datatable(c,
rownames = FALSE,
extensions = 'Buttons',
options = list(pageLength = 100,
dom = 'B',
buttons = c('excel'))) %>%
DT::formatStyle(columns = colnames(c), `fontSize` = '80%')4 Spatial distribution of data
5 Data lookup
tabletoshow <- mydatatogo %>% mutate_if(is.numeric, ~round(., 3))
DT::datatable(tabletoshow,
rownames = FALSE,
extensions = 'Buttons',
options = list(scroller = TRUE,
pageLength = 500,
scrollX = TRUE,
scrollY=500,
buttons = c('excel'))) %>%
DT::formatStyle(columns = colnames(tabletoshow), `fontSize` = '80%')6 Species lookup
6.1 Insecta
fig <- plot_ly(mydatatogo %>%
filter(Class=="Insecta"),
x = ~GD, y = ~vel_abs, z = ~SHIFT_abs, color = ~Order,
text = ~paste('Species:', spp, '<br>Family:', Family, '<br>Order:', Order,'<br>ECO:', ECO)) %>%
add_markers() %>%
layout(scene = list(xaxis = list(title = 'GD'),
yaxis = list(title = 'Isotherm shift'),
zaxis = list(title = 'Range shift')))
fig## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors6.2 Aves
fig <- plot_ly(mydatatogo %>%
filter(Class=="Aves"),
x = ~GD, y = ~vel_abs, z = ~SHIFT_abs,
text = ~paste('Species:', spp, '<br>Family:', Family, '<br>Order:', Order,'<br>ECO:', ECO)) %>%
add_markers() %>%
layout(scene = list(xaxis = list(title = 'GD'),
yaxis = list(title = 'Isotherm shift'),
zaxis = list(title = 'Range shift')))
fig6.3 Plants
fig <- plot_ly(mydatatogo %>%
filter(Class=="Magnoliopsida" | Class=="Liliopsida"),
x = ~GD, y = ~vel_abs, z = ~SHIFT_abs, color = ~Class,
text = ~paste('Species:', spp, '<br>Family:', Family, '<br>Order:', Order,'<br>ECO:', ECO)) %>%
add_markers() %>%
layout(scene = list(xaxis = list(title = 'GD'),
yaxis = list(title = 'Isotherm shift'),
zaxis = list(title = 'Range shift')))
fig6.4 Arachnida
fig <- plot_ly(mydatatogo %>%
filter(Class=="Arachnida"),
x = ~GD, y = ~vel_abs, z = ~SHIFT_abs,
text = ~paste('Species:', spp, '<br>Family:', Family, '<br>Order:', Order,'<br>ECO:', ECO)) %>%
add_markers() %>%
layout(scene = list(xaxis = list(title = 'GD'),
yaxis = list(title = 'Isotherm shift'),
zaxis = list(title = 'Range shift')))
fig6.5 Fish
fig <- plot_ly(mydatatogo %>%
filter(Class=="Actinopterygii"),
x = ~GD, y = ~vel_abs, z = ~SHIFT_abs,
text = ~paste('Species:', spp, '<br>Family:', Family, '<br>Order:', Order,'<br>ECO:', ECO)) %>%
add_markers() %>%
layout(scene = list(xaxis = list(title = 'GD'),
yaxis = list(title = 'Isotherm shift'),
zaxis = list(title = 'Range shift')))
fig6.6 Example species
##############################
# species with low GD and low range shift
tmp <- mydatatogo %>%
filter(GD >0) %>%
filter(GD < quantile(GD, .1) & SHIFT_abs < quantile(SHIFT_abs, .1)) %>%
select(spp,Class,GD,SHIFT_abs) %>%
group_by(spp,Class) %>%
summarise(GD = round(mean(GD),5),
SHIFT_abs = round(mean(SHIFT_abs),3)) %>%
arrange(GD)
# these species
tmp## # A tibble: 33 × 4
## # Groups: spp [33]
## spp Class GD SHIFT_abs
## <fct> <fct> <dbl> <dbl>
## 1 Rumex obtusifolius Magnoliopsida 0.00002 0.061
## 2 Cosmia trapezina Insecta 0.00016 0.048
## 3 Platycheirus granditarsis Insecta 0.00017 0.161
## 4 Micranthes nivalis Magnoliopsida 0.00019 0.021
## 5 Alliaria petiolata Magnoliopsida 0.00023 0.1
## 6 Persicaria maculosa Magnoliopsida 0.00025 0.079
## 7 Hypericum perforatum Magnoliopsida 0.00027 0.04
## 8 Epilobium palustre Magnoliopsida 0.00028 0.104
## 9 Calidris canutus Aves 0.0003 0.067
## 10 Accipiter striatus Aves 0.00031 0.18
## # ℹ 23 more rows# from these classes
table(tmp$Class) # mostly plants##
## Insecta Aves Magnoliopsida Polypodiopsida Arachnida
## 3 7 20 0 0
## Liliopsida Actinopterygii
## 2 1# do for each class
my_classes <- unique(mydatatogo$Class)
for(i in 1:length(my_classes)){
tmp <- mydatatogo %>%
filter(GD >0) %>%
filter(Class == my_classes[i]) %>%
filter(GD < quantile(GD, .1) & vel_abs < quantile(vel_abs, .1)) %>%
select(spp,Class) %>%
unique()
print(tmp)
}## spp Class
## 1 Araschnia levana Insecta
## 2 Argynnis paphia Insecta
## 3 Bombus hortorum Insecta
## 5 Bombus huntii Insecta
## 6 Bombus terricola Insecta
## 7 Bombus wurflenii Insecta
## 9 Cordulegaster boltonii Insecta
## 10 Libellula depressa Insecta
## spp Class
## 1 Acrocephalus schoenobaenus Aves
## 2 Amaurornis phoenicurus Aves
## 4 Aquila chrysaetos Aves
## 5 Ardea cinerea Aves
## 6 Chloris chloris Aves
## 8 Gallinago gallinago Aves
## 9 Loxia pytyopsittacus Aves
## 10 Numenius arquata Aves
## 13 Pipilo maculatus Aves
## 14 Podiceps cristatus Aves
## 15 Vanellus vanellus Aves
## [1] spp Class
## <0 rows> (or 0-length row.names)
## [1] spp Class
## <0 rows> (or 0-length row.names)
## [1] spp Class
## <0 rows> (or 0-length row.names)
## spp Class
## 1 Poa glauca Liliopsida
## spp Class
## 1 Sebastes borealis Actinopterygii# there are no plants or arachnid in areas of slow velocity##############################
# species with high GD and fast range shift
tmp <- mydatatogo %>%
filter(GD >0) %>%
filter(GD > quantile(GD, .9) & SHIFT_abs > quantile(SHIFT_abs, .9)) %>%
select(spp,Class,GD,SHIFT_abs) %>%
group_by(spp,Class) %>%
summarise(GD = round(mean(GD),5),
SHIFT_abs = round(mean(SHIFT_abs),3)) %>%
arrange(GD)
# these species
tmp## # A tibble: 45 × 4
## # Groups: spp [45]
## spp Class GD SHIFT_abs
## <fct> <fct> <dbl> <dbl>
## 1 Corvus corax Aves 0.0150 14.4
## 2 Chesias rufata Insecta 0.0152 9.62
## 3 Xysticus ulmi Arachnida 0.0152 8.61
## 4 Platycheirus angustatus Insecta 0.0156 11.2
## 5 Pachygnatha clercki Arachnida 0.0162 10.9
## 6 Vireo bellii Aves 0.0164 13.9
## 7 Maso sundevalli Arachnida 0.0164 9.38
## 8 Walckenaeria vigilax Arachnida 0.0170 10.4
## 9 Dysdera crocata Arachnida 0.0172 8.87
## 10 Salticus scenicus Arachnida 0.0173 10.8
## # ℹ 35 more rows# from these classes
table(tmp$Class) # mostly plants##
## Insecta Aves Magnoliopsida Polypodiopsida Arachnida
## 17 12 0 0 15
## Liliopsida Actinopterygii
## 0 1# do for each class
my_classes <- unique(mydatatogo$Class)
for(i in 1:length(my_classes)){
tmp <- mydatatogo %>%
filter(GD >0) %>%
filter(Class == my_classes[i]) %>%
filter(GD > quantile(GD, .9) & SHIFT_abs > quantile(SHIFT_abs, .9)) %>%
select(spp,Class) %>%
unique()
print(tmp)
}## spp Class
## 1 Andrena bicolor Insecta
## 3 Chesias rufata Insecta
## 4 Coenonympha tullia Insecta
## 5 Cordulia aenea Insecta
## 6 Eupeodes latifasciatus Insecta
## 7 Helophilus hybridus Insecta
## 8 Hydrobius fuscipes Insecta
## 9 Hypomecis punctinalis Insecta
## 10 Lasioglossum lativentre Insecta
## 11 Lasius niger Insecta
## 12 Myrmica scabrinodis Insecta
## 13 Nycterosea obstipata Insecta
## 14 Platycheirus angustatus Insecta
## 15 Trichoplusia ni Insecta
## spp Class
## 1 Ardea alba Aves
## 3 Athene cunicularia Aves
## 4 Garrulus glandarius Aves
## 5 Meleagris gallopavo Aves
## 10 Motacilla flava Aves
## 12 Numenius phaeopus Aves
## 14 Nycticorax nycticorax Aves
## 15 Phoenicurus phoenicurus Aves
## 16 Phylloscopus borealis Aves
## 17 Phylloscopus trochiloides Aves
## 18 Pinicola enucleator Aves
## spp Class
## 1 Rumex conglomeratus Magnoliopsida
## 2 Sagina procumbens Magnoliopsida
## 3 Vaccinium oxycoccos Magnoliopsida
## [1] spp Class
## <0 rows> (or 0-length row.names)
## spp Class
## 1 Haplodrassus signifer Arachnida
## 2 Larinioides cornutus Arachnida
## 3 Mitopus morio Arachnida
## 4 Steatoda bipunctata Arachnida
## [1] spp Class
## <0 rows> (or 0-length row.names)
## spp Class
## 1 Anarhichas denticulatus Actinopterygii
## 2 Mallotus villosus Actinopterygii
## 3 Zeus faber Actinopterygii# there are no plants or arachnid in areas of slow velocity##############################
# species with low genetic diversity and fast range shift
tmp <- mydatatogo %>%
filter(GD >0) %>%
filter(GD < quantile(GD, .1) & SHIFT_abs > quantile(SHIFT_abs, .9)) %>%
select(spp,Class,GD,SHIFT_abs) %>%
group_by(spp,Class) %>%
summarise(GD = round(mean(GD),5),
SHIFT_abs = round(mean(SHIFT_abs),3))
# these species
tmp## # A tibble: 35 × 4
## # Groups: spp [35]
## spp Class GD SHIFT_abs
## <fct> <fct> <dbl> <dbl>
## 1 Accipiter striatus Aves 0.00031 8.71
## 2 Acrocephalus schoenobaenus Aves 0.00087 12.6
## 3 Amaurornis phoenicurus Aves 0.00025 23.0
## 4 Andrena subopaca Insecta 0.00053 10.8
## 5 Araschnia levana Insecta 0.00012 20.5
## 6 Argynnis paphia Insecta 0.0005 39.2
## 7 Bombus hortorum Insecta 0.00049 13.4
## 8 Bombus terrestris Insecta 0.00085 13.6
## 9 Catocala fulminea Insecta 0.0008 12.0
## 10 Colias croceus Insecta 0.00016 14.8
## # ℹ 25 more rows# from these classes
table(tmp$Class) # mostly insects and some birds##
## Insecta Aves Magnoliopsida Polypodiopsida Arachnida
## 24 8 0 0 2
## Liliopsida Actinopterygii
## 0 1##############################
# species with low genetic diversity and fast range expansion at the leading edge
tmp <- mydatatogo %>%
filter(GD >0 & Param=="LE") %>%
filter(GD < quantile(GD, .1) & SHIFT_abs > quantile(SHIFT_abs, .9)) %>%
select(spp,Class,GD,SHIFT_abs) %>%
group_by(spp,Class) %>%
summarise(GD = round(mean(GD),5),
SHIFT_abs = round(mean(SHIFT_abs),3))
# these species
tmp## # A tibble: 21 × 4
## # Groups: spp [21]
## spp Class GD SHIFT_abs
## <fct> <fct> <dbl> <dbl>
## 1 Amaurornis phoenicurus Aves 0.00025 23.0
## 2 Andrena subopaca Insecta 0.00053 10.8
## 3 Araschnia levana Insecta 0.00012 32.5
## 4 Argynnis paphia Insecta 0.0005 39.2
## 5 Bombus hortorum Insecta 0.00049 13.4
## 6 Bombus terrestris Insecta 0.00085 16.2
## 7 Catocala fulminea Insecta 0.0008 12.0
## 8 Colias croceus Insecta 0.00016 14.8
## 9 Lithosia quadra Insecta 0.00075 15.1
## 10 Lomographa bimaculata Insecta 0.00076 13.1
## # ℹ 11 more rows# from these classes
table(tmp$Class) # mostly insects##
## Insecta Aves Magnoliopsida Polypodiopsida Arachnida
## 19 1 0 0 0
## Liliopsida Actinopterygii
## 0 16.7 Histograms
6.7.1 Per Class
p1 <- ggplot(mydatatogo, aes(SHIFT_abs))+
geom_density(aes(fill=Class),alpha=.5)+
scale_x_log10()+
scale_y_continuous(trans="log1p",labels = function(x)round(x,1))+
theme_classic() +
labs(y = "Density", x = bquote('Absolute velocity of range shift '(km.yr^1)))
p2 <- ggplot(mydatatogo, aes(GD))+
geom_density(aes(fill=Class),alpha=.5)+
scale_x_log10()+
scale_y_continuous(trans="log1p",labels = function(x)round(x,1))+
theme_classic() +
labs(y = "", x = "Genetic diversity")
p3 <- ggplot(mydatatogo, aes(vel_abs))+
# geom_histogram(binwidth = .5,aes(fill=Class),alpha=.5,color="black") +
geom_density(aes(fill=Class),alpha=.5,) +
scale_x_log10()+
scale_y_continuous(trans="log1p",labels = function(x)round(x,1))+
theme_classic() +
labs(y = "", x = bquote('Absolute velocity of isotherm shifts '(km.yr^1)))
legend <- cowplot::get_legend(
# create some space to the left of the legend
p1 + theme(legend.box.margin = margin(0, 0, 0, 12))
)## Warning in get_plot_component(plot, "guide-box"): Multiple components found;
## returning the first one. To return all, use `return_all = TRUE`.cowplot::plot_grid(plotlist = list(p1 + theme(legend.position = "none"),
p2 + theme(legend.position = "none"),
p3 + theme(legend.position = "none"),
legend),
rel_widths = c(3,3,3,2), align = "v", nrow=1)## Warning in scale_x_log10(): log-10 transformation introduced infinite values.## Warning: Removed 430 rows containing non-finite outside the scale range
## (`stat_density()`).## Warning: Graphs cannot be vertically aligned unless the axis parameter is set.
## Placing graphs unaligned.
6.7.2 Shift
ggplot(mydatatogo, aes(SHIFT_abs))+
geom_density(aes(fill = Param), color = NA, alpha = .3)+
theme_classic() +
labs(y = "Density", x = "Range shift")
6.7.3 GD
hist(mydatatogo$GD)
### GD vs Class
ggplot(mydatatogo, aes(GD))+
geom_density(aes(fill = Class), color = NA, alpha = .3)+
theme_classic() +
labs(y = "Density", x = "Nucleotide diversity (log)")
m <- aov(GD ~ Class, mydatatogo)
summary(m)## Df Sum Sq Mean Sq F value Pr(>F)
## Class 6 0.01983 0.003305 70.68 <2e-16 ***
## Residuals 4681 0.21891 0.000047
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tres <- TukeyHSD(m)
round(tres$Class,3)## diff lwr upr p adj
## Aves-Insecta 0.000 0.000 0.001 0.978
## Magnoliopsida-Insecta -0.004 -0.005 -0.003 0.000
## Polypodiopsida-Insecta -0.004 -0.009 0.000 0.038
## Arachnida-Insecta 0.004 0.002 0.005 0.000
## Liliopsida-Insecta -0.005 -0.007 -0.004 0.000
## Actinopterygii-Insecta -0.001 -0.003 0.000 0.058
## Magnoliopsida-Aves -0.005 -0.005 -0.004 0.000
## Polypodiopsida-Aves -0.005 -0.009 0.000 0.025
## Arachnida-Aves 0.004 0.002 0.005 0.000
## Liliopsida-Aves -0.006 -0.007 -0.004 0.000
## Actinopterygii-Aves -0.002 -0.003 0.000 0.016
## Polypodiopsida-Magnoliopsida 0.000 -0.004 0.004 1.000
## Arachnida-Magnoliopsida 0.008 0.007 0.010 0.000
## Liliopsida-Magnoliopsida -0.001 -0.003 0.001 0.573
## Actinopterygii-Magnoliopsida 0.003 0.001 0.005 0.000
## Arachnida-Polypodiopsida 0.008 0.004 0.013 0.000
## Liliopsida-Polypodiopsida -0.001 -0.006 0.003 0.992
## Actinopterygii-Polypodiopsida 0.003 -0.001 0.007 0.413
## Liliopsida-Arachnida -0.009 -0.011 -0.007 0.000
## Actinopterygii-Arachnida -0.005 -0.007 -0.003 0.000
## Actinopterygii-Liliopsida 0.004 0.002 0.006 0.000Genetic diversity differs among classes.
6.7.4 Velocity
hist(mydatatogo$vel_abs)
hist(mydatatogo$vel_abs_log)
We should use SHIFT_abs with a Gamma family or SHIFT_abs_log in a gaussian family.
6.7.5 Distance of GD data and range shift data
summary(mydatatogo$dist_cent/1000)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 55.14 1185.10 2162.06 3602.37 5633.14 19510.91hist(mydatatogo$dist_cent/1000,
main = "Distance from study area centroid",
xlab = "Distance in km")
summary(mydatatogo$dist_edge/1000)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 84.41 1243.71 2227.49 3648.08 5610.20 18961.05hist(mydatatogo$dist_edge/1000,
main = "Distance from the range edge of study area",
xlab = "Distance in km")
7 Scatter plots
7.1 SHIFT vs Velocity
ggplot(mydatatogo,
aes(x=vel_abs, y=SHIFT_abs, color=sqrt(GD)*100))+
scale_x_continuous(
trans="log",
labels=function(x) round(x,3),
breaks = exp(round(seq(min(mydatatogo$vel_abs_log),
max(mydatatogo$vel_abs_log),
length.out = 5),3)))+
scale_y_continuous(
trans="log",
labels=function(x) round(x,3),
breaks = exp(round(seq(min(mydatatogo$SHIFT_abs_log),
max(mydatatogo$SHIFT_abs_log),
length.out = 5),3)))+
scale_color_viridis_c(option = "D", direction = -1) +
geom_point(alpha = .3, size = 3)+
theme_classic()+
theme(legend.justification = 1,
legend.position="top",
aspect.ratio=1)+
labs(x = bquote('Absolute velocity of isotherm shifts '(km.yr^1)),
y = bquote('Absolute velocity of range shift '(km.yr^1)),
color = "Nucleotide diversity x100") +
geom_smooth(method = "lm", color = "red")
ggplot(mydatatogo,
aes(x=vel_abs, y=SHIFT_abs, color=sqrt(GD)*100))+
scale_x_continuous(
trans="log",
labels=function(x) round(x,3),
breaks = exp(round(seq(min(mydatatogo$vel_abs_log),
max(mydatatogo$vel_abs_log),
length.out = 5),3)))+
scale_y_continuous(
trans="log",
labels=function(x) round(x,3),
breaks = exp(round(seq(min(mydatatogo$SHIFT_abs_log),
max(mydatatogo$SHIFT_abs_log),
length.out = 5),3)))+
scale_color_viridis_c(option = "D", direction = -1) +
geom_point(alpha = .3, size = 3)+
theme_classic()+
theme(legend.justification = 1,
legend.position="top",
aspect.ratio=1)+
labs(x = bquote('Absolute velocity of isotherm shifts '(km.yr^1)),
y = bquote('Absolute velocity of range shift '(km.yr^1)),
color = "Nucleotide diversity x100") +
geom_smooth(method = "lm", color = "red") +
facet_wrap(.~Param)
ggplot(mydatatogo,
aes(x=vel_abs, y=SHIFT_abs, color=sqrt(GD)*100))+
scale_x_continuous(
trans="log",
labels=function(x) round(x,3),
breaks = exp(round(seq(min(mydatatogo$vel_abs_log),
max(mydatatogo$vel_abs_log),
length.out = 5),3)))+
scale_y_continuous(
trans="log",
labels=function(x) round(x,3),
breaks = exp(round(seq(min(mydatatogo$SHIFT_abs_log),
max(mydatatogo$SHIFT_abs_log),
length.out = 5),3)))+
scale_color_viridis_c(option = "D", direction = -1) +
geom_point(alpha = .3, size = 3)+
theme_classic()+
theme(legend.justification = 1,
legend.position="top",
aspect.ratio=1)+
labs(x = bquote('Absolute velocity of isotherm shifts '(km.yr^1)),
y = bquote('Absolute velocity of range shift '(km.yr^1)),
color = "Nucleotide diversity x100") +
geom_smooth(method = "lm", color = "red") +
facet_wrap(.~Class,scales = "free_x")
7.2 SHIFT vs GD
ggplot(mydatatogo,
aes(x=GD, y=SHIFT_abs,color=vel_abs))+
scale_y_continuous(
trans="log",
labels=function(x) round(x,3),
breaks = exp(round(seq(min(mydatatogo$SHIFT_abs_log),
max(mydatatogo$SHIFT_abs_log),
length.out = 6),3)))+
scale_color_viridis_c(option = "inferno",direction=-1)+
geom_point(alpha = .5, size = 3)+
theme_classic()+
theme(legend.justification = 1,
legend.position="top",
aspect.ratio=1)+
labs(x = "Nucleotide diversity",
y = bquote('Absolute velocity of range shift '(km.yr^1)),
color = bquote('Absolute velocity of isotherm shifts '(km.yr^1))) +
geom_smooth(method = "lm", color = "red")
ggplot(mydatatogo,
aes(x=GD, y=SHIFT_abs,color=vel_abs))+
scale_y_continuous(
trans="log",
labels=function(x) round(x,3),
breaks = exp(round(seq(min(mydatatogo$SHIFT_abs_log),
max(mydatatogo$SHIFT_abs_log),
length.out = 6),3)))+
scale_color_viridis_c(option = "inferno",direction=-1)+
geom_point(alpha = .5, size = 3)+
theme_classic()+
theme(legend.justification = 1,
legend.position="top",
aspect.ratio=1)+
labs(x = "Nucleotide diversity",
y = bquote('Absolute velocity of range shift '(km.yr^1)),
color = bquote('Absolute velocity of isotherm shifts '(km.yr^1))) +
geom_smooth(method = "lm", color = "red") +
facet_wrap(.~Param)
ggplot(mydatatogo,
aes(x=GD, y=SHIFT_abs, color=vel_abs))+
scale_y_continuous(
trans="log",
labels=function(x) round(x,3),
breaks = exp(round(seq(min(mydatatogo$SHIFT_abs_log),
max(mydatatogo$SHIFT_abs_log),
length.out = 6),3)))+
scale_color_viridis_c(option = "inferno",direction=-1)+
geom_point(alpha = .3, size = 3)+
theme_classic()+
theme(legend.justification = 1,
legend.position="top",
aspect.ratio=1)+
labs(x = "Nucleotide diversity",
y = bquote('Absolute velocity of range shift '(km.yr^1)),
color = bquote('Absolute velocity of isotherm shifts '(km.yr^1))) +
geom_smooth(method = "lm", color = "red") +
facet_wrap(Class~., scales = "free")
7.3 GD vs Velocity
ggplot(mydatatogo,
aes(x=vel_abs, y=GD))+
geom_point(alpha = .3, size = 3, shape=21, fill = "black")+
scale_x_continuous(
trans="log",
labels=function(x) round(x,3),
breaks = exp(round(seq(min(mydatatogo$vel_abs_log),
max(mydatatogo$vel_abs_log),
length.out = 5),3)))+
theme_classic()+
theme(aspect.ratio=1)+
labs(x = bquote('Absolute velocity of isotherm shifts '(km.yr^1)),
y = "Nucleotide diversity") +
geom_smooth(method = "lm", color = "red")
ggplot(mydatatogo,
aes(x=vel_abs, y=GD))+
geom_point(alpha = .3, size = 3, shape=21, fill = "black")+
scale_x_continuous(
trans="log",
labels=function(x) round(x,3),
breaks = exp(round(seq(min(mydatatogo$vel_abs_log),
max(mydatatogo$vel_abs_log),
length.out = 5),3)))+
theme_classic()+
theme(aspect.ratio=1)+
labs(x = bquote('Absolute velocity of isotherm shifts '(km.yr^1)),
y = "Nucleotide diversity") +
geom_smooth(method = "lm", color = "red") +
facet_wrap(.~Param)
7.4 GD vs Latitude
ggplot(mydatatogo,
aes(x=LatGen, y=GD))+
geom_point(alpha = .3, size = 3, shape=21, fill = "black")+
theme_classic()+
theme(aspect.ratio=1)+
labs(x = 'Absolute latitude from average GD centroid',
y = "Nucleotide diversity") +
geom_smooth(method = "lm", color = "red")
summary(lm(GD~LatGen,mydatatogo))##
## Call:
## lm(formula = GD ~ LatGen, data = mydatatogo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.006665 -0.004226 -0.002627 0.001114 0.047280
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.010e-03 5.102e-04 15.698 < 2e-16 ***
## LatGen -5.438e-05 1.064e-05 -5.112 3.32e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.007118 on 4686 degrees of freedom
## Multiple R-squared: 0.005545, Adjusted R-squared: 0.005333
## F-statistic: 26.13 on 1 and 4686 DF, p-value: 3.318e-07ggplot(mydatatogo,
aes(x=abs(Lat), y=GD))+
geom_point(alpha = .3, size = 3, shape=21, fill = "black")+
theme_classic()+
theme(aspect.ratio=1)+
labs(x = 'Absolute latitude from study area centroid',
y = "Nucleotide diversity") +
geom_smooth(method = "lm", color = "red")
summary(lm(GD~abs(Lat),mydatatogo))##
## Call:
## lm(formula = GD ~ abs(Lat), data = mydatatogo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.005495 -0.004249 -0.002711 0.001257 0.048839
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.267e-03 7.136e-04 7.381 1.85e-13 ***
## abs(Lat) 3.539e-06 1.320e-05 0.268 0.789
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.007138 on 4686 degrees of freedom
## Multiple R-squared: 1.534e-05, Adjusted R-squared: -0.0001981
## F-statistic: 0.07188 on 1 and 4686 DF, p-value: 0.78867.5 Velocity vs Latitude
ggplot(mydatatogo,
aes(x=abs(Lat), y=vel_abs))+
geom_point(alpha = .3, size = 3, shape=21, fill = "black")+
theme_classic()+
theme(aspect.ratio=1)+
labs(x = 'Absolute latitude from study area centroid',
y = bquote('Absolute velocity of isotherm shifts '(km.yr^1))) +
geom_smooth(method = "lm", color = "red")
summary(lm(vel_abs~abs(Lat),mydatatogo))##
## Call:
## lm(formula = vel_abs ~ abs(Lat), data = mydatatogo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.4155 -1.0096 -0.5917 0.8004 5.4244
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.394117 0.135111 2.917 0.00355 **
## abs(Lat) 0.017222 0.002499 6.890 6.31e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.351 on 4686 degrees of freedom
## Multiple R-squared: 0.01003, Adjusted R-squared: 0.009818
## F-statistic: 47.47 on 1 and 4686 DF, p-value: 6.31e-127.6 SHIFT vs TajimasD
ggplot(mydatatogo,
aes(x=TajimasD, y=SHIFT_abs,color=vel_abs))+
scale_color_viridis_c(option = "inferno",direction=-1)+
geom_point(alpha = .5, size = 3)+
scale_y_continuous(
trans="log",
labels=function(x) round(x,3),
breaks = exp(round(seq(min(mydatatogo$SHIFT_abs_log),
max(mydatatogo$SHIFT_abs_log),
length.out = 5),3)))+
theme_classic()+
theme(legend.justification = 1,
legend.position="top",
aspect.ratio=1)+
labs(x = "Tajima's D",
y = "Range shift velocity (abs)\nkm/year",
color = bquote('Absolute velocity of isotherm shifts '(km.yr^1))) +
geom_smooth(method = "lm", color = "red")
ggplot(mydatatogo,
aes(x=TajimasD, y=SHIFT_abs, color=vel_abs))+
scale_y_continuous(
trans="log",
labels=function(x) round(x,3),
breaks = exp(round(seq(min(mydatatogo$SHIFT_abs_log),
max(mydatatogo$SHIFT_abs_log),
length.out = 5),3)))+
scale_color_viridis_c(option = "inferno",direction=-1)+
geom_point(alpha = .5, size = 3)+
theme_classic()+
theme(legend.justification = 1,
legend.position="top",
aspect.ratio=1)+
labs(x = "Tajima's D",
y = "Range shift velocity (abs)\nkm/year",
color = bquote('Absolute velocity of isotherm shifts '(km.yr^1))) +
geom_smooth(method = "lm", color = "red") +
facet_wrap(.~Param)
ggplot(mydatatogo,
aes(x=TajimasD, y=SHIFT_abs, color=vel_abs))+
scale_y_continuous(
trans="log",
labels=function(x) round(x,3),
breaks = exp(round(seq(min(mydatatogo$SHIFT_abs_log),
max(mydatatogo$SHIFT_abs_log),
length.out = 5),3)))+
scale_color_viridis_c(option = "inferno",direction=-1)+
geom_point(alpha = .5, size = 3)+
theme_classic()+
theme(legend.justification = 1,
legend.position="top",
aspect.ratio=1)+
labs(x = "Tajima's D",
y = "Range shift velocity (abs)\nkm/year",
color = bquote('Absolute velocity of isotherm shifts '(km.yr^1))) +
geom_smooth(method = "lm", color = "red") +
facet_wrap(Class~., scales = "free")
7.7 Tajima’s D vs Velocity
ggplot(mydatatogo,
aes(x=vel_abs, y=TajimasD))+
scale_x_continuous(
trans="log",
labels=function(x) round(x,3),
breaks = exp(round(seq(min(mydatatogo$vel_abs_log),
max(mydatatogo$vel_abs_log),
length.out = 5),3)))+
geom_point(alpha = .3, size = 3, shape=21, fill = "black")+
theme_classic()+
theme(aspect.ratio=1)+
labs(x = bquote('Absolute velocity of isotherm shifts '(km.yr^1)),
y = "Tajima's D") +
geom_smooth(method = "lm", color = "red")
ggplot(mydatatogo,
aes(x=vel_abs, y=TajimasD))+
scale_x_continuous(
trans="log",
labels=function(x) round(x,3),
breaks = exp(round(seq(min(mydatatogo$vel_abs_log),
max(mydatatogo$vel_abs_log),
length.out = 5),3)))+
geom_point(alpha = .3, size = 3, shape=21, fill = "black")+
theme_classic()+
theme(aspect.ratio=1)+
labs(x = bquote('Absolute velocity of isotherm shifts '(km.yr^1)),
y = "Tajima's D") +
geom_smooth(method = "lm", color = "red") +
facet_wrap(.~Param)
7.8 N obs for modeling
7.9 N obs for modeling II
7.10 North America vs Europe
7.10.1 N data per Continent
## # A tibble: 5 × 3
## Continent StudyAreas RangeShifts
## <fct> <int> <int>
## 1 Europe 38 3494
## 2 NorthAmerica 29 1136
## 3 Australian 2 28
## 4 Asian 5 24
## 5 African 2 6Most data comes from Europe or North America. But also note Europe has about 3x more data than North America.
7.10.2 SHIFT vs Velocity
Most of the global pattern seems to be driven by Europe. North America shows not relationship, or a tendency of inverse relationship, between Range shift and Climate velocity (more lags in North America?).
Looking across range edges we can more clearly see that North America lags more than Europe, specially at the leading and trailing edges.
7.10.3 SHIFT vs GD
Also for GD, effects seem to be larger in Europe than North America.
Same across range edges…
7.10.4 Distribution of shift values
##
## Welch Two Sample t-test
##
## data: log(SHIFT_abs) by Continent
## t = -1.0793, df = 2080.7, p-value = 0.2806
## alternative hypothesis: true difference in means between group Europe and group NorthAmerica is not equal to 0
## 95 percent confidence interval:
## -0.16648677 0.04828797
## sample estimates:
## mean in group Europe mean in group NorthAmerica
## 0.2183777 0.27747718 Statistics
We explore how species range shifts (km/year) (response variable) varies in function of climate change velocity (km/year) and GD (Nucleotide diversity) using a general linear mixed-effect modelling (GLMM) framework. Our goal is to estimate whereas GD increases (or decreases) species range shift in response to climate change velocity. To account for methodological differences in range shift estimates, we included the methodological variables from BioShifts described above as fixed-effect terms in our models. As our dataset comprises species across multiple taxonomic classes, with different ecology and sensitivity to climate variability, observations are considered as non-independent from each other. Flight capacity would confer birds greater dispersal potential than plants, ectotherms and short lived like insects should be more sensitive to climate variability than endotherms like fishes. We control for taxonomy and ecological variability by including Class as a random-intercept in our model.
We model absolute values for range shifts and climate velocity as we expect GD to impact the magnitude of the coupling between range shifts and climate velocity, not the direction of this relationship. As we are modeling positive values, we fitted our GLMMs using a Gamma family with a log link. In addition to these models, we also fitted Gaussian models using log-transformed values of range shift and climate velocity. Gamma models performed better (in terms of R2 and AIC) than Gaussian models, therefore we discuss results from Gamma models only.
We test whereas the modulating effect of GD varies across the multiple range positions. Different levels of exposure are expected at different positions of a species range, thereby leading to influencing ecological dynamics. At the warm edge species leave closer to their upper thermal tolerance limit and thus should be more sensitive to warming. In contrast, a warming at the cold edge may represent colonization opportunity. Because GD is at the intraspecific scale, a single GD value can be matched to multiple estimates of range shift estimates across different positions of the range for the same species. To avoid false duplicate data, we run separate analyses to the different range positions.
In summary, three models are fitted:
allW: a model that considers the whole dataset (the position is accounting for in the model as a fixed and additive effect only due to the imbalance in data among positions) weighting for the imbalance in he number of obs among species (i.e., in the model every species will have the same weight whatever the number of obs)
CEW, LEW and TEW: three independent models that model centroid, leading and trailing edges range shifts, respectively. Like in the allW model, observations used by each model are weighted to control for imbalance in he number of obs among species (i.e., in the model every species will have the same weight whatever the number of obs)
Here is the formula for models W: \[ Range Shift = Climate Velocity * Genetic Diversity + RangePosition + Methods + (1|Class) \]
- allW2: differs from allW in two ways:
- the obs are weighted in order to control for the imbalance in both the number of obs among species (i.e., in the model every species will have the same weight whatever the number of obs) and the number of obs among the positions.
- as imbalance in data is controlled the position is accounted for an additive and interaction effect (with climate velocity and genetic diversity in the model).
Here is the formula for models W2: \[ Range Shift = Climate Velocity * Genetic Diversity * RangePosition + Methods + (1|Class) \]
We first test if the species range shift varies in function of both climate change velocity, genetic diversity and the interaction between these two variables, while account for methodology and phylogeny.
8.1 Gaussian model
gau_velxGD_edge1 <- as.formula(
"SHIFT_abs_log ~ scale(vel_abs_log) * scale(GD) + Param +
LogNtempUnits + LogExtent + ContinuousGrain + PrAb + Quality +
(1 | Class)")
gau_velxGD_edge2 <- as.formula(
"SHIFT_abs_log ~ scale(vel_abs_log) * scale(GD) + Param +
(1 | Class + NtempUnitsF + ExtentF + ContinuousGrain + PrAb + Quality)")# Model with zero GD + fixed-effect methods
gau1_allW <- glmmTMB(gau_velxGD_edge1,
weights = weight_obs_spp,
data = mydatatogo)
summary(gau1_allW)## Family: gaussian ( identity )
## Formula:
## SHIFT_abs_log ~ scale(vel_abs_log) * scale(GD) + Param + LogNtempUnits +
## LogExtent + ContinuousGrain + PrAb + Quality + (1 | Class)
## Data: mydatatogo
## Weights: weight_obs_spp
##
## AIC BIC logLik deviance df.resid
## 31.3 121.7 -1.7 3.3 4674
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Class (Intercept) 5.654e-10 2.378e-05
## Residual 1.647e+00 1.283e+00
## Number of obs: 4688, groups: Class, 7
##
## Dispersion estimate for gaussian family (sigma^2): 1.65
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.63155 16.48208 -0.281 0.779
## scale(vel_abs_log) 0.26049 1.61659 0.161 0.872
## scale(GD) 0.09111 1.23484 0.074 0.941
## ParamLE 1.32396 3.97813 0.333 0.739
## ParamTE 0.38772 7.34213 0.053 0.958
## LogNtempUnits -0.08744 2.23778 -0.039 0.969
## LogExtent 0.55938 2.68848 0.208 0.835
## ContinuousGrain 0.14666 1.84439 0.080 0.937
## PrAbabundance 0.80590 5.11069 0.158 0.875
## QualityRAW 0.23119 3.29523 0.070 0.944
## QualityRESURVEYED -0.52969 4.98264 -0.106 0.915
## scale(vel_abs_log):scale(GD) -0.06014 1.40771 -0.043 0.966round(MuMIn::r.squaredGLMM(gau1_allW),2)## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.## Warning in r.squaredGLMM.glmmTMB(gau1_allW): the effects of zero-inflation and
## dispersion model are ignored## R2m R2c
## [1,] 0.35 0.35AIC(gau1_allW)## [1] 31.33657performance::check_collinearity(gau1_allW)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## scale(vel_abs_log) 1.34 [1.30, 1.39] 1.16 0.75
## scale(GD) 1.04 [1.02, 1.09] 1.02 0.96
## Param 2.61 [2.49, 2.74] 1.62 0.38
## LogNtempUnits 1.79 [1.72, 1.86] 1.34 0.56
## LogExtent 1.75 [1.68, 1.83] 1.32 0.57
## ContinuousGrain 1.35 [1.31, 1.40] 1.16 0.74
## PrAb 2.25 [2.15, 2.35] 1.50 0.45
## Quality 1.85 [1.78, 1.93] 1.36 0.54
## scale(vel_abs_log):scale(GD) 1.07 [1.05, 1.11] 1.04 0.93
## Tolerance 95% CI
## [0.72, 0.77]
## [0.92, 0.98]
## [0.37, 0.40]
## [0.54, 0.58]
## [0.55, 0.59]
## [0.71, 0.77]
## [0.43, 0.47]
## [0.52, 0.56]
## [0.90, 0.95]# Model without zero GD + fixed-effect methods
gau2_allW <- glmmTMB(gau_velxGD_edge1,
weights = weight_obs_spp,
data = mydatatogo %>%
filter_all(all_vars(!is.infinite(.))))
summary(gau2_allW)## Family: gaussian ( identity )
## Formula:
## SHIFT_abs_log ~ scale(vel_abs_log) * scale(GD) + Param + LogNtempUnits +
## LogExtent + ContinuousGrain + PrAb + Quality + (1 | Class)
## Data: mydatatogo %>% filter_all(all_vars(!is.infinite(.)))
## Weights: weight_obs_spp
##
## AIC BIC logLik deviance df.resid
## 31.3 121.7 -1.7 3.3 4674
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Class (Intercept) 5.654e-10 2.378e-05
## Residual 1.647e+00 1.283e+00
## Number of obs: 4688, groups: Class, 7
##
## Dispersion estimate for gaussian family (sigma^2): 1.65
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.63155 16.48208 -0.281 0.779
## scale(vel_abs_log) 0.26049 1.61659 0.161 0.872
## scale(GD) 0.09111 1.23484 0.074 0.941
## ParamLE 1.32396 3.97813 0.333 0.739
## ParamTE 0.38772 7.34213 0.053 0.958
## LogNtempUnits -0.08744 2.23778 -0.039 0.969
## LogExtent 0.55938 2.68848 0.208 0.835
## ContinuousGrain 0.14666 1.84439 0.080 0.937
## PrAbabundance 0.80590 5.11069 0.158 0.875
## QualityRAW 0.23119 3.29523 0.070 0.944
## QualityRESURVEYED -0.52969 4.98264 -0.106 0.915
## scale(vel_abs_log):scale(GD) -0.06014 1.40771 -0.043 0.966round(MuMIn::r.squaredGLMM(gau2_allW),2)## Warning in r.squaredGLMM.glmmTMB(gau2_allW): the effects of zero-inflation and
## dispersion model are ignored## R2m R2c
## [1,] 0.35 0.35AIC(gau2_allW)## [1] 31.33657performance::check_collinearity(gau2_allW)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## scale(vel_abs_log) 1.34 [1.30, 1.39] 1.16 0.75
## scale(GD) 1.04 [1.02, 1.09] 1.02 0.96
## Param 2.61 [2.49, 2.74] 1.62 0.38
## LogNtempUnits 1.79 [1.72, 1.86] 1.34 0.56
## LogExtent 1.75 [1.68, 1.83] 1.32 0.57
## ContinuousGrain 1.35 [1.31, 1.40] 1.16 0.74
## PrAb 2.25 [2.15, 2.35] 1.50 0.45
## Quality 1.85 [1.78, 1.93] 1.36 0.54
## scale(vel_abs_log):scale(GD) 1.07 [1.05, 1.11] 1.04 0.93
## Tolerance 95% CI
## [0.72, 0.77]
## [0.92, 0.98]
## [0.37, 0.40]
## [0.54, 0.58]
## [0.55, 0.59]
## [0.71, 0.77]
## [0.43, 0.47]
## [0.52, 0.56]
## [0.90, 0.95]# Model with zero GD + random-effect methods
gau3_allW <- glmmTMB(gau_velxGD_edge2,
weights = weight_obs_spp,
data = mydatatogo)
summary(gau3_allW)## Family: gaussian ( identity )
## Formula:
## SHIFT_abs_log ~ scale(vel_abs_log) * scale(GD) + Param + (1 |
## Class + NtempUnitsF + ExtentF + ContinuousGrain + PrAb + Quality)
## Data: mydatatogo
## Weights: weight_obs_spp
##
## AIC BIC logLik deviance df.resid
## 29.5 113.4 -1.7 3.5 4675
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Class (Intercept) 2.959e-11 5.439e-06
## NtempUnitsF (Intercept) 1.179e-11 3.433e-06
## ExtentF (Intercept) 1.549e-11 3.936e-06
## ContinuousGrain (Intercept) 2.826e-11 5.316e-06
## PrAb (Intercept) 1.433e-11 3.785e-06
## Quality (Intercept) 2.799e-11 5.291e-06
## Residual 1.929e+00 1.389e+00
## Number of obs: 4688, groups:
## Class, 7; NtempUnitsF, 9; ExtentF, 7; ContinuousGrain, 4; PrAb, 2; Quality, 3
##
## Dispersion estimate for gaussian family (sigma^2): 1.93
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.45721 2.41040 -0.190 0.850
## scale(vel_abs_log) 0.36842 1.70571 0.216 0.829
## scale(GD) 0.09615 1.33201 0.072 0.942
## ParamLE 1.52535 3.33898 0.457 0.648
## ParamTE 0.67865 7.24845 0.094 0.925
## scale(vel_abs_log):scale(GD) -0.11118 1.51143 -0.074 0.941round(MuMIn::r.squaredGLMM(gau3_allW),2)## Warning in r.squaredGLMM.glmmTMB(gau3_allW): the effects of zero-inflation and
## dispersion model are ignored## R2m R2c
## [1,] 0.22 0.22AIC(gau3_allW)## [1] 29.4947performance::check_collinearity(gau3_allW)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## scale(vel_abs_log) 1.27 [1.23, 1.32] 1.13 0.78
## scale(GD) 1.03 [1.01, 1.08] 1.02 0.97
## Param 1.33 [1.29, 1.38] 1.15 0.75
## scale(vel_abs_log):scale(GD) 1.06 [1.03, 1.10] 1.03 0.95
## Tolerance 95% CI
## [0.76, 0.81]
## [0.92, 0.99]
## [0.72, 0.78]
## [0.91, 0.97]# Model without zero GD + random-effect methods
gau4_allW <- glmmTMB(gau_velxGD_edge2,
weights = weight_obs_spp,
data = mydatatogo %>%
filter_all(all_vars(!is.infinite(.))))
summary(gau4_allW)## Family: gaussian ( identity )
## Formula:
## SHIFT_abs_log ~ scale(vel_abs_log) * scale(GD) + Param + (1 |
## Class + NtempUnitsF + ExtentF + ContinuousGrain + PrAb + Quality)
## Data: mydatatogo %>% filter_all(all_vars(!is.infinite(.)))
## Weights: weight_obs_spp
##
## AIC BIC logLik deviance df.resid
## 29.5 113.4 -1.7 3.5 4675
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Class (Intercept) 2.959e-11 5.439e-06
## NtempUnitsF (Intercept) 1.179e-11 3.433e-06
## ExtentF (Intercept) 1.549e-11 3.936e-06
## ContinuousGrain (Intercept) 2.826e-11 5.316e-06
## PrAb (Intercept) 1.433e-11 3.785e-06
## Quality (Intercept) 2.799e-11 5.291e-06
## Residual 1.929e+00 1.389e+00
## Number of obs: 4688, groups:
## Class, 7; NtempUnitsF, 9; ExtentF, 7; ContinuousGrain, 4; PrAb, 2; Quality, 3
##
## Dispersion estimate for gaussian family (sigma^2): 1.93
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.45721 2.41040 -0.190 0.850
## scale(vel_abs_log) 0.36842 1.70571 0.216 0.829
## scale(GD) 0.09615 1.33201 0.072 0.942
## ParamLE 1.52535 3.33898 0.457 0.648
## ParamTE 0.67865 7.24845 0.094 0.925
## scale(vel_abs_log):scale(GD) -0.11118 1.51143 -0.074 0.941round(MuMIn::r.squaredGLMM(gau4_allW),2)## Warning in r.squaredGLMM.glmmTMB(gau4_allW): the effects of zero-inflation and
## dispersion model are ignored## R2m R2c
## [1,] 0.22 0.22AIC(gau4_allW)## [1] 29.4947performance::check_collinearity(gau4_allW)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## scale(vel_abs_log) 1.27 [1.23, 1.32] 1.13 0.78
## scale(GD) 1.03 [1.01, 1.08] 1.02 0.97
## Param 1.33 [1.29, 1.38] 1.15 0.75
## scale(vel_abs_log):scale(GD) 1.06 [1.03, 1.10] 1.03 0.95
## Tolerance 95% CI
## [0.76, 0.81]
## [0.92, 0.99]
## [0.72, 0.78]
## [0.91, 0.97]8.2 Gamma model
gam_velxGD_edge1 <- as.formula(
"SHIFT_abs ~ scale(vel_abs) * scale(GD) + Param +
LogNtempUnits + LogExtent + ContinuousGrain + PrAb + Quality +
(1 | Class)")
gam_velxGD_edge2 <- as.formula(
"SHIFT_abs ~ scale(vel_abs) * scale(GD) + Param +
(1 | Class + NtempUnitsF + ExtentF + ContinuousGrain + PrAb + Quality)")
gam_velxGDxedge1 <- as.formula(
"SHIFT_abs ~ scale(vel_abs) * scale(GD) * Param +
LogNtempUnits + LogExtent + ContinuousGrain + PrAb + Quality +
(1 | Class)")
gam_velxGDxedge2 <- as.formula(
"SHIFT_abs ~ scale(vel_abs) * scale(GD) * Param +
(1 | Class + NtempUnitsF + ExtentF + ContinuousGrain + PrAb + Quality)")
gam_velxGD1 <- as.formula(
"SHIFT_abs ~ scale(vel_abs) * scale(GD) +
LogNtempUnits + LogExtent + ContinuousGrain + PrAb + Quality +
(1 | Class)")
gam_velxGD2 <- as.formula(
"SHIFT_abs ~ scale(vel_abs) * scale(GD) +
(1 | Class + NtempUnitsF + ExtentF + ContinuousGrain + PrAb + Quality)")# Model with zero GD + fixed-effect methods
gam1_allW <- glmmTMB(gam_velxGD_edge1,
weights = weight_obs_spp,
family = Gamma(link = "log"),
data = mydatatogo)
summary(gam1_allW)## Family: Gamma ( log )
## Formula:
## SHIFT_abs ~ scale(vel_abs) * scale(GD) + Param + LogNtempUnits +
## LogExtent + ContinuousGrain + PrAb + Quality + (1 | Class)
## Data: mydatatogo
## Weights: weight_obs_spp
##
## AIC BIC logLik deviance df.resid
## 32.0 122.4 -2.0 4.0 4674
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Class (Intercept) 3.997e-10 1.999e-05
## Number of obs: 4688, groups: Class, 7
##
## Dispersion estimate for Gamma family (sigma^2): 1.01
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.41138 13.23528 -0.333 0.739
## scale(vel_abs) 0.25962 1.25963 0.206 0.837
## scale(GD) 0.05690 0.97526 0.058 0.953
## ParamLE 1.37732 3.25643 0.423 0.672
## ParamTE 0.48468 5.58047 0.087 0.931
## LogNtempUnits -0.06834 1.65348 -0.041 0.967
## LogExtent 0.63347 2.13100 0.297 0.766
## ContinuousGrain 0.07478 1.42256 0.053 0.958
## PrAbabundance 0.74783 4.05792 0.184 0.854
## QualityRAW 0.21596 2.50755 0.086 0.931
## QualityRESURVEYED -0.34747 4.07881 -0.085 0.932
## scale(vel_abs):scale(GD) -0.03098 0.95395 -0.032 0.974round(MuMIn::r.squaredGLMM(gam1_allW),2)## Warning in r.squaredGLMM.glmmTMB(gam1_allW): the effects of zero-inflation and
## dispersion model are ignored## R2m R2c
## delta 0.44 0.44
## lognormal 0.53 0.53
## trigamma 0.32 0.32AIC(gam1_allW)## [1] 32.01896performance::check_collinearity(gam1_allW)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## scale(vel_abs) 1.51 [1.45, 1.57] 1.23 0.66
## scale(GD) 1.03 [1.01, 1.08] 1.01 0.97
## Param 2.72 [2.60, 2.85] 1.65 0.37
## LogNtempUnits 1.63 [1.57, 1.70] 1.28 0.61
## LogExtent 1.57 [1.51, 1.63] 1.25 0.64
## ContinuousGrain 1.33 [1.28, 1.38] 1.15 0.75
## PrAb 2.31 [2.21, 2.42] 1.52 0.43
## Quality 1.94 [1.86, 2.02] 1.39 0.52
## scale(vel_abs):scale(GD) 1.04 [1.02, 1.09] 1.02 0.96
## Tolerance 95% CI
## [0.64, 0.69]
## [0.92, 0.99]
## [0.35, 0.38]
## [0.59, 0.64]
## [0.61, 0.66]
## [0.73, 0.78]
## [0.41, 0.45]
## [0.49, 0.54]
## [0.92, 0.98]# Model without zero GD + fixed-effect methods
gam2_allW <- glmmTMB(gam_velxGD_edge1,
weights = weight_obs_spp,
family = Gamma(link = "log"),
data = mydatatogo %>%
filter_all(all_vars(!is.infinite(.))))
summary(gam2_allW)## Family: Gamma ( log )
## Formula:
## SHIFT_abs ~ scale(vel_abs) * scale(GD) + Param + LogNtempUnits +
## LogExtent + ContinuousGrain + PrAb + Quality + (1 | Class)
## Data: mydatatogo %>% filter_all(all_vars(!is.infinite(.)))
## Weights: weight_obs_spp
##
## AIC BIC logLik deviance df.resid
## 32.0 122.4 -2.0 4.0 4674
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Class (Intercept) 3.997e-10 1.999e-05
## Number of obs: 4688, groups: Class, 7
##
## Dispersion estimate for Gamma family (sigma^2): 1.01
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.41138 13.23528 -0.333 0.739
## scale(vel_abs) 0.25962 1.25963 0.206 0.837
## scale(GD) 0.05690 0.97526 0.058 0.953
## ParamLE 1.37732 3.25643 0.423 0.672
## ParamTE 0.48468 5.58047 0.087 0.931
## LogNtempUnits -0.06834 1.65348 -0.041 0.967
## LogExtent 0.63347 2.13100 0.297 0.766
## ContinuousGrain 0.07478 1.42256 0.053 0.958
## PrAbabundance 0.74783 4.05792 0.184 0.854
## QualityRAW 0.21596 2.50755 0.086 0.931
## QualityRESURVEYED -0.34747 4.07881 -0.085 0.932
## scale(vel_abs):scale(GD) -0.03098 0.95395 -0.032 0.974round(MuMIn::r.squaredGLMM(gam2_allW),2)## Warning in r.squaredGLMM.glmmTMB(gam2_allW): the effects of zero-inflation and
## dispersion model are ignored## R2m R2c
## delta 0.44 0.44
## lognormal 0.53 0.53
## trigamma 0.32 0.32AIC(gam2_allW)## [1] 32.01896performance::check_collinearity(gam2_allW)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## scale(vel_abs) 1.51 [1.45, 1.57] 1.23 0.66
## scale(GD) 1.03 [1.01, 1.08] 1.01 0.97
## Param 2.72 [2.60, 2.85] 1.65 0.37
## LogNtempUnits 1.63 [1.57, 1.70] 1.28 0.61
## LogExtent 1.57 [1.51, 1.63] 1.25 0.64
## ContinuousGrain 1.33 [1.28, 1.38] 1.15 0.75
## PrAb 2.31 [2.21, 2.42] 1.52 0.43
## Quality 1.94 [1.86, 2.02] 1.39 0.52
## scale(vel_abs):scale(GD) 1.04 [1.02, 1.09] 1.02 0.96
## Tolerance 95% CI
## [0.64, 0.69]
## [0.92, 0.99]
## [0.35, 0.38]
## [0.59, 0.64]
## [0.61, 0.66]
## [0.73, 0.78]
## [0.41, 0.45]
## [0.49, 0.54]
## [0.92, 0.98]# Model with zero GD + random-effect methods
gam3_allW <- glmmTMB(gam_velxGD_edge2,
weights = weight_obs_spp,
family = Gamma(link = "log"),
data = mydatatogo)
summary(gam3_allW)## Family: Gamma ( log )
## Formula:
## SHIFT_abs ~ scale(vel_abs) * scale(GD) + Param + (1 | Class +
## NtempUnitsF + ExtentF + ContinuousGrain + PrAb + Quality)
## Data: mydatatogo
## Weights: weight_obs_spp
##
## AIC BIC logLik deviance df.resid
## 30.3 114.1 -2.1 4.3 4675
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Class (Intercept) 2.435e-10 1.560e-05
## NtempUnitsF (Intercept) 2.983e-11 5.462e-06
## ExtentF (Intercept) 1.745e-10 1.321e-05
## ContinuousGrain (Intercept) 1.397e-10 1.182e-05
## PrAb (Intercept) 2.158e-11 4.645e-06
## Quality (Intercept) 9.462e-11 9.727e-06
## Number of obs: 4688, groups:
## Class, 7; NtempUnitsF, 9; ExtentF, 7; ContinuousGrain, 4; PrAb, 2; Quality, 3
##
## Dispersion estimate for Gamma family (sigma^2): 1.21
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.44945 1.93175 0.233 0.816
## scale(vel_abs) 0.27806 1.34445 0.207 0.836
## scale(GD) 0.05929 1.08378 0.055 0.956
## ParamLE 1.18622 2.71805 0.436 0.663
## ParamTE 0.70271 5.73465 0.122 0.902
## scale(vel_abs):scale(GD) -0.06090 1.02973 -0.059 0.953round(MuMIn::r.squaredGLMM(gam3_allW),2)## Warning in r.squaredGLMM.glmmTMB(gam3_allW): the effects of zero-inflation and
## dispersion model are ignored## R2m R2c
## delta 0.20 0.20
## lognormal 0.28 0.28
## trigamma 0.12 0.12AIC(gam3_allW)## [1] 30.25481performance::check_collinearity(gam3_allW)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## scale(vel_abs) 1.41 [1.36, 1.46] 1.19 0.71
## scale(GD) 1.01 [1.00, 1.14] 1.01 0.99
## Param 1.41 [1.36, 1.46] 1.19 0.71
## scale(vel_abs):scale(GD) 1.03 [1.01, 1.08] 1.01 0.97
## Tolerance 95% CI
## [0.68, 0.74]
## [0.87, 1.00]
## [0.68, 0.74]
## [0.92, 0.99]# Model without zero GD + random-effect methods
gam4_allW <- glmmTMB(gam_velxGD_edge2,
weights = weight_obs_spp,
family = Gamma(link = "log"),
data = mydatatogo %>%
filter_all(all_vars(!is.infinite(.))))
summary(gam4_allW)## Family: Gamma ( log )
## Formula:
## SHIFT_abs ~ scale(vel_abs) * scale(GD) + Param + (1 | Class +
## NtempUnitsF + ExtentF + ContinuousGrain + PrAb + Quality)
## Data: mydatatogo %>% filter_all(all_vars(!is.infinite(.)))
## Weights: weight_obs_spp
##
## AIC BIC logLik deviance df.resid
## 30.3 114.1 -2.1 4.3 4675
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Class (Intercept) 2.435e-10 1.560e-05
## NtempUnitsF (Intercept) 2.983e-11 5.462e-06
## ExtentF (Intercept) 1.745e-10 1.321e-05
## ContinuousGrain (Intercept) 1.397e-10 1.182e-05
## PrAb (Intercept) 2.158e-11 4.645e-06
## Quality (Intercept) 9.462e-11 9.727e-06
## Number of obs: 4688, groups:
## Class, 7; NtempUnitsF, 9; ExtentF, 7; ContinuousGrain, 4; PrAb, 2; Quality, 3
##
## Dispersion estimate for Gamma family (sigma^2): 1.21
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.44945 1.93175 0.233 0.816
## scale(vel_abs) 0.27806 1.34445 0.207 0.836
## scale(GD) 0.05929 1.08378 0.055 0.956
## ParamLE 1.18622 2.71805 0.436 0.663
## ParamTE 0.70271 5.73465 0.122 0.902
## scale(vel_abs):scale(GD) -0.06090 1.02973 -0.059 0.953round(MuMIn::r.squaredGLMM(gam4_allW),2)## Warning in r.squaredGLMM.glmmTMB(gam4_allW): the effects of zero-inflation and
## dispersion model are ignored## R2m R2c
## delta 0.20 0.20
## lognormal 0.28 0.28
## trigamma 0.12 0.12AIC(gam4_allW)## [1] 30.25481performance::check_collinearity(gam4_allW)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI Increased SE Tolerance
## scale(vel_abs) 1.41 [1.36, 1.46] 1.19 0.71
## scale(GD) 1.01 [1.00, 1.14] 1.01 0.99
## Param 1.41 [1.36, 1.46] 1.19 0.71
## scale(vel_abs):scale(GD) 1.03 [1.01, 1.08] 1.01 0.97
## Tolerance 95% CI
## [0.68, 0.74]
## [0.87, 1.00]
## [0.68, 0.74]
## [0.92, 0.99]8.3 Summary
gau_models <- list(gau1_allW,gau2_allW,gau3_allW,gau4_allW)
gam_models_R2 <- sapply(gau_models, function(x) round(MuMIn::r.squaredGLMM(x),2))## Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
## dispersion model are ignored
## Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
## dispersion model are ignored
## Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
## dispersion model are ignored
## Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
## dispersion model are ignoredresults_gau <- data.frame(model = c("gau1_allW","gau2_allW","gau3_allW","gau4_allW"),
formula = c("gau_velxGD_edge1", "gau_velxGD_edge1", "gau_velxGD_edge2", "gau_velxGD_edge2"),
GD = rep(c("with zeros","without zeros")),
R2m = gam_models_R2[1,],
R2c = gam_models_R2[2,],
AIC = round(sapply(gau_models,AIC),1))
DT::datatable(results_gau,
rownames = FALSE,
extensions = 'Buttons',
options = list(pageLength = 100,
dom = 'B',
buttons = c('excel'))) %>%
DT::formatStyle(columns = colnames(results_gau), `fontSize` = '80%')gam_models <- list(gam1_allW=gam1_allW,gam2_allW=gam2_allW,gam3_allW=gam3_allW,gam4_allW=gam4_allW)
gam_models_R2 <- sapply(gam_models, function(x) round(MuMIn::r.squaredGLMM(x),2))## Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
## dispersion model are ignored
## Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
## dispersion model are ignored
## Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
## dispersion model are ignored
## Warning in r.squaredGLMM.glmmTMB(x): the effects of zero-inflation and
## dispersion model are ignoredresults_gam <- data.frame(model = c("gam1_allW","gam2_allW","gam3_allW","gam4_allW"),
formula = c("gam_velxGD_edge1", "gam_velxGD_edge1", "gam_velxGD_edge2", "gam_velxGD_edge2"),
GD = rep(c("with zeros","without zeros")),
R2m_delta = gam_models_R2[1,],
R2m_lognormal = gam_models_R2[2,],
R2m_trigamma = gam_models_R2[3,],
R2c_delta = gam_models_R2[4,],
R2c_lognormal = gam_models_R2[5,],
R2c_trigamma = gam_models_R2[6,],
AIC = round(sapply(gam_models,AIC),1))
DT::datatable(results_gam,
rownames = FALSE,
extensions = 'Buttons',
options = list(pageLength = 100,
dom = 'B',
buttons = c('excel'))) %>%
DT::formatStyle(columns = colnames(results_gam), `fontSize` = '80%')8.4 Bootstrap models
# Global parameters
no=2 # CPUs saved -- not used for parallel computing
nCPU=parallel::detectCores()-no # CPUs for parallel computing
s1=10 # the seed to make reproducible random stuff
nB=5000 # the number of bootstrap
#bootstrap test: H0= coefficient equal to 0; H1= coefficient higher than 0
test1=function(x,mu=0){
return(1-(length(x[x>mu])/length(x)))
}
#bootstrap test: H0= coefficient equal to 0; H1= coefficient lower than 0
test2=function(x,mu=0){
return(1-(length(x[x<mu])/length(x)))
}8.4.1 allW
my.cluster2 <- parallel::makeCluster(
nCPU,
type = "PSOCK"
)
clusterExport(cl = my.cluster2, varlist = c("s1","nB","glmmTMB","mydatatogo","gam_velxGD_edge1","ranef"))
#register cluster
doParallel::registerDoParallel(cl = my.cluster2)
ex=pblapply(1:nB, function(i) {
gc(reset=T)
print(i)
set.seed(s1+i)
n1=sample(1:nrow(mydatatogo),rep=T,size=nrow(mydatatogo))
mydatatogo2=mydatatogo[n1,]
x1=data.frame(table(mydatatogo2$spp))
x1=subset(x1,Freq>0)
x1$weight_obs_spp=(1/x1$Freq)*(1/nrow(x1))
mydatatogo2=merge(mydatatogo2[,1:48],x1[,c(1,3)],by.x="spp",by.y="Var1")
gamX1 <- glmmTMB(gam_velxGD_edge1,
family = Gamma(link = "log"),
weights = weight_obs_spp,
data = mydatatogo2)
r2=round(MuMIn::r.squaredGLMM(gamX1),2)
r2=data.frame(r2m=r2[1,1],
r2c=r2[1,2],
moy_GD=mean(mydatatogo2$GD),
sd_GD=sd(mydatatogo2$GD),
moy_VC=mean(mydatatogo2$vel_abs),
sd_VC=sd(mydatatogo2$vel_abs),
nB=i)
coeff=data.frame(summary(gamX1)$coeff$cond,
var=row.names(summary(gamX1)$coeff$cond),
nB=i)
randEff=as.data.frame(ranef(gamX1,condVar=T))
randEff$nB=i
return(list(r2,coeff,randEff))
}, cl = my.cluster2)
parallel::stopCluster(cl = my.cluster2)
nom="allW"
randEff=rlist::list.rbind(lapply(ex,"[[",3))
coeff=rlist::list.rbind(lapply(ex,"[[",2))
r2=rlist::list.rbind(lapply(ex,"[[",1))
write.table(r2,
here::here("Output",paste0("R2_",nom,".csv")),
sep=";",dec=".",row=F)
write.table(coeff,
here::here("Output",paste0("coeff_",nom,".csv")),
sep=";",dec=".",row=F)
write.table(randEff,here::here("Output",paste0("randEff_",nom,".csv")),
sep=";",dec=".",row=F)8.4.2 allW2
my.cluster2 <- parallel::makeCluster(
nCPU,
type = "PSOCK"
)
clusterExport(cl = my.cluster2, varlist = c("s1","nB","glmmTMB","mydatatogo","gam_velxGDxedge1","ranef"))
#register cluster
doParallel::registerDoParallel(cl = my.cluster2)
ex=pblapply(1:nB, function(i) {
gc(reset=T)
print(i)
set.seed(s1+i)
n1=sample(1:nrow(mydatatogo),rep=T,size=nrow(mydatatogo))
mydatatogo2=mydatatogo[n1,]
mydatatogoCE=subset(mydatatogo2,Param=="O")
x1=data.frame(table(mydatatogoCE$spp))
x1=subset(x1,Freq>0)
x1$weight_obs_spp=(1/x1$Freq)*(1/nrow(x1))
mydatatogoCE=merge(mydatatogoCE[,1:48],x1[,c(1,3)],by.x="spp",by.y="Var1")
mydatatogoLE=subset(mydatatogo2,Param=="LE")
x1=data.frame(table(mydatatogoLE$spp))
x1=subset(x1,Freq>0)
x1$weight_obs_spp=(1/x1$Freq)*(1/nrow(x1))
mydatatogoLE=merge(mydatatogoLE[,1:48],x1[,c(1,3)],by.x="spp",by.y="Var1")
mydatatogoTE=subset(mydatatogo2,Param=="TE")
x1=data.frame(table(mydatatogoTE$spp))
x1=subset(x1,Freq>0)
x1$weight_obs_spp=(1/x1$Freq)*(1/nrow(x1))
mydatatogoTE=merge(mydatatogoTE[,1:48],x1[,c(1,3)],by.x="spp",by.y="Var1")
mydatatogoTE$weight_obs_spp=mydatatogoTE$weight_obs_spp*(1/3)
mydatatogoCE$weight_obs_spp=mydatatogoCE$weight_obs_spp*(1/3)
mydatatogoLE$weight_obs_spp=mydatatogoLE$weight_obs_spp*(1/3)
mydatatogo2=rbind(mydatatogoTE,mydatatogoCE,mydatatogoLE)
gamX1 <- glmmTMB(gam_velxGDxedge1,
family = Gamma(link = "log"),
weights = weight_obs_spp,
data = mydatatogo2)
r2=round(MuMIn::r.squaredGLMM(gamX1),2)
r2=data.frame(r2m=r2[1,1],
r2c=r2[1,2],
moy_GD=mean(mydatatogo2$GD),
sd_GD=sd(mydatatogo2$GD),
moy_VC=mean(mydatatogo2$vel_abs),
sd_VC=sd(mydatatogo2$vel_abs),
nB=i)
coeff=data.frame(summary(gamX1)$coeff$cond,
var=row.names(summary(gamX1)$coeff$cond),
nB=i)
randEff=as.data.frame(ranef(gamX1,condVar=T))
randEff$nB=i
return(list(r2,coeff,randEff))
}, cl = my.cluster2)
parallel::stopCluster(cl = my.cluster2)
nom="allW2"
randEff=rlist::list.rbind(lapply(ex,"[[",3))
coeff=rlist::list.rbind(lapply(ex,"[[",2))
r2=rlist::list.rbind(lapply(ex,"[[",1))
write.table(r2,
here::here("Output",paste0("R2_",nom,".csv")),
sep=";",dec=".",row=F)
write.table(coeff,
here::here("Output",paste0("coeff_",nom,".csv")),
sep=";",dec=".",row=F)
write.table(randEff,here::here("Output",paste0("randEff_",nom,".csv")),
sep=";",dec=".",row=F)8.4.3 LEW
mydatatogoS=subset(mydatatogo,Param=="LE")
x1=data.frame(table(mydatatogoS$spp))
x1=subset(x1,Freq>0)
x1$weight_obs_spp=(1/x1$Freq)*(1/nrow(x1))
mydatatogoS=merge(mydatatogoS[,1:48],x1[,c(1,3)],by.x="spp",by.y="Var1")
my.cluster2 <- parallel::makeCluster(
nCPU,
type = "PSOCK"
)
clusterExport(cl = my.cluster2, varlist = c("s1","nB","glmmTMB","mydatatogoS","gam_velxGD1","ranef"))
#register cluster
doParallel::registerDoParallel(cl = my.cluster2)
ex=pblapply(1:nB, function(i) {
gc(reset=T)
print(i)
set.seed(s1+i)
n1=sample(1:nrow(mydatatogoS),rep=T,size=nrow(mydatatogoS))
mydatatogo2=mydatatogoS[n1,]
x1=data.frame(table(mydatatogo2$spp))
x1=subset(x1,Freq>0)
x1$weight_obs_spp=(1/x1$Freq)*(1/nrow(x1))
mydatatogo2=merge(mydatatogo2[,1:48],x1[,c(1,3)],by.x="spp",by.y="Var1")
gamX1 <- glmmTMB(gam_velxGD1,
family = Gamma(link = "log"),
weights = weight_obs_spp,
data = mydatatogo2)
r2=round(MuMIn::r.squaredGLMM(gamX1),2)
r2=data.frame(r2m=r2[1,1],
r2c=r2[1,2],
moy_GD=mean(mydatatogo2$GD),
sd_GD=sd(mydatatogo2$GD),
moy_VC=mean(mydatatogo2$vel_abs),
sd_VC=sd(mydatatogo2$vel_abs),
nB=i)
coeff=data.frame(summary(gamX1)$coeff$cond,
var=row.names(summary(gamX1)$coeff$cond),
nB=i)
randEff=as.data.frame(ranef(gamX1,condVar=T))
randEff$nB=i
return(list(r2,coeff,randEff))
}, cl = my.cluster2)
parallel::stopCluster(cl = my.cluster2)
nom="LEW"
randEff=rlist::list.rbind(lapply(ex,"[[",3))
coeff=rlist::list.rbind(lapply(ex,"[[",2))
r2=rlist::list.rbind(lapply(ex,"[[",1))
write.table(r2,
here::here("Output",paste0("R2_",nom,".csv")),
sep=";",dec=".",row=F)
write.table(coeff,
here::here("Output",paste0("coeff_",nom,".csv")),
sep=";",dec=".",row=F)
write.table(randEff,here::here("Output",paste0("randEff_",nom,".csv")),
sep=";",dec=".",row=F)8.4.4 TEW
mydatatogoS=subset(mydatatogo,Param=="TE")
x1=data.frame(table(mydatatogoS$spp))
x1=subset(x1,Freq>0)
x1$weight_obs_spp=(1/x1$Freq)*(1/nrow(x1))
mydatatogoS=merge(mydatatogoS[,1:48],x1[,c(1,3)],by.x="spp",by.y="Var1")
my.cluster2 <- parallel::makeCluster(
nCPU,
type = "PSOCK"
)
clusterExport(cl = my.cluster2, varlist = c("s1","nB","glmmTMB","mydatatogoS","gam_velxGD1","ranef"))
#register cluster
doParallel::registerDoParallel(cl = my.cluster2)
ex=pblapply(1:nB, function(i) {
gc(reset=T)
print(i)
set.seed(s1+i)
n1=sample(1:nrow(mydatatogoS),rep=T,size=nrow(mydatatogoS))
mydatatogo2=mydatatogoS[n1,]
x1=data.frame(table(mydatatogo2$spp))
x1=subset(x1,Freq>0)
x1$weight_obs_spp=(1/x1$Freq)*(1/nrow(x1))
mydatatogo2=merge(mydatatogo2[,1:48],x1[,c(1,3)],by.x="spp",by.y="Var1")
gamX1 <- glmmTMB(gam_velxGD1,
family = Gamma(link = "log"),
weights = weight_obs_spp,
data = mydatatogo2)
r2=round(MuMIn::r.squaredGLMM(gamX1),2)
r2=data.frame(r2m=r2[1,1],
r2c=r2[1,2],
moy_GD=mean(mydatatogo2$GD),
sd_GD=sd(mydatatogo2$GD),
moy_VC=mean(mydatatogo2$vel_abs),
sd_VC=sd(mydatatogo2$vel_abs),
nB=i)
coeff=data.frame(summary(gamX1)$coeff$cond,
var=row.names(summary(gamX1)$coeff$cond),
nB=i)
randEff=as.data.frame(ranef(gamX1,condVar=T))
randEff$nB=i
return(list(r2,coeff,randEff))
}, cl = my.cluster2)
parallel::stopCluster(cl = my.cluster2)
nom="TEW"
randEff=rlist::list.rbind(lapply(ex,"[[",3))
coeff=rlist::list.rbind(lapply(ex,"[[",2))
r2=rlist::list.rbind(lapply(ex,"[[",1))
write.table(r2,
here::here("Output",paste0("R2_",nom,".csv")),
sep=";",dec=".",row=F)
write.table(coeff,
here::here("Output",paste0("coeff_",nom,".csv")),
sep=";",dec=".",row=F)
write.table(randEff,here::here("Output",paste0("randEff_",nom,".csv")),
sep=";",dec=".",row=F)8.4.5 CEW
mydatatogoS=subset(mydatatogo,Param=="O")
x1=data.frame(table(mydatatogoS$spp))
x1=subset(x1,Freq>0)
x1$weight_obs_spp=(1/x1$Freq)*(1/nrow(x1))
mydatatogoS=merge(mydatatogoS[,1:48],x1[,c(1,3)],by.x="spp",by.y="Var1")
my.cluster2 <- parallel::makeCluster(
nCPU,
type = "PSOCK"
)
clusterExport(cl = my.cluster2, varlist = c("s1","nB","glmmTMB","mydatatogoS","gam_velxGD1","ranef"))
#register cluster
doParallel::registerDoParallel(cl = my.cluster2)
ex=pblapply(1:nB, function(i) {
gc(reset=T)
print(i)
set.seed(s1+i)
n1=sample(1:nrow(mydatatogoS),rep=T,size=nrow(mydatatogoS))
mydatatogo2=mydatatogoS[n1,]
x1=data.frame(table(mydatatogo2$spp))
x1=subset(x1,Freq>0)
x1$weight_obs_spp=(1/x1$Freq)*(1/nrow(x1))
mydatatogo2=merge(mydatatogo2[,1:48],x1[,c(1,3)],by.x="spp",by.y="Var1")
gamX1 <- glmmTMB(gam_velxGD1,
family = Gamma(link = "log"),
weights = weight_obs_spp,
data = mydatatogo2)
r2=round(MuMIn::r.squaredGLMM(gamX1),2)
r2=data.frame(r2m=r2[1,1],
r2c=r2[1,2],
moy_GD=mean(mydatatogo2$GD),
sd_GD=sd(mydatatogo2$GD),
moy_VC=mean(mydatatogo2$vel_abs),
sd_VC=sd(mydatatogo2$vel_abs),
nB=i)
coeff=data.frame(summary(gamX1)$coeff$cond,
var=row.names(summary(gamX1)$coeff$cond),
nB=i)
randEff=as.data.frame(ranef(gamX1,condVar=T))
randEff$nB=i
return(list(r2,coeff,randEff))
}, cl = my.cluster2)
parallel::stopCluster(cl = my.cluster2)
nom="CEW"
randEff=rlist::list.rbind(lapply(ex,"[[",3))
coeff=rlist::list.rbind(lapply(ex,"[[",2))
r2=rlist::list.rbind(lapply(ex,"[[",1))
write.table(r2,
here::here("Output",paste0("R2_",nom,".csv")),
sep=";",dec=".",row=F)
write.table(coeff,
here::here("Output",paste0("coeff_",nom,".csv")),
sep=";",dec=".",row=F)
write.table(randEff,here::here("Output",paste0("randEff_",nom,".csv")),
sep=";",dec=".",row=F)8.5 Collect data
resR2_ok_file <- here::here("Output","summary_R2.csv")
res_ok_file <- here::here("Output","summary_coeff.csv")
if(all(file.exists(resR2_ok_file,res_ok_file))){
resR2_ok_file <- read.table(resR2_ok_file, sep=";", dec=".")
res_ok_file <- read.table(res_ok_file, sep=";", dec=".")
} else {
mod=c("allW","CEW","TEW","LEW","allW2")
a=1
s1=10 #the seed to make reproducible random staff
for(i in 1:length(mod)){
print(a)
r2=read.csv2(here::here("Output",paste0("R2_",mod[i],".csv")),
sep=";",dec=".")
#r2$nB=1:nrow(r2)
r2a=subset(r2,is.na(r2m)==F)
print(paste(mod[i],': ',nrow(r2a),sep=""))
if(nrow(r2a)>=nB){ #random selection of nB bootstrap model (criteria: model convergence)
set.seed(s1)
r2a=r2a[sample(1:nrow(r2a),size=nB,replace=F),]
}
if(nrow(r2a)<nB){
print(paste0("sampling size is less than ",nB," bootstraps"))
}else{
resR2=data.frame(type="r2m",moy=mean(r2a$r2m),sd=sd(r2a$r2m),median=median(r2a$r2m),q025=quantile(r2a$r2m,probs=0.025),q975=quantile(r2a$r2m,probs=0.975))
resR2=rbind(resR2,data.frame(type="r2c",moy=mean(r2a$r2c),sd=sd(r2a$r2c),median=median(r2a$r2c),q025=quantile(r2a$r2c,probs=0.025),q975=quantile(r2a$r2c,probs=0.975)))
#resR2$n_sing=nrow(subset(r2a,sing==T))
resR2$model=mod[i]
c1=read.csv2(here::here("Output",paste("coeff_",mod[i],".csv",sep="")),
sep=";",dec=".",h=T)
#c1$nB=rep(1:nrow(r2),each=nrow(c1)/6000)
c1=merge(c1,data.frame(nB=r2a$nB),by.x="nB",by.y="nB")
rr3=c1
res=data.frame(var=names(tapply(rr3$Estimate,rr3$var,mean)),moy=tapply(rr3$Estimate,rr3$var,mean),sd=tapply(rr3$Estimate,rr3$var,sd),median=tapply(rr3$Estimate,rr3$var,median),q025=tapply(rr3$Estimate,rr3$var,quantile,probs=0.025),
p975=tapply(rr3$Estimate,rr3$var,quantile,probs=0.975), pv.inf0=tapply(rr3$Estimate,rr3$var,test2,mu=0),pv.sup0=tapply(rr3$Estimate,rr3$var,test1,mu=0))
res$model=mod[i]
if(a==1){
resR2_ok=resR2
res_ok=res
}
else{
resR2_ok=rbind(resR2_ok,resR2)
res_ok=rbind(res_ok,res)
}
a=a+1
}
}
write.table(resR2_ok,
here::here("Output","summary_R2.csv"),
sep=";",dec=".",row=F)
#R2 statistic computed from 5000 boostrap model:
#type= R2 type (r2m= marginal R2; r2c= conditional R2)
#moy= mean R2 value
#sd= R2 standard deviation value
#median= median R2 value
#q025= 2.5% quantile of the bootstrap distribution
#q975= 97.5% quantile of the bootstrap distribution
#n_sing= number of singular model
#model= model name abbreviation
write.table(res_ok,
here::here("Output","summary_coeff.csv"),
sep=";",dec=".",row=F)
#coefficient statistics
#Var = Term of the model
#moy = mean effect computed from the 5000 bootstrap model
#sd = standard deviation computed from the 5000 bootstrap model
#med = median effect computed from the 5000 bootstrap model
#q025 = 2.5% quantile of the bootstrap coefficient distribution
#q975 = 97.5% quantile of the bootstrap coefficient distribution
#pv.sup0 = p-value of the bootstrap test testing the significance of the effect (H0= coefficient equal to 0; H1= coefficient greater than 0)
#pv.inf0 = p-value of the bootstrap test testing the significance of the effect (H0= coefficient equal to 0; H1= coefficient lower than 0)
#model= model name abbreviation
# #significant positive effect
# subset(res_ok,model=="allW" & pv.sup0<0.05)
# subset(res_ok,model=="CEW" & pv.sup0<0.05)
# subset(res_ok,model=="LEW" & pv.sup0<0.05)
# subset(res_ok,model=="TEW" & pv.sup0<0.05)
# subset(res_ok,model=="allW2" & pv.sup0<0.05)
#
# #significant negative effect
# subset(res_ok,model=="allW" & pv.inf0<0.05)
# subset(res_ok,model=="CEW" & pv.inf0<0.05)
# subset(res_ok,model=="LEW" & pv.inf0<0.05)
# subset(res_ok,model=="TEW" & pv.inf0<0.05)
# subset(res_ok,model=="allW2" & pv.inf0<0.05)
#
# #non-significant negative effect
# subset(res_ok,model=="allW" & pv.inf0>=0.05 & pv.sup0>=0.05)
# subset(res_ok,model=="CEW" & pv.inf0>=0.05 & pv.sup0>=0.05)
# subset(res_ok,model=="LEW" & pv.inf0>=0.05 & pv.sup0>=0.05)
# subset(res_ok,model=="TEW" & pv.inf0>=0.05 & pv.sup0>=0.05)
# subset(res_ok,model=="allW2" & pv.inf0>=0.05 & pv.sup0>=0.05)
}8.5.1 Test climate velocity value for which GD is significant
8.5.1.1 allW model
res <- here::here("Output","summary_GDeff_allW.csv")
if(!file.exists(res)){
c1=read.csv2(here::here("Output","coeff_allW.csv"),
sep=";",dec=".",h=T)
dsel=mydatatogo
v1=seq(round(min(dsel$vel_abs),1),round(max(dsel$vel_abs),1),by=0.1)
my.cluster2 <- parallel::makeCluster(
nCPU,
type = "PSOCK"
)
clusterExport(cl = my.cluster2, varlist = c("c1","dsel","v1","test1","test2"))
#register cluster
doParallel::registerDoParallel(cl = my.cluster2)
ex=pblapply(1:nB,function(i){
gc(reset=T)
print(i)
c1a=subset(c1,nB==i)
eff=c1a$Estimate[c1a$var=="scale(GD)"]+(c1a$Estimate[c1a$var=="scale(vel_abs):scale(GD)"]*((v1-mean(dsel$vel_abs))/sd(dsel$vel_abs)))
tmp=data.frame(vel_abs=v1,GDeff=eff,nB=i)
return(tmp)
}, cl = my.cluster2)
parallel::stopCluster(cl = my.cluster2)
ex <- do.call("rbind",ex)
res=data.frame(vel_abs=names(tapply(ex$GDeff,ex$vel_abs,mean)),
moy=tapply(ex$GDeff,ex$vel_abs,mean),
sd=tapply(ex$GDeff,ex$vel_abs,sd),
median=tapply(ex$GDeff,ex$vel_abs,median),
q025=tapply(ex$GDeff,ex$vel_abs,quantile,probs=0.025),
p975=tapply(ex$GDeff,ex$vel_abs,quantile,probs=0.975),
pv.inf0=tapply(ex$GDeff,ex$vel_abs,test2,mu=0),
pv.sup0=tapply(ex$GDeff,ex$vel_abs,test1,mu=0))
write.table(res,
here::here("Output","summary_GDeff_allW.csv"),
sep=";",dec=".",row=F)
# #significant positive effect
# subset(res,pv.sup0<0.05)
# #significant negative effect
# subset(res,pv.inf0<0.05)
# #non-significant negative effect
# subset(res,pv.inf0>=0.05 & pv.sup0>=0.05)
}8.5.1.2 CEW model
res <- here::here("Output","summary_GDeff_CEW.csv")
if(!file.exists(res)){
c1=read.csv2(here::here("Output","coeff_CEW.csv"),
sep=";",dec=".",h=T)
dsel=subset(mydatatogo,Param=="O")
v1=seq(round(min(dsel$vel_abs),1),round(max(dsel$vel_abs),1),by=0.1)
my.cluster2 <- parallel::makeCluster(
nCPU,
type = "PSOCK"
)
clusterExport(cl = my.cluster2, varlist = c("c1","dsel","v1","test1","test2"))
#register cluster
doParallel::registerDoParallel(cl = my.cluster2)
ex=pblapply(1:nB,function(i){
gc(reset=T)
print(i)
c1a=subset(c1,nB==i)
eff=c1a$Estimate[c1a$var=="scale(GD)"]+(c1a$Estimate[c1a$var=="scale(vel_abs):scale(GD)"]*((v1-mean(dsel$vel_abs))/sd(dsel$vel_abs)))
tmp=data.frame(vel_abs=v1,GDeff=eff,nB=i)
return(tmp)
}, cl = my.cluster2)
parallel::stopCluster(cl = my.cluster2)
ex <- do.call("rbind",ex)
res=data.frame(vel_abs=names(tapply(ex$GDeff,ex$vel_abs,mean)),
moy=tapply(ex$GDeff,ex$vel_abs,mean),
sd=tapply(ex$GDeff,ex$vel_abs,sd),
median=tapply(ex$GDeff,ex$vel_abs,median),
q025=tapply(ex$GDeff,ex$vel_abs,quantile,probs=0.025),
p975=tapply(ex$GDeff,ex$vel_abs,quantile,probs=0.975),
pv.inf0=tapply(ex$GDeff,ex$vel_abs,test2,mu=0),
pv.sup0=tapply(ex$GDeff,ex$vel_abs,test1,mu=0))
write.table(res,
here::here("Output","summary_GDeff_CEW.csv"),
sep=";",dec=".",row=F)
# #significant positive effect
# subset(res,pv.sup0<0.05)
# #significant negative effect
# subset(res,pv.inf0<0.05)
# #non-significant negative effect
# subset(res,pv.inf0>=0.05 & pv.sup0>=0.05)
}8.5.1.3 LEW model
res <- here::here("Output","summary_GDeff_LEW.csv")
if(!file.exists(res)){
c1=read.csv2(here::here("Output","coeff_LEW.csv"),
sep=";",dec=".",h=T)
dsel=subset(mydatatogo,Param=="LE")
v1=seq(round(min(dsel$vel_abs),1),round(max(dsel$vel_abs),1),by=0.1)
my.cluster2 <- parallel::makeCluster(
nCPU,
type = "PSOCK"
)
clusterExport(cl = my.cluster2, varlist = c("c1","dsel","v1","test1","test2"))
#register cluster
doParallel::registerDoParallel(cl = my.cluster2)
ex=pblapply(1:nB,function(i){
gc(reset=T)
print(i)
c1a=subset(c1,nB==i)
eff=c1a$Estimate[c1a$var=="scale(GD)"]+(c1a$Estimate[c1a$var=="scale(vel_abs):scale(GD)"]*((v1-mean(dsel$vel_abs))/sd(dsel$vel_abs)))
tmp=data.frame(vel_abs=v1,GDeff=eff,nB=i)
return(tmp)
}, cl = my.cluster2)
parallel::stopCluster(cl = my.cluster2)
ex <- do.call("rbind",ex)
res=data.frame(vel_abs=names(tapply(ex$GDeff,ex$vel_abs,mean)),
moy=tapply(ex$GDeff,ex$vel_abs,mean),
sd=tapply(ex$GDeff,ex$vel_abs,sd),
median=tapply(ex$GDeff,ex$vel_abs,median),
q025=tapply(ex$GDeff,ex$vel_abs,quantile,probs=0.025),
p975=tapply(ex$GDeff,ex$vel_abs,quantile,probs=0.975),
pv.inf0=tapply(ex$GDeff,ex$vel_abs,test2,mu=0),
pv.sup0=tapply(ex$GDeff,ex$vel_abs,test1,mu=0))
write.table(res,
here::here("Output","summary_GDeff_LEW.csv"),
sep=";",dec=".",row=F)
# #significant positive effect
# subset(res,pv.sup0<0.05)
# #significant negative effect
# subset(res,pv.inf0<0.05)
# #non-significant negative effect
# subset(res,pv.inf0>=0.05 & pv.sup0>=0.05)
}8.5.1.4 TEW model
res <- here::here("Output","summary_GDeff_TEW.csv")
if(!file.exists(res)){
c1=read.csv2(here::here("Output","coeff_TEW.csv"),
sep=";",dec=".",h=T)
dsel=subset(mydatatogo,Param=="TE")
v1=seq(round(min(dsel$vel_abs),1),round(max(dsel$vel_abs),1),by=0.1)
my.cluster2 <- parallel::makeCluster(
nCPU,
type = "PSOCK"
)
clusterExport(cl = my.cluster2, varlist = c("c1","dsel","v1","test1","test2"))
#register cluster
doParallel::registerDoParallel(cl = my.cluster2)
ex=pblapply(1:nB,function(i){
gc(reset=T)
print(i)
c1a=subset(c1,nB==i)
eff=c1a$Estimate[c1a$var=="scale(GD)"]+(c1a$Estimate[c1a$var=="scale(vel_abs):scale(GD)"]*((v1-mean(dsel$vel_abs))/sd(dsel$vel_abs)))
tmp=data.frame(vel_abs=v1,GDeff=eff,nB=i)
return(tmp)
}, cl = my.cluster2)
parallel::stopCluster(cl = my.cluster2)
ex <- do.call("rbind",ex)
res=data.frame(vel_abs=names(tapply(ex$GDeff,ex$vel_abs,mean)),
moy=tapply(ex$GDeff,ex$vel_abs,mean),
sd=tapply(ex$GDeff,ex$vel_abs,sd),
median=tapply(ex$GDeff,ex$vel_abs,median),
q025=tapply(ex$GDeff,ex$vel_abs,quantile,probs=0.025),
p975=tapply(ex$GDeff,ex$vel_abs,quantile,probs=0.975),
pv.inf0=tapply(ex$GDeff,ex$vel_abs,test2,mu=0),
pv.sup0=tapply(ex$GDeff,ex$vel_abs,test1,mu=0))
write.table(res,
here::here("Output","summary_GDeff_TEW.csv"),
sep=";",dec=".",row=F)
# #significant positive effect
# subset(res,pv.sup0<0.05)
# #significant negative effect
# subset(res,pv.inf0<0.05)
# #non-significant negative effect
# subset(res,pv.inf0>=0.05 & pv.sup0>=0.05)
}8.5.1.5 allW2 model
resTE <- here::here("Output","summary_GDeff_allW2_TE.csv")
resCE <- here::here("Output","summary_GDeff_allW2_CE.csv")
resLE <- here::here("Output","summary_GDeff_allW2_LE.csv")
if(!all(file.exists(resTE,resCE,resLE))){
c1=read.csv2(here::here("Output","coeff_allW2.csv"),
sep=";",dec=".",h=T)
dsel=mydatatogo
v1=seq(round(min(dsel$vel_abs),1),round(max(dsel$vel_abs),1),by=0.1)
my.cluster2 <- parallel::makeCluster(
nCPU,
type = "PSOCK"
)
clusterExport(cl = my.cluster2, varlist = c("c1","dsel","v1","test1","test2"))
#register cluster
doParallel::registerDoParallel(cl = my.cluster2)
ex=pblapply(1:nB,function(i){
gc(reset=T)
print(i)
c1a=subset(c1,nB==i)
eff=c1a$Estimate[c1a$var=="scale(GD)"]+(c1a$Estimate[c1a$var=="scale(vel_abs):scale(GD)"]*((v1-mean(dsel$vel_abs))/sd(dsel$vel_abs)))
tmp=data.frame(vel_abs=v1,GDeff_CE=eff)
eff=c1a$Estimate[c1a$var=="scale(GD)"]+c1a$Estimate[c1a$var=="scale(GD):ParamLE"]+((c1a$Estimate[c1a$var=="scale(vel_abs):scale(GD)"]+c1a$Estimate[c1a$var=="scale(vel_abs):scale(GD):ParamLE"])*((v1-mean(dsel$vel_abs))/sd(dsel$vel_abs)))
tmp$GDeff_LE=eff
eff=c1a$Estimate[c1a$var=="scale(GD)"]+c1a$Estimate[c1a$var=="scale(GD):ParamTE"]+((c1a$Estimate[c1a$var=="scale(vel_abs):scale(GD)"]+c1a$Estimate[c1a$var=="scale(vel_abs):scale(GD):ParamTE"])*((v1-mean(dsel$vel_abs))/sd(dsel$vel_abs)))
tmp$GDeff_TE=eff
tmp$nB=i
return(tmp)
}, cl = my.cluster2)
parallel::stopCluster(cl = my.cluster2)
ex <- do.call("rbind",ex)
#summary statistics
resCE=data.frame(vel_abs=names(tapply(ex$GDeff_CE,ex$vel_abs,mean)),
moy=tapply(ex$GDeff_CE,ex$vel_abs,mean),
sd=tapply(ex$GDeff_CE,ex$vel_abs,sd),
median=tapply(ex$GDeff_CE,ex$vel_abs,median),
q025=tapply(ex$GDeff_CE,ex$vel_abs,quantile,probs=0.025),
p975=tapply(ex$GDeff_CE,ex$vel_abs,quantile,probs=0.975),
pv.inf0=tapply(ex$GDeff_CE,ex$vel_abs,test2,mu=0),
pv.sup0=tapply(ex$GDeff_CE,ex$vel_abs,test1,mu=0))
resLE=data.frame(vel_abs=names(tapply(ex$GDeff_LE,ex$vel_abs,mean)),
moy=tapply(ex$GDeff_LE,ex$vel_abs,mean),
sd=tapply(ex$GDeff_LE,ex$vel_abs,sd),
median=tapply(ex$GDeff_LE,ex$vel_abs,median),
q025=tapply(ex$GDeff_LE,ex$vel_abs,quantile,probs=0.025),
p975=tapply(ex$GDeff_LE,ex$vel_abs,quantile,probs=0.975),
pv.inf0=tapply(ex$GDeff_LE,ex$vel_abs,test2,mu=0),
pv.sup0=tapply(ex$GDeff_LE,ex$vel_abs,test1,mu=0))
resTE=data.frame(vel_abs=names(tapply(ex$GDeff_TE,ex$vel_abs,mean)),
moy=tapply(ex$GDeff_TE,ex$vel_abs,mean),
sd=tapply(ex$GDeff_TE,ex$vel_abs,sd),
median=tapply(ex$GDeff_TE,ex$vel_abs,median),
q025=tapply(ex$GDeff_TE,ex$vel_abs,quantile,probs=0.025),
p975=tapply(ex$GDeff_TE,ex$vel_abs,quantile,probs=0.975),
pv.inf0=tapply(ex$GDeff_TE,ex$vel_abs,test2,mu=0),
pv.sup0=tapply(ex$GDeff_TE,ex$vel_abs,test1,mu=0))
write.table(resTE,
here::here("Output","summary_GDeff_allW2_TE.csv"),
sep=";",dec=".",row=F)
write.table(resCE,
here::here("Output","summary_GDeff_allW2_CE.csv"),
sep=";",dec=".",row=F)
write.table(resLE,
here::here("Output","summary_GDeff_allW2_LE.csv"),
sep=";",dec=".",row=F)
# #significant positive effect
# subset(resCE,pv.sup0<0.05)
# #significant negative effect
# subset(resCE,pv.inf0<0.05)
# #non-significant negative effect
# subset(resCE,pv.inf0>=0.05 & pv.sup0>=0.05)
# #significant positive effect
# subset(resLE,pv.sup0<0.05)
# #significant negative effect
# subset(resLE,pv.inf0<0.05)
# #non-significant negative effect
# subset(resLE,pv.inf0>=0.05 & pv.sup0>=0.05)
# #significant positive effect
# subset(resTE,pv.sup0<0.05)
# #significant negative effect
# subset(resTE,pv.inf0<0.05)
# #non-significant negative effect
# subset(resTE,pv.inf0>=0.05 & pv.sup0>=0.05)
}8.5.2 Test GD value for which the climate velocity effect is significant
8.5.2.1 allW model
res <- here::here("Output","summary_VAeff_allW.csv")
if(!file.exists(res)){
c1=read.csv2(here::here("Output","coeff_allW.csv"),
sep=";",dec=".",h=T)
dsel=mydatatogo
v1=seq(round(min(dsel$GD),2),round(max(dsel$GD),2),by=0.01)
my.cluster2 <- parallel::makeCluster(
nCPU,
type = "PSOCK"
)
clusterExport(cl = my.cluster2, varlist = c("c1","dsel","v1","test1","test2"))
ex=pblapply(1:nB,function(i){
gc(reset=T)
print(i)
c1a=subset(c1,nB==i)
eff=c1a$Estimate[c1a$var=="scale(vel_abs)"]+(c1a$Estimate[c1a$var=="scale(vel_abs):scale(GD)"]*((v1-mean(dsel$GD))/sd(dsel$GD)))
tmp=data.frame(GD=v1,VAeff=eff,nB=i)
return(tmp)
}, cl = my.cluster2)
parallel::stopCluster(cl = my.cluster2)
ex <- do.call("rbind",ex)
res=data.frame(GD=names(tapply(ex$VAeff,ex$GD,mean)),
moy=tapply(ex$VAeff,ex$GD,mean),
sd=tapply(ex$VAeff,ex$GD,sd),
median=tapply(ex$VAeff,ex$GD,median),
q025=tapply(ex$VAeff,ex$GD,quantile,probs=0.025),
p975=tapply(ex$VAeff,ex$GD,quantile,probs=0.975),
pv.inf0=tapply(ex$VAeff,ex$GD,test2,mu=0),
pv.sup0=tapply(ex$VAeff,ex$GD,test1,mu=0))
write.table(res,
here::here("Output","summary_VAeff_allW.csv"),
sep=";",dec=".",row=F)
# #significant positive effect
# subset(res,pv.sup0<0.05)
# #significant negative effect
# subset(res,pv.inf0<0.05)
# #non-significant negative effect
# subset(res,pv.inf0>=0.05 & pv.sup0>=0.05)
}8.5.2.2 CEW model
res <- here::here("Output","summary_VAeff_CEW.csv")
if(!file.exists(res)){
c1=read.csv2(here::here("Output","coeff_CEW.csv"),
sep=";",dec=".",h=T)
dsel=subset(mydatatogo,Param=="O")
v1=seq(round(min(dsel$GD),2),round(max(dsel$GD),2),by=0.01)
my.cluster2 <- parallel::makeCluster(
nCPU,
type = "PSOCK"
)
clusterExport(cl = my.cluster2, varlist = c("c1","dsel","v1","test1","test2"))
ex=pblapply(1:nB,function(i){
gc(reset=T)
print(i)
c1a=subset(c1,nB==i)
eff=c1a$Estimate[c1a$var=="scale(vel_abs)"]+(c1a$Estimate[c1a$var=="scale(vel_abs):scale(GD)"]*((v1-mean(dsel$GD))/sd(dsel$GD)))
tmp=data.frame(GD=v1,VAeff=eff,nB=i)
return(tmp)
}, cl = my.cluster2)
parallel::stopCluster(cl = my.cluster2)
ex <- do.call("rbind",ex)
res=data.frame(GD=names(tapply(ex$VAeff,ex$GD,mean)),
moy=tapply(ex$VAeff,ex$GD,mean),
sd=tapply(ex$VAeff,ex$GD,sd),
median=tapply(ex$VAeff,ex$GD,median),
q025=tapply(ex$VAeff,ex$GD,quantile,probs=0.025),
p975=tapply(ex$VAeff,ex$GD,quantile,probs=0.975),
pv.inf0=tapply(ex$VAeff,ex$GD,test2,mu=0),
pv.sup0=tapply(ex$VAeff,ex$GD,test1,mu=0))
write.table(res,
here::here("Output","summary_VAeff_CEW.csv"),
sep=";",dec=".",row=F)
# #significant positive effect
# subset(res,pv.sup0<0.05)
# #significant negative effect
# subset(res,pv.inf0<0.05)
# #non-significant negative effect
# subset(res,pv.inf0>=0.05 & pv.sup0>=0.05)
}8.5.2.3 LEW model
res <- here::here("Output","summary_VAeff_LEW.csv")
if(!file.exists(res)){
c1=read.csv2(here::here("Output","coeff_LEW.csv"),
sep=";",dec=".",h=T)
dsel=subset(mydatatogo,Param=="LE")
v1=seq(round(min(dsel$GD),2),round(max(dsel$GD),2),by=0.01)
my.cluster2 <- parallel::makeCluster(
nCPU,
type = "PSOCK"
)
clusterExport(cl = my.cluster2, varlist = c("c1","dsel","v1","test1","test2"))
ex=pblapply(1:nB,function(i){
gc(reset=T)
print(i)
c1a=subset(c1,nB==i)
eff=c1a$Estimate[c1a$var=="scale(vel_abs)"]+(c1a$Estimate[c1a$var=="scale(vel_abs):scale(GD)"]*((v1-mean(dsel$GD))/sd(dsel$GD)))
tmp=data.frame(GD=v1,VAeff=eff,nB=i)
return(tmp)
}, cl = my.cluster2)
parallel::stopCluster(cl = my.cluster2)
ex <- do.call("rbind",ex)
res=data.frame(GD=names(tapply(ex$VAeff,ex$GD,mean)),
moy=tapply(ex$VAeff,ex$GD,mean),
sd=tapply(ex$VAeff,ex$GD,sd),
median=tapply(ex$VAeff,ex$GD,median),
q025=tapply(ex$VAeff,ex$GD,quantile,probs=0.025),
p975=tapply(ex$VAeff,ex$GD,quantile,probs=0.975),
pv.inf0=tapply(ex$VAeff,ex$GD,test2,mu=0),
pv.sup0=tapply(ex$VAeff,ex$GD,test1,mu=0))
write.table(res,
here::here("Output","summary_VAeff_LEW.csv"),
sep=";",dec=".",row=F)
# #significant positive effect
# subset(res,pv.sup0<0.05)
# #significant negative effect
# subset(res,pv.inf0<0.05)
# #non-significant negative effect
# subset(res,pv.inf0>=0.05 & pv.sup0>=0.05)
}8.5.2.4 TEW model
res <- here::here("Output","summary_VAeff_TEW.csv")
if(!file.exists(res)){
c1=read.csv2(here::here("Output","coeff_TEW.csv"),sep=";",dec=".",h=T)
dsel=subset(mydatatogo,Param=="TE")
v1=seq(round(min(dsel$GD),2),round(max(dsel$GD),2),by=0.01)
my.cluster2 <- parallel::makeCluster(
nCPU,
type = "PSOCK"
)
clusterExport(cl = my.cluster2, varlist = c("c1","dsel","v1","test1","test2"))
ex=pblapply(1:nB,function(i){
gc(reset=T)
print(i)
c1a=subset(c1,nB==i)
eff=c1a$Estimate[c1a$var=="scale(vel_abs)"]+(c1a$Estimate[c1a$var=="scale(vel_abs):scale(GD)"]*((v1-mean(dsel$GD))/sd(dsel$GD)))
tmp=data.frame(GD=v1,VAeff=eff,nB=i)
return(tmp)
}, cl = my.cluster2)
parallel::stopCluster(cl = my.cluster2)
ex <- do.call("rbind",ex)
res=data.frame(GD=names(tapply(ex$VAeff,ex$GD,mean)),
moy=tapply(ex$VAeff,ex$GD,mean),
sd=tapply(ex$VAeff,ex$GD,sd),
median=tapply(ex$VAeff,ex$GD,median),
q025=tapply(ex$VAeff,ex$GD,quantile,probs=0.025),
p975=tapply(ex$VAeff,ex$GD,quantile,probs=0.975),
pv.inf0=tapply(ex$VAeff,ex$GD,test2,mu=0),
pv.sup0=tapply(ex$VAeff,ex$GD,test1,mu=0))
write.table(res,
here::here("Output","summary_VAeff_TEW.csv"),
sep=";",dec=".",row=F)
# #significant positive effect
# subset(res,pv.sup0<0.05)
# #significant negative effect
# subset(res,pv.inf0<0.05)
# #non-significant negative effect
# subset(res,pv.inf0>=0.05 & pv.sup0>=0.05)
}8.5.2.5 allW2 model
resTE <- here::here("Output","summary_VAeff_allW2_TE.csv")
resCE <- here::here("Output","summary_VAeff_allW2_CE.csv")
resLE <- here::here("Output","summary_VAeff_allW2_LE.csv")
if(!all(file.exists(resTE,resCE,resLE))){
c1=read.csv2(here::here("Output","coeff_allW2.csv"),sep=";",dec=".",h=T)
dsel=mydatatogo
v1=seq(round(min(dsel$GD),3),round(max(dsel$GD),3),by=0.001)
my.cluster2 <- parallel::makeCluster(
nCPU,
type = "PSOCK"
)
clusterExport(cl = my.cluster2, varlist = c("c1","dsel","v1","test1","test2"))
ex=pblapply(1:nB,function(i){
gc(reset=T)
print(i)
c1a=subset(c1,nB==i)
eff=c1a$Estimate[c1a$var=="scale(vel_abs)"]+(c1a$Estimate[c1a$var=="scale(vel_abs):scale(GD)"]*((v1-mean(dsel$GD))/sd(dsel$GD)))
tmp=data.frame(GD=v1,VAeff_CE=eff)
eff=c1a$Estimate[c1a$var=="scale(vel_abs)"]+c1a$Estimate[c1a$var=="scale(vel_abs):ParamLE"]+((c1a$Estimate[c1a$var=="scale(vel_abs):scale(GD)"]+c1a$Estimate[c1a$var=="scale(vel_abs):scale(GD):ParamLE"])*((v1-mean(dsel$GD))/sd(dsel$GD)))
tmp$VAeff_LE=eff
eff=c1a$Estimate[c1a$var=="scale(vel_abs)"]+c1a$Estimate[c1a$var=="scale(vel_abs):ParamTE"]+((c1a$Estimate[c1a$var=="scale(vel_abs):scale(GD)"]+c1a$Estimate[c1a$var=="scale(vel_abs):scale(GD):ParamTE"])*((v1-mean(dsel$GD))/sd(dsel$GD)))
tmp$VAeff_TE=eff
tmp$nB=i
return(tmp)
}, cl = my.cluster2)
parallel::stopCluster(cl = my.cluster2)
ex <- do.call("rbind",ex)
#summary statistictics
resCE=data.frame(GD=names(tapply(ex$VAeff_CE,ex$GD,mean)),
moy=tapply(ex$VAeff_CE,ex$GD,mean),
sd=tapply(ex$VAeff_CE,ex$GD,sd),
median=tapply(ex$VAeff_CE,ex$GD,median),
q025=tapply(ex$VAeff_CE,ex$GD,quantile,probs=0.025),
p975=tapply(ex$VAeff_CE,ex$GD,quantile,probs=0.975),
pv.inf0=tapply(ex$VAeff_CE,ex$GD,test2,mu=0),
pv.sup0=tapply(ex$VAeff_CE,ex$GD,test1,mu=0))
resLE=data.frame(GD=names(tapply(ex$VAeff_LE,ex$GD,mean)),
moy=tapply(ex$VAeff_LE,ex$GD,mean),
sd=tapply(ex$VAeff_LE,ex$GD,sd),
median=tapply(ex$VAeff_LE,ex$GD,median),
q025=tapply(ex$VAeff_LE,ex$GD,quantile,probs=0.025),
p975=tapply(ex$VAeff_LE,ex$GD,quantile,probs=0.975),
pv.inf0=tapply(ex$VAeff_LE,ex$GD,test2,mu=0),
pv.sup0=tapply(ex$VAeff_LE,ex$GD,test1,mu=0))
resTE=data.frame(GD=names(tapply(ex$VAeff_TE,ex$GD,mean)),
moy=tapply(ex$VAeff_TE,ex$GD,mean),
sd=tapply(ex$VAeff_TE,ex$GD,sd),
median=tapply(ex$VAeff_TE,ex$GD,median),
q025=tapply(ex$VAeff_TE,ex$GD,quantile,probs=0.025),
p975=tapply(ex$VAeff_TE,ex$GD,quantile,probs=0.975),
pv.inf0=tapply(ex$VAeff_TE,ex$GD,test2,mu=0),
pv.sup0=tapply(ex$VAeff_TE,ex$GD,test1,mu=0))
write.table(resTE,
here::here("Output","summary_VAeff_allW2_TE.csv"),
sep=";",dec=".",row=F)
write.table(resLE,
here::here("Output","summary_VAeff_allW2_LE.csv"),
sep=";",dec=".",row=F)
write.table(resCE,
here::here("Output","summary_VAeff_allW2_CE.csv"),
sep=";",dec=".",row=F)
# #significant positive effect
# subset(resCE,pv.sup0<0.05)
# #significant negative effect
# subset(resCE,pv.inf0<0.05)
# #non-significant negative effect
# subset(resCE,pv.inf0>=0.05 & pv.sup0>=0.05)
# #significant positive effect
# subset(resLE,pv.sup0<0.05)
# #significant negative effect
# subset(resLE,pv.inf0<0.05)
# #non-significant negative effect
# subset(resLE,pv.inf0>=0.05 & pv.sup0>=0.05)
# #significant positive effect
# subset(resTE,pv.sup0<0.05)
# #significant negative effect
# subset(resTE,pv.inf0<0.05)
# #non-significant negative effect
# subset(resTE,pv.inf0>=0.05 & pv.sup0>=0.05)
}9 Tables goodness-of-fit
c1=read.csv2(here::here("Output","summary_R2.csv"),
sep=";",dec=".",h=T)
c1_m <- c1 %>%
select(moy, type, model) %>%
spread(type,moy)
c1_sd <- c1 %>%
select(sd, type, model) %>%
spread(type,sd)
c1_m$r2m <- paste0(round(c1_m$r2m,2), " (", round(c1_sd$r2m,2), ")")
c1_m$r2c <- paste0(round(c1_m$r2c,2), " (", round(c1_sd$r2c,2), ")")
DT::datatable(c1_m,
rownames = FALSE,
extensions = 'Buttons',
options = list(pageLength = 100,
dom = 'B',
buttons = c('excel'))) %>%
DT::formatStyle(columns = colnames(c1_m), `fontSize` = '80%')10 Tables coefficients
c1=read.csv2(here::here("Output","summary_coeff.csv"),
sep=";",dec=".",h=T)10.1 allW
mytab=subset(c1,model=="allW")
mytab <- mytab %>%
mutate(moy = paste(round(moy,2), paste0("(",round(sd,2),")")),
q025 = round(q025,3),
p975 = round(p975,3))
mytab <- mytab[-1,c(-3,-4,-7:-9)]
mytab$var[grep("scale(vel_abs)",mytab$var,fixed = TRUE)] <- gsub("scale(vel_abs)","Climate velocity",mytab$var[grep("scale(vel_abs)",mytab$var,fixed = TRUE)],fixed = TRUE)
mytab$var[grep("scale(GD)",mytab$var,fixed = TRUE)] <- gsub("scale(GD)","Genetic diversity",mytab$var[grep("scale(GD)",mytab$var,fixed = TRUE)],fixed = TRUE)
mytab$var[grep("Param",mytab$var,fixed = TRUE)] <- gsub("Param","RangePosition",mytab$var[grep("Param",mytab$var,fixed = TRUE)],fixed = TRUE)
names(mytab) <- c("Predictor variable", "Avg. coeff. (SD)", "Quant. 2.5%", "Quant. 9.75%")
mytab <- mytab[c(10,9,11,1:3,6:8,4,5),]
rownames(mytab) <- NULL
allWtab <- mytab
DT::datatable(allWtab,
rownames = FALSE,
extensions = 'Buttons',
options = list(pageLength = 100,
dom = 'B',
buttons = c('excel'))) %>%
DT::formatStyle(columns = colnames(allWtab), `fontSize` = '80%')10.2 TECELEW
# TE
mytab=subset(c1,model=="TEW")
mytab <- mytab %>%
mutate(moy = paste(round(moy,2), paste0("(",round(sd,3),")")),
q025 = round(q025,3),
p975 = round(p975,3),
`Range position` = "Trailing edge")
mytab <- mytab[-1,c(-3,-4,-7:-9)]
mytab$var[grep("scale(vel_abs)",mytab$var,fixed = TRUE)] <- gsub("scale(vel_abs)","Climate velocity",mytab$var[grep("scale(vel_abs)",mytab$var,fixed = TRUE)],fixed = TRUE)
mytab$var[grep("scale(GD)",mytab$var,fixed = TRUE)] <- gsub("scale(GD)","Genetic diversity",mytab$var[grep("scale(GD)",mytab$var,fixed = TRUE)],fixed = TRUE)
names(mytab) <- c("Predictor variable", "Avg. coeff. (SD)", "Quant. 2.5%", "Quant. 9.75%", "Range position")
mytab <- mytab[c(8,7,9,1:6),c(5,1:4)]
rownames(mytab) <- NULL
mytabTE <- mytab
# CE
mytab=subset(c1,model=="CEW")
mytab <- mytab %>%
mutate(moy = paste(round(moy,2), paste0("(",round(sd,3),")")),
q025 = round(q025,3),
p975 = round(p975,3),
`Range position` = "Centroid")
mytab <- mytab[-1,c(-3,-4,-7:-9)]
mytab$var[grep("scale(vel_abs)",mytab$var,fixed = TRUE)] <- gsub("scale(vel_abs)","Climate velocity",mytab$var[grep("scale(vel_abs)",mytab$var,fixed = TRUE)],fixed = TRUE)
mytab$var[grep("scale(GD)",mytab$var,fixed = TRUE)] <- gsub("scale(GD)","Genetic diversity",mytab$var[grep("scale(GD)",mytab$var,fixed = TRUE)],fixed = TRUE)
names(mytab) <- c("Predictor variable", "Avg. coeff. (SD)", "Quant. 2.5%", "Quant. 9.75%", "Range position")
mytab <- mytab[c(8,7,9,1:6),c(5,1:4)]
rownames(mytab) <- NULL
mytabCE <- mytab
# LE
mytab=subset(c1,model=="LEW")
mytab <- mytab %>%
mutate(moy = paste(round(moy,2), paste0("(",round(sd,3),")")),
q025 = round(q025,3),
p975 = round(p975,3),
`Range position` = "Leading edge")
mytab <- mytab[-1,c(-3,-4,-7:-9)]
mytab$var[grep("scale(vel_abs)",mytab$var,fixed = TRUE)] <- gsub("scale(vel_abs)","Climate velocity",mytab$var[grep("scale(vel_abs)",mytab$var,fixed = TRUE)],fixed = TRUE)
mytab$var[grep("scale(GD)",mytab$var,fixed = TRUE)] <- gsub("scale(GD)","Genetic diversity",mytab$var[grep("scale(GD)",mytab$var,fixed = TRUE)],fixed = TRUE)
names(mytab) <- c("Predictor variable", "Avg. coeff. (SD)", "Quant. 2.5%", "Quant. 9.75%", "Range position")
mytab <- mytab[c(8,7,9,1:6),c(5,1:4)]
rownames(mytab) <- NULL
mytabLE <- mytab
# Combine
TECELEWtab <- rbind(mytabTE,mytabCE,mytabLE)
DT::datatable(TECELEWtab,
rownames = FALSE,
extensions = 'Buttons',
options = list(pageLength = 100,
dom = 'B',
buttons = c('excel'))) %>%
DT::formatStyle(columns = colnames(TECELEWtab), `fontSize` = '80%')10.3 TECELEW2
mytab=subset(c1,model=="allW2")
mytab <- mytab %>%
mutate(moy = paste(round(moy,2), paste0("(",round(sd,2),")")),
q025 = round(q025,3),
p975 = round(p975,3))
mytab <- mytab[-1,c(-3,-4,-7:-9)]
mytab$var[grep("scale(vel_abs)",mytab$var,fixed = TRUE)] <- gsub("scale(vel_abs)","Climate velocity",mytab$var[grep("scale(vel_abs)",mytab$var,fixed = TRUE)],fixed = TRUE)
mytab$var[grep("scale(GD)",mytab$var,fixed = TRUE)] <- gsub("scale(GD)","Genetic diversity",mytab$var[grep("scale(GD)",mytab$var,fixed = TRUE)],fixed = TRUE)
mytab$var[grep("Param",mytab$var,fixed = TRUE)] <- gsub("Param","RangePosition",mytab$var[grep("Param",mytab$var,fixed = TRUE)],fixed = TRUE)
names(mytab) <- c("Predictor variable", "Avg. coeff. (SD)", "Quant. 2.5%", "Quant. 9.75%")
mytab <- mytab[c(9:14,15:17,1:8),]
rownames(mytab) <- NULL
TECELEW2tab <- mytab
DT::datatable(TECELEW2tab,
rownames = FALSE,
extensions = 'Buttons',
options = list(pageLength = 100,
dom = 'B',
buttons = c('excel'))) %>%
DT::formatStyle(columns = colnames(TECELEW2tab), `fontSize` = '80%')11 Graphics coefficients
11.1 allW
ggplot(allWtab %>%
rowwise() %>%
mutate(`Avg. coeff.` = as.numeric(strsplit(`Avg. coeff. (SD)`, " ")[[1]][1]),
`Predictor variable` = factor(`Predictor variable`, levels = `Predictor variable`),
sig = ifelse((`Quant. 2.5%` > 0 & `Quant. 9.75%` < 0) | (`Quant. 2.5%` < 0 & `Quant. 9.75%` > 0),"n","y")),
aes(x = `Avg. coeff.`, y = `Predictor variable`, color = sig)) +
scale_color_manual(values = c("y" = "black", "n" = "gray")) +
ylim(rev(allWtab$`Predictor variable`)) +
labs(x="Avg. Estimate\n2.5% - 9.75% quantiles",y="")+
geom_point(size=2) +
geom_vline(xintercept = 0, linetype = "dashed") +
geom_linerange(aes(xmin = `Quant. 2.5%`, xmax = `Quant. 9.75%`)) +
theme_classic() +
theme(legend.position = "none")
# focus on ecology
selected_vars <- allWtab$`Predictor variable`[1:3]
ggplot(allWtab %>%
rowwise() %>%
mutate(`Avg. coeff.` = as.numeric(strsplit(`Avg. coeff. (SD)`, " ")[[1]][1]),
`Predictor variable` = factor(`Predictor variable`, levels = `Predictor variable`),
sig = ifelse((`Quant. 2.5%` > 0 & `Quant. 9.75%` < 0) | (`Quant. 2.5%` < 0 & `Quant. 9.75%` > 0),"n","y")) %>%
filter(`Predictor variable` %in% selected_vars),
aes(x = `Avg. coeff.`, y = `Predictor variable`, color = sig)) +
scale_color_manual(values = c("y" = "black", "n" = "gray")) +
labs(x="Avg. Estimate\n2.5% - 9.75% quantiles",y="")+
ylim(rev(selected_vars)) +
geom_point(size=2) +
geom_vline(xintercept = 0, linetype = "dashed") +
geom_linerange(aes(xmin = `Quant. 2.5%`, xmax = `Quant. 9.75%`)) +
theme_classic() +
theme(legend.position = "none")
11.2 TECELEW
ggplot(TECELEWtab %>%
rowwise() %>%
mutate(`Avg. coeff.` = as.numeric(strsplit(`Avg. coeff. (SD)`, " ")[[1]][1]),
`Predictor variable` = factor(`Predictor variable`, levels = `Predictor variable`),
sig = ifelse((`Quant. 2.5%` > 0 & `Quant. 9.75%` < 0) | (`Quant. 2.5%` < 0 & `Quant. 9.75%` > 0),"n","y")),
aes(x = `Avg. coeff.`, y = `Predictor variable`, color = sig)) +
scale_color_manual(values = c("y" = "black", "n" = "gray")) +
ylim(unique(rev(TECELEWtab$`Predictor variable`))) +
labs(x="Avg. Estimate\n2.5% - 9.75% quantiles",y="")+
geom_point(size=2) +
geom_vline(xintercept = 0, linetype = "dashed") +
geom_linerange(aes(xmin = `Quant. 2.5%`, xmax = `Quant. 9.75%`)) +
theme_classic()+
facet_grid(.~`Range position`)+
theme(legend.position = "none")
# focus on ecology
selected_vars <- TECELEWtab$`Predictor variable`[1:3]
ggplot(TECELEWtab %>%
rowwise() %>%
mutate(`Avg. coeff.` = as.numeric(strsplit(`Avg. coeff. (SD)`, " ")[[1]][1]),
`Predictor variable` = factor(`Predictor variable`, levels = `Predictor variable`),
sig = ifelse((`Quant. 2.5%` > 0 & `Quant. 9.75%` < 0) | (`Quant. 2.5%` < 0 & `Quant. 9.75%` > 0),"n","y")) %>%
filter(`Predictor variable` %in% selected_vars),
aes(x = `Avg. coeff.`, y = `Predictor variable`, color = sig)) +
scale_color_manual(values = c("y" = "black", "n" = "gray")) +
ylim(rev(selected_vars)) +
labs(x="Avg. Estimate\n2.5% - 9.75% quantiles",y="")+
geom_point(size=2) +
geom_vline(xintercept = 0, linetype = "dashed") +
geom_linerange(aes(xmin = `Quant. 2.5%`, xmax = `Quant. 9.75%`)) +
theme_classic()+
theme(strip.background = element_blank())+
facet_wrap(.~`Range position`,ncol = 1)+
theme(legend.position = "none")
11.3 TECELEW2
ggplot(TECELEW2tab %>%
rowwise() %>%
mutate(`Avg. coeff.` = as.numeric(strsplit(`Avg. coeff. (SD)`, " ")[[1]][1]),
`Predictor variable` = factor(`Predictor variable`, levels = `Predictor variable`),
sig = ifelse((`Quant. 2.5%` > 0 & `Quant. 9.75%` < 0) | (`Quant. 2.5%` < 0 & `Quant. 9.75%` > 0),"n","y")),
aes(x = `Avg. coeff.`, y = `Predictor variable`, color = sig)) +
scale_color_manual(values = c("y" = "black", "n" = "gray")) +
ylim(unique(rev(TECELEW2tab$`Predictor variable`))) +
labs(x="Avg. Estimate\n2.5% - 9.75% quantiles",y="")+
geom_point(size=2) +
geom_vline(xintercept = 0, linetype = "dashed") +
geom_linerange(aes(xmin = `Quant. 2.5%`, xmax = `Quant. 9.75%`)) +
theme_classic() +
theme(legend.position = "none")
# focus on ecology
selected_vars <- TECELEW2tab$`Predictor variable`[1:9]
ggplot(TECELEW2tab %>%
rowwise() %>%
mutate(`Avg. coeff.` = as.numeric(strsplit(`Avg. coeff. (SD)`, " ")[[1]][1]),
`Predictor variable` = factor(`Predictor variable`, levels = `Predictor variable`),
sig = ifelse((`Quant. 2.5%` > 0 & `Quant. 9.75%` < 0) | (`Quant. 2.5%` < 0 & `Quant. 9.75%` > 0),"n","y")) %>%
filter(`Predictor variable` %in% selected_vars),
aes(x = `Avg. coeff.`, y = `Predictor variable`, color = sig)) +
scale_color_manual(values = c("y" = "black", "n" = "gray")) +
ylim(rev(selected_vars)) +
labs(x="Avg. Estimate\n2.5% - 9.75% quantiles",y="")+
geom_point(size=2) +
geom_vline(xintercept = 0, linetype = "dashed") +
geom_linerange(aes(xmin = `Quant. 2.5%`, xmax = `Quant. 9.75%`)) +
theme_classic()+
theme(strip.background = element_blank())+
theme(legend.position = "none")
12 Graphics effect GD
# GD_graph_params <-
# list(geom_line(),
# labs(x = "Genetic diversity",
# y = bquote('Absolute velocity of range shift '(km.yr^1)),
# color = bquote('Absolute velocity of isotherm shifts '(km.yr^1))),
# scale_color_gradientn(colors = c("#00AFBB", "#E7B800", "#FC4E07"),
# guide = guide_colourbar(title.position = "right",
# barwidth = 1, barheight = 10),
# limits=c(0, 7),
# breaks=c(seq(0,7,by=1))),
# theme_classic(),
# scale_x_continuous(breaks=c(c(seq(0,0.05,by=.01))),
# labels=c("0",as.character(seq(0.01,0.05,by=.01))),
# limits=c(0,0.055), expand=c(0,0)),
# scale_y_continuous(breaks=c(seq(0,6,by=2)), limits=c(0,6), expand=c(0,0)),
# theme(legend.title = element_text(angle = -90),
# legend.title.align = 0.5))
rampPallFunGD <- function(x){
tmp <- colorRampPalette(colors = c("#00AFBB", "#E7B800", "#FC4E07"))(length(x))
pos <- which(x==0)
tmp[pos] <- "lightgray"
return(tmp)
}
GD_graph_params <-
list(geom_line(size = 1),
labs(x = "Genetic diversity",
y = bquote('Absolute velocity of range shift '(km.yr^1)),
color = bquote('Absolute velocity of isotherm shifts '(km.yr^1))),
theme_classic(),
scale_x_continuous(n.breaks = 5,
limits=c(0,0.055), expand=c(0,0),
labels=c("0",as.character(seq(0.01,0.05,length.out=5)))),
scale_y_continuous(n.breaks = 4,
limits=c(0,6),
expand=c(0,0)),
theme(legend.title = element_text(angle = -90),
legend.title.align = 0.5))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: The `legend.title.align` argument of `theme()` is deprecated as of ggplot2
## 3.5.0.
## ℹ Please use theme(legend.title = element_text(hjust)) instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.12.1 AllW
res=read.csv2(here::here("Output","summary_GDeff_allW.csv"),
sep=";",dec=".",h=T)
c1=read.csv2(here::here("Output","summary_coeff.csv"),
sep=";",dec=".",h=T)
c1a=subset(c1,model=="allW")
int=c1a$moy[c1a$var=="(Intercept)"]
dsel=mydatatogo
v1=seq(round(min(dsel$vel_abs),1),round(max(dsel$vel_abs),1),by=0.1)
v2=seq(round(min(dsel$GD),4),round(max(dsel$GD),4),length.out=nrow(mydatatogo))
for(i in 1:nrow(res)){
pred1=exp(int+res$moy[i]*((v2-mean(mydatatogo$GD))/sd(mydatatogo$GD))+c1a$moy[c1a$var=="LogExtent"]*mean(mydatatogo$LogExtent)+c1a$moy[c1a$var=="LogNtempUnits"]*mean(mydatatogo$LogNtempUnits))
pred1=data.frame(GD=v2,pred1,
vel_abs=res$vel_abs[i],
sig = ifelse(res$pv.inf0[i]<.05 | res$pv.sup0[i]<.05,1,0),
lower = res$q025[i],
upper = res$p975[i])
if(i==1){
pred2=pred1
}else{
pred2=rbind(pred2,pred1)
}
}
res$ID=as.character(round(res$vel_abs,1))
pred2$ID=as.character(round(pred2$vel_abs,1))
resX=subset(res,pv.inf0<0.05)
if(nrow(resX>0)){
selN=round(seq(min(resX$vel_abs),max(resX$vel_abs),le=5),1)
predX=merge(pred2,data.frame(ID=as.character(selN),signifN=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX=pred2
predX$signifN=NA
}
resX=subset(res,pv.sup0<0.05)
if(nrow(resX>0)){
selP=round(seq(min(resX$vel_abs),max(resX$vel_abs),le=5),1)
predX2=merge(predX,data.frame(ID=as.character(selP),signifP=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX2=predX
predX2$signifP=NA
}
predX1=subset(predX2,signifN==1 | signifP==1)
sigcolor <- predX2 %>% group_by(vel_abs) %>% summarise(x = mean(sig))
gg1 <- ggplot(data=predX1, aes(GD, pred1, color = vel_abs, group=ID)) +
GD_graph_params +
scale_color_gradientn(colors = rampPallFunGD(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 8),
limits=c(0, 7),
breaks=c(seq(0,7,by=1)))
legend <- cowplot::get_legend(
gg1 +
scale_color_gradientn(colors = rampPallFunGD(x=rep(1,nrow(sigcolor))),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 15),
limits=c(0, 7),
breaks=c(seq(0,7,by=1))))## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.## Warning in get_plot_component(plot, "guide-box"): Multiple components found;
## returning the first one. To return all, use `return_all = TRUE`.gg1a <- gg1 + theme(legend.position="none") + ggtitle("") + theme(plot.title = element_text(hjust = 0.5))res=read.csv2(here::here("Output","summary_GDeff_TEW.csv"),
sep=";",dec=".",h=T)
c1=read.csv2(here::here("Output","summary_coeff.csv"),
sep=";",dec=".",h=T)
c1a=subset(c1,model=="allW")
int=c1a$moy[c1a$var=="(Intercept)"]
dsel=subset(mydatatogo,Param=="TE")
v1=seq(round(min(dsel$vel_abs),1),round(max(dsel$vel_abs),1),by=0.1)
v2=seq(round(min(dsel$GD),4),round(max(dsel$GD),4),length.out=nrow(mydatatogo))
for(i in 1:nrow(res)){
pred1=exp(int+res$moy[i]*((v2-mean(mydatatogo$GD))/sd(mydatatogo$GD))+c1a$moy[c1a$var=="LogExtent"]*mean(mydatatogo$LogExtent)+c1a$moy[c1a$var=="LogNtempUnits"]*mean(mydatatogo$LogNtempUnits))
pred1=data.frame(GD=v2,pred1,
vel_abs=res$vel_abs[i],
sig = ifelse(res$pv.inf0[i]<.05 | res$pv.sup0[i]<.05,1,0),
lower = res$q025[i],
upper = res$p975[i])
if(i==1){
pred2=pred1
}else{
pred2=rbind(pred2,pred1)
}
}
res$ID=as.character(round(res$vel_abs,1))
pred2$ID=as.character(round(pred2$vel_abs,1))
resX=subset(res,pv.inf0<0.05)
if(nrow(resX>0)){
selN=round(seq(min(resX$vel_abs),max(resX$vel_abs),le=5),1)
predX=merge(pred2,data.frame(ID=as.character(selN),signifN=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX=pred2
predX$signifN=NA
}
resX=subset(res,pv.sup0<0.05)
if(nrow(resX>0)){
selP=round(seq(min(resX$vel_abs),max(resX$vel_abs),le=5),1)
predX2=merge(predX,data.frame(ID=as.character(selP),signifP=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX2=predX
predX2$signifP=NA
}
predX1=subset(predX2,signifN==1 | signifP==1)
sigcolor <- predX2 %>% group_by(vel_abs) %>% summarise(x = mean(sig))
gg1 <- ggplot(data=predX1, aes(GD, pred1, color = vel_abs, group=ID)) +
GD_graph_params +
scale_color_gradientn(colors = rampPallFunGD(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 8),
limits=c(0, 7),
breaks=c(seq(0,7,by=1))) res=read.csv2(here::here("Output","summary_GDeff_CEW.csv"),
sep=";",dec=".",h=T)
c1=read.csv2(here::here("Output","summary_coeff.csv"),
sep=";",dec=".",h=T)
c1a=subset(c1,model=="allW")
int=c1a$moy[c1a$var=="(Intercept)"]
dsel=subset(mydatatogo,Param=="O")
v1=seq(round(min(dsel$vel_abs),1),round(max(dsel$vel_abs),1),by=0.1)
v2=seq(round(min(dsel$GD),4),round(max(dsel$GD),4),length.out=nrow(mydatatogo))
for(i in 1:nrow(res)){
pred1=exp(int+res$moy[i]*((v2-mean(mydatatogo$GD))/sd(mydatatogo$GD))+c1a$moy[c1a$var=="LogExtent"]*mean(mydatatogo$LogExtent)+c1a$moy[c1a$var=="LogNtempUnits"]*mean(mydatatogo$LogNtempUnits))
pred1=data.frame(GD=v2,pred1,
vel_abs=res$vel_abs[i],
sig = ifelse(res$pv.inf0[i]<.05 | res$pv.sup0[i]<.05,1,0),
lower = res$q025[i],
upper = res$p975[i])
if(i==1){
pred2=pred1
}else{
pred2=rbind(pred2,pred1)
}
}
res$ID=as.character(round(res$vel_abs,1))
pred2$ID=as.character(round(pred2$vel_abs,1))
resX=subset(res,pv.inf0<0.05)
if(nrow(resX>0)){
selN=round(seq(min(resX$vel_abs),max(resX$vel_abs),le=5),1)
predX=merge(pred2,data.frame(ID=as.character(selN),signifN=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX=pred2
predX$signifN=NA
}
resX=subset(res,pv.sup0<0.05)
if(nrow(resX>0)){
selP=round(seq(min(resX$vel_abs),max(resX$vel_abs),le=5),1)
predX2=merge(predX,data.frame(ID=as.character(selP),signifP=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX2=predX
predX2$signifP=NA
}
predX1=subset(predX2,signifN==1 | signifP==1)
sigcolor <- predX2 %>% group_by(vel_abs) %>% summarise(x = mean(sig))
gg2 <- ggplot(data=predX1, aes(GD, pred1, color = vel_abs, group=ID)) +
GD_graph_params +
scale_color_gradientn(colors = rampPallFunGD(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 8),
limits=c(0, 7),
breaks=c(seq(0,7,by=1))) res=read.csv2(here::here("Output","summary_GDeff_LEW.csv"),
sep=";",dec=".",h=T)
c1=read.csv2(here::here("Output","summary_coeff.csv"),
sep=";",dec=".",h=T)
c1a=subset(c1,model=="allW")
int=c1a$moy[c1a$var=="(Intercept)"]
dsel=subset(mydatatogo,Param=="LE")
v1=seq(round(min(dsel$vel_abs),1),round(max(dsel$vel_abs),1),by=0.1)
v2=seq(round(min(dsel$GD),4),round(max(dsel$GD),4),length.out=nrow(mydatatogo))
for(i in 1:nrow(res)){
pred1=exp(int+res$moy[i]*((v2-mean(mydatatogo$GD))/sd(mydatatogo$GD))+c1a$moy[c1a$var=="LogExtent"]*mean(mydatatogo$LogExtent)+c1a$moy[c1a$var=="LogNtempUnits"]*mean(mydatatogo$LogNtempUnits))
pred1=data.frame(GD=v2,pred1,
vel_abs=res$vel_abs[i],
sig = ifelse(res$pv.inf0[i]<.05 | res$pv.sup0[i]<.05,1,0),
lower = res$q025[i],
upper = res$p975[i])
if(i==1){
pred2=pred1
}else{
pred2=rbind(pred2,pred1)
}
}
res$ID=as.character(round(res$vel_abs,1))
pred2$ID=as.character(round(pred2$vel_abs,1))
resX=subset(res,pv.inf0<0.05)
if(nrow(resX>0)){
selN=round(seq(min(resX$vel_abs),max(resX$vel_abs),le=5),1)
predX=merge(pred2,data.frame(ID=as.character(selN),signifN=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX=pred2
predX$signifN=NA
}
resX=subset(res,pv.sup0<0.05)
if(nrow(resX>0)){
selP=round(seq(min(resX$vel_abs),max(resX$vel_abs),le=5),1)
predX2=merge(predX,data.frame(ID=as.character(selP),signifP=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX2=predX
predX2$signifP=NA
}
predX1=subset(predX2,signifN==1 | signifP==1)
sigcolor <- predX2 %>% group_by(vel_abs) %>% summarise(x = mean(sig))
gg3 <- ggplot(data=predX1, aes(GD, pred1, color = vel_abs, group=ID)) +
GD_graph_params +
scale_color_gradientn(colors = rampPallFunGD(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 8),
limits=c(0, 7),
breaks=c(seq(0,7,by=1))) 12.1.1 Fig across edges
(gg1a + legend) +
plot_layout(nrow = 1, widths = c(1,.3))
12.1.2 Group fig v1
p1 <- gg1 + theme(legend.position="none") + ggtitle("Trailing edge") + theme(plot.title = element_text(hjust = 0.5))
p2 <- gg2 + theme(legend.position="none") + ggtitle("Centroid") + theme(plot.title = element_text(hjust = 0.5))
p3 <- gg3 + theme(legend.position="none") + ggtitle("Leading edge") + theme(plot.title = element_text(hjust = 0.5))
legend <- cowplot::get_legend(
gg3 +
scale_color_gradientn(colors = rampPallFunGD(x=rep(1,nrow(sigcolor))),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 15),
limits=c(0, 7),
breaks=c(seq(0,7,by=1))))## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.## Warning in get_plot_component(plot, "guide-box"): Multiple components found;
## returning the first one. To return all, use `return_all = TRUE`.(p1 + p2 + p3 + legend) +
plot_layout(axis_titles = "collect", nrow = 1, widths = c(1,1,1,.3)) + # common axes
plot_annotation(tag_levels = list(c('(a)','(b)','(c)',''))) # figure tags
12.1.3 Group fig v2
ranges <- ggplot_build(p1) %>%
pluck("layout", "panel_params", 1) %>%
`[`(c("x.range", "y.range"))
x <- (0.004 - ranges$x.range[1]) / (ranges$x.range[2] - ranges$x.range[1])
y <- (3 - ranges$y.range[1]) / (ranges$y.range[2] - ranges$y.range[1])
p1 <- gg1 +
theme(legend.title=element_blank(),
plot.title = element_text(hjust = 0.5),
legend.position = c(x, y), legend.justification = c(0, 0)) +
ggtitle("Trailing edge") ## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.p2 <- gg2 +
theme(legend.title=element_blank(),
plot.title = element_text(hjust = 0.5),
legend.position = c(x, y), legend.justification = c(0, 0)) +
ggtitle("Centroid")
p3 <- gg3 +
theme(legend.title=element_blank(),
plot.title = element_text(hjust = 0.5),
legend.position = c(x, y), legend.justification = c(0, 0)) +
ggtitle("Leading edge")
legend <- cowplot::get_legend(
gg3 +
scale_color_gradientn(colors = rampPallFunGD(x=rep(1,nrow(sigcolor))),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 15),
limits=c(0, 7),
breaks=c(seq(0,7,by=1))))## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.## Warning in get_plot_component(plot, "guide-box"): Multiple components found;
## returning the first one. To return all, use `return_all = TRUE`.(p1 + p2 + p3 + legend) +
plot_layout(axis_titles = "collect", nrow = 1, widths = c(1,1,1,.3)) + # common axes
plot_annotation(tag_levels = list(c('(a)','(b)','(c)',''))) # figure tags
12.2 allW2
Model with all edges and looking at the effect of GD at each velocity level
## allW2_TE GDeff
res=read.csv2(here::here("Output","summary_GDeff_allW2_TE.csv"),
sep=";",dec=".",h=T)
c1=read.csv2(here::here("Output","summary_coeff.csv"),
sep=";",dec=".",h=T)
c1a=subset(c1,model=="allW2")
int=c1a$moy[c1a$var=="(Intercept)"]+c1a$moy[c1a$var=="ParamTE"]
dsel=subset(mydatatogo,Param=="TE")
v1=seq(round(min(dsel$vel_abs),1),round(max(dsel$vel_abs),1),by=0.1)
v2=seq(round(min(dsel$GD),4),round(max(dsel$GD),4),length.out=nrow(mydatatogo))
for(i in 1:nrow(res)){
pred1=exp(int+res$moy[i]*((v2-mean(mydatatogo$GD))/sd(mydatatogo$GD))+c1a$moy[c1a$var=="LogExtent"]*mean(mydatatogo$LogExtent)+c1a$moy[c1a$var=="LogNtempUnits"]*mean(mydatatogo$LogNtempUnits))
pred1=data.frame(GD=v2,pred1,
vel_abs=res$vel_abs[i],
sig = ifelse(res$pv.inf0[i]<.05 | res$pv.sup0[i]<.05,1,0),
lower = res$q025[i],
upper = res$p975[i])
if(i==1){
pred2=pred1
}else{
pred2=rbind(pred2,pred1)
}
}
res$ID=as.character(round(res$vel_abs,1))
pred2$ID=as.character(round(pred2$vel_abs,1))
resX=subset(res,pv.inf0<0.05)
if(nrow(resX>0)){
selN=round(seq(min(resX$vel_abs),max(resX$vel_abs),le=5),1)
predX=merge(pred2,data.frame(ID=as.character(selN),signifN=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX=pred2
predX$signifN=NA
}
resX=subset(res,pv.sup0<0.05)
if(nrow(resX>0)){
selP=round(seq(min(resX$vel_abs),max(resX$vel_abs),le=5),1)
predX2=merge(predX,data.frame(ID=as.character(selP),signifP=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX2=predX
predX2$signifP=NA
}
predX1=subset(predX2,signifN==1 | signifP==1)
sigcolor <- predX2 %>% group_by(vel_abs) %>% summarise(x = mean(sig))
gg1 <- ggplot(data=predX1, aes(GD, pred1, color = vel_abs, group=vel_abs)) +
GD_graph_params +
scale_color_gradientn(colors = rampPallFunGD(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 8),
limits=c(0, 7),
breaks=c(seq(0,7,by=1)))
gg12 <- ggplot(data=predX1, aes(GD, pred1, color = vel_abs, group=vel_abs)) +
GD_graph_params +
scale_color_gradientn(colors = rampPallFunGD(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 15),
limits=c(0, 7),
breaks=c(seq(0,7,by=1)))
gg13 <- ggplot(data=predX1, aes(GD, pred1, color = vel_abs, group=vel_abs)) +
stat_density_2d(data = mydatatogo %>% filter(Param == "TE") %>% mutate(pred1 = vel_abs, ID = vel_abs),
aes(GD, SHIFT_abs, fill = ..density.., group = NA), bins = 20, geom = "raster", contour = FALSE) +
GD_graph_params +
scale_color_gradientn(colors = rampPallFunGD(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 8),
limits=c(0, 7),
breaks=c(seq(0,7,by=1)))+
scale_fill_viridis_c(name = "N obs",option = "B", direction = -1, limits = c(.001, 65), na.value = "white")res=read.csv2(here::here("Output","summary_GDeff_allW2_CE.csv"),
sep=";",dec=".",h=T)
c1=read.csv2(here::here("Output","summary_coeff.csv"),
sep=";",dec=".",h=T)
c1a=subset(c1,model=="allW2")
int=c1a$moy[c1a$var=="(Intercept)"]
dsel=subset(mydatatogo,Param=="O")
v1=seq(round(min(dsel$vel_abs),1),round(max(dsel$vel_abs),1),by=0.1)
v2=seq(round(min(dsel$GD),4),round(max(dsel$GD),4),length.out=nrow(mydatatogo))
for(i in 1:nrow(res)){
pred1=exp(int+res$moy[i]*((v2-mean(mydatatogo$GD))/sd(mydatatogo$GD))+c1a$moy[c1a$var=="LogExtent"]*mean(mydatatogo$LogExtent)+c1a$moy[c1a$var=="LogNtempUnits"]*mean(mydatatogo$LogNtempUnits))
pred1=data.frame(GD=v2,pred1,
vel_abs=res$vel_abs[i],
sig = ifelse(res$pv.inf0[i]<.05 | res$pv.sup0[i]<.05,1,0),
lower = res$q025[i],
upper = res$p975[i])
if(i==1){
pred2=pred1
}else{
pred2=rbind(pred2,pred1)
}
}
res$ID=as.character(round(res$vel_abs,1))
pred2$ID=as.character(round(pred2$vel_abs,1))
resX=subset(res,pv.inf0<0.05)
if(nrow(resX>0)){
selN=round(seq(min(resX$vel_abs),max(resX$vel_abs),le=5),1)
predX=merge(pred2,data.frame(ID=as.character(selN),signifN=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX=pred2
predX$signifN=NA
}
resX=subset(res,pv.sup0<0.05)
if(nrow(resX>0)){
selP=round(seq(min(resX$vel_abs),max(resX$vel_abs),le=5),1)
predX2=merge(predX,data.frame(ID=as.character(selP),signifP=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX2=predX
predX2$signifP=NA
}
predX1=subset(predX2,signifN==1 | signifP==1)
sigcolor <- predX2 %>% group_by(vel_abs) %>% summarise(x = mean(sig))
gg2 <- ggplot(data=predX1, aes(GD, pred1, color = vel_abs, group=vel_abs)) +
GD_graph_params +
scale_color_gradientn(colors = rampPallFunGD(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 8),
limits=c(0, 7),
breaks=c(seq(0,7,by=1)))
gg22 <- ggplot(data=predX1, aes(GD, pred1, color = vel_abs, group=vel_abs)) +
GD_graph_params +
scale_color_gradientn(colors = rampPallFunGD(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 15),
limits=c(0, 7),
breaks=c(seq(0,7,by=1)))
gg23 <- ggplot(data=predX1, aes(GD, pred1, color = vel_abs, group=vel_abs)) +
stat_density_2d(data = mydatatogo %>% filter(Param == "O") %>% mutate(pred1 = vel_abs, ID = vel_abs),
aes(GD, SHIFT_abs, fill = ..density.., group = NA), bins = 20, geom = "raster", contour = FALSE) +
GD_graph_params +
scale_color_gradientn(colors = rampPallFunGD(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 8),
limits=c(0, 7),
breaks=c(seq(0,7,by=1)))+
scale_fill_viridis_c(name = "N obs",option = "B", direction = -1, limits = c(.001, 65), na.value = "white")## allW2_LE GDeff
res=read.csv2(here::here("Output","summary_GDeff_allW2_LE.csv"),
sep=";",dec=".",h=T)
c1=read.csv2(here::here("Output","summary_coeff.csv"),
sep=";",dec=".",h=T)
c1a=subset(c1,model=="allW2")
int=c1a$moy[c1a$var=="(Intercept)"]+c1a$moy[c1a$var=="ParamLE"]
dsel=subset(mydatatogo,Param=="LE")
v1=seq(round(min(dsel$vel_abs),1),round(max(dsel$vel_abs),1),by=0.1)
v2=seq(round(min(dsel$GD),4),round(max(dsel$GD),4),length.out=nrow(mydatatogo))
for(i in 1:nrow(res)){
pred1=exp(int+res$moy[i]*((v2-mean(mydatatogo$GD))/sd(mydatatogo$GD))+c1a$moy[c1a$var=="LogExtent"]*mean(mydatatogo$LogExtent)+c1a$moy[c1a$var=="LogNtempUnits"]*mean(mydatatogo$LogNtempUnits))
pred1=data.frame(GD=v2,pred1,
vel_abs=res$vel_abs[i],
sig = ifelse(res$pv.inf0[i]<.05 | res$pv.sup0[i]<.05,1,0),
lower = res$q025[i],
upper = res$p975[i])
if(i==1){
pred2=pred1
}else{
pred2=rbind(pred2,pred1)
}
}
res$ID=as.character(round(res$vel_abs,1))
pred2$ID=as.character(round(pred2$vel_abs,1))
resX=subset(res,pv.inf0<0.05)
if(nrow(resX>0)){
selN=round(seq(min(resX$vel_abs),max(resX$vel_abs),le=5),1)
predX=merge(pred2,data.frame(ID=as.character(selN),signifN=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX=pred2
predX$signifN=NA
}
resX=subset(res,pv.sup0<0.05)
if(nrow(resX>0)){
selP=round(seq(min(resX$vel_abs),max(resX$vel_abs),le=5),1)
predX2=merge(predX,data.frame(ID=as.character(selP),signifP=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX2=predX
predX2$signifP=NA
}
predX1=subset(predX2,signifN==1 | signifP==1)
sigcolor <- predX2 %>% group_by(vel_abs) %>% summarise(x = mean(sig))
gg3 <- ggplot(data=predX1, aes(GD, pred1, color = vel_abs, group=ID)) +
GD_graph_params +
scale_color_gradientn(colors = rampPallFunGD(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 8),
limits=c(0, 7),
breaks=c(seq(0,7,by=1)))
gg32 <- ggplot(data=predX1, aes(GD, pred1, color = vel_abs, group=vel_abs)) +
GD_graph_params +
scale_color_gradientn(colors = rampPallFunGD(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 15),
limits=c(0, 7),
breaks=c(seq(0,7,by=1)))
gg33 <- ggplot(data=predX1, aes(GD, pred1, color = vel_abs, group=vel_abs)) +
stat_density_2d(data = mydatatogo %>% filter(Param == "LE") %>% mutate(pred1 = vel_abs, ID = vel_abs),
aes(GD, SHIFT_abs, fill = ..density.., group = NA), bins = 20, geom = "raster", contour = FALSE) +
GD_graph_params +
scale_color_gradientn(colors = rampPallFunGD(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 8),
limits=c(0, 7),
breaks=c(seq(0,7,by=1)))+
scale_fill_viridis_c(name = "N obs",option = "B", direction = -1, limits = c(.001, 65), na.value = "white")12.2.1 Group fig v1
p1 <- gg1 + theme(legend.position="none") + ggtitle("Trailing edge") + theme(plot.title = element_text(hjust = 0.5))
p2 <- gg2 + theme(legend.position="none") + ggtitle("Centroid") + theme(plot.title = element_text(hjust = 0.5))
p3 <- gg3 + theme(legend.position="none") + ggtitle("Leading edge") + theme(plot.title = element_text(hjust = 0.5))
legend <- cowplot::get_legend(
gg3 +
scale_color_gradientn(colors = rampPallFunGD(x=rep(1,nrow(sigcolor))),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 15),
limits=c(0, 7),
breaks=c(seq(0,7,by=1))))## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.## Warning in get_plot_component(plot, "guide-box"): Multiple components found;
## returning the first one. To return all, use `return_all = TRUE`.(p1 + p2 + p3 + legend) +
plot_layout(axis_titles = "collect", nrow = 1, widths = c(1,1,1,.3)) + # common axes
plot_annotation(tag_levels = list(c('(a)','(b)','(c)',''))) # figure tags
p1 <- gg1 + theme(legend.position="none") + ggtitle("Trailing edge") + theme(plot.title = element_text(hjust = 0.5)) p1 <- gg13 + theme(legend.position="none") + ggtitle("Trailing edge") + theme(plot.title = element_text(hjust = 0.5))
p2 <- gg23 + theme(legend.position="none") + ggtitle("Centroid") + theme(plot.title = element_text(hjust = 0.5))
p3 <- gg33 + theme(legend.position="none") + ggtitle("Leading edge") + theme(plot.title = element_text(hjust = 0.5))
legend <- cowplot::get_legend(
gg33 +
scale_color_gradientn(colors = rampPallFunGD(x=rep(1,nrow(sigcolor))),
name = bquote(Absolute~velocity~of~range~shift~(km.yr^1)),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 10),
limits=c(0, 7),
breaks=c(seq(0,7,by=1)))+
scale_fill_viridis_c(name = "Number of observations", option = "A",
direction = -1, limits = c(.001, 65), na.value = "white",
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 8)))## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.(p1 + p2 + p3 + legend) +
plot_layout(axis_titles = "collect", nrow = 1, widths = c(1,1,1,.3)) + # common axes
plot_annotation(tag_levels = list(c('(a)','(b)','(c)',''))) # figure tags
12.2.2 Group fig v2
ranges <- ggplot_build(p1) %>%
pluck("layout", "panel_params", 1) %>%
`[`(c("x.range", "y.range"))## Warning: Removed 35 rows containing non-finite outside the scale range
## (`stat_density2d()`).## Warning: The following aesthetics were dropped during statistical transformation:
## colour.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?x <- (0.004 - ranges$x.range[1]) / (ranges$x.range[2] - ranges$x.range[1])
y <- (3 - ranges$y.range[1]) / (ranges$y.range[2] - ranges$y.range[1])
p1 <- gg12 +
ggtitle("Trailing edge")
p2 <- gg22 +
ggtitle("Centroid")
p3 <- gg32 +
ggtitle("Leading edge")
(p1 + p2 + p3) +
plot_layout(axis_titles = "collect", nrow = 1, widths = c(1,1,1)) + # common axes
plot_annotation(tag_levels = list(c('(a)','(b)','(c)'))) # figure tags
12.2.3 Group fig v3
ranges <- ggplot_build(p1) %>%
pluck("layout", "panel_params", 1) %>%
`[`(c("x.range", "y.range"))
x <- (0.004 - ranges$x.range[1]) / (ranges$x.range[2] - ranges$x.range[1])
y <- (3.3 - ranges$y.range[1]) / (ranges$y.range[2] - ranges$y.range[1])
p1 <- gg1 +
theme(legend.title=element_blank(),
plot.title = element_text(hjust = 0.5),
legend.position = c(x, y), legend.justification = c(0, 0)) +
ggtitle("Trailing edge")
p2 <- gg2 +
theme(legend.title=element_blank(),
plot.title = element_text(hjust = 0.5),
legend.position = c(x, y), legend.justification = c(0, 0)) +
ggtitle("Centroid")
p3 <- gg3 +
theme(legend.title=element_blank(),
plot.title = element_text(hjust = 0.5),
legend.position = c(x, y), legend.justification = c(0, 0)) +
ggtitle("Leading edge")
legend <- cowplot::get_legend(
gg3 +
scale_color_gradientn(colors = rampPallFunGD(x=rep(1,nrow(sigcolor))),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 15),
limits=c(0, 7),
breaks=c(seq(0,7,by=1))))## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.## Warning in get_plot_component(plot, "guide-box"): Multiple components found;
## returning the first one. To return all, use `return_all = TRUE`.(p1 + p2 + p3 + legend) +
plot_layout(axis_titles = "collect", nrow = 1, widths = c(1,1,1,.3)) + # common axes
plot_annotation(tag_levels = list(c('(a)','(b)','(c)',''))) # figure tags
In general, GD helps species to stay in place as climate velocities increase across all positions of a species range. At the trailing edge, the effect of GD of range shifts is significantly negative when climate velocity is high (>3km.yr^1), indicating that species stay in place (or lag more from climate change) when exposure is high.
As climate velocity increases, low GD species move faster than high GD species to keep up with climate change. This corroborates the hypotheses that GD increases species’ ability to stay in place. The fact that low GD species show a positive relationship with climate velocity indicates that climate exposure poses higher pressure on them relative to high GD species. As a consequence, low GD species have to move to track changing climates and stay under their thermal niches, otherwise they would go extinct. As GD increases, the effect of climate exposure on species range shift is relaxed – species track climate change less closely.
Across all range edges, low GD species move faster as climate velocity increases, whereas high GD species move slower.
12.2.4 Significant velocity values
res=lapply(c("LE","CE","TE"), function(x){
tmp <- read.csv2(here::here("Output",paste0("summary_GDeff_allW2_","TE",".csv")),
sep=";",dec=".",h=T)
tmp$edge = x
return(tmp)
})
res <- rbindlist(res)
res %>%
group_by(edge) %>%
summarise(sig_inf = paste(vel_abs), collapse = " ,")## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## `summarise()` has grouped output by 'edge'. You can override using the
## `.groups` argument.## # A tibble: 210 × 3
## # Groups: edge [3]
## edge sig_inf collapse
## <chr> <chr> <chr>
## 1 CE 0 " ,"
## 2 CE 0.1 " ,"
## 3 CE 0.2 " ,"
## 4 CE 0.3 " ,"
## 5 CE 0.4 " ,"
## 6 CE 0.5 " ,"
## 7 CE 0.6 " ,"
## 8 CE 0.7 " ,"
## 9 CE 0.8 " ,"
## 10 CE 0.9 " ,"
## # ℹ 200 more rows13 Graphics effect velocity
# GD_graph_params <-
# list(geom_line(),
# labs(x = "Genetic diversity",
# y = bquote('Absolute velocity of range shift '(km.yr^1)),
# color = bquote('Absolute velocity of isotherm shifts '(km.yr^1))),
# scale_color_gradientn(colors = c("#00AFBB", "#E7B800", "#FC4E07"),
# guide = guide_colourbar(title.position = "right",
# barwidth = 1, barheight = 10),
# limits=c(0, 7),
# breaks=c(seq(0,7,by=1))),
# theme_classic(),
# scale_x_continuous(breaks=c(c(seq(0,0.05,by=.01))),
# labels=c("0",as.character(seq(0.01,0.05,by=.01))),
# limits=c(0,0.055), expand=c(0,0)),
# scale_y_continuous(breaks=c(seq(0,6,by=2)), limits=c(0,6), expand=c(0,0)),
# theme(legend.title = element_text(angle = -90),
# legend.title.align = 0.5))
rampPallFunVel <- function(x){
tmp <- viridis_pal(option = "D")(length(x))
pos <- which(x==0)
tmp[pos] <- "lightgray"
return(tmp)
}
vel_graph_params <-
list(geom_line(size = 1),
labs(x = bquote('Absolute velocity of isotherm shifts '(km.yr^1)),
y = bquote('Absolute velocity of range shift '(km.yr^1)),
color = "Genetic diversity"),
theme_classic(),
scale_x_continuous(n.breaks = 8,
limits=c(0,7),
expand=c(0,0)),
scale_y_continuous(n.breaks = 6,
limits=c(0,10),
expand=c(0,0)),
theme(legend.title = element_text(angle = -90),
legend.title.align = 0.5))13.1 AllW
## AllW
res=read.csv2(here::here("Output","summary_VAeff_allW.csv"),
sep=";",dec=".",h=T)
c1=read.csv2(here::here("Output","summary_coeff.csv"),
sep=";",dec=".",h=T)
c1a=subset(c1,model=="allW")
int=c1a$moy[c1a$var=="(Intercept)"]
dsel=mydatatogo
v2=seq(round(min(dsel$vel_abs),1),round(max(dsel$vel_abs),1),by=0.1)
for(i in 1:nrow(res)){
pred1=exp(int+res$moy[i]*((v2-mean(mydatatogo$vel_abs))/sd(mydatatogo$vel_abs))+c1a$moy[c1a$var=="LogExtent"]*mean(mydatatogo$LogExtent)+c1a$moy[c1a$var=="LogNtempUnits"]*mean(mydatatogo$LogNtempUnits))
pred1=data.frame(GD=res$GD[i],
pred1,
vel_abs=v2,
sig = ifelse(res$pv.inf0[i]<.05 | res$pv.sup0[i]<.05,1,0),
lower = res$q025[i],
upper = res$p975[i])
if(i==1){
pred2=pred1
}else{
pred2=rbind(pred2,pred1)
}
}
res$ID=as.character(round(res$GD,3))
pred2$ID=as.character(round(pred2$GD,3))
resX=subset(res,pv.inf0<0.05)
if(nrow(resX>0)){
if(nrow(resX)<5){
selN=round(seq(min(resX$GD),max(resX$GD),le=nrow(resX)),3)
}else{
selN=round(seq(min(resX$GD),max(resX$GD),le=5),3)
}
predX=merge(pred2,data.frame(ID=as.character(selN),signifN=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX=pred2
predX$signifN=NA
}
resX=subset(res,pv.sup0<0.05)
if(nrow(resX>0)){
if(nrow(resX)<5){
selP=round(seq(min(resX$GD),max(resX$GD),le=nrow(resX)),3)
}else{
selP=round(seq(min(resX$GD),max(resX$GD),le=5),3)
}
predX2=merge(predX,data.frame(ID=as.character(selP),signifP=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX2=predX
predX2$signifP=NA
}
predX1=subset(predX2,signifN==1 | signifP==1)
sigcolor <- predX2 %>% group_by(vel_abs) %>% summarise(x = mean(sig))
gg1 <- ggplot(data=predX1, aes(vel_abs, pred1, color = GD, group=ID)) +
vel_graph_params +
scale_color_gradientn(colors = rampPallFunVel(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 8),
limits=c(0, 0.05),
breaks=c(seq(0,0.05,by=.01)))
legend <- cowplot::get_legend(
gg1 +
scale_color_gradientn(colors = rampPallFunVel(x=rep(1,nrow(sigcolor))),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 15),
limits=c(0, 0.05),
breaks=c(seq(0,0.05,by=.01))))## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.## Warning in get_plot_component(plot, "guide-box"): Multiple components found;
## returning the first one. To return all, use `return_all = TRUE`.gg1a <- gg1 + theme(legend.position="none") + ggtitle("") + theme(plot.title = element_text(hjust = 0.5))
### TEW
res=read.csv2(here::here("Output","summary_VAeff_TEW.csv"),
sep=";",dec=".",h=T)
c1=read.csv2(here::here("Output","summary_coeff.csv"),
sep=";",dec=".",h=T)
c1a=subset(c1,model=="allW")
int=c1a$moy[c1a$var=="(Intercept)"]
dsel=subset(mydatatogo,Param=="TE")
v2=seq(round(min(dsel$vel_abs),1),round(max(dsel$vel_abs),1),by=0.1)
for(i in 1:nrow(res)){
pred1=exp(int+res$moy[i]*((v2-mean(mydatatogo$vel_abs))/sd(mydatatogo$vel_abs))+c1a$moy[c1a$var=="LogExtent"]*mean(mydatatogo$LogExtent)+c1a$moy[c1a$var=="LogNtempUnits"]*mean(mydatatogo$LogNtempUnits))
pred1=data.frame(GD=res$GD[i],
pred1,
vel_abs=v2,
sig = ifelse(res$pv.inf0[i]<.05 | res$pv.sup0[i]<.05,1,0),
lower = res$q025[i],
upper = res$p975[i])
if(i==1){
pred2=pred1
}else{
pred2=rbind(pred2,pred1)
}
}
res$ID=as.character(round(res$GD,3))
pred2$ID=as.character(round(pred2$GD,3))
resX=subset(res,pv.inf0<0.05)
if(nrow(resX>0)){
if(nrow(resX)<5){
selN=round(seq(min(resX$GD),max(resX$GD),le=nrow(resX)),3)
}else{
selN=round(seq(min(resX$GD),max(resX$GD),le=5),3)
}
predX=merge(pred2,data.frame(ID=as.character(selN),signifN=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX=pred2
predX$signifN=NA
}
resX=subset(res,pv.sup0<0.05)
if(nrow(resX>0)){
if(nrow(resX)<5){
selP=round(seq(min(resX$GD),max(resX$GD),le=nrow(resX)),3)
}else{
selP=round(seq(min(resX$GD),max(resX$GD),le=5),3)
}
predX2=merge(predX,data.frame(ID=as.character(selP),signifP=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX2=predX
predX2$signifP=NA
}
predX1=subset(predX2,signifN==1 | signifP==1)
sigcolor <- predX2 %>% group_by(vel_abs) %>% summarise(x = mean(sig))
gg1 <- ggplot(data=predX1, aes(vel_abs, pred1, color = GD, group=ID)) +
vel_graph_params +
scale_color_gradientn(colors = rampPallFunVel(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 8),
limits=c(0, 0.05),
breaks=c(seq(0,0.05,by=.01))) ### CEW
res=read.csv2(here::here("Output","summary_VAeff_CEW.csv"),
sep=";",dec=".",h=T)
c1=read.csv2(here::here("Output","summary_coeff.csv"),
sep=";",dec=".",h=T)
c1a=subset(c1,model=="allW")
int=c1a$moy[c1a$var=="(Intercept)"]
dsel=subset(mydatatogo,Param=="O")
v2=seq(round(min(dsel$vel_abs),1),round(max(dsel$vel_abs),1),by=0.1)
for(i in 1:nrow(res)){
pred1=exp(int+res$moy[i]*((v2-mean(mydatatogo$vel_abs))/sd(mydatatogo$vel_abs))+c1a$moy[c1a$var=="LogExtent"]*mean(mydatatogo$LogExtent)+c1a$moy[c1a$var=="LogNtempUnits"]*mean(mydatatogo$LogNtempUnits))
pred1=data.frame(GD=res$GD[i],
pred1,
vel_abs=v2,
sig = ifelse(res$pv.inf0[i]<.05 | res$pv.sup0[i]<.05,1,0),
lower = res$q025[i],
upper = res$p975[i])
if(i==1){
pred2=pred1
}else{
pred2=rbind(pred2,pred1)
}
}
res$ID=as.character(round(res$GD,3))
pred2$ID=as.character(round(pred2$GD,3))
resX=subset(res,pv.inf0<0.05)
if(nrow(resX>0)){
if(nrow(resX)<5){
selN=round(seq(min(resX$GD),max(resX$GD),le=nrow(resX)),3)
}else{
selN=round(seq(min(resX$GD),max(resX$GD),le=5),3)
}
predX=merge(pred2,data.frame(ID=as.character(selN),signifN=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX=pred2
predX$signifN=NA
}
resX=subset(res,pv.sup0<0.05)
if(nrow(resX>0)){
if(nrow(resX)<5){
selP=round(seq(min(resX$GD),max(resX$GD),le=nrow(resX)),3)
}else{
selP=round(seq(min(resX$GD),max(resX$GD),le=5),3)
}
predX2=merge(predX,data.frame(ID=as.character(selP),signifP=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX2=predX
predX2$signifP=NA
}
predX1=subset(predX2,signifN==1 | signifP==1)
sigcolor <- predX2 %>% group_by(vel_abs) %>% summarise(x = mean(sig))
gg2 <- ggplot(data=predX1, aes(vel_abs, pred1, color = GD, group=ID)) +
vel_graph_params +
scale_color_gradientn(colors = rampPallFunVel(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 8),
limits=c(0, 0.05),
breaks=c(seq(0,0.05,by=.01))) ### LEW
res=read.csv2(here::here("Output","summary_VAeff_LEW.csv"),
sep=";",dec=".",h=T)
c1=read.csv2(here::here("Output","summary_coeff.csv"),
sep=";",dec=".",h=T)
c1a=subset(c1,model=="allW")
int=c1a$moy[c1a$var=="(Intercept)"]
dsel=subset(mydatatogo,Param=="LE")
v2=seq(round(min(dsel$vel_abs),1),round(max(dsel$vel_abs),1),by=0.1)
for(i in 1:nrow(res)){
pred1=exp(int+res$moy[i]*((v2-mean(mydatatogo$vel_abs))/sd(mydatatogo$vel_abs))+c1a$moy[c1a$var=="LogExtent"]*mean(mydatatogo$LogExtent)+c1a$moy[c1a$var=="LogNtempUnits"]*mean(mydatatogo$LogNtempUnits))
pred1=data.frame(GD=res$GD[i],
pred1,
vel_abs=v2,
sig = ifelse(res$pv.inf0[i]<.05 | res$pv.sup0[i]<.05,1,0),
lower = res$q025[i],
upper = res$p975[i])
if(i==1){
pred2=pred1
}else{
pred2=rbind(pred2,pred1)
}
}
res$ID=as.character(round(res$GD,3))
pred2$ID=as.character(round(pred2$GD,3))
resX=subset(res,pv.inf0<0.05)
if(nrow(resX>0)){
if(nrow(resX)<5){
selN=round(seq(min(resX$GD),max(resX$GD),le=nrow(resX)),3)
}else{
selN=round(seq(min(resX$GD),max(resX$GD),le=5),3)
}
predX=merge(pred2,data.frame(ID=as.character(selN),signifN=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX=pred2
predX$signifN=NA
}
resX=subset(res,pv.sup0<0.05)
if(nrow(resX>0)){
if(nrow(resX)<5){
selP=round(seq(min(resX$GD),max(resX$GD),le=nrow(resX)),3)
}else{
selP=round(seq(min(resX$GD),max(resX$GD),le=5),3)
}
predX2=merge(predX,data.frame(ID=as.character(selP),signifP=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX2=predX
predX2$signifP=NA
}
predX1=subset(predX2,signifN==1 | signifP==1)
sigcolor <- predX2 %>% group_by(vel_abs) %>% summarise(x = mean(sig))
gg3 <- ggplot(data=predX1, aes(vel_abs, pred1, color = GD, group=ID)) +
vel_graph_params +
scale_color_gradientn(colors = rampPallFunVel(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 8),
limits=c(0, 0.05),
breaks=c(seq(0,0.05,by=.01))) 13.1.1 Fig across edges
(gg1a + legend) +
plot_layout(nrow = 1, widths = c(1,.3))
13.1.2 Group fig v1
p1 <- gg1 + theme(legend.position="none") + ggtitle("Trailing edge") + theme(plot.title = element_text(hjust = 0.5))
p2 <- gg2 + theme(legend.position="none") + ggtitle("Centroid") + theme(plot.title = element_text(hjust = 0.5))
p3 <- gg3 + theme(legend.position="none") + ggtitle("Leading edge") + theme(plot.title = element_text(hjust = 0.5))
legend <- cowplot::get_legend(
gg3 +
scale_color_gradientn(colors = rampPallFunVel(x=rep(1,nrow(sigcolor))),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 15),
limits=c(0, 7),
breaks=c(seq(0,7,by=1))))## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.## Warning in get_plot_component(plot, "guide-box"): Multiple components found;
## returning the first one. To return all, use `return_all = TRUE`.(p1 + p2 + p3 + legend) +
plot_layout(axis_titles = "collect", nrow = 1, widths = c(1,1,1,.3)) + # common axes
plot_annotation(tag_levels = list(c('(a)','(b)','(c)',''))) # figure tags
13.1.3 Group fig v2
ranges <- ggplot_build(p1) %>%
pluck("layout", "panel_params", 1) %>%
`[`(c("x.range", "y.range"))
x <- (.5 - ranges$x.range[1]) / (ranges$x.range[2] - ranges$x.range[1])
y <- (5 - ranges$y.range[1]) / (ranges$y.range[2] - ranges$y.range[1])
p1 <- gg1 +
theme(legend.title=element_blank(),
plot.title = element_text(hjust = 0.5),
legend.position = c(x, y), legend.justification = c(0, 0)) +
ggtitle("Trailing edge")
p2 <- gg2 +
theme(legend.title=element_blank(),
plot.title = element_text(hjust = 0.5),
legend.position = c(x, y), legend.justification = c(0, 0)) +
ggtitle("Centroid")
p3 <- gg3 +
theme(legend.title=element_blank(),
plot.title = element_text(hjust = 0.5),
legend.position = c(x, y), legend.justification = c(0, 0)) +
ggtitle("Leading edge")
p4 <- gg3 +
theme(legend.title=element_blank(),
plot.title = element_text(hjust = 0.5),
legend.position = c(x, y), legend.justification = c(0, 0)) +
ggtitle("Leading edge")
legend <- cowplot::get_legend(
gg3 +
scale_color_gradientn(colors = rampPallFunVel(x=rep(1,nrow(sigcolor))),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 15),
limits=c(0, 7),
breaks=c(seq(0,7,by=1))))## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.## Warning in get_plot_component(plot, "guide-box"): Multiple components found;
## returning the first one. To return all, use `return_all = TRUE`.(p1 + p2 + p3 + legend) +
plot_layout(axis_titles = "collect", nrow = 1, widths = c(1,1,1,.3)) + # common axes
plot_annotation(tag_levels = list(c('(a)','(b)','(c)',''))) # figure tags
13.2 allW2
Model with all edges and looking at the effect of GD at each velocity level
## allW2_TE VAeff
res=read.csv2(here::here("Output","summary_VAeff_allW2_TE.csv"),
sep=";",dec=".",h=T)
c1=read.csv2(here::here("Output","summary_coeff.csv"),
sep=";",dec=".",h=T)
c1a=subset(c1,model=="allW2")
int=c1a$moy[c1a$var=="(Intercept)"]+c1a$moy[c1a$var=="ParamTE"]
dsel=subset(mydatatogo,Param=="TE")
v2=seq(round(min(dsel$vel_abs),1),round(max(dsel$vel_abs),1),by=0.1)
for(i in 1:nrow(res)){
pred1=exp(int+res$moy[i]*((v2-mean(mydatatogo$vel_abs))/sd(mydatatogo$vel_abs))+c1a$moy[c1a$var=="LogExtent"]*mean(mydatatogo$LogExtent)+c1a$moy[c1a$var=="LogNtempUnits"]*mean(mydatatogo$LogNtempUnits))
pred1=data.frame(GD=res$GD[i],
pred1,
vel_abs=v2,
sig = ifelse(res$pv.inf0[i]<.05 | res$pv.sup0[i]<.05,1,0),
lower = res$q025[i],
upper = res$p975[i])
if(i==1){
pred2=pred1
}else{
pred2=rbind(pred2,pred1)
}
}
res$ID=as.character(round(res$GD,3))
pred2$ID=as.character(round(pred2$GD,3))
resX=subset(res,pv.inf0<0.05)
if(nrow(resX>0)){
if(nrow(resX)<5){
selN=round(seq(min(resX$GD),max(resX$GD),le=nrow(resX)),3)
}else{
selN=round(seq(min(resX$GD),max(resX$GD),le=5),3)
}
predX=merge(pred2,data.frame(ID=as.character(selN),signifN=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX=pred2
predX$signifN=NA
}
resX=subset(res,pv.sup0<0.05)
if(nrow(resX>0)){
if(nrow(resX)<5){
selP=round(seq(min(resX$GD),max(resX$GD),le=nrow(resX)),3)
}else{
selP=round(seq(min(resX$GD),max(resX$GD),le=5),3)
}
predX2=merge(predX,data.frame(ID=as.character(selP),signifP=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX2=predX
predX2$signifP=NA
}
predX1=subset(predX2,signifN==1 | signifP==1)
sigcolor <- predX2 %>% group_by(GD) %>% summarise(x = mean(sig))
gg1 <- ggplot(data=predX1, aes(vel_abs, pred1, color = GD, group=ID)) +
vel_graph_params +
scale_color_gradientn(colors = rampPallFunVel(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 8),
limits=c(0, 0.055),
breaks=c(seq(0,0.055,by=.01)))
gg12 <- ggplot(data=predX1, aes(vel_abs, pred1, color = GD, group=ID)) +
vel_graph_params +
scale_color_gradientn(colors = rampPallFunVel(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 15),
limits=c(0, 0.055),
breaks=c(seq(0,0.055,by=.01))) ## allW2_CE VAeff
res=read.csv2(here::here("Output","summary_VAeff_allW2_CE.csv"),
sep=";",dec=".",h=T)
c1=read.csv2(here::here("Output","summary_coeff.csv"),
sep=";",dec=".",h=T)
c1a=subset(c1,model=="allW2")
int=c1a$moy[c1a$var=="(Intercept)"]
dsel=subset(mydatatogo,Param=="O")
v2=seq(round(min(dsel$vel_abs),1),round(max(dsel$vel_abs),1),by=0.1)
for(i in 1:nrow(res)){
pred1=exp(int+res$moy[i]*((v2-mean(mydatatogo$vel_abs))/sd(mydatatogo$vel_abs))+c1a$moy[c1a$var=="LogExtent"]*mean(mydatatogo$LogExtent)+c1a$moy[c1a$var=="LogNtempUnits"]*mean(mydatatogo$LogNtempUnits))
pred1=data.frame(GD=res$GD[i],
pred1,
vel_abs=v2,
sig = ifelse(res$pv.inf0[i]<.05 | res$pv.sup0[i]<.05,1,0),
lower = res$q025[i],
upper = res$p975[i])
if(i==1){
pred2=pred1
}else{
pred2=rbind(pred2,pred1)
}
}
res$ID=as.character(round(res$GD,3))
pred2$ID=as.character(round(pred2$GD,3))
resX=subset(res,pv.inf0<0.05)
if(nrow(resX>0)){
if(nrow(resX)<5){
selN=round(seq(min(resX$GD),max(resX$GD),le=nrow(resX)),3)
}else{
selN=round(seq(min(resX$GD),max(resX$GD),le=5),3)
}
predX=merge(pred2,data.frame(ID=as.character(selN),signifN=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX=pred2
predX$signifN=NA
}
resX=subset(res,pv.sup0<0.05)
if(nrow(resX>0)){
if(nrow(resX)<5){
selP=round(seq(min(resX$GD),max(resX$GD),le=nrow(resX)),3)
}else{
selP=round(seq(min(resX$GD),max(resX$GD),le=5),3)
}
predX2=merge(predX,data.frame(ID=as.character(selP),signifP=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX2=predX
predX2$signifP=NA
}
predX1=subset(predX2,signifN==1 | signifP==1)
sigcolor <- predX2 %>% group_by(GD) %>% summarise(x = mean(sig))
gg2 <- ggplot(data=predX1, aes(vel_abs, pred1, color = GD, group=ID)) +
vel_graph_params +
scale_color_gradientn(colors = rampPallFunVel(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 8),
limits=c(0, 0.055),
breaks=c(seq(0,0.055,by=.01)))
gg22 <- ggplot(data=predX1, aes(vel_abs, pred1, color = GD, group=ID)) +
vel_graph_params +
scale_color_gradientn(colors = rampPallFunVel(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 15),
limits=c(0, 0.055),
breaks=c(seq(0,0.055,by=.01))) ## allW2_LE VAeff
res=read.csv2(here::here("Output","summary_VAeff_allW2_LE.csv"),
sep=";",dec=".",h=T)
c1=read.csv2(here::here("Output","summary_coeff.csv"),
sep=";",dec=".",h=T)
c1a=subset(c1,model=="allW2")
int=c1a$moy[c1a$var=="(Intercept)"]+c1a$moy[c1a$var=="ParamLE"]
dsel=subset(mydatatogo,Param=="LE")
v2=seq(round(min(dsel$vel_abs),1),round(max(dsel$vel_abs),1),by=0.1)
for(i in 1:nrow(res)){
pred1=exp(int+res$moy[i]*((v2-mean(mydatatogo$vel_abs))/sd(mydatatogo$vel_abs))+c1a$moy[c1a$var=="LogExtent"]*mean(mydatatogo$LogExtent)+c1a$moy[c1a$var=="LogNtempUnits"]*mean(mydatatogo$LogNtempUnits))
pred1=data.frame(GD=res$GD[i],
pred1,
vel_abs=v2,
sig = ifelse(res$pv.inf0[i]<.05 | res$pv.sup0[i]<.05,1,0),
lower = res$q025[i],
upper = res$p975[i])
if(i==1){
pred2=pred1
}else{
pred2=rbind(pred2,pred1)
}
}
res$ID=as.character(round(res$GD,3))
pred2$ID=as.character(round(pred2$GD,3))
resX=subset(res,pv.inf0<0.05)
if(nrow(resX>0)){
if(nrow(resX)<5){
selN=round(seq(min(resX$GD),max(resX$GD),le=nrow(resX)),3)
}else{
selN=round(seq(min(resX$GD),max(resX$GD),le=5),3)
}
predX=merge(pred2,data.frame(ID=as.character(selN),signifN=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX=pred2
predX$signifN=NA
}
resX=subset(res,pv.sup0<0.05)
if(nrow(resX>0)){
if(nrow(resX)<5){
selP=round(seq(min(resX$GD),max(resX$GD),le=nrow(resX)),3)
}else{
selP=round(seq(min(resX$GD),max(resX$GD),le=5),3)
}
predX2=merge(predX,data.frame(ID=as.character(selP),signifP=1),by.x="ID",by.y="ID",all.x=T)
}else{
predX2=predX
predX2$signifP=NA
}
predX1=subset(predX2,signifN==1 | signifP==1)
sigcolor <- predX2 %>% group_by(GD) %>% summarise(x = mean(sig))
gg3 <- ggplot(data=predX1, aes(vel_abs, pred1, color = GD, group=ID)) +
vel_graph_params +
scale_color_gradientn(colors = rampPallFunVel(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 8),
limits=c(0, 0.055),
breaks=c(seq(0,0.055,by=.01)))
gg32 <- ggplot(data=predX1, aes(vel_abs, pred1, color = GD, group=ID)) +
vel_graph_params +
scale_color_gradientn(colors = rampPallFunVel(sigcolor$x),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 15),
limits=c(0, 0.055),
breaks=c(seq(0,0.055,by=.01))) 13.2.1 Group fig v1
p1 <- gg1 + theme(legend.position="none") + ggtitle("Trailing edge") + theme(plot.title = element_text(hjust = 0.5))
p2 <- gg2 + theme(legend.position="none") + ggtitle("Centroid") + theme(plot.title = element_text(hjust = 0.5))
p3 <- gg3 + theme(legend.position="none") + ggtitle("Leading edge") + theme(plot.title = element_text(hjust = 0.5))
legend <- cowplot::get_legend(
gg3 +
scale_color_gradientn(colors = rampPallFunVel(x=rep(1,nrow(sigcolor))),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 15),
limits=c(0, 0.055),
breaks=c(seq(0,0.055,by=.01))))## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.## Warning in get_plot_component(plot, "guide-box"): Multiple components found;
## returning the first one. To return all, use `return_all = TRUE`.(p1 + p2 + p3 + legend) +
plot_layout(axis_titles = "collect", nrow = 1, widths = c(1,1,1,.3)) + # common axes
plot_annotation(tag_levels = list(c('(a)','(b)','(c)',''))) # figure tags
13.2.2 Group fig v2
ranges <- ggplot_build(p1) %>%
pluck("layout", "panel_params", 1) %>%
`[`(c("x.range", "y.range"))
x <- (0.004 - ranges$x.range[1]) / (ranges$x.range[2] - ranges$x.range[1])
y <- (3 - ranges$y.range[1]) / (ranges$y.range[2] - ranges$y.range[1])
p1 <- gg12 +
ggtitle("Trailing edge")
p2 <- gg22 +
ggtitle("Centroid")
p3 <- gg32 +
ggtitle("Leading edge")
(p1 + p2 + p3) +
plot_layout(axis_titles = "collect", nrow = 1, widths = c(1,1,1)) + # common axes
plot_annotation(tag_levels = list(c('(a)','(b)','(c)'))) # figure tags
13.2.3 Group fig v3
ranges <- ggplot_build(p1) %>%
pluck("layout", "panel_params", 1) %>%
`[`(c("x.range", "y.range"))
x <- (0.5 - ranges$x.range[1]) / (ranges$x.range[2] - ranges$x.range[1])
y <- (5.5 - ranges$y.range[1]) / (ranges$y.range[2] - ranges$y.range[1])
p1 <- gg1 +
theme(legend.title=element_blank(),
plot.title = element_text(hjust = 0.5),
legend.position = c(x, y), legend.justification = c(0, 0)) +
ggtitle("Trailing edge")
p2 <- gg2 +
theme(legend.title=element_blank(),
plot.title = element_text(hjust = 0.5),
legend.position = c(x, y), legend.justification = c(0, 0)) +
ggtitle("Centroid")
p3 <- gg3 +
theme(legend.title=element_blank(),
plot.title = element_text(hjust = 0.5),
legend.position = c(x, y), legend.justification = c(0, 0)) +
ggtitle("Leading edge")
legend <- cowplot::get_legend(
gg3 +
scale_color_gradientn(colors = rampPallFunVel(x=rep(1,nrow(sigcolor))),
guide = guide_colourbar(title.position = "right",
barwidth = 1, barheight = 15),
limits=c(0, 7),
breaks=c(seq(0,7,by=1))))## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.## Warning in get_plot_component(plot, "guide-box"): Multiple components found;
## returning the first one. To return all, use `return_all = TRUE`.(p1 + p2 + p3 + legend) +
plot_layout(axis_titles = "collect", nrow = 1, widths = c(1,1,1,.3)) + # common axes
plot_annotation(tag_levels = list(c('(a)','(b)','(c)',''))) # figure tags