**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")
<- c("BEMA","LEMA","SOAS")
three_sp <- which(names(neigh_list) %in% three_sp)
sp.pos
<- neigh_list[sp.pos]
data # keep only fitness and neighbours columns
for(i in 1:length(data)){
<- data[[i]][,2:length(data[[i]])]
data[[i]]
}<- names(data)
focal_column
# Beverton-Holt model
<- "BH"
model_family
# the "bobyqa" algorithm works best for these species
<- "bobyqa"
optimization_method
# pairwise alphas, but no covariate effects
<- "pairwise"
alpha_form <- "none"
lambda_cov_form <- "none"
alpha_cov_form
# no fixed terms, i.e. we fit both lambdas and alphas
<- NULL
fixed_terms
# for demonstration purposes
<- 3
bootstrap_samples
# a limited number of timesteps.
<- 10
timesteps
# for demonstration purposes
<- c(10,10,10)
initial_abundances names(initial_abundances) <- three_sp
# standard initial values,
# not allowing for intraspecific facilitation
<- list(lambda = 10,
initial_values alpha_intra = 0.1,
alpha_inter = 0.1)
# lambda_cov = 0.1,
# alpha_cov = 0.1)
<- list(lambda = 1,
lower_bounds alpha_intra = 0,
alpha_inter = -1)
# lambda_cov = 0,
# alpha_cov = 0)
<- list(lambda = 100,
upper_bounds alpha_intra = 1,
alpha_inter = 1)
# lambda_cov = 1,
# alpha_cov = 1)
```

With all initial values set, we can fit the parameters.

```
<- cxr_pm_multifit(data = data,
cxr_fit 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).

```
<- abundance_projection(cxr_fit = cxr_fit,
ab # 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.

```
<- as.data.frame(ab)
ab.df $timestep <- 1:nrow(ab.df)
ab.df<- tidyr::gather(ab.df,key = "sp",value = "abund",-timestep)
ab.df <- ggplot2::ggplot(ab.df,
abund.plot ::aes(x = timestep,
ggplot2y = abund, group = sp)) +
::geom_line(ggplot2::aes(color = sp)) +
ggplot2::ylab("number of individuals") + ggplot2::xlab("time") +
ggplot2::ggtitle("Projected abundances of three plant species")+
ggplot2NULL
```

` 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
<- function(lambda,
BH_project_alpha_pairwise_lambdacov_global_alphacov_pairwise
alpha_intra,
alpha_inter,
lambda_cov,
alpha_cov,
abundance,
covariates){
# put together intra and inter coefficients,
# be sure names match
<- names(abundance)
spnames
<- c(alpha_intra,alpha_inter)
alpha <- alpha[spnames]
alpha <- list()
alpha_covs for(ia in 1:length(alpha_cov)){
<- alpha_cov[[ia]][spnames]
alpha_covs[[ia]]
}
<- length(abundance)
numsp <- NA_real_
expected_abund
# model
= 1
num <- as.matrix(covariates)
focal.cov.matrix for(v in 1:ncol(focal.cov.matrix)){
<- num + lambda_cov[v]*focal.cov.matrix[,v]
num
}<- list()
cov_term_x for(v in 1:ncol(focal.cov.matrix)){
<- focal.cov.matrix[,v]
cov_temp for(z in 1:length(abundance)){
#create alpha_cov_i*cov_i vector
+(length(abundance)*(v-1))]] <-
cov_term_x[[z# alpha_cov[z+(ncol(abund)*(v-1))]
* cov_temp
alpha_cov[[v]][z]
}
}<- list()
cov_term for(z in 0:(length(abundance)-1)){
<- cov_term_x[[z+1]]
cov_term_x_sum if(ncol(focal.cov.matrix) > 1){
for(v in 2:ncol(focal.cov.matrix)){
<- cov_term_x_sum +
cov_term_x_sum + length(abundance)]]
cov_term_x[[v
}
}+1]] <- cov_term_x_sum
cov_term[[z
}<- 1 #create the denominator term for the model
term for(z in 1:length(abundance)){
<- term + (alpha[z] + cov_term[[z]]) * abundance[z]
term
}<- (lambda * (num) / term) * abundance[names(lambda)]
expected_abund
expected_abund }
```