Allow a continuous size (dispersion) in nbinom and nbinomMu

[hiddevandebeek]

Summary

llikNbinom() and llikNbinomMu() require an integer size. In the negative
binomial's mean/dispersion (NB2) parameterisation size is a real dispersion
parameter, not a count -- the pmf uses gamma(size + x) / (gamma(size) * x!), and
stats::dnbinom(x, size, mu) has always accepted continuous size. Continuous
overdispersion is the norm in count models, so this blocks nbinomMu() models in
nlmixr2.

The R-level checkmate::assertIntegerish(size) is only the visible half. size was
also stored in an Eigen::VectorXi and assigned with N(0) = (int)(size), so the
C-API entry points (rxLlikNbinom / rxLlikNbinomMu) -- which rxode2 model solves and
focei call directly, bypassing the R assertion -- silently truncated it.

Current behaviour

R entry point, any non-integer size:

rxode2ll::llikNbinomMu(2L, 0.7, 5)
#> Error: Assertion on 'size' failed: Must be of type 'integerish',
#>        but element 1 is not close to an integer.

Model/solve path, which never reaches that assertion:

m <- function() { model({ ll <- llikNbinomMu(time, size, mu) }) }
et <- et(2:3); et$size <- 3.9; et$mu <- 5
rxSolve(m, et)$ll
#> -2.090736 -2.049914      # == dnbinom(2:3, size = 3, mu = 5, log = TRUE)
dnbinom(2:3, size = 3.9, mu = 5, log = TRUE)
#> -2.113954 -2.014227      # what it should have returned

No warning -- the returned value is the log-likelihood at trunc(size).

For 0 < size < 1 the truncation gives size = 0, Stan's neg_binomial_2_lpmf
throws, and the exception crosses the extern "C" boundary uncaught:

et$size <- 0.5
rxSolve(m, et)
#> terminate called without an active exception
#> [R process aborts]

Because the truncation quantises size, the log-likelihood is a step function of the
dispersion -- its derivative is zero almost everywhere:

k <- seq(2, 4, by = 0.25); et <- et(rep(3, 9)); et$size <- k; et$mu <- 5
rxSolve(m, et)$ll
#> -2.128648 -2.128648 -2.128648 -2.128648 -2.049914
#> -2.049914 -2.049914 -2.049914 -2.011349
# 9 distinct k -> 3 distinct log-likelihoods; 6 of 8 slopes exactly 0

This is why the dispersion cannot be estimated even when its true value is an
integer: any optimiser has to move k through non-integer values, and there it sees
no gradient.

Downstream, a plain nlmixr2 count model with a dispersion parameter cannot be fit:

mod <- function() {
  ini({ lmu <- 1.2; lk <- 0; etaMu ~ 0.2 })
  model({
    mu <- exp(lmu + etaMu)
    k <- exp(lk)
    dv ~ nbinomMu(k, mu)
  })
}
nlmixr2(mod, dat, "focei")
#> Error: Could not fit data
#>   Could not find the best eta even hessian reset and eta reset for ID 1.

Across a sweep of simulated data sets (100 subjects, 10 observations each; true k
in {0.3, 0.7, 1.5, 3}, four replicates each, plus three nbinom() size/prob fits),
0 of 19 fits succeed on the current release -- including every replicate whose
true k is exactly 3, for the reason above.

Changes

• src/llikNbinom.cpp, src/llikNbinom2.cpp: store size in an
  Eigen::VectorXd instead of an Eigen::VectorXi (struct field, constructor and
  llik_nbinom() / llik_nbinomMu() signatures), and assign N(0) = size without
  the (int) cast. Stan's neg_binomial_2_lpmf already accepts a real precision
  argument, so nothing else in the likelihood or its derivative changes.
• Same files: the bounds guard size < 0.0 || size > INT_MAX becomes size <= 0.0.
  The INT_MAX bound is meaningless once size is a double, and NB2 requires
  size > 0; non-positive values return NA_REAL rather than reaching Stan and
  aborting the process. The guard on x is unchanged -- x is a genuine count.
• R/llik.R: assertIntegerish(size, ...) becomes
  assertNumeric(size, lower = 0, finite = TRUE) in llikNbinom() and
  llikNbinomMu(). assertIntegerish(x) is kept.
• llikBinom() is deliberately untouched: there size is a number of trials and
  the integer restriction is correct.

Tests

Two existing tests asserted the old behaviour and are updated:

• llikNbinomInternal returns NA for size > INT_MAX
• llikNbinomMuInternal returns NA for size > INT_MAX

size = 2^31 is a valid dispersion; it was only out of range because it had to fit
in an int. Both now assert the value equals stats::dnbinom(...), which it does
exactly. They are joined by new size <= 0 returns NA tests.

New regression tests cover continuous size for both parameterisations: values
against stats::dnbinom(), and dMu / dProb against central finite differences.

Verification

• fx matches stats::dnbinom(x, size, mu = mu, log = TRUE) to a maximum absolute
  error of 3.6e-15 over a 112-point grid of x in 0:6, size in
  {0.3, 0.7, 1.5, 4.25} and mu in {0.5, 2, 5, 20}; the size/prob form matches
  stats::dnbinom(x, size, prob = prob, log = TRUE) equally well.

