# Coexistence metrics

#### David Garcia-Callejas and cxr team

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")
    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 <- 0

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

    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.4542872

Thus, 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_SA

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_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.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_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.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.
    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.

    spfitness <- species_fitness(er.fit)
    spfitness

    ## fitness_BEMA fitness_HOMA fitness_LEMA
    ##     4.695324    15.291338    23.486368

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.