pars_randomEffect.qmd

# Random effects {#sec-parsRandomEffects}

In this section, I discuss the plot random effects affecting growth, survival, and recruitment intercepts.
Random effects were used to control for unknown variances grouped within the plot level.
For each vital rate and species, the variance among plots was defined by the parameter $\sigma_{plot}$ where random effects were generated with $N(0, \sigma_{plot})$.

```{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)
library(corrr)
library(sf)
```

```{r re_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 parameters
map_dfr(
  spIds$species_id_old,
  ~ readRDS(
    paste0(
      pars_path, '/growth/', models['growth'], '/posterior_pop_', .x, '.RDS'
    )
  ) |>
  filter(par %in% c('r', 'sigma_plot')) |>
  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', 'sigma_plot')) |>
  pivot_wider(names_from = par) |>
  bind_cols(species_id = .x)
) ->
pars_survival

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

pars_growth |>
  select(species_id, sigma_plot) |>
  bind_cols(vr = 'growth') |>
  bind_rows(
    pars_survival |>
      select(species_id, sigma_plot) |>
      bind_cols(vr = 'survival')
  ) |>
  bind_rows(
    pars_rec |>
      select(species_id, sigma_plot) |>
      bind_cols(vr = 'recruitment')
  ) ->
pars_all

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

## Variance among plots

```{r re_sigma,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-reSigma
#| fig-width: 8.5
#| fig-height: 6
#| fig-cap: "Posterior distribution of variance in plot random effects across species and vital rates. Note that the parameter is in the scale of the intercept, so it is not comparable between demographic models."

pars_all |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  ggplot() +
  aes(sigma_plot, fct_reorder(species_name, sigma_plot)) +
  aes(fill = vr) +
  ggridges::geom_density_ridges2(color = NA, alpha = 0.7) +
  theme_classic() +
  xlab(expression(sigma[plot])) +
  scale_fill_manual(
    values = c('#5ab4ac', '#8856a7', '#d8b365')
  ) +
  ylab('') +
  theme(axis.text.y = element_text(face = "italic")) +
  labs(fill = 'Vital rate') +
  theme(legend.position = 'none') ->
p1

pars_all |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  ggplot() +
  aes(vr, sigma_plot) +
  aes(fill = vr) +
  geom_boxplot() +
  theme_classic() +
  scale_fill_manual(
    values = c('#5ab4ac', '#8856a7', '#d8b365')
  ) +
  theme(legend.position = 'none') +
  xlab('') +
  ylab(expression(sigma[plot])) ->
p2

ggarrange(p1, p2, ncol = 2)

```


## Intercept + random effects

To better understand how plot variation affects growth, survival, and recruitment, the following figures show the intercept distribution for each vital rate and the distribution of a randomly generated offset for plots based on their $\sigma_{plot}$ parameters.


#### Growth

```{r re_intGrowth,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-intGrowth
#| fig-width: 7
#| fig-height: 6
#| fig-cap: "Posterior distribution of the growth rate intercept (red) and intercept with plot random effects (blue). Plot random effects were randomly generated using a normal distribution with zero mean and the posterior distribution of the among plot variance."

pars_growth |>
  mutate(
    re = rnorm(n(), mean = 0, sd = sigma_plot),
    `Intercept +\nrandom effects` = r + re,
    Intercept = r
  ) |>
  pivot_longer(cols = c(Intercept, `Intercept +\nrandom effects`)) |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  ggplot() +
  aes(value, fct_reorder(species_name, value)) + 
  aes(color = name) +
  stat_pointinterval(position = position_dodge(width = 0.75)) +
  theme_classic() +
  scale_color_manual(
    values = c('#fc8d59', '#91bfdb')
  ) +
  theme(axis.text.y = element_text(face = "italic")) +
  labs(
    x = expression(Gamma),
    y = '',
    color = ''
  )
```


#### Survival

```{r re_intSurv,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-intSurv
#| fig-width: 7
#| fig-height: 6
#| fig-cap: "Posterior distribution of the survival intercept (red) and the intercept with plot random effects (blue). Plot random effects were randomly generated using a normal distribution with zero mean and the posterior distribution of the among plot variance."

pars_survival |>
  mutate(
    re = rnorm(n(), mean = 0, sd = sigma_plot),
    `Intercept +\nrandom effects` = psi + re,
    Intercept = psi
  ) |>
  pivot_longer(cols = c(Intercept, `Intercept +\nrandom effects`)) |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  ggplot() +
  aes(value, fct_reorder(species_name, value)) + 
  aes(color = name) +
  stat_pointinterval(position = position_dodge(width = 0.75)) +
  theme_classic() +
  scale_color_manual(
    values = c('#fc8d59', '#91bfdb')
  ) +
  theme(axis.text.y = element_text(face = "italic")) +
  labs(
    x = expression(psi),
    y = '',
    color = ''
  )
```


#### Recruitment

```{r re_intRec,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
#| label: fig-intRec
#| fig-width: 7
#| fig-height: 6
#| fig-cap: "Posterior distribution of the survival intercept (red) and the intercept with plot random effects (blue). Plot random effects were randomly generated using a normal distribution with zero mean and the posterior distribution of the among plot variance."

pars_rec |>
  mutate(
    re = rnorm(n(), mean = 0, sd = sigma_plot),
    `Intercept +\nrandom effects` = mPop_log + re,
    Intercept = mPop_log
  ) |>
  pivot_longer(cols = c(Intercept, `Intercept +\nrandom effects`)) |>
  left_join(
    spIds,
    by = c('species_id' = 'species_id_old')
  ) |>
  ggplot() +
  aes(value, fct_reorder(species_name, value)) + 
  aes(color = name) +
  stat_pointinterval(position = position_dodge(width = 0.75)) +
  theme_classic() +
  scale_color_manual(
    values = c('#fc8d59', '#91bfdb')
  ) +
  theme(axis.text.y = element_text(face = "italic")) +
  labs(
    x = expression(phi),
    y = '',
    color = ''
  )
```


## Correlation in Plot Random Effects Among Species

In this section, we compare the correlation among species regarding the plot offset of each demographic model.

```{r prepareRandomEffects,echo=Echo,eval=FALSE,cache=Cache,warning=Warng,message=Msg}
#setting eval=FALSE to run only once
# code to extract the average posterior distribution for each species-plot-demographic model

plot_mean_ls <- list()
for(vr in c('growth', 'mort', 'recruit'))
{
  spIds$species_id_old |>
    map_dfr(
      ~readRDS(
        paste0(pars_path, '/', vr, '/', models[vr], '/posterior_plot_', .x, '.RDS')
      ) |>
      group_by(plot_id) |> 
      reframe(value = mean(value)) |>
      bind_cols(species_id = .x)
    ) ->
  plot_mean_ls[[vr]]
}

saveRDS(plot_mean_ls, file.path('data', 'plot_mean_ls.RDS'))
```

```{r loadData,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
plot_mean_ls <- readRDS(file.path('data', 'plot_mean_ls.RDS'))

spIds_new <- spIds |>
  select(species_id_old, species_name) |>
  rename(species_id = species_id_old) #|>
  # mutate(
  #   names = str_split(species_name, ' '),
  #   genus = map_chr(names, \(x) unlist(x)[1]),
  #   sp = map_chr(names, \(x) unlist(x)[2]),
  #   genus = paste0(str_extract(genus, '^.{1}'), '.'),
  #   species_name = paste(genus, sp)
  # ) |>
  # select(!c(names, genus, sp))
```

#### Growth

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

plot_mean_ls[['growth']] |>
  left_join(
    spIds_new
  ) |>
  select(!species_id) |> 
  pivot_wider(
    names_from = species_name,
    values_from = value
  ) |>
  select(!plot_id) |>
  correlate(quiet = TRUE) |>
  shave() |> # Remove the upper triangle
  rplot() +
  scale_color_gradient2(limits = c(-1,1), low = rgb(0.57, 0.75, .86, 1), high = rgb(0.99, 0.55, 0.35, 1)) +
  theme(
    legend.position = 'bottom',
    axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1, face = 'italic'),
    axis.text.y = element_text(face = 'italic')
  )
```


#### Survival

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

plot_mean_ls[['mort']] |>
  left_join(
    spIds_new
  ) |>
  select(!species_id) |> 
  pivot_wider(
    names_from = species_name,
    values_from = value
  ) |>
  select(!plot_id) |>
  correlate(quiet = TRUE) |>
  shave() |> # Remove the upper triangle
  rplot() +
  scale_color_gradient2(limits = c(-1,1), low = rgb(0.57, 0.75, .86, 1), high = rgb(0.99, 0.55, 0.35, 1)) +
  theme(
    legend.position = 'bottom',
    axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1, face = 'italic'),
    axis.text.y = element_text(face = 'italic')
  )
```


#### Recruitment

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

plot_mean_ls[['recruit']] |>
  left_join(
    spIds_new
  ) |>
  select(!species_id) |> 
  pivot_wider(
    names_from = species_name,
    values_from = value
  ) |>
  select(!plot_id) |>
  correlate(quiet = TRUE) |>
  shave() |> # Remove the upper triangle
  rplot() +
  scale_color_gradient2(limits = c(-1,1), low = rgb(0.57, 0.75, .86, 1), high = rgb(0.99, 0.55, 0.35, 1)) +
  theme(
    legend.position = 'bottom',
    axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1, face = 'italic'),
    axis.text.y = element_text(face = 'italic')
  )
```

#### Distribution of the correlation among species par

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

plot_mean_ls |>
  map(
    ~ .x |>
      left_join(
        spIds_new
      ) |>
      select(!species_id) |> 
      pivot_wider(
        names_from = species_name,
        values_from = value
      ) |>
      select(!plot_id) |>
      correlate(quiet = TRUE) |>
      shave() |>
      pivot_longer(cols = !term) |>
      drop_na()
  ) |>
  bind_rows(.id = 'vr') |>
  ggplot() +
  aes(value) +
  aes(fill = vr) +
  geom_density(color = NA, alpha = 0.7) +
  theme_classic() +
  labs(
    x = 'Correlation',
    y = '',
    fill = ''
  ) +
  scale_fill_manual(
    values = c('#018571', '#a6611a', '#dfc27d'),
    labels = c('Growth', 'Survival', 'Recruitment')
  )
```

## Spatial Distribution of Plot Random Effects

Our modeling approach used a simple hierarchical model that did not explicitly account for the plot identity.
However, this section allows us to explore potential spatial structure by visualizing the spatial distribution of the plot's random effects.
As illustrated in the figure below, we categorized the continuous offset values into nine distinct bin classes to aid visualization. 
To address the issue of overlapping points within the small figures, we use a transparency gradient that allows us to assign more emphasis to plots with higher absolute values.

```{r prepareSpatialPlot,echo=Echo,eval=Eval,cache=Cache,warning=Warng,message=Msg}
readRDS(file.path(data_path, 'treeData.RDS')) |>
  select(species_id, plot_id, longitude, latitude) |>
  bind_rows(
    readRDS(file.path(data_path, 'fec_dt.RDS')) |>
      select(species_id, plot_id, longitude, latitude)
  ) |>
  group_by(plot_id) |>
  reframe(
    longitude = unique(longitude),
    latitude = unique(latitude)
  ) ->
plot_location

plot_mean_ls |>
  bind_rows(.id = 'vr') |>
  pivot_wider(
    names_from = vr,
    values_from = value
  ) |>
  left_join(
    plot_location
  ) |>
  pivot_longer(cols = c(growth, mort, recruit)) |>
  mutate(
    value = cut(value, c(-2.5, -1.5, -.5, -.25, 0, .25, .5, 1.5, 2.5), include.lowest = TRUE)
  ) |>
  pivot_wider() |>
  st_as_sf(
    coords = c('longitude', 'latitude'),
    crs = 4326
  ) |>
  left_join(spIds_new) ->
spatial_re

plot_spatial_rf <- function(sp)
{ 
  # growth
  spatial_re |>
      filter(species_id == sp) |>
      filter(!is.na(growth)) |>
      ggplot() +
      geom_sf(aes(color = growth), size = 0.1) +
      labs(
        subtitle = 'Growth'
      ) +
      theme_minimal() +
      scale_color_manual(values = col_vec) +
      theme(
        legend.position = 'none',
        axis.text.y = element_blank(),
        axis.text.x = element_blank()
      ) ->
  p1
  
  # survival
  spatial_re |>
      filter(species_id == sp) |>
      filter(!is.na(mort)) |>
      ggplot() +
      geom_sf(aes(color = mort), size = 0.1) +
      labs(
        subtitle = 'Survival'
      ) +
      theme_minimal() +
      scale_color_manual(values = col_vec) +
      theme(
        legend.position = 'none',
        axis.text.y = element_blank(),
        axis.text.x = element_blank()
      ) ->
  p2

  # survival
  spatial_re |>
      filter(species_id == sp) |>
      filter(!is.na(recruit)) |>
      ggplot() +
      geom_sf(aes(color = recruit), size = 0.1) +
      labs(
        subtitle = 'Recruitment'
      ) +
      theme_minimal() +
      scale_color_manual(values = col_vec) +
      theme(
        legend.position = 'none',
        axis.text.y = element_blank(),
        axis.text.x = element_blank()
      ) ->
  p3

  print(
      annotate_figure(
        ggarrange(p1, p2, p3, nrow = 1),
        top = text_grob(
          spIds_new$species_name[spIds_new$species_id == sp],
          face = 'italic', size = 10
        )
      )
  )
}

col_vec <- c(
    rgb(0.57, 0.75, .86, 1),
    rgb(0.57, 0.75, .86, 0.8),
    rgb(0.57, 0.75, .86, 0.5),
    rgb(0.57, 0.75, .86, 0.1),
    rgb(0.99, 0.55, 0.35, 0.1),
    rgb(0.99, 0.55, 0.35, 0.5),
    rgb(0.99, 0.55, 0.35, 0.8),
    rgb(0.99, 0.55, 0.35, 1)
  )
```

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

lvs <- levels(cut(rnorm(100, 0, 3), c(-2.5, -1.5, -.5, -.25, 0, .25, .5, 1.5, 2.5), include.lowest = TRUE))
tibble(lvs) |>
  mutate(
    x = row_number(),
    y = 1,
    col = col_vec
  ) |>
  ggplot() +
  aes(x, y, label = lvs, color = col) +
  geom_point(size = 5) +
  geom_text(aes(color = NULL), nudge_y = 0.25, size = 3) +
  ylim(0.7, 1.3) +
  theme_void() +
  theme(legend.position = 'none') +
  scale_color_manual(values = col_vec[c(4:1, 5:9)])
```

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

for(sp in spIds$species_id_old)
  plot_spatial_rf(sp)
```