Check log-likelihood calculations and interpretation

[ArtPoon]

I obtained these results on a simulated tree with 500 tips where there was assumed to be two risk groups with a transmission rate discordance by a factor of 3.

> res <- clmp(t1, nrates=1)
log likelihood for 1 state model is 898.696965
rates: 185.511109 
Q: [    *    ]
> res <- clmp(t1, nrates=2)
log likelihood for 2 state model is 903.303757
rates: 170.256368 481.467315 
Q: [    *   2.990365 ]
   [ 6.715851   *    ]
> table(res)
res
  0   1 
904  79 
> res <- clmp(t1, nrates=3)
ConditionNumber: maximal condition number 9.01e+12 reached. maxEW=1.87e-02,minEW=1.93e-15,maxdiagC=1.16e-02,mindiagC=2.48e-15
log likelihood for 3 state model is 903.801291
rates: 74.498138 180.094874 484.310980 
Q: [    *   0.00000046.299770 ]
   [ 7.564244   *   0.000000 ]
   [ 0.0000008.826223   *    ]
> table(res)
res
  1 
983 

[ArtPoon]

• I would expect log-likelihood values to be negative.
• The likelihood for a model a greater number of parameters should be tend to be higher, thus the log-likelihood should be a smaller negative value (closer to 0).

[ArtPoon]

Using the examples/test.nwk tree, I get the following:

> clmp(t1)
log likelihood for 2 state model is 2238.543290
rates: 495.368077 1305.115690 
Q: [    *   2.526692 ]
   [ 23.309447   *    ]
> res0 <- clmp(t1, nrates=1)
log likelihood for 1 state model is 2227.311795
rates: 531.823491 
Q: [    *    ]

[ArtPoon]

Log-likelihood calculation is being rescaled in mmpp.c:likelihood():

        new_scale = get_scale(&w->L[i * nrates], nrates);
        for (pstate = 0; pstate < nrates; ++pstate) {
            w->L[i * nrates + pstate] /= pow(10, new_scale);
        }
        w->scale[i] += new_scale;
    }
    // from utils.c
    lik = (reconstruct ? max_doubles : sum_doubles)(&w->L[rt * nrates], nrates);
    return log10(lik) + w->scale[rt];

Modify to return scaling constant?

[ArtPoon]

It looks like scaling constant is always set to 0 if nrates=1, so I'm not sure the above is a satisfactory answer...