Handling outputs from the same model with different likelihoods / score functions

[ben18785]

PKPD models will likely need outputs from the same model with different likelihoods. E.g. one output follows a normal model; another follows a log-normal. At the moment, this sort of thing is not possible in PINTS without lots of manual hacking. Similarly, it is not possible to natively handle an inference problem with multiple outputs where those outputs are measured at different times.

I wanted to open the forum on this. Here's a potential proposal.

We create a MasterProblemCollection with the following methods:

• initialised with a forward model and values as is usual. It is also (crucially) initialised with some way of indexing which outputs correspond to which SubProblem (see below). At a simple level, this could be a list [1,1,1,2,2,2,2,3,3,3,3,3] which would indicate that the first three forward model outputs correspond to subproblem1, the second four to subproblem2 and the last five to subproblem3. It is also initiated with a list of times lists: one list of times for each SubProblem.
• evaluate takes an argument index (in the above example, this would be either 1, 2 or 3) which specifies which output set to return. The first time evaluate is called, the result is cached so that the model is only solved once -- not once for each output set.
• n_outputs also takes index and returns the relevant number of outputs for the SubProblem
• times takes index and returns the appropriate time list

We also create a SubOutputProblem:

• intialised with a MasterProblemCollection and an index
• evaluate calls MasterProblemCollection.evaluate(index)
• and so on for n_outputs etc

Example use case

User has a forward model with three outputs and would to model the first two using a GaussianLogLikelihood and the third with LogNormalLogLikelihood.

They would do:

times_1 = [1, 2, 4, 6]
times_2 = [3, 5]
times = [times_1, times_2]
index_list = [1, 1, 2]
master = MasterProblemCollection(forward_model, times, value, index_list)

subproblem_1 = SubProblem(master, 1)
subproblem_2 = SubProblem(master, 2)

loglikelihood_1 = GaussianLogLikelihood(subproblem_1)
loglikelihood_2 = LogNormalLogLikelihood(subproblem_2)
logprior_1 = ...
logprior_2 = ...

logposterior_1 = LogPosterior(loglikelihood_1, logprior_1)
logposterior_2 = LogPosterior(loglikelihood_2, logprior_2)

The user could then create a wrapper to create a LogPosterior class that just returns the sum of logposterior_1 and logposterior_2 which they could use for inference.

@martinjrobins @MichaelClerx @chonlei (and anyone else!) interested to hear thoughts on this

[MichaelClerx]

PKPD models will likely need outputs from the same model with different likelihoods. E.g. one output follows a normal model; another follows a log-normal.

I don't understand what you mean here. The PKPD model would be an ODE model, right? So where do the distributions come in?

Similarly, it is not possible to natively handle an inference problem with multiple outputs where those outputs are measured at different times.

Just for my intuition, can you give me examples of why you'd want or need to do this?

[martinjrobins]

just a note to say we will definitely need this for the pkpdapp inference. Different outputs of the ODE model will have different noise models, hence the need for different distributions on different outputs. The different outputs might also be measured at different times, e.g. if you are measuring the concentration of the drug as well as the size of a tumour, these wouldn't neccessarily be at the same time.

[ben18785]

Thanks @MichaelClerx

Yes, PKPD models are just ODE models. What I mean is that, for these types of model, when you want to do inference, it seems to often be the case that some of the outputs would be modelling using Gaussian distributions; others with log-normal distributions.

Re: different outputs having different measured times, I can give you a good example from epidemiology but there are lots of cases in PKPD modelling it seems due to non-regular measurement of different things (e.g. tumour size and inflammation markers). For SEIR-type modelling, it is often the case that some outputs (e.g. cases and deaths) are measured regularly (e.g. daily) whereas others are measured less regularly (e.g. the proportion of population who are seropositive).

[MichaelClerx]

Thanks both! Worth keeping #1372 in mind too I guess

[MichaelClerx]

Caching sounds useful whether or not we have a split analysis of the model output, so... separate ticket?
I'm assuming you want to evaluate all times in a single simulation? Then I'd say there's two obvious ways:

1. Like you suggest, construct some kind of SuperProblem from SubProblems (let's shy away from "master" for now?). Then this problem would combine the requirements of individual subproblems, run a single sim (or retrieve a cached output), and then split the output up again. E.g. it would have to combine your times_1 = [1, 2, 4, 6] and times_2 = [3, 5] into times=[1,2,3,4,5,6] and then be clever enough to split the resulting data sets again.

2. Instead, we could also consider filtering. E.g. have a class ProblemFilter(problem, filter) where filter is something like a boolean array of the same size as the expected output, or a lambda that the user provides. (Or to make it friendlier you could have explicit filters for times and for outputs, but that wouldn't be quite as flexible?) Then you'd (1) Create a MultiSeriesProblem as usual, (but set caching=True) (2) Create a bunch of filters, (3) Create a composed likelihood (as in your example)

[ben18785]

Yep, spot on with what I was thinking re: 1 and, yes, happy to shy away from "master"!

Re: filters, I'm not sure I get ya? Also, currently an issue with MultiSeriesProblem is that this doesn't allow different time arrays for each output, so I guess that'd need to be re-engineered?

[ben18785]

@MichaelClerx I'm going to give my suggestion a go here as need to make progress for PKPD app. But, look forward to your thoughts on the PR!

[MichaelClerx]

Yep, spot on with what I was thinking re: 1 and, yes, happy to shy away from "master"!

Re: filters, I'm not sure I get ya? Also, currently an issue with MultiSeriesProblem is that this doesn't allow different time arrays for each output, so I guess that'd need to be re-engineered?

Hold on! Can I try explaining the filter thing a bit?

[ben18785]

Sure!

[MichaelClerx]

Basically, in your proposal you start by creating a superproblem that then needs to combine some sub problems, run a sim, and split its output again.
I'm wondering why you don't cut out the first step, and just create methods to split a problem's output into sub outputs?

So you'd have an ordinary MultiSeriesProblem, but then you'd make sub problems based on its output, and finally combine those sub problems into a single likelihood (although you could hide output as well, using this method)

[MichaelClerx]

All the splitting code would need to do would be e.g.

def evaluate(parameters):
    all_values = parent_problem.evaluate(parameters)
    selected_values = all_values[self._numpy_filter]
    return selected_values

or

def evaluate(parameters):
    all_values = parent_problem.evaluate(parameters)
    selected_values = self._very_flexible_filter(all_values)
    return selected_values

[MichaelClerx]

The thinking behind this is:

• we've already got ways to combine problems into a single function
• having it work with any old multi (or single) series problem is more flexible than requiring a special type
• this system would potentially allow for more complicated filters than in your original proposal (the user just passes in either a numpy filter or a lambda function)
• the efficiency comes from using caching in the problem class separately, so the code we write to do that will be useful to others as well

[ben18785]

Ok, thanks. So would the MultiSeriesProblem take all the times where all outputs were evaluated? E.g.:

times_1 = [1, 2, 4, 6] # output 1
times_2 = [3, 5] # output 2
times_all = [1,2,3,4,5,6]
problem = MultiOutputProblem(model, times_all, values_all)

# construct sub problems which are then used to construct individual log-likelihoods
subproblem_1 = SubProblem(problem, filter_1)
...

How would you handle values_all? In the current MultiOutputProblem this needs to be a numpy array of size n_times by n_outputs. But in the inference problems I'm thinking about we won't have observations for all times for all outputs.