pars_intercept.qmd

# Intercept {#sec-parsIntercept}

This chapter evaluates the intercept of the growth, survival, and recruitment models.
We use trait information extracted from the literature to compare the average estimation for each demographic model.
The traits of growth classes, maximum observed size, maximum observed age, and shade tolerance are extracted from @burns1990silvics, while the seed mass comes from @diaz2022.

```{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'
)

# Using parameters from the `forest-IPM` repo
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)

map_dfr(
  spIds$species_id_old,
  ~ readRDS(
    paste0(
      pars_path, '/growth/', models['growth'], '/posterior_pop_', .x, '.RDS'
    )
  ) |>
  filter(par %in% c('r', 'Lmax')) |>
  pivot_wider(names_from = par) |>
  bind_cols(species_id = .x)
) ->
pars_growth

map_dfr(
  spIds$species_id_old,
  ~ readRDS(
    paste0(
      pars_path, '/mort/', models['mort'], '/posterior_pop_', .x, '.RDS'
    )
  ) |>
  filter(par %in% c('psi')) |>
  pivot_wider(names_from = par) |>
  bind_cols(species_id = .x)
) ->
pars_mort

map_dfr(
  spIds$species_id_old,
  ~ readRDS(
    paste0(
      pars_path, '/recruit/', models['recruit'], '/posterior_pop_', .x, '.RDS'
    )
  ) |>
  filter(par %in% c('mPop_log', 'p_log')) |>
  pivot_wider(names_from = par) |>
  bind_cols(species_id = .x)
) ->
pars_rec

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


## Growth rate

Because the growth rate intercept decreases non-linearly with size and is governed by the interaction of two parameters ($\Gamma$ and $\zeta_{\infty}$), we computed the average growth rate over a 10-year interval.
Also, because the size of the individual is integrated into the Von Bertalanffy growth model, we computed the 10-year growth average starting from the lower size threshold of 12.7 cm, where growth is optimal.

```{r intcerpt_growth,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-intGrowth
#| fig-width: 8.5
#| fig-height: 5.5
#| fig-cap: "Intercept of the growth model using the 10-year average growth rate. Species are classified by their general growth trait following @burns1990silvics"

pars_growth |>
  mutate(
    size_t0 = 127,
    deltaTime = 10,
    rPlotTime = exp(-exp(r) * deltaTime),
    size_t1 = size_t0 * rPlotTime + Lmax * (1 - rPlotTime),
    growth = (size_t1 - size_t0)/deltaTime
  ) |>
  # # This is the formula for the intercept (growth at zero close to zero)
  # mutate(
  #   Int = Lmax * exp(-exp(r))
  # ) |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  ggplot() +
  aes(growth, fct_reorder(species_name, growth)) +
  aes(color = growth_rate) +
  stat_pointinterval() +
  ylab('') +
  xlab('10 years average growth rate (mm/year)') +
  scale_color_manual(
    values = c("#87bc45", "#edbf33", "#ea5545")
  ) +
  theme_classic() +
  theme(
    axis.text.y = element_text(face = "italic"),
    legend.position = 'none'
  ) ->
p1

pars_growth |>
  group_by(species_id) |>
  reframe(
    r = mean(r),
    Lmax = mean(Lmax)
  ) |>
  mutate(
    size_t0 = 127,
    deltaTime = 10,
    rPlotTime = exp(-exp(r) * deltaTime),
    size_t1 = size_t0 * rPlotTime + Lmax * (1 - rPlotTime),
    growth = (size_t1 - size_t0)/deltaTime
  ) |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  ggplot() +
  aes(growth_rate, growth) +
  aes(fill = growth_rate) +
  geom_boxplot() +
  scale_fill_manual(
    values = c("#87bc45", "#edbf33", "#ea5545")
    # values = c("#20bc45", "#87bc45", "#edbf33", "#ea5545", "#ba0000")
  ) +
  xlab('') +
  ylab('Average 10-years growth rate (mm/year)') +
  theme_classic() +
  theme(legend.position = 'none') +
  labs(fill = '') ->
p2

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

The following figure compares the maximum observed size in the literature with the asymptotic size ($\zeta_{\infty}$), denoting the size where the growth rate converges to zero.

```{r intcerpt_maxObsLmax,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-efszGrowth2
#| fig-width: 6.5
#| fig-height: 5
#| fig-cap: "Correlation between maximum observed size in literature following @burns1990silvics and estimated parameter $\\zeta_{\\infty}$"

parsGrowth_mean <- pars_growth |>
  group_by(species_id) |>
  reframe(
    Lmax = mean(Lmax),
    r = mean(r)
  ) |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) 

pars_growth |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  ggplot() +
  aes(max_size/10, Lmax/10, group = species_id) +
  # aes(color = shade_sylvics) +
  stat_pointinterval(alpha = 0.4) +
  geom_abline(slope = 1, intercept = 0, alpha = 0.3) +
  geom_text_repel(
    data = parsGrowth_mean,
    aes(x = max_size/10, y = Lmax/10, label = species_name),
    alpha = 0.8,
    size = 2.6,
    fontface = 'italic'
  ) +
  theme_classic() +
  xlab('Maximum observed dbh (cm)') +
  ylab('Maximum predicted dbh (cm)')
```


## Mortality probability

Similar to the growth rate, the survival model has a temporal component in which survival probability ($\psi$) decreases exponentially with time.
So, we computed the 10-year average mortality probability as the intercept for the survival model.

```{r intcerpt_survival,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-efszMort
#| fig-width: 8.5
#| fig-height: 5.5
#| fig-cap: "Intercept of the survival model using a 10-year mortality probability. Species are classified by their shade tolerance trait following @burns1990silvics"

pars_mort |>
  mutate(psi = 1 - 1/(1+exp(-psi))^10) |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  ggplot() +
  aes(psi, fct_reorder(species_name, psi)) +
  aes(color = shade_sylvics) +
  stat_pointinterval() +
  ylab('') +
  xlab('10 years mortality probability') +
  scale_color_manual(
    # values = c("#87bc45", "#edbf33", "#ea5545")
    values = c("#20bc45", "#87bc45", "#edbf33", "#ea5545", "#ba0000")
  ) +
  theme_classic() +
  theme(
    axis.text.y = element_text(face = "italic"),
    legend.position = 'none'
  ) ->
p1

pars_mort |>
  group_by(species_id) |>
  reframe(psi = mean(psi)) |>
  mutate(psi = 1 - (1/(1+exp(-psi))^10)) |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  ggplot() +
  aes(shade_sylvics, psi, fill = shade_sylvics) +
  geom_boxplot() +
  scale_fill_manual(
    # values = c("#87bc45", "#edbf33", "#ea5545")
    values = c("#20bc45", "#87bc45", "#edbf33", "#ea5545", "#ba0000")
  ) +
  xlab('Shade tolerance') +
  ylab('10 years mortality probability') +
  theme_classic() +
  theme(legend.position = 'none') +
  labs(fill = '') ->
p2

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

Alternatively, we can also use the parameter $\psi$ to derive the expected longevity ($L$), which can be determined as $L = e^{\psi}$.
The @fig-efszMort2 correlated the maximum size observed in the literature [@burns1990silvics] with longevity ($L$).

```{r intcerpt_maxAgevsPsi,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-efszMort2
#| fig-width: 7.5
#| fig-height: 5
#| fig-cap: "Correlation between maximum size observed in the literature [@burns1990silvics] and the predicted longevity posterior distribution."

parsMort_mean <- pars_mort |>
  group_by(species_id) |>
  reframe(psi = mean(psi)) |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) 

pars_mort |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  ggplot() +
  aes(max_age, exp(psi), group = species_id) +
  # aes(color = shade_sylvics) +
  stat_pointinterval(alpha = 0.4) +
  geom_abline(slope = 1, intercept = 0, alpha = 0.3) +
  geom_text_repel(
    data = parsMort_mean,
    aes(x = max_age, y = exp(psi), label = species_name),
    alpha = 0.8,
    size = 2.6,
    fontface = 'italic'
  ) +
  theme_classic() +
  xlab('Maximum observed age') +
  ylab('Expected longevity')
```


## Ingrowth rate

The ingrowth model comprises two intercept components: the annual ingrowth rate per square meter ($\phi$) and the annual survival probability $\rho$.
In the context of the survival submodel for the ingrowth rate, we computed the 10-year survival probability of ingrowth individuals just like in the survival model.

```{r intcerpt_ingrowth,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-intcerpt_ingrowth
#| fig-width: 8.5 
#| fig-height: 5.5
#| fig-cap: "Intercept of the ingrwoth model for the number of individuals that ingress the population per year per m$^2$. Species are classified by their successional status following @burns1990silvics."

# Fig 3: m ingrowth in function of species and shade tolerance
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  pars_rec |>
    filter(mPop_log < quantile(mPop_log, 0.975)) |>
    left_join(
      spIds,
      by = c('species_id' = 'species_id_old')
    ) |>
    ggplot() +
    aes(mPop_log, fct_reorder(species_name, exp(mPop_log))) +
    aes(color = succession) +
    stat_pointinterval() +
    ylab('') +
    xlab('Annual ingrowth rate per m² (log)') +
    scale_color_manual(
      values = c('#c2e699', '#fb6a4a', '#cb181d')
      # values = c("#87bc45", "#edbf33", "#ea5545")
      # values = c("#20bc45", "#87bc45", "#edbf33", "#ea5545", "#ba0000")
    ) +
    theme_classic() +
    theme(
      axis.text.y = element_text(face = "italic"),
      legend.position = 'none'
    ) ->
  p1

  pars_rec |>
    filter(mPop_log < quantile(mPop_log, 0.975)) |>
    left_join(
      spIds,
      by = c('species_id' = 'species_id_old')
    ) |>
    filter(!is.na(succession)) |>
    ggplot() +
    aes(log(sm), mPop_log, color = succession) +
    stat_pointinterval() +
    scale_color_manual(
      values = c('#fcae91', '#fb6a4a', '#cb181d')
      # values = c("#87bc45", "#edbf33", "#ea5545")
      # values = c("#20bc45", "#87bc45", "#edbf33", "#ea5545", "#ba0000")
    ) +
    theme_classic() +
    labs(
      x = 'Seed mass (ln g)',
      y = 'Annual ingrowth rate per m² (log)',
      color = ''
    ) +
    theme(legend.position = 'top') ->
  p2

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

```{r intcerpt_ingP,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-intcerpt_ingP
#| fig-width: 8.5
#| fig-height: 5.5
#| fig-cap: "Intercept of the annual survival probability for the ingrowth model. Species are classified by their shade tolerance trait following @burns1990silvics."

pars_rec |>
  mutate(p = exp(-exp(p_log))^10) |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  ggplot() +
  aes(p, fct_reorder(species_name, p)) +
  aes(color = shade_sylvics) +
  stat_pointinterval() +
  ylab('') +
  xlab('10 years survival probability') +
  scale_color_manual(
    # values = c("#87bc45", "#edbf33", "#ea5545")
    values = c("#20bc45", "#87bc45", "#edbf33", "#ea5545", "#ba0000")
  ) +
  theme_classic() +
  theme(
    axis.text.y = element_text(face = "italic"),
    legend.position = 'none'
  ) ->
p1

pars_rec |>
  group_by(species_id) |>
  reframe(
    p_log = mean(p_log)
  ) |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  ggplot() +
  aes(shade_sylvics, exp(-exp(p_log))^10) +
  aes(fill = shade_sylvics) +
  geom_boxplot() +
  scale_fill_manual(
    # values = c("#87bc45", "#edbf33", "#ea5545")
    values = c("#20bc45", "#87bc45", "#edbf33", "#ea5545", "#ba0000")
  ) +
  ylab('10 years survival probability') +
  xlab('') +
  theme_classic() +
  theme(legend.position = 'none') ->
p2

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