R/hmscH.formula.R

In HMSC: Hierarchical modelling of species community

hmscH.formula <-
function(formula,data,start,priors=NULL,niter=12000,nburn=2000,nrot=100,thin=1,likelihoodtype="logit",details=TRUE,verbose=TRUE,...){
#### F. Guillaume Blanchet - Janvier 2013, March 2013, April 2013, June 2013
##########################################################################################
	call<-match.call()
	
	### Basic objects
	nparamX<-ncol(start$paramX)
	nsp<-nrow(start$paramX)
	nmeans<-length(start$means)
	
	### General checks
	if(niter < nburn){
		stop("'niter' should be equal or larger than 'burning'")
	}
	
	if(is.null(priors)){
		priors<-list(means0=rep(0,length(start$mean)),
					 kappa0=0,
					 nu0=rep(0,nparamX),
					 Lambda0=diag(nparamX))
	}
	
	#=====================================================================================================
	#### Construct the model using Single Component Adaptative Metropolis (SCAM) with eigenvector rotation
	#=====================================================================================================
	#--------------------
	#### Starting objects
	#--------------------
	#### Rotation objects (to use when no rotation is carried out)
	rotRes<-vector("list",length=4)
	names(rotRes)<-c("values","vectors","wParamX","wOuterParamX")
	rotRes[[1]]<-matrix(1,nrow=nsp,ncol=nparamX)
	rotRes[[2]]<-array(diag(nparamX),dim=c(nparamX,nparamX,nsp))
	rotRes[[3]]<-matrix(0,nrow=nsp,ncol=nparamX)
	rotRes[[4]]<-array(0,dim=c(nparamX,nparamX,nsp))

	accept<-array(0,dim=c(nsp,nparamX,2))
	sdvalue<-matrix(1,nrow=nsp,ncol=nparamX)
	scamRes<-list(start=start,rotation=rotRes,sdvalue=sdvalue,accept=accept)

	#### Store burning results
	burningObj<-vector("list",length=3)
	names(burningObj)<-c("paramX","means","sigma")
	burningObj[[1]]<-array(dim=c(nsp,nparamX,nburn))
	burningObj[[2]]<-matrix(NA,nrow=nburn,ncol=nmeans)
	burningObj[[3]]<-array(dim=c(nparamX,nparamX,nburn))
	
	#### Store detailed information
	if(details){
		detailsObj<-vector("list",length=2)
		names(detailsObj)<-c("sdvalue","acceptRatio")
		detailsObj[[1]]<-array(dim=c(nsp,nparamX,nburn))
		detailsObj[[1]][,,1]<-sdvalue
		detailsObj[[2]]<-array(dim=c(nsp,nparamX,nburn))
		detailsObj[[2]][,,1]<-matrix(0,nrow=nsp,ncol=nparamX)
	}
	
	#### Store results after burning
	modelres<-vector("list",length=5)
	names(modelres)<-c("paramX","means","sigma","sdvalue","acceptRatio")
	nvalues<-floor((niter-nburn)/thin)
	modelres[[1]]<-array(dim=c(nsp,nparamX,nvalues))
	modelres[[2]]<-matrix(NA,nrow=nvalues,ncol=nmeans)
	modelres[[3]]<-array(dim=c(nparamX,nparamX,nvalues))
	
	### Generate print counter
	if(verbose){
		printCounter<-round(seq(niter/5,niter,length=5))
	}
	#==================================
	#### SCAM with eigenvector rotation
	#==================================
	itsave<-1
	for(i in 1:niter){
		scamRes<-scamH(formula,data,scamRes$start,priors,scamRes$rotation$values,scamRes$rotation$vectors,scamRes$rotation$wParamX,scamRes$rotation$wOuterParamX,scamRes$sdvalue,scamRes$accept,likelihoodtype=likelihoodtype,burning=(i < nburn),rotation=(i > nrot & i < nburn),iter=i) # Pourrait changer
		#-------------------------
		#### Store burning results
		#-------------------------
		if(i < nburn){
			burningObj[[1]][,,i]<-scamRes$start$paramX
			burningObj[[2]][i,]<-scamRes$start$means
			burningObj[[3]][,,i]<-scamRes$start$sigma
			
			if(details){
				detailsObj[[1]][,,i+1]<-scamRes$sdvalue
				detailsObj[[2]][,,i+1]<-scamRes$accept[,,1]/scamRes$accept[,,2]
			}
		}else{
			it<-i-nburn
			if(it==thin*itsave){
				modelres[[1]][,,itsave]<-scamRes$start$paramX
				modelres[[2]][itsave,]<-scamRes$start$means
				modelres[[3]][,,itsave]<-scamRes$start$sigma
				itsave<-itsave+1
			}
		}
		
		### Print number of iterations completed
		if(verbose){
			if(any(i==printCounter)){
				print(i)
			}
		}
	}
	
	modelres[[4]]<-scamRes$sdvalue
	modelres[[5]]<-scamRes$accept[,,1]/scamRes$accept[,,2]
	
	#===================
	#### Name the pieces
	#===================
	### paramX
	dimnames(modelres[[1]])[[1]]<-rownames(start$paramX)
	dimnames(modelres[[1]])[[2]]<-colnames(start$paramX)

	### means
	dimnames(modelres[[2]])[[2]]<-colnames(start$paramX)
	
	### sigma
	dimnames(modelres[[3]])[[1]]<-colnames(start$paramX)
	dimnames(modelres[[3]])[[2]]<-colnames(start$paramX)
	
	#==========
	### Results
	#==========
	if(details){
		res<-list(call=call,formula=formula,model=modelres,burning=burningObj,details=detailsObj,data=data)
	}else{
		res<-list(call=call,formula=formula,model=modelres,burning=burningObj,data=data)
	}
	
	attr(res,"likelihood")<-likelihoodtype
	class(res)<-c("hmsc","HMSCformula")
	return(res)
}