pars_climate.qmd

# Climate effect {#sec-parsClimate}

In this section, I discuss the effect of climate on growth, survival, and recruitment vital rates.
We employ the bioclimatic variables of mean annual temperature (MAT) and mean annual precipitation (MAP), averaging these variables over the years, covering the time intervals spanning two measurement periods.
The effect of these two variables on growth, survival, and recruitment is modeled using an unimodal function characterized by two key parameters: the optimal climate ($\xi$), representing the climatic conditions where growth and survival reach their peak, and the climate breadth ($\sigma$), which quantifies the rate of change in the climate's impact as one moves away from this optimal point.
In the fit process, we transformed the climate breadth parameters as $\sigma = \tau^{-1}$ for sampling efficiency.

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

```{r intercept_loadPars,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
data_path <- readLines('_data.path')
pars_path <- file.path(data_path, 'output_sim_processed')
models <- c(
  'growth' = 'intcpt_plot_comp_clim',
  'mort' = 'intcpt_plot_comp_clim',
  'recruit' = 'intcpt_plot_comp_clim'
)

## Species if info
spIds <- read_csv(
  file.path(
    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)


# load data to unscale climate variables and function to scale climate back to original range
clim_scaleRange <- readRDS(file.path(data_path, 'climate_scaleRange.RDS')) |>
  bind_cols()

unscaleClim <- function(
  value,
  # if value is temp or prec
  clim,
  # temp and prec range are vectors of format [min, max]
  temp_rg = clim_scaleRange$bio_01_mean,
  prec_rg = clim_scaleRange$bio_12_mean
) {
  if(clim == 'temp') {
    return( value * (temp_rg[2] - temp_rg[1]) + temp_rg[1] )
  }else if(clim == 'prec') {
    return( value * (prec_rg[2] - prec_rg[1]) + prec_rg[1] )
  } else{
    stop('Clim must be either `temp` or `prec` character.')
  }
}


# Load parameters
map_dfr(
    spIds$species_id_old,
    ~ readRDS(
        paste0(
          pars_path, '/growth/', models['growth'], '/posterior_pop_', .x, '.RDS'
        )
      ) |>
    pivot_wider(names_from = par) |>
    select(iter, c(contains('optimal'), contains('tau'))) |>
    bind_cols(species_id = .x)
  ) ->
pars_growth

map_dfr(
    spIds$species_id_old,
    ~ readRDS(
      paste0(
        pars_path, '/mort/', models['mort'], '/posterior_pop_', .x, '.RDS'
      )
    ) |>
    pivot_wider(names_from = par) |>
    select(iter, c(contains('optimal'), contains('tau'))) |>
    bind_cols(species_id = .x)
  ) ->
  pars_survival

map_dfr(
    spIds$species_id_old,
    ~ readRDS(
      paste0(
        pars_path, '/recruit/', models['recruit'], '/posterior_pop_', .x, '.RDS'
      )
    ) |>
    pivot_wider(names_from = par) |>
    select(iter, c(contains('optimal'), contains('tau'))) |>
    bind_cols(species_id = .x)
  ) ->
  pars_recruit

treeData <- readRDS(file.path(data_path, 'treeData.RDS')) |>
  filter(species_id %in% spIds$species_id_old)
```


## Optimal climate

@fig-matDist and @fig-mapDist show the distribution of the optimal climate parameter, meaning the climate location where the vital rate is at its maximum, across the 31 tree species.
The parameter mean and their uncertainty are displayed with the dot interval for growth (green), recruitment (yellow), and survival (brown).
The density plot in light gray is the distribution of the climate variable among all observed trees across space and time.

#### Mean annual temperature

