R/scamH.formula.R

In HMSC: Hierarchical modelling of species community

scamH.formula <-
function(formula,data,start,priors,rotationVal,rotationVect,wParamX,wOuterParamX,sdvalue,accept,likelihoodtype="logit",burning=TRUE,rotation=TRUE,iter,...){
#### F. Guillaume Blanchet - Janvier 2013, March 2013, April 2013
##########################################################################################
	#========================
	#### Initiation of object
	#========================
	nparamX<-ncol(start$paramX)
	nsp<-ncol(data$Y)
	
	paramXCurrent<-start$paramX
	means<-start$means
	sigma<-start$sigma
	
	#### Calculate log-likelihood - mean
	likeCurrentMean<-numeric()
	for(i in 1:nsp){
		dataLL<-data.frame(data$Y[,i],data$X)
		colnames(dataLL)[1]<-attributes(data)$Ypattern
		likeCurrentMean[i]<-loglikelihood.formula(formula,paramXCurrent[i,],dataLL,likelihoodtype=likelihoodtype)
	}
	
	#### Calculate log-likelihood - sigma
	likeCurrentSigma<-numeric()
	for(i in 1:nsp){
		likeCurrentSigma[i]<-loglikelihood(sigma,paramXCurrent[i,],means,likelihoodtype="sigma")
	}
	
	#=============
	#### Algorithm
	#=============
	#----------------------------------------
	### Update means (each one independently)
	#----------------------------------------
	paramXNew<-matrix(NA,nrow=nsp,ncol=nparamX)
	colnames(paramXNew)<-colnames(paramXCurrent)
	
	for(k in 1:nsp){
		for(j in 1:nparamX){
			### Sample from a normal distribution and weight the sample by the eigenvalue of the covariance matrix
			paramXNew[k,]<-paramXCurrent[k,]+rnorm(1,mean=0,sd=sdvalue[k,j])*sqrt(rotationVal[k,j])*rotationVect[,j,k]
			
			#### Calculate log-likelihood - mean
			dataLL<-data.frame(data$Y[,k],data$X)
			colnames(dataLL)[1]<-attributes(data)$Ypattern
			likeNewMean<-loglikelihood.formula(formula,paramXNew[k,],dataLL,likelihoodtype=likelihoodtype)
			
			#### Calculate log-likelihood - sigma
			likeNewSigma <- loglikelihood(sigma,paramXNew[k,],means,likelihoodtype="sigma")
			
			accept[k,j,2]<-accept[k,j,2]+1
			
			### Reject if the likelihood is -Inf
			if(!(is.infinite(likeNewMean) & sign(likeNewMean)==-1)){
				### Accept or reject the new paramX
				if(runif(1) <= exp(likeNewMean+likeNewSigma-likeCurrentMean[k]-likeCurrentSigma[k])){
					paramXCurrent[k,]<-paramXNew[k,]
					likeCurrentMean[k]<-likeNewMean
					likeCurrentSigma[k]<-likeNewSigma
					accept[k,j,1]<-accept[k,j,1]+1
				}
			}
		}
	}
	
	#--------------------------
	#### means and sigma update
	#--------------------------
	meansparamX<-colMeans(paramXCurrent)
	SigmaparamX<-crossprod(sweep(paramXCurrent,2,meansparamX,FUN="-"))
	
	MeansSp<-(priors$kappa0/(priors$kappa0+nsp))*priors$means0+(nsp/(priors$kappa0+nsp))*meansparamX
	kappaSp<-priors$kappa0+nsp
	nuSp<-priors$nu0+nsp
	
	LambdaSp<-priors$Lambda0 + SigmaparamX + (priors$kappa0*nsp)/(priors$kappa0+nsp)*outer(meansparamX-priors$means0,meansparamX-priors$means0)
	iLambdaSp<-solve(LambdaSp)
	
	isigma <- rWishart(1,nuSp,iLambdaSp)[,,1]
	sigma <- solve(isigma)
	sigmaKappa<-sigma/kappaSp
	
	means<-mvrnorm(1,MeansSp,sigmaKappa)
	
	if(burning){
		wAccept<-1-0.1*exp(-iter/500) # w
		### Greediness during the adaptation period
		greedAdapt<-1+exp(-iter/500) # q
		
		### Calculate the accept ratio
		acceptRatio<-accept[,,1]/accept[,,2]
		
		### Modify the standard deviation
		sdvalue<-sdvalue*greedAdapt^(acceptRatio-0.44)
		
		### Weight the accept ratio
		accept<-accept*wAccept
		
		### Weight the current ParamX
		wParamX<-wParamX+wAccept*paramXCurrent # s1
		### Weight the outer product of the current ParamX
		for(j in 1:nsp){
			wOuterParamX[,,j]<-wOuterParamX[,,j]+wAccept*outer(paramXCurrent[j,],paramXCurrent[j,]) # s2
		}
		
		if(rotation){
			### sum the weight
			weight<-sum(1-0.1*exp(-(1:iter)/500))
			
			for(j in 1:nsp){
				### Construct a weighted covariance matrix
				covMat<-(wOuterParamX[,,j]-outer(wParamX[j,],wParamX[j,])/weight)/(weight-1)+diag(sqrt(.Machine$double.eps),nparamX)
				rotationBase<-eigen(covMat)
				
				rotationVal[j,]<-rotationBase$values
				rotationVect[,,j]<-rotationBase$vectors
			}
		}
	}
	
	#### Result object
	colnames(sigma)<-names(means)
	rownames(sigma)<-names(means)
	startNew<-list(paramX=paramXCurrent,sigma=sigma,means=means)
	rotNew<-list(values=rotationVal,vectors=rotationVect,wParamX=wParamX,wOuterParamX=wOuterParamX)
	res<-list(start=startNew,rotation=rotNew,sdvalue=sdvalue,accept=accept)
	
	return(res)
}