R/species_mix-class.R
species_mix.RdFits a mixture-of-regressions (SAM) to identify species archetypes.
species_mix(
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(),
inits = NULL,
titbits = TRUE
)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 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 model, family can be a function from family, such as binomial(). Or it can be a character string of one of the following: "bernoulli", "binomial", "poisson","negative.binomial", "tweedie" and "gaussian". An error will be thrown if the family is not recognized. The family you use is generally specific to the type of data you have. Users can also specify a link function within the family object, eg: binomial(link="logit"). Currently, ecomix can implement "logit", "cloglog", "log" and "identity" link functions.
a numeric vector of length nrow(data) (n sites) that represents a known amount of effort for models with a log or cloglog link function. 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. These are weights used on the species log-likelihoods. We tend to use these for doing the Bayesian Bootstrap, where we pass a vector S long of Dirichlet weights. 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.
species_mix is used to fit mixtures of glms to multivariate species data. The function uses BFGS to optimize the mixture likelihood. There is the option to use ECM algorithm to get appropriate starting parameters. `species_mix` acts as a wrapper for species_mix.fit that allows for easier data input. The data frames are merged into the appropriate format for the use in species_mix.fit. Minima is found using vmmin (BFGS). Currently 'bernoulli', 'binomial', 'poisson', 'negative.binomial', 'tweedie' and 'gaussian' distributions can be fitted using the species_mix.
species_mix is used to fit mixtures of glms to multivariate species data. The function uses BFGS to optimize the mixture likelihood. There is the option to use ECM algorithm to get appropriate starting parameters. `species_mix` acts as a wrapper for species_mix.fit that allows for easier data input. The data frames are merged into the appropriate format for the use in species_mix.fit. Minima is found using vmmin (BFGS).
Control arguments for species_mix are as follows:
The maximum number of iterations. Default is 1000.
Should any reporting be performed? Default is FALSE, for reporting.
Non-negative integer. If positive, tracing information on the progress of the optimization is produced. Higher values may produce more tracing information. If quiet=TRUE, then trace will be set to zero.
Use glmnet to get initial starting values, default is TRUE, if FALSE, glm2 will be used.
Which initialsation approach to use for covariates, default is "random2" which is uses k-means on the groups and adds some random noise. Alternatives are 'kmeans' and 'kmed'
The amount of gaussian variation to add to starting values, default is NA, which adds variantion in based on the standard deviation of the parameters. Otherwise should be a double like 0.3
The number of presences required for each species to be included when generaling the starting values, the default is to remove all species <= 30 occurrences during inititalisation, these species will still be included in the model fit.
Use ECM to generate starting values? If FALSE, the model will use starting values without the ECM.
The number of ECM refits to run, the default is 2, users could increase this to potentially get better starting values.
The numer of ECM steps/iterations to run per ECM fit, the default is 5.
The relative tolerance for the ECM, default is quiet large 1e-2, makin this smaller will make the convergance tolerance stickter, you'll have to also adjust the "ecm.steps", to be larger if you make this stickter.
Print the intialisation values being passed to c++. Handy for debugging especially if 'vmmin' error call
The frequency of reports for optimisation. Default is 10 a report for 10th iteration.
The abstol for vmmin, the default is -Inf. See optim for further details
Relative convergence tolerance. The algorithm stops if it is unable to reduce the value by a factor of reltol * (abs(val) + reltol) at a step. Defaults to sqrt(.Machine$double.eps), typically about 1e-8.
Should optimisation for estimation occur? If TRUE (default) optimisation will occur. If FALSE no optimisation is performed.
Should we return the scores (derivates) when doing optimsation? If FALSE (default) scores will not be returned. If TRUE scores will be returned.
Should the log-likelihood be caulcated? If TRUE (default) then log-likelihood is calculated and returned. If FALSE then the log-likelihood is not calculated for return.
Should the scores be evaluated at the (final) parameter values. If TRUE (default) then they are calculated. If FALSE then they are not calculated.
A boolean 0 or 1. Default is 0 (false) and no penalities will be applied within the optimsation.
A numeric scalar. This is the penalty for the mixing coefs. The penality is an from a Dirichlet distribution.
A numeric scalar. This is the penalty for the alpha parameters in the species model. They are assumed to come from a normal distribution with standard deviation given as this parameter (default is 10).
A numeric scalar. This is the penalty for the beta parameters in the group model. They are assumed to come from a normal distribution with standard deviation given as this parameter (default is 10).
A numeric scalar. This is the penalty for the gamma parameters in the species model. They are assumed to come from a normal distribution with standard deviation given as this parameter (default is 10).
# \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.705982
#> final value 819.703422
#> converged
# }