R/communitySimulHT.default.R

In HMSC: Hierarchical modelling of species community

communitySimulHT.default <-
function(X,Tr,Random=NULL,datatype="binomial",paramXDist,paramX=NULL,paramTr=NULL,paramRandom=NULL,sigma=NULL,RandomVar=NULL,RandomCommSp=NULL,...){
### F. Guillaume Blanchet - July 2013
##########################################################################################
	### Make sure X is a matrix
	X<-as.matrix(X)
	
	### Make sure Tr is a matrix
	if(is.null(dim(Tr))){
		Tr<-t(as.matrix(Tr))
	}else{
		Tr<-as.matrix(Tr)
	}
	
	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")
		}
	}
	
	### matrices dimensions
	nsite<-nrow(X)
	nexp<-ncol(X)
	nsp<-ncol(Tr)
	ntr<-nrow(Tr)
	
	if(is.null(colnames(X))){
		colnames(X)<-paste("x",1:nexp,sep="")
	}
	if(is.null(rownames(X))){
		rownames(X)<-paste("site",1:nsite,sep="")
	}
	if(is.null(colnames(Tr))){
		colnames(Tr)<-paste("sp",1:nsp,sep="")
	}
	if(is.null(rownames(Tr))){
		rownames(Tr)<-paste("t",1:ntr,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<-mvrnorm(n=nlev,mu=rep(0,nsp),Sigma=RandomCov)
			colnames(paramRandom)<-paste("sp",1:nsp,sep="")
			rownames(paramRandom)<-lev
			paramRandom<-t(paramRandom)
		}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)
	}
	
	### Trait parameters
	if(is.null(paramTr)){
		paramTr<-matrix(rnorm(ntr*nexp),nrow=nexp,ncol=ntr)
		rownames(paramTr)<-paste("p",1:nexp,sep="")
		colnames(paramTr)<-paste("t",1:ntr,sep="")
	}
	
	### Community Mean (projector of traits and parameters)
	meanSp<-matrix(NA,ncol=nsp,nrow=nexp)
	for(i in 1:nsp){
		meanSp[,i]<-tcrossprod(Tr[,i],paramTr)
	}
	colnames(meanSp)<-paste("sp",1:nsp,sep="")
	rownames(meanSp)<-paste("p",1:nexp,sep="")

	paramXNULL<-is.null(paramX)
	### Mean of the community paramXeters
	if(paramXNULL){
		### 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<-matrix(NA,nrow=nsp,ncol=nexp)
		### paramX for each species
		for(i in 1:nsp){
			paramX[i,]<-mvrnorm(1,mu=meanSp[,i],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
					if(datatype=="poisson"){
						for(i in 1:nsptochange){
							paramXtochange<-mvrnorm(nsptochange,mu=meanSp[,sptochange[i]],Sigma=sigma)
							LinearModeltochange<-tcrossprod(X,matrix(paramXtochange,ncol=ncol(X)))
							probMattochange<-ilog(LinearModeltochange)
							suppressWarnings(dataMat[,sptochange[i]]<-rpois(nsite,probMattochange))
						}
					}
					if(datatype=="nbinomial"){
						for(i in 1:nsptochange){
							paramXtochange<-mvrnorm(nsptochange,mu=meanSp[,sptochange[i]],Sigma=sigma)
							LinearModeltochange<-tcrossprod(X,matrix(paramXtochange,ncol=ncol(X)))
							meanMattochange<-ilog(LinearModeltochange)
							probMattochange<-paramXDist/(paramXDist+meanMattochange)
							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 the model always yield lambda too large for rpois(), change the model")
					}
					if(datatype=="nbinomial"){
						stop("A community could not be simulated because 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,Tr=Tr)
		attributes(data)<-list(names=c("Y","X","Tr"),Ypattern="sp")
		class(data)<-c("HMSCdata","HMSCdataTrait")
	}else{
		data<-list(Y=dataMat,X=X,Tr=Tr,Random=RandomMat)
		attributes(data)<-list(names=c("Y","X","Tr","Random"),Ypattern="sp")
		class(data)<-c("HMSCdata","HMSCdataTrait","HMSCdataRandom")
	}
	#### Format parameters
	if(paramXNULL){
		if(is.null(Random)){
			allparam<-list(paramX=paramX,paramTr=paramTr,means=meanSp,sigma=sigma)
		}else{
			allparam<-list(paramX=paramX,paramTr=paramTr,paramRandom=paramRandom,means=meanSp,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)
}