• dMu and dProb match central finite differences of the continuous-size log
  density to ~1e-9.

• The documented integer-size examples (llikNbinomMu(46:54, 100, 40),
  llikNbinom(46:54, 100, 0.5)) are unchanged.

• The solve path no longer aborts for 0 < size < 1, returns the continuous-size
  value to 8.9e-16, and no longer agrees with the truncated-size reference.

• tests/testthat/test-llik.R: 44 blocks, 68 expectations, 0 failures, 0 errors.

• The log-likelihood is once again strictly monotone in the dispersion: over the same
  k grid that produced 3 distinct values, all 9 are distinct, no slope is zero, and
  each matches stats::dnbinom() exactly.

• The focei sweep described above -- same seeds and data that fail 19 out of 19 on
  the current release -- converges 19 out of 19 against the patched build, all with
  relative convergence (4). Mean estimated k over four replicates:

  true k	mean estimate	range	relative bias
  0.3	0.308	0.284 - 0.330	+2.5%
  0.7	0.746	0.695 - 0.797	+6.6%
  1.5	1.501	1.359 - 1.566	+0.1%
  3.0	3.038	2.730 - 3.398	+1.3%

  The three nbinom() size/prob fits recover 0.611, 2.399 and 2.959 against true
  values of 0.7, 2.5 and 3.0.

Built and tested on Windows, R 4.5.3.

Note on NEWS.md

The NEWS.md change is purely additive, under the current 2.0.16 section; the
released 2.0.15 section is untouched. Its INT_MAX bullet lists size alongside
x for all five discrete distributions, which after this change is accurate only for
rxLlikBinom -- rather than edit a shipped entry, the new bullet states where the
INT_MAX bound still applies.

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[mattfidler]

Review summary

Reviewed by Claude and independently by Antigravity (Gemini 3.1 Pro). Both rounds were run against a real build — this environment needed RcppParallel 6.0.0 installed from source first, since 2.0.16 now requires it.

The core change is right. size in NB2 is a real dispersion parameter, storing it in an Eigen::VectorXi was truncating it, and neg_binomial_2_lpmf() already takes a real precision argument, so VectorXd is the correct fix and nothing downstream in the likelihood or its derivative needs to change. I confirmed the values against stats::dnbinom() independently.

Findings

Both reviewers converged on one theme: the mean handed to Stan can still be out of domain, and since rxLlikNbinom() / rxLlikNbinomMu() are extern "C" with no handler between them and the solver, the resulting exception doesn't become an R error — it aborts the process. That's the same failure mode this PR set out to fix for size < 1.

1. New in this PR. Dropping the size > INT_MAX bound opened a case the old bound had been incidentally covering: in the size/prob form mu = size*(1-prob)/prob can now overflow to Inf for a large finite size with a small prob. llikNbinom(1L, 1e300, 1e-10) reached Stan and threw; on main the INT_MAX bound returned NA.

2. Pre-existing, same functions. prob == 0 (mu = Inf), prob == 1 (mu = 0), prob outside [0, 1] via the C API (mu < 0), and a non-positive mu in the mean/dispersion form all threw as well. All four are reachable — the R-level assertions use lower = 0, upper = 1, so llikNbinomMu(1L, 10, 0) and llikNbinom(1L, 10, 0) both get through. Worth noting these are the outliers: llikBinom/llikGeom/llikPois all sanitise their degenerate parameter values.

Fixed in b00650d by guarding the derived mean rather than re-bounding size — one check covers the overflow and every boundary case, and keeps size genuinely unbounded, which was the point of the PR. Non-positive mu guarded in the same way. All now return NA with tests for each.

Notes, not blocking

• Returning NA at prob == 1 / mu == 0 deviates from stats::dnbinom(), which is defined there (1 at x == 0, 0 otherwise). These are degenerate points with no usable derivative and Stan can't evaluate them, so NA beats aborting. Recorded in NEWS.md.
• Agreement with stats::dnbinom() degrades with size — ~5e-14 absolute at size = 1e4, ~4.7e-10 at size = 2^31 — from Stan's large-phi Poisson branch, not from this change. At size = 2^31 that's 1.9e-10 relative, and a dispersion that large is Poisson anyway. Not a problem, but it's why the new size > INT_MAX tests want a tolerance rather than exact equality.
• The NEWS.md handling of the shipped 2.0.15 INT_MAX bullet is the right call — clarifying in a new entry rather than rewriting a released one.

Verification

• testthat: 74 passing, 0 failures, 1 skip (the ~103 GB R_xlen_t test).
• R CMD check: 0 errors, 0 warnings. The single NOTE is a stray untracked directory in my working tree, not from this PR.
• Round-two review by both reviewers over the cumulative diff: no findings.

Looks good to merge from my side. The one thing I'd flag for @mattfidler is timing rather than code: CRAN-SUBMISSION points at 3c6339b, so this is a behaviour change landing on top of an already-submitted 2.0.16.