```{r climate_optClim,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-matDist
#| fig-width: 8.5
#| fig-height: 7
#| fig-cap: "Distribution for the optimal annual mean temperature ($\\xi_{MAT}$) for growth (green), recruitment (yellow), and survival (brown). The gray density plot is the annual mean temperature distribution among all observed trees across space and time."

treeData |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  filter(!is.na(bio_01_mean_scl)) |>
  ggplot() +
  aes(unscaleClim(bio_01_mean_scl, 'temp'), fct_reorder(species_name, bio_01_mean_scl)) +
  ggridges::geom_density_ridges2(color = NA, alpha = 0.5) +
  stat_pointinterval(
    data = pars_growth |>
      group_by(species_id) |>
      mutate(growth = optimal_temp, iter = 1:n()) |>
      select(species_id, iter, growth) |>
      ungroup() |>
      left_join(
        pars_survival |>
        group_by(species_id) |>
        mutate(survival = optimal_temp, iter = 1:n()) |>
        select(species_id, iter, survival) |>
        ungroup()
      ) |>
      left_join(
        pars_recruit |>
        group_by(species_id) |>
        mutate(recruit = optimal_temp, iter = 1:n()) |>
        select(species_id, iter, recruit) |>
        ungroup()
      ) |>
      pivot_longer(
        cols = c(growth, survival, recruit)
      ) |>
      left_join(
        spIds,
        by = c('species_id' = 'species_id_old')
      ),
    aes(unscaleClim(value, 'temp'), species_name, color = name),
    alpha = 0.95,
    position = position_dodge(width = 0.4) # dodge to avoid overlap
  ) +
  theme_classic() +
  scale_color_manual(
    values = c('#018571', '#dfc27d', '#a6611a'),
    labels = c('Growth', 'Recruitment', 'Survival')
  ) +
  theme(
    axis.text.y = element_text(face = "italic"),
    panel.grid.major.y = element_line(colour = rgb(0,0,0,.1))
  ) +
  labs(
    color = expression(xi[MAT]),
    x = expression(paste("Mean annual temperature (", degree, "C)")),
    y = ''
  )
```

#### Mean annual precipitation

```{r climate_optPrec,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-mapDist
#| fig-width: 8.5
#| fig-height: 7
#| fig-cap: "Distribution for the optimal mean annual precipitation ($\\xi_{MAP}$) for growth (green), recruitment (yellow), and survival (brown). The density plot in gray is the distribution of annual precipitation variable among all observed trees across space and time."

treeData |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  ggplot() +
  aes(unscaleClim(bio_12_mean_scl, 'prec'), fct_reorder(species_name, bio_12_mean_scl)) +
  ggridges::geom_density_ridges2(color = NA, alpha = 0.5) +
  stat_pointinterval(
    data = pars_growth |>
      group_by(species_id) |>
      mutate(growth = optimal_prec, iter = 1:n()) |>
      select(species_id, iter, growth) |>
      ungroup() |>
      left_join(
        pars_survival |>
        group_by(species_id) |>
        mutate(survival = optimal_prec, iter = 1:n()) |>
        select(species_id, iter, survival) |>
        ungroup()
      ) |>
      left_join(
        pars_recruit |>
        group_by(species_id) |>
        mutate(recruit = optimal_prec, iter = 1:n()) |>
        select(species_id, iter, recruit) |>
        ungroup()
      ) |>
      pivot_longer(
        cols = c(growth, survival, recruit)
      ) |>
      left_join(
        spIds,
        by = c('species_id' = 'species_id_old')
      ),
    aes(unscaleClim(value, 'prec'), species_name, color = name),
    alpha = 0.95,
    position = position_dodge(width = 0.4) # dodge to avoid overlap
  ) +
  theme_classic() +
  scale_color_manual(
    values = c('#018571', '#dfc27d', '#a6611a'),
    labels = c('Growth', 'Recruitment', 'Survival')
  ) +
  xlab('Annual Precipitation (mm)') +
  ylab('') +
  theme(
    axis.text.y = element_text(face = "italic"),
    panel.grid.major.y = element_line(colour = rgb(0,0,0,.1))
  ) +
  labs(color = expression(xi[MAP]),)
```

#### Optimal temperature vs temperature mid range

