Using make.simmap with previously inferred transition rates

[Belinda333]

Dear all,

I am working with a binary trait for which I already used a BISSE analyses, which indicated not only that the speciation rate of my taxa of interested is state-dependent but also yielded the transition rates between trait states. I was wondering if I can use that inferred transition rates from the BISSE analyses for ancestral trait reconstruction with the make.simmap function? I tried inputting it under the model but I received the following error: mod<-matrix(c(0,0.005,0.216,0),2) mtrees<-make.simmap(tree,data, model=mod, nsim=500) Error in nlminb(q.init, function(p) lik(makeQ(m, p, index.matrix), pi = pi), : 'd' must be a nonempty numeric vector Is that a bug or am I inputting my transition rates incorrectly? Also, since I read in the manual that the values in the diagonal of the transition rate matrix will be ignored I just inserted a zero value- is that correct?

I would highly appreciate some help/feedback! Thanks a lot!

Best wishes, Belinda

[liamrevell]

Hi. I posted an explanation of how to do this on my blog, here.

[Belinda333]

Great!!! Thanks a lot. I will directly check it out!

[Belinda333]

I have just one more question; the as.Qmatrix doesn't work for any matrix one creates by oneself, so how can I modify it? You see, my problem is that the transition rates inferred with BISSE are quite different from the ones suggested using the ARD model and stochastic character mapping. Hence also the trait mapping results differ quite a bit between asr.marginal (using a model that accounts for state-dependent speciation) and make.simmap. Therefore I wanted to check how the results from the stochastic character mapping look like when I use the same transition rates as inferred by the BISSE anlyses ... Thanks a lot again!

[Belinda333]

Ok, I just used a trick and simply overwrote the values for the transition rates given under the fitted ARD model and now it seems to run.

[liamrevell]

Hi Belinda. Just create a matrix in the normal way. So, for instance, if your Q[1,2] = 0.5 and your Q[2,1] = 1.2 for states a and b then you could do:

Q<-matrix(c(0,0.5,1.2,0),2,2,byrow=TRUE)
diag(Q)<--rowSums(Q)
rownames(Q)<-colnames(Q)<-c("a","b")
Q

Does this make sense?

(Addendum: make.simmap will take a k by k matrix, for k states, as input -- not just an object of class "Qmatrix".)

[Belinda333]

Dear Liam,

ok. Thank you. Will try that too but before I always received an error when using a regular matrix... Anyways, my trick also worked and make.simmap was running using my specified transition rates. Thanks again for helping!

Best wishes, Belinda

[Belinda333]

Short update: no, when I it as regular matrix I get again an error: Error in EXPM(Q * tree$edge.length[j])[NN[j, 1], ] : subscript out of bounds

Nevertheless, the trick to modify the rates in the matrix inferred under ARD worked fine; I basically just did the following: fitARD[["rates"]]<-c(0.216,0.005) estQ<-as.Qmatrix(fitARD) class(estQ)<-c("Qmatrix", "matrix") isSymmetric(estQ) So no worries!