R/quantile.hmsc.R

In HMSC: Hierarchical modelling of species community

quantile.hmsc<-function(x,probs=c(0.025,0.5,0.975),type=8,...){
####
#### Calculates quantiles of the probabilities for an x of class hmsc
####
#### Argument :
####
#### x : x of class hmsc
#### probs : numeric vector of probabilities with values ranging from 0 to 1.  (Values up to 2e-14 outside that range are accepted and moved to the nearby endpoint.)
#### type : an integer between 1 and 9 selecting one of the nine quantile algorithms detailed below to be used. Default is 8 as suggested by Hyndman and Fan (1996).
#### ... : Arguments passed from quantile.default
####
#### F. Guillaume Blanchet - June 2013
##########################################################################################
	X<-as.matrix(x$data$X)
	
	if(any(class(x)=="HMSCdataRandom")){
		Random<-as.matrix(x$data$Random)
	}
	
	nsp<-dim(x$model$paramX)[1]
	nparam<-dim(x$model$paramX)[2]
	niter<-dim(x$model$paramX)[3]
	nsites<-nrow(X)
	nprobs<-length(probs)
	
	quant<-array(dim=c(nsites,nsp,nprobs))
	dimnames(quant)[[1]]<-rownames(x$data$Y)
	dimnames(quant)[[2]]<-colnames(x$data$Y)
	dimnames(quant)[[3]]<-as.character(probs)
	
	if(class(x)[2]=="HMSCdata"){
		for(i in 1:nsp){
			if(any(class(x)=="HMSCdataRandom")){
				modelsp<-X%*%x$model$paramX[i,,]+Random%*%x$model$paramRandom[i,,]
			}else{
				modelsp<-X%*%x$model$paramX[i,,]
			}
			
			#### Apply the inverse logit transformation (this approach is quicker than using binomial(link="logit")$linkinv)
			if(attr(x,"likelihood")=="logit"){
				probsp<-ilogit(modelsp)
			}
			if(attr(x,"likelihood")=="poisson"){
				probsp<-ilog(modelsp)
			}
			quant[,i,]<-t(apply(probsp,1,quantile,probs=probs,type=type))
		}
	}
	if(class(x)[2]=="HMSCformula"){
		for(i in 1:nsp){
			modelsp<-matrix(NA,nrow=nsites,ncol=niter)
			for(j in 1:niter){
				modelsp[,j]<-parseformula(x$formula,x$model$paramX[i,,j],x$data$X)$model
			}
			
			#### Apply the inverse logit transformation (this approach is quicker than using binomial(link="logit")$linkinv)
			if(attr(x,"likelihood")=="logit"){
				probsp<-ilogit(modelsp)
			}
			if(attr(x,"likelihood")=="poisson"){
				probsp<-ilog(modelsp)
			}
			quant[,i,]<-t(apply(probsp,1,quantile,probs=probs,type=type))
		}
	}
	quant
}