still confused about setting up a negative binomial likelihood

[bbolker]

Searching the pkgdown site for "neg_bin" shows me:

• "specifying likelihood and prior components", which tells me I need to use mp_neg_bin, and points me to the SHIVER example. Searching this page for "nbinom" and "neg_bin" doesn't find anything. I also can't find anything in the tmb.R model definition file
• the Distributions man page, which says

  All distributional parameter arguments can be specified either as a numeric value, a character string giving the parameter name, or a distributional parameter object (See fit_distr_params).

• the fit_distr_params page also doesn't say anything explicit about neg_bin or nbinom: the only examples there are using mp_normal() to specify Gaussian priors on the transmission parameter, which is a sufficiently different use case that I'm having a hard trouble adapting it.

If I use traj = list(I = mp_neg_bin(disp=phi) I get "object 'phi' not found". If I use disp="phi" (or disp="log_phi") instead I get

Error in $<-.data.frame(*tmp*, "time", value = numeric(0)) :
replacement has 0 rows, data has 67

(I used to use options(macpan2_default_loss = "neg_bin") but that's deprecated now ...

grepping through the macpan2 code base I see the following types of usage:

./tests/testthat/test-calibrator-traj.R:        I = mp_neg_bin(disp = 0.5),

(makes sense, but I want to estimate rather than use a fixed value)

./tests/testthat/test-distributions.R:    , traj = list(I = mp_neg_bin(disp = c(1,5)))

(OK, this is an error-checker, makes sense)

./tests/testthat/test-distributions.R:    , traj = list(I = mp_neg_bin(disp = mp_fit(2)))

I guess this is what I want to do, but what does mp_fit(2) do ...??

From ?mp_fit:

x: numeric starting value of the distributional parameter to fit, or character name of an existing variable in the model with a default starting value to use.

OK, it makes sense to use a starting value of 2, but how is the name of this parameter decided? Magic? I still get the "replacement has 0 rows, data has 67" error, but maybe there's something else wrong with my code?

I know I could try to start from something closer to the standard calibration example to make sure that I'm not doing something else wrong (minimal changes), but right now I'm out of steam.

---

My example so far:

## mapping/profiling
library(macpan2)

sir <- (mp_tmb_library("starter_models", "sir", package = "macpan2")
  |> mp_tmb_insert(default = list(beta=70/52, gamma=60/52,
                                  N=40000, I=1, phi=6))
)

## convert data to macpan2 format
sirdat <- fitode::SierraLeone2014 |>
        dplyr::transmute(
          matrix = "I",
          row = 0,
          col = 0,
          value = confirmed)

sir_calibrator = mp_tmb_calibrator(sir
  , data = sirdat,
  , traj = list(I = mp_neg_bin(disp=mp_fit(2))),
## also tried other variants here like   par = c("beta", "gamma") ...
  , par = c("log_beta", "log_gamma", "log_I")
)

[RichardSichengZhao]

This might be related, but seems outdated:

https://canmod.github.io/macpan2/articles/calibration_advanced.html?q=binom#a-slightly-harder-more-realistic-example

[RichardSichengZhao]

Adding time = c(1:67) in the dataframe seems fixed the "replacement has 0 rows, data has 67" error.

[bbolker]

thanks @RichardSichengZhao :

• unfortunately your first reference is indeed outdated (I wrote that machinery originally; it's been superseded by other stuff now, but I can't figure out the new stuff ...)
• your second hint was very valuable. I've added some checks to produce a more useful error message in this case ... mp_fit(2) may indeed have worked
• my point about this being hard to find and understand still stands ...

[stevencarlislewalker]

thanks @RichardSichengZhao :

* unfortunately your first reference is indeed outdated (I wrote that machinery originally; it's been superseded by other stuff now, but I can't figure out the new stuff ...)

I'll work more on docs.

* your second hint was very valuable. I've added some checks to produce a more useful error message in this case ... `mp_fit(2)` may indeed have worked

It does FWIW.

* my point about this being hard to find and understand still stands ...

Agreed. I'll work on docs.