GP prior for growth rate

[sbfnk]

I had a look at whether a model with a GP prior on the growth rate (instead of R) would perform better (in particular: easier/faster to fit as it avoids calculations around the generation time).

It seems to work but no discernible improvement in sampling - not sure if worth keeping:

> ease_of_fit
# A tibble: 5 x 3
  gp_on      mean_time divergent_transitions
  <chr>          <dbl>                 <int>
1 infections      7.37                     0
2 gr             55.7                      0
3 R0             63.6                      0
4 R_t-1         121.                      19
5 gr_t-1        179.                      96

[seabbs]

I am a bit lost in the changes here but I thought I had already nearly done this by placing a prior on I_t = GP*I_(t-1)? I guess the major difference is the log scale and that above has no clear interpretation for the GP effects.

(

EpiNow2/inst/stan/functions/infections.stan

Line 73 in d558305

infections[i] = infections[i - 1] * exp_noise[i];

)

Just looking at the branch.

Looked at the branch and I see more what your getting at now. I think it makes sense and is worth exploring. When I looked at the model above I was quite surprised there wasn't more speed up as it seems like there really should be.

My plan (#206) was to generalise update_Rt to be for covariates and then use it for both R -> infectiousness and back calculation which I think might give you the changes here but perhaps in a slightly cleaner way.

[sbfnk]

The way I thought about this is that all the models have the same relationship between I, Rt and r:

where the models differ only in the choices of GP prior and where it applies (in very simple terms, what is being smoothed) - available models being (where C is the number of cases and me the mean shift):

1. I[t] = C(t + m) * GP
2. I[t] = GP
3. I[1] = C(1 + m) * GP[1]
   I[t>1] = I(t - 1) * GP[t]
4. R[t] = R0 * log(GP) * log(sum of bp)
5. R[1] = R0 * log(GP[1]) * log(sum of bp)
   R[t>1] = R[t - 1] * log(GP[t])

to which I thought it would make sense to add (in order to avoid calculations with / dependency on the generation time):

6. r[t] = r0 + GP + (sum of bp)
7. r[1] = r0 * GP[1]
   r[t>1] = r[t-1] * GP[t]

What I hadn't realised is that model (3) was already implemented - which is indeed very similar to (6), differences being the centering of the prior around r0 based on given R0 prior, exponentiation of the GP and availability of break points (which are trickier to interpret than with R i think).

Completely agree that adding covariates would be a great extension (in a flexible way, i.e. including arbitrary/nonlinear interactions as in brms), and that backcalculation and generating infectiousness are really the same model (given by the equations above + observation model) so it should be possible to streamline/simplify.

[seabbs]

Yes totally agree. Really nice to have the model family written down. The other extension would be to add more links. An obvious one being the logit (so that R0 has an upper bound). I was thinking of doing this in two stages.

1. Expanding update_rt to be used for Rt, back calculation and also growth rates (this would involved doing some other standardisation bits and bobs + use code from this PR). It would be good to do this in order to make the relationships you have outlined above obvious and to reduce to amount of code that is copied.
2. Adding covariates. I am still looking to what is the best thing to do for this. An easy option is to implement from scratch (well using the stan manual) but that would just give you linear/polynomial relationships. Looking at rstanarm, brms , etc to see if any bits and pieces can be copied to give more advanced features (splines basically and if ever extend to multiple types of observations other nice to haves).

[sbfnk]

Great - only thought re 2 is that it would be nice if an arbitrary relationship could be used / passed in (e.g. as Stan code) as with the S gene dropout analysis using brms - but I appreciate this might require quite a re-vamp and be beyond the scope of this package (i.e. easier done just in Stan).

[seabbs]

The issue with that is you would need to recompile the model (as in brms) rather than supply it compiled.

I agree would be nice though. Could maybe set it up to work as compiled and optionally add other options you can only use if you recompile.

I was also thinking you could just focus on the custom model family and therefore have something that reproduced much of what EpiNow2 does but as a plug in for brms.