Issue 118: First pass at missing data support

[seabbs]

This PR adds support for missing data in the observation model by moving from a vectorised to loop approach. It also renames the negative binomial observation model to match the distribution naming scheme.

[seabbs]

@SamuelBrand1 any thoughts here? The issue I had here was in trying to give both a vectorised and non-vectorised likelihood option but perhaps just accepting it can't be vectorised is the way forward for now?

It all seems surprisingly clunky to me as well so perhaps I have just got the wrong end of the stick/.

[SamuelBrand1]

@SamuelBrand1 any thoughts here? The issue I had here was in trying to give both a vectorised and non-vectorised likelihood option but perhaps just accepting it can't be vectorised is the way forward for now?

It all seems surprisingly clunky to me as well so perhaps I have just got the wrong end of the stick/.

Yeah, the underlying maths here is that vectorised == multivariate distribution and non-vectorised == conditionally independent.

Having a vector with element type Union{Missing, Integer} doesn't make much sense to also use arraydist, whereas we are asking it to be multivariate and conditionally independent at the same time. The reason we want that is because we believe that the performance of the ReverseDiff AD is better that way.

So I'd suggest just trying out a non-vectorised approach and checking if it really does have a big performance hit. The AD systems evolve independent of the PPL and so what might have been true when they wrote their performance tips might not be true now. If it doesn't then a loop is more flexible, missing elements will just work as expected.

[SamuelBrand1]

Another option is to make the concept of missing part of the generative process; e.g. there is a probability model for observed cases (we have that) and another that filters that and unobserved cases on a day are given value -1 (for example, and maintaining common eltype in what I imagine is an Integer array).

Then the usage would be missing -> -1 in data wrangling.

[seabbs]

Another option is to make the concept of missing part of the generative process; e.g. there is a probability model for observed cases (we have that) and another that filters that and unobserved cases on a day are given value -1 (for example, and maintaining common eltype in what I imagine is an Integer array).

Then the usage would be missing -> -1 in data wrangling.

I really don't love this and it feels very clunky but I see it would work and enable vectorisation.

[seabbs]

Note I have had to update the getting started examples to use generated_y_t vs y_t which I think is an internal Turing approach to account for having both observed data (i.e. y_t) and generated data (i.e generated_y_t) which can now happen due to the support for partial missingness?

[SamuelBrand1]

Do you have any sense that this is noticeably bad for inference time?

[seabbs]

the only real evidence I have for this is running the getting started example. Comparing CI runs there is perhaps a small slowdown (https://github.com/CDCgov/Rt-without-renewal/actions) but it doesn't seem massive (though hard to know how that would scale given so little of that CI check is from inference).

For interest in epinowcast we have this nice benchmark CI that really helps with these kinds of questions: epinowcast/epinowcast#442 (comment)

[SamuelBrand1]

the only real evidence I have for this is running the getting started example. Comparing CI runs there is perhaps a small slowdown (https://github.com/CDCgov/Rt-without-renewal/actions) but it doesn't seem massive (though hard to know how that would scale given so little of that CI check is from inference).

For interest in epinowcast we have this nice benchmark CI that really helps with these kinds of questions: epinowcast/epinowcast#442 (comment)

Very nice. We can do that with BenchmarkTools I think

[SamuelBrand1]

Another option is to make the concept of missing part of the generative process; e.g. there is a probability model for observed cases (we have that) and another that filters that and unobserved cases on a day are given value -1 (for example, and maintaining common eltype in what I imagine is an Integer array).

Then the usage would be missing -> -1 in data wrangling.

I really don't love this and it feels very clunky but I see it would work and enable vectorisation.

We could have the not missing (or missing indices) as input-able data?

[SamuelBrand1]

@seabbs The example script doesn't use a missing / Int mix.

I've tried that out here, but then generated quantities starts to chuck an error. So sadly this doesn't work yet :-(.

[seabbs]

The example script doesn't use a missing / Int mix.

That is why I was using it for benchmarking as it should show the difference in speed for complete data (which is what we care about for the benchmark).

I tested the generated quantities in tests but yes ideally we would have a small example to had partially complete data at some point (an issue and not a priority for now I think).

[seabbs]

e tried that out here, but then generated quantities starts to chuck an error. So sadly this doesn't work yet :-(.

Where did you try this? Shouldn't it be here:

Rt-without-renewal/EpiAware/docs/src/examples/getting_started.jl

Line 240 in 02fcf74

truth_data = generated_obs

[SamuelBrand1]

e tried that out here, but then generated quantities starts to chuck an error. So sadly this doesn't work yet :-(.

Where did you try this? Shouldn't it be here:

Rt-without-renewal/EpiAware/docs/src/examples/getting_started.jl

Line 240 in 02fcf74

truth_data = generated_obs

Ooops I hadn't committed my change

[seabbs]

I added a small example here as part of the getting started example: e28cc8c

[SamuelBrand1]

I added a small example here as part of the getting started example: e28cc8c

Nice one. LGTM now.

[SamuelBrand1]

Very nice.