extinction_risk.qmd

# Suitable probability {#sec-suitableProbability}

```{r,include=FALSE,echo=FALSE}
Echo=FALSE
Eval=TRUE
Cache=TRUE
Warng=FALSE
Msg=FALSE
library(tidyverse)
library(ggdist)
library(ggpubr)
library(ggrepel)
library(ggdensity)
```

In this section, we discuss the concept and application of suitable probability derived from the stochastic population growth rate in the IPM.
The variability in the population growth rate ($\lambda$) arises from both environmental stochasticity and parametric uncertainty.
The former combines spatial and temporal variation in the climate (mean annual temperature and precipitation) and variation in the strength of competition for light.
The latter combines model uncertainty at the individual level for the growth, survival, and recruitment models, along with the spatial heterogeneity determined by the variance of these demographic models among local populations (plot level variability).
The inherent demographic uncertainty and environmental stochasticity introduce variance into $\lambda$, and the strength of this variance directly influences the suitable probability, especially for populations with lower performance or density [@Holt2005].
To illustrate the concept, consider two populations with equal mean growth rates ($\overline{\lambda}$) but varying standard deviations ($\sigma_\lambda$).
We define suitable probability as the area under the distribution for $\lambda \geq 1$ [@fig-lambda_sigma].

```{r re_sigma,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-lambda_sigma
#| fig-width: 6
#| fig-height: 4
#| fig-cap: "Distribution of 2000 draws of population growth rate $\\lambda$ with a mean of 1.05 and low (green) and high (yellow) standard deviations. The vertical dashed line represents the threshold between viable and non-viable populations."

# generate random population growth rates
set.seed(0.0)
expand_grid(
    mean = c(1.05),
    sd = c(0.05, 0.2)
  ) |>
  group_by(sd) |>
  expand_grid(draw = 1:2000) |>
  mutate(
    lambda = rnorm(n(), mean, sd)
  ) |>
  ggplot() +
  aes(lambda) +
  aes(fill = factor(sd)) +
  geom_density(alpha = 0.7, color = 'transparent') +
  geom_vline(xintercept = 1, linetype = 2) +
  scale_fill_manual(
    values = c('#5ab4ac', '#d8b365'),
    labels = c(
      expression(bar(lambda)~'='~'1.05; '~sigma[lambda]~'='~'0.05'),
      expression(bar(lambda)~'='~'1.05; '~sigma[lambda]~'='~'0.2')
    )
  ) +
  theme_classic() +
  labs(
    x = expression(lambda),
    y = '',
    fill = NULL,
    # subtitle = expression('Distribution of 2000 draws of'~lambda~'~'~plain(N(1.05, sd)))
  ) +
  scale_x_continuous(breaks = seq(0.5, 1.5, 0.5)) +
  theme(
    axis.ticks.y = element_blank(),
    axis.text.y = element_blank()
  )

```

To estimate suitable probability, we first determine the cumulative distribution function, $F(x)$, from the generic probability density function, $log(\lambda) = f(t)$, as follows:

