R/marginal_lik.R

In bite: Bayesian Integrative Models of Trait Evolution

#' @title Calculate marginal likelihood by thermodynamic integration (LTI)
#' @description  Calculate the marginal likelihood from a logfile generated by \code{\link{mcmc_bite}} with thermodynamic integration (Lartillot and Philippe, 2006) or stepping stone (Xie et al., 2011).
#' @param mcmc.log the output file of a \code{\link{mcmc_bite}} run
#' @param burnin number or proportion of iteration to delete
#' @param method one of "TI" for thermodynamic integration and "SS" for stepping stone integration (the default)
#' @export
#' @return a length one numeric double giving the marginal likelihood of the model.
#' @author Theo Gaboriau and Simon Joly 
#' @examples 
#' ## Load test data
#' data(Anolis_traits)
#' data(Anolis_tree)
#' data(Anolis_map)
#' 
#' ## Run a MCMC chain with thermodynamic Integration
#' set.seed(300)
#' my.jive <- make_jive(Anolis_tree, Anolis_traits[,-3],
#'   model.priors = list(mean="BM", logvar="OU"))
#' bite_ex <- tempdir()
#' logfile <- sprintf("%s/my.jive_mcmc_TI.log", bite_ex)
#' mcmc_bite(my.jive, log.file=logfile, ncat=10, sampling.freq=10,
#'  print.freq=100, ngen=1000, burnin=0) 
#'  
#' ## import the results in R
#' res <- read.csv(logfile, header = TRUE, sep = "\t")
#'  
#' mlikTI <- marginal_lik(res, burnin = 0.1, method = "TI")
#' mlikTI
#' 
#' mlikSS <- marginal_lik(res, burnin = 0.1, method = "SS")
#' mlikSS
#' @encoding UTF-8


marginal_lik <- function(mcmc.log, burnin = 0, method = "SS") {

	if (!("temperature" %in% colnames(mcmc.log))) stop('No \'temperature\' column in the log file')
  
  temp <- unique(mcmc.log$temperature)
  
  if(burnin>0){
    if(burnin < 1) burnin <- sapply(temp, function(t) quantile(mcmc.log$iter[mcmc.log$temperature == t], burnin))
    else burnin <- sapply(temp, function(t) min(mcmc.log$iter[mcmc.log$temperature == t])) + burnin
    burn <- unlist(lapply(1:length(temp), function(t) mcmc.log[mcmc.log$temperature == temp[t],"iter"] <= burnin[t]))
  } else {
    burn <- rep(FALSE, nrow(mcmc.log))
  }
    
	
	# The marginal likelihood
	lik.ti = 0
	
	# Calculate the intervals in temperature
	intervals <- abs(diff(unique(mcmc.log$temperature)))
	
	
	if(method == "TI"){
	  
	  # Get mean marginal likelihood for each temperature
	  lik <- sapply(temp, function(t) mean(rowSums(mcmc.log[mcmc.log$temperature == t & !burn, grepl("prior", colnames(mcmc.log))], na.rm = TRUE))) 
	  
	  if (length(lik)<=1) warning('Only one temerature for the thermodynamic integration')
	  
	  
	  # Estimate the integral
	  for (i in 1:length(intervals)) {
	    lik.ti = lik.ti + ((lik[i]+lik[i+1])/2 * intervals[i])
	  }
	  
	 } else if(method == "SS"){
	   
	   ## Estimate RSS for each k and sum it
	   for(i in 1:length(intervals)) {
	     betaliks <- rowSums(mcmc.log[mcmc.log$temperature == temp[i] & !burn, grepl("prior", colnames(mcmc.log))])
	     Lmax <- max(betaliks)
	     n <- length(betaliks)
	     lik.ti = lik.ti + intervals[i]*Lmax + log(1/n * sum(exp(intervals[i]*(betaliks - Lmax)))) 
	   }
	   
	 } else {
	    stop(sprintf("%s is not a supported method"), method)
	 }
	  
	return(lik.ti)

}