R/communitySimulH.default.R

In HMSC: Hierarchical modelling of species community

communitySimulH.default <-
function(X,nsp,Random=NULL,datatype="binomial",paramXDist,paramX=NULL,paramRandom=NULL,means=NULL,sigma=NULL,RandomVar=NULL,RandomCommSp=NULL,...){
### F. Guillaume Blanchet - Février 2013, Avril 2013, June 2013
##########################################################################################
	type<-c("binomial","poisson","nbinomial")
	datatype<-match.arg(datatype,type)
	
	### Basic check
	if(datatype=="nbinomial"){
		if(paramXDist==0){
			stop("'paramXDist' needs to be larger than 0")
		}
	}
	
	### X matrix
	nsite<-nrow(X)
	nexp<-ncol(X)
	
	if(is.null(colnames(X))){
		colnames(X)<-paste("x",1:nexp,sep="")
	}
	if(is.null(rownames(X))){
		rownames(X)<-paste("site",1:nsite,sep="")
	}
	
	#================
	### Random effect
	#================
	if(!is.null(Random)){
		if(!is.factor(Random)){
			stop("'Random' should be a factor")
		}
		if(length(Random)!=nsite){
			stop("'Random' should have a length equal to the number of rows of 'Y'")
		}
		### Transform Random to a 0-1 table
		nlev<-nlevels(Random)
		lev<-levels(Random)
		RandomMat<-matrix(0, nsite, nlev)
		RandomMat[(1:nsite)+nsite * (unclass(Random) - 1)]<-1
		colnames(RandomMat)<-lev
		rownames(RandomMat)<-paste("site",1:nsite,sep="")
		
		### Construct a covariance matrix for the random effect
		if(is.null(RandomVar)){
			RandomVar<-1
		}
		if(is.null(RandomCommSp)){
			RandomCommSp<-runif(1,0,RandomVar-0.0001)
		}
		
		if(RandomVar<RandomCommSp){
			stop("'RandomVar' should be larger than 'RandomCommSp'")
		}
		RandomCov<-matrix(RandomCommSp,nsp,nsp)
		diag(RandomCov)<-RandomVar
		
		### Construct random effect
		if(is.null(paramRandom)){
			paramRandom<-matrix(mvrnorm(n=nlev,mu=rep(0,nsp),Sigma=RandomCov),ncol=nlev)
			rownames(paramRandom)<-paste("sp",1:nsp,sep="")
			colnames(paramRandom)<-lev
		}else{
			if(is.null(rownames(paramRandom))){
				rownames(paramRandom)<-paste("sp",1:nsp,sep="")
			}
			if(is.null(colnames(paramRandom))){
				colnames(paramRandom)<-lev
			}
		}
		randomEffect<-tcrossprod(RandomMat,paramRandom)
	}
	
	paramXNULL<-is.null(paramX)
	### Mean of the community paramX
	if(is.null(paramX)){
		if(is.null(means)){
			means<-rnorm(nexp)
		}
		if(is.null(names(means))){
			names(means)<-paste("p",1:nexp,sep="")
		}
		
		### Sigma matrix of the community
		if(is.null(sigma)){
			sigma<-as.matrix(rWishart(1,nexp+1,diag(nexp))[,,1])
		}else{
			if(!isSymmetric(sigma)){
				stop("'sigma' is not a symmetric matrix")
			}
		}
		if(is.null(colnames(sigma))){
			colnames(sigma)<-paste("p",1:nexp,sep="")
		}
		if(is.null(rownames(sigma))){
			rownames(sigma)<-paste("p",1:nexp,sep="")
		}
		### paramX for each species
		paramX<-mvrnorm(nsp,mu=means,Sigma=sigma)
		colnames(paramX)<-paste("p",1:nexp,sep="")
		rownames(paramX)<-paste("sp",1:nsp,sep="")
	}
	
	#=================================
	### Model (Site by species matrix)
	#=================================
	LinearModel<-tcrossprod(X,paramX)
	
	if(!is.null(Random)){
		LinearModel<-LinearModel+randomEffect
	}
	
	if(datatype=="binomial"){
		probMat<-ilogit(LinearModel)
	}
	if(datatype=="poisson"){
		probMat<-ilog(LinearModel)
	}
	if(datatype=="nbinomial"){
		meanMat<-ilog(LinearModel)
		probMat<-paramXDist/(paramXDist+meanMat)
	}
	
	colnames(probMat)<-paste("sp",1:nsp,sep="")
	rownames(probMat)<-paste("site",1:nsite,sep="")
	
	#===================================
	### Sampled (Site by species matrix)
	#===================================
	dataMat<-matrix(NA,nsite,nsp)
	
	if(datatype=="binomial"){
		for(i in 1:nsite){
			dataMat[i,]<-rbinom(nsp,1,probMat[i,])
		}
	}
	
	if(datatype=="poisson"){
		for(i in 1:nsite){
			suppressWarnings(dataMat[i,]<-rpois(nsp,probMat[i,]))
		}
	}
	if(datatype=="nbinomial"){
		for(i in 1:nsp){
			suppressWarnings(dataMat[,i]<-rnbinom(nsite,paramXDist,probMat[,i]))
		}
	}
	
	colnames(dataMat)<-paste("sp",1:nsp,sep="")
	rownames(dataMat)<-paste("site",1:nsite,sep="")
	
	### Check if NAs were generated
	if(datatype=="poisson" | datatype=="nbinomial"){
		if(paramXNULL){
			dataMatNA<-is.na(dataMat)
			### Replace species generated with NAs
			if(any(dataMatNA)){
				sptochange<-unique(which(dataMatNA,arr.ind=TRUE)[,2])
				nsptochange<-length(sptochange)
				
				counter<-1
				while(nsptochange!=0 & counter <= 200){
					### Generate new paramXeters, probabilities and data
					paramXtochange<-mvrnorm(nsptochange,mu=means,Sigma=sigma)
					LinearModeltochange<-tcrossprod(X,matrix(paramXtochange,ncol=ncol(X)))
					if(datatype=="poisson"){
						probMattochange<-ilog(LinearModeltochange)
						for(i in 1:nsptochange){
							suppressWarnings(dataMat[,sptochange[i]]<-rpois(nsite,probMattochange[,i]))
						}
					}
					if(datatype=="nbinomial"){
						meanMattochange<-ilog(LinearModeltochange)
						probMattochange<-paramXDist/(paramXDist+meanMattochange)
						for(i in 1:nsptochange){
							suppressWarnings(dataMat[,sptochange[i]]<-rnbinom(nsite,paramXDist,probMattochange[,i]))
						}
					}
					paramX[sptochange,]<-paramXtochange
					probMat[sptochange,]
					dataMatNA<-is.na(dataMat)
					sptochange<-unique(which(dataMatNA,arr.ind=TRUE)[,2])
					nsptochange<-length(sptochange)
					counter<-counter+1
				}
				if(counter > 200){
					if(datatype=="poisson"){
						stop("A community could not be simulated because it the model always yield lambda too large for rpois(), change the model")
					}
					if(datatype=="nbinomial"){
						stop("A community could not be simulated because it the model always yield means too large for rnbinomial(), change the model or paramXDist")
					}
				}
			}
		}
	}
	
	#### Format data
	if(is.null(Random)){
		data<-list(Y=dataMat,X=X)
		attributes(data)<-list(names=c("Y","X"),Ypattern="sp")
		class(data)<-"HMSCdata"
	}else{
		data<-list(Y=dataMat,X=X,Random=RandomMat)
		attributes(data)<-list(names=c("Y","X","Random"),Ypattern="sp")
		class(data)<-c("HMSCdata","HMSCdataRandom")
	}
	#### Format parameters
	if(paramXNULL){
		if(is.null(Random)){
			allparam<-list(paramX=paramX,means=means,sigma=sigma)
		}else{
			allparam<-list(paramX=paramX,paramRandom=paramRandom,means=means,sigma=sigma,RandomVar=RandomVar,RandomCommSp=RandomCommSp)
		}
		class(allparam)<-"HMSCparam"
	}else{
		allparam<-paramX
	}
	
	if(datatype=="binomial"){
		res<-list(data=data,param=allparam,probMat=probMat)
	}
	if(datatype=="poisson"){
		res<-list(data=data,param=allparam,lambdaMat=probMat)
	}
	if(datatype=="nbinomial"){
		res<-list(data=data,param=allparam,probMat=probMat,paramXDist=paramXDist)
	}
	class(res)<-"communitySimul"
	
	return(res)
}