knitr::opts_chunk$set(message = FALSE, warning = FALSE)
Projecting species abundances
Population dynamics models, like the ones included in cxr
, are used for understanding the influence of different factors in shaping species densities, and also for inferring future dynamics based on the expected factors. We have included in cxr
basic functionality for projecting these dynamics. For that, we can use any of the default models included in the package, or -again- implement any user-defined model.
In this vignette, we show how to project the abundances of three species from our dataset a number of timesteps. For that, we first estimate model coefficients with the cxr_pm_multifit
function, as in previous vignettes.
library(cxr) data("neigh_list") three_sp <- c("BEMA","LEMA","SOAS") sp.pos <- which(names(neigh_list) %in% three_sp) data <- neigh_list[sp.pos] # keep only fitness and neighbours columns for(i in 1:length(data)){ data[[i]] <- data[[i]][,2:length(data[[i]])] } focal_column <- names(data) # Beverton-Holt model model_family <- "BH" # the "bobyqa" algorithm works best for these species optimization_method <- "bobyqa" # pairwise alphas, but no covariate effects alpha_form <- "pairwise" lambda_cov_form <- "none" alpha_cov_form <- "none" # no fixed terms, i.e. we fit both lambdas and alphas fixed_terms <- NULL # for demonstration purposes bootstrap_samples <- 3 # a limited number of timesteps. timesteps <- 10 # for demonstration purposes initial_abundances <- c(10,10,10) names(initial_abundances) <- three_sp # standard initial values, # not allowing for intraspecific facilitation initial_values <- list(lambda = 10, alpha_intra = 0.1, alpha_inter = 0.1) # lambda_cov = 0.1, # alpha_cov = 0.1) lower_bounds <- list(lambda = 1, alpha_intra = 0, alpha_inter = -1) # lambda_cov = 0, # alpha_cov = 0) upper_bounds <- list(lambda = 100, alpha_intra = 1, alpha_inter = 1) # lambda_cov = 1, # alpha_cov = 1)
With all initial values set, we can fit the parameters.
cxr_fit <- cxr_pm_multifit(data = data, focal_column = focal_column, model_family = model_family, # covariates = salinity, optimization_method = optimization_method, alpha_form = alpha_form, lambda_cov_form = lambda_cov_form, alpha_cov_form = alpha_cov_form, initial_values = initial_values, lower_bounds = lower_bounds, upper_bounds = upper_bounds, fixed_terms = fixed_terms, bootstrap_samples = bootstrap_samples)
Projecting abundances from a cxr
object is straightforward: just call the function abundance_projection
with the cxr
fit and the number of timesteps to project as well as the initial abundances. In this example, our fit did not include the effect of covariates, but if this was the case, we would also need to specify the value of each covariate in the projected timesteps (see the help of abundance_projection
for details).
ab <- abundance_projection(cxr_fit = cxr_fit, # covariates = covariates_proj, timesteps = timesteps, initial_abundances = initial_abundances)
This function returns a simple numeric matrix with the projected abundances. Lastly, here is a basic plot showing the projection. This is a pedagogical example and it does not make sense to interpret it ecologically, but it shows how one of the three species selected tends to become more abundant in the short term.
ab.df <- as.data.frame(ab) ab.df$timestep <- 1:nrow(ab.df) ab.df <- tidyr::gather(ab.df,key = "sp",value = "abund",-timestep) abund.plot <- ggplot2::ggplot(ab.df, ggplot2::aes(x = timestep, y = abund, group = sp)) + ggplot2::geom_line(ggplot2::aes(color = sp)) + ggplot2::ylab("number of individuals") + ggplot2::xlab("time") + ggplot2::ggtitle("Projected abundances of three plant species")+ NULL
abund.plot
Including your own projection model
As with other features of cxr
, you can implement your own model for projecting species abundances. If you have gone through vignette 4, you should already be familiar with the format of cxr
models, and how to make them available to the package.
Here, for reference, we show the complete code of the most complex Beverton-Holt model included. You can use this example as a template for the function name and arguments. The instructions of vignette 4 are also applicable here.
#' Beverton-Holt model for projecting abundances, #' with specific alpha values and global covariate effects on alpha and lambda #' #' @param lambda named numeric lambda value. #' @param alpha_intra single numeric value. #' @param alpha_inter numeric vector with interspecific alpha values. #' @param lambda_cov numeric vector with effects of covariates over lambda. #' @param alpha_cov named list of named numeric vectors #' with effects of each covariate over alpha values. #' @param abundance named numeric vector of abundances in the previous timestep. #' @param covariates matrix with observations in rows and covariates in named columns. #' Each cell is the value of a covariate in a given observation. #' #' @return numeric abundance projected one timestep #' @export BH_project_alpha_pairwise_lambdacov_global_alphacov_pairwise <- function(lambda, alpha_intra, alpha_inter, lambda_cov, alpha_cov, abundance, covariates){ # put together intra and inter coefficients, # be sure names match spnames <- names(abundance) alpha <- c(alpha_intra,alpha_inter) alpha <- alpha[spnames] alpha_covs <- list() for(ia in 1:length(alpha_cov)){ alpha_covs[[ia]] <- alpha_cov[[ia]][spnames] } numsp <- length(abundance) expected_abund <- NA_real_ # model num = 1 focal.cov.matrix <- as.matrix(covariates) for(v in 1:ncol(focal.cov.matrix)){ num <- num + lambda_cov[v]*focal.cov.matrix[,v] } cov_term_x <- list() for(v in 1:ncol(focal.cov.matrix)){ cov_temp <- focal.cov.matrix[,v] for(z in 1:length(abundance)){ #create alpha_cov_i*cov_i vector cov_term_x[[z+(length(abundance)*(v-1))]] <- # alpha_cov[z+(ncol(abund)*(v-1))] alpha_cov[[v]][z] * cov_temp } } cov_term <- list() for(z in 0:(length(abundance)-1)){ cov_term_x_sum <- cov_term_x[[z+1]] if(ncol(focal.cov.matrix) > 1){ for(v in 2:ncol(focal.cov.matrix)){ cov_term_x_sum <- cov_term_x_sum + cov_term_x[[v + length(abundance)]] } } cov_term[[z+1]] <- cov_term_x_sum } term <- 1 #create the denominator term for the model for(z in 1:length(abundance)){ term <- term + (alpha[z] + cov_term[[z]]) * abundance[z] } expected_abund <- (lambda * (num) / term) * abundance[names(lambda)] expected_abund }
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.