Using your own models

One of the design principles of cxr is to allow users to estimate parameters from their own models. cxr already includes four families of well-known population dynamic models: Beverton-Holt (BH), Lotka-Volterra (LV), Ricker (RK), and Law-Watkinson (LW), all in their discrete time formulation (see the accompanying publication of cxr for formulation of these general families). Users can choose between these default families by setting the model_family argument to the acronym of the chosen model. In this vignette we show, step by step, how other user-defined models can be integrated in the package.

First of all, a common feature of the families included is that they rely on a common set of parameters, namely lambda,alpha (separated in alpha_intra and alpha_inter), lambda_cov, and alpha_cov. For now, user-defined models within cxr must be restricted to these parameters. A detailed description of these parameters can be found in previous vignettes. This does not preclude users to hard-code other model parameters inside the model functions, and future releases will aim to relax this assumption.

At the end of this vignette we include the full R code of a model template that can be used directly into cxr. This template is also included as a stand-alone R file in the source files of the package (cxr_model_template.R), and here we first explain each section of it.

First, model functions should be defined with the common set of arguments recognized by cxr.

family_pm_alpha_form_lambdacov_form_alphacov_form <- function(par,
                                                              fitness,
                                                              neigh_intra_matrix = NULL,
                                                              neigh_inter_matrix,
                                                              covariates,
                                                              fixed_parameters)

Assuming data for n observations of a focal species and s neighbour species (including itself), the function arguments are:

par: one-dimensional numeric vector including all parameters to fit.
fitness: one-dimensional vector of length n with fitness observations in log scale.
neigh_intra_matrix: matrix of n rows and 1 column, with observations of intraspecific neighbours.
neigh_inter_matrix: matrix of n rows and s − 1 columns, with observations of interspecific neighbours.
covariates: dataframe with n rows and as many columns as covariates.
fixed_parameters: list with values of parameters that are not to be fitted.

the name of the function should be adapted to the parameters you want to fit:

‘family’ is the two-letter acronym of your model family.
‘alpha_form’ is one of ‘alpha_none’,‘alpha_global’, or ‘alpha_pairwise’, depending on whether no alphas, a single alpha, or pairwise alphas are included in your model.
‘lambda_cov_form’ is one of ‘lambda_cov_none’ or ‘lambda_cov_global’, depending on whether you include the effect of covariates over lambda.
‘alpha_cov_form’ is one of ‘alpha_cov_none’, ‘alpha_cov_global’, or ‘alpha_cov_pairwise’, depending on whether and how to include the effect of covariates over alpha values: no inclusion, a global parameter for each covariate on the alpha values, or interaction-specific parameters for each covariate.

Thus, if you code a model from your User Model (UM) family that includes pairwise alphas, but no effects of covariates over lambda or alpha, your function must be named as:

UM_pm_alpha_pairwise_lambdacov_none_alphacov_none <- function(par,
                                                              fitness,
                                                              neigh_intra_matrix = NULL,
                                                              neigh_inter_matrix,
                                                              covariates,
                                                              fixed_parameters)

The first step inside the function is to retrieve the different parameters from the one-dimensional par vector. You should not have to change this section except for commenting or commenting out the section where parameters are retrieved. For example, if your model does not include lambda_cov or alpha_cov, you should comment their sections. Note that it includes a final parameter, sigma, that should always be there, as it keeps track of the error associated to the fit.

This is how this part in your ‘UM_pm_alpha_pairwise_lambdacov_none_alphacov_none’ model would look like:

# retrieve parameters -----------------------------------------------------
# parameters to fit are all in the "par" vector,
# so we need to retrieve them one by one
# order is {lambda,lambda_cov,alpha,alpha_cov,sigma}

# comment or uncomment sections for the different parameters
# depending on whether your model includes them
# note that the section on alpha_inter includes two
# possibilities, depending on whether a single alpha is 
# fitted for all interactions (global) or each pairwise alpha is 
# different (pairwise)
# both are commented, you need to uncomment the appropriate one

