Index Notation in ipmr • ipmr

Overview

This vignette describes how to write formulas using ipmr’s notation for models derived from parameter sets. Parameter sets may be derived from mixed models, or fixed effects that are discrete. The notation is meant to mirror the standard mathematical notation used to represent the models. The examples here will primarily use lme4 to illustrate how to go from vital rate models to functional forms that you can use in ipmr. You are of course free to use whatever model fitting software you like, but the number of different possibilities are far too great to cover in this vignette.

The ipmr notation is meant to make implementing models more concise, as it allows you to avoid re-typing/copy+pasting the same kernel formats and just substituting, for example, plot name or transition year. It works by appending an index suffix to kernels, vital rates, and parameters that can take multiple values (e.g. P becomes P_yr). We then add a named list where the names correspond to the index and the entries correspond to the values it can take. We pass pass the list to par_set_indices (i.e. par_set_indices = list(yr = 1:5)), and set uses_par_sets = TRUE.

NB: because of the way these indices are handled internally, they should always appear at the end of the name they modify. For example, use s_slope_site instead of s_site_slope. Multiple indices are fine as long as they chained together at the end (e.g. s_slope_site_year). Models may still work if formatted with indices in the middle of a name, but they also may not (often with obscure error messages).

Quick start guide

Some quick translations between mathematical notation, R model syntax, and ipmr style notation. In the ipmr column, index suffixes that are appended are highlighted in bold.

Suffix names are not restricted to the examples below, and can be anything you want them to be.

Notation guide

The mathematical notation in this vignette will take the following form:

Coefficients/parameters: $\alpha$ - intercept, $\beta$ - slope, $\mu$ - mean, $\sigma$ standard deviation.

Vital rates: subscripts on the coefficients/parameters. For example, $\alpha_G$ for the intercept of a growth model.