```{r optTempVsMid,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-optTempVsMid
#| fig-width: 8
#| fig-height: 4
#| fig-cap: "Optimal Mean annual temperature for the growth, survival, and recruitment models as function of each species midpoint temperature range."

treeData |>
  # mutate(
  #   bio_01_mean_nat = unscaleClim(bio_01_mean_scl, 'temp'),
  #   bio_12_mean_nat = unscaleClim(bio_12_mean_scl, 'prec')
  # ) |>
  group_by(species_id) |>
  reframe(
    mid_temp = (max(bio_01_mean, na.rm = TRUE)+min(bio_01_mean, na.rm = TRUE))/2,
    mid_prec = (max(bio_12_mean, na.rm = TRUE)+min(bio_12_mean, na.rm = TRUE))/2
  ) ->
sp_mid

pars_growth |>
  bind_cols(vr = 'Growth') |>
  bind_rows(
    pars_survival |>
      bind_cols(vr = 'Survival')
  ) |>
  bind_rows(
    pars_recruit |>
      bind_cols(vr = 'Recruitment')
  ) |>
  left_join(sp_mid) |>
  ggplot() +
  aes(mid_temp, unscaleClim(optimal_temp, 'temp')) +
  aes(color = vr) +
  stat_pointinterval(position = position_dodgejust(0.5)) +
  geom_smooth(method = 'lm') +
  scale_color_manual(
    values = c('#018571', '#dfc27d', '#a6611a'),
    labels = c('Growth', 'Recruitment', 'Survival')
  ) +
  theme_classic() +
  labs(
    x = 'Annual Mean Temperature midpoint (°C)',
    y = expression(xi[MAT]),
    color = ''
  )
```


#### Optimal precipitaiton vs precipitation mid range

```{r optPrecVsMid,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-optPrecVsMid
#| fig-width: 8
#| fig-height: 4
#| fig-cap: "Optimal Mean annual temperature for the growth, survival, and recruitment models as function of each species midpoint temperature range."

treeData |>
  # mutate(
  #   bio_01_mean_nat = unscaleClim(bio_01_mean_scl, 'temp'),
  #   bio_12_mean_nat = unscaleClim(bio_12_mean_scl, 'prec')
  # ) |>
  group_by(species_id) |>
  reframe(
    mid_temp = (max(bio_01_mean, na.rm = TRUE)+min(bio_01_mean, na.rm = TRUE))/2,
    mid_prec = (max(bio_12_mean, na.rm = TRUE)+min(bio_12_mean, na.rm = TRUE))/2
  ) ->
sp_mid

pars_growth |>
  bind_cols(vr = 'Growth') |>
  bind_rows(
    pars_survival |>
      bind_cols(vr = 'Survival')
  ) |>
  bind_rows(
    pars_recruit |>
      bind_cols(vr = 'Recruitment')
  ) |>
  left_join(sp_mid) |>
  ggplot() +
  aes(mid_prec, unscaleClim(optimal_prec, 'prec')) +
  aes(color = vr) +
  stat_pointinterval(position = position_dodgejust(0.5)) +
  geom_smooth(method = 'lm') +
  scale_color_manual(
    values = c('#018571', '#dfc27d', '#a6611a'),
    labels = c('Growth', 'Recruitment', 'Survival')
  ) +
  theme_classic() +
  labs(
    x = 'Annual Mean Precipitation midpoint (°C)',
    y = expression(xi[MAP]),
    color = ''
  )
```

## Climate breadth

The climate breadth is the second parameter to account for the climate effect on the growth, survival, and recruitment functions.
We define climate breadth as the precision parameter, which is the inverse of the variance of the bell shape ($\tau = \frac{1}{\sigma^2}$).
The closer this value is to zero, the higher is the breadth around the mean.
In other words, when climate breadth is zero, the bell-shaped unimodal function becomes an almost flat line.

#### Climate breadth 