# likewise for the section on alpha_cov

# --------------------------------------------------------------------------

pos <- 1

# if a parameter is passed within the "par" vector,
# it should be NULL in the "fixed_parameters" list

# lambda
if(is.null(fixed_parameters$lambda)){
  lambda <- par[pos]
  pos <- pos + 1
}else{
  lambda <- fixed_parameters[["lambda"]]
}

# lambda_cov
# if(is.null(fixed_parameters$lambda_cov)){
#   lambda_cov <- par[pos:(pos+ncol(covariates)-1)]
#   pos <- pos + ncol(covariates)
# }else{
#   lambda_cov <- fixed_parameters[["lambda_cov"]]
# }

# alpha_intra
if(!is.null(neigh_intra_matrix)){
  # intra
  if(is.null(fixed_parameters[["alpha_intra"]])){
    alpha_intra <- par[pos]
    pos <- pos + 1
  }else{
    alpha_intra <- fixed_parameters[["alpha_intra"]]
  }
}else{
  alpha_intra <- NULL
}

# alpha_inter
if(is.null(fixed_parameters[["alpha_inter"]])){
  # uncomment for alpha_global
  # alpha_inter <- par[pos]
  # pos <- pos + 1
  
  # uncomment for alpha_pairwise
  alpha_inter <- par[pos:(pos+ncol(neigh_inter_matrix)-1)]
  pos <- pos + ncol(neigh_inter_matrix)
}else{
  alpha_inter <- fixed_parameters[["alpha_inter"]]
}

# alpha_cov
# if(is.null(fixed_parameters$alpha_cov)){
#   # uncomment for alpha_cov_global
#   # alpha_cov <- par[pos:(pos+ncol(covariates)-1)]
#   # pos <- pos + ncol(covariates)
#   
#   # uncomment for alpha_cov_pairwise
#   # alpha_cov <- par[pos:(pos+(ncol(covariates)*
#   # (ncol(neigh_inter_matrix)+ncol(neigh_intra_matrix)))-1)]
#   # pos <- pos + (ncol(covariates)*(ncol(neigh_inter_matrix)+ncol(neigh_intra_matrix)))
# }else{
#   alpha_cov <- fixed_parameters[["alpha_cov"]]
# }

# sigma - this is always necessary
sigma <- par[length(par)]

At this point, you can code your model. This is how the equivalent Beverton-Holt model looks like, for reference:

# MODEL CODE HERE ---------------------------------------------------------

# the model should return a "pred" value
# a function of lambda, alpha_intra, alpha_inter, lambda_cov, alpha_cov 
# and neigh_intra_matrix, neigh_inter_matrix, and covariates

# we do not differentiate alpha_intra from alpha_inter in this model
# so, put together alpha_intra and alpha_inter, and the observations
# with intraspecific ones at the beginning
if(!is.null(alpha_intra)){
  alpha <- c(alpha_intra,alpha_inter)
  all_neigh_matrix <- cbind(neigh_intra_matrix,neigh_inter_matrix)
}else{
  alpha <- alpha_inter
  all_neigh_matrix <- neigh_inter_matrix
}

term = 1 #create the denominator term for the model
for(z in 1:ncol(all_neigh_matrix)){
  term <- term + alpha[z]*all_neigh_matrix[,z] 
}
pred <- lambda/ term

# MODEL CODE ENDS HERE ----------------------------------------------------

Note how there are no hard-coded checks for ensuring that alpha values and neigh_matrix values are appropriately sorted. These checks are performed in the cxr_pm_fit function that internally sorts data and call this model for optimization.

The last part of the model function is returning the negative log-likelihood value for that particular parameter combination. This is the value that the numerical optimization algorithms aim to minimize. In general, you should not need to change this code.

