paul-buerkner / brms

brms R package for Bayesian generalized multivariate non-linear multilevel models using Stan
https://paul-buerkner.github.io/brms/
GNU General Public License v2.0
1.25k stars 177 forks source link

Option to not loop over observations when using custom_family #1084

Closed seabbs closed 3 years ago

seabbs commented 3 years ago

Relates to #1075 and this discussion.

Firstly, thanks for a great package. Its been a great help in what has been a very busy year!

We have been using stan directly for our Covid-19 related work estimating the effective reproduction number (methods) but have recently been exploring using a wrapper around brms for some related work estimating the case fatality ratio. This currently uses the following model:

D_{i,t} \sim \mathrm{NB}\left(c_i\sum_{\tau = 0}^{30} \xi_{i}(\tau | \mu, \sigma) C_{i, t-\tau}, \phi \right)

where,

c_{i,t} = \mathrm{logit}^{-1}\left(\alpha f_{it} + \beta C_{i,t}+ \gamma_{i}\right)

and c_{i, t} is defined to be mu as estimated by brms. This is all still very preliminary and this particular model may be able to be replaced with something more standard but the issue we have been running into is because of the loop over observations the code is a little convoluted (as the convolution needs to be calculated over the vector of cases). We have also been unable to define the log_lik within the brms framework which is a shame.

I am a big fan of reducing the amount of custom code I and others write and would like to build more on your work if possible (for example I have been looking at implementing a regression module into our reproduction number estimation package which I would really like to avoid!). I've looked at the code for brms and not made any progress with working out what changes would be needed but am happy to have a go with some guidance. More generally, if the idea of a wrapper package around brms with some support for these kinds of infectious disease models is something you would be interested in supporting with some guidance on how best to integrate with brms that would be excellent. Whilst the large scale changes suggested in #1075 might help for more complex models I think quite a lot can probably be achieved with the current implementation + potentially a few tweaks.

Unrelated but of use is ideally we would be able to output mu as a tracking tool. Any thoughts on that would be appreciated (I think this is already very doable but I have not managed to think about it correctly yet).

paul-buerkner commented 3 years ago

I will think of a good option.

Can you elaborate on your last sentence please? I didn't quite get it.

paul-buerkner commented 3 years ago

This is now supported in the github version via argument loop of custom_family.

seabbs commented 3 years ago

Amazing thanks for resolving this (and sorry for the lack of reply - sleep deprivation!).

To clarify the last sentence - the issue I was having is reconstructing mu as an output where y ~ mu * cases, mu ~ beta * x. If mu is something we are interested in (so here for example the case fatality ratio), x is some set of covariates and beta are their effect modifiers. I can open a new issue for this but I assumed I just hadn't thought about it enough.

For interest, some of the work we have been doing using the custom_family support in brms is now published (though this CFR extension just ended up as gov reports for now but 🤞🏻): https://science.sciencemag.org/content/372/6538/eabg3055

paul-buerkner commented 3 years ago

Thanks for linking to the paper. I am happy that brms was useful for your study!

For your application with regard to mu, I recommend using brms' non-linear syntax to separate terms.