```{r climate_tauSp,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-climBreadth
#| fig-width: 8.5
#| fig-height: 6
#| fig-cap: "Mean annual temperature and mean annual precipitation effect size across the growth (green), recruitment (yellow), and survival (brown) vital rates. The higher the value, the stronger the negative effect of climate when moving from the optimal climate."

pars_growth |>
  mutate(temp_growth = tau_temp,
          prec_growth = tau_prec) |>
  group_by(species_id) |>
  mutate(iter = 1:n()) |>
  select(species_id, iter, contains('growth')) |>
  ungroup() |>
  left_join(
    pars_survival |>
      mutate(temp_survival = tau_temp,
            prec_survival = tau_prec) |>
      group_by(species_id) |>
      mutate(iter = 1:n()) |>
      select(species_id, iter, contains('survival')) |>
      ungroup()
  ) |>
  left_join(
    pars_recruit |>
      mutate(temp_recruit = tau_temp,
            prec_recruit = tau_prec) |>
      group_by(species_id) |>
      mutate(iter = 1:n()) |>
      select(species_id, iter, contains('recruit')) |>
      ungroup()
  ) |>
  pivot_longer(
    cols = c(contains('temp_'), contains('prec_')),
    names_to = c('clim', 'vr'),
    names_sep = '_'
  ) |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  mutate(
    clim = case_match(
      clim,
      "temp" ~ "Mean annual temperature",
      "prec" ~ "Mean annual precipitation"
    )
  ) |>
  ggplot() +
  aes(value, fct_reorder(species_name, value)) +
  aes(color = vr) +
  facet_wrap(~clim) +
  stat_pointinterval() +
  theme_classic() +
  scale_color_manual(
    values = c('#018571', '#dfc27d', '#a6611a'),
    labels = c('Growth', 'Recruitment', 'Survival')
  ) +
  xlab(expression(tau)) +
  ylab('') +
  theme(axis.text.y = element_text(face = "italic")) +
  labs(color = '')
```


```{r climate_tauSumm,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-climBreadthSumm
#| fig-width: 8
#| fig-height: 5
#| fig-cap: "Mean annual temperature and mean annual precipitation effect size across the growth (green), recruitment (yellow), and survival (brown) vital rates. The higher the value, the stronger the negative effect of climate when moving from the optimal climate."

pars_growth |>
  group_by(species_id) |>
  reframe(
    temp_growth = mean(tau_temp),
    prec_growth = mean(tau_prec)
  ) |>
  select(species_id, contains('growth')) |>
  ungroup() |>
  left_join(
    pars_survival |>
      group_by(species_id) |>
      reframe(
        temp_survival = mean(tau_temp),
        prec_survival = mean(tau_prec)
      ) |>
      select(species_id, contains('survival')) |>
      ungroup()
  ) |>
  left_join(
    pars_recruit |>
      group_by(species_id) |>
      reframe(
        temp_recruit = mean(tau_temp),
        prec_recruit = mean(tau_prec)
      ) |>
      select(species_id, contains('recruit')) |>
      ungroup()
  ) |>
  pivot_longer(
    cols = c(contains('temp_'), contains('prec_')),
    names_to = c('clim', 'vr'),
    names_sep = '_'
  ) |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  mutate(
    clim = case_match(
      clim,
      "temp" ~ "Mean annual temperature",
      "prec" ~ "Mean annual precipitation"
    ),
    vr = case_match(
      vr,
      'growth' ~ 'Growth',
      'recruit' ~ 'Recruitment',
      'survival' ~ 'Survival'
    )
  ) |>
  ggplot() +
  aes(vr, value) +
  aes(fill = clim) +
  geom_boxplot(alpha = 0.8) +
  theme_classic() +
  scale_fill_manual(
    values = c('#386cb0', '#fdc086')
  ) +
  ylab(expression(tau)) +
  xlab('') +
  theme(axis.text.y = element_text(face = "italic")) +
  labs(fill = '') +
  theme(legend.position = 'bottom')
```

<!-- 
#### Climate breadth vs optimal climate

We have seen with @fig-matDist and @fig-mapDist that the optimal temperature and precipitation are often at the border of the distribution. 
In @fig-optClimVsTau, we plot the correlation between the optimal climate's position with its distribution (0 representing the lower limit and 1 indicating the upper limit) and the climate breadth parameter. -->

