R/hmscHT.HMSCdata.R

In HMSC: Hierarchical modelling of species community

hmscHT.HMSCdata <-
function(data,start,priors=NULL,niter=12000,nburn=2000,nrot=100,thin=1,likelihoodtype="logit",details=TRUE,verbose=TRUE,...){
#### F. Guillaume Blanchet - July 2013
##########################################################################################
	call<-match.call()
	
	### Basic objects
	nparamX<-ncol(start$paramX)
	nsp<-nrow(start$paramX)
	nparamTr<-ncol(start$paramTr)
	nexp<-nrow(start$paramTr)
	
	if(is.null(start$means)){
		start$means<-matrix(NA,ncol=nsp,nrow=nexp)
		for(i in 1:nsp){
			start$means[,i]<-tcrossprod(data$Tr[,i],start$paramTr)
		}
	}
	
	### General checks
	if(niter < nburn){
		stop("'niter' should be equal or larger than 'burning'")
	}
	
	if(is.null(priors)){
		if(!any(class(data)=="HMSCdataRandom")){
			priors<-list(means0=rep(0,nparamX),kappa0=0,nu0=rep(0,nparamX),Lambda0=diag(nparamX))
		}else{
			priors<-list(means0=rep(0,nparamX),kappa0=0,nu0=rep(0,nparamX),Lambda0=diag(nparamX),priorRandom=c(-1,1))
			nparamRandom<-ncol(start$paramRandom)
		}
	}
	
	#=====================================================================================================
	#### Construct the model using Single Component Adaptative Metropolis (SCAM) with eigenvector rotation
	#=====================================================================================================
	#--------------------
	#### Starting objects
	#--------------------
	#### Rotation objects (to use when no rotation is carried out)
	if(!any(class(data)=="HMSCdataRandom")){
		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))
	}else{
		rotRes<-vector("list",length=9)
		names(rotRes)<-c("values","vectors","wParamX","wOuterParamX","valuesRandom","vectorsRandom","wParamRandom","wOuterParamRandom","wRandomCov")
		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))

		rotRes[[5]]<-matrix(1,nrow=nsp,ncol=nparamRandom)
		rotRes[[6]]<-array(diag(nparamRandom),dim=c(nparamRandom,nparamRandom,nsp))
		rotRes[[7]]<-matrix(0,nrow=nsp,ncol=nparamRandom)
		rotRes[[8]]<-array(0,dim=c(nparamRandom,nparamRandom,nsp))

		rotRes[[9]]<-1
	}
	
	if(!any(class(data)=="HMSCdataRandom")){
		accept<-array(0,dim=c(nsp,nparamX,2))
		sdvalue<-matrix(1,nrow=nsp,ncol=nparamX)
	}else{
		acceptparamX<-array(0,dim=c(nsp,nparamX,2))
		acceptparamRandom<-array(0,dim=c(nsp,nparamRandom,2))
		acceptRandomCov<-c(0,0)
		accept<-list(accept=acceptparamX,acceptRandom=acceptparamRandom,acceptRandomCov=acceptRandomCov)
		
		sdvalueParamX<-matrix(1,nrow=nsp,ncol=nparamX)
		sdvalueParamRandom<-matrix(1,nrow=nsp,ncol=nparamRandom)
		sdvalueRandomCov<-1
		sdvalue<-list(sdvalue=sdvalueParamX,sdvalueRandom=sdvalueParamRandom,sdvalueRandomCov=sdvalueRandomCov)
	}
	
	scamRes<-list(start=start,rotation=rotRes,sdvalue=sdvalue,accept=accept)
	
	#### Store burning results
	if(!any(class(data)=="HMSCdataRandom")){
		burningObj<-vector("list",length=4)
		names(burningObj)<-c("paramX","paramTr","means","sigma")
		burningObj[[1]]<-array(dim=c(nsp,nparamX,nburn))
		burningObj[[2]]<-array(dim=c(nparamX,nparamTr,nburn))
		burningObj[[3]]<-array(dim=c(nparamX,nsp,nburn))
		burningObj[[4]]<-array(dim=c(nparamX,nparamX,nburn))
	}else{
		burningObj<-vector("list",length=6)
		names(burningObj)<-c("paramX","paramTr","means","sigma","paramRandom","RandomVar")
		burningObj[[1]]<-array(dim=c(nsp,nparamX,nburn))
		burningObj[[2]]<-array(dim=c(nparamX,nparamTr,nburn))
		burningObj[[3]]<-array(dim=c(nparamX,nsp,nburn))
		burningObj[[4]]<-array(dim=c(nparamX,nparamX,nburn))
		burningObj[[5]]<-array(dim=c(nsp,nparamRandom,nburn))
		burningObj[[6]]<-vector("numeric",length=nburn)
	}
	#### Store detailed information
	if(details){
		if(!any(class(data)=="HMSCdataRandom")){
			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)
		}else{
			detailsObj<-vector("list",length=6)
			names(detailsObj)<-c("sdvalue","acceptRatio","acceptRatioRandomCov","acceptRatioRandom","sdvalueRandom","sdvalueRandomCov")
			detailsObj[[1]]<-array(dim=c(nsp,nparamX,nburn))
			detailsObj[[1]][,,1]<-sdvalueParamX
			detailsObj[[2]]<-array(dim=c(nsp,nparamX,nburn))
			detailsObj[[2]][,,1]<-matrix(0,nrow=nsp,ncol=nparamX)
			detailsObj[[3]]<-numeric()
			detailsObj[[3]][1]<-0
			detailsObj[[4]]<-array(dim=c(nsp,nparamRandom,nburn))
			detailsObj[[4]][,,1]<-matrix(0,nrow=nsp,ncol=nparamRandom)
			detailsObj[[5]]<-array(dim=c(nsp,nparamRandom,nburn))
			detailsObj[[5]][,,1]<-sdvalueParamRandom
			detailsObj[[6]]<-numeric()
			detailsObj[[6]][1]<-0
		}
	}
	
	#### Store results after burning
	if(!any(class(data)=="HMSCdataRandom")){
		modelres<-vector("list",length=6)
		names(modelres)<-c("paramX","paramTr","means","sigma","sdvalue","acceptRatio")
		nvalues<-floor((niter-nburn)/thin)
		modelres[[1]]<-array(dim=c(nsp,nparamX,nvalues))
		modelres[[2]]<-array(dim=c(nparamX,nparamTr,nvalues))
		modelres[[3]]<-array(dim=c(nparamX,nsp,nvalues))
		modelres[[4]]<-array(dim=c(nparamX,nparamX,nvalues))
	}else{
		modelres<-vector("list",length=11)
		names(modelres)<-c("paramX","paramTr","means","sigma","sdvalue","acceptRatio","paramRandom","sdvalueRandomCov","sdvalueRandom","acceptRatioRandomCov","acceptRatioRandom")
		nvalues<-floor((niter-nburn)/thin)
		modelres$paramX<-array(dim=c(nsp,nparamX,nvalues))
		modelres$paramTr<-array(dim=c(nparamX,nparamTr,nvalues))
		modelres$means<-array(dim=c(nparamX,nsp,nvalues))
		modelres$sigma<-array(dim=c(nparamX,nparamX,nvalues))
		modelres$paramRandom<-array(dim=c(nsp,nparamRandom,nvalues))
	}
	
	### Generate print counter
	if(verbose){
		printCounter<-round(seq(niter/5,niter,length=5))
	}
	#==================================
	#### SCAM with eigenvector rotation
	#==================================
	itsave<-1
	if(!any(class(data)=="HMSCdataRandom")){
		for(i in 1:niter){
			scamRes<-scamHT(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$paramTr
				burningObj[[3]][,,i]<-scamRes$start$means
				burningObj[[4]][,,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$paramTr
					modelres[[3]][,,itsave]<-scamRes$start$means
					modelres[[4]][,,itsave]<-scamRes$start$sigma
					itsave<-itsave+1
				}
			}
			
			### Print number of iterations completed
			if(verbose){
				if(any(i==printCounter)){
					print(i)
				}
			}
		}
	}else{
		for(i in 1:niter){
			#---------------
			### Fixed effect
			#---------------
			scamResHT<-scamHT(data,scamRes$start,priors,scamRes$rotation$values,scamRes$rotation$vectors,scamRes$rotation$wParamX,scamRes$rotation$wOuterParamX,scamRes$sdvalue$sdvalue,scamRes$accept$accept,likelihoodtype=likelihoodtype,burning=(i < nburn),rotation=(i > nrot & i < nburn),iter=i) # Pourrait changer 
			
			### Assign the new estimated pieces of scamResH to scamRes
			scamRes$start$paramX<-scamResHT$start$paramX
			scamRes$start$sigma<-scamResHT$start$sigma
			scamRes$start$means<-scamResHT$start$means
			scamRes$start$paramTr<-scamResHT$start$paramTr
			
			scamRes$rotation$values<-scamResHT$rotation$values
			scamRes$rotation$vectors<-scamResHT$rotation$vectors
			scamRes$rotation$wParamX<-scamResHT$rotation$wParamX
			scamRes$rotation$wOuterParamX<-scamResHT$rotation$wOuterParamX
			
			scamRes$sdvalue$sdvalue<-scamResHT$sdvalue
			scamRes$accept$accept<-scamResHT$accept
			
			#----------------
			### Random effect
			#----------------
			scamResRandom<-scamRandom(data,start=scamRes$start,priors=priors,rotationVal=scamRes$rotation$valuesRandom,rotationVect=scamRes$rotation$vectorsRandom,wParamRandom=scamRes$rotation$wParamRandom,wOuterParamRandom=scamRes$rotation$wOuterParamRandom,wRandomCov=scamRes$rotation$wRandomCov,sdvalue=scamRes$sdvalue,accept=scamRes$accept,likelihoodtype=likelihoodtype,burning=(i < nburn),rotation=(i > nrot & i < nburn),iter=i)
			
			### Assign the new estimated pieces of scamResH to scamRes
			scamRes$start$paramRandom<-scamResRandom$start$paramRandom
			scamRes$start$RandomVar<-scamResRandom$start$RandomVar
			
			scamRes$rotation$valuesRandom<-scamResRandom$rotation$valuesRandom
			scamRes$rotation$vectorsRandom<-scamResRandom$rotation$vectorsRandom
			scamRes$rotation$wParamRandom<-scamResRandom$rotation$wParamRandom
			scamRes$rotation$wOuterParamRandom<-scamResRandom$rotation$wOuterParamRandom
			scamRes$rotation$wRandomCov<-scamResRandom$rotation$wRandomCov
			
			scamRes$sdvalue<-scamResRandom$sdvalue
			scamRes$accept<-scamResRandom$accept
			
			#-------------------------
			#### Store burning results
			#-------------------------
			if(i < nburn){
				burningObj$paramX[,,i]<-scamRes$start$paramX
				burningObj$paramTr[,,i]<-scamRes$start$paramTr
				burningObj$means[,,i]<-scamRes$start$means
				burningObj$sigma[,,i]<-scamRes$start$sigma
				burningObj$paramRandom[,,i]<-scamRes$start$paramRandom
				burningObj$RandomVar[i]<-scamRes$start$RandomVar
				
				if(details){
					detailsObj$sdvalue[,,i+1]<-scamRes$sdvalue$sdvalue
					detailsObj$acceptRatio[,,i+1]<-scamRes$accept$accept[,,1]/scamRes$accept$accept[,,2]
					
					detailsObj$sdvalueRandomCov[i+1]<-scamRes$sdvalue$sdvalueRandomCov[1]
					detailsObj$acceptRatioRandomCov[i+1]<-scamRes$accept$acceptRandomCov[1]/scamRes$accept$acceptRandomCov[2]
					
					detailsObj$sdvalueRandom[,,i+1]<-scamRes$sdvalue$sdvalueRandom
					detailsObj$acceptRatioRandom[,,i+1]<-scamRes$accept$acceptRandom[,,1]/scamRes$accept$acceptRandom[,,2]
				}
			}else{
				it<-i-nburn
				if(it==thin*itsave){
					modelres$paramX[,,itsave]<-scamRes$start$paramX
					modelres$paramTr[,,itsave]<-scamRes$start$paramTr
					modelres$means[,,itsave]<-scamRes$start$means
					modelres$sigma[,,itsave]<-scamRes$start$sigma
					modelres$paramRandom[,,itsave]<-scamRes$start$paramRandom
					modelres$RandomVar[itsave]<-scamRes$start$RandomVar
					itsave<-itsave+1
				}
			}
			
			### Print number of iterations completed
			if(verbose){
				if(any(i==printCounter)){
					print(i)
				}
			}
		}
	}
	
	if(!any(class(data)=="HMSCdataRandom")){
		modelres$sdvalue<-scamRes$sdvalue
		modelres$acceptRatio<-scamRes$accept[,,1]/scamRes$accept[,,2]
	}else{
		modelres$sdvalue<-scamRes$sdvalue$sdvalue
		modelres$acceptRatio<-scamRes$accept$accept[,,1]/scamRes$accept$accept[,,2]
		
		modelres$sdvalueRandomCov<-scamRes$sdvalue$sdvalueRandomCov
		modelres$acceptRatioRandomCov<-scamRes$accept$acceptRandomCov[1]/scamRes$accept$acceptRandomCov[2]
		
		modelres$sdvalueRandom<-scamRes$sdvalue$sdvalueRandom
		modelres$acceptRatioRandom<-scamRes$accept$acceptRandom[,,1]/scamRes$accept$acceptRandom[,,2]
	}
	
	#===================
	#### Name the pieces
	#===================
	### paramX
	dimnames(modelres$paramX)[[3]]<-paste("iter",1:dim(modelres[[1]])[3],sep="")
	dimnames(modelres$paramX)[[1]]<-rownames(start$paramX)
	dimnames(modelres$paramX)[[2]]<-colnames(start$paramX)
	
	### paramTr
	dimnames(modelres$paramTr)[[3]]<-paste("iter",1:dim(modelres[[1]])[3],sep="")
	dimnames(modelres$paramTr)[[1]]<-rownames(start$paramTr)
	dimnames(modelres$paramTr)[[2]]<-colnames(start$paramTr)

	### means
	dimnames(modelres$means)[[3]]<-paste("iter",1:dim(modelres[[1]])[3],sep="")
	dimnames(modelres$means)[[1]]<-colnames(start$paramX)
	dimnames(modelres$means)[[2]]<-rownames(start$paramX)
	
	### sigma
	dimnames(modelres$sigma)[[3]]<-paste("iter",1:dim(modelres[[1]])[3],sep="")
	dimnames(modelres$sigma)[[1]]<-colnames(start$paramX)
	dimnames(modelres$sigma)[[2]]<-colnames(start$paramX)
	
	### paramRandom
	if(any(class(data)=="HMSCdataRandom")){
		dimnames(modelres$paramRandom)[[1]]<-rownames(start$paramRandom)
		dimnames(modelres$paramRandom)[[2]]<-colnames(start$paramRandom)
	}
	
	#==========
	### Results
	#==========
	if(details){
		res<-list(call=call,model=modelres,burning=burningObj,details=detailsObj,data=data)
	}else{
		res<-list(call=call,model=modelres,burning=burningObj,data=data)
	}
	
	attr(res,"likelihood")<-likelihoodtype
	
	if(!any(class(data)=="HMSCdataRandom")){
		class(res)<-c("hmsc","HMSCdata","HMSCdataTrait")
	}else{
		class(res)<-c("hmsc","HMSCdata","HMSCdataTrait","HMSCdataRandom")
	}
	
	return(res)
}