Added censored log-likelihood and some tests for it

[codecov[bot]]

Codecov Report

All modified and coverable lines are covered by tests :white_check_mark:

Project coverage is 100.00%. Comparing base (abcc8e5) to head (9be1b06).

Additional details and impacted files

```diff @@ Coverage Diff @@ ## main #1500 +/- ## ========================================= Coverage 100.00% 100.00% ========================================= Files 97 97 Lines 9577 9665 +88 ========================================= + Hits 9577 9665 +88 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[k-shep]

@DavAug and @rccreswell

I’ve had a go at the censored log-likelihood, which we discussed last year. I appreciate that there are quite a lot of tests, so I thought I’d spell out the cases I’m testing (for a variety of objects):

1. none of the data are censored;
2. some data are lower/upper censored and some aren’t censored;
3. some of the data are lower censored, some are upper censored and some are not censored;
4. all of the data are lower/upper censored.

I’m testing the likelihood and S1 are correct for these cases and the following object types (for the values):

1. list;
2. one dimensional array;
3. two dimensional array single (i.e. a 1D array in a 2D array format);
4. two dimensional array multi (i.e. a ‘proper’ 2D array).

[DavAug]

Ok, just to follow up on my comment before: I had a quick look at the PK modelling reference and your implementation seems to be consistent with method M3 (except, perhaps, a factor of 2, but I am assuming this comes from the difference in the Phi and CDF definitions).

I am still not quite convinced that this is the best method to handle censored data, because it will effectively drive the model to go as far into the censored region as possible to maximise the likelihood. So, for a model with sufficient flexibility, you might get strange results where the model shoots off to minus infinity whenever there is a censored measurement below the censoring threshold (or to infinity for censored measurements above the censoring threshold). That doesn't really seem to be what we want.

Their M4 suggestion would do a little bit better by bounding the range of censored measurements from below for the lower censored measurements (and above for the upper censored measurements) but then you can already anticipate that the Gaussian likelihood will just be maximised when the model positions itself right in the middle between the lower bound and the censoring threshold whenever a censored measurement is present.

I think these shortcomings are very exciting and leave quite a bit of room for some methodological work 😊

[DavAug]

Also one more methodological comment:

Given that a censored measurement contributes with the CDF to the likelihood, as opposed to uncensored measurements contributing with the PDF, the likelihood will be larger the more censored measurements are present (by construction because the CDF is the integral over the PDF).

For the maximum likelihood approach in the PK modelling reference this matters maybe less (although I still don't think it's ideal), for Bayesian inference this will have the consequence that posterior distributions will be narrower than they should be (higher likelihood contributions have a similar effect to having more data). This feels very counterintuitive, as posteriors should probably be wider for uncensored measurements. Just something to think about. 😊

[k-shep]

Thanks @DavAug for the review!

I think I’ve made all the changes you suggested, and the code is much simpler now.

Just as a general comment, I had previously tried to copy the code or coding style from elsewhere in the file, so I think some of your suggestions might also be relevant for that code.

[k-shep]

Ok, just to follow up on my comment before: I had a quick look at the PK modelling reference and your implementation seems to be consistent with method M3 (except, perhaps, a factor of 2, but I am assuming this comes from the difference in the Phi and CDF definitions).

I am still not quite convinced that this is the best method to handle censored data, because it will effectively drive the model to go as far into the censored region as possible to maximise the likelihood. So, for a model with sufficient flexibility, you might get strange results where the model shoots off to minus infinity whenever there is a censored measurement below the censoring threshold (or to infinity for censored measurements above the censoring threshold). That doesn't really seem to be what we want.

Their M4 suggestion would do a little bit better by bounding the range of censored measurements from below for the lower censored measurements (and above for the upper censored measurements) but then you can already anticipate that the Gaussian likelihood will just be maximised when the model positions itself right in the middle between the lower bound and the censoring threshold whenever a censored measurement is present.

I think these shortcomings are very exciting and leave quite a bit of room for some methodological work 😊

Yes – my industrial supervisor recommended using the M3 method, so this is why I started implementing this censored log-likelihood. I think in the reference, they’re multiplying the whole likelihood (both censored and uncensored terms) by -2, although I’m not sure why.

Thank you for your thoughts on this – I think it is potentially an interesting thing to investigate and, if I have time, I’m hoping I might be able to do this later in my PhD. In the papers that I came across when I was looking at this before, I think people dealt with the censored data in different ways – I think I came across the use of the M3 method and also removing the censored data, and I believe when I tried to compare these for simulated data they gave different results in terms of parameter inference (via finding the MLE). There didn’t seem to be a consensus for the best way to deal with censored data, at least in the field I was looking at – I haven’t, however, looked too much into the general theory of this.

From looking at the experiments I tried before, it looks like your thoughts on censoring ‘narrowing’ the likelihood (and posterior) are also correct. The likelihood surface does seem to look flatter when values are not censored, compared to when they are.

[k-shep]

Thanks @DavAug for reviewing this again. I’ve made the changes you’ve suggested. I have changed the documentation, although I’ve written it differently to how you suggested. Would you be OK just checking that you agree what I’ve written makes sense now?

[DavAug]

Thanks @DavAug for reviewing this again. I’ve made the changes you’ve suggested. I have changed the documentation, although I’ve written it differently to how you suggested. Would you be OK just checking that you agree what I’ve written makes sense now?

Hi @k-shep

Of course, no problem! Yes the docstring makes sense and reads very well! Good job 😊

[k-shep]

Thanks @DavAug :) - should I now just click the merge button?

[DavAug]

Thanks @DavAug :) - should I now just click the merge button?

Yes 😊

[k-shep]

Hi @DavAug, I just noticed that I didn't put this log-likelihood class in the correct position for it to be in alphabetical order in that file. Is the best way for me to correct this by opening a new issue and pull request, and making this change?

[DavAug]

Hi,

Good catch! Yes, create an issue for it and a PR when you have changed the order 😊

Thank you!

Best, David

[k-shep]

Thanks @DavAug - I've created the PR for this here: https://github.com/pints-team/pints/pull/1526