Hi Sergei,

I am trying to compare to simple models but am getting some strange results. Was wondering if you could suggest why...

The null model is a standard GY94 model - all parameters except omega are estimated once-off over the entire alignment and omega is optimised separately at each site, holding the other parameters fixed. The alternative model defines omega = omega1*f[i] + omega2*(1-f[i]), where f[i] is a weight specific to branch i, and optimises omega1 and omega2 separately at each site. For each site, I compute the null and alternative likelihoods and perform a LRT.

I would expect omega (under the null model) to lie somewhere in between omega1 and omega2 (under the alt model); however, I am not getting this. I typically get omega < omega1 < omega2. (e.g. site 9) Surely a larger omega would then produce a better likelihood under the null model? I've tried larger starting values for omega, but seem to always end up with the same estimate.

I also get some sites where the likelihood actually gets worse under the alt model! At site 45, for example, I start the optimisation with omega1=omega2=0.7515262 (the value of omega under the null model) which has a likelihood of -88.87692. After optimising the alt model, I get omega1=408.816, omega2=2596.717 and a likelihood of -1e+100. 

Any suggestions? I've attached my code and data.

Thanks!
Miguel

Thanks Sergei! I should have picked that up