R/communitySimulH.formula.R

In HMSC: Hierarchical modelling of species community

communitySimulH.formula <-
function(formula,data,nsp,datatype="binomial",paramXDist,paramX=NULL,means=NULL,sigma=NULL,...){
### F. Guillaume Blanchet - Avril 2013 - June 2013
##########################################################################################
	### Make sure any variation of datatype works
	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")
		}
	}
	
	if(length(formula)==2){
		spNames<-"sp"
		
		#### Extract the model in the formula
		explanaFormula<-formula[[2L]]
	}else{
		spNames<-as.character(formula[[2L]])
		
		#### Extract the model in the formula
		explanaFormula<-formula[[3L]]
	}
	
	nsite<-nrow(data)
	#### Check if the names in the formula and are associated to variables in the data
	varNames<-colnames(data)
	if(is.null(varNames)){
		stop("variables in 'data' need to be named")
	}
	formulaNames<-all.vars(explanaFormula)
	
	varInFormula<-which(formulaNames%in%varNames)
	nvarInFormula<-length(varInFormula)
	
	paramXNULL<-is.null(paramX)
	if(paramXNULL){
		paramXInFormula<-which(!(formulaNames%in%varNames))
		nparamXInFormula<-length(paramXInFormula)
		paramXNames<-formulaNames[which(!(formulaNames%in%varNames))]
	}else{
		paramXNames<-colnames(paramX)
		paramXInFormula<-which(formulaNames%in%paramXNames)
		nparamXInFormula<-length(paramXInFormula)
	}
	
	if(paramXNULL){
		#### Generate 'means'
		if(is.null(means)){
			means<-rnorm(nparamXInFormula)
		}
		
		if(is.null(names(means))){
			names(means)<-paramXNames
		}
		
		#### Sigma matrix of the community
		if(is.null(sigma)){
			sigma<-as.matrix(rWishart(1,nparamXInFormula+1,diag(nparamXInFormula))[,,1])
		}else{
			if(!isSymmetric(sigma)){
				stop("'sigma' is not a symmetric matrix")
			}
		}
		
		if(is.null(colnames(sigma))){
			colnames(sigma)<-paramXNames
		}
		if(is.null(rownames(sigma))){
			rownames(sigma)<-paramXNames
		}
		
		#### paramX for each species
		paramX<-mvrnorm(nsp,mu=means,Sigma=sigma)
		colnames(paramX)<-paramXNames
		rownames(paramX)<-paste(spNames,1:nsp,sep="")
	}
	
	#==================================================
	#### Construct a list including data and paramXeters
	#==================================================
	nVarparamX<-length(formulaNames)
	listData<-vector("list",length=nVarparamX)
	
	names(listData)[1:nvarInFormula]<-formulaNames[varInFormula]
	names(listData)[(nvarInFormula+1):(nparamXInFormula+nvarInFormula)]<-formulaNames[paramXInFormula]
	
	counter<-1
	
	#### Data
	for(i in which(varNames%in%formulaNames[varInFormula])){
		listData[[counter]]<-data[,i]
		counter<-counter+1
	}
	
	#### paramXeters
	modelMat<-matrix(NA,nrow=nrow(data),ncol=nsp)
	
	for(h in 1:nsp){
		counterparamX<-counter
		for(i in 1:ncol(paramX)){
			listData[[counterparamX]]<-paramX[h,i]
			counterparamX<-counterparamX+1
		}
		modelMat[,h]<-eval(explanaFormula,listData)
	}
	
	#=================================
	### Model (Site by species matrix)
	#=================================
	if(datatype=="binomial"){
		probMat<-ilogit(modelMat)
	}
	if(datatype=="poisson"){
		probMat<-ilog(modelMat)
	}
	if(datatype=="nbinomial"){
		meanMat<-ilog(modelMat)
		probMat<-paramXDist/(paramXDist+meanMat)
	}
	
	colnames(probMat)<-paste(spNames,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(spNames,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)
				
				counterCheck<-1
				while(nsptochange!=0 & counterCheck <= 200){
					### Generate new paramXeters, probabilities and data
					for(h in sptochange){
						counterparamX<-counter
						paramX[h,]<-mvrnorm(1,mu=means,Sigma=sigma)
						for(i in 1:ncol(paramX)){
							listData[[counterparamX]]<-paramX[h,i]
							counterparamX<-counterparamX+1
						}
						modelMat[,h]<-eval(explanaFormula,listData)
						
						if(datatype=="poisson"){
							probMat[,h]<-ilog(modelMat[,h])
							suppressWarnings(dataMat[,h]<-rpois(nsite,probMat[,h]))
						}
						if(datatype=="nbinomial"){
							meanMat[,h]<-ilog(modelMat[,h])
							probMat[,h]<-paramXDist/(paramXDist+meanMat[,h])
							suppressWarnings(dataMat[,h]<-rnbinom(nsite,paramXDist,probMat[,h]))
						}
					}
					dataMatNA<-is.na(dataMat)
					sptochange<-unique(which(dataMatNA,arr.ind=TRUE)[,2])
					nsptochange<-length(sptochange)
					counterCheck<-counterCheck+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
	rownames(data)<-paste("site",1:nsite,sep="")
	
	data<-list(Y=dataMat,X=data[,which(colnames(data)%in%formulaNames[varInFormula])])
	attributes(data)<-list(names=c("Y","X"),Ypattern=spNames)
	class(data)<-"HMSCdata"
	
	#### Format paramXeters
	if(paramXNULL){
		allparamX<-list(paramX=paramX,means=means,sigma=sigma)
		class(allparamX)<-"HMSCparamX"
	}else{
		allparamX<-paramX
	}
	
	if(datatype=="binomial"){
		res<-list(data=data,param=allparamX,probMat=probMat)
	}
	if(datatype=="poisson"){
		res<-list(data=data,param=allparamX,lambdaMat=probMat)
	}
	if(datatype=="nbinomial"){
		res<-list(data=data,param=allparamX,probMat=probMat,paramXDist=paramXDist)
	}
	class(res)<-"communitySimul"
	
	return(res)
}