R/basicStatHMSC.R

In HMSC: Hierarchical modelling of species community

basicStatHMSC<-function(x,stats=c("mean","sd")){
#### F. Guillaume Blanchet - June 2013
##########################################################################################
	statsbase<-c("mean","sd")
	stats<-match.arg(stats,statsbase,several.ok =TRUE)
	
	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)
	
	### Means
	if(any(stats=="mean")){
		means<-matrix(NA,nrow=nsites,ncol=nsp)
		rownames(means)<-rownames(x$data$Y)
		colnames(means)<-colnames(x$data$Y)
	}
	
	### standard deviation
	if(any(stats=="sd")){
		sds<-matrix(NA,nrow=nsites,ncol=nsp)
		rownames(sds)<-rownames(x$data$Y)
		colnames(sds)<-colnames(x$data$Y)
	}
	
	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)
			}
			
			### Means
			if(any(stats=="mean")){
				means[,i]<-apply(probsp,1,mean)
			}
			### Standard deviations
			if(any(stats=="sd")){
				sds[,i]<-apply(probsp,1,sd)
			}
		}
	}
	
	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)
			}
			### Means
			if(any(stats=="mean")){
				means[,i]<-apply(probsp,1,mean)
			}
			### Standard deviations
			if(any(stats=="sd")){
				sds[,i]<-apply(probsp,1,sd)
			}
		}
	}
	
	### Result object 
	if(any(stats=="mean") & !any(stats=="sd")){
		res<-means
	}
	if(!any(stats=="mean") & any(stats=="sd")){
		res<-sds
	}
	if(any(stats=="mean") & any(stats=="sd")){
		res<-array(dim=c(nsites,nsp,2))
		res[,,1]<-means
		res[,,2]<-sds
		dimnames(res)[[1]]<-rownames(x$data$Y)
		dimnames(res)[[2]]<-colnames(x$data$Y)
		dimnames(res)[[3]]<-c("mean","sd")
	}

	res
}