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:
and
\(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:
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)
param_list <- 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
)
inv_logit <- function(x) {
return( 1 / (1 + exp(-x)) )
}
f_fun <- function(age, s_age, pb_age, pr_age, recr) {
if(age == 0) return(0)
s_age * pb_age * pr_age * recr * 0.5
}
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.
age_size_ipm <- init_ipm(sim_gen = "general",
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
)
lam <- lambda(age_size_ipm)
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.
v_a <- left_ev(age_size_ipm, iterations = 100)
w_a <- right_ev(age_size_ipm, iterations = 100)
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]])
}