# the routine returns the sum of negative log-likelihoods of the data and model:
# DO NOT CHANGE THIS
llik <- dnorm(fitness, mean = (log(pred)), sd = (sigma), log=TRUE)
return(sum(-1*llik))

When you have a functioning model, you should be able to use it within cxr as soon as you load your model into the global environment, which can be done simply by ‘sourcing’ your R file named with the name of your model, i.e. ‘UM_pm_alpha_pairwise_lambdacov_none_alphacov_none.R’ in our example.

Thus, for fitting a certain dataset with your custom model, you would call the cxr_pm_fit function with the appropriate parameters:

# load your model into the global environment
source("./pm_UM_alpha_pairwise_lambdacov_none_alphacov_none.R")
# fit your data
custom_fit <- cxr_pm_fit(data = custom_data, # assuming custom_data is already set...
                         focal_column = my_focal, # assuming my_focal is already set...
                         model_family = "UM",
                         covariates = NULL, # as we have no covariate effects
                         alpha_form = "pairwise",
                         lambda_cov_form = "none",
                         alpha_cov_form = "none")

If the function cannot find a model matching the provided ‘model_family’,‘alpha_form’,‘lambda_cov_form’, and ‘alpha_cov_form’, it will output a message and return NULL. On a final note, if you think your bright new model can be useful for fellow scientists, we encourage you to make a pull request for integrating it in the set of default model families included in cxr. An interesting update to the package would be, in addition to include other model families, to implement non-linear relationships of covariate effects over lambda and alpha values.

There are several issues that you should have in mind when developing custom models. First, recall that numerical optimization methods may be bounded or unbounded. If your model does not return meaningful values for part of the parameter space defined by the hypervolume of parameter boundaries, bounded optimization methods (such as ‘L-BFGS-B’ or ‘bobyqa’) may fail. This is well exemplified by the discrete Lotka-Volterra model:

λ_i − α_iiN_i − α_ijNj

which can easily produce negative values and thus undefined log-likelihoods. The second point to note is that user-defined population dynamics models are sufficient to obtain tailored lambda and alpha values, but importantly, several MCT metrics are also model-specific.

Model-specific coexistence metrics

Among the coexistence metrics that can be calculated with cxr, five of them are model-specific, i.e. have different formulations depending on the underlying population dynamics model from which lambda and alpha values are obtained. These are:

demographic ratio (and therefore, average fitness differences in the MCT formulation)
competitive ability
competitive effects and responses
species fitness in the absence of niche differences

All these have different formulations for different model families. Therefore, if you want to estimate coexistence metrics based on custom models, you should provide a full set of formulations for these metrics. Here we show their formulation for Beverton-Holt models, for reference. The supplementary material of Godoy and Levine (2014) and that of Hart et al. (2018) are good places to explore the mathematical underpinnings of these particularities.

demographic ratio:

BH_demographic_ratio <- function(pair_lambdas){
  (pair_lambdas[1]-1)/(pair_lambdas[2]-1)
}

competitive ability:

BH_competitive_ability <- function(lambda, pair_matrix){
  if(all(pair_matrix >= 0)){
    (lambda - 1)/sqrt(pair_matrix[1,1] * pair_matrix[1,2])
  }else{
    NA_real_
  }
}

competitive effects and responses: these are obtained with an optimization function cxr_er_fit, which is similar to cxr_pm_fit (see vignette 1). We include at the end of this vignette the complete model template for a custom effect and response model. The formulation of the Beverton-Holt effect and response model is simply

pred <- lambda.part/ (1+ e.part*r.part )

where each part potentially depends on neighbour densities and covariates.

species fitness in the absence of niche differences:

BH_species_fitness <- function(lambda, competitive_response){
  (lambda-1)/competitive_response
}

Therefore, in order to compute coexistence metrics based on your custom population dynamics model, you should code your own versions of these metrics, save them with appropriate names (i.e. changing the model family acronym at the beginning of the functions), and source them in the global environment, in order for the cxr interface to be able to find and use them.

