comparison of phylogeographic models

[Arley Camargo]

Hi Beasties,
I am running phylogeographic analysis in Beast using the standard
random walk and the relaxed random walk sensu Lemey et al. (2009,
2010), and I would like now to compare the fit of these alternative
diffusion models.

1) Can they be compared using the MLE based on path and stepping-stone
sampling of Baele et al. 2012? Based on the tutorial, and some
preliminary runs, looks like the MLE sampling is performed after the
mcmc is complete, is this correct? I was wondering if it is possible
to perform the MLE in a separate Beast run instead of combining the
mcmc and MLE sampling in the same XML file. I have some mcmc runs that
took very long (~a week), and it would be very nice to perform MLE
without running the mcmc again.

2) I also noticed the the path and stepping-stone analyses provide a
marginal likelihood value but no indication of a variance or standard
deviation. Is there any way of estimating variance of the marginal
likelihood? should I run the sampling step several times for obtaining
such an estimate?

3) I would also like to obtain a plot similar to Fig. 3B of Lemey et
al. 2010. Skyride analysis shows a steady increase in population size
and I'd like to see if the diffusion rate also increases over time.
How can this joint plot of pop. size and diffusion rate be made?

thanks for your advice

arley

[Guy Baele]

Hi Arley,

I will answer some of your questions concerning MLE and rely on other BEAST users to answer the remaining questions.

1) Typically you would perform the MLE after your mcmc run. That way, the mcmc run functions as a (unnecessarily large) burn-in to your marginal likelihood estimation. It is possible to do the marginal likelihood estimation in parallel to your mcmc run, but I would strongly recommend doing an mcmc run first that is run long enough to pass the burn-in phase of your mcmc chain (but without running the full-length chain of course).

2) If you run the MLE using PS/SS twice and you get the same order of models, you should be fine. On the other hand, if you do one run but with enough path steps / ratios and corresponding chain lengths (but I don't know what the minimal values for the settings are), that should be fine as well. SS should converge towards the marginal likelihood value faster than does PS, so I guess you could say that if the difference between these two values is small enough, then your settings are probably stringent enough.

Cheers,

Guy

Op woensdag 16 mei 2012 14:42:50 UTC+1 schreef Arley Camargo het volgende:

[Arley Camargo]

Thanks Guy!
your reply is very helpful!
arley

2012/5/16 Guy Baele <bael...@gmail.com>

--
You received this message because you are subscribed to the Google Groups "beast-users" group.

To view this discussion on the web visit https://groups.google.com/d/msg/beast-users/-/hYqJ1bQAVh8J.

To post to this group, send email to beast...@googlegroups.com.
To unsubscribe from this group, send email to beast-users...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/beast-users?hl=en.

[Guy Baele]

Hi Arley,

Coming back to your first question, yes you can compare the standard random walk and the relaxed random walk sensu Lemey et al. (2009, 2010).

I wasn't sure about the details, so I asked Philippe Lemey (I'm working in his lab).

Sorry for taking so long to respond.

Cheers,

Guy

Op donderdag 17 mei 2012 13:44:01 UTC+2 schreef Arley Camargo het volgende:

2012/5/16 Guy Baele <bael...@gmail.com>

To unsubscribe from this group, send email to beast-users+unsubscribe@googlegroups.com.

[Arley Camargo]

Thanks Guy for the follow-up
In fact, I obtained MLE for two separate runs for each model as follows:

	SRW1	SRW2	RRW1	RRW2
AIC	15228.7145	15094.8925	22080.161	23678.6418
PS	-4230.362678	-4273.998253	-3865.656574	-3886.36949
SS	-4233.538871	-4277.119099	-3867.424673	-3886.198907

the results seem well consistent between runs and based on them, I believe a bayes factor can be calculated as the quotient between the MLE of the models, correct?:

ln BF = MLE SRW / MLE SRW

If this is correct, I am obtaining ln BF ~ 1, which is BF ~ 3 for SRW in comparison with RRW. As far as I know, this value represents an acceptable support for the SRW model in this case, correct?

Now, this suggest that the diffusion rate has remained constant through the tree, and that's why I wanted to plot dispersal rate as a function of time as in Fig. 3B of Lemey et al. 2010. I think I can summarize this parameter for each branch of the MCC tree but not across branches. I think Lemey used a custom script for obtaining these summary estimates but I have not been able to find the code in Lemey's lab software website.

I wanted to look at the change in the dispersal rate over time because BSP plot suggest recent population growth, and I was expecting a concomitant increase in dispersal rate like in Fig. 3B of Lemey et al. 2010. However, I think it is possible that a constant dispersal rate could be consistent with a steady population growth, do you think this intepretation is correct?

finally, I' am obtaining a dispersal rate estimate of ~1,000. I calibrated the substitution rate to million of years. This means that the dispersal rate is 1,000 km per million years? thus 1 meter per year?

thanks for your feedback
arley

2012/5/30 Guy Baele <bael...@gmail.com>

To view this discussion on the web visit https://groups.google.com/d/msg/beast-users/-/Fx0-btCIxY0J.

To post to this group, send email to beast...@googlegroups.com.

To unsubscribe from this group, send email to beast-users...@googlegroups.com.

[Guy Baele]

Hi Arley,

I agree that, given the large difference between the SRW and RRW models, the results are consistent between runs. There is still some variation on your outcomes, but that's ok. Could you post the settings you used for the MLE (i.e. chain length and number of path steps)?

To calculate the log Bayes factor, you can simply subtract (and not divide!!) the two log marginal likelihood estimates, for example: ln BF = -3865.656574 + 4230.362678 = 364.706104 which means you have very strong support in favor of the RRW model.

Why would you want to compare model fit between two independent MLE runs of the same model? In other words, calculating MLE SRW / MLE SRW does not make any sense to me ...

I suggest you have another close look at your results and what your conclusions should be. For the interpretation of (log) Bayes factors, I recommend reading Kass & Raftery (1995).

Cheers,

Guy

Op woensdag 30 mei 2012 18:06:02 UTC+2 schreef Arley Camargo het volgende:

[Arley Camargo]

Thanks Guy!
now the calculations makes more sense and I'm glad that the evidence for the RRW model is very decisive. that makes sense with the population growth. I'll check now with Lemey's code if I can see a change in dispersal rate consistent with the change in pop. size. and also with paleo-distribution models.

I was somewhat surprised that the AICM was always lower for the SRW model, but I guess this might be an example of the lower reliability of this method for model comparison based on the performance evaluations of your paper.

here are the settings I used for MLE:

    <marginalLikelihoodEstimator chainLength="500000" pathSteps="100" pathScheme="betaquantile" alpha="0.30">
        <samplers>
            <mcmc idref="mcmc"/>
        </samplers>
        <pathLikelihood id="pathLikelihood">
            <source>
                <posterior idref="posterior"/>
            </source>
            <destination>
                <prior idref="prior"/>
            </destination>
        </pathLikelihood>
        <log id="MLE.SRW" logEvery="1000" fileName="MLE.SRW.log">
            <pathLikelihood idref="pathLikelihood"/>
        </log>
    </marginalLikelihoodEstimator>

cheers,
a

2012/5/31 Guy Baele <bael...@gmail.com>

To view this discussion on the web visit https://groups.google.com/d/msg/beast-users/-/Rkhh_EQ-XIwJ.

To post to this group, send email to beast...@googlegroups.com.

To unsubscribe from this group, send email to beast-users...@googlegroups.com.