[R-sig-phylo] transition matrix in geiger (sim.char)

Emmanuel Paradis Emmanuel.Paradis at mpl.ird.fr
Wed Feb 3 18:39:15 CET 2010


As a side note, ape 2.5 (available on CRAN since yesterday) has a 
function rTraitDisc that does the same job. The diagonal of the rate 
matrix, if given, is always ignored. It's possible to specify 
Alexandre's model as:

rTraitDisc(tree, k = 4, rate = 0.5)

It's also possible to specify the equilibrium frequencies (assumed 
balanced by default). And if you're not happy with the traditional 
Markovian models, you can make your own model with an R function.

Emmanuel

Liam J. Revell wrote on 03/02/2010 17:39:
> Hi Alex,
> The columns in your matrix should sum to 0.0, so if you want the 
> transitions to any of the four states to be equiprobable, try:
> 
>  > q<-list(rbind(c(-1, 1/3, 1/3, 1/3),
>        c(1/3, -1, 1/3, 1/3),
>        c(1/3, 1/3, -1, 1/3),
>        c(1/3, 1/3, 1/3, -1)));
> 
> Then type:
> 
>  > results<-sim.char(tree,q,nsims=1,model="discrete",root.state=1);
> 
> Note that you don't have to set the root.state, but if not assigned it 
> will assume a root.state of 1.
> 
> If you want to simulate many characters at once under the same model, 
> you can do that in one of two ways, depending on how you want the output 
> to be stored.  You can increase nsims, e.g.:
> 
>  > results<-sim.char(tree,q,nsims=10,model="discrete",root.state=1);
> 
> for 10 simulations.  Or, you can add additional matrices to the list of 
> transition matrices, q, e.g.:
> 
>  > for (i in 2:10){
>    q[[i]]=q[[1]];
> }
> 
> The latter gives your results in a more convenient form, from which:
> 
>  > results<-as.matrix(results[,,1]);
> 
> puts all the results from all 10 simulations into one numspecies x 10 
> matrix.
> 
> Good luck!
> 
> - Liam
> 
> Liam J. Revell
> NESCent, Duke University
> web: http://anolis.oeb.harvard.edu/~liam/
> NEW email: lrevell at nescent.org
> 
> 
> 
> Alexandre Antonelli wrote:
>> Hi,
>>
>>    I am trying to use the function sim.char in geiger without 
>> success.  I would like to let a character (in this case 'species 
>> distribution',  with four states) evolve along a tree, and try 
>> different transition  costs in the matrix.
>>
>>   However, I don't manage to get more than one character with 2  
>> states. Any ideas what I am doing wrong? It is very unclear to me how  
>> the transition matrix should be compiled (if columns/rows have to sum  
>> up to one, etc). The manual indicated that the number of states per  
>> character would be dependent on the size of the matrix, but that does  
>> not seem to be the case.
>>
>>   These are the commands I used, and the transition matrix I intended  
>> to apply:
>>
>> ## 0: Area A, 1:Area B, 2:Area C, 3:Area D
>> ##
>> ##        From     0    1    2    3
>> ## To    0        -.5    .5    .5    .5
>> ##        1        .5    -.5    .5    .5
>> ##        2        .5    .5    -.5    .5   
>> ##        3        .5    .5    .5    -.5
>>
>>
>> require (geiger)
>>
>> tree <- read.nexus(file="input.tre")
>>
>> q<-list (rbind(    c(-.5, .5, .5, .5),
>>         c(.5, -.5, .5, .5),
>>         c(.5, .5, -.5, .5),
>>         c(.5, .5, .5, -.5)))
>>
>> results<-sim.char(tree, q, model="discrete",n=2)
>>
>> ------
>> Thanks for any help!
>>
>> Best wishes,
>>
>> Alex
>>
>>
>> Dr Alexandre Antonelli
>> Institute of Systematic Botany
>> Zollikerstrasse 107
>> CH 8008  Zürich
>> Switzerland
>> Phone: +41 (0)44 634 8416
>> Mobile: + 41 (0) 76 7345116
>> alexandre.antonelli at systbot.uzh.ch
>> http://www.systbot.uzh.ch/Personen/PostdoktorandInnen/AlexandreAntonelli.html 
>>
>> ******************************************************
>>
>>
>>     [[alternative HTML version deleted]]
>>
>>   
>> ------------------------------------------------------------------------
>>
>> _______________________________________________
>> R-sig-phylo mailing list
>> R-sig-phylo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
>>
> 
> _______________________________________________
> R-sig-phylo mailing list
> R-sig-phylo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
> 

-- 
Emmanuel Paradis
IRD, Montpellier, France
   ph: +33 (0)4 67 16 64 47
  fax: +33 (0)4 67 16 64 40
http://ape.mpl.ird.fr/



More information about the R-sig-phylo mailing list