R summary estimate using bpmodels

[adamkucharski]

Given data on early transmission clusters, would be useful to have an initial estimate of R with uncertainty, as well as an estimate of P(R>1). Below gives a simple implementation using a profile likelihood with quickfit (there may be a more stable efficient package for this, however).

Calculating P(R>1) would require doing the inverse inference, i.e. updating quickfit (or another package) to find the level of confidence that gives an upper bound corresponding to R=1.

#devtools::install_github("adamkucharski/quickfit")
#devtools::install_github('sbfnk/bpmodels')
library(quickfit)
library(stats4)
library(bpmodels)

# Early clusters
cluster_sizes <- c(rep(1,10),5,10)

# Estimate R and k with quickfit
log_l <- function(x,a,b){
  bpmodels::chain_ll(x,"nbinom",stat="size",
                     mu=a,size=b) 
}

estim <- quickfit::calculate_profile(log_l,
                            data_in=cluster_sizes,
                            n_param=2,
                            a_initial = 0.5,
                            b_initial = 1,
                            precision=0.01)

# Reformat and output CI for R:
estim <- data.frame(t(unlist(estim)))

c(estim$estimate.a,estim$profile_out.a1,estim$profile_out.a2)

[adamkucharski]

This is perhaps better placed as a vignette with {bpmodels}? Or an extension example linked to issue #16 ?

[jamesmbaazam]

It does seem like the core engine here is from {bpmodels}, so better linked to the original package. #16 is a full pipeline and that's where I get confused about whether it belongs in {episoap}? Might also be useful to use the early/middle/late pipelines categorizations to determine where vignettes should live.

[joshwlambert]

I agree that this example would be a better fit for a vignette in {quickfit} or {bpmodels}.

@jamesmbaazam you make a good point about the pipeline vignette and {episoap}. I think we should aim to port a version of the pipeline over to {episoap} once we're happy with it.