Hi all,

I am trying to examine a correlation between two discrete traits, using stochastic mapping.

I tried two options and ran into problems with both:

1. Using stochastic mapping of Phytools -

   make.simmap - The problem is that I need to use constraints - for one character I need to set the state of the root to 0 (I figured how to do that), while for the second character I need to set the reverse rate to 0. Couldn't figure out this second constraints. Tried this

   mtrees<-make.simmap(tree,mtree$states,nsim=10, model=matrix(c(0,1,0,0),2))

   But model=matrix(c(0,1,0,0)) leads to the following error:

   BLAS/LAPACK routine 'DGEBAL' gave error code -3
   In addition: Warning message:
   In matrix(XX$rates[II], m, m, dimnames = list(lvls, lvls)) :
   data length [3] is not a sub-multiple or multiple of the number of rows [2]
   The same command works if I use , model=matrix(c(0,1,0,0),2))
   

   In addition, I wonder if anyone has examined the distribution of two characters - i.e. how can I calculate the percentage of time spent in (0,0), (0,1), etc. for the classes generated by phytools.

2. Using Diversitree -

   I used asr.stoch with a likelihood function and parameters. I managed to make constraints in the two needed ways (i.e. on the root, and on the rates)

   pars <- c(.1, .1, .03, .03, .03, .06)
   set.seed(1)
   phy <- trees(pars, "bisse", max.taxa=4, max.t=Inf, x0=0)[[1]]
   h <- history.from.sim.discrete(phy, 0:1)
   plot(h, phy, main="True history")
   lik <- make.mk2(phy, phy$tip.state)
   #lik<-constrain(q13~q12,q21~q12,q23~q12,q31~q12,q32~q12)
   argnames(lik)
   lik2<-constrain(lik, q10~0)
   argnames(lik2)
   fit = find.mle(lik2,c(.1))
   st.s1 <- asr.stoch(lik2, fit$par)
   plot(st.s1, phy)
   st.s2 <- asr.stoch(lik2, fit$par)
   

   What I do not understand is how to simulate stochastic mappings on my trees randomly - I want to use the same rate parameters, but with no constraints on the tips. These simulations will help me determine whether the statistic I am examining is significantly large for my data.

   In addition, I wonder if anyone has examined the distribution of two characters - i.e. how can I calculate the percentage of time spent in (0,0), (0,1), etc. for the classes generated by Diversitree.

Thanks a lot for your help!