Hi Sergei,

Please could you advise on the following problem:

I have specified a substitution model with a 3 category GDD on omega and fitted this in HyPhy (likelihood = -6772.52). I now wish to compare this null model to an alternative model which allows for the largest omega value to differ between cases and controls (different models are allocated to different branches in the tree file).

As a check, I constrained all of the parameters in the alternative model to the values obtained under the null model and computed the likelihood. However, I now get a likelihood of -6793.75 (constrained parameters, no optimising) if I specify:

[code]
Model C = (SubsMatrixC, EFV, 0);
Model V = (SubsMatrixV, EFV, 0);
[/code]

and set all the parameters of SubsMatrixC = SubsMatrixV. On the other hand, if I replace "SubsMatrixV" with "SubMatrixC" in the above code, I get the previous likelihood of -6772.52! What could be going wrong? And which likelihood is correct?

After optimising the unconstrained omegas with all other parameters fixed, I get a likelihood of -6793.42 for the alternative model, which is actually worse than the first likelihood I got with the omegas constrained to be equal for cases and controls!! (but very slightly better than the second likelihood I got for the null model).

I have attached my code so that you can check for any obvious errors. I recall having a similar problem previously where I got different likelihoods from the same constrained Q matrix: Multimedia File Viewing and Clickable Links are available for Registered Members only!!  You need to Login

Thanks!
Miguel

---

Hi Miguel,

If you are using a HyPhy build that is not than ~8 months old, then the calculation for the mixed (along) branches model should be correct. That said, it is possible that by adding a new REL model to the subset of branches can decrease the likelihood of a single REL model.

Consider the calculation of single site likelihood at a site under a two omega REL model, i.e.


Code (]
Pr[omega = A):

= p Pr[omega = B] = 1-p

Hence,


Code (]
L[site):

= L[site | omega = A] * p + L[site | site | omega = B] * (1-p)

Now consider splitting the branches in a tree into two (arbitrary, non-empty) classes, and assigning a separate copy of the same REL model to both classes (hence creating 4 rate classes). This seems like what you are doing with models C and V initially. The likelihood calculation is now


Code (]
L[site):

= L[site | omega = A, omega2 = A] * p^2 + L[site | site | omega = B, omega2 = B] * (1-p)^2 + L[site | omega = A, omega2 = B] * p (1-p) + L[site | omega = B, omega2 = A] * p* (1-p)

Notice that the new likelihood can be LOWER than the previous one, because the two terms corresponding to the single model likelihood have weights adding to LESS than one (p^2 + (1-p)^2), and the other two likelihood terms could DECREASE the likelihood.

HTH,
Sergei

Makes sense! Thanks... Miguel