R/loglikelihood.default.R

In HMSC: Hierarchical modelling of species community

loglikelihood.default <-
function(y,X,paramX,Random=NULL,paramRandom=NULL,likelihoodtype="logit",...,prior){
#### Check to make code faster
####
#### F. Guillaume Blanchet - Avril 2013
##########################################################################################
	#### Transform all objects into matrix
	y<-as.matrix(y)
	X<-as.matrix(X)
	paramX<-matrix(paramX,ncol=ncol(X))
	
	if(!is.null(Random)){
		Random<-as.matrix(Random)
		paramRandom<-matrix(paramRandom,ncol=ncol(Random))
	}
	
	if(is.null(Random)){
		### Without random effect
		rawModel<-tcrossprod(X,paramX)
	}else{
		### With random effect
		rawModel<-tcrossprod(X,paramX)+tcrossprod(Random,paramRandom)
	}
	
	if(likelihoodtype=="logit"){
		model<-ilogit(rawModel)
		
		like<-(model*y)+((1-model)*(1-y))
#		like<-ifelse(y==1,model,1-model)
		if(any(like==0)){
			### Small correction for values exactly equal to 0
			like[which(like==0,arr.ind=TRUE)]<-.Machine$double.eps
		}
		resFinal<-sum(log(like))
	}
	if(likelihoodtype=="poisson"){
		model<-ilog(rawModel)
		resFinal<-sum(y*rawModel-model)
	}
	if(likelihoodtype=="normal"){
		rawModel
		resFinal<- -1/(2*prior) * sum((y-rawModel)^2)
	}
		return(resFinal)
}