Variable parameters in rootFrequencies and parameter combinations in substModel

[jjmccollum]

This is more a question about existing features than a new feature request; if these features are not already native to BEAST, then I imagine they'll have to be implemented in a plugin.

I'm using BEAST for a text-critical phylogenetic application. In this setting, it is common for the user to make broad, but not precisely quantified judgments of relative likelihood between different variant readings (= states) at the root of the tree. These judgments would normally correspond to the numerical entries of the @value attribute of a frequency element at each site, as in the following example:

<rootFrequencies id="rootfreqs.character1" spec="Frequencies">
    <frequencies id="rootfrequencies.character1" spec="parameter.RealParameter" value="0.2 0.8" estimate="false"/>
</rootFrequencies>

But to accommodate the more qualitative nature of the user's judgments, I'd like to employ common parameters that represent odds ratios between root frequency values and can be estimated as part of the model. For instance, if we had one parameter A representing "much more likely" and one parameter B representing "slightly more likely," we would have the constraint 1 < B < A on these parameters, and instead of value="0.2 0.8" in the above example, we would somehow encode that the frequency of state 1 is A times the frequency of state 0. Can this be done using existing mechanisms in BEAST XML?

Similarly, for substitution models, I would like to use custom rate parameters that correspond to different mechanisms of scribal error and estimate the values of these rate parameters. As is common with text-critical phylogenetic data, the substitution models vary from site to site. Certain mechanisms may appear for some state transitions and not others, and more than one mechanism may explain the same state transition. (For the latter scenario, different mechanisms are assumed to be mutually exclusive.) Now, for two different rate parameters α and β, I know that I can already encode a simple transition rate matrix like

-α  α
 β -β

with something like

<substModel id="SubstModel.character1" spec="GeneralSubstitutionModel">
    <frequencies id="freqs.character1" spec="Frequencies">
        <frequencies id="frequencies.character1" spec="parameter.RealParameter" value="0.5 0.5" estimate="false"/>
    </frequencies>
    <parameter id="rates.character1" spec="parameter.CompoundValuable" name="rates">
        <var idref="alpha_rate"/>
        <var idref="beta_rate"/>
    </parameter>
</substModel>

for each site. But if I want something like

-α    α
α+β -α-β

where the rate parameters combine in certain ways, is there a way I can encode that in BEAST XML using existing mechanisms?

[rbouckaert]

There is no convenient way to do this as far as I am aware off other than write a few classes and put them in a package. I am happy to assist if that seems too complicated.

One thing to be aware of in the example you gave is that there is only a single degree of freedom in the binary "general" substitution model: rates are normalised to give an expected rate of 1, so α and β are not identifiable and perhaps parameterising as

-α      α
1-α   -1+α

and restricting α to be in the range 0 to 0.5 could be an alternative to ensure the α<β.

[jjmccollum]

@rbouckaert What I have in mind is that the rate parameters α, β, etc. would be common over all sites. (Ideally, we would have one copy of each such parameter per branch, as we do for branch length, but this may be too complex, so making the rate variables global over all sites and branches would also be acceptable.) In this case, if I have at least one site whose substitution matrix contains a 1, would this allow us to determine the rate parameters even at sites whose matrices contain only parameters?

For example, say at site 1 we have the matrix

-α    α
α+β -α-β

at site 2, we have

-α    α
 α   -α

and at site 3 we have

-β-1  β   1
 β   -2β  β
 1    β  -β-1

Since the values assigned to α and β are assumed to hold over all sites, will the first two substitution models work as they are, or will normalization still pose a problem?

[rbouckaert]

The way you have things set up looks OK to me, and indeed leaving at least one rate fixed (to 1) makes everything identifiable.

What I tried to point out is that in a rate matrix for N states, there are N*N entries in the rate matrix, but N+1 of them are not free: N on the diagonal and 1 because of the normalisation. So, for binary data, N=2, leaving just 1 degree of freedom. But since you have more sources of information for α and β this restriction does not apply in your case: α and β are identifiably for their individual sites, so having another rate matrix as function of α and β should not be a problem.

Anyway, let me know if you need any help implementing these substitution models or the constraints you want to put on parameters.

[jjmccollum]

I may need your help in implementing these changes. I'm working with @rbturnbull on this, and he's more experienced in adding features to BEAST than I am, but it would also be a valuable skill for me to learn, so any help would be appreciated!

[rbturnbull]

Hi @jjmccollum and @rbouckaert - this all looks good to me. I'm happy to work with @jjmccollum on this and we'll see how we go. We might have questions in the future.

[jjmccollum]

It looks like we might already be able to encode rate parameter combinations in substitution model matrices using RPNcalculator. For the matrix

                 -1                      1
        Clar_rate + Byz_rate  -Clar_rate - Byz_rate

the XML

<parameter spec="parameter.CompoundValuable" id="rates.character1" name="rates">
    <!-- Start rate vars -->
    <var idref="default_rate"/>
    <var spec="RPNcalculator" expression="Clar_rate Byz_rate +">
        <parameter idref="Clar_rate"/>
        <parameter idref="Byz_rate"/>
    </var>
    <!-- End rate vars -->
</parameter>

appears to work; at the very least, it is not raising any exceptions when I input the XML to BEAST.

[jjmccollum]

Just adding some thoughts on this process in case anyone has any relevant remarks. Substitutional model matrices with arbitrary combinations of variable rates work for short runs (~100,000 samples), but in longer runs (~20,000,000 samples), the process eventually reaches an assignment of rates that makes at least one site's substitution matrix singular (which happens if the site substitution models have @spec="GeneralSubstitutionModel") or ill-conditioned, resulting in repeated "SVD is not converged" warnings and eventually a java.lang.ArrayIndexOutOfBoundsException in the RobustEigenDecomposition.tql2 method (which happens if the site substitution models have @spec="ComplexSubstitutionModel"). It makes sense that this happens given how arbitrary the rate matrices are allowed to be, but it would be nice if I could force a perturbation of the rates whenever they result in a singular or ill-conditioned matrix or if I could adjust the rates together with an operator in a way that avoids this outcome (although I'm not sure if that's possible).

[rbouckaert]

@jjmccollum would putting a prior on rates to ensure there are not many (almost) zero rates in the matrix solve the problem? Could that be a feasible solution for your analysis?

[jjmccollum]

@rbouckaert It looks like that did the trick! With the priors on rates, I'm no longer seeing the same errors. Thanks for the suggestion!