Dear Matt,

BEB is really just a hack to attempt to approximate the results of a fully Bayesian inference without actually doing it; it allows one to incorporate sampling errors in parameter estimates.

Generally speaking, if your data set is small and parameter estimate errors are lare (the most likely scenario where NEB and BEB might disagree), there is very little power to detect individual sites under selection, regardless of which method you use. The fact that LRT is significant (for omega>1) does not guarantee that you will be able to find at least one site under selection, e.g. the LRT could be due to small contributions from a number of sites, but individually, no single site is reliably significant.

If you are really interested in pinning down the significance of your findings, do some parametric simulations; generate 100 data sets using the omega distribution (and all other parameters, such as base frequencies and branch lengths) inferred from your alignment, and see what type of false positive error in detecting individual sites under selection you get. Then you can tweak the bayes factor cutoff to make that small.

Cheers,
Sergei

Dear Matt,

Ah yes - because likelihood function states are attached to data panels, they are not available for reloaded data (even though parameter estimate tables load just fine). This is a fairly fundamental design assumption and will not be easy to change, I am afraid (it has to do with the need for the GUI to know how to build a copy of the likelihood function on simulated data, and likelihood functions created in the batch language directly are next to impossible to clone runtime).

Do this however:

1). Load your model fit (e.g. REL); simulate 100 copies from it using the little gear wheels icon from the console window (you should get a Simulate LF option).

2). Write a little wrapper script following Multimedia File Viewing and Clickable Links are available for Registered Members only!!  You need to Login to run the same (e.g. REL fit) on every data set.

You would effectively make a loop

Code:

for (file=0; file<100; file=file+1)
{
     set up analysis options
     ExecuteAFile ("analysis file", options);
     Record relevant results
     (or better yet, just save a copy of the entire
	likelihood function for each replicate, and
	then write a separate processor script to
	tabulate what you need, say omega or
	sites under selection)
}

 

Cheers,
Sergei

Dear Matt,

Instead of "d", place an empty string "" in the stdinArray structure to terminate the variable length selection.

HTH,
Sergei

Dear Andrew,

If you want general statistical books on bootstrap, then any by the inventor of the procedure (Bradley Efron) would do. In your case, the application is rather simple

0). Fit your data using some model M, with estimated parameter vector R (either direct model parameters, like branch lengths, or functions of model parameters and the data, such as the Bayes factor for dN/dS at site 10).
1). The objective of the bootstrap is to estimate sampling properties of elements of R (e.g. mean, variance etc) the the function of random data. Ideally, one would prefer to repeat the same analysis on a large number of independent data samples which were generated using the same process as that which produced the observed data. Since such independent replicates are not available, we do the next best thing and assume that the model we estimated was what actually generated the data, and then simulate under it. Efron showed that this actually works really well asymptotically.

In your case, for each simulated data, record the actual dN/dS for each site.
2). Run the same inference procedure as the one you used in step 0.
3). Tabulate the proportions of (p_F) false positives (sites with dN<dS which were inferred to be under selection) and of true positives (p_T) (sites with dN>dS which were inferred to be under selection) as a function of the significance level of your test (e.g. the Bayes Factor = 50, 100, 200 etc). If you plot p_T vs p_F for a fixed value of the Bayes factor, you will get what is known as an ROC curve.
4). Decide what Bayes factor gives you decent performance (low p_F and high p_T, decreasing the BF will increase both, increasing the BF will decrease both)
5). Now use that Bayes factor cutoff to reparse your original results and claim that the sites you found under selection are reasonably robust based on bootstrapped sampling properties of the estimator.

Anisimova and Yang had two papers on this in MBE in 2001/2002, and our Not So Different ... paper in MBE 2005 uses the same procedure.

HTH,
Sergei

P.S. Take a look at Multimedia File Viewing and Clickable Links are available for Registered Members only!!  You need to Login for a canned set of scripts we used for the Not So Different paper.