Kernels: bivariate kernels are denoted with capital letters (e.g. a growth kernel $G(z',z)$ ) and univariate functions are denoted with lower case letters (e.g. a survival function $s(z)$ ).

Indices: superscripts on coefficients. For example, $\beta_g^{site}$ for a site specific growth slope.

Probability density functions: these are denoted with $f_p$ , where $_p$ is replaced by the vital rate. For example, $f_G(z',\mu_G, \sigma_G)$ would denote a probability density function for growth.

Math Formula corresponds to the mathematical notation for the model. R Formula corresponds to how one could write the model using R’s model formula syntax (think `lm`, `glm`, `lmer`). ipmr shows an equivalent way to write this model in kernel formula or vital rate expression. Model Number indicates which rows are grouped together, as some models have too many terms to cleanly fit on a line, or have multiple equivalent representations.
Math Formula	R Formula	ipmr	Model Number
$\mu_G^{yr} = \alpha_G + \alpha_G^{yr} * \beta_G * z$	`size_2 ~ size_1 + (1 \| year), family = gaussian()`	mu_g_yr = g_int + g_int_yr + g_slope * z	1
$G(z',z) = f_G(z', \mu_g^{yr}, \sigma_G)$		g = dnorm(z_2, mu_g_yr, sd_g)	1
$Logit(\mu_{s}^{plot,yr}) = \alpha_s + \alpha_{s}^{plot,yr} + \beta * z$	`surv ~ size_1 + (1 \| year) + (1 \| plot), family = binomial()`	s_pl_yr = s_int + s_int_pl_yr + s_slope * z	2
	`surv ~ size_1 + (1 \| year) + (1 \| plot), family = binomial()`	s_pl_yr = s_int + s_int_pl + s_int_yr + s_slope * z	2
		s = inv_logit(s_pl_yr)	2
$log(\mu_{r}^{yr}) = \alpha_{r} + \alpha_r^{yr} + (\beta_r + \beta_r^{yr}) * z$	`fec ~ size_1 + (size_1 \| year), family = poisson()`	r = exp(r_int + r_int_yr + (r_slope + r_slope_yr) * z)	3

Internally, ipmr generates the various levels using expand.grid. This means that every combination of the values in par_set_indices must exist in the data_list, or the model will fail with an error along the lines of (Error in eval_tidy: object 'x_yz' not found). In order to exclude levels that do not exist in your data, you can add a vector to the list in par_set_indices called drop_levels. This should contain the values you wish to exclude as a character vector. For example, if data are collected for sites = c("a", "b", "c"), and years = c(2005:2008), but there is no data from site "a" in 2007, we can use par_set_indices = list(site = c("a", "b", "c"), year = c(2005:2008), drop_levels = c("a_2007")).

When picking values to appear in par_set_indices, it is best to avoid longer words or phrases where words are separated by underscores (e.g. treatment = c("h_removal", "c_addition")), as this may cause problems when values are substituted into expressions. CamelCase is generally safer in these cases (e.g. use treatment = c("hRemoval", "cAddition") instead of treatment = c("h_removal", "c_addition")).

Models with a single indexed variable

These models are pretty common in the IPM literature. Examples include IPMs where a single site has been sampled over many years, or sampled many sites across one transition.

This example will use a hypothetical study that spans 5 transitions and has 6 consecutive censuses. The temporal variation will be denoted with a index/superscript appended to each term that it modifies ( $x^{yr}$ /x_yr). To limit complexity, let’s pretend our exploratory analysis and model selection indicated the best overall models included random year-specifc intercepts for growth ( $g^{yr}$ ) and survival ( $s^{yr}$ ), but not probability of reproduction ( $r_p$ ), recruit production ( $r_r$ ), or the recruit size distribution ( $r_d$ ). We will work with a simple model to start (i.e. 1 continuous state variable, no discrete states). This yields the following IPM ( $z,z'$ represent size/height/weight at time $t$ and $t+1$ , respectively):

$n(z', t + 1) = \int_L^U [P^{yr}(z',z) + F(z',z)]n(z,t)dz$
$P^{yr}(z', z) = s^{yr}(z) * g^{yr}(z',z)$
$F(z',z) = r_p(z) + r_r(z) + r_d(z')$

Note that only the $P$ kernel get the $yr$ subscript appended to them - this is the only one that our model selection process decided had substantial year to year variation. Thus, we can actually write the $F$ kernel the exact same way as we would in a simple IPM when we implement the model in ipmr.

The vital rate models could be written on paper as follows:

Growth
- $\mu_G^{yr}(z) = (\alpha_G + \alpha_G^{yr})+ \beta_G * z$
- $G^{yr}(z',z) = f_G(z', \mu_G^{yr}(z), \sigma_G)$
Survival
- $Logit(s^{yr}(z)) = (\alpha_s + \alpha_s^{yr}) + \beta_s * z$
Probability of reproducing
- $Logit(r_p(z)) = \alpha_{r_p} + \beta_{r_p} * z$
Recruit production
- $log(r_r(z)) = \alpha_{r_r} + \beta_{r_r} * z$
Recruit size distribution
- $r_d(z') = f_{r_d}(z', \mu_{r_d}, \sigma_{r_d})$

and modeled in R as follows:

library(lme4)
library(ipmr)

grow_mod <- lmer(size_2 ~ size_1 + (1 | year), data = grow_data)

surv_mod <- glmer(surv ~ size_1 + (1 | year), data = surv_data, family = binomial)

repr_mod <- glm(repr ~ size_1, data = repr_data, family = binomial)

recr_mod <- glm(recr ~ size_1, data = recr_data, family = poisson)

rcsz_mu  <- mean(rcsz_data$size_2)
rcsz_sd  <- sd(rcsz_data$size_2)

We could then extract the fixed coefficients with the fixef method and the conditional modes of the random effects via the ranef method:

beta_gs     <- fixef(grow_mod)

alpha_g_yrs <- ranef(grow_mod)

beta_ss     <- fixef(surv_mod)

alpha_s_yrs <- ranef(surv_mod)

The non-mixed models can use the coef method to extract a vector of coefficients:

repr_coef   <- coef(repr_mod)

recr_coef   <- coef(recr_mod)

The next step is to transform these values into something ipmr can use. This must always be a named list where the names of the list components are the names used in the vital expressions or kernel formula. Transforming the $\beta$ s is easier - the output from fixef is always a named vector, so we just use as.list() and change the names to whatever we want to use. We then bundle it into one big list with all the fixed coefficients.

g_list <- as.list(beta_gs)
names(g_list) <- c("alpha_g", "beta_g")

s_list <- as.list(beta_ss)
names(s_list) <- c("alpha_s", "beta_s")

r_p_list <- as.list(recr_coef)
names(f_p_list ) <- c("alpha_r_p", "beta_r_p")

r_r_list <- as.list(repr_coef)
names(r_r_list) <- c("alpha_r_r", "beta_r_r")

r_d_list <- list(mu_r_d = rcsz_mu, sigma_r_d = rcsz_sd)

all_fixed_params <- c(g_list, s_list, r_p_list, r_r_list, r_d_list)

For the random effects, we probably need to do some a bit more processing, though this really depends on which modeling framework you use and the format of the outputs. For lme4, it looks like this:

g_alpha_list <- as.list(unlist(alpha_g_yrs))
names(g_alpha_list) <- paste("alpha_g_", 2001:2006, sep = "")

s_alpha_list <- as.list(unlist(alpha_s_yrs))
names(s_alpha_list) <- paste("alpha_s_", 2001:2006, sep = "")

First, we strip away all of lme4-related attributes with unlist(), then convert it to the flat list format that we need. Next, we set the names. Now, we can combine it with our all_fixed_params list and we have all of the parameters we need for the kernel functions:

all_params <- c(all_fixed_params, g_alpha_list, s_alpha_list)

We are now ready to begin implementing the kernels. This model will be a simple, density independent, stochastic kernel-resampled model. The parts where the index syntax is used are highlighted by comments in the code block:

ex_ipm <-init_ipm(sim_gen    = "simple",
                  di_dd      = "di", 
                  det_stoch  = "stoch",
                  kern_param = "kern") %>%
  define_kernel(
    
    # The _yr index is appended as a suffix with an underscore. Notice how the formula
    # from Eq 2 is translated into the formula argument and the kernel name
    
    name             = "P_yr",
    family           = "CC",
    formula          = s_yr * g_yr,
    
    # We also modify each parameter name is indexed as well.
    # Here, I've split out the linear predictor and the inverse logit
    # transformation into separate steps to avoid over cluttering
    
    s_lin_p          = alpha_s + alpha_s_yr + beta_s * z_1,
    s_yr             = 1 / (1 + exp(- s_lin_p)),
    
    # We do the same with growth as we did with survival and the P_yr formula
    
    g_yr             = dnorm(z_2, mu_g_yr, sigma_g),
    mu_g_yr          = alpha_g + alpha_g_yr + beta_g * z_1,
    
    data_list        = all_params,
    states           = list(c("z")),
    
    # We signal that the kernel has parameter sets
    
    uses_par_sets    = TRUE,
    
    # And provide the values that each index can take as a named list. 
    # The name(s) in this list MUST match the index used in the expressions.
    
    par_set_indices = list(yr = 2001:2006),
    evict_cor        = TRUE,
    evict_fun        = truncated_distributions("norm", "g_yr")
  ) %>%
  
  # This kernel doesn't get anything special because there are no
  # time-varying parameter values.
  
  define_kernel(
    name          = "F",
    family        = "CC",
    formula       = r_p * r_r * r_d,
    r_p_lin_p     = alpha_r_p + beta_r_p * z_1,
    r_p           = 1 / (1 + exp( -r_p_lin_p)),
    r_r           = exp(alpha_r_r + beta_r_r * z_1),
    r_d           = dnorm(z_2, mu_r_d, sigma_r_d),
    data_list     = all_params,
    states        = list(c("z")),
    uses_par_sets = FALSE,
    evict_cor     = TRUE,
    evict_fun     = truncated_distributions("norm", "r_d")
  ) %>%
  define_impl(
    make_impl_args_list(
      kernel_names = c("P_yr", "F"),
      int_rule     = rep("midpoint", 2),
      state_start    = rep("z", 2),
      state_end      = rep("z", 2)
    )
  ) %>%
  define_domains(
    z = c(0.1, 10, 50)
  ) %>%
  define_pop_state(
    n_z = runif(50)
  ) %>%
  make_ipm(
    iterate    = TRUE,
    iterations = 100,
    kernel_seq = sample(2001:2006, size = 100, replace = TRUE),
    normalize_pop_size = TRUE
  )

IPMs from discrete parameter sets and continuously varying environments

Sometimes, models will have discretely and continuously varying components as functions of different variables. For example, a study could consider multiple sites (discrete component) and spatial and/or temporal fluctuations in the environment at each site (continuous component). We’ll use a hypothetical study that models demography as a function of environment across multiple sites to illustrate this.

We’ll implement a general IPM of a perennial plant species. Our data will come from 6 sites, "A" - "F", and the vital rate intercepts for the plants will vary across each one. Its seeds can either enter a seed bank or recruit immediately, producing recruits at the next time step. We’ll pretend that data from the recruit size distribution was pooled across sites, and so does not vary with space. The seed bank will be age structured and its parameters will be space- and time-invariant. Seeds in the first year can either survive ( $s_{{sb}_1}$ /s_sb1) and recruit ( $r_{{sb}_1}$ ), or survive and become two year old seeds. Two year old seeds can either survive ( $s_{{sb}_2}$ /s_sb2) and recruit ( $r_{{sb}_2}$ /r_sb2) or die.

The vital rate information fed into our hypothetical regression models are ( $s^{site}$ /s_site), growth ( $g^{site}$ /g_site), probability of reproduction ( $c_{p}^{site}$ /c_p_site), and seed production ( $c_{r}^{site}$ /c_r_site). They are all functions of size (z) and environmental covariates ( $\theta^{t}$ and $\theta^{p}$ , denoted together in the kernel formulas as just $\theta$ ).

The full IPM looks like this:

$$ n(z, t + 1) = \int_L^U [P^{site}(z', z,\theta) + r_i * F^{site}(z', z, \theta)]n(z, t)dz +\\ s_{{sb}_1} * r_{{sb}_1} * sb_1(t) * c_d(z') + \\s_{{sb}_2} * r_{{sb}_2} * sb_2(t) * c_d(z') $$
$sb_1(t + 1) = (1-r_i) * \int_L^U [c_p^{site}(z, \theta) * c_r^{site}(z,\theta)] n(z, t)dz$
$sb_2(t + 1) = s_{{sb}_1} * (1-r_{{sb}_1}) * sb_1(t)$

The sub kernels are:

$P^{site}(z',z, \theta) = s^{site}(z, \theta) * G^{site}(z',z, \theta)$
$F^{site}(z',z, \theta) = c_{{p}}^{site}(z, \theta) * c_{r}^{site}(z, \theta) * c_d(z')$
$\theta^{t} \sim Norm(\mu^{t}, \sigma^{t})$
$\theta^{p} \sim Gamma(\alpha^{p}, \beta^{p})$

And the size/climate dependent vital rate functions and example R code are:

$Logit(s^{site}(z, \theta)) = (\alpha_s + \alpha_{s}^{site}) + \beta_s^{z} * z + \beta_s^{t} * \theta^{t} + \beta_s^{p} * \theta^{p}$
- example R code: glmer(survival ~ size + temp + precip + (1 | site), data = my_surv_data, family = binomial())
$G^{site}(z',z,\theta) = f_G(z', \mu_G^{site}(z, \theta), \sigma_G)$
$\mu_G^{site}(z, \theta) = (\alpha_G + \alpha_G^{site}) + \beta_G^{z} * z + \beta_G^{t} * \theta^{t} + \beta_G^{p} * \theta^{p}$
- example R code: lmer(size_next ~ size + temp + precip + (1 | site), data = my_grow_data)
$Logit(c_{p}^{site}(z, \theta)) = (\alpha_{c_p} + \alpha_{c_p}^{site}) + \beta_{c_p}^{z} * z + \beta_{c_p}^{t} * \theta^{t} + \beta_{c_p}^{p} * \theta^{p}$

example R code: glmer(repro ~ size + temp + precip + (1 | site), data = my_repro_data, family = binomial())

$log(c_{r}^{site}(z, \theta)) = (\alpha_{c_r} + \alpha_{c_r}^{site}) + \beta_{c_r}^{z} * z + \beta_{c_r}^{t} * \theta_{t} + \beta_{c_r}^{p} * \theta_{p}$
- example R code: glmer(flower_n ~ size + temp + precip + (1 | site), data = my_flower_data, family = poisson())
$c_d(z') = f_{c_d}(z', \mu_{c_d}, \sigma_{c_d})$
- example R code: mean(my_recr_data$size_next); sd(my_recr_data$size_next)

Simulating parameters

Since this example doesn’t work with actual data, we need to simulate some parameters for each site. Additionally, we’ll specify the distribution parameters for each $\theta$ and write a function that samples them once.

library(ipmr)

set.seed(51)

# Set up the sites, and the parameters

sites <- LETTERS[1:6]

all_pars <- list(
  r_i       = 0.6,
  r_sb1     = 0.2,
  r_sb2     = 0.3,
  s_sb1     = 0.4,
  s_sb2     = 0.1,
  mu_c_d    = 4,
  sigma_c_d = 0.9,
  sigma_g   = 3, 
  beta_s_z    = 2.5,
  beta_s_prec = 0.3,
  beta_s_temp = -0.5,
  beta_g_z    = 0.99,
  beta_g_prec = 0.01,
  beta_g_temp = 0.1,
  beta_c_p_z    = 0.4,
  beta_c_p_prec = 0.6,
  beta_c_p_temp = -0.3,
  beta_c_r_z    = 0.005,
  beta_c_r_prec = -0.3,
  beta_c_r_temp = 0.9
)

g_alphs <- rnorm(6, mean = 4.5, sd = 1) %>%
  as.list() %>%
  setNames(
    paste('alpha_g_', sites, sep = "")
  )

s_alphs <- rnorm(6, mean = -7, sd = 1.5) %>%
  as.list() %>%
  setNames(
    paste('alpha_s_', sites, sep = "")
  )

c_p_alphs <- rnorm(6, mean = -10, sd = 2)  %>%
  as.list() %>%
  setNames(
    paste('alpha_c_p_', sites, sep = "")
  )

c_r_alphs <- runif(6, 0.01, 0.1) %>%
  as.list() %>%
  setNames(
    paste('alpha_c_r_', sites, sep = "")
  )


all_pars <- c(all_pars,
              g_alphs, 
              s_alphs, 
              c_p_alphs,
              c_r_alphs)

# This list contains parameters that define the distributions for environmental
# covariates

env_params <- list(
  temp_mu = 12,
  temp_sd = 2,
  prec_shape = 1000,
  prec_rate  = 2
)

# This function samples the environmental covariate distributions and returns
# the values in a named list. We can reference the names in the list in our
# vital rate expressions.

sample_fun <- function(env_params) {
  
  temp <- rnorm(1, env_params$temp_mu, env_params$temp_sd)
  
  prec <- rgamma(1, env_params$prec_shape, env_params$prec_rate)
  
  out  <- list(temp = temp,
               prec = prec)
  
  return(out)
  
}

Implementing the model

We’ve written out our model - now we can implement it. We’ll use the general, density independent, stochastic parameter resampled machinery for this.

ex_ipm <- init_ipm(sim_gen    = "general",
                   di_dd      = "di", 
                   det_stoch  = "stoch",
                   kern_param = "param") %>%
  define_kernel(
    
    name = "P_site",
    family = "CC",
    
    # site is appended so that we don't have to write out s_A, s_B, etc.
    
    formula          = s_site * g_site * d_z,
    
    s_site_lin       = alpha_s_site + 
                       beta_s_z * z_1 + 
                       beta_s_prec * prec + 
                       beta_s_temp * temp,
    s_site           = 1 / (1 + exp(-s_site_lin)),
    mu_g_site        = alpha_g_site +
                       beta_g_z * z_1 + 
                       beta_g_prec * prec + 
                       beta_g_temp * temp,
    g_site           = dnorm(z_2, mu_g_site, sigma_g),
    data_list        = all_pars,
    states           = list(c("z")),
    uses_par_sets    = TRUE,
    par_set_indices = list(site = LETTERS[1:6]),
    evict_cor        = TRUE,
    evict_fun        = truncated_distributions("norm", "g_site")
  ) %>%
  define_kernel(
    name             = "F_site",
    family           = "CC",
    formula          = c_p_site * c_r_site * c_d * d_z,
    c_p_site_lin     = alpha_c_p_site + 
                       beta_c_p_z * z_1 + 
                       beta_c_p_prec * prec + 
                       beta_c_p_temp * temp,
    c_p_site         = 1 / (1 + exp(-c_p_site_lin)),
    c_r_site         = exp(alpha_c_r_site + 
                           beta_c_r_z * z_1 +
                           beta_c_r_prec * prec +
                           beta_c_r_temp * temp),
    c_d              = dnorm(z_2, mu_c_d, sigma_c_d),
    data_list        = all_pars,
    states           = list(c("z", "sb1", "sb2")),
    uses_par_sets    = TRUE,
    par_set_indices = list(site = LETTERS[1:6]),
    evict_cor        = TRUE,
    evict_fun        = truncated_distributions("norm", "c_d")
    ) %>%
  define_kernel(
    name             = "go_sb_1_site",
    family           = "CD",
    formula          = (1 - r_i) * c_p_site * c_r_site * d_z,
    c_p_site_lin     = alpha_c_p_site + 
                       beta_c_p_z * z_1 + 
                       beta_c_p_prec * prec +
                       beta_c_p_temp * temp,
    c_p_site         = 1 / (1 + exp(-c_p_site_lin)),
    c_r_site         = exp(alpha_c_r_site + 
                           beta_c_r_z * z_1 + 
                           beta_c_r_prec * prec + 
                           beta_c_r_temp * temp),
    data_list        = all_pars,
    states           = list(c("z", "sb1")),
    uses_par_sets    = TRUE,
    par_set_indices = list(site = LETTERS[1:6]),
    evict_cor        = FALSE
  ) %>%
  define_kernel(
    name          = "go_sb_2",
    family        = "DD",
    formula       = s_sb1 * (1 - r_sb1),
    data_list     = all_pars,
    states        = list(c("sb1", "sb2")),
    uses_par_sets = FALSE,
    evict_cor     = FALSE
  ) %>%
  define_kernel(
    name          = "leave_sb_1",
    family        = "DC",
    formula       = s_sb1 * r_sb1 * c_d,
    c_d           = dnorm(z_2, mu_c_d, sigma_c_d),
    data_list     = all_pars,
    states        = list(c("z", "sb1")),
    uses_par_sets = FALSE,
    evict_cor     = TRUE,
    evict_fun     = truncated_distributions("norm", "c_d")
  ) %>%
  define_kernel(
    name          = "leave_sb_2",
    family        = "DC",
    formula       = s_sb2 * r_sb2 * c_d,
    c_d           = dnorm(z_2, mu_c_d, sigma_c_d),
    data_list     = all_pars,
    states        = list(c("z", "sb2")),
    uses_par_sets = FALSE,
    evict_cor     = TRUE,
    evict_fun     = truncated_distributions("norm", "c_d")
  ) %>%
  define_impl(
    make_impl_args_list(
      kernel_names = c("P_site",
                       "F_site",
                       "go_sb_1_site",
                       "go_sb_2", 
                       "leave_sb_1", 
                       "leave_sb_2"),
      int_rule     = rep("midpoint", 6),
      state_start    = c("z", "z", "z", "sb1", "sb1", "sb2"),
      state_end      = c("z", "z", "sb1", "sb2", "z", "z")
    )
  ) %>%
  define_domains(
    z = c(0.1, 100, 200)
  ) %>%
  define_pop_state(
    n_z   = runif(200),
    n_sb1 = 20,
    n_sb2 = 20
  ) %>%
  define_env_state(
    env_state = sample_fun(env_params),
    data_list = list(env_params = env_params,
                     sample_fun = sample_fun)
  ) %>%  
  make_ipm(
    iterations         = 20,
    kernel_seq         = sample(LETTERS[1:6], 20, replace = TRUE),
    normalize_pop_size = TRUE
  )

lambda(ex_ipm)