$$
F_{\lambda}(x) = P(\lambda \le x) = \int_{-\infty}^{0} f(t)dt
$${#eq-cdf}

This function represents the cumulative distribution from $-\infty$ to $x$.
Subsequently, we define the suitable probability ($\Lambda$) as the inverse of the cumulative distribution function for $x = 1$:

$$
  \Lambda = 1 - F_{\lambda}(1)
$${#eq-sp}

@fig-lambda_sigma2 illustrates the relationship between $\lambda$ variability and $\Lambda$.
The higher the variability around $\lambda$, the more blurred the mean difference between two populations will become.
Increasing variance tends to reduce suitable probability when the population growth rate is below 1.

```{r re_sigma,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-lambda_sigma2
#| fig-width: 6
#| fig-height: 4
#| fig-cap: "Cumulative density function $F(\\lambda)$ defines suitable probability as a function of $\\lambda$ variability for different population growth rates. The smooth lines represent 2000 replications for varying mean ($\\lambda$) and standard deviation ($\\sigma$)."

set.seed(0.0)
expand_grid(
    mean = c(.6, 0.9, 1, 1.1, 1.4),
    sd = seq(0.15, 2, 0.01)
  ) |>
  mutate(group = row_number()) |>
  group_by(group) |>
  expand_grid(draw = 1:2000) |>
  mutate(
    lambda = rnorm(n(), mean, sd)
  ) ->
random_lambdas


# function to get the area under the curve from Threshold
area_thr <- function(x, threshold = 1, normalized = TRUE) {
  dd = density(x)
  den_weight =  dd$y
  if(normalized) {
    return( sum( (den_weight/sum(den_weight))[dd$x >= threshold] ) )
  }else{
    return( sum( den_weight[dd$x >= threshold] ) )
  }
}

random_lambdas |>
  group_by(mean, sd) |>
  reframe(ext = area_thr(lambda)) |>
  ggplot() +
  aes(sd, ext) +
  aes(color = factor(mean)) +
  geom_smooth() +
  scale_colour_viridis_d() +
  theme_classic() +
  labs(
    x = expression(sigma),
    y = expression('Suitable probability ('~Lambda~')'),
    color = expression(bar(lambda))
  )

# # using the normal distribution Z-score
# expand_grid(
#     mean = seq(0.9, 1.1, 0.001),
#     sd = seq(0.05, 2.5, 0.001)
#   ) |>
#   mutate(prop_above = 1 - pnorm((1-mean)/sd)) |>
#   ggplot() +
#   aes(mean, sd) +
#   aes(fill = prop_above) +
#   geom_raster() +
#   scale_fill_gradient2(midpoint = 0.5) +
#   theme_classic()
```


## Simulating the effect of demographic uncertainty and environmental stochasticity

```{r loadSim,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
lambdas <- readRDS(gzcon(url('https://github.com/willvieira/forest-IPM/raw/master/simulations/uncertainty_sim/output_processed/lambdas.RDS')))
```

In this section, we investigate how increasing the uncertainty of parameters and the variability of covariates affects the suitable probability.
We selected two species, *Quercus prinus* and *Quercus rubra*, with an average $\lambda$ close to 1, and were sensitive to competition and climate covariates.
We set these species to conditions with population growth rates closest to 1 but with different average values (lower for *Q. rubra* and higher for *Q. prinus*) for comparison [@fig-lambdaDistsp].

```{r lambdaDistsp,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-lambdaDistsp
#| fig-width: 6
#| fig-height: 4
#| fig-cap: "Population growth rate under demographic uncertainty from the 500 posterior distributions."

lambdas |>
  filter(noise_size == 'noise1') |>
  ggplot() +
  aes(log(lambda), species_name) +
  aes(fill = species_name) +
  stat_halfeye() +
  geom_vline(xintercept = 0, linetype = 2, alpha = 0.4) +
  theme_classic() +
  scale_fill_manual(
    values = c('#d8b365', '#5ab4ac')
  ) +
  labs(
    x = expression('ln('~lambda~')'),
    y = NULL
  ) +
  theme(
    legend.position = 'none',
    axis.text.y = element_text(face = 'italic')
  )

```

We used the same 500 parameter draws from the posterior distribution defined in the sensitivity analysis section (@sec-sensAnalysis).
We conducted four simulations to understand the cumulative impact of demographic uncertainty, environmental, and environmentally induced competition stochasticity.
First, we quantified suitable probability under demographic uncertainty using the posterior parameter distribution while keeping the covariates at their mean conditions, effectively simulating zero noise from the covariates.
Subsequently, we introduced environmental stochasticity and demographic uncertainty by adding noise to the mean temperature conditions.
Similarly, we introduced noise to the mean basal area of heterospecific individuals while keeping temperature constant at the mean.
Finally, we added all sources of variability together, encompassing demographic, environmental, and environmentally induced competition variability. 
We increased the standard deviation (SD) for temperature from 0.1 to 1.5, and for basal area, SD increased from 1 to 8.

@fig-lambdaVsNoise shows the changes in $\lambda$ variability with increased variability of the covariates.
The panels illustrate the effect of climate stochasticity, competition stochasticity, and both combined factors.

```{r lambdaVsNoise,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-lambdaVsNoise
#| fig-width: 8
#| fig-height: 4
#| fig-cap: "Distribution of population growth rate for different noise intensities applied to the temperature and basal area abundance covariates. We have omitted the values of $\\sigma$ as they differ between the covariates. However, the noise increases from left to right."

lambdas |>
  filter(lambda > quantile(lambda, 0.001)) |>
  ggplot() +
  aes(noise_size, log(lambda)) +
  aes(fill = species_name) +
  facet_wrap(~sim) +
  stat_slab(alpha = 0.6) +
  geom_hline(yintercept = 0, alpha = 0.4, linetype = 2) +
  scale_fill_manual(
    values = c('#d8b365', '#5ab4ac')
  ) +
  theme_classic() +
  labs(
    y = expression('ln('~lambda~')'),
    x = expression('Noise intensity ('~sigma~') increasing from left to right'),
    fill = ''
  ) +
  theme(
    legend.text = element_text(face = 'italic'),
    axis.text.x = element_blank(),
    panel.grid.major.x = element_line(color = rgb(0, 0, 0, 0.1)),
    legend.position = 'top'
  )
```

@fig-lambdaVsNoise and @fig-sdVsNoise highlight that increased noise around the average covariate conditions leads to higher variability in $\lambda$.
It is important to note that the y-scale on the panels in @fig-sdVsNoise is standardized, showing a higher $\lambda$ variability for *Q. rubra* compared to *Q. prinus*.
This suggests that climate and competition are more significant predictors for all vital rates of *Q. rubra* compared to *Q. prinus*.
In summary, higher variability in covariates will likely increase the variability in $\lambda$, proportional to the covariate sensitivity.

```{r sdVsNoise,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-sdVsNoise
#| fig-width: 8
#| fig-height: 4
#| fig-cap: "Correlation between $\\lambda$ standard deviation and variability intensity. The axis values of $\\sigma$ have been omitted as they differ between the covariates. However, the noise increases from left to right. The line represents a smoothed approximation from the point data."

lambdas |>
  group_by(species_name, sim, noise_size) |>
  reframe(
    lambda_mean = mean(lambda),
    lambda_sd = sd(lambda)
  ) |>
  # pivot_longer(cols = contains('lambda')) |>
  ggplot() +
  aes(noise_size, lambda_sd) +
  aes(group = sim, color = sim) +
  # geom_path(linewidth = .8, alpha = 0.3) +
  # geom_smooth(level = NA, alpha = 0) +
  geom_point(size = 2) +
  stat_smooth(geom = 'line', alpha = 0.4, size = 1) +
  facet_wrap(~species_name, scales = 'free') +
  scale_color_manual(
    values = c('#7fc97f', '#fb8072', '#fdc086')
  ) +
  theme_classic() +
  labs(
    y = expression(sigma[lambda]),
    x = expression('Noise intensity ('~sigma~') increasing from left to right'),
    color = NULL
  ) +
  theme(
    strip.text= element_text(face = 'italic'),
    strip.background = element_blank(),
    axis.text.x = element_blank(),
    legend.position = 'top'
  )
```

Once we determine the uncertainty in $\lambda$, we can calculate population-level suitable probability using the cumulative distribution function $F(\lambda)$.
@fig-noiseVsExt illustrates how suitable probability increases for *Q. prinus* and decreases for *Q. rubra* with increased variability of covariates.

```{r noiseVsExt,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-noiseVsExt
#| fig-width: 8
#| fig-height: 4
#| fig-cap: "Correlation between suitable probability and variability intensity. The axis values of $\\sigma$ have been omitted as they differ between the covariates. However, the noise increases from left to right. The line represents a smoothed approximation from the point data."

lambdas |>
  # filter(species_name != 'Picea mariana') |>
  group_by(species_name, sim, noise_size) |>
  reframe(A = area_thr(lambda)) |>
  ggplot() +
  aes(noise_size, A) +
  aes(group = sim, color = sim) +
  # geom_path(linewidth = .8, alpha = 0.3) +
  # geom_smooth(level = NA, alpha = 0) +
  geom_point(size = 2) +
  stat_smooth(geom = 'line', alpha = 0.4, size = 1) +
  facet_wrap(~species_name, scales = 'free') +
  scale_color_manual(
    values = c('#7fc97f', '#fb8072', '#fdc086')
  ) +
  theme_classic() +
  labs(
    y = 'Suitable probability',
    x = expression('Noise intensity ('~sigma~') increasing from left to right'),
    color = NULL
  ) +
  theme(
    strip.text= element_text(face = 'italic'),
    strip.background = element_blank(),
    axis.text.x = element_blank(),
    legend.position = 'top'
  )
```


```{r sigmaVsExt,echo=Echo,eval=FALSE,cache=Cache,warning=Warng,message=Msg}
#| label: fig-sigmaVsExt
#| fig-width: 8
#| fig-height: 4
#| fig-cap: ""

lambdas |>
  group_by(species_name, sim, noise_size) |>
  reframe(
    lambda_mean = mean(lambda),
    lambda_sd = sd(lambda),
    A = area_thr(lambda)
  ) |>
  ggplot() +
  aes(lambda_sd, A) +
  aes(color = sim) +
  stat_smooth(geom = 'line', alpha = 0.4, size = 1) +
  geom_point(size = 2) +
  facet_wrap(~species_name, scales = 'free') +
  scale_color_manual(
    values = c('#7fc97f', '#fb8072', '#fdc086')
  ) +
  theme_classic() +
  labs(
    y = 'Suitable probability',
    x = expression(sigma[lambda]),
    color = NULL
  ) +
  theme(
    strip.text= element_text(face = 'italic'),
    strip.background = element_blank(),
    legend.position = 'top'
  )
```

## Objective

This exploration shows that parametric uncertainty and environmental stochasticity can increase $\lambda$ variability and change suitable probability.
In this chapter, we want to assess the effect of spatio-temporal environmental variability and demographic uncertainty on the suitable probability of each species.
Specifically, we aim to model how the suitable probability changes from the center toward the cold and hot temperature range distribution.
By having the suitable probability varying across the temperature range, we can predict the effect of climate and competition, along with environmental variability and model uncertainty, on the range limits of forest trees.
This new approach allows us to test whether species interaction (here competition) can set range limits while accounting for the uncertainty of the estimations.
We conclude this study with a novel theory that uses suitable probability to establish the link between individual demographic rates and metapopulation dynamics.

## Methods

### Sources of population growth rate ($\lambda$) variability

The estimation of $\lambda$ occurs at the local population scale, specifically at the plot level in our study.
Within a given geographic location, such as a specific latitude where several plots are located, the variance of $\lambda$ among those plots arises from spatio-temporal variations in both climate and competition covariates.
For instance, climate stochasticity introduces noise in annual temperature and precipitation, leading to environmental variations. 
Similarly, even with identical climate conditions, two locations can exhibit different community abundance and composition, resulting in variability in the strength of competition. 
Beyond these spatio-temporal environmentally-induced variations, $\lambda$ can still vary due to model uncertainty.

In our study, demographic models track uncertainty at two complementary scales: individual and plot levels.
At the individual level, plots with the same climate and competition conditions may have different $\lambda$ values due to the uncertainty in the demographic models. 
Similarly, even with the same environmental conditions and averaged parameter values (eliminating demographic uncertainty at the individual level), two plots can still yield different $\lambda$ values due to the spatial uncertainty of each demographic model assigned among plots through the random effects.
Therefore, variability in the population growth rate can arise from spatio-temporal variations in both the environment and the parameters.

### Modeling suitable probability

We can evaluate the suitable probability of a species at various scales, ranging from a single local plot up to several plots in a region.
At the plot level, sources of variability in $\lambda$ stem from parametric uncertainty at the individual level and temporal environmental stochasticity in climate and competition.
When assessing suitable probability at the metapopulation scale by considering multiple plots simultaneously, we can additionally account for spatial variability in the environment and spatial uncertainty in plot-level parameters.

With the exception of parameter uncertainty at the individual level, all other sources of variability exhibit spatial dependence.
This implies that environmental variability (from climate, competition, or both) and parameter uncertainty at the plot level can vary based on their spatial location.
For instance, plots at the borders of the species distribution may experience more temperature variability than those at the center.
Additionally, plot-level parametric uncertainty can be spatially structured, capturing potential features of demographic variability beyond the climatic and competition covariates, such as historical factors or local edaphic conditions. 

At the metapopulation scale, we can model how suitable probability changes across the species' range distribution, considering both fixed climate and competition effects and the underlying spatio-temporal variability.
In our study, we are particularly interested in examining how suitable probability evolves from the center toward the cold and hot ranges. 
For that, we categorized all plot-year observations of a species based on the gradient of mean annual temperature (MAT), divided into cold and hot ranges using the MAT centroid among all plots ($\frac{max(MAT) + min(MAT)}{2}$).
We chose to use MAT instead of latitude because we are interested in the species' climatic niche, although the two variables are highly correlated.
The choice of MAT over latitude is motivated by our interest in understanding the species' climatic niche.
Despite our choice, MAT and latitude are highly correlated for all species.

We assessed suitable probability separately for the cold and hot ranges, employing a linear model to determine the relationship between $\lambda$ and MAT.
The spatio-temporal variability of $\lambda$ arising from environmental stochasticity and parameter uncertainty influences the variance of the linear model.
As this variance may change depending on the range position, we introduce a submodel for the variance of the linear model to be dependent on MAT.
To accommodate potential asymmetry in this variance, we use a Skew Normal Distribution ($\mathcal{SN}$) incorporating an additional parameter ($\alpha$) that can introduce right or left-skewed tails to the variance:

$$
\begin{split}
&log(\lambda) \sim \mathcal{SN}(\xi, \omega, \alpha) \\[2pt]
&\xi = \beta_{1, \xi} \times MAT + \beta_{0, \xi} \\[2pt]
&\omega = e^{\beta_{1, \omega} \times MAT + \beta_{0, \omega}}
\end{split}
$$ {#eq-metamodel}

Here, $\xi$ is the location parameter or the $\lambda$ average, and $\omega$ is the scale representing the variance around the mean.

In @fig-metamodel_example, we show some examples of the model's flexibility.
At each panel, two sets of parameter combinations illustrate the different behavior of the model.
Each panel illustrates an element of the model attributed to accommodate the data variability.
In panel 1, the most straightforward case shows how suitable probability changes linearly with MAT and constant variance.
Panel 2 shows two cases with larger and smaller variances derived from the spatiotemporal variability of the conditions and model uncertainty.
Panel 3 demonstrates how this variability can change across the range distribution.
Finally, panel 4 shows how the asymmetry of the variability can bias the variance towards lower (yellow) or higher (green) values.
Note that in panel 4, the illustration does not fully represent the asymmetry in the model, as the tool to create the curves does not allow asymmetric variance.
Therefore, while the figure illustrates a shift in the average and variance, we should instead see a shift only in the variance.

```{r example_of_how_nice_and_flexible_is_my_model,echo=FALSE,eval=FALSE}
library(tidyverse)

n = 50
tibble(
  beta_mean = rnorm(n, -0.5, 0.3),
  beta_sigma = rnorm(n, -.1, .1),
  inter = rnorm(n, 2, .5),
  sigma_inter = rnorm(n, -1, .8),
  alpha = rnorm(n, 6, 1)
) |>
mutate(
  row_id = row_number()
) |>
group_by(row_id) |>
expand_grid(temp = seq(0, 10)) |>
group_by(temp, row_id) |>
expand_grid(rep = 1:100) |>
mutate(
  y = sn::rsn(
    n(),
    xi = inter + beta_mean * temp,
    omega = exp(sigma_inter + beta_sigma * temp),
    alpha = alpha
  )
) ->
pred_dat

pred_dat |>
  pull(y) |>
  quantile(probs = c(0.01, .99)) ->
yLim

pred_dat |>
  group_split(row_id) |>
  map(
    ~ .x |>
      ggplot() +
      aes(temp, y) +
      ggdist::stat_lineribbon(.width = 0.9) +
      geom_hline(yintercept = 0, alpha = 0.7, linetype = 2) +
      theme_classic() +
      theme(legend.position = 'none') +
      ylim(yLim)
  ) ->
gg_plots

animation::saveGIF(
  expr = {
    walk(
      gg_plots,
      ~ print(.)
    )
  },
  movie.name = 'test.gif',
  interval = .5,
  ani.width = 600,
  ani.height = 400
)

# TODO:
# - each idea should be introduced incrementally
# - First the slope and intercept variation
# - Then the intercept to the variance and what it represents (spatio-temporal...)
# - Then show that there is also a slope to the variance in function of MAT
# - Finally, the skewness can introduce possible bias in the variance
```

```{r metamodel_example,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-metamodel_example
#| fig-width: 8
#| fig-height: 6
#| fig-cap: "The linear model from the correlation between $\\lambda$ and Mean annual temperature allows us to assess how suitable probability changes across the climate range. At each panel, two sets of parameter combination illustrate the different behaviour of the model. Each panel illustrates an element of the model attributed to accomodate the data variability."

library(tidyverse)
tibble(
  beta_mean = c(-.05, -.1,-.05, -.1,-.05, -.1, -.05, -.1),
  beta_sigma = c(0, 0, 0, 0, .05, -.15, .05, -.15),
  inter = c(.8, 1.4, .8, 1.4, .8, 1.4, .8, 1.4),
  sigma_inter = c(-1.8, -1.8, -2, -1, -2, -1, -2, -1),
  alpha = c(0, 0, 0, 0, 0, 0, -1, 1),
  group = factor(rep(1:2, 4)),
  plot_n = factor(rep(1:4, each = 2))
) |>
mutate(row_id = row_number()) |>
group_by(row_id) |>
expand_grid(temp = seq(0, 18, 0.1)) |>
group_by(temp, row_id) |>
expand_grid(replication = 1:2000) |>
mutate(
  y = sn::rsn(
    n(),
    xi = inter + beta_mean * temp,
    omega = exp(sigma_inter + beta_sigma * temp),
    alpha = alpha
  )
) |>
ggplot() +
aes(temp, y) +
aes(group = group, color = group, fill = group) +
facet_wrap(~plot_n) +
ggdist::stat_lineribbon(.width = 0.9, alpha = 0.7, size = 2) +
# geom_point(alpha = 0.5, size = .3) +
theme_classic() +
scale_color_manual(values = c('#d8b365', '#5ab4ac')) +
scale_fill_manual(values = c('#d8b365', '#5ab4ac')) +
geom_hline(yintercept = 0, linetype = 2, alpha = 0.8) +
labs(
  x = 'Mean annual temperature (°C)',
  y = expression('ln('~lambda~')'),
  color = NULL, fill = NULL
) +
theme(legend.position = 'none')
```

### Simulations

We computed $\lambda$ for each species based on the plot-year observations in the dataset, considering both environmentally induced variability and parameter uncertainty.
For every observed species-plot-year combination, we incorporated temporal stochasticity in climate conditions by using the mean and standard deviation of mean annual temperature and precipitation calculated from the years between measurements.
For instance, in the case of a plot observed twice, we calculated $\lambda$ for the second observation with climate conditions drawn randomly from a normal distribution with mean and standard deviation defined from climate observations within the year interval.
Similarly, temporal stochasticity in competition arises from variation in abundance and composition between measured years.
By iteratively performing this calculation, drawing parameter values randomly from the posterior distribution, we introduced demographic uncertainty at the individual level.
For each species-plot-year measurement, we replicated the calculation of $\lambda$ 100 times. 
By applying this approach across all plots, we naturally incorporate spatial variation in climate and competition conditions and spatial uncertainty in plot-level parameters.

For each species-plot-year-replication combination, we calculated $\lambda$ under two simulated conditions.
The first condition assumed no competition, with heterospecific competition defined at zero and conspecific competition at a minimal proportion (total population size ($N$) = 0.1).
This condition simulated the population growth rate under the fundamental niche.
The second condition simulated the invasion growth rate, representing the population growth rate when the focal species is rare ($N = 0.1$) and heterospecific competition set to the observed abundance of the competitive species.
This condition simulated the population growth rate under the realized niche.
The code for these computations on High-performance computing is available on GitHub: [https://github.com/willvieira/forest-IPM/tree/master/simulations/lambda_plot](https://github.com/willvieira/forest-IPM/tree/master/simulations/lambda_plot).

After calculating $\lambda$ for each species-plot-year-replication combination, we used these estimates to fit the linear model describing how $\lambda$ and its variability change across the mean annual temperature gradient.
We fitted the species-specific linear model for the hot and cold ranges using the Hamiltonian Monte Carlo (HMC) algorithm via the Stan software [version 2.30.1 @stan2022stan] and the `cmdstandr` R package [version 0.5.3 @cmdstanr].
To reduce the sampling time, we used a sample of 5000 plots for each species to fit the model.
This sample was necessary only for 6 out of the 31 species.
The R and Stan scripts for these computations on High-performance computing are accessible on GitHub: [https://github.com/willvieira/forest-IPM/tree/master/simulations/model_lambdaPlot/](https://github.com/willvieira/forest-IPM/tree/master/simulations/model_lambdaPlot/).

With the fitted models, we leveraged the posterior distribution to estimate the suitable probability of a species for any value of MAT under fundamental or realized niche for the cold and hot ranges.
Specifically, we estimated suitable probability under four different MAT conditions encountered by the species: at the border and the center of each cold and hot range.
We defined the border of the cold range as the minimum observed MAT for the focal species in the dataset, while the hot range was defined as the maximum observed MAT.
The center location is defined as the centroid of MAT for the focal species.
Although the center location has the same MAT for the cold and hot ranges, both are retained because the model is fitted separately for the cold and hot ranges.
Finally, we estimated suitable probability for each location under no competition (fundamental niche) and heterospecific competition (realized niche) conditions, using the empirical cumulative distribution function over 1000 predictive draws.

```{r check_for_spatial_aggreation_in_the_temperature_variability, echo=FALSE,eval=FALSE}
library(sf)

treeData = readRDS(file.path(readLines('_data.path'), 'growth_dt.RDS'))

treeData |>
  group_by(plot_id, year1) |>
  reframe(
    delta_year = unique(deltaYear),
    temp_var = unique(bio_01_sd)/delta_year,
    prec_var = unique(bio_12_sd)/delta_year,
    latitude = unique(latitude),
    longitude = unique(longitude)
  ) |>
  drop_na() |>
  st_as_sf(
    coords = c('longitude', 'latitude'),
    crs = 4326
  ) |>
  ggplot() + 
  aes(color = log(temp_var), fill = log(temp_var)) +
  geom_sf(size = 0.5, alpha = 0.7)

# define which species have border matching with the dataset border
# meaning it is not necessarily the species border
db_range <- treeData |>
  pull(bio_01_mean) |>
  # drop_na() |>
  range()

treeData |>
  filter(species_id %in% spIds$species_id_old) |>
  group_by(species_id) |>
  reframe(
    cold_diff = abs(min(bio_01_mean) - db_range[1]),
    hot_diff = abs(max(bio_01_mean) - db_range[2])
  ) |>
  ggplot() +
  aes(cold_diff, hot_diff) +
  geom_hline(yintercept = 0, linetype = 2, alpha = 0.7) +
  geom_vline(xintercept = 0, linetype = 2, alpha = 0.7) +
  geom_point() +
  ggrepel::geom_label_repel(aes(label = species_id), size = 2) +
  theme_classic()
```

## Results

```{r download_processed_data,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
out_pars <- readRDS(gzcon(url("https://github.com/willvieira/forest-IPM/raw/master/simulations/model_lambdaPlot/param_posterior.RDS"))) |>
  filter(par != 'lp__')
sp_dt <- readRDS(gzcon(url('https://github.com/willvieira/forest-IPM/raw/master/simulations/model_lambdaPlot/suitable_prob.RDS')))
sim_pars <- readRDS(gzcon(url('https://github.com/willvieira/forest-IPM/raw/master/simulations/model_lambdaPlot/simulation_pars.RDS')))

out_pars |>
  left_join(sim_pars) |>
  mutate(
    sim = factor(sim),
    cond = factor(cond)
  ) ->
out_pars
```

### Example using *Abies balsemea*

We used a linear model to summarize the evolution of the population growth rate ($\lambda$) and its variability across the cold and warm ranges (@eq-metamodel).
The @fig-res_example shows the observed $\lambda$ distribution and the fit of the underlying model on the mean annual temperature gradient for the species *Abies balsamea*.
Each point represents a plot-year-replication encompassing the complete spatio-temporal sources of variability arising from the stochastic environment and parameter uncertainty.
The black line represents the average fitted model (how $\lambda$ changes with MAT), and the shaded area is the 90th quantile of model uncertainty (how the variability of $\lambda$ changes with MAT).
From this uncertainty, we can deduce the suitable probability.
In this example, in the cold range, the mean and variance of $\lambda$ decrease towards the cold boundary.
By comparing the model under heterospecific competition with that without competition, we can observe a slight decrease in the average at the cold border but a larger decrease in uncertainty (@fig-res_example at the bottom left).

```{r res_example,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-res_example
#| fig-width: 8
#| fig-height: 6
#| fig-cap: "Distribution of stochastic population growth rate ($\\lambda$) for *Abies balsamea* over the mean annual temperature gradient for cold (left panels) and warm (right panels) ranges under no competition (fundamental niche) and heterospecific competition (realized niche). The dots represent $\\lambda$ over the plot-year-replication combinations. The model's average line and 90% prediction intervals are estimated using 500 draws from the posterior distribution."

Sp = '18032ABIBAL'

readRDS(file.path('data', 'lambda_MAT_dt.RDS')) |>
  list2env(envir = environment()) |>
  invisible()

out_pars |>
  filter(species_id == Sp & border == 'cold' & sim == 2 & cond != 3) |>
  select(cond, iter, par, value) |>
  pivot_wider(names_from = par, values_from = value) |>
  group_by(cond, iter) |>
  expand_grid(
    temp = seq(temp_range_cold[1], temp_range_cold[2], length.out = 200)
  ) |>
  mutate(
    lambda = sn::rsn(n(), beta_mean * temp + inter, exp(sigma_inter + beta_sigma * temp), alpha)
  ) |>
  ggplot() +
  aes(temp, lambda) +
  aes(fill = cond) +
  facet_wrap(
    ~cond,
    ncol = 1,
    labeller = as_labeller(setNames(c('No competition', 'Heterospecific competition'), 1:2))
  ) +
  geom_point(
    data = db_sim_cold |>
      mutate(cond = factor(cond)) |>
      group_by(cond) |>
      slice_sample(n = 1e4),
    aes(temp, log(lambda)),
    alpha = 0.4, size = 0.3
  ) +
  stat_lineribbon(.width = 0.9, alpha = 0.6) +
  # scale_color_manual(values = c('#5ab4ac', '#d8b365')) +
  scale_fill_manual(values = c('#5ab4ac', '#d8b365')) +
  geom_hline(yintercept = 0, linetype = 2) +
  theme_classic() +
  theme(
    plot.title = element_text(face = 'italic'),
    legend.position = 'none',
    strip.background = element_blank()
  ) +
  labs(
    subtitle = 'Cold',
    x = 'Mean annual temperature (°C)',
    y = expression('ln('~lambda~')')
  ) +
  xlim(temp_range_cold) +
  ylim(yLim) ->
p1

out_pars |>
  filter(species_id == Sp & border == 'hot' & sim == 2 & cond != 3) |>
  select(cond, iter, par, value) |>
  pivot_wider(names_from = par, values_from = value) |>
  group_by(cond, iter) |>
  expand_grid(
    temp = seq(temp_range_hot[1], temp_range_hot[2], length.out = 200)
  ) |>
  mutate(
    lambda = sn::rsn(n(), beta_mean * temp + inter, exp(sigma_inter + beta_sigma * temp), alpha)
  ) |>
  ggplot() +
  aes(temp, lambda) +
  aes(fill = cond) +
  facet_wrap(
    ~cond,
    ncol = 1,
    labeller = as_labeller(setNames(c('No competition', 'Heterospecific competition'), 1:2))
  ) +
  geom_point(
    data = db_sim_hot |>
      mutate(cond = factor(cond)) |>
      group_by(cond) |>
      slice_sample(n = 1e4),
    aes(temp, log(lambda)),
    alpha = 0.4, size = 0.3
  ) +
  stat_lineribbon(.width = 0.9, alpha = 0.6) +
  # scale_color_manual(values = c('#5ab4ac', '#d8b365')) +
  scale_fill_manual(values = c('#5ab4ac', '#d8b365')) +
  geom_hline(yintercept = 0, linetype = 2) +
  theme_classic() +
  theme(
    strip.background = element_blank(),
    legend.position = 'none'
  ) +
  labs(
    subtitle = 'Hot',
    x = 'Mean annual temperature (°C)',
    y = NULL
  ) +
  xlim(temp_range_hot) +
  ylim(yLim) ->
p2

ggarrange(p1, p2, ncol = 2)
```

We can estimate the suitable probability using the empirical cumulative distribution approach (@eq-sp) from the linear model predictions.
The @fig-sp-example shows the suitable probability expected over the mean annual temperature of the same species.
As expected, the suitable probability was reduced toward the cold border and was lower under heterospecific competition (yellow curve).
We can also observe that the decrease in suitable probability toward the border is nonlinear, becoming more substantial for heterospecific competition than for the no-competition condition.

```{r sp-example,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-sp-example
#| fig-width: 8
#| fig-height: 3.5
#| fig-cap: "Suitable probability of *Abies balsamea* over the mean annual temperature gradient for cold and warm ranges under no competition (green) and heterospecific (yellow). The vertical dotted line represents the range limits of the MAT observed in the dataset."

# equation to estimate suitable probability using the ecdf
get_ext <- function(mean_y, var_y, alpha) {
  P = ecdf(sn::rsn(1000, mean_y, var_y, alpha))
  return(1 - P(0))
}

Sp = '18032ABIBAL'

out_pars |>
  filter(species_id == Sp & border == 'cold' & sim == 2 & cond != 3) |>
  mutate(
    cond = case_match(
      as.character(cond),
      '1' ~ 'No competition',
      '2' ~ 'Heterospecific competition'
    )
  ) |>
  group_by(cond, par) |>
  reframe(value = mean(value)) |>
  pivot_wider(names_from = par, values_from = value) |>
  group_by(cond) |>
  expand_grid(
    temp = seq(species_range[1] - abs(species_range[1] * 0.2), temp_range_cold[2], length.out = 1000)
  ) |>
  mutate(
    mean = beta_mean * temp + inter,
    sd = exp(sigma_inter + beta_sigma * temp)
  ) |>
  rowwise() |>
  mutate(
    ext = get_ext(mean, sd, alpha)
  ) |>
  ggplot() +
  aes(temp, ext) +
  aes(color = cond) +
  geom_smooth(size = 1.3) +
  geom_vline(xintercept = species_range[1], linetype = 2, alpha = 0.6) +
  theme_classic() +
  scale_color_manual(values = c('#d8b365', '#5ab4ac')) +
  theme(
    plot.title = element_text(face = 'italic'),
    legend.position = 'top'
  ) +
  labs(
    subtitle = 'Cold',
    x = 'Mean annual temperature (°C)',
    y = 'Suitable probability',
    color = NULL
  ) +
  ylim(0, 1) ->
p1

out_pars |>
  filter(species_id == Sp & border == 'hot' & sim == 2 & cond != 3) |>
  mutate(
    cond = case_match(
      as.character(cond),
      '1' ~ 'No competition',
      '2' ~ 'Heterospecific competition'
    )
  ) |>
  group_by(cond, par) |>
  reframe(value = mean(value)) |>
  pivot_wider(names_from = par, values_from = value) |>
  group_by(cond) |>
  expand_grid(
    temp = seq(temp_range_hot[1], species_range[2] + species_range[2] * 0.2, length.out = 1000)
  ) |>
  mutate(
    mean_y = beta_mean * temp + inter,
    sd_y = exp(sigma_inter + beta_sigma * temp)
  ) |>
  rowwise() |>
  mutate(
    ext = get_ext(mean_y, sd_y, alpha)
  ) |>
  ggplot() +
  aes(temp, ext) +
  aes(color = cond) +
  geom_smooth(size = 1.3) +
  geom_vline(xintercept = species_range[2], linetype = 2, alpha = 0.6) +
  theme_classic() +
  scale_color_manual(values = c('#d8b365', '#5ab4ac')) +
  theme(
    axis.text.y = element_blank(),
    legend.position = 'top'
  ) +
  labs(
    subtitle = 'Hot',
    x = 'Mean annual temperature (°C)',
    y = NULL,
    color = NULL
  ) +
  ylim(0, 1) ->
p2

ggarrange(
  p1, p2,
  ncol = 2,
  common.legend = TRUE,
  legend = 'top'
)
```

### Effect of climate and competition on suitable probability

In the final step, we analyze the effect of competition on the suitable probability at the border and center of the temperature range distribution for all species.
The @fig-sp_comp_vs_nocomp2 compares the suitable probability under no competition to those under heterospecific competition for four locations from the mean annual temperature gradient.
Overall, suitable probability was high among the species, with an average of 0.78.
Among the four locations, species presented a lower suitable probability at the border of the hot range, with an average of 0.67.
Almost all species, in all conditions, are distributed below the identity line (1:1), meaning that heterospecific competition reduces the suitable probability.
In the cold range, temperature decrease toward the border had little effect on the suitable probability for both competitive conditions.
However, at the hot range, however, we can observe a significant shift in suitable probability from the center to the border.
Across the temperature range, there is a monotonic decrease in suitable probability from the cold border toward the hot border.

```{r sp_comp_vs_nocomp,echo=Echo,eval=FALSE,cache=Cache,warning=Warng,message=Msg}
#| label: fig-sp_comp_vs_nocomp
#| fig-width: 8
#| fig-height: 4
#| fig-cap: "The relationship in the suitable probability between no competition and with heterospecific competition. Each point represents the 31 forest species at their border (red) or center (blue) locations. Species points are grouped by a Multivariate Normal Density function with 75% probability. Note that we omitted the parameter uncertainty in this figure to avoid overlap and increase clarity."

sp_dt |>
  select(species_id, iter, cond, range_pos, border, ext) |>
  pivot_wider(
    names_from = cond,
    values_from = ext,
    names_prefix = 'cond_'
  ) |>
  mutate(
    sp_diff = cond_2-cond_1
  ) |>
  group_by(species_id, range_pos, border) |>
  reframe(
    cond_1 = mean(cond_1),
    cond_2 = mean(cond_2)
  ) |>
  ggplot() +
  aes(cond_1, cond_2) +
  aes(color = range_pos, fill = range_pos) +
  facet_wrap(~border) +
  geom_hdr(
    aes(color = NULL),
    probs = .75, alpha = .4,
    method = 'mvnorm'
  ) +
  geom_point() +
  geom_abline(slope = 1, intercept = 0, alpha = 0.7) +
  scale_color_manual(values = c('#fc8d59', '#91bfdb')) +
  scale_fill_manual(values = c('#fc8d59', '#91bfdb')) +
  theme_classic() +
  labs(
    x = 'Suitable probability (no competition)',
    y = 'Suitable probability (heter competition)',
    color = NULL,
    fill = NULL
  ) +
  theme(
    legend.position = 'top',
    strip.background = element_blank()
  )

```

```{r sp_comp_vs_nocomp2,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-sp_comp_vs_nocomp2
#| fig-width: 8
#| fig-height: 4
#| fig-cap: "Suitable probability for the 31 forest species in their distribution from the cold to the hot border of the mean annual temperature gradient. Note that we omitted the parameter uncertainty of each species in this figure to avoid overlap and increase clarity."

xLabel_order = c('Cold.border', 'Cold.center', 'Hot.center', 'Hot.border')

sp_dt |>
  select(species_id, iter, cond, range_pos, border, ext) |>
  mutate(
    cond = factor(cond, levels = c(1, 2), labels = c('No competition', 'Heterospecific competition'))
  ) |>
  group_by(species_id, border, cond, range_pos) |>
  reframe(
    ext = mean(ext)
  ) |>
  mutate(
    gp = factor(interaction(border, range_pos), levels = xLabel_order)
  ) ->
plot_dt

# get ABIBAL example to help undertand the figure
plot_dt |>
  filter(species_id == '18032ABIBAL') ->
abibal
 
plot_dt |>
  ggplot() +
  aes(gp, ext) +
  # aes(color = factor(cond, levels = c('No competition', 'Heterospecific competition'))) +
  aes(color = cond) +
  geom_point(position = position_jitterdodge()) +
  geom_boxplot(fill = 'transparent', outlier.color = 'transparent') +
  scale_color_manual(values = c('#5ab4ac', '#d8b365')) +
  geom_vline(xintercept = 2.5, linetype = 2, alpha = 0.7) +
  # add abibal
  geom_point(
    data = abibal,
    aes(gp, ext, fill = cond),
    color = 'red',
    position = position_jitterdodge()
  ) +
  theme_classic() +
  scale_x_discrete(
    breaks = xLabel_order,
    labels = c('Border', 'Center', 'Center', 'Border')
  ) +
  labs(
    y = expression('Suitable probability ('~Lambda~')'),
    x = "Species' relative position in the mean annual temperature gradient",
    color = NULL, fill = NULL
  ) +
  theme(
    legend.position = 'top',
    strip.background = element_blank()
  ) +
  annotate(geom = 'text', x = 1.5, y = 1.08, label = 'Cold') +
  annotate(geom = 'text', x = 3.5, y = 1.08, label = 'Hot') +
  coord_cartesian(ylim = c(0, 1), clip = 'off')
```

While the negative effect of climate on suitable probability is evident in the above figure, distinguishing the individual effects of competition remains challenging.
Hence, in @fig-effect_of_comp, we assessed the impact of competition on reducing suitable probability by subtracting the suitable probability under heterospecific competition from the suitable probability with no competition.
Even when isolating the effect of climate, the negative impact of heterospecific competition on suitable probability increases toward the border of the hot range.

```{r effect_of_comp,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-effect_of_comp
#| fig-width: 7
#| fig-height: 4
#| fig-cap: "The effect of competition on suitable probability across the four range positions. Each species point is color-coded based on its shade tolerance following @burns1990silvics."

spIds <- read_csv(
  file.path(
    readLines('_data.path'), 'species_id.csv'
  )
) |>
mutate(
  shade = factor(shade, levels = c('tolerant', 'intermediate', 'intolerant')),
  shade_sylvics = factor(shade_sylvics, levels = c('very-tolerant',
                                            'tolerant',
                                            'intermediate',
                                            'intolerant',
                                            'very-intolerant'
                                            )),
  succession = factor(succession, levels = c('pioneer',
                                              'intermediate',
                                              'climax'))
) |>
filter(sp_to_analyze)

sp_dt |>
  select(species_id, iter, cond, border, range_pos, ext) |>
  group_by(species_id, border, cond, range_pos) |>
  reframe(
    ext = mean(ext)
  ) |>
  pivot_wider(
    names_from = cond,
    values_from = ext,
    names_prefix = 'cond_'
  ) |>
   mutate(
    diff = cond_2 - cond_1,
    y = paste(border, range_pos, sep = ' '),
    y = factor(y, levels = rev(c('Hot border', 'Hot center', 'Cold center', 'Cold border')))
  ) |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  ggplot() +
  aes(y, diff) +
  geom_hline(yintercept = 0, alpha = 0.7) +
  geom_boxplot(fill = 'transparent', outlier.shape = NA) +
  geom_point(
    aes(fill = shade_sylvics),
    shape = 21, size = 2, alpha = .8, color = 'grey', position = position_jitterdodge(seed = 15)
  ) +
  scale_fill_manual(values = c("#20bc45", "#87bc45", "#edbf33", "#ea5545", "#ba0000")) +
  geom_vline(xintercept = 2.5, linetype = 2, alpha = 0.7) +
  theme_classic() +
  theme(
    legend.position = 'top'
  ) +
  labs(
    y = 'Effect of competition on suitable probability',
    x = '',
    fill = NULL
  )
```

For each species, we estimated the difference in suitable probability between the border and center (for climate) and between heterospecific and non-competition.
In @fig-diff_sp_raw, we compare the suitable probability for each species under the climate or competition conditions, where each line and their respective color determines the amount of change in suitable probability.

```{r diff_sp_raw,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-diff_sp_raw
#| fig-width: 8
#| fig-height: 7
#| fig-cap: "Difference in suitable probability between the center and border of the distribution (left panels) and no competition and heterospecific competition (right) for the cold (top) and hot (bottom) ranges. The color of the line connecting each species' conditions represents the intensity of change in suitable probability—the more intense the color, the greater the shift between conditions."

sp_dt |>
  filter(border == 'Cold') |>
  select(species_id, iter, cond, range_pos, ext) |>
  mutate(
    cond = factor(cond, levels = c(1, 2), labels = c('No competition', 'Heterospecific competition'))
  ) |>
  group_by(species_id, cond, range_pos) |>
  reframe(
    ext = mean(ext)
  ) |>
  pivot_wider(
    names_from = range_pos,
    values_from = ext
  ) |> 
   mutate(
    diff = border - center
  ) |>
  pivot_longer(
    cols = c(center, border),
    names_to = 'range_pos',
    values_to = 'ext'
  ) |>
  mutate(range_pos = factor(range_pos, levels = c('center', 'border'))) |>
  ggplot() +
  aes(range_pos, ext) +
  geom_jitter(width = 0.01, alpha = 0.8) +
  geom_line(aes(color = diff, group = species_id), alpha = 0.8) +
  facet_wrap(~cond) +
  scale_color_gradient2() +
  theme_classic() +
  labs(
    title = 'Cold',
    subtitle = 'Center vs border',
    y = expression('Suitable probability ('~Lambda~')'),
    x = NULL,
    color = NULL,
  ) +
  ylim(0.05, 1) +
  theme(
    strip.background = element_blank(),
    axis.text.x = element_blank()
  ) ->
p1

sp_dt |>
  filter(border == 'Cold') |>
  select(species_id, iter, cond, range_pos, ext) |>
  group_by(species_id, cond, range_pos) |>
  reframe(
    ext = mean(ext)
  ) |>
  pivot_wider(
    names_from = cond,
    values_from = ext,
    names_prefix = 'cond_'
  ) |>
   mutate(
    diff = cond_2 - cond_1
  ) |>
  pivot_longer(
    cols = contains('cond_'),
    names_to = 'cond',
    values_to = 'ext'
  ) |>
  mutate(
    range_pos = factor(range_pos, levels = c('center', 'border')),
    cond = factor(cond, levels = paste0('cond_', c(1, 2)), labels = c('No comp', 'Heter competition'))
  ) |>
  ggplot() +
  aes(cond, ext) +
  geom_jitter(width = 0.01, alpha = 0.8) +
  geom_line(aes(color = diff, group = species_id), alpha = 0.8) +
  facet_wrap(~range_pos) +
  scale_color_gradient2() +
  theme_classic() +
  labs(
    title = '',
    subtitle = 'No competition vs heterospecific competition',
    y = NULL,
    x = NULL,
    color = NULL,
  ) +
  ylim(0.05, 1) +
  theme(
    strip.background = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank()
  ) ->
p2

sp_dt |>
  filter(border == 'Hot') |>
  select(species_id, iter, cond, range_pos, ext) |>
  mutate(
    cond = factor(cond, levels = c(1, 2), labels = c('No competition', 'Heterospecific competition'))
  ) |>
  group_by(species_id, cond, range_pos) |>
  reframe(
    ext = mean(ext)
  ) |>
  pivot_wider(
    names_from = range_pos,
    values_from = ext
  ) |> 
   mutate(
    diff = border - center
  ) |>
  pivot_longer(
    cols = c(center, border),
    names_to = 'range_pos',
    values_to = 'ext'
  ) |>
  mutate(range_pos = factor(range_pos, levels = c('center', 'border'))) |>
  ggplot() +
  aes(range_pos, ext) +
  geom_jitter(width = 0.01, alpha = 0.8) +
  geom_line(aes(color = diff, group = species_id), alpha = 0.8) +
  facet_wrap(~cond) +
  scale_color_gradient2() +
  theme_classic() +
  labs(
    title = 'Hot',
    subtitle = 'Center vs border',
    y = expression('Suitable probability ('~Lambda~')'),
    x = NULL,
    color = NULL,
  ) +
  ylim(0.05, 1) +
  theme(
    strip.background = element_blank()
  ) ->
p3

sp_dt |>
  filter(border == 'Hot') |>
  select(species_id, iter, cond, range_pos, ext) |>
  group_by(species_id, cond, range_pos) |>
  reframe(
    ext = mean(ext)
  ) |>
  pivot_wider(
    names_from = cond,
    values_from = ext,
    names_prefix = 'cond_'
  ) |>
   mutate(
    diff = cond_2 - cond_1
  ) |>
  pivot_longer(
    cols = contains('cond_'),
    names_to = 'cond',
    values_to = 'ext'
  ) |>
  mutate(
    range_pos = factor(range_pos, levels = c('center', 'border')),
    cond = factor(cond, levels = paste0('cond_', c(1, 2)), labels = c('No comp', 'Heter competition'))
  ) |>
  ggplot() +
  aes(cond, ext) +
  geom_jitter(width = 0.01, alpha = 0.8) +
  geom_line(aes(color = diff, group = species_id), alpha = 0.8) +
  facet_wrap(~range_pos) +
  scale_color_gradient2() +
  theme_classic() +
  labs(
    title = '',
    subtitle = 'No competition vs heterospecific competition',
    y = NULL,
    x = NULL,
    color = NULL,
  ) +
  ylim(0.05, 1) +
  theme(
    strip.background = element_blank(),
    axis.text.y = element_blank()
  ) ->
p4

ggarrange(
  p1, p2, p3, p4,
  nrow = 2, ncol = 2,
  common.legend = TRUE,
  legend = 'bottom'
)
```

We further estimated the relative effect of climate and competition on changing suitable probability from the center to the border of the species distribution (@fig-diff_sp_hist).
Specifically, we quantified the distinct contributions of temperature and competition to changes in suitable probability.
Positive values of the relative difference signify an increase in suitable probability from the center towards the border, negative values indicate a decrease, and zero implies no net change.
In the hot range, most species exhibited negative values of the relative difference in suitable probability, indicating a decline in suitability due to higher temperature and intensified competition toward the hot border of the species distribution.
At the cold range, most species showed positive values for the competition effect, meaning a reduction in the effect of competition toward the cold border.
However, the climate effect at the cold range showed a less clear pattern, with some species experiencing an increase and others a decrease in suitable probability as temperature decreased toward the cold range.
Overall, the relative difference in suitable probability from the center toward the border was more influenced by climate rather than competition.

```{r diff_sp_hist,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-diff_sp_hist
#| fig-width: 8
#| fig-height: 3.5
#| fig-cap: "Difference in suitable probability for climate and competition effects over the cold and hot ranges. Negative values denote a decrease in species suitable probability from the center towards the distribution border, while positive values indicate an increase. Specifically, a negative value for climate at the hot (or cold) range signifies a reduction in suitable probability as temperature rises (or falls) towards the border."

sp_dt |>
  select(species_id, iter, border, cond, range_pos, ext) |>
  pivot_wider(
    names_from = c(border, cond, range_pos),
    values_from = ext
  ) |>
  mutate(
    diff_Cold_centerBorder_nc = Cold_1_border - Cold_1_center,
    diff_Cold_comp_center = Cold_2_center - Cold_1_center,
    diff_Cold_comp_border = Cold_2_border - Cold_1_border,
    diff_Cold_comp_nc = diff_Cold_comp_border - diff_Cold_comp_center,
    diff_Hot_centerBorder_nc = Hot_1_border - Hot_1_center,
    diff_Hot_comp_center = Hot_2_center - Hot_1_center,
    diff_Hot_comp_border = Hot_2_border - Hot_1_border,
    diff_Hot_comp_nc = diff_Hot_comp_border - diff_Hot_comp_center,
  ) |>
  select(
    species_id, iter, contains('_nc')
  ) |>
  pivot_longer(cols = contains('_nc')) |>
  mutate(
    name = gsub('diff_', '', name),
    species = species_id,
    species_id = NULL
  ) |>
  pivot_wider() |>
  pivot_longer(
    cols = contains('_'),
    names_sep = '_',
    names_to = c('border', 'sim', 'x')
  ) |>
  mutate(
    sim = factor(sim, levels = c('comp', 'centerBorder'), labels = c('Competition', 'Climate'))
  ) |>
  group_by(species, border, sim, x) |>
  reframe(value = mean(value)) |>
  ggplot() +
  aes(x = value, y = sim) +
  aes(fill = value) +
  facet_wrap(~border) +
  # stat_eye(alpha = 0.5, .width = 0, size = 7) +
  geom_point(size = 2.5, shape = 21, color = 'grey', position = position_jitterdodge(seed = 15)) +
  geom_boxplot(fill = 'transparent', outlier.color = 'transparent', color = 'black') +
  # geom_violin(fill = 'transparent')+
  scale_fill_gradient2(
    midpoint = 0,
    name = ''
  ) +
  geom_vline(xintercept = 0, linetype = 2, alpha = 0.3) +
  theme_classic() +
  theme(
    legend.position = 'none',
    strip.background = element_blank()
  ) +
  labs(
    x = expression('Relative difference in suitable probability between center and border ('~Delta*Lambda~')'),
    y = NULL,
    fill = NULL
  )
```

In our investigation of whether suitable probability decreases simultaneously at both borders (indicating a unimodal shape) or exhibits a linear pattern from one border to the other (@fig-diff_sp_hot_vs_cold), we found distinct patterns under the climate and competition effects.
For the climate effect, the majority of species demonstrated a decrease in suitable probability at one border while the other remained unchanged.
Additionally, a few species displayed a clear linear pattern of decreasing suitable probability from the cold to the hot border, with only one species (*Betula papyrifera*) showing a pronounced unimodal shape.
Conversely, under the competition effect, most species exhibited a decrease in suitable probability at the hot border and an increase at the cold border, indicating a linear rise in the impact of competition from the cold to the hot border of the distribution.

```{r diff_sp_hot_vs_cold,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-diff_sp_hot_vs_cold
#| fig-width: 8
#| fig-height: 4.5
#| fig-cap: "The relationship in the difference in suitable probability from the center to the border ($\\Delta\\Lambda$) between the cold and hot borders. A species at the bottom-left area (both $\\Delta\\Lambda$ are negative) indicates an unimodal decrease in suitable probability at the borders. Conversely, species at the top-right square exhibit an inverse unimodal shape.  In the top-left or bottom-right areas, suitable probability linearly decreases or increases from the cold toward the hot border, respectively. For the climate effect (left panel), species are colored following the centroid mean annual temperature among all observed plots. In the climate effect (left panel), species are color-coded based on the centroid of mean annual temperature among all observed plots. In the competition effect (right panel), species are classified by shade tolerance following @burns1990silvics."

## Species info
spIds <- read_csv(
  file.path(
    readLines('_data.path'), 'species_id.csv'
  )
) |>
mutate(
  shade = factor(shade, levels = c('tolerant', 'intermediate', 'intolerant')),
  shade_sylvics = factor(shade_sylvics, levels = c('very-tolerant',
                                            'tolerant',
                                            'intermediate',
                                            'intolerant',
                                            'very-intolerant'
                                            )),
  succession = factor(succession, levels = c('pioneer',
                                              'intermediate',
                                              'climax'))
) |>
filter(sp_to_analyze)

sp_dt |>
  select(species_id, iter, border, cond, range_pos, ext) |>
  pivot_wider(
    names_from = c(border, cond, range_pos),
    values_from = ext
  ) |>
  mutate(
    diff_Cold_centerBorder_nc = Cold_1_border - Cold_1_center,
    diff_Cold_comp_center = Cold_2_center - Cold_1_center,
    diff_Cold_comp_border = Cold_2_border - Cold_1_border,
    diff_Cold_comp_nc = diff_Cold_comp_border - diff_Cold_comp_center,
    diff_Hot_centerBorder_nc = Hot_1_border - Hot_1_center,
    diff_Hot_comp_center = Hot_2_center - Hot_1_center,
    diff_Hot_comp_border = Hot_2_border - Hot_1_border,
    diff_Hot_comp_nc = diff_Hot_comp_border - diff_Hot_comp_center,
  ) |>
  select(
    species_id, iter, contains('_nc')
  ) |>
  pivot_longer(cols = contains('_nc')) |>
  mutate(
    name = gsub('diff_', '', name),
    species = species_id,
    species_id = NULL
  ) |>
  pivot_wider() |>
  pivot_longer(
    cols = contains('_'),
    names_sep = '_',
    names_to = c('border', 'sim', 'x')
  ) |>
  mutate(
    sim = factor(sim, levels = c('centerBorder', 'comp'), labels = c('Climate', 'Competition'))
  ) |>
  group_by(species, border, sim) |>
  reframe(value = mean(value)) |>
  pivot_wider(names_from = border, values_from = value) |>
  left_join(
    sp_dt |>
      filter(range_pos == 'center' & border == 'Cold') |>
      select(species_id, temp) |>
      distinct() |>
      rename(species = species_id)
  ) ->
plot_dt

plot_dt |>
  filter(sim == 'Climate') |>
  ggplot() +
  aes(Hot, Cold) +
  aes(color = temp) +
  geom_point(size = 2) +
  scale_color_viridis_c() +
  geom_vline(xintercept = 0, linetype = 2, alpha = .5) +
  geom_hline(yintercept = 0, linetype = 2, alpha = 0.5) +
  theme_classic() +
  theme(
    legend.position = 'bottom',
    legend.key.height= unit(.8, 'cm'),
    strip.background = element_blank()
  ) +
  labs(
    subtitle = 'Climate',
    x = expression(Delta*Lambda~'- Hot'),
    y = expression(Delta*Lambda~'- Cold'),
    color = 'MAT (°C)'
  ) +
  xlim(-.5, .3) +
  ylim(-.5, .3) ->
p1

plot_dt |>
  left_join(
    spIds,
    by = c('species' = 'species_id_old')
  ) |>
  filter(sim == 'Competition') |>
  ggplot() +
  aes(Hot, Cold) +
  aes(color = shade_sylvics) +
  geom_point(size = 2) +
  scale_color_manual(values = c("#20bc45", "#87bc45", "#edbf33", "#ea5545", "#ba0000")) +
  geom_vline(xintercept = 0, linetype = 2, alpha = .5) +
  geom_hline(yintercept = 0, linetype = 2, alpha = 0.5) +
  theme_classic() +
  theme(
    legend.position = 'bottom',
    strip.background = element_blank()
  ) +
  labs(
    subtitle = 'Competition',
    x = expression(Delta*Lambda~'- Hot'),
    y = NULL,
    color = NULL
  ) +
  guides(color = guide_legend(nrow = 2)) +
  xlim(-.5, .3) +
  ylim(-.5, .3) ->
p2

ggarrange(p1, p2, nrow = 1)
```

### Suitable probability across the mean annual temperature gradient

We further investigated the potential dependency of climate and competition effects on the geographical location of the species.
Specifically, we hypothesized that the impact of climate might be more pronounced for species distributed in hotter temperatures than those in colder temperatures.
<!-- The @fig-sp_over_MAT1 shows the suitable probability distribution at the border for the cold and hot ranges over the MAT gradient.
We can observe that heterospecific competition decreases the suitable probability almost similarly for all species across the temperature gradient.
The only exception is for the cold range species located toward higher temperatures, which showed a relatively higher effect of competition than species located in cold temperatures (@fig-sp_over_MAT1 left panel). -->
In @fig-sp_diff_over_MAT, we evaluated the difference in suitable probability for climate and competition effects along the mean annual temperature gradient.
Notably, the only condition showing a significant trend was the climate effect in the cold range (blue group at the top panel).
Not surprisingly, suitable probability decreased as temperature declined toward the cold range, but this pattern was observed only for species located in colder conditions.
Conversely, the suitable probability increased when species were situated in hotter conditions as temperature decreased toward the cold range.

```{r sp_over_MAT1,echo=Echo,eval=FALSE,cache=Cache,warning=Warng,message=Msg}
#| label: fig-sp_over_MAT1
#| fig-width: 8
#| fig-height: 4
#| fig-cap: "Suitable probability at the cold and hot border for 31 species located over the mean annual temperature gradient. Species points are grouped by a Multivariate Normal Density function with 75% probability."

sp_dt |>
  filter(range_pos == 'border') |>
  mutate(
    cond = case_match(
      as.character(cond),
      '1' ~ 'No competition',
      '2' ~ 'Heterospecific competition'
    )
  ) |>
  ggplot() +
  aes(temp, ext) +
  aes(color = cond, fill = cond) +
  facet_wrap(~border, scales = 'free_x') +
  geom_hdr(
    aes(color = NULL),
    probs = .9, alpha = .4,
    method = 'mvnorm'
  ) +
  ylim(0, 1.1) +
  stat_pointinterval() +
  theme_classic() +
  theme(
    legend.position = 'top',
    strip.background = element_blank()
  ) +
  scale_color_manual(values = c('#d8b365', '#5ab4ac')) +
  scale_fill_manual(values = c('#d8b365', '#5ab4ac')) +
  labs(
    x = 'Mean annual temperature (°C)',
    y = 'Suitable probability at the border',
    color = NULL,
    fill = NULL
  )
```

```{r sp_diff_over_MAT,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-sp_diff_over_MAT
#| fig-width: 8
#| fig-height: 5.5
#| fig-cap: "Relative difference in suitable probability between the center and border for climate and competition for 31 species located over the mean annual temperature gradient. The difference in suitable probability for the climate effect is marginalized under the no-competition condition. Species points are grouped by a Multivariate Normal Density function with 75% probability."

sp_dt |>
  # filter(cond == 1) |>
  select(species_id, iter, cond, range_pos, border, ext) |>
  pivot_wider(
    names_from = c(range_pos, cond),
    values_from = ext,
    names_prefix = 'range_'
  ) |>
  mutate(
    comp_border = range_border_2 - range_border_1,
    comp_center = range_center_2 - range_center_1,
    diff_comp = comp_border - comp_center,
    diff_clim = range_border_1 - range_center_1
  ) |>
  select(species_id, border, contains('diff')) |>
  left_join(
    sp_dt |>
      select(species_id, border, range_pos, temp) |>
      rename(brd = border) |>
      distinct() |>
      pivot_wider(names_from = range_pos, values_from = temp) |>
      mutate(
        temp = (border + center)/2,
        size = border - center,
        size = if_else(brd == 'Cold', abs(size), size)
      ) |>
      select(species_id, brd, temp, size) |>
      rename(border = brd)
  ) |>
  pivot_longer(
    cols = contains('diff')
  ) |>
  mutate(
    name = factor(name, levels = c('diff_clim', 'diff_comp'), labels = c('Climate', 'Competition'))
  ) |>
  ggplot() +
  aes(temp, value) +
  aes(color = border, fill = border) +
  facet_wrap(~name, nrow = 2) +
  geom_hline(yintercept = 0, alpha = 0.7) +
  geom_hdr(
    aes(color = NULL),
    probs = .75, alpha = .4,
    method = 'mvnorm'
  ) +
  stat_pointinterval() +
  theme_classic() +
  scale_color_manual(values = c('#91bfdb', '#fc8d59')) +
  scale_fill_manual(values = c('#91bfdb', '#fc8d59')) +
  theme(
    # legend.position = 'top',
    strip.background = element_blank()
  ) +
  labs(
    x = 'Mean annual temperature centroid between center and border (°C)',
    y = expression('Relative difference in suitable probability between center and border ('~Delta*Lambda~')'),
    color = NULL, fill = NULL
  )

# sp_dt |>
#   filter(range_pos == 'border') |>
#   select(species_id, temp, iter, cond, range_pos, border, ext) |>
#   pivot_wider(
#     names_from = cond,
#     values_from = ext,
#     names_prefix = 'cond_'
#   ) |>
#   mutate(
#     sp_diff = cond_2-cond_1
#   ) |>
#   # filter(species_id != '19462FAGGRA') |>
#   ggplot() +
#   aes(temp, sp_diff) +
#   aes(color = border, fill = border) +
#   # facet_wrap(~border, scales = 'free_x') +
#   geom_hdr(
#     aes(color = NULL),
#     probs = .75, alpha = .4,
#     method = 'mvnorm'
#   ) +
#   geom_hline(yintercept = 0, linetype = 2, alpha = 0.7) +
#   stat_pointinterval() +
#   theme_classic() +
#   theme(
#     legend.position = 'top',
#     strip.background = element_blank()
#   ) +
#   scale_color_manual(values = c('#91bfdb', '#fc8d59')) +
#   scale_fill_manual(values = c('#91bfdb', '#fc8d59')) +
#   labs(
#     x = 'Mean annual temperature (°C)',
#     y = 'Difference in suitable probability',
#     color = NULL,
#     fill = NULL
#   )
```

```{r sp_diff_over_rangeSize,echo=Echo,eval=FALSE,cache=Cache,warning=Warng,message=Msg}

sp_dt |>
  # filter(cond == 1) |>
  select(species_id, iter, cond, range_pos, border, ext) |>
  pivot_wider(
    names_from = c(range_pos, cond),
    values_from = ext,
    names_prefix = 'range_'
  ) |>
  mutate(
    comp_border = range_border_2 - range_border_1,
    comp_center = range_center_2 - range_center_1,
    diff_comp = comp_border - comp_center,
    diff_clim = range_border_1 - range_center_1
  ) |>
  select(species_id, border, contains('diff')) |>
  left_join(
    sp_dt |>
      select(species_id, border, range_pos, temp) |>
      rename(brd = border) |>
      distinct() |>
      pivot_wider(names_from = range_pos, values_from = temp) |>
      mutate(
        temp = (border + center)/2,
        size = border - center,
        size = if_else(brd == 'Cold', abs(size), size)
      ) |>
      select(species_id, brd, size) |>
      rename(border = brd)
  ) |>
  pivot_longer(
    cols = contains('diff')
  ) |>
  mutate(
    name = factor(name, levels = c('diff_clim', 'diff_comp'), labels = c('Climate', 'Competition'))
  ) |>
  ggplot() +
  aes(size, value) +
  aes(color = border, fill = border) +
  facet_wrap(~name, nrow = 2) +
  geom_hline(yintercept = 0, alpha = 0.7) +
  geom_hdr(
    aes(color = NULL),
    probs = .75, alpha = .4,
    method = 'mvnorm'
  ) +
  geom_smooth(method = 'lm') +
  stat_pointinterval() +
  theme_classic() +
  scale_color_manual(values = c('#91bfdb', '#fc8d59')) +
  scale_fill_manual(values = c('#91bfdb', '#fc8d59')) +
  theme(
    # legend.position = 'top',
    strip.background = element_blank()
  ) +
  labs(
    x = 'Mean annual temperature centroid between center and border (°C)',
    y = 'Relative difference in suitable probability between center and border ',
    color = NULL, fill = NULL
  )

```

## Appendices

```{r appen,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| fig-width: 8
#| fig-height: 8

# data with temperature range limits of each species db and cold/hot ranges
tempRange_ls <- readRDS(file.path('data', 'temperatureRange_sp.RDS'))

# species id to use latin name
spIds <- read_csv(file.path(readLines('_data.path'), 'species_id.csv'))

for(Sp in unique(out_pars$species_id))
{
  # define temperature ranges
  temp_cold = c(
    tempRange_ls[[Sp]][['species_range']][1] - abs(tempRange_ls[[Sp]][['species_range']][1] * 0.2),
    tempRange_ls[[Sp]][['temp_range_cold']][2]
  )
  temp_hot = c(
    tempRange_ls[[Sp]][['temp_range_hot']][1],
    tempRange_ls[[Sp]][['species_range']][2] + tempRange_ls[[Sp]][['species_range']][2] * 0.2
  )

  out_pars |>
    filter(species_id == Sp & border == 'cold' & sim == 2 & cond != 3) |>
    select(cond, iter, par, value) |>
    pivot_wider(names_from = par, values_from = value) |>
    group_by(cond, iter) |>
    expand_grid(
      temp = seq(temp_cold[1], temp_cold[2], length.out = 25)
    ) |>
    mutate(
      lambda = sn::rsn(n(), beta_mean * temp + inter, exp(sigma_inter + beta_sigma * temp), alpha)
    ) ->
  dat_cold

  out_pars |>
    filter(species_id == Sp & border == 'hot' & sim == 2 & cond != 3) |>
    select(cond, iter, par, value) |>
    pivot_wider(names_from = par, values_from = value) |>
    group_by(cond, iter) |>
    expand_grid(
      temp = seq(temp_hot[1], temp_hot[2], length.out = 25)
    ) |>
    mutate(
      lambda = sn::rsn(n(), beta_mean * temp + inter, exp(sigma_inter + beta_sigma * temp), alpha)
    ) ->
  dat_hot

  yLim = quantile(c(dat_cold$lambda, dat_hot$lambda), probs = c(0.01, 0.999))
  
  dat_cold |>
    ggplot() +
    aes(temp, lambda) +
    aes(color = cond, fill = cond) +
    facet_wrap(
      ~cond,
      ncol = 1,
      scales = 'free_x',
      labeller = as_labeller(setNames(c('No competition', 'Heterospecific competition'), 1:2))
    ) +
    stat_lineribbon(.width = 0.9, alpha = 0.6) +
    scale_color_manual(values = c('#5ab4ac', '#d8b365')) +
    scale_fill_manual(values = c('#5ab4ac', '#d8b365')) +
    geom_hline(yintercept = 0, linetype = 2) +
    geom_vline(xintercept = tempRange_ls[[Sp]][['species_range']][1], linetype = 2, alpha = 0.6) +
    theme_classic() +
    theme(
      plot.title = element_text(face = 'italic'),
      legend.position = 'none',
      strip.background = element_blank(),
      axis.text.x = element_blank()
    ) +
    labs(
      title = spIds |> filter(species_id_old == Sp) |> pull(species_name),
      subtitle = 'Cold',
      x = NULL,
      y = expression('ln('~lambda~')')
    ) +
    scale_y_continuous(
      limits = yLim,
      labels = scales::number_format(accuracy = 0.01)
    ) ->
  p1

  dat_hot |>
    ggplot() +
    aes(temp, lambda) +
    aes(color = cond, fill = cond) +
    facet_wrap(
      ~cond,
      ncol = 1,
      scales = 'free_x',
      labeller = as_labeller(setNames(c('No competition', 'Heterospecific competition'), 1:2))
    ) +
    stat_lineribbon(.width = 0.9, alpha = 0.6) +
    geom_hline(yintercept = 0, linetype = 2) +
    geom_vline(xintercept = tempRange_ls[[Sp]][['species_range']][2], linetype = 2, alpha = 0.6) +
    scale_color_manual(values = c('#5ab4ac', '#d8b365')) +
    scale_fill_manual(values = c('#5ab4ac', '#d8b365')) +
    theme_classic() +
    theme(
      strip.background = element_blank(),
      legend.position = 'none',
      axis.text.y = element_blank(),
      axis.text.x = element_blank()
    ) +
    labs(
      title = '',
      subtitle = 'Hot',
      x = NULL,
      y = NULL
    ) +
    scale_y_continuous(
      limits = yLim,
      labels = scales::number_format(accuracy = 0.01)
    ) ->
  p2

  out_pars |>
    filter(species_id == Sp & border == 'cold' & sim == 2 & cond != 3) |>
    group_by(cond, par) |>
    reframe(value = mean(value)) |>
    pivot_wider(names_from = par, values_from = value) |>
    group_by(cond) |>
    expand_grid(
      temp = seq(temp_cold[1], temp_cold[2], length.out = 1000)
    ) |>
    mutate(
      mean = beta_mean * temp + inter,
      sd = exp(sigma_inter + beta_sigma * temp)
    ) |>
    rowwise() |>
    mutate(
      ext = get_ext(mean, sd, alpha)
    ) |>
    ggplot() +
    aes(temp, ext) +
    aes(color = cond) +
    geom_smooth(size = 1.3) +
    geom_vline(xintercept = tempRange_ls[[Sp]][['species_range']][1], linetype = 2, alpha = 0.6) +
    theme_classic() +
    scale_color_manual(values = c('#5ab4ac', '#d8b365')) +
    theme(
      plot.title = element_text(face = 'italic'),
      legend.position = 'none'
    ) +
    labs(
      x = 'Mean annual temperature (°C)',
      y = 'Suitable probability'
    ) +
    scale_y_continuous(
      limits = c(0, 1),
      labels = scales::number_format(accuracy = 0.01)
    ) ->
  p3

  out_pars |>
    filter(species_id == Sp & border == 'hot' & sim == 2 & cond != 3) |>
    group_by(cond, par) |>
    reframe(value = mean(value)) |>
    pivot_wider(names_from = par, values_from = value) |>
    group_by(cond) |>
    expand_grid(
      temp = seq(temp_hot[1], temp_hot[2], length.out = 1000)
    ) |>
    mutate(
      mean_y = beta_mean * temp + inter,
      sd_y = exp(sigma_inter + beta_sigma * temp)
    ) |>
    rowwise() |>
    mutate(
      ext = get_ext(mean_y, sd_y, alpha)
    ) |>
    ggplot() +
    aes(temp, ext) +
    aes(color = cond) +
    geom_smooth(size = 1.3) +
    geom_vline(xintercept = tempRange_ls[[Sp]][['species_range']][2], linetype = 2, alpha = 0.6) +
    theme_classic() +
    scale_color_manual(values = c('#5ab4ac', '#d8b365')) +
    theme(
      legend.position = 'none',
      axis.text.y = element_blank()
    ) +
    labs(
      x = 'Mean annual temperature (°C)',
      y = NULL
    ) +
    scale_y_continuous(
      limits = c(0, 1),
      labels = scales::number_format(accuracy = 0.01)
    ) ->
  p4

  print(
    ggarrange(
      p1, p2, p3, p4,
      ncol = 2, nrow = 2,
      heights = c(1, 0.5)
    )
  )
}
```