make_jive: Create a list that can be used as an input to mcmc_bite

Description

This function creates a jive object from a matrix of intraspecific observations and species phylogeny. The obtained jive object is a list that can than be used as an input to mcmc_bite function Intraspecific observations should be stored as matrix, where lines are vector of observations for each species, with NA for no data. Phylogenetic tree can be either a simmap object (make.simmap) or phylo object (as.phylo)

Usage

make_jive(phy = NULL, traits, map = NULL, model.priors = list(mean =
  "BM", logvar = "OU"), scale = FALSE, nreg = NULL, lik.f = NULL,
  init = NULL)

Arguments

phy

phylogenetic tree provided as either a simmap or a phylo object

traits

matrix of traits value for every species of phy (see details)

map

matrix mapping regimes on every edge of phy (see details)

model.priors

list giving model specification for trait evolution preferably given along with variable names. Supported models are "OU", "BM", "WN". The user can also specify if the assumptions of the model should be relaxed and can also enter a function (see details)

scale

boolean indicating whether the tree should be scaled to unit length for the model fitting

nreg

integer giving the number of regimes for a Beast analysis. Only evaluated if phy == NULL

lik.f

alternative likelihood function of the form function(pars.lik, traits, counts) to model intraspecific variation (see details)

init

matrix giving initial values for parameters with the variables in rows and the species in columns (see examples)

Value

A list of functions and tuning parameters (of class "JIVE" and "list") representing the plan of the hierarchical model to parse into mcmc_bite.

Details

This function creates a jive object needed for mcmc_bite function. Trait values must be stored as a matrix, where lines are vectors of observations for each species, with NA for no data. Rownames are species names that should match exactly tip labels of the phylogenetic tree.

Phylogenetic tree must be provided as either simmap object or as a phylo object. If the phylogenetic tree is a phylo object but model specification indicates multiple regimes, user must provide a mapping of the regime in map. If you keep the phy = NULL options the JIVE object can only be parsed to the xml_bite function.

map is a matrix giving the mapping of regimes on phy edges. Each row correspond to an edge in phy and each column correspond to a regime. If map is provided the map from the simmap object is ignored.

trait evolution can be modeled with Ornstein-Uhlenbeck (OU), Brownian Motion (BM) or White Noise (WN) processes. Multiple regimes can be defined for both models and will apply on thetas: c("OU", "theta"), sigmas: c("OU", "sigma") or alphas: c("OU", "alpha") for OU and on sigmas only for WN: c("WN", "sigma") and BM: c("BM", "sigma"). While using the OU model, the user can also relax the stationarity of the root: c("OU", "root") and relax several assumptions at the same time c("OU", "root", "theta") Species-specific distributions are modeled as multivariate normal distributions. User defined functions of trait evolution can be used in model.priors. The function should be of the form: function(tree, x, pars) and return a loglikelihood value with "tree" being the phylogenetic tree, x being a vector of trait value of size equal to the number of species and ordered as tree$tip.label and pars should be a vector of model parameters (see examples)

parameters used in the different pre-defined models:

White Noise model (WN):

root: root value
sigma_sq: evolutionary rate, n regimes if "sigma" is specified in model.priors

Brownian Motion model (BM):

root: root value
sigma_sq: evolutionary rate, n regimes if "sigma" is specified in model.priors

Ornstein Uhlenbeck model (OU):

root: root value. Only used if "root" is specified in model.priors
sigma_sq: evolutionary rate, n regimes if "sigma" is specified in model.priors
theta: optimal value, n regimes if "theta" is specified in model.priors
alpha: strength of selection, n regimes if "alpha" is specified in model.priors

Examples

Run this code

# NOT RUN {
## Load test data
data(Anolis_traits)
data(Anolis_tree)
data(Anolis_map)

## JIVE object to run jive with single regimes
my.jive <- make_jive(phy = Anolis_tree, traits = Anolis_traits[,-3],
 model.priors = list(mean = "BM", logvar= c("OU", "root")))

## JIVE object to run jive with multiple regimes
my.jive <- make_jive(Anolis_tree, Anolis_traits[,-3], map = Anolis_map,
 model.priors =list(mean = "BM", logvar = c("OU", "theta", "alpha")))

## JIVE object to run jive from an ancestral state reconstruction (stochastic mapping)
# First generate simmap object
library(phytools)
n= length(Anolis_tree$tip.label)
trait = rep(0,n)
trait[c(4,3,14,16, 6,5)] = 1
names(trait) =  Anolis_tree$tip.label

mapped_tree=make.simmap(Anolis_tree, trait, model='SYM')
plotSimmap(mapped_tree)

my.jive <- make_jive(mapped_tree, Anolis_traits[,-3]
, model.priors = list(mean = "OU" , logvar = c("OU", "theta")))
 
 ## Jive object using another model of trait evolution (EB from mvMORPH)
 library(mvMORPH)
 early_burst <- function(tree, x, pars){
  suppressMessages(mvEB(tree, x, method = "inverse", optimization = "fixed", 
   echo = FALSE)$llik(pars, root.mle = FALSE))
 }
 
 my.jive <- make_jive(phy = Anolis_tree, traits = Anolis_traits[,-3]
, model.priors = list(mean = early_burst , logvar = c("OU", "root")))
 initial.values <- c(0.1, 1, 50)
 window.size <- c(0.1, 0.2, 1)
 proposals <- list("slidingWin", "slidingWin", "slidingWin")
 hyperprior <- list(hpfun("Gamma", hp.pars = c(1.1, 5)), hpfun("Gamma", hp.pars = c(3, 5)),
                     hpfun("Uniform", hp.pars = c(30, 80)))
 names(initial.values) <- names(window.size) <- c("sigma_sq", "beta", "root")
 names(proposals) <- names(hyperprior) <- c("sigma_sq", "beta", "root")
 my.jive <- control_jive(jive = my.jive, level = "prior", intvar = "mean",
  pars = names(initial.values), window.size = window.size,
  initial.values = initial.values, proposals = proposals, hyperprior = hyperprior)
 
 ## Jive object using another model of intraspecific variation (uniform model)
 lik_unif <- function(pars.lik, traits, counts){
   if(!"mid" %in% names(pars.lik)) stop("'mid' parameter cannot be found in model.priors")
   if(!"logrange" %in% names(pars.lik)){
    stop("'logrange' parameter cannot be found in model.priors")
   }

   min.sp <- pars.lik$mid - 1/2*exp(pars.lik$logrange)
   max.sp <- pars.lik$mid + 1/2*exp(pars.lik$logrange)
   
   log.lik.U <- sapply(1:length(traits), function(i){
   sum(dunif(traits[[i]], min.sp[i], max.sp[i], log = TRUE))
   })
   
   if (is.na(sum(log.lik.U))) {
     return(-Inf)
   } else {
     return(log.lik.U)
   }
 }
 
 init_unif <- sapply(Anolis_tree$tip.label, function(sp){
  logrange <- log(diff(range(Anolis_traits[Anolis_traits[,1] == sp, 3])) + 2)
  mid <- mean(range(Anolis_traits[Anolis_traits[,1] == sp, 3]))
  c(mid = mid, logrange = logrange)
 })
 
 my.jive <- make_jive(phy = Anolis_tree, traits = Anolis_traits[,-2],  
 model.priors = list(mid = "BM" , logrange = c("OU", "root")),
 lik.f = lik_unif, init = init_unif)
 
# }