[R-sig-phylo] phangorn package; pml tree

Klaus Schliep klaus.schliep at gmail.com
Tue Jul 20 14:23:43 CEST 2010


Dear Erwan,

there are two different functions 'pml' which estimates the likelihood for
the given parameters and returns an object of class 'pml' and 'optim.pml'
which performs optimisation on an object of class 'pml'.

In your case you first to construct an object of class 'pml'

assume you have the base frequencies given in a vector bf like
bf = c(0.1, 0.2, 0.3, 0.4)

fit = pml(mytree, mydata, k=4, bf=bf)
fit        #tells you a little bit about your model so far.

# if you want to use a Gamma model you need to specify the
# number of rate classes via k (here k=4)
# update.pml(fit) is a convenience function if you want to change just
# one or two parameters e.g.
fit = update(fit, bf = c(.25, .25, .25, .25))
# changes the base frequencies of your 'pml' object

#The next step is to optimise the model:

fit1 = optim.pml(fit, optNni=TRUE, optGamma=TRUE, model='GTR')
fit1


# this optimises a GTR model (base frequencies and rate matrix),
# performs NNI rearrangements and optimises the shape parameter.
# If you want to have the base frequencies fixed you would use

fit2 = optim.pml(fit, optNni=TRUE, optQ=TRUE, optGamma=TRUE)
# or by calling the appropriate  model: 'SYM'
# fit2 = optim.pml(fit, optNni=TRUE, optGamma=TRUE, model='SYM')
fit2

anova(fit2, fit1)

I hope this helped you a bit. If you have further question, do not bother
to send me another mail.

On 7/19/10, Erwan DELRIEU-TROTTIN <erwan.delrieu-trottin at etu.upmc.fr> wrote:
>
> Hi everyone,
>
> Phangorn package allows to perform ML trees after having chosen a model.
> I've run jmodeltest which has told me that a model GTR+G would be
> appropriate for my data.
> I would like to enter the different parameters in the pml function in
> order to draw a tree.
> There's an example for pml but I don't fully understand it.
>
> To perform a JC + Gamma + I - model, they do:
>    fitJC_GI <- update(fitJC, k=4, inv=.2)
> I is the the proportion of invariable site, but I don't understand why
> .2 is chosen in that exemple.
> What would it be if you wanted to perform GTR+G? or +F?
What is an GTR + F model??

>
> We can fix bf, the base frequencies, and I would like to, but I do not
> find how to code it: should I put the different values between brackets?
>

> I hope that my approach is appropriate.
>
> Many thanks in advance for your help.
>
> Erwan Delrieu-Trottin
>
> --
> Erwan Delrieu-Trottin, PhD student
> USR 3278 CNRS - EPHE
> Centre de Biologie et Ecologie Tropicale et Méditerranéenne
> Université de Perpignan, 52 Av. Paul Alduy
> 66860 Perpignan cedex, France
> e-mail : erwan.delrieu-trottin at etu.upmc.fr
>
> _______________________________________________
> R-sig-phylo mailing list
> R-sig-phylo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
>

Cheers,
Klaus



Dr. Klaus Schliep
Postdoctoral Fellow
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris



More information about the R-sig-phylo mailing list