SDMfit is used to fit a single species SDM, what we call a 'local model' of trophicSDM. It returns an object of class 'SDMfit'. Requires basically the same inputs of trophicSDM, with the requirement to specify with the parameter 'focal' the species that is modeled by the SDMfit.
SDMfit(
focal,
Y,
X,
G,
formula.foc,
sp.formula = NULL,
sp.partition = NULL,
mode = "prey",
method = "stan_glm",
family,
penal = NULL,
iter = 1000,
chains = 2,
verbose = TRUE
)the name of the species to be modeled
The sites x species matrix containing observed species distribution (e.g. presence-absence).
The design matrix, i.e. sites x predictor matrix containing the value of each explanatory variable (e.g. the environmental conditions) at each site.
The species interaction network (aka metaweb). Needs to be an igraph object. Links must go from predator to preys. It needs to be a directed acyclic graph.
The formula for the abiotic part of the species distribution model.
(optional) It allows to specify a particular definition of the biotic part of the model, e.g., using composite variables (e.g., richness), or an interaction of the biotic and abiotic component. More details in 'Details'.
(optional) a list to specify groups of species that are used to compute composite variables, e.g., a species can be modeled as a function of the richness of each group of preys. It has to be a list, each element is a vector containing the names of species in the group.
"prey" if bottom-up control (default), "predators" otherwise. Notice that G needs to be such that links point from predators to prey.
which SDM method to use. For now the available choices are: "glm" (frequentist) or "stan_glm" (full Bayesian MCMC, default). Notice that using "glm" does not allow error propagation when predicting.
the family parameter of the glm function (see glm). family=gaussian(link ="identity") for gaussian data or family=binomial(link = "logit") or binomial(link = "probit") for presence-absence data.
(optional, default to NULL) Penalisation method to shrink regression coefficients.If NULL (default), the model does not penalise the regression coefficient. For now, available penalisation method are "horshoe" for stan_glm, "elasticnet" for glm and "coeff.signs" (prey coefficients are set to positive and predator coefficients to negative) for glm and stan_glm.
(for method="stan_glm" only) Number of iterations for each MCMC chain if stan_glm is used
(for method="stan_glm" only) Number of MCMC chains (default to 2)
Whether to print algorithm progresses
A list containing 'm', a "SDMfit" object and 'form.all', a string describing the formula of the SDMfit object. The "SDM" fit object contains:
The output of the function used to fit the SDM. E.g., an object of class "glm" is method = "glm", an object of class "stanreg" if method = "stan_glm".
A numeric vector of standard errors on parameters
The formula used to fit the SDM (both abiotic and biotic terms)
The input parameters used to fit the SDM.
The name of the species modeled
The model.frame data.frame used to fit the model
The inferred coefficients (with credible intervals or p-values when available)
The AIC of the local model
The log.likelihood of the local model
"sp.formula" and "sp.partition" can be combined to define any kind of composite variables for the biotic part of the formula. "sp.formula" can be :
A string defining a formula as function of "richness". E.g., sp.formula="richness+I(richness)^2" (species are modeled as a function of a quadratic polynomial of their prey richness), "I(richness>0)" (species are modeled as a function of a dummy variable that is equal to 1 when at least one species is present). Importantly, when group of preys (or predators) are specified by "sp.partition", species are modeled as a function of the composite variable specified by "sp.formula" for each of their prey groups.
A more flexible option is to specify sp.formula as a list (whose names are species' names) that contains for each species the definition biotic part of the model. Notice that, in this case, the function does not check that the model is a DAG. This allow to define any kind of composite variable, or to model interactions between environmental covariates and preys (or predators).
data(Y,X,G)
# Run a local model (i.e. a SDM) for species Y6
mySDM = SDMfit("Y6", Y, X, G, "~X_1 + X_2", mode = "prey",
method = "stan_glm", family = binomial(link = "logit"))
#>
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 5.7e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.57 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup)
#> Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling)
#> Chain 1: Iteration: 1000 / 1000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 0.346 seconds (Warm-up)
#> Chain 1: 0.285 seconds (Sampling)
#> Chain 1: 0.631 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 2).
#> Chain 2:
#> Chain 2: Gradient evaluation took 8.9e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.89 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 2: Iteration: 1 / 1000 [ 0%] (Warmup)
#> Chain 2: Iteration: 501 / 1000 [ 50%] (Sampling)
#> Chain 2: Iteration: 1000 / 1000 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 0.213 seconds (Warm-up)
#> Chain 2: 0.228 seconds (Sampling)
#> Chain 2: 0.441 seconds (Total)
#> Chain 2:
mySDM$m
#> ==================================================================
#> Local SDMfit for species Y6 with no penalty , fitted using stan_glm
#> ==================================================================
#> * Useful S3 methods
#> print(), coef(), plot(), predict()
#> $model gives the stanreg class object
#> ==================================================================