Age-Size IPMs

Introduction to age \(\times\) size classified models

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).

Mathematical overview

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

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

and

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

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

The sub-kernels have the following form:

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

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

  1. 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\)

  2. 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\)

  3. 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)\)

  4. 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\)

  5. 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\)

Model parameterization

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:

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

  2. 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:

  1. define_domains() is exactly the same.

  2. 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.

  3. Optionally, define_env_state() (though not here because it’s a deterministic model)

  4. 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

Further analyses

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]]) 
  
}