Introduction
The cxr package facilitates the estimation of key
metrics from modern coexistence theory (MCT, Chesson 2000), by obtaining
species vital rates and interactions coefficients from population
dynamics models. The metrics that can be obtained with cxr
relate to the idea pioneered by Chesson that the degree to which two
species can coexist depends both on their stabilizing niche differences
and on their average fitness differences. In this vignette we review the
metrics that can be computed with cxr, starting with those
related to niche differences and moving on to the more complex situation
of average fitness differences and its different demographic and
density-dependent components.
Stabilizing niche differences
cxr allows the calculation of niche overlap (ρ) between pairs of species, from
which niche differences can be easily obtained as 1 − ρ. Conceptually, if species
limit themselves much more than they limit their competitors niche
overlap will be very low and niche differences will be close to the
maximum (i.e. 1). Conversely, if species limit themselves and other
species similarly niche overlap will be very high and niche differences
will be close to zero. In the absence of niche differences, the species
with higher fitness (next section) is defined as the superior
competitor. For obtaining niche overlap, thus, we need interaction
coefficients among pairs of species.
First, we specify the data to use and the starting parameter values and bounds.
library(cxr)
data("neigh_list")
data <- neigh_list
# keep only fitness and neighbours columns
for(i in 1:length(data)){
data[[i]] <- data[[i]][,2:length(data[[i]])]
}
focal_column <- names(data)
model_family <- "RK"
optimization_method <- "bobyqa" # we use a bounded method
alpha_form <- "pairwise"
lambda_cov_form <- "none"
alpha_cov_form <- "none"
initial_values = list(lambda = 1,
alpha_intra = 0.1,
alpha_inter = 0.1)
lower_bounds = list(lambda = 0,
alpha_intra = 0.01,
alpha_inter = 0.01)
upper_bounds = list(lambda = 100,
alpha_intra = 1,
alpha_inter = 1)
fixed_terms <- NULL
bootstrap_samples <- 0Now we obtain, for each focal species, values for lambda (offspring production in the absence of interactions, in our dataset average viable seed production per species) and alpha (intra and interspecific pairwise interaction coefficients).
all.sp.fit <- cxr_pm_multifit(data = data,
model_family = model_family,
focal_column = focal_column,
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)With these parameters, we can compute pairwise niche differences as 1 - niche overlap. Niche overlap following Chesson’s (2013) definition can only include interaction coefficients with positive values (i.e. competition). Niche overlap can range from 0 to infinity, which implies that niche differences can range from minus infinity to zero. Negative niche differences can be interpreted as a signature of priority effects, meaning in a very broad sense that the first species of a pair to arrive to a community is the winner. Niche differences between 0 and 1 reflect the degree to which species A and B limit each other compared to themselves. Within this range of values, species might coexist or be competitive excluded according to how stabilizing niche differences offset average fitness differences (see Adler et al. (2007) for more details).
Finally, it is worth noting that for those cases in which interaction
coefficients are negative (i.e. negative interaction coefficients in
Beverton-Holt or Ricker models imply facilitative interactions), the
cxr package provides another definition of niche
differences, developed by Saavedra et al. (2017). This definition of
niche differences uses another set of tools based on the Structural
Approach (SA) to species coexistence, and while it is qualitatively
coherent with the MCT framework, niche differences from both definitions
do not exactly match (see Saavedra et al. (2017) for more details).
The function niche_overlap computes both definitions
(MCT and SA). If the input is a cxr_pm_multifit object it
will compute the niche overlap between all pairs of focal species (but
it can also accept other arguments, see the help of the function for
details).
niche_overlap_all_pairs <- niche_overlap(cxr_multifit = all.sp.fit)
# as not all species combinations occur,
# several interaction coefficients are NA
# which, in turn, return NA values for niche overlap
# so, remove them and check the first computed values
niche_overlap_all_pairs <- niche_overlap_all_pairs[complete.cases(niche_overlap_all_pairs),]
head(niche_overlap_all_pairs)## sp1 sp2 niche_overlap_MCT niche_overlap_SA
## 4 BEMA HOMA 0.5907839 0.7137810
## 5 BEMA LEMA 0.5750315 0.6965845
## 6 BEMA MEEL 1.3114612 0.9136775
## 7 BEMA MESU 0.7446739 0.8561565
## 8 BEMA PAIN 2.4908512 0.8219943
## 9 BEMA PLCO 0.3721694 0.4542872Thus, niche differences are simply
niche_overlap_all_pairs$MCT_niche_diff <- 1 - niche_overlap_all_pairs$niche_overlap_MCT
niche_overlap_all_pairs$SA_niche_diff <- 1 - niche_overlap_all_pairs$niche_overlap_SAAverage fitness differences between pair of species
The second type of pairwise differences relevant to species coexistence are the average fitness differences. These are, properly speaking, a ratio, and reflect the degree to which one species is a superior competitor over another. Conceptually, they range from 1 to infinity, so that when fitness differences are equal to 1, both species are equivalent competitors. Ratios lower than 1 mean than the denominator species is the superior competitor, so that the standard ratio is calculated with the superior competitor in the numerator.
Average fitness differences are defined by two components, the ‘demographic ratio’ and the ‘competitive response ratio’. The ‘demographic ratio’ is a density independent term that describes the degree to which one species (species j) has higher offspring production than another species (species i). The ‘competitive response ratio’ is a density dependent term, which describes the degree to which species i is more sensitive to both intra and interspecific competition than species j. Because both ratios define the average fitness differences, this means that a species can be a superior competitor because it produces high number of for instance seeds/eggs or because its offspring production is little reduced in the presence of competitors. This formulation is explained in greater detail by Godoy and Levine (2014).
Both components of average fitness differences can be computed with
cxr. As it is the case with the niche overlap calculations,
average fitness differences in the context of MCT are only defined for
negative interactions, i.e. positive alpha values. Saavedra et
al. (2017) developed a structural analog of this metric, that we also
include in our package. Thus, the function avg_fitness_diff
returns the demographic ratio, competitive response ratio, and average
fitness differences for each species pair in the context of MCT, as well
as the structural analog of these differences. It accepts the same
arguments as the niche_overlap function (see help for
details), and here we use the multispecies fit obtained above to
calculate fitness differences across all pairs of focal species.
avg_fitness_diff_all_pairs <- avg_fitness_diff(cxr_multifit = all.sp.fit)
# as with niche overlap, remove NA values
avg_fitness_diff_all_pairs <- avg_fitness_diff_all_pairs[complete.cases(
avg_fitness_diff_all_pairs),]
# average fitness ratio of sp1 over sp2
# if < 1, sp2 is the superior competitor,
# and the average fitness difference is the inverse ratio,
# i.e. sp2 over sp1.
head(avg_fitness_diff_all_pairs)## sp1 sp2 demographic_ratio competitive_response_ratio average_fitness_ratio
## 1 BEMA BEMA 1.0000000 1.000000 1.000000
## 5 HOMA BEMA 0.6990714 1.692666 1.183295
## 6 LEMA BEMA 1.0253884 1.647534 1.689362
## 7 MEEL BEMA 0.6350817 3.757493 2.386315
## 8 MESU BEMA 1.0470844 2.133580 2.234038
## 9 PAIN BEMA 0.7667031 7.136586 5.471643
## structural_fitness_diff
## 1 0.000000
## 5 17.213741
## 6 15.298219
## 7 3.996220
## 8 25.177242
## 9 8.488462Estimation of species’ competitive ability
From average fitness differences, we can compute species’ competitive ability. The formulation of species’ competitive ability changes according to the population model selected to describe the dynamics of interacting species. The mathematical procedure to obtain the definition of species competitive ability was firstly described in Godoy and Levine (2014) for the annual plant model, and expanded to a wider family of models by Hart et al. (2018).
Competitive ability can be calculated given a set of species
parameters and model family (e.g. as specified in cxr
objects). Again, the set of arguments accepted by the
competitive_ability function is the same as the arguments
passed to the niche_overlap and
avg_fitness_diff functions. The easiest way is to provide a
cxr_pm_multifit object with all necessary information.
competitive_ability_all_pairs <- competitive_ability(cxr_multifit = all.sp.fit)
competitive_ability_all_pairs <- competitive_ability_all_pairs[complete.cases(
competitive_ability_all_pairs),]
head(competitive_ability_all_pairs)## sp1 sp2 competitive_ability_sp1
## 5 HOMA BEMA 307.4578
## 6 LEMA BEMA 438.9503
## 7 MEEL BEMA 279.3146
## 8 MESU BEMA 460.5170
## 9 PAIN BEMA 337.2028
## 10 PLCO BEMA 238.7581Estimation of species fitness in the absence of niche differences
In the previous steps, average fitness differences and therefore species competitive ability are computed combining density independent and density dependent effects. This means that these metrics are estimated in the presence of stabilizing niche differences. However in some cases, it can be interesting for the user to estimate species fitness in the absence of niche differences. This is the case for instance if we want to list species according to a competitive hierarchy. With this procedure, we would be able to tease apart which species would be the first superior competitor if we remove niche differences, which would be the second superior competitor and so on, up to the weakest competitor. In order to define this species fitness against all other competitors, we need to collapse pairwise interaction coefficients into two components: how species respond to overall competition (competitive response) and how species affect other species (competitive effect). Both components are defined in Godoy et al. (2014).
In cxr, we first need to obtain these components from
observational data, and this calculation is done with the
cxr_er_fit function. This function is similar to
cxr_pm_fit, with some caveats. It accepts a list of
observational dataframes, but in this case the number of observations of
each focal species must match (this is in order to compute balanced
parameters). Furthermore, the set of focal species needs to be the same
as the set of neighbours, unlike in cxr_pm_fit.
We first set the initial data and values
data("neigh_list")
# For obtaining effect and responses, all species
# need to have the same number of observations.
# We selct 3 species that have >250 observations
names(neigh_list)## [1] "BEMA" "CETE" "CHFU" "CHMI" "HOMA" "LEMA" "MEEL" "MESU" "PAIN" "PLCO"
## [11] "POMA" "POMO" "PUPA" "SASO" "SCLA" "SOAS" "SPRU"## BEMA CETE CHFU CHMI HOMA LEMA MEEL MESU PAIN PLCO POMA POMO PUPA SASO SCLA SOAS
## 287 10 160 5 288 273 76 229 108 169 90 16 98 290 100 87
## SPRU
## 44# BEMA, HOMA, LEMA, SASO, have > 250 observations.
example_sp <- c(1,5,6) #corresponds to c("BEMA","HOMA","LEMA")
n.obs <- 250
data <- neigh_list[example_sp]
# use a bounded optimization method
optimization_method <- "L-BFGS-B"
# no fixed terms, i.e. we fit all parameters
fixed_terms <- NULL
# according to a Ricker model (for consistency with previous examples)
model_family <- "RK"
# no standard error calculation in this example
bootstrap_samples <- 0
# keep only fitness and neighbours columns
# and subset to 'n.obs' rows
for(i in 1:length(data)){
data[[i]] <- data[[i]][1:n.obs,c(2,example_sp+2)]#2:length(data[[i]])]
}
# set initial values and bounds
initial_values_er = list(lambda = 10,
effect = 1,
response = 1)
lower_bounds_er = list(lambda = 1,
effect = 0.1,
response = 0.1)
upper_bounds_er = list(lambda = 100,
effect = 10,
response = 10)and obtain the maximum-likelihood estimation of competitive effects and responses.
er.fit <- cxr_er_fit(data = data,
model_family = model_family,
optimization_method = optimization_method,
initial_values = initial_values_er,
lower_bounds = lower_bounds_er,
upper_bounds = upper_bounds_er,
fixed_terms = fixed_terms,
bootstrap_samples = bootstrap_samples)With this information, we can readily obtain specific species fitness
by passing the cxr_er_fit object as argument to the
species_fitness function.
## fitness_BEMA fitness_HOMA fitness_LEMA
## 4.695324 15.291338 23.486368References
Adler, P. B., HilleRisLambers, J., & Levine, J. M. (2007). A niche for neutrality. Ecology letters, 10(2), 95-104.
Chesson, P. (2000). Mechanisms of maintenance of species diversity. Annual review of Ecology and Systematics, 31(1), 343-366.
Chesson, P. (2013). Species competition and predation. In Ecological systems (pp. 223-256). Springer, New York, NY.
Godoy, O., & Levine, J. M. (2014). Phenology effects on invasion success: insights from coupling field experiments to coexistence theory. Ecology, 95(3), 726-736.
Godoy, O., Kraft, N. J., & Levine, J. M. (2014). Phylogenetic relatedness and the determinants of competitive outcomes. Ecology Letters, 17(7), 836-844.
Hart, S. P., Freckleton, R. P., & Levine, J. M. (2018). How to quantify competitive ability. Journal of Ecology, 106(5), 1902-1909.
Saavedra, S., Rohr, R. P., Bascompte, J., Godoy, O., Kraft, N. J., & Levine, J. M. (2017). A structural approach for understanding multispecies coexistence. Ecological Monographs, 87(3), 470-486.