The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
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:
the name of the function should be adapted to the parameters you want to fit:
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:
\(\lambda_i - \alpha_{ii}N_i - \alpha_{ij}Nj\)
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:
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.
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_
}
}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
simplywhere each part potentially depends on neighbour densities and covariates.
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.
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.