muhisse_example.R

# muhisse example
# from https://speciationextinction.info/articles/MuHiSSE-vignette.html

# load packages
library(hisse)
library(diversitree)

# simulate simple four state model
pars <- c(.1,  .15,  .2, .1, # lambda 1, 2, 3, 4
          .03, .045, .06, 0.03, # mu 1, 2, 3, 4
          .05, .05, .00,        # q12, q13, q14
          .05, .00, .05,     # q21, q23, q24
          .05, .00, .05,     # q31, q32, q34
          .00, .05, .05)
set.seed(2)
phy <- tree.musse(pars, 30, x0=1)
states <- phy$tip.state

lik <- make.musse(phy, states, 4)
#lik <- make.musse(phy, states, 3)
diversitree.free = lik(pars)
print(diversitree.free)

#create states dataframe
states <- data.frame(phy$tip.state, phy$tip.state, 
                     row.names=names(phy$tip.state))
states <- states[phy$tip.label,]
states.trans <- states
for(i in 1:Ntip(phy)){
  if(states[i,1] == 1){
    states.trans[i,1] = 0
    states.trans[i,2] = 0
  }
  if(states[i,1] == 2){
    states.trans[i,1] = 0
    states.trans[i,2] = 1
  }
  if(states[i,1] == 3){
    states.trans[i,1] = 1
    states.trans[i,2] = 0
  }
  if(states[i,1] == 4){
    states.trans[i,1] = 1
    states.trans[i,2] = 1
  }
}

# set up and run a MuSSE model
turnover <- c(1,1,1,1)
extinction.fraction <- c(1,1,1,1)
f = c(1,1,1,1)

trans.rate <- TransMatMakerMuHiSSE(hidden.traits=0)
print(trans.rate)

states.trans <- cbind(phy$tip.label, states.trans)
dull.null <- MuHiSSE(phy=phy, data=states.trans, f=f, turnover=turnover, 
                     eps=extinction.fraction, hidden.states=FALSE, 
                     trans.rate=trans.rate)


# set up and run a MuHiSSE
turnover <- c(1,2,3,4,5,6,7,8)
extinction.fraction <- rep(1, 8) 
f = c(1,1,1,1)

trans.rate <- TransMatMakerMuHiSSE(hidden.traits=1)
print(trans.rate)

MuHiSSE <- MuHiSSE(phy=phy, data=states.trans, f=f, turnover=turnover, 
                   eps=extinction.fraction, hidden.states=TRUE, 
                   trans.rate=trans.rate)