```{r climate_tauvsRangeSize,echo=Echo,eval=FALSE,cache=Cache,warning=Warng,message=Msg}
#| label: fig-optClimVsTau
#| fig-width: 9
#| fig-height: 6
#| fig-cap: "Climate breadth in function of mean optimal climate relative to the species climate range in which 0 means lower climate limit and 1 means upper climate limit."

treeData |>
  # mutate(
  #   bio_01_mean_nat = unscaleClim(bio_01_mean_scl, 'temp'),
  #   bio_12_mean_nat = unscaleClim(bio_12_mean_scl, 'prec')
  # ) |>
  group_by(species_id) |>
  reframe(
    max_temp = quantile(bio_01_mean_scl, 1, na.rm = TRUE),
    min_temp = quantile(bio_01_mean_scl, 0, na.rm = TRUE),
    max_prec = quantile(bio_12_mean_scl, 1, na.rm = TRUE),
    min_prec = quantile(bio_12_mean_scl, 0, na.rm = TRUE)
  ) ->
sp_range

pars_growth |>
  bind_cols(vr = 'growth') |>
  bind_rows(
    pars_survival |>
      bind_cols(vr = 'survival')
  ) |>
  bind_rows(
    pars_recruit |>
      bind_cols(vr = 'recruitment')
  ) |>
  left_join(sp_range) |>
  mutate(
    optimal_temp_scl = (optimal_temp - min_temp)/(max_temp - min_temp),
    optimal_prec_scl = (optimal_prec - min_prec)/(max_prec - min_prec)
  ) |>
  group_by(species_id, vr) |>
  reframe(
    temp = mean(optimal_temp_scl),
    prec = mean(optimal_prec_scl)
  ) |>
  pivot_longer(
    cols = !c(species_id, vr),
    names_to = 'clim',
    values_to = 'clim_value'
  ) ->
meanOptimalClim_scl


pars_growth |>
  group_by(species_id) |>
  mutate(
    iter = 1:n(),
    temp_growth = tau_temp,
    prec_growth = tau_prec
  ) |>
  ungroup() |>
  select(species_id, iter, temp_growth, prec_growth) |>
  left_join(
    pars_survival |>
      group_by(species_id) |>
      mutate(
        iter = 1:n(),
        temp_survival = tau_temp,
        prec_survival = tau_prec
      ) |>
      ungroup() |>
      select(species_id, iter, temp_survival, prec_survival)
  ) |> 
  left_join(
    pars_recruit |>
      group_by(species_id) |>
      mutate(
        iter = 1:n(),
        temp_recruitment = tau_temp,
        prec_recruitment = tau_prec
      ) |>
      ungroup() |>
      select(species_id, iter, temp_recruitment, prec_recruitment)
  ) |>
  pivot_longer(
    cols = c(contains('temp'), contains('prec')),
    names_to = c('clim', 'vr'),
    names_sep = '_'
  ) |>
  left_join(meanOptimalClim_scl) |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  mutate(
    clim = case_match(clim, 'temp' ~ 'Mean annual temperature',
                            'prec' ~ 'Mean annual precipitation')
  ) |>
  ggplot() +
  aes(x = clim_value, y = value) +
  aes(color = vr) +
  facet_grid(clim~vr) +
  stat_pointinterval(alpha = 0.8) +
  scale_color_manual(
    values = c('#018571', '#dfc27d', '#a6611a')
  ) +
  geom_smooth(method = 'lm') +
  theme_classic() +
  xlab('Mean optimal climate relative to the species range') +
  ylab(expression(tau)) +
  labs(color = '') +
  theme(legend.position = 'top')
```


#### Climate breadth vs observed climate range size

In @fig-climRangeVsTau, we plot the climate range size of a species with its climate breadth.

