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

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:

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

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

\(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:

\(P(z', z, N) = s(z, N) * G(z', z)\)

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

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

\(\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:

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

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

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

\(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`

.

```
library(ipmr)
= list(
data_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
<- rnorm(5, 0, 0.3)
g_r_int <- rnorm(5, 0, 0.7)
s_r_int <- rnorm(5, 0, 0.2)
r_s_r_int
<- paste("r_", 1:5, sep = "")
nms
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 = "")
<- c(data_list, g_r_int, s_r_int, r_s_r_int) params
```

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.

```
<- init_ipm(sim_gen = "simple",
dd_ipm 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.

```
<- define_kernel(
dd_ipm 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.

```
<- define_kernel(
dd_ipm 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.

```
<- 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`

.

```
<- lambda(dd_ipm, type_lambda = "all")
time_step_lams <- lambda(dd_ipm, type_lambda = "stochastic", burn_in = 0.15)
stoch_lam
<- colSums(dd_ipm$pop_state$n_size)
pop_sizes
plot(pop_sizes, type = "l")
```