---
title: "Density Dependent IPMs"
output:
  rmarkdown::html_vignette: 
    toc: yes
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{Density Dependent IPMs}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

# Density dependent models

Density dependent model classes are now implemented. This vignette will get more details shortly. For now, see the example below:

## Example of a simple, stochastic, kernel-resampled model with density dependence

This example assumes that density dependence is modeled as a fixed effect in survival and recruit production models, and assumes there is no density dependence in growth or probability of reproducing models. The survival ($s(z, N)$/`s_yr`), growth ($G_{yr}(z',z)$/ `g_yr`), and number of recruit models ($r_{s,yr}(z, N)$/`r_s_yr`) have year-specific intercepts as well.

The mathematical form for the IPM is below:

1. $n(z', t+1) = K_{yr}(z', z, N)n(z, t)dz$

2. $N = \int_L^Un(z,t)dz$

3. $K_{yr}(z', z, N) = P_{yr}(z', z, N) + F_{yr}(z', z, N)$

Here, $N$ represents the total population size. The kernel values fluctuate as a function of $N$ at each iteration of the model.

The $P_{yr}(z', z, N)$ kernel is comprised of a density independent function for growth (Eq 6-7) and a density dependent function for survival (Eq 5). $f_G$ denotes a Gaussian probability density function: 

4. $P(z', z, N) = s(z, N) * G(z', z)$

5. $Logit(s(z, N)) = \alpha_s + \alpha_{s,yr} + \beta_s^z * z + \beta_s^{N} * N$

6. $G(z', z, \theta) = f_G(z', \mu_{G,yr}(z), \sigma_G)$

7. $\mu_{G,yr}(z) = \alpha_G + \alpha_{G,yr} + \beta_G^z * z$

The $F_{yr}(z',z, N)$ kernel is comprised of a density independent function for recruit size (Eq 10) and probability of reproducing (Eq 9), and a density dependent function for number of recruits produced by parents (Eq 11). $f_{r_d}$ denotes a Gaussian probability density function: 

8. $F_{yr}(z', z, N) = r_r(z) * r_{s,yr}(z, N) + r_d(z')$

9. $Logit(r_r(z)) = \alpha_{r_r} + \beta_{r_r}^z * z$

10. $r_d(z') = f_{r_d}(z', \mu_{r_d}, \sigma_{r_d})$

11. $Log(r_{s,yr}(z, N)) = \alpha_{r_s} + \alpha_{{r_s},yr} + \beta_{r_s}^z * z + \beta_{r_s}^N * N$

We'll simulate a 50 year time series using hypothetical parameter values. The fixed parameter values are created as with a density independent model. The difference is that we now have two more parameters: `s_dd`, and `r_s_dd`. These are the coefficients that correspond to $\beta_s^N$ and $\beta_{r_s}^N$, respectively. The chunk below initializes the data list object, which we name `params`.


```{r eval = FALSE}

library(ipmr)

data_list = list(
  s_int     = 1.03,
  s_slope   = 2.2,
  s_dd      = -0.7,
  g_int     = 8,
  g_slope   = 0.92,
  sd_g      = 0.9,
  r_r_int   = 0.09,
  r_r_slope = 0.05,
  r_s_int   = 0.1,
  r_s_slope = 0.005,
  r_s_dd    = -0.03,
  mu_rd     = 9,
  sd_rd     = 2
)

# Now, simulate some random intercepts for growth, survival, and offspring production

g_r_int   <- rnorm(5, 0, 0.3)
s_r_int   <- rnorm(5, 0, 0.7)
r_s_r_int <- rnorm(5, 0, 0.2)

nms <- paste("r_", 1:5, sep = "")

names(g_r_int) <- paste("g_", nms, sep = "")
names(s_r_int) <- paste("s_", nms, sep = "")
names(r_s_r_int) <- paste("r_s_", nms, sep = "")

params     <- c(data_list, g_r_int, s_r_int, r_s_r_int)

```


Next, we initialize the model using `init_ipm`. The difference is that the second argument is now changed to `"dd"` to denote that this is a density dependent model.

```{r eval = FALSE}

dd_ipm <- init_ipm(sim_gen = "simple",
                   di_dd = "dd", 
                   det_stoch = "stoch", 
                   kern_param = "kern")
```

Once we've done that, we're ready to begin specifying the kernel forms. One previously not mentioned aspect of `define_pop_state()` is that, in addition to defining initial conditions, 2 additional helper variables are generated: `n_stateVariable_t` and `n_stateVariable_t_1`. These can be used to reference the population states in vital rate and/or kernel expressions.

These will look very similar to the ones we specified for density-independent models, except that we now include the term `s_dd * sum(n_size_t)` in the survival expression. `sum(n_size_t)` is the syntax `ipmr` uses to denote total population size. Further down, there is an example of how to use subsets of the trait distribution. 

```{r eval = FALSE}
dd_ipm <- define_kernel(
  proto_ipm        = dd_ipm,
  name             = "P_yr",
  formula          = s_yr * g_yr,
  family           = "CC",
  s_yr             = plogis(s_int + s_r_yr + s_slope * size_1 + s_dd * sum(n_size_t)),
  g_yr             = dnorm(size_2, g_mu_yr, sd_g),
  g_mu_yr          = g_int + g_r_yr + g_slope * size_1,
  data_list        = params,
  states           = list(c("size")),
  uses_par_sets    = TRUE,
  par_set_indices = list(yr = 1:5),
  evict_cor        = TRUE,
  evict_fun        = truncated_distributions("norm", "g_yr")
) 
  
```

Other than the inclusion of the density dependent term in the survival expression, this should look quite similar to the density-independent kernel-resampled models from the Introduction vignette. We are now ready to continue defining the $F_{yr}(z',z,N)$ kernel.


```{r eval = FALSE}

dd_ipm <- define_kernel(
  proto_ipm        = dd_ipm,
  name             = "F_yr",
  formula          = r_r * r_s_yr * r_d,
  family           = "CC",
  r_r              = plogis(r_r_int + r_r_slope * size_1),
  r_s_yr           = exp(r_s_int + r_s_r_yr + r_s_slope * size_1 + r_s_dd * sum(n_size_t)),
  r_d              = dnorm(size_2, mu_rd, sd_rd),
  data_list        = params,
  states           = list(c("size")),
  uses_par_sets    = TRUE,
  par_set_indices = list(yr = 1:5),
  evict_cor        = TRUE,
  evict_fun        = truncated_distributions("norm", "r_d")
  ) 

```

Again, we've add the `f_s_dd * sum(n_size_t)` to the expression for `f_s_yr`, but otherwise, not much is different from how we've defined density independent models. The rest of the model definition process is unchanged.

```{r eval = FALSE}
 dd_ipm <-  dd_ipm %>%
  define_impl(
    make_impl_args_list(
      kernel_names = c("P_yr", "F_yr"),
      int_rule     = rep("midpoint", 2),
      state_start    = rep("size", 2),
      state_end      = rep("size", 2)
    )
  ) %>%
  define_domains(
    size = c(0, 50, 200)
  ) %>%
  define_pop_state(
    n_size = runif(200)
  ) %>%
  make_ipm(
    iterate = TRUE,
    iterations = 50,
    kernel_seq = sample(1:5, 50, replace = TRUE)
  )


```

`lambda` methods are defined for all density-dependent models as well. It is fairly straightforward to plot population sizes for these models by extracting the column sums of the arrays in `pop_state`.  

```{r eval = FALSE}
time_step_lams <- lambda(dd_ipm, type_lambda = "all")
stoch_lam      <- lambda(dd_ipm, type_lambda = "stochastic", burn_in = 0.15)

pop_sizes <- colSums(dd_ipm$pop_state$n_size)

plot(pop_sizes, type = "l")

```