```{r climate_tauvsRangeSize,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-climRangeVsTau
#| fig-width: 9
#| fig-height: 6
#| fig-cap: "Climate breadth in function of climate range size. The higher the climate range size, the more climate conditions the species experienced."

# function to get range distribution from quantile extremes
get_climateRange <- function(
  var,
  prob_min = 0.01,
  prob_max = 0.99,
  narm = TRUE
){
  min_max = quantile(var, probs = c(prob_min, prob_max), na.rm = narm)
  return( min_max[2] - min_max[1] )
}

treeData |>
  # mutate(
  #   bio_01_mean_nat = unscaleClim(bio_01_mean_scl, 'temp'),
  #   bio_12_mean_nat = unscaleClim(bio_12_mean_scl, 'prec')
  # ) |>
  group_by(species_id) |>
  summarise(
    across(contains('_scl'), get_climateRange)
  ) |>
  rename(
    temp = bio_01_mean_scl,
    prec = bio_12_mean_scl
  ) |>
  pivot_longer(
    cols = c(temp, prec),
    names_to = 'clim',
    values_to = 'clim_val'
  ) ->
sp_range

pars_growth |>
  group_by(species_id) |>
  mutate(
    iter = 1:n(),
    temp_growth = tau_temp,
    prec_growth = tau_prec
  ) |>
  ungroup() |>
  select(species_id, iter, temp_growth, prec_growth) |>
  left_join(
    pars_survival |>
      group_by(species_id) |>
      mutate(
        iter = 1:n(),
        temp_survival = tau_temp,
        prec_survival = tau_prec
      ) |>
      ungroup() |>
      select(species_id, iter, temp_survival, prec_survival)
  ) |> 
  left_join(
    pars_recruit |>
      group_by(species_id) |>
      mutate(
        iter = 1:n(),
        temp_recruitment = tau_temp,
        prec_recruitment = tau_prec
      ) |>
      ungroup() |>
      select(species_id, iter, temp_recruitment, prec_recruitment)
  ) |>
  pivot_longer(
    cols = c(contains('temp'), contains('prec')),
    names_to = c('clim', 'vr'),
    names_sep = '_'
  ) |>
  left_join(sp_range) |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  mutate(
    clim = case_match(clim, 'temp' ~ 'Mean annual temperature',
                            'prec' ~ 'Mean annual precipitation')
  ) |>
  ggplot() +
  aes(clim_val, value) +
  aes(color = vr) +
  facet_grid(clim~vr) +
  # aes(tau_temp, fct_reorder(species_name, tau_temp)) +
  stat_pointinterval(
    alpha = 0.8
  ) +
  scale_color_manual(
    values = c('#018571', '#dfc27d', '#a6611a')
  ) +
  geom_smooth(method = 'lm') +
  theme_classic() +
  xlab('Climate range') +
  ylab(expression(tau)) +
  labs(color = '') +
  theme(legend.position = 'top')
```


## Conditional effect at the lower and upper climate range

To have an integral understanding of the climate effect across the range of the species, we computed the effect size of each climate variable on each vital rate for both the minimum and maximum observed climate conditions.
Effect size ranges from zero (no effect) to 1 (maximum negative effect).
We defined each species' minimum and maximum climate conditions separately using the 0.1 and 99% quantile probability distributions.
In @fig-climEffectSize, we changed the signal of the lower minimum condition to negative to better visualize both lower and upper ranges but consider the absolute value.

```{r climate_prepareRangePerform,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}

treeData |>
  group_by(species_id) |>
  summarise(
    temp_max = quantile(bio_01_mean_scl, prob = 0.99, na.rm = TRUE),
    temp_min = quantile(bio_01_mean_scl, prob = 0.01, na.rm = TRUE),
    prec_max = quantile(bio_12_mean_scl, prob = 0.99, na.rm = TRUE),
    prec_min = quantile(bio_12_mean_scl, prob = 0.01, na.rm = TRUE)
  ) ->
sp_range

# climate effect:
# 0 means no effect of climate
# 1 means that vital rate was reduced to zero at this climate
clim_effect <- function(x_pos, optimal, tau)
  return( 1 - exp(-tau * (x_pos - optimal)^2) )

pars_growth |>
  group_by(species_id) |>
  mutate(iter = row_number()) |>
  ungroup() |>
  left_join(sp_range) |>
  mutate(
    across(contains('temp_'), ~ clim_effect(.x, optimal_temp, tau_temp)),
    across(contains('prec_'), ~ clim_effect(.x, optimal_prec, tau_prec))
  ) |>
  select(species_id, contains('temp_'), contains('prec_')) ->
growth_borderEffect

pars_survival |>
  group_by(species_id) |>
  mutate(iter = row_number()) |>
  ungroup() |>
  left_join(sp_range) |>
  mutate(
    across(contains('temp_'), ~ clim_effect(.x, optimal_temp, tau_temp)),
    across(contains('prec_'), ~ clim_effect(.x, optimal_prec, tau_prec))
  ) |>
  select(species_id, contains('temp_'), contains('prec_')) ->
survival_borderEffect

pars_recruit |>
  group_by(species_id) |>
  mutate(iter = row_number()) |>
  ungroup() |>
  left_join(sp_range) |>
  mutate(
    across(contains('temp_'), ~ clim_effect(.x, optimal_temp, tau_temp)),
    across(contains('prec_'), ~ clim_effect(.x, optimal_prec, tau_prec))
  ) |>
  select(species_id, contains('temp_'), contains('prec_')) ->
recruitment_borderEffect
```

