**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 (\(\rho\)) between pairs of species, from which niche differences can be easily obtained as \(1 - \rho\). 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")
<- neigh_list
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 <- "RK"
model_family <- "bobyqa" # we use a bounded method
optimization_method <- "pairwise"
alpha_form <- "none"
lambda_cov_form <- "none"
alpha_cov_form = list(lambda = 1,
initial_values alpha_intra = 0.1,
alpha_inter = 0.1)
= list(lambda = 0,
lower_bounds alpha_intra = 0.01,
alpha_inter = 0.01)
= list(lambda = 100,
upper_bounds alpha_intra = 1,
alpha_inter = 1)
<- NULL
fixed_terms <- 0 bootstrap_samples
```

Now 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).

```
<- cxr_pm_multifit(data = data,
all.sp.fit 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(cxr_multifit = all.sp.fit)
niche_overlap_all_pairs # 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[complete.cases(niche_overlap_all_pairs),]
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.4542872
```

Thus, niche differences are simply

```
$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_SA niche_overlap_all_pairs
```

**Average 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(cxr_multifit = all.sp.fit)
avg_fitness_diff_all_pairs
# as with niche overlap, remove NA values
<- avg_fitness_diff_all_pairs[complete.cases(
avg_fitness_diff_all_pairs
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.488462
```

**Estimation 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(cxr_multifit = all.sp.fit)
competitive_ability_all_pairs <- competitive_ability_all_pairs[complete.cases(
competitive_ability_all_pairs
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.7581
```

**Estimation 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"
```

`sapply(neigh_list,nrow)`

```
## 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.
<- c(1,5,6) #corresponds to c("BEMA","HOMA","LEMA")
example_sp <- 250
n.obs <- neigh_list[example_sp]
data
# use a bounded optimization method
<- "L-BFGS-B"
optimization_method
# no fixed terms, i.e. we fit all parameters
<- NULL
fixed_terms
# according to a Ricker model (for consistency with previous examples)
<- "RK"
model_family
# no standard error calculation in this example
<- 0
bootstrap_samples
# keep only fitness and neighbours columns
# and subset to 'n.obs' rows
for(i in 1:length(data)){
<- data[[i]][1:n.obs,c(2,example_sp+2)]#2:length(data[[i]])]
data[[i]]
}
# set initial values and bounds
= list(lambda = 10,
initial_values_er effect = 1,
response = 1)
= list(lambda = 1,
lower_bounds_er effect = 0.1,
response = 0.1)
= list(lambda = 100,
upper_bounds_er effect = 10,
response = 10)
```

and obtain the maximum-likelihood estimation of competitive effects and responses.

```
<- cxr_er_fit(data = data,
er.fit 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.

```
<- species_fitness(er.fit)
spfitness spfitness
```

```
## fitness_BEMA fitness_HOMA fitness_LEMA
## 44.63158 21.14203 42.08254
```

**References**

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.