[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