communitySimul: Simulate a species community
In HMSC: Hierarchical modelling of species community

Description Usage Arguments Details Value Note Author(s) Examples

Function communitySimulH simulates a species community on a transect assuming that the species stem from linear or non-linear (when using formula) relationships with the proposed descriptors. Function communitySimulHT also simulates a species community but it assumes that species stem from linear relationships with descriptors and species traits. A plot of the probability of occurrence of the model and of the simulated data can also be drawn with the function plot.

## Default S3 method:
communitySimulH(X, nsp, Random=NULL, datatype = "binomial",paramXDist,paramX=NULL,paramRandom=NULL, means = NULL, sigma = NULL,RandomVar=NULL,RandomCommSp=NULL,...)
## S3 method for class 'formula'
communitySimulH(formula,data,nsp,datatype="binomial",paramXDist,paramX=NULL,means=NULL,sigma=NULL,...)

## Default S3 method:
communitySimulHT(X, Tr, Random=NULL, datatype = "binomial",paramXDist,paramX=NULL,paramTr=NULL,paramRandom=NULL, sigma = NULL,RandomVar=NULL,RandomCommSp=NULL,...)

## S3 method for class 'communitySimul'
plot(x, ...)

`formula`	An object of class formula describing the model to be fitted. Variables and parameters need to be specified. (See details)
`X`	Matrix of descriptors defining the species. Rows are sites and the columns are descriptors.
`Tr`	Matrix defining species traits. Rows are traits and columns are species.
`Random`	A factor defining a random effect on the samples.
`data`	A data.frame containing the variables in the model.
`x`	Object of class communitySimul.
`nsp`	Numeric. Number of species to be simulated in the community.
`paramX`	A matrix of model parameters defining how each species (rows) is characterized by each descriptors (columns).
`paramTr`	A matrix of model parameters defining how each descriptors (rows) characterizes species traits (columns).
`paramRandom`	A matrix of model parameters defining how each level of the random effect (rows) characterizes species (columns).
`means`	Vector of means used to generate the model parameters. (See details)
`sigma`	Covariance matrix used to generate the model parameters. (See details)
`RandomVar`	Scalar defining the within group variance in `Random`. (See details)
`RandomCommSp`	A scalar defining how the random effect influence species independently (0) and on as a community (1). (See details)
`datatype`	Character string defining the data type in the simulated community. Any unambiguous variation of wording used in this argument is accepted. (See details)
`paramXDist`	A value that defines the free parameter of distributions (i.e. shape, size, range) from which species are sampled. (See details)
`...`	Parameters passed to other functions.

The left hand-side of formula (the response) is not necessary. If it is included, the string of character that defines the response variable will be used as the common character in the species's name. Also, since the parameters and the variables are required when defining the formula, the function can generate data for linear as well as non-linear regression models.

The descriptors in X are used without any modifications (or additions) to simulate the species community. As such, a column of 1 should be included in X for the model used to construct the community to include an intercept. As is the case for data.

The values in means and sigma are used as parameter of a multivariate normal distribution to generate the model's parameters (mvrnorm is used in the function). When the paramX is set to NULL, the means and the sigma will be randomly generated if they are also set to NULL. When values are given to paramX the values of the means and the sigma are not used and if either is set to NULL, no data will be generated for either set of parameter. When generated, the values for means are randomly sampled from rnorm and the values for sigma are randomly sampled from a Wishart distribution (rWishart).

When using the communitySimulHT, means are not necessary because they are calculated from the species traits, that is from Tr and paramTr

The argument paramXDist is only active when datatype = "nbinomial".

If the arguments RandomVar is NULL, the value given to RandomVar is 1. As for the argument RandomCommSp, if it is NULL, the value given to RandomCommSp is sampled from a Uniform distribution runif with a minimum of 0 and a maximum of RandomVar-0.001.

The communitySimul functions can be used for prediction if the parameters of a model (paramX and paramTr) are given with the proper variables given.

Currently, this function simulates different types of community:

`binomial`	For presence/absence of species (using a logit link).
`poisson`	For count data. The current implementation checks for NAs in the simulated data if parameters were not supplied to the function. The function will try 200 times to remove NAs by generating a different set of paramters. If NAs are still produced, an error message will be sent. NAs are generate because the parameters randomly selected in `paramX` yield a lambda so large that `rpois` cannot make a proper sample.
`nbinomial`	For count data. The current implementation checks for NAs as for `poisson`. `paramXDist` is used to define the dispersion parameter (`size` argument in `rnbinom`) of the negative binomial distribution, for this reason it must be strictly positive. This values will be returned by the functions.

The functions communitySimulH and communitySimulHT return an object of class communitySimul with the following components:

`data`	an object of class HMSCdata
`param`	an object of class HMSCparam
`probMat`	Matrix. If `datatype` is `binomial`, `probMat` defines the occurrence probability of each species for each site. If `datatype` is `nbinomial`, `probMat` defines the probability of finding a species with a predefine dispersion parameter.
`lambdaMat`	A matrix defining the `lambda` of the Poisson distribution.

Currently, the plot function is only valid for occurrence data simulated on a transect.

F. Guillaume Blanchet

### Construct some descriptors
#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
### Simulate community using the linear algorithm
#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
desc<-cbind(1,scale(1:50),scale(1:50)^2)

#=======================================
### Simulate communities with 20 species
#=======================================
### Using only a single descriptors
comm1<-communitySimulH(matrix(desc[,1],ncol=1),nsp=20)
plot(comm1)

comm2<-communitySimulH(matrix(desc[,2],ncol=1),nsp=20)
plot(comm2)

comm3<-communitySimulH(matrix(desc[,3],ncol=1),nsp=20)
plot(comm3)

#------------------------
### Using all descriptors
#------------------------
comm<-communitySimulH(desc,nsp=20)
plot(comm)

### Defining the means and sigma
mn<-rnorm(3)
covMat<-rWishart(1,4,diag(3))[,,1]

commDes1<-communitySimulH(desc,nsp=20,means=mn,sigma=covMat)
commDes2<-communitySimulH(desc,nsp=20,means=mn,sigma=covMat)
commDes3<-communitySimulH(desc,nsp=20,means=mn,sigma=covMat)

dev.new()
plot(commDes1)
dev.new()
plot(commDes2)
dev.new()
plot(commDes3)

#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
### Simulate community using the formula algorithm
#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
desc<-cbind(1,scale(1:50),scale(1:50)^2)
colnames(desc)<-paste("x",1:3,sep="")
formule <- ~p1*x1+x2/p2+exp(p3*x3)

comm<-communitySimulH(formule,data=desc,nsp=20)
plot(comm)

#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
### Simulate community using descriptors and traits
#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
desc <- cbind(1,scale(1:50),scale(1:50)^2)
trait <- t(cbind(1,scale(1:20)))
commTrait<-communitySimulHT(desc,trait)
plot(commTrait)

#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
### Simulate community with a random effect
#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
desc <- cbind(1,scale(1:50),scale(1:50)^2)
arbi<-as.factor(rep(1:2,each=25))
commRandomEff<-communitySimulH(desc,Random=arbi,nsp=20)
plot(commRandomEff)

#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
### Simulate community with a random effect where traits are considered
#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
desc<-cbind(1,scale(1:50),scale(1:50)^2)
trait <- t(cbind(1,scale(1:10)*2))
arbi<-as.factor(c(rep(1:2,each=16),rep(3,18)))
commTraitRandomEff<-communitySimulHT(desc,trait,Random=arbi)
plot(commTraitRandomEff)