Setting priors for snapper in BEAUti

[bmichanderson]

Hi,

I'm fairly inexperienced setting priors in BEAUti, and I was confused trying to interpret what is stated in the SNAPPER paper (Stoltz et al. 2021) and in the SPEEDEMON paper (Douglas & Bouckaert 2022) with the options available in BEAUti for the SNAPPER template and YSC.

1. In Stoltz et al. (2021), they model theta with a "CIR" prior in SNAPP and SNAPPER. I can see this as an option in SNAPP, but not in the SNAPPER template in BEAUti. I don't understand exactly what the prior is doing (linking adjacent pop sizes?), nor what settings would make the most sense, but I can't seem to use it in SNAPPER.

2. Under the Yule Skyline Collapse prior, the birth rate appears to be drawn from a gamma distribution with a set shape and a rate that is sampled during the run. I want to set this sensibly for realistic values of lambda (calculated from my own data based on substitutions per site and a range of expected species numbers). In Douglas & Bouckaert 2022, it seems they set the birth rate to come from a gamma with alpha = 2 and beta from a Log-normal of mean = 1.59 and sd = 0.2. Since the two values are described as "shape" and "rate" in the BEAUti form, I assume that means the mean of the gamma distribution is 2 / 5 = 0.4 (assuming the mean of the log-normal is in log space, so real space = 5?). But comparing that to typical lambda values in the delimitation tutorial, it seems too low. My own data suggest a lambda from about 6 to 30. On the other hand, if the "rate" is instead "scale", then mean would be 2 * 5 = 10, which seems more sensible to me. Or maybe I'm just not understanding what is going on. I guess I would like to know how to set the prior for these to match a lambda I infer from my data.

I realize that SPEEDEMON has just come out, so it may be a bit before there is a tutorial or some hints about setting priors properly (e.g. how many groups, whether to check equal epochs and linked means and what that means biologically). Thanks for any help clearing this up!

Cheers,

Ben

[rbouckaert]

Hi Ben,

1. The CIR prior is only available in SNAPPER by editing the XML: just copy the XML fragment generated by BEAUti through the SNAPP template into the XML with the SNAPPER template. It is indeed a prior where child branches have a mean that is a function of the population of the parent branch.

2. With the Yule Skyline Collapse model, the parameters are indeed shape and rate, so the mean is shape/rate. I do not know how sensitive the result is to this rate prior, but found that the BICEPS prior is fairly insensitive, so perhaps it is the same here. Regardless, you can run for a few different values of the prior to see whether it makes any difference.

For SPEEDEMON, you want to check equal epochs, since the collapsed nodes will otherwise interfere with the groups. Linking epochs means that the mean of the previous epoch is used as prior mean for the next epoch. The alternative implies each epoch independent of other epochs draw from the same prior distribution. The number of groups is determined by how much information there is in the timing of internal nodes in the tree. Not having enough nodes means there is no information to inform the birth rate, too many nodes can be too restrictive since a constant birth rate is assumed throughout the epoch. When linking, it seems at least 6 internal nodes per epoch is enough, though we are still busy gaining experience in what makes most sense there.

Cheers,

Remco

[rbouckaert]

Just confirming that the Yule Skyline appears not very sensitive to the prior on birth rate rate: running a log normal with mean in real space for the values of mean in 0.1, 1.0, 10 and 100 results in nearly identical birth rate estimates for the HCV example that comes with BICEPS. It uses 10 epochs for 63 taxa, so the prior on the birth rate rate is already overtaken by the data in the first epoch. Of course, your data may be more sensitive to the prior than is the case of the HCV data.

[bmichanderson]

Thanks for looking into that Remco!
Just to clarify setting the CIR prior, I have these lines in my SNAPPER xml:

            <prior id="SnapperCoalescentRatePrior.t:speedemon_rand" name="distribution" x="@snapperCoalescentRate.t:speedemon_rand">
                <Gamma id="Gamma.2" name="distr">
                    <parameter id="RealParameter.4" spec="parameter.RealParameter" estimate="false" name="alpha">2.0</parameter>
                    <parameter id="parameter.hyperGamma-beta-SnapperCoalescentRatePrior.t:speedemon_rand" spec="parameter.RealParameter" estimate="false" name="beta">0.02</parameter>
                </Gamma>

and these lines in my SNAPP xml:

<distribution id="snapprior.speedemon_rand" spec="snap.likelihood.SnAPPrior" coalescenceRate="@coalescenceRate" lambda="@lambda" rateprior="CIR" tree="@Tree.speedemon_rand">
                <parameter id="alpha" spec="parameter.RealParameter" estimate="false" lower="0.0" name="alpha">2.0</parameter>
                <parameter id="beta" spec="parameter.RealParameter" estimate="false" lower="0.0" name="beta">0.02</parameter>
                <parameter id="kappa" spec="parameter.RealParameter" estimate="false" lower="0.0" name="kappa">1.0</parameter>
            </distribution>

Which parts of the distribution definition need to be transferred? The entries don't look equivalent to me (e.g. "spec" and "coalescenceRate" are not in the SNAPPER entry). Is it just the rateprior="CIR" bit plus the alpha, beta and kappa? Do I need to change SnAPPrior to something that works with SNAPPER?

Cheers,

Ben