```{r climate_figBorderGrowthMortality,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-climEffectSize
#| fig-width: 9
#| fig-height: 11
#| fig-cap: "Mean annual temperature and mean annual precipitation effect size on growth, recruitment, and survival vital rates between the minimum and maximum observed climate condition. Note that the effect size ranges from 0 (no effect) to 1 (maximum negative effect), and the effect at the lower climate condition was converted to negative for visual clarity."

# fig overall performance across climate and vital rates
growth_borderEffect |>
  pivot_longer(
    cols = !species_id,
    names_to = c('clim', 'range'),
    names_sep = '_'
  ) |> 
  bind_cols(vr = 'growth') |>
  bind_rows(
    survival_borderEffect |>
      pivot_longer(
        cols = !species_id,
        names_to = c('clim', 'range'),
        names_sep = '_'
      ) |> 
      bind_cols(vr = 'survival')
  ) |>
  bind_rows(
    recruitment_borderEffect |>
      pivot_longer(
        cols = !species_id,
        names_to = c('clim', 'range'),
        names_sep = '_'
      ) |> 
      bind_cols(vr = 'recruitment')
  ) |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  mutate(
    value = ifelse(range == 'min', -value, value),
    range = case_match(
      range,
      'max' ~ 'Upper',
      'min' ~ 'Lower'
    ),
    clim = case_match(
      clim,
      'temp' ~ 'Temperature',
      'prec' ~ 'Precipitation'
    )
  ) |>
  ggplot() +
  aes(value, species_name) +
  aes(color = as.factor(range)) +
  facet_grid(clim~vr) +
  stat_pointinterval() +
  theme_classic() +
  scale_color_manual(
    values = c('#fc8d59', '#91bfdb')
  ) +
  geom_vline(xintercept = 0, alpha = 0.8) +
  theme(
    axis.text.y = element_text(face = "italic"),
    legend.position = 'top'
  ) +
  xlab('Marginal efect at the lower (left) and upper (rigth) climate range') +
  ylab('') +
  labs(color = 'Climate\nrange position')
```

```{r climate_figBorderGrowthMortalitySumm,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-climEffectSizeSunn
#| fig-width: 8
#| fig-height: 5
#| fig-cap: "Summary between the species of the climate effect size for each vital rate between the minimum and maximum observed climate condition."

growth_borderEffect |>
  pivot_longer(
    cols = !species_id,
    names_to = c('clim', 'range'),
    names_sep = '_'
  ) |> 
  bind_cols(vr = 'growth') |>
  bind_rows(
    survival_borderEffect |>
      pivot_longer(
        cols = !species_id,
        names_to = c('clim', 'range'),
        names_sep = '_'
      ) |> 
      bind_cols(vr = 'survival')
  ) |>
  bind_rows(
    recruitment_borderEffect |>
      pivot_longer(
        cols = !species_id,
        names_to = c('clim', 'range'),
        names_sep = '_'
      ) |> 
      bind_cols(vr = 'recruitment')
  ) |>
  mutate(
    clim = case_match(clim, 'temp' ~ 'Mean annual temperature', 'prec' ~ 'Mean annual precipitation'),
    range = case_match(range, 'max' ~ 'Upper limit', 'min' ~ 'Lower limit')
  ) |>
  group_by(species_id, vr, clim, range) |>
  reframe(value = mean(value)) |>
  ggplot() +
  aes(range, -value) +
  aes(fill = vr) +
  facet_wrap(~clim) +
  geom_boxplot() +
  scale_fill_manual(
    values = c('#018571', '#dfc27d', '#a6611a')
  ) +
  theme_classic() +
  geom_hline(yintercept = 0, alpha = 0.8) +
  xlab('') +
  ylab('Effect size at climate ranges') +
  labs(fill = '')
```

