What kinds of attributes might we consider alongside a typical attribute such as 'size'?

• size and age (e.g. D. Z. Childs et al. 2003; Moore et al. 2016)
• size and maturation status (e.g. D. Z. Childs et al. 2011)
• size dormant stages (e.g. Rees et al. 2006; Metcalf et al. 2008; Jongejans et al. 2011; Eager et al. 2013)
• size and 'quality' (e.g. Metcalf et al. 2009)

Of course, we don't have to only work with 'size'.

• age and lay date (Dylan Z Childs, Sheldon, and Rees 2016) (not really an IPM…)
• stage and stage duration (Valpine 2009)
• infection status and area of infected tissue (Bruno et al. 2011)
• infection status and parasite number (C. J. E. Metcalf et al. 2016; Wilber et al. 2016)

The space of individual states \(\textbf{Z}\) can include a set of discrete points \(\mathbf{D} = \{ z_1 , \cdots z_D \}\) and a set of continuous domains \(\mathbf{C} = \{ \textbf{Z} _{D + 1} ,\textbf{Z} _{D + 2} , \cdots \textbf{Z} _{D + C} \}\).

• Each continuous domain is either a closed interval \([L_j,U_j]\) or a closed finite rectangle in \(d\)-dimensional space.

• To simplify life, each set in \(\mathbf{D}\) or \(\mathbf{C}\) will be called a component and denoted \(\textbf{Z}_j ,j = 1,2, \cdots , M = D + C\).

For example, \(\textbf{Z}_1\) and \(\textbf{Z}_2\) might represent different sexes, different colour morphs, maturation states (juvenile vs adult), or environmental states, within which individuals are classified by size.

The population state \(n(z,t)\) also consists of a set of \(M\) components:

1. Discrete values \(n_j(t), j = 1, 2, · · · , D\) giving the total number of individuals in each discrete component.

2. Continuous functions \(n_j(z_j , t), j = D+1, D+2,···, M\) giving the state distribution on the continuous components.

The total number of individuals in continuous component \(j\) is given by:

\[ \int_{Z_j} n_j (z_j , t) dz_j \]

Transitions within and among components are described by a set of kernel components \(K_{ij} , 1 ≤ i, j ≤ M\). Four kinds of transition:

1. Continuous to Continuous: \(K_{ij} = K_{ij}(z′,z)\) is a bivariate kernel, just like the ones you've already seen. These are typically fecundity or survival-growth kernels.

2. Discrete to Discrete: \(K_{ij}\) is just a number, such as the proportion of seeds in the seed bank that fail to survive and fail to germinate.

3. Discrete to Continuous: \(K_{ij} = k_{ij}(z′)\) is a state distribution, such as the size distribution of seedlings derived from germinating seeds.

4. Continuous to Discrete: \(K_{ij} = k_{ij}(z)\), is a function that describes the state-dependent contribution individuals to a discrete state, such as the number of seeds produced by an adult plant of size \(z\).

• Many parameters to build one of these…

• IPM version requires someting like 10-20 parameters

The general form of a deterministic age-size IPM can be written:

\[ \begin{eqnarray} n_0(z',t + 1) &=& \sum\limits_{a = 0}^M {\int\limits_{Z} F_a(z',z)n_a(z,t) \, dz} \\ n_a(z',t + 1) &=& \int\limits_{Z} P_{a - 1}(z',z)n_{a - 1}(z,t) dz \quad a = 1,2, \ldots ,M \end{eqnarray} \] This assumes a maximum age beyond which no individual lives. If we want to include a maximum 'greybeard' class this adds one more equation to the model:

\[ n_{M+1}(z',t + 1) = \int\limits_Z{\left[P_{M}(z',z)n_{M}(z,t) + P_{M+1}(z',z)n_{M+1}(z,t)\right] \, dz} \] It's often a good idea to include this extra class because organisms don't really have a maximum age, but data from older individuals is limited.

Next, is a bit of simple house keeping. We need to store the model coefficients in a vector.

m_par <- 
  c(
    surv    = coef(mod_surv),
    grow    = coef(mod_grow),
    grow_sd = summary(mod_grow)$sigma,
    repr    = coef(mod_repr),
    recr    = coef(mod_recr),
    rcsz    = coef(mod_rcsz),
    rcsz_sd = summary(mod_rcsz)$sigma
  )

# clean up the names -- make the intercept labels shorter
names(m_par) <- 
  sub(pattern = "(Intercept)", replace = "int", 
      x = names(m_par), fixed = TRUE)
# clean up the names -- replace . with underscores
names(m_par) <- 
  sub(pattern = ".", replace = "_", 
      x = names(m_par), fixed = TRUE)

Over to you. You need to define five functions:

# Growth kernel
g_z1z <- function(z1, z, a, m_par){
  <define me> 
}
# Survival function, logistic regression
s_z <- function(z, a, m_par){
  <define me> 
}
# Reproduction function, logistic regression
pb_z <- function(z, a, m_par){
  <define me> 
}
# Recruitment function, logistic regression
pr_z <- function(a, m_par) {
  <define me> 
}
# Recruit size kernel
c_z1z <- function(z1, z, a, m_par){
  <define me> 
}

How do we implement an age-size model? Recall the model:

\[ \begin{eqnarray} n_0(z',t + 1) &=& \sum\limits_{a = 0}^M {\int\limits_{Z} F_a(z',z)n_a(z,t) \, dz} \\ n_a(z',t + 1) &=& \int\limits_{Z} P_{a - 1}(z',z)n_{a - 1}(z,t) dz \quad a = 1,2, \ldots ,M \end{eqnarray} \]

The integrals are all one-dimensional, and can be handled in the usual way by the midpoint rule. How do we "put it all together" so that we can iterate the model? There are several options:

1. Embed the iteration matrices for each \(F_a(z',z)\) and \(P_a(z',z)\) into one big, ordinary matrix. This is not efficient because most elements of this matrix will be zero.

2. Embed the iteration matrices for each \(F_a(z',z)\) and \(P_a(z',z)\) into one block-sparse matrix. This is more efficient, but taking care of the indexing can be error prone.

3. Work with a set of age-specific matrices stored in a list. This is more computationally efficient, it's easier and safer to code (and works well for other kinds of general IPM)

We're going to use #3.

# Calculate the 'stuff' we need for integration
mk_intpar <- function(m, L, U, M) {
  <body>
}
# make initial state distribution 
mk_init_nt0 <- function(i_par) {
  <body>
}
# Build the list of age/process specific iteration matrices
mk_age_IPM <- function(i_par, m_par) {
  <body>
}

I'll give you these. Copy them to the end of your script. Then add the next two lines, run everything, and then spend some time looking at i_par, IPM_sys and x:

i_par <- mk_intpar(m = 200, M = 20, L = 1.6, U = 3.7)
IPM_sys <- mk_age_IPM(i_par, m_par)
x <- mk_init_nt0(i_par)

If we iterate an age-size IPM it will converge to a stable age-size distribution, \(w_a(z)\), with the population growing at a constant rate \(\lambda\). This means we can always calculate \(w_a(z)\) and \(\lambda\) by iteration once we know how to implement the model.

What about reproductive value? We need 'the transpose iteration'. With a finite maximum age \(M\), the transpose iteration is \[ v_a(t + 1) = P_a^T v_{a + 1}(t) + F_a^T v_0 (t), \quad a=0,1,2,\cdots,M \] Here \(K_a^T\) represents the transpose kernel \(K_a^T(z',z) = K_a(z,z')\). In practice this is done with the transposes of the iteration matrices for \(P\) and \(F\).

For the greybeard age-class \(M+1\), if the model includes one, \[ v_{M+1}(t + 1) = P_{M+1}^T v_{M + 1}(t) + F_a^T v_0 (t). \]

See if you can calculate lambda by iteration and produce the following plots of \(w_a(z)\) and \(v_a(z)\).

Hint: You need to write code to iterate the model using mk_intpar, mk_age_IPM and mk_init_nt0. You'll need two versions, one that uses r_iter (\(w_a(z)\)) and another that uses l_iter (\(v_a(z)\)).

The long-run, generation-to-generation growth rate \(R_0\) for an age-size IPM is calculated in the same way as for a simple IPM. We add up the expected offspring production in each year of life from a cohort with size distribution \(n_0(z)\) at birth using the 'next generation kernel' \(R\). If individuals do not live beyond age \(M\), then \(R\) is \[ (F_0+F_1P_0+F_2P_1P_0+F_3P_2P_1P_0+ \cdots +F_MP_{M-1} \cdots P_0) = \left(F_0 + \sum\limits_{a=1}^M F_a\prod\limits_{m=0}^{a-1}P_m\right) \] If individuals beyond age \(M\) are lumped into an "\(M+1\) or older" class, then \(R\) includes an extra term: \[ F_{M+1}(I+P_{M+1}+P_{M+1}^2+ \cdots )\prod\limits_{m=0}^{M}P_m = F_{M+1}(I-P_{M+1})^{-1}\prod\limits_{m=0}^{M}P_m. \]

\(R_0\) is the dominant eigenvalue of the \(R\) kernel. How do we do this?

We need a new function to do generation-to-generation iteration:

matmulR <- function(x, IPM_sys, n_iter=100) {
  with(IPM_sys, {
    xtot <- rep(0, length(x))
    for (ia in seq_len(n_iter)) {
      if (ia < na) {
        Pnow <- P[[ia]]
        Fnow <- F[[ia]]
      } else {
        Pnow <- P[[na]]
        Fnow <- F[[na]]
      }
      xtot <- Fnow %*% x + xtot
      x <- Pnow %*% x 
    }
    xtot
  })
}

Write some R code that uses the matmulR function to compute the 'next generation' \(n_0(z)\) and \(R_0\) by iteration.