Template for population dynamics models

family_pm_alpha_form_lambdacov_form_alphacov_form <- function(par,
                                                              fitness,
                                                              neigh_intra_matrix = NULL,
                                                              neigh_inter_matrix,
                                                              covariates,
                                                              fixed_parameters){
  
  
  # retrieve parameters -----------------------------------------------------
  # parameters to fit are all in the "par" vector,
  # so we need to retrieve them one by one
  # order is {lambda,lambda_cov,alpha_intra,alpha_inter,alpha_cov,sigma}
  
  # comment or uncomment sections for the different parameters
  # depending on whether your model includes them
  # note that the section on alpha_inter includes two
  # possibilities, depending on whether a single alpha is 
  # fitted for all interactions (global) or each pairwise alpha is 
  # different (pairwise)
  # both are commented, you need to uncomment the appropriate one
  
  # likewise for the section on alpha_cov
  
  # --------------------------------------------------------------------------
  
  pos <- 1
  
  # if a parameter is passed within the "par" vector,
  # it should be NULL in the "fixed_parameters" list
  
  # lambda
  if(is.null(fixed_parameters$lambda)){
    lambda <- par[pos]
    pos <- pos + 1
  }else{
    lambda <- fixed_parameters[["lambda"]]
  }
  
  # lambda_cov
  if(is.null(fixed_parameters$lambda_cov)){
    lambda_cov <- par[pos:(pos+ncol(covariates)-1)]
    pos <- pos + ncol(covariates)
  }else{
    lambda_cov <- fixed_parameters[["lambda_cov"]]
  }
  
  # alpha_intra
  if(!is.null(neigh_intra_matrix)){
    # intra
    if(is.null(fixed_parameters[["alpha_intra"]])){
      alpha_intra <- par[pos]
      pos <- pos + 1
    }else{
      alpha_intra <- fixed_parameters[["alpha_intra"]]
    }
  }else{
    alpha_intra <- NULL
  }
  
  # alpha_inter
  if(is.null(fixed_parameters[["alpha_inter"]])){
    # uncomment for alpha_global
    # alpha_inter <- par[pos]
    # pos <- pos + 1
    
    # uncomment for alpha_pairwise
    # alpha_inter <- par[pos:(pos+ncol(neigh_inter_matrix)-1)]
    # pos <- pos + ncol(neigh_inter_matrix)
  }else{
    alpha_inter <- fixed_parameters[["alpha_inter"]]
  }
  
  # alpha_cov
  if(is.null(fixed_parameters$alpha_cov)){
    # uncomment for alpha_cov_global
    # alpha_cov <- par[pos:(pos+ncol(covariates)-1)]
    # pos <- pos + ncol(covariates)
    
    # uncomment for alpha_cov_pairwise
    # alpha_cov <- par[pos:(pos+(ncol(covariates)*
    # (ncol(neigh_inter_matrix)+ncol(neigh_intra_matrix)))-1)]
    # pos <- pos + (ncol(covariates)*(ncol(neigh_inter_matrix)+ncol(neigh_intra_matrix)))
  }else{
    alpha_cov <- fixed_parameters[["alpha_cov"]]
  }
  
  # sigma - this is always necessary
  sigma <- par[length(par)]
  
  # now, parameters have appropriate values (or NULL)
  # next section is where your model is coded
  
  # MODEL CODE HERE ---------------------------------------------------------
  
  # the model should return a "pred" value
  # a function of lambda, alpha_intra, alpha_inter, lambda_cov, alpha_cov 
  # and neigh_intra_matrix, neigh_inter_matrix, and covariates
  pred <- 0
  
  # MODEL CODE ENDS HERE ----------------------------------------------------
  
  # the routine returns the sum of log-likelihoods of the data and model:
  # DO NOT CHANGE THIS
  llik <- dnorm(fitness, mean = (log(pred)), sd = (sigma), log=TRUE)
  return(sum(-1*llik))
}

Template for effect and response models

family_er_lambdacov_form_effectcov_form_responsecov_form <- function(par,
                                                                     fitness,
                                                                     target,
                                                                     density,
                                                                     covariates,
                                                                     fixed_parameters){
  
  num.sp <- nrow(target)
  
  # parameters to fit are all in the "par" vector,
  # so we need to retrieve them one by one
  # order is {lambda,lambda_cov,effect,effect_cov,response,response_cov,sigma}
  
  # comment or uncomment sections for the different parameters
  # depending on whether your model includes them
  # note that effect and response models must always include
  # lambda, effect, and response, at least.
  
  pos <- 1
  
  # if a parameter is passed within the "par" vector,
  # it should be NULL in the "fixed_parameters" list
  
  # lambda
  if(is.null(fixed_parameters[["lambda"]])){
    lambda <- par[pos:(pos + num.sp - 1)]
    pos <- pos + num.sp
  }else{
    lambda <- fixed_parameters[["lambda"]]
  }
  
  # lambda_cov
  if(is.null(fixed_parameters$lambda_cov)){
    # the covariate effects are more efficient in a matrix form
    # with species in rows (hence byrow = T, because by default
    # the vector is sorted first by covariates)
    lambda_cov <- matrix(par[pos:(pos+(ncol(covariates)*num.sp)-1)],
                         nrow = num.sp,
                         byrow = TRUE)
    pos <- pos + ncol(covariates)*num.sp
  }else{
    lambda_cov <- fixed_parameters[["lambda_cov"]]
  }
  
  # effect
  if(is.null(fixed_parameters[["effect"]])){
    effect <- par[pos:(pos + num.sp - 1)]
    pos <- pos + num.sp
  }else{
    effect <- fixed_parameters[["effect"]]
  }
  
  # effect_cov
  if(is.null(fixed_parameters$effect_cov)){
    effect_cov <- matrix(par[pos:(pos+(ncol(covariates)*num.sp)-1)],
                         nrow = num.sp,
                         byrow = TRUE)
    pos <- pos + ncol(covariates)*num.sp
  }else{
    effect_cov <- fixed_parameters[["effect_cov"]]
  }
  
  # response
  if(is.null(fixed_parameters[["response"]])){
    response <- par[pos:(pos + num.sp - 1)]
    pos <- pos + num.sp
  }else{
    response <- fixed_parameters[["response"]]
  }
  
  # response_cov
  if(is.null(fixed_parameters[["response_cov"]])){
    response_cov <- matrix(par[pos:(pos+(ncol(covariates)*num.sp)-1)],
                           nrow = num.sp,
                           byrow = TRUE)
    pos <- pos + ncol(covariates)*num.sp
  }else{
    response_cov <- fixed_parameters[["response_cov"]]
  }
  
  sigma <- par[length(par)]
  
  # now, parameters have appropriate values (or NULL)
  # next section is where your model is coded
  
  # MODEL CODE HERE ---------------------------------------------------------
  
  # the model should return a "pred" value
  # a function of lambda, effect, response, lambda_cov, effect_cov, response_cov
  pred <- 0
  
  # MODEL CODE ENDS HERE ----------------------------------------------------
  
  # the routine returns the sum of log-likelihoods of the data and model:
  # DO NOT CHANGE THIS
  llik <- dnorm(fitness, mean = (log(pred)), sd = (sigma), log=TRUE)
  return(sum(-1*llik))
  
}

References

Godoy, O., & Levine, J. M. (2014). Phenology effects on invasion success: insights from coupling field experiments to coexistence theory. Ecology, 95(3), 726-736.

Hart, S. P., Freckleton, R. P., & Levine, J. M. (2018). How to quantify competitive ability. Journal of Ecology, 106(5), 1902-1909.