<!-- 
In @fig-latVsEffect, we compare the climate effect at the lower observed range with the average latitude distribution among all observations of the species.
We have tested multiple combinations of lower/upper and temperature/precipitation. This was the one where there was a clear pattern: species more distributed towards the poles are the ones whose growth rate is more affected by temperature at the lower border. -->

```{r climate_effctSizeRangeSize2,echo=Echo,eval=FALSE,cache=Cache,warning=Warng,message=Msg}
#| label: fig-latVsEffect
#| fig-width: 9
#| fig-height: 5
#| fig-cap: "Effect size of annual mean temperature at lower observed range in function of the median latitudinal position of species."

survival_borderEffect |>
  pivot_longer(
    cols = !species_id,
    names_to = c('clim', 'range'),
    names_sep = '_'
  ) |>
  filter(clim == 'temp' & range == 'min') |>
  group_by(species_id) |>
  reframe(value = mean(value)) |>
  left_join(
    treeData |>
      filter(!is.na(latitude)) |>
      group_by(species_id) |>
      reframe(
        temp = median(latitude)
      )
  ) |>
  bind_cols(vr = 'survival') |>
  bind_rows(
    growth_borderEffect |>
      pivot_longer(
        cols = !species_id,
        names_to = c('clim', 'range'),
        names_sep = '_'
      ) |>
      filter(clim == 'temp' & range == 'min') |>
      group_by(species_id) |>
      reframe(value = mean(value)) |>
      left_join(
        treeData |>
          filter(!is.na(latitude)) |>
          group_by(species_id) |>
          reframe(
            temp = median(latitude)
          )
      ) |>
      bind_cols(vr = 'growth')
  ) |>
  bind_rows(
    recruitment_borderEffect |>
      pivot_longer(
        cols = !species_id,
        names_to = c('clim', 'range'),
        names_sep = '_'
      ) |>
      filter(clim == 'temp' & range == 'min') |>
      group_by(species_id) |>
      reframe(value = mean(value)) |>
      left_join(
        treeData |>
          filter(!is.na(latitude)) |>
          group_by(species_id) |>
          reframe(
            temp = median(latitude)
          )
      ) |>
      bind_cols(vr = 'recruitment')
  ) |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) ->
parsGrowthSurv_mean

survival_borderEffect |>
  pivot_longer(
    cols = !species_id,
    names_to = c('clim', 'range'),
    names_sep = '_'
  ) |>
  bind_cols(vr = 'survival') |>
  bind_rows(
    growth_borderEffect |>
      pivot_longer(
        cols = !species_id,
        names_to = c('clim', 'range'),
        names_sep = '_'
      ) |>
      bind_cols(vr = 'growth')
  ) |>
  bind_rows(
    recruitment_borderEffect |>
      pivot_longer(
        cols = !species_id,
        names_to = c('clim', 'range'),
        names_sep = '_'
      ) |>
      bind_cols(vr = 'recruitment')
  ) |>
  filter(clim == 'temp' & range == 'min') |>
  left_join(
    treeData |>
      filter(!is.na(latitude)) |>
      group_by(species_id) |>
      reframe(
        temp = median(latitude)
      )
  ) |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  ggplot() +
  aes(temp, value) +
  facet_wrap(~vr) +
  stat_pointinterval(alpha = 0.6) +
  geom_text_repel(
    data = parsGrowthSurv_mean,
    aes(x = temp, y = value, label = species_name),
    alpha = 0.8,
    size = 2,
    fontface = 'italic'
  ) +
  geom_smooth(method = 'lm') +
  theme_classic() +
  xlab('Median latitude position') +
  ylab('Effect of annual mean temperature at lower range')
```