Parameterization of Negative Binomial for Dynamical System Model Calibration/Fitting

Question

I'm studying about the applications of bayesian inference to fitting dynamical systems to observations.

So the model itself is a deterministic SIR model:

$$ f(R_0,D_{inf})=\begin{cases} \frac{dS}{t}&=-\beta S\frac{I}{N} \\ \frac{dI}{t}&=\beta S\frac{I}{N}-\nu I \\ \frac{dR}{t}&=\nu I \\ \end{cases}\tag{1} $$ where $\beta=R_0/D_{inf}$, $\nu=1/D_{inf}$. The parameters we will calibrate/estimate will be $R_0$ and $D_{inf}$.

So the concept is finding the $R_0$ and $D_{inf}$ makes the model outputs, specifically the $I_t$ best fit the data, say $d_t$ and the uncertainty in the fit is captured by some probability distribution, i.e, Poisson:

$$d_t\sim Poisson(I_t)\tag{2}$$

meaning the actual observed time series $d_t$ follows some probability distribution around $I_t$.

---

Problem

In the case of $Poisson$ distribution, it is relatively easy to construct the likelihood for $R_0$ and $D_{inf}$ since all you have to do is to find $R_0$ and $D_{inf}$ that maximizes the likelihood generated from $(1)$ and $(2)$:

\begin{align*} L(R_0,D_{inf}|d_1,d_2,\dots,d_n) &= \prod_{t=1}^nL(I_t|d_t) \\ &= \prod_{t=1}^n\frac{\exp(-I_t)I_t^{d_t}}{d_t!} \end{align*}

then one could just use the following function for getting the likelihood

dpois(
   x = d_t,
   lambda = I_t
)

But how is this done for negative binomial? Generating likelihood we would have

\begin{align*} L(R_0,D_{inf}|d_1,d_2,\dots,d_n) &=\prod_{t=1}^n\frac{\Gamma (d_t+r_t)}{d_t!\Gamma (r_t)}p^{r_t}(1-p)^{d_t} \end{align*} where $I_t=\frac{r_t(1-p)}{p}$ would be the mean of the Negative binomial distribution.

how could I use dnbinom() for in this method of parameterization?

dnbinom(
 x=d_t,
 mu = I_t,
 size = ???
)

I'm slightly confused by the documentation and would like help on what should size be.

Please comment as well if this belongs to R Collective (StackOverflow) rather than CrossValidate rather than a straight downvote.

Answer 1

Checking ?dnbinom, the equivalent parameterisation is: your $d$ is R's x, your $p$ is R's prob, your $r$ is R's size, and your $I$ is R's mu.

mu is the mean, $I_t$ in your case. Then prob = size/(size+mu) (see docs). So

$$ \texttt{size} = r/(r + I_t) $$

according to your parameterisation.