Many species exhibit age-dependent demography in addition to some
other continuous measures of size. Long term, age classified data sets
aren’t nearly as common as non-age classified data sets, but can be
exceptionally insightful when available. `ipmr`

includes some
methods to deal with specifying them.

While we may be tempted to think of age as a continuous quantity, it must be a discrete state in a discrete time model. This makes the age \(\times\) size IPM a special case of the general IPM. This also uses the parameter index syntax, and it is largely the same as for other, non-age structured IPMs, but includes a couple extra bits that must be specified.

The rest of this vignette assumes you are familiar with the suffix
notation that `ipmr`

uses. If you haven’t already covered
this, it would be best to read at least the Intro
and General
IPM vignettes before continuing. We’ll model Soay sheep (*Ovis
aries*) from St. Kilda using parameters from Chapter 6 of Ellner,
Childs, & Rees (2016).

The deterministic form of an age \(\times\) size model is generally:

- \(n_0(z',t+1) = \sum\limits_{a=0}^M \int_L^UF_a(z',z)n_a(z,t)dz\)

and

- \(n_a(z',t+1) = \int_L^U P_{a-1}(z',z)n_{a-1}(z,t)dz\) for \(a = 1,2,...,M\)

\(M\) indicates the maximum age beyond which we assume no individuals can survive. If we aren’t comfortable with that assumption, then we can add a third equation to specify the number of individuals in the “\(M+1\) or older” age class:

- \(n_{M+1}(z',t+1) = \int_L^U[P_M(z',z)n_M(z,t) + P_{M+1}(z',z)n_{M+1}(z,t)]dz\)

If we do this, then the upper limit in the sum of Eq 1 must be modified to include \(M+1\). The Soay sheep model we are about to work on includes this \(M+1\) age group.

The sub-kernels in this model are formed using regression models for:

survival (

`s_age`

/\(s(z,a)\))growth (

`g_age`

/\(G(z',z,a)\))probability of adults reproducing (

`pb_age`

/\(p_b(z,a)\))probability that a recruit survives to the first census (

`pr_age`

/\(p_r(a)\))recruit size distribution (

`rcsz`

/\(C_d(z',z)\)). The recruit size distribution has a maternal effect of parental weight on recruit weight.

The sub-kernels have the following form:

\(P_a(z', z, a) = s(z, a) * G(z', z, a)\)

\(F_a(z', z, a) = s(z, a) * p_b(z, a) * p_r(a) * C_d(z', z) * 0.5\)

\(F_0 \equiv 0\) so that age-0 recruits cannot reproduce.

This model only tracks females. Assuming an equal sex ratio, we multiply the fecundity kernels by 0.5. We could change the weighting based on observed data.

The vital rates are as follows:

Survival (

`s_age`

/\(s(z,a)\)): A logistic regression with size and age as fixed effects.Example code:

`glm(surv ~ size + age, data = survival_data, family = binomial())`

Mathematical form: \(Logit(s(z,a)) = \alpha_s + \beta_s^z * z + \beta_s^a *age\)

Growth (

`g_age`

/\(G(z',z,a)\)): A linear model with size and age as fixed effects. \(f_G\) denotes a normal probability density function.Example code:

`lm(size_next ~ size + age, data = growth data)`

Mathematical form:

\(G(z',z,a) = f_G(\mu_G(z, a), \sigma_G)\)

\(\mu_G(z, a) = \alpha_G + \beta_G^z * z + \beta_G^a * age\)

Probability of reproduction (

`pb_age`

/\(p_b(z,a)\)): A logistic regression with size and age as fixed effects.Example code:

`glm(repr ~ size + age, data = repr_data, family = binomial())`

Mathematical form: \(Logit(p_b(z,a)) = \alpha_{p_b} + \beta_{p_b}^z * z + \beta_{p_b}^a * age)\)

Probability of recruitment (

`pr_age`

/\(p_r(a)\)): A logistic regression with age as a fixed effect.Example code:

`glm(recr ~ age, data = recr_data, family = binomial())`

Mathematical form: \(Logit(p_r(a)) = \alpha_{p_r} + \beta_{p_r}^a * age\)

Recruit size distribution (

`rcsz`

/\(C_d(z',z)\)): A linear model with parent size as a fixed effect. \(f_{C_d}\) denotes a normal probability density function.Example code:

`lm(rcsz ~ size, data = rcsz_data)`

Mathematical form:

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

\(\mu_{C_d}(z) = \alpha_{C_d} + \beta_{C_d}^z * z\)

This example directly reproduces the model found in Ellner, Rees
& Childs (2016), chapter 6.2. The code that implements that version
can be found here.
This example will skip the IBM simulation and model fitting and just
focus on the new syntax (`ipmr`

doesn’t .

First, we set up our parameter list and define a couple functions to
help out. The `f_fun`

is used to wrap `formula`

argument in `"F_age"`

kernel so that we can concisely express
that age-0 individuals do not reproduce. Note that it is not possible to
use an `ifelse()`

statement instead, because `age`

will be a single number, and `ifelse()`

always returns a
value that is the same length as its input (and we want the return value
to be a \(100\times100\) matrix).

```
library(ipmr)
<- list(
param_list surv_int = -17,
surv_z = 6.68,
surv_a = -0.334,
grow_int = 1.27,
grow_z = 0.612,
grow_a = -0.00724,
grow_sd = 0.0787,
repr_int = -7.88,
repr_z = 3.11,
repr_a = -0.078,
recr_int = 1.11,
recr_a = 0.184,
rcsz_int = 0.362,
rcsz_z = 0.709,
rcsz_sd = 0.159
)
<- function(x) {
inv_logit
return( 1 / (1 + exp(-x)) )
}
<- function(age, s_age, pb_age, pr_age, recr) {
f_fun
if(age == 0) return(0)
* pb_age * pr_age * recr * 0.5
s_age
}
```

Next, we begin to initialize our kernels. `init_ipm`

now
has a fifth argument - `uses_age`

. This is a logical and
denotes that we are specifying a model with individuals cross-classified
by both age and size. The `sim_gen`

argument is
`"general`

“, because age-size models are always general IPMs,
and `det_stoch = "det"`

because we are only working on a
deterministic version of this model for now. We’ll append the
`_age`

index to every variable that is age dependent.
Additionally, we can use `age`

as a standalone term - these
will be substituted during the model building as well.

```
<- init_ipm(sim_gen = "general",
age_size_ipm di_dd = "di",
det_stoch = "det",
uses_age = TRUE) %>%
define_kernel(
name = "P_age",
family = "CC",
formula = s_age * g_age * d_wt,
s_age = inv_logit(surv_int + surv_z * wt_1 + surv_a * age),
g_age = dnorm(wt_2, mu_g_age, grow_sd),
mu_g_age = grow_int + grow_z * wt_1 + grow_a * age,
data_list = param_list,
states = list(c("wt")),
uses_par_sets = FALSE,
age_indices = list(age = c(0:20), max_age = 21),
evict_cor = FALSE
)
```

The primary difference in defining this kernel vs any other indexed
model is that we now specify `uses_par_sets = FALSE`

, and
pass a list to `age_indices`

instead of
`par_set_indices`

. `age_indices`

takes a list with
at least one, but possibly two components:

`age`

: This is the age range for the model. It should always start from 0 and be a sequence of integers.Optionally,

`max_age`

: This is used to denote that while individuals may get older in reality, this value will be the highest that we model. In effect, it creates “very old, but not quite dead” group which can survive and remain in the same age class.

Next, we continue defining the models as we did before.

```
<-
age_size_ipm define_kernel(
proto_ipm = age_size_ipm,
name = "F_age",
family = "CC",
formula = f_fun(age, s_age, pb_age, pr_age, rcsz) * d_wt,
s_age = inv_logit(surv_int + surv_z * wt_1 + surv_a * age),
pb_age = inv_logit(repr_int + repr_z * wt_1 + repr_a * age),
pr_age = inv_logit(recr_int + recr_a * age),
rcsz = dnorm(wt_2, rcsz_mu, rcsz_sd),
rcsz_mu = rcsz_int + rcsz_z * wt_1,
data_list = param_list,
states = list(c("wt")),
uses_par_sets = FALSE,
age_indices = list(age = c(0:20), max_age = 21),
evict_cor = FALSE
)
```

The `"F_age"`

kernel has a custom function in the
`formula`

slot that allows us to always set fecundity for
age-0 individuals to 0. We also add `_age`

suffixes and
`age`

terms to the appropriate equations.

Our call to `define_impl()`

will look a little different.
In examples in the other vignettes using the index syntax, we never
added indices to `state_start`

/`state_end`

.
However, we need to append them here. This is because we have to make
sure that our `P_age`

kernels produce `wt_age`

individuals, whereas our `F_age`

kernels must produce
`wt_0`

individuals (i.e. only age-0 lambs). We do this using
the `state_start`

and `state_end`

arguments.

```
<- age_size_ipm %>%
age_size_ipm define_impl(
make_impl_args_list(
kernel_names = c("P_age", "F_age"),
int_rule = rep("midpoint", 2),
state_start = c("wt_age", "wt_age"),
state_end = c("wt_age", "wt_0")
) )
```

The rest of the steps are largely similar to previous examples:

`define_domains()`

is exactly the same.`define_pop_state()`

: we append`_age`

suffix again here, because we want to track individual weights within age groups. This will create 22 copies of the`wt`

domain, 1 for each`age_indices`

.Optionally,

`define_env_state()`

(though not here because it’s a deterministic model)`make_ipm()`

is exactly the same.

```
<- age_size_ipm %>%
age_size_ipm define_domains(
wt = c(1.6, 3.7, 100)
%>%
) define_pop_state(
n_wt_age = runif(100)
%>%
) make_ipm(
usr_funs = list(inv_logit = inv_logit,
f_fun = f_fun),
iterate = TRUE,
iterations = 100
)
<- lambda(age_size_ipm)
lam lam
```

There are a number of other calculations that we can perform using
functions provided by `ipmr`

. `left_ev()`

and
`right_ev`

also work for age \(\times\) size models. We’ll extract them
and plot their values.

```
<- left_ev(age_size_ipm, iterations = 100)
v_a <- right_ev(age_size_ipm, iterations = 100)
w_a
par(mfrow = c(1, 2))
plot(1:100, seq(0, max(unlist(w_a)), length.out = 100), type = "n",
ylab = expression(paste("w"[a],"(z)")),
xlab = "Size bin")
for(i in seq_along(w_a)) {
lines(w_a[[i]])
}
plot(1:100,
seq(0, max(unlist(v_a)), length.out = 100),
type = "n",
ylab = expression(paste("v"[a], "(z)")),
xlab = "Size bin")
for(i in seq_along(v_a)) {
lines(v_a[[i]])
}
```