R/lik_wn.R

In bite: Bayesian Integrative Models of Trait Evolution

#' @import ape

# input: n - number of species, par.n - which parameter has been updated, pars - c(sig1, ..., sigN, the0), tree and map
# does: calculate log-likelihood; 
lik_wn <- function(x, n, pars, Pi, par.n, data, map){
  
  sigma <- pars[Pi[1,]==1]
  root <- pars[Pi[2,]==1]
  upd <- apply(Pi, 2,function(x) which(x == 1))[par.n]
  
  S <- map$S ##Mapping of the start and end of each epoch nepochs x 2
  gamma <- map$gamma ##Indicator mapping of the epochs lived by each species nepochs x n
  beta <- map$beta ##Indicator mapping of the regimes on each epoch nepochs x nreg
  
  ## Expectation
  if(any(upd == 2)){
    data$E  <- matrix(1, ncol(gamma), 1)
    data$E[,] <- root # ancestral mean
    rownames(data$E) <- colnames(gamma)
  }
  
  # Variance Matrix (SIGMA)
  if(any(upd == 1)){
    var.reg <- sigma * t((S[,2]-S[,1]) * beta)
    V <-  matrix(0, ncol(gamma), ncol(gamma))
    diag.elts <- 1 + 0:(ncol(gamma) - 1) * (ncol(gamma) + 1)
    V[diag.elts] <- colSums(var.reg) %*% gamma
    
    # determinant
    data$det <- as.numeric(determinant(V)$modulus)
    # inverse
    data$inv <- solve(V)
  }
  
  ## Log likelihood (multinormal)
  loglik <- (-n/2 * log(2 * pi) - 1/2 * data$det - 1/2 * (crossprod(x - data$E,data$inv)%*%(x - data$E)))
  
  if (is.na(loglik) | (class(loglik) == "try-error" )) {
    return(list(loglik = -Inf, data = data))
  } else {
    return(list(loglik = loglik, data = data))
  }
  
}