R/species_mix-class.R, R/species_mix_s3-class.R
species_mix.multifit.RdThis version of species mix is useful for fitting models which have complex likelihoods. The multiple starts will enable optimization of the log-likelihood using multiple starts.
species_mix.multifit(
archetype_formula = NULL,
species_formula = stats::as.formula(~1),
data,
nArchetypes = 3,
family = "bernoulli",
offset = NULL,
weights = NULL,
spp_weights = NULL,
size = NULL,
power = 1.6,
control = list(ecm.prefit = FALSE),
inits = NULL,
titbits = FALSE,
nstart = 10,
mc.cores = 1
)
# S3 method for species_mix.multifit
print(x, ...)an object of class "formula" (or an object that can be coerced to that class). The response variable (left hand side of the formula) needs to be either 'occurrence', 'abundance', 'biomass' or 'quantity' data. The type of response data will help specify the type of error distribution to be used. The dependent variables (the right hind side) of this formula specifies the dependence of the species archetype probabilities on covariates. For all model the basic formula structure follows something like this: cbind(spp1,spp2,spp3)~1+temperature+rainfall
an object of class "formula" (or an object that can be
coerced to that class). The right hand side of this formula specifies the
dependence of the species"'" data on covariates (typically different
covariates to archetype_formula to avoid confusing confounding).
Current the formula is set at ~ 1 by default for species-specific intercepts
for the archetype models. If you include a species specific formula which has
more than an intercept you will be fitting a partial species archetype model
which has species specific covariates and archetype specific covariates.
a matrix or data.frame which contains the 'species_data' matrix, a const and the covariates in the structure of spp1, spp2, spp3, const, temperature, rainfall. dims of matrix should be nsites*(nspecies+const+covariates).
The number of archetypes (mixing components/groups) to estimate from the data. This need to be explicitly declared. By default it is set to 3.
The family of statistical family to use within the ecomix models. a choice between "bernoulli", "binomial", "poisson", "negative.binomial", "tweedie" and "gaussian" families are possible and applicable to specific types of data.
a numeric vector of length nrow(data) (n sites) that is included into the model as an offset. It is included into the conditional part of the model where conditioning is performed on the SAM.
a numeric vector of length ncol(Y) (n species) that is used as weights in the log-likelihood calculations. If NULL (default) then all weights are assumed to be identically 1. Because we are estimating the log-likelihood over species (rather than sites), the weights should be a vector n species long.
a numeric vector of n species long. This is used for undertaking a Bayesian Bootstrap. See 'vcov.species_mix' for more details.
The size of the sample for a binomial model (defaults to 1).
The power parameter for 'tweedie' family. Default is 1.6, and this is assigned to all species
a list of control parameters for optimization and calculation. See details.
NULL a numeric vector that provides approximate starting values for species_mix coefficients. These are family specific, but at a minimum you will need pi (additive_logistic transformed), alpha (intercepts) and beta (mixing coefs).
either a boolean or a vector of characters. If TRUE (default for species_mix(qv)), then some objects used in the estimation of the model"'"s parameters are returned in a list entitled "titbits" in the model object. Some functions, for example plot.species_mix(qv) and predict.species_mix(qv), will require some or all of these pieces of information. If titbits=FALSE (default for species_mix.multifit(qv)), then an empty list is returned. If a character vector, then just those objects are returned. Possible values are:"Y" for the outcome matrix, "X" for the model matrix for the SAM model, "offset" for the offset in the model, "site_spp_weights" for the model weights, "archetype_formula" for the formula for the SAMs, "species_formula" for the formula for the species-specific model, "control" for the control arguments used in model fitting, "family" for the conditional distribution of the species data. Care needs to be taken when using titbits=TRUE in species_mix.multifit(qv) calls as titbits is created for EACH OF THE MODEL FITS. If the data is large or if nstart is large, then setting titbits=TRUE may give users problems with memory.
for species_mix.multifit only. The number of random starts to perform for re-fitting. Default is 10, which will need increasing for serious use.
for species_mix.multifit only. The number of cores to spread the re-fitting over.
A species mix multifit object
Ignored
species_mix.multifit is used to fit species archetype models with multiple fits. This can be useful for complex models where there are many local maximum in the likelihood. The multiple restarts using in `species_mix.multifit` can help search the log-likelihood using starting values with some random variance (noise).
# \donttest{
library(ecomix)
set.seed(42)
sam_form <- stats::as.formula(paste0('cbind(',paste(paste0('spp',1:20),
collapse = ','),")~x1+x2"))
sp_form <- ~ 1
beta <- matrix(c(-2.9,-3.6,-0.9,1,.9,1.9),3,2,byrow=TRUE)
dat <- data.frame(y=rep(1,100),x1=stats::runif(100,0,2.5),
x2=stats::rnorm(100,0,2.5))
dat[,-1] <- scale(dat[,-1])
simulated_data <- species_mix.simulate(archetype_formula = sam_form,
species_formula = sp_form,
data = dat,beta=beta,family="bernoulli")
#> Random alpha from normal (-1,0.5) distribution
fm1 <- species_mix(archetype_formula = sam_form,species_formula = sp_form,
data = simulated_data, family = 'bernoulli', nArchetypes=3)
#> SAM modelling
#> There are 3 archetypes
#> There are 100 site observations for 20 species
#> The model for the archetype (grouping) is ~x1 + x2
#> The model for the species is ~1
#> You are implementing a bernoulli Species Archetype Model.
#> This model uses a logit link function.
#> Using ECM algorithm to find starting values; using 1 refits
#> ECM restart 1 of 1
#> Initialising starting values
#> Initial groups parameter estimates by K-means clustering
#> Iteration: 1 | New loglik -950.728 | Ratio loglik 0
#> Iteration: 2 | New loglik -829.843 | Ratio loglik 0.8728
#> Iteration: 3 | New loglik -820.246 | Ratio loglik 0.9884
#> Iteration: 4 | New loglik -819.839 | Ratio loglik 0.9995
#> initial value 819.837517
#> iter 10 value 819.706573
#> final value 819.703416
#> converged
fmods <- species_mix.multifit(archetype_formula = sam_form,
species_formula = sp_form, data=simulated_data, family = 'bernoulli',
nstart = 10, nArchetypes=3)
#> SAM modelling
#> Running species_mix.multifit using 10 starts to fit 3 Archetypes to the data.
#> There are 3 archetypes
#> There are 100 site observations for 20 species
#> The model for the archetype (grouping) is ~x1 + x2
#> The model for the species is ~1
#> You are implementing a bernoulli Species Archetype Model.
#> This model uses a logit link function.
#>
|
| | 0%.
|
|======== | 11%.
|
|================ | 22%.
|
|======================= | 33%.
|
|=============================== | 44%.
|
|======================================= | 56%.
|
|=============================================== | 67%.
|
|====================================================== | 78%.
|
|============================================================== | 89%.
|
|======================================================================| 100%.
# }
#Print information about a species_mix model
# \donttest{
print(fmods)
#> A multiple fit bernoulli species_mix model object
#>
#> You fitted 10 random starts using 3 Archetypes.
#>
#> The best model based on BIC was start 7
#> species_mix model
#>
#> Pi
#> Archetype1 Archetype2 Archetype3
#> 0.40 0.25 0.35
#>
#> Coefficients
#> $alpha
#> spp1 spp2 spp3 spp4 spp5 spp6 spp7
#> -0.8015555 -1.6402804 -0.5754542 -1.0843063 -1.0207942 -1.1778805 -0.7333331
#> spp8 spp9 spp10 spp11 spp12 spp13 spp14
#> -1.2033406 -1.6419186 -0.1858125 -0.8569556 -1.3912914 -0.7040136 -0.7040533
#> spp15 spp16 spp17 spp18 spp19 spp20
#> -0.8931633 -0.3781974 -1.3860213 -1.0552188 -1.2037696 -0.4297331
#>
#> $beta
#> x1 x2
#> Archetype1 -0.8300914 0.9550959
#> Archetype2 -2.9107038 -3.5526441
#> Archetype3 0.9916798 1.9873880
#>
#> You can access more elements of this model using x$multiple_fits[[ 7 ]]
# }