Getting started with cxr
Introduction
The cxr package provides a general interface to obtain
estimates of species vital rates and interaction coefficients between
species pairs from empirical data. These estimations are critical to
parameterize population models describing the dynamics of interacting
species. They also allow computing a series of metrics associated with
modern coexistence theory that inform about the likelihood of species to
coexist. These metrics are 1) niche differences that stabilize
coexistence between competing species and 2) average fitness
differences that drive competitive dominance and, in the absence of
niche differences, determine the superior competitor.
The package also allows exploring how environmental variation modifies both the intrinsic ability of species to produce offspring and the interaction coefficients between pairs of species (including with itself). This feature opens the possibility of exploring how stabilizing niche differences and average fitness differences vary across environmental gradients, and therefore, it allows analyzing whether the likelihood of species to coexist changes across environmental conditions.
Here we demonstrate the basic functionality of the package using a published observational dataset (see Lanuza et al. (2018) and vignette 2 (Data formats) for a description of the dataset). With this example, we will specifically estimate seed production in the absence of neighbors (lambda) and the strength and sign of species interactions between species pairs (alpha matrix). These values are the basis for estimating the degree of niche overlap (i.e 1- niche differences) and average fitness differences between species pairs, which are covered in vignette 3 (Coexistence metrics).
Finally, these estimations of lambda and alpha parameters are the basis for analyzing more complex dynamics such as the stability of the dynamics of multispecies communities (see Godoy et al. 2017) and multitrophic coexistence (see Godoy et al. 2018).
Fitting a single species
First, we load the package and the associated data. The included dataset contains, for each individual, its reproductive success and the number of neighbors per species in a 7.5 cm buffer (see vignette 2).
First, we draw the values of a single focal species.
my.sp <- "HOMA"
# get data from the list
obs_homa <- neigh_list[[my.sp]]
# no need for ID column
obs_homa <- subset(obs_homa,select = -c(obs_ID))
# For each observation, we need the individual plant fitness
# and the number of neighbours per species (in columns).
head(obs_homa)## fitness BEMA CETE CHFU CHMI HOMA LEMA MEEL MESU PAIN PLCO POMA POMO PUPA SASO
## 1 12 2 0 0 0 35 0 0 0 0 0 0 0 0 0
## 2 12 0 0 0 0 34 2 0 0 0 0 0 0 0 0
## 3 12 0 0 0 0 42 1 0 0 0 0 0 0 0 0
## 4 12 0 0 0 0 42 2 0 0 0 0 0 0 0 0
## 5 12 1 0 0 0 38 0 0 0 0 0 0 0 0 0
## 6 12 1 0 0 0 52 1 0 0 0 0 0 0 0 0
## SCLA SOAS SPRU
## 1 0 0 0
## 2 0 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 0 0Next, we estimate both the reproduction success in the absence of
neighbors (lambda) and the competition matrix (alpha). This is done by
fitting a model that mathematically relates the reproductive success to
the number of neighbors observed. In this first example, we fit the
selected species with a Ricker model (‘RK’ model family, see vignette 4)
and fairly standard initial values. The default optimization method
(Nelder-Mead) does not allow for lower or upper bounds in
model parameters, so these arguments are commented out. In ecological
terms, this optimization process allows estimating the strength of both
competitive and facilitative interactions, yet bounded optimization
algorithms can be used to restrict the analysis to either competition or
facilitation (i.e. positive or negative alpha values). We can also
specify whether we want to compute standard errors numerically, by
setting the argument bootstrap_samples to the number of
bootstrap replications for the calculation.
#?cxr_pm_fit #check the help file for a description of the arguments
fit_homa <- cxr_pm_fit(data = obs_homa,
focal_column = my.sp,
model_family = "BH",
covariates = NULL,
optimization_method = "Nelder-Mead",
alpha_form = "pairwise",
lambda_cov_form = "none",
alpha_cov_form = "none",
initial_values = list(lambda = 45,
alpha_intra = .1,
alpha_inter = .1),
#not aplicable to this optimazation method
# lower_bounds = list(lambda = 0,
# alpha_intra = 0,
# alpha_inter = 0),
# upper_bounds = list(lambda = 10,
# alpha_intra = 1,
# alpha_inter = 1),
fixed_terms = NULL,
# a low number of bootstrap samples
# for demonstration purposes,
# increase it for robust results.
bootstrap_samples = 3) For a quick summary of the fit, we can run a summary on the resulting object.
##
## model: 'BH_pm_alpha_pairwise_lambdacov_none_alphacov_none'
## optimization method: 'Nelder-Mead'
## ----------
## focal taxa ID: HOMA
## observations: 286
## neighbours: 17
## covariates: 0
## ----------
## focal lambda: 45.91924
## alpha_intra: 0.06096145
## mean alpha_inter: NA
## mean lambda_cov: - not fit -
## mean alpha_cov: - not fit -
## negative log-likelihood of the fit: 161.8463
## ----------This object is actually a list with several elements. We can thus access these elements as usual:
## [1] "model_name" "model"
## [3] "data" "focal_ID"
## [5] "optimization_method" "initial_values"
## [7] "fixed_terms" "lambda"
## [9] "alpha_intra" "alpha_inter"
## [11] "lambda_cov" "alpha_cov"
## [13] "lambda_standard_error" "alpha_intra_standard_error"
## [15] "alpha_inter_standard_error" "lambda_cov_standard_error"
## [17] "alpha_cov_standard_error" "log_likelihood"## lambda
## 45.91924## HOMA
## 0.06096145## BEMA CETE CHFU CHMI LEMA MEEL
## 0.111167394 NA NA NA 0.135454916 0.446891371
## MESU PAIN PLCO POMA POMO PUPA
## 0.247409952 0.376771962 0.268328939 -0.009558564 0.386357248 0.162581693
## SASO SCLA SOAS SPRU
## 0.097926694 0.565216139 NA 0.556207593Note that some interaction coefficients are set to NA because species do not co-occur but are nevertheless listed as neighbours with densities equal to zero in all focal observations. To check whether model fit reproduces well observed data, represent the fit according to obtained values.
require(ggplot2)
ggplot(obs_homa, aes(HOMA , fitness)) +
geom_point() +
stat_function(fun = function(x) fit_homa$lambda/(1+fit_homa$alpha_intra*x), lwd = 1.5, colour = "blue")
Fitting several species at once
Most likely users will want to fit model parameters to data from two
or more focal species. In order to do that with a single call, we
provide the function cxr_pm_multifit, which has a very
similar interface to cxr_pm_fit. Here we show how multiple
species can be fit using this function. For this multispecies case, rows
in the alpha element of the returning list correspond to species i and
columns to species j for each αij
coefficient. The diagonal corresponds to intraspecific coefficients. In
order to showcase other capabilities of the package, we include in this
example the effect of a covariate over the fitted lambda and alpha
parameters. This covariate, soil salinity, is also included as a dataset
in the package (see vignette 2). We consider that the covariate has a
linear effect in both the modification of lambda and alpha
parameters.
my.sp <- c("BEMA","CETE","LEMA")
obs_3sp <- neigh_list[my.sp]
# discard ID column
for(i in 1:length(obs_3sp)){
obs_3sp[[i]] <- obs_3sp[[i]][,2:length(obs_3sp[[i]])]
}
# load covariates: salinity
data("salinity_list")
salinity <- salinity_list[my.sp]
# keep only salinity column
for(i in 1:length(salinity)){
salinity[[i]] <- as.matrix(salinity[[i]][,2:length(salinity[[i]])])
colnames(salinity[[i]]) <- "salinity"
}Note how the data is passed in a list with as many elements as focal species. Each element is a dataframe with observations of the corresponding focal species. Same for the covariates data, it must be a list with as many elements as focal species. Each element is a dataframe (or a matrix) with a column for each covariate (one, in this case) and the same number of observations as its associated species data.
## [1] "BEMA" "CETE" "LEMA"## # A tibble: 6 × 18
## fitness BEMA CETE CHFU CHMI HOMA LEMA MEEL MESU PAIN PLCO POMA
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 116 1 0 0 0 0 0 0 0 0 0 34
## 2 68 0 0 0 0 0 1 0 0 0 0 47
## 3 36 0 0 0 0 0 0 0 0 0 0 38
## 4 64 0 0 0 0 0 5 0 0 0 0 21
## 5 144 2 0 0 0 0 0 0 0 0 0 55
## 6 56 1 0 0 0 0 4 0 0 0 0 17
## # ℹ 6 more variables: POMO <dbl>, PUPA <dbl>, SASO <dbl>, SCLA <dbl>,
## # SOAS <dbl>, SPRU <dbl>## [1] 287## salinity
## [1,] 0.9860
## [2,] 1.0280
## [3,] 1.0000
## [4,] 0.8830
## [5,] 0.8420
## [6,] 0.9121## [1] 287We fit the model as above, but using the cxr_pm_multifit
function. To allow proper convergence we fit each species with different
starting lambda values
fit_3sp <- list()
all_sp <- names(obs_3sp)
for (i in 1:length(obs_3sp)){
obs <- obs_3sp[i]
my.sp <- all_sp[i]
salinity.sp <- salinity[i]
lambda <- mean(obs[[1]]$fitness)
upper_lambda <- as.numeric(max(obs[[1]]$fitness))
fit <- cxr_pm_multifit(data = obs,
focal_column = my.sp,
model_family = "BH",
# here we use a bounded method for demonstration purposes
optimization_method = "bobyqa",
covariates = salinity.sp,
alpha_form = "pairwise",
lambda_cov_form = "global", # effect of covariates over lambda
alpha_cov_form = "global", # effect of covariates over alpha
initial_values = list(lambda = lambda,
alpha_intra = 1,
alpha_inter = 1,
lambda_cov = 0.1,
alpha_cov = 0.1),
lower_bounds = list(lambda = 0,
alpha_intra = -1,
alpha_inter = -1,
lambda_cov = -2,
alpha_cov = -2),
upper_bounds = list(lambda = upper_lambda,
alpha_intra = 3,
alpha_inter = 3,
lambda_cov = 2,
alpha_cov = 2),
# no standard errors
bootstrap_samples = 0)
fit_3sp[[i]] <- fit
}
names(fit_3sp) <- all_spWe can also have a glimpse of this multispecies fit for each species separately. This is an example for the first species.
## $model_name
## [1] "BH_pm_alpha_pairwise_lambdacov_global_alphacov_global"
##
## $data
## $data$BEMA
## # A tibble: 287 × 18
## fitness BEMA CETE CHFU CHMI HOMA LEMA MEEL MESU PAIN PLCO POMA
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 116 1 0 0 0 0 0 0 0 0 0 34
## 2 68 0 0 0 0 0 1 0 0 0 0 47
## 3 36 0 0 0 0 0 0 0 0 0 0 38
## 4 64 0 0 0 0 0 5 0 0 0 0 21
## 5 144 2 0 0 0 0 0 0 0 0 0 55
## 6 56 1 0 0 0 0 4 0 0 0 0 17
## 7 32 0 0 0 0 0 1 0 0 0 0 60
## 8 112 0 0 0 0 0 1 0 0 0 0 35
## 9 36 0 0 0 0 65 0 0 0 0 0 0
## 10 68 1 0 0 0 36 1 0 0 0 0 0
## # ℹ 277 more rows
## # ℹ 6 more variables: POMO <dbl>, PUPA <dbl>, SASO <dbl>, SCLA <dbl>,
## # SOAS <dbl>, SPRU <dbl>
##
##
## $taxa
## [1] "BEMA"
##
## $covariates
## $covariates$BEMA
## salinity
## [1,] 0.9860
## [2,] 1.0280
## [3,] 1.0000
## [4,] 0.8830
## [5,] 0.8420
## [6,] 0.9121
## [7,] 0.7450
## [8,] 0.8690
## [9,] 0.8380
## [10,] 0.7910
## [11,] 0.6830
## [12,] 0.9010
## [13,] 0.8750
## [14,] 0.8440
## [15,] 0.7600
## [16,] 0.6940
## [17,] 1.0220
## [18,] 0.8850
## [19,] 0.7680
## [20,] 0.9310
## [21,] 0.7770
## [22,] 0.7580
## [23,] 0.6410
## [24,] 0.6570
## [25,] 0.6820
## [26,] 0.7310
## [27,] 0.6400
## [28,] 0.5450
## [29,] 0.6180
## [30,] 0.5670
## [31,] 0.8070
## [32,] 0.5730
## [33,] 0.5160
## [34,] 0.6580
## [35,] 0.8520
## [36,] 0.5550
## [37,] 0.7780
## [38,] 0.5500
## [39,] 0.5750
## [40,] 0.7350
## [41,] 0.4780
## [42,] 0.4320
## [43,] 0.5840
## [44,] 0.4960
## [45,] 0.5890
## [46,] 0.5300
## [47,] 0.4820
## [48,] 0.6700
## [49,] 0.4120
## [50,] 0.3330
## [51,] 0.5700
## [52,] 0.5430
## [53,] 0.4460
## [54,] 0.5500
## [55,] 0.4460
## [56,] 0.5430
## [57,] 0.7520
## [58,] 0.4350
## [59,] 0.5950
## [60,] 0.7120
## [61,] 0.5580
## [62,] 0.6410
## [63,] 0.6180
## [64,] 0.4950
## [65,] 0.3450
## [66,] 0.3990
## [67,] 0.6940
## [68,] 0.8300
## [69,] 0.5260
## [70,] 0.4220
## [71,] 0.4990
## [72,] 0.5680
## [73,] 0.6400
## [74,] 0.6660
## [75,] 0.5740
## [76,] 0.6450
## [77,] 0.6160
## [78,] 0.5160
## [79,] 0.7340
## [80,] 0.6520
## [81,] 0.7900
## [82,] 0.7700
## [83,] 0.5970
## [84,] 0.8320
## [85,] 0.9160
## [86,] 0.5880
## [87,] 0.7190
## [88,] 0.9450
## [89,] 0.9820
## [90,] 0.7590
## [91,] 0.8740
## [92,] 0.8030
## [93,] 0.8330
## [94,] 0.6110
## [95,] 0.6060
## [96,] 0.6130
## [97,] 0.8390
## [98,] 0.7650
## [99,] 0.7770
## [100,] 0.4490
## [101,] 0.5760
## [102,] 0.5730
## [103,] 1.0510
## [104,] 1.0930
## [105,] 0.6910
## [106,] 0.5760
## [107,] 0.5790
## [108,] 0.5850
## [109,] 0.6500
## [110,] 0.6800
## [111,] 1.1320
## [112,] 0.5930
## [113,] 0.4930
## [114,] 0.6100
## [115,] 0.4660
## [116,] 0.5640
## [117,] 0.5780
## [118,] 0.5990
## [119,] 0.5240
## [120,] 0.4790
## [121,] 0.4800
## [122,] 0.7380
## [123,] 0.7820
## [124,] 0.7310
## [125,] 0.6230
## [126,] 0.6230
## [127,] 0.4760
## [128,] 0.7030
## [129,] 0.5390
## [130,] 0.4840
## [131,] 0.5710
## [132,] 0.5990
## [133,] 0.5160
## [134,] 0.4670
## [135,] 0.6530
## [136,] 0.7170
## [137,] 0.5560
## [138,] 0.5420
## [139,] 0.4400
## [140,] 0.6300
## [141,] 0.6800
## [142,] 0.9670
## [143,] 0.9320
## [144,] 0.6370
## [145,] 0.6330
## [146,] 0.6590
## [147,] 0.6420
## [148,] 0.7240
## [149,] 0.7560
## [150,] 0.4730
## [151,] 0.7330
## [152,] 0.9350
## [153,] 0.7970
## [154,] 1.0070
## [155,] 0.8100
## [156,] 0.9150
## [157,] 0.9080
## [158,] 0.8260
## [159,] 0.9330
## [160,] 0.7550
## [161,] 0.7800
## [162,] 0.9210
## [163,] 0.9690
## [164,] 0.7510
## [165,] 0.5980
## [166,] 0.6920
## [167,] 0.6490
## [168,] 0.6940
## [169,] 0.7460
## [170,] 0.6330
## [171,] 0.6420
## [172,] 0.6830
## [173,] 0.7570
## [174,] 0.7950
## [175,] 0.7470
## [176,] 0.6750
## [177,] 0.7630
## [178,] 0.7110
## [179,] 0.7550
## [180,] 0.7510
## [181,] 0.7860
## [182,] 0.9480
## [183,] 0.7310
## [184,] 0.6110
## [185,] 0.6170
## [186,] 0.6600
## [187,] 0.9440
## [188,] 0.7990
## [189,] 0.7130
## [190,] 0.5270
## [191,] 0.6250
## [192,] 0.8190
## [193,] 0.8790
## [194,] 0.6780
## [195,] 0.8260
## [196,] 0.6800
## [197,] 1.0160
## [198,] 0.8380
## [199,] 0.8990
## [200,] 0.9000
## [201,] 0.7830
## [202,] 0.7470
## [203,] 1.2460
## [204,] 1.1330
## [205,] 0.7870
## [206,] 0.6990
## [207,] 0.7900
## [208,] 0.8630
## [209,] 0.7320
## [210,] 0.7580
## [211,] 0.7380
## [212,] 0.8030
## [213,] 0.7610
## [214,] 0.8350
## [215,] 0.8430
## [216,] 0.9950
## [217,] 0.9260
## [218,] 0.7390
## [219,] 0.7720
## [220,] 0.9240
## [221,] 1.0160
## [222,] 1.1540
## [223,] 1.2620
## [224,] 0.8670
## [225,] 0.9310
## [226,] 1.0230
## [227,] 1.1400
## [228,] 1.2500
## [229,] 0.9640
## [230,] 1.0150
## [231,] 1.0250
## [232,] 0.9830
## [233,] 1.3380
## [234,] 1.0040
## [235,] 1.0440
## [236,] 1.0330
## [237,] 0.8630
## [238,] 1.0070
## [239,] 1.1030
## [240,] 0.7300
## [241,] 0.8030
## [242,] 0.9990
## [243,] 0.8900
## [244,] 0.8800
## [245,] 0.9570
## [246,] 0.7620
## [247,] 0.8340
## [248,] 0.9110
## [249,] 1.1080
## [250,] 0.6650
## [251,] 1.1050
## [252,] 1.1460
## [253,] 0.8250
## [254,] 0.6990
## [255,] 0.8270
## [256,] 1.0260
## [257,] 1.4740
## [258,] 1.4890
## [259,] 1.5450
## [260,] 1.6300
## [261,] 1.6520
## [262,] 1.3220
## [263,] 1.3940
## [264,] 1.3480
## [265,] 1.4130
## [266,] 1.5310
## [267,] 1.5360
## [268,] 1.2260
## [269,] 1.1680
## [270,] 1.2860
## [271,] 1.3910
## [272,] 1.1920
## [273,] 1.4640
## [274,] 1.4150
## [275,] 1.2810
## [276,] 1.5940
## [277,] 1.4100
## [278,] 1.4350
## [279,] 1.0920
## [280,] 1.0370
## [281,] 1.2750
## [282,] 1.0240
## [283,] 1.2080
## [284,] 1.6270
## [285,] 1.2850
## [286,] 1.2800
## [287,] 1.5360
##
##
## $optimization_method
## [1] "bobyqa"
##
## $initial_values
## $initial_values$lambda
## [1] 82.9547
##
## $initial_values$alpha_intra
## [1] 1
##
## $initial_values$alpha_inter
## [1] 1
##
## $initial_values$lambda_cov
## [1] 0.1
##
## $initial_values$alpha_cov
## [1] 0.1
##
##
## $fixed_terms
## NULL
##
## $lambda
## BEMA
## 83.15007
##
## $alpha_matrix
## BEMA CETE CHFU CHMI HOMA LEMA MEEL MESU PAIN
## BEMA 0.8296157 NA NA NA 0.2511788 0.9555702 1.263203 0.2699469 1.23746
## PLCO POMA POMO PUPA SASO SCLA SOAS
## BEMA 0.2058185 0.03381329 1.22791 1.261585 0.4169347 1.410967 1.260929
## SPRU
## BEMA 0.2401413
##
## $lambda_cov
## salinity
## BEMA 1.984515
##
## $alpha_cov
## $alpha_cov$salinity
## BEMA CETE CHFU CHMI HOMA LEMA MEEL
## BEMA 0.3886128 0.3886128 0.3886128 0.3886128 0.3886128 0.3886128 0.3886128
## MESU PAIN PLCO POMA POMO PUPA SASO
## BEMA 0.3886128 0.3886128 0.3886128 0.3886128 0.3886128 0.3886128 0.3886128
## SCLA SOAS SPRU
## BEMA 0.3886128 0.3886128 0.3886128
##
##
## $lambda_standard_error
## NULL
##
## $alpha_standard_error
## NULL
##
## $lambda_cov_standard_error
## NULL
##
## $alpha_cov_standard_error
## NULL
##
## $log_likelihood
## BEMA
## 506.9143
##
## attr(,"class")
## [1] "cxr_pm_multifit"The numerical estimation of parameters depends on the model with
which to estimate fitness values, the optimization method, and the
underlying data. In our example dataset, some species are better
represented than others, and the function will raise warnings if the
estimation can be improved or, for example, if any fitted parameter is
equal to the lower or upper bounds provided. The
cxr_pm_multifit function will behave conservatively and
return NULL if parameter fitting fails for any of
the focal species passed as arguments, printing an informative message
about which species failed to fit. In such cases, users may either fit
each species separately or call cxr_pm_multifit without the
problematic species.
Importantly, bounded methods can be very sensitive to the initial
values and bounds, so, as in any numerical optimization problem, you
should double-check the values obtained, at least by computing standard
errors on your parameters with an adequate number of boostrap samples,
and preferably by performing sensitivity analyses on the different
parameters (in this and other vignettes, we have not included any of
these checks, as our examples are merely for demonstration purposes).
Aside from these recommendations, in the cxr objects
returned from our functions, the negative log-likelihood of the fit is
also included, which may be useful in helping users choose a certain
optimization algorithm for a particular dataset. Remember that the
convention is to present negative log-likelihood, and more negative
values are better. In general, it is recommended to test different
optimization algorithms, as they may produce fairly different results
(Mullen 2014). In this example, the method used (“bobyqa”, a
well-established bounded optimization algorithm) returns the following
negative log-likelihood values:
## BEMA
## 506.9143Including environmental variability
In the above example for fitting multiple species at once, we have
already offered a glimpse on how to include the effect of environmental
covariates over lambda and alpha values. The relevant arguments in
cxr_pm_fit and cxr_pm_multifit are
‘lambda_cov_form’ and ‘alpha_cov_form’. If these are set to ‘none’, no
effect of covariates is considered. Otherwise, there are a couple
options to model this effect (all the options consider linear effects,
for now). First, the effect of covariates over lambda can be ‘global’,
meaning that each covariate has a global parameter affecting lambda.
This is formulated as
$\lambda (1 + \sum_{k=1}^{s} \theta_k c_k)$
where s is the number of
environmental covariates, ck is the
observed value of the i-th covariate and θk is the
‘lambda_cov’ parameter of the cxr functions.
The effect over alpha values can likewise be ‘global’, so that, for focal species i, each covariate affects the alpha values αij equally through a global parameter:
$\alpha_* + \sum_{k=1}^{s} \psi_k c_k$
where * represents the set of all
pairwise interaction for a given focal species, and ψk is the global
parameter relating covariate k
to the alpha values. The last possibility included in cxr
is for the effect of covariates over alpha values to be
interaction-specific, which is coded by specifying ‘alpha_cov_form’ as
‘pairwise’:
$\alpha_{ij} + \sum_{k=1}^{s} \psi_{ijk} c_k$
In this case, each covariate ck will have a specific parameter ψijk for its effect over each interaction αij.
We can retrieve the fitted values for ‘lambda_cov’ (θ) and ‘alpha_cov’ (ψ) simply from the output of the function:
## salinity
## BEMA 1.984515## $salinity
## BEMA CETE CHFU CHMI HOMA LEMA MEEL
## BEMA 0.3886128 0.3886128 0.3886128 0.3886128 0.3886128 0.3886128 0.3886128
## MESU PAIN PLCO POMA POMO PUPA SASO
## BEMA 0.3886128 0.3886128 0.3886128 0.3886128 0.3886128 0.3886128 0.3886128
## SCLA SOAS SPRU
## BEMA 0.3886128 0.3886128 0.3886128Parameter boundaries
cxr includes the possibility of providing boundaries for
your parameters, so that the fitted values fall inside a certain
interval. This is done via the arguments lower_bounds and
upper_bounds of the fitting functions. These arguments
accept a list whose elements must be named and correspond to one or more
of the function parameters, i.e. lambda,
alpha_intra, alpha_inter,
lambda_cov, and alpha_cov. If these boundaries
are of length 1, as in this example:
lower_bounds = list(lambda = 0,
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)the same boundaries will be used for all values in the appropriate
set, e.g. each interspecific alpha term (alpha_inter) will
have the same boundaries. Otherwise, if you are interested in fitting
parameters for more than one covariate, you may want to specify varying
boundaries for lambda_cov and alpha_cov
depending on the different covariates. In that case, the functions
accept, for these elements, vectors with as many values as covariates.
For example, if you have two covariates for which you want to calculate
associated lambda_cov and alpha_cov values,
but these covariates have very different magnitudes, you may provide two
lambda_cov and alpha_cov boundaries.
# fit three species, as in the previous example
data <- neigh_list[1:3]
# keep only fitness and neighbours columns
for(i in 1:length(data)){
data[[i]] <- data[[i]][,2:length(data[[i]])]
}
focal_column <- names(data)
# covariates
# in this example, we use salinity
# and two more covariates, randomly generated
data("salinity_list")
# observations for the first three species
cov_list <- salinity_list[1:3]
for(i in 1:length(cov_list)){
# keep only salinity column
cov_list[[i]] <- as.matrix(cov_list[[i]][,2:length(cov_list[[i]])])
# add two random covariates
cov_list[[i]] <- cbind(cov_list[[i]],
runif(nrow(cov_list[[i]]),1,10),
runif(nrow(cov_list[[i]]),10,100))
colnames(cov_list[[i]]) <- c("salinity","cov2","cov3")
}
# this is how each element of the covariates list looks like
head(cov_list[[1]])## salinity cov2 cov3
## [1,] 0.9860 6.012198 83.14349
## [2,] 1.0280 3.503852 79.85387
## [3,] 1.0000 5.377400 84.40572
## [4,] 0.8830 2.721512 91.69667
## [5,] 0.8420 7.907873 55.78722
## [6,] 0.9121 7.482548 90.60003# function parameters
model_family <- "BH"
covariates <- cov_list
# bobyqa is generally more robust than other bounded methods
optimization_method <- "bobyqa"
alpha_form <- "pairwise"
lambda_cov_form <- "global"
alpha_cov_form <- "pairwise"
# note how lambda_cov and alpha_cov
# have different initial values for each covariate effect
# the commented assignations are also possible,
# giving equal initial values to all parameters
initial_values = list(lambda = 1,
alpha_intra = 0.1,
alpha_inter = 0.1,
lambda_cov = c(0.1,0.2,0.1),
alpha_cov = c(0.1,0.2,0.1))
# lambda_cov = c(0.1),
# alpha_cov = c(0.1))
# same with boundaries
lower_bounds = list(lambda = 0,
alpha_intra = 0,
alpha_inter = -1,
lambda_cov = c(-1,0,-1),
alpha_cov = c(-1,0,-1))
# lambda_cov = c(-1),
# alpha_cov = c(-1))
upper_bounds = list(lambda = 100,
alpha_intra = 1,
alpha_inter = 1,
lambda_cov = c(1,2,1),
alpha_cov = c(1,2,1))
# lambda_cov = c(1),
# alpha_cov = c(1))
fixed_terms <- NULL
bootstrap_samples <- 3# this fit is fairly complex, it may take a while,
# it also raises warnings, suggesting that either the data,
# model, or initial values/boundaries can be improved.
# This is consistent with having observational data
# and, furthermore, random covariates
fit_multi_cov <- cxr_pm_multifit(data = data,
focal_column = focal_column,
model_family = model_family,
covariates = covariates,
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)Keeping parameters fixed
There is the option to provide fixed values for model parameters,
using the fixed_terms argument in the functions. As an
example, if you want to keep the lambda parameter fixed,
for a single species fit, you would modify the function arguments
accordingly:
fixed_terms <- list(lambda = 1)
# now lambda does not appear in 'initial_values'
initial_values <- list(alpha_intra = 0,
alpha_inter = 0,
lambda_cov = 0,
alpha_cov = 0)
# lower and upper bounds should, likewise,
# not contain lambdaThis is extended to the cxr_pm_multifit function, by
having fixed_terms be a list with as many elements as focal
species. For example, having three focal species, where we want to keep
lambda fixed:
fixed_terms <- list(list(lambda = 1), # focal sp 1
list(lambda = 1.2), # focal sp 2
list(lambda = 1.3)) # focal sp 3In this release of cxr, you must specify the same
initial values for every focal species when using the
cxr_pm_multifit option. This means, in practice, that if
you want to fit species with different sets of fixed parameters,
(e.g. only lambda for one species, but lambda
and alpha_intra for another), or with different initial
values, you should fit them separately, using
cxr_pm_fit.
References
Godoy, O., Stouffer, D. B., Kraft, N. J., & Levine, J. M. (2017). Intransitivity is infrequent and fails to promote annual plant coexistence without pairwise niche differences. Ecology, 98(5), 1193-1200.
Lanuza, J. B., Bartomeus, I., & Godoy, O. (2018). Opposing effects of floral visitors and soil conditions on the determinants of competitive outcomes maintain species diversity in heterogeneous landscapes. Ecology letters, 21(6), 865-874.
Godoy, O., Bartomeus, I., Rohr, R. P., & Saavedra, S. (2018). Towards the integration of niche and network theories. Trends in ecology & evolution, 33(4), 287-300.
Mullen, K. M. (2014). Continuous global optimization in R. Journal of Statistical Software, 60(6), 1-45.