Allow a continuous size (dispersion) in nbinom and nbinomMu

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: rxode2ll
-Version: 2.0.16
+Version: 2.0.16.9000
 Title: Log-Likelihood Functions for 'rxode2'
 Authors@R: c(
    person("Matthew L.","Fidler",
@@ -14,7 +14,7 @@ Suggests:
 Imports:
     Rcpp (>= 1.0.8),
     checkmate,
-    RcppParallel (>= 6.0.0)
+    RcppParallel
 Description: Provides the log-likelihoods with gradients from 'stan'
     (Carpenter et al (2015), <doi:10.48550/arXiv.1509.07164>) needed
     for generalized log-likelihood estimation in 'nlmixr2'

NEWS.md

@@ -1,3 +1,38 @@
+# rxode2ll (development version)
+
+* `llikNbinom()` and `llikNbinomMu()` now accept a continuous (non-integer)
+  `size`.  In the negative binomial's mean/dispersion parameterisation `size`
+  is a real dispersion parameter, not a count, and `stats::dnbinom()` has
+  always allowed it.  Previously `size` was stored in an integer vector and
+  silently truncated, so a `size` above 1 returned the log-likelihood at
+  `trunc(size)` and a `size` between 0 and 1 truncated to 0 and aborted; the
+  R-level `assertIntegerish(size)` rejected such values outright.  Truncation
+  also made the log-likelihood a step function of `size`, so the dispersion
+  could not be estimated even when its true value was an integer.  This
+  blocked `nbinomMu()` models with continuous overdispersion in `nlmixr2`.
+  `size` must now be strictly positive; non-positive values return `NA`.  The
+  `INT_MAX` bound added in 2.0.15 still applies to `x` everywhere, and to
+  `size` for `rxLlikBinom` where `size` is a number of trials.
+
+* `llikNbinom()` and `llikNbinomMu()` now return `NA` instead of aborting the R
+  process when the mean handed to the underlying likelihood is out of domain.
+  This covers a `prob` of 0 (and, through the C API, a `prob` outside `[0, 1]`
+  or a negative `mu`), and -- because `size` is no longer bounded by `INT_MAX`
+  -- a large `size` with a small `prob` whose product overflows.  Previously
+  these threw a C++ exception that escaped the `rxLlikNbinom()` /
+  `rxLlikNbinomMu()` entry points used by `rxode2` solves and `focei`, killing
+  the session.
+
+* `llikNbinom()` at `prob == 1` and `llikNbinomMu()` at `mu == 0` now match
+  `stats::dnbinom()`.  The distribution is degenerate at zero there, so `fx` is
+  `0` (a likelihood of 1) at `x == 0` and `-Inf` otherwise.  These are the only
+  points where `stats::dnbinom()` defines a value that the underlying Stan
+  likelihood cannot compute -- it requires a strictly positive mean.  No
+  derivative exists at a point mass, so `dProb` / `dMu` are `NA`; an unusable
+  gradient beats a `NA` log-likelihood, and beats the abort these used to
+  cause.  `prob == 0` remains `NA`, matching the `NaN` R returns there.
+
+
 # rxode2ll 2.0.16
 
 * Fix ABI issues in rxode2 and the nlmixr2 ecosystem by requiring

R/llik.R

@@ -88,9 +88,9 @@ llikBinom <- function(x, size, prob, full=FALSE) {
 }
 
 #' Calculate the log likelihood of the negative binomial function (and its derivatives)
-#' 
+#'
 #' @param x  Number of successes
-#' @param size Size of trial
+#' @param size Dispersion parameter; may be non-integer
 #' @param prob probability of success
 #' 
 #' @inheritParams llikNorm
@@ -108,7 +108,7 @@ llikBinom <- function(x, size, prob, full=FALSE) {
 #'
 llikNbinom <- function(x, size, prob, full=FALSE) {
   checkmate::assertIntegerish(x, min.len=0, lower=0, any.missing=FALSE)
-  checkmate::assertIntegerish(size, min.len=0, lower=0, any.missing=FALSE)
+  checkmate::assertNumeric(size, min.len=0, lower=0, any.missing=FALSE, finite=TRUE)
   checkmate::assertNumeric(prob, min.len=0, lower=0, upper=1, any.missing=FALSE, finite=TRUE)
   .df <- try(data.frame(x=x, size=size, prob=prob), silent=TRUE)
   if (inherits(.df, "try-error")) {
@@ -123,8 +123,8 @@ llikNbinom <- function(x, size, prob, full=FALSE) {
 #' 
 #' @param x  Number of successes
 #' 
-#' @param size Size of trial
-#' 
+#' @param size Dispersion parameter; may be non-integer
+#'
 #' @param mu mu parameter for negative binomial
 #' 
 #' @inheritParams llikNorm
@@ -144,7 +144,7 @@ llikNbinom <- function(x, size, prob, full=FALSE) {
 #' 
 llikNbinomMu <- function(x, size, mu, full=FALSE) {
   checkmate::assertIntegerish(x, min.len=0, lower=0, any.missing=FALSE)
-  checkmate::assertIntegerish(size, min.len=0, lower=0, any.missing=FALSE)
+  checkmate::assertNumeric(size, min.len=0, lower=0, any.missing=FALSE, finite=TRUE)
   checkmate::assertNumeric(mu, min.len=0, lower=0, any.missing=FALSE, finite=TRUE)
   .df <- try(data.frame(x=x, size=size, mu=mu), silent=TRUE)
   if (inherits(.df, "try-error")) {

src/llikNbinom.cpp

@@ -9,25 +9,25 @@
 // R , sigma=scale
 struct nbinom_llik {
   const Eigen::VectorXi y_;
-  const Eigen::VectorXi N_;
-  
-  nbinom_llik(const Eigen::VectorXi& y, Eigen::VectorXi& N) : y_(y), N_(N) { }
+  const Eigen::VectorXd N_; // dispersion; real-valued, not a count
+
+  nbinom_llik(const Eigen::VectorXi& y, Eigen::VectorXd& N) : y_(y), N_(N) { }
 
   template <typename T>
   Eigen::Matrix<T, -1, 1> operator()(const Eigen::Matrix<T, -1, 1>& theta) const {
     T p = theta[0]; // prob = size/(size+mu); size/prob-size = mu
 
     Eigen::Matrix<T, -1, 1> lp(y_.size());
     for (Eigen::Index i = 0; i < y_.size(); ++i) {
-      T mu = (double)(N_[i])*(1.0-p)/p;
+      T mu = N_[i]*(1.0-p)/p;
       lp[i] = stan::math::neg_binomial_2_lpmf(y_[i], mu, N_[i]);
     }
     return lp;
   }
 };
 
 
-stanLl llik_nbinom(Eigen::VectorXi& y, Eigen::VectorXi& N, Eigen::VectorXd& params) {
+stanLl llik_nbinom(Eigen::VectorXi& y, Eigen::VectorXd& N, Eigen::VectorXd& params) {
   rx_stan_math_thread_init_rev_autodiff();
   nbinom_llik f(y, N);
   Eigen::VectorXd fx;
@@ -69,8 +69,41 @@ static inline void llikNbinomFull(double* ret, double x, double size, double pro
     ret[5] = NA_REAL;
     return;
   }
-  if (x < 0.0 || x > static_cast<double>(INT_MAX) ||
-      size < 0.0 || size > static_cast<double>(INT_MAX)) {
+  if (x < 0.0 || x > static_cast<double>(INT_MAX) || size <= 0.0) {
+    ret[0] = isNbinom;
+    ret[1] = x;
+    ret[2] = size;
+    ret[3] = prob;
+    ret[4] = NA_REAL;
+    ret[5] = NA_REAL;
+    return;
+  }
+  if (prob == 1.0) {
+    // prob == 1 makes the distribution degenerate at zero, so the likelihood is
+    // 1 at x == 0 and 0 elsewhere, matching stats::dnbinom().  Stan cannot
+    // evaluate this point at all (it needs a strictly positive mean) and there
+    // is no usable derivative here, so report the value R would and leave the
+    // gradient NA rather than returning NA for both.  Test prob rather than the
+    // derived mean below: that mean can also underflow to 0 for a denormal size
+    // with a prob well under 1, where the true log-likelihood is finite and
+    // large-negative rather than -Inf.
+    ret[0] = isNbinom;
+    ret[1] = x;
+    ret[2] = size;
+    ret[3] = prob;
+    ret[4] = (x == 0.0) ? 0.0 : R_NegInf;
+    ret[5] = NA_REAL;
+    return;
+  }
+  // neg_binomial_2_lpmf() needs a positive finite mean, which this
+  // parameterisation reaches through mu = size*(1-prob)/prob: that is Inf at
+  // prob == 0, negative outside [0, 1], and -- now that size is no longer
+  // bounded by INT_MAX -- can overflow for a large size with a small prob, or
+  // underflow to 0 for a denormal size.  Stan throws on all of those, and the
+  // exception would escape the extern "C" entry points below and abort the R
+  // process, so return NA here.
+  double meanPar = size*(1.0-prob)/prob;
+  if (!R_finite(meanPar) || meanPar <= 0.0) {
     ret[0] = isNbinom;
     ret[1] = x;
     ret[2] = size;
@@ -80,10 +113,10 @@ static inline void llikNbinomFull(double* ret, double x, double size, double pro
     return;
   }
   Eigen::VectorXi y(1);
-  Eigen::VectorXi N(1);
+  Eigen::VectorXd N(1);
   Eigen::VectorXd params(1);
   y(0) = (int)(x);
-  N(0) = (int)(size);
+  N(0) = size;
   params(0) = prob;
   stanLl ll = llik_nbinom(y, N, params);
   ret[0] = isNbinom;

src/llikNbinom2.cpp

@@ -10,8 +10,8 @@
 
 struct nbinomMu_llik {
   const Eigen::VectorXi y_;
-  const Eigen::VectorXi N_;
-  nbinomMu_llik(const Eigen::VectorXi& y, Eigen::VectorXi& N) : y_(y), N_(N) { }
+  const Eigen::VectorXd N_; // dispersion; real-valued, not a count
+  nbinomMu_llik(const Eigen::VectorXi& y, Eigen::VectorXd& N) : y_(y), N_(N) { }
 
   template <typename T>
   Eigen::Matrix<T, -1, 1> operator()(const Eigen::Matrix<T, -1, 1>& theta) const {
@@ -24,7 +24,7 @@ struct nbinomMu_llik {
   }
 };
 
-stanLl llik_nbinomMu(Eigen::VectorXi& y, Eigen::VectorXi& N, Eigen::VectorXd& params) {
+stanLl llik_nbinomMu(Eigen::VectorXi& y, Eigen::VectorXd& N, Eigen::VectorXd& params) {
   rx_stan_math_thread_init_rev_autodiff();
   nbinomMu_llik f(y, N);
   Eigen::VectorXd fx;
@@ -67,8 +67,13 @@ static inline void llikNbinomMuFull(double* ret, double x, double size, double m
     ret[5] = NA_REAL;
     return;
   }
+  // neg_binomial_2_lpmf() needs a positive finite mean; a negative mu (reachable
+  // through the C API, which bypasses the R-level assertions) makes it throw,
+  // and that exception would escape the extern "C" entry points below and abort
+  // the R process, so return NA here instead.  mu == 0 is degenerate rather than
+  // out of domain and is handled just below.
   if (x < 0.0 || x > static_cast<double>(INT_MAX) ||
-      size < 0.0 || size > static_cast<double>(INT_MAX)) {
+      size <= 0.0 || mu < 0.0) {
     ret[0] = isNbinomMu;
     ret[1] = x;
     ret[2] = size;
@@ -77,11 +82,25 @@ static inline void llikNbinomMuFull(double* ret, double x, double size, double m
     ret[5] = NA_REAL;
     return;
   }
+  if (mu == 0.0) {
+    // mu == 0 makes the distribution degenerate at zero, so the likelihood is 1
+    // at x == 0 and 0 elsewhere, matching stats::dnbinom().  Stan cannot
+    // evaluate this point at all (it needs a strictly positive mean) and there
+    // is no usable derivative here, so report the value R would and leave the
+    // gradient NA rather than returning NA for both.
+    ret[0] = isNbinomMu;
+    ret[1] = x;
+    ret[2] = size;
+    ret[3] = mu;
+    ret[4] = (x == 0.0) ? 0.0 : R_NegInf;
+    ret[5] = NA_REAL;
+    return;
+  }
   Eigen::VectorXi y(1);
-  Eigen::VectorXi N(1);
+  Eigen::VectorXd N(1);
   Eigen::VectorXd params(1);
   y(0) = (int)(x);
-  N(0) = (int)(size);
+  N(0) = size;
   params(0) = mu;
   stanLl ll = llik_nbinomMu(y, N, params);
   ret[0] = isNbinomMu;

tests/testthat/test-llik.R

@@ -39,6 +39,31 @@ test_that("log-liklihood tests for NbinomMu (including derivatives)", {
   expect_equal(fromR$fx, dnbinom(et$time, size=100, mu=40, log=TRUE))
 })
 
+test_that("nbinom size may be continuous (non-integer dispersion)", {
+  # size used to be truncated to an int, silently returning the value at
+  # trunc(size) for size > 1 and aborting the process for 0 < size < 1.
+  size <- c(0.05, 0.3, 0.7, 0.99, 1.5, 2.5, 3.9)
+  x <- rep(2L, length(size))
+  expect_equal(llikNbinomMu(x, size, rep(5, length(size)))$fx,
+               dnbinom(2L, size=size, mu=5, log=TRUE))
+  expect_equal(llikNbinom(x, size, rep(0.4, length(size)))$fx,
+               dnbinom(2L, size=size, prob=0.4, log=TRUE))
+})
+
+test_that("nbinom derivatives are correct for continuous size", {
+  h <- 1e-5
+  size <- c(0.3, 0.7, 1.5, 4.25)
+  x <- rep(2L, length(size))
+  dMuNum <- (dnbinom(2L, size=size, mu=5 + h, log=TRUE) -
+             dnbinom(2L, size=size, mu=5 - h, log=TRUE)) / (2 * h)
+  expect_equal(llikNbinomMu(x, size, rep(5, length(size)))$dMu, dMuNum,
+               tolerance=1e-5)
+  dProbNum <- (dnbinom(2L, size=size, prob=0.4 + h, log=TRUE) -
+               dnbinom(2L, size=size, prob=0.4 - h, log=TRUE)) / (2 * h)
+  expect_equal(llikNbinom(x, size, rep(0.4, length(size)))$dProb, dProbNum,
+               tolerance=1e-5)
+})
+
 test_that("log-liklihood tests for beta (including derivatives)", {
   et <- data.frame(time=seq(1e-4, 1-1e-4, length.out=21))
   et$shape1 <- 0.5
@@ -293,11 +318,55 @@ test_that("llikNbinomInternal returns NA for x > INT_MAX", {
   expect_true(is.na(res$dProb))
 })
 
-test_that("llikNbinomInternal returns NA for size > INT_MAX", {
+test_that("llikNbinomInternal accepts size > INT_MAX (size is a real dispersion)", {
   big <- 2^31
   res <- llikNbinomInternal(5, big, 0.5)
-  expect_true(is.na(res$fx))
-  expect_true(is.na(res$dProb))
+  expect_equal(res$fx, stats::dnbinom(5, size = big, prob = 0.5, log = TRUE))
+})
+
+test_that("llikNbinomInternal returns NA for size <= 0", {
+  expect_true(is.na(llikNbinomInternal(5, 0, 0.5)$fx))
+  expect_true(is.na(llikNbinomInternal(5, -1, 0.5)$fx))
+})
+
+test_that("llikNbinomInternal returns NA for an out-of-domain prob", {
+  # mu = size*(1-prob)/prob is Inf at prob == 0 and negative outside [0, 1];
+  # each makes neg_binomial_2_lpmf() throw, which aborts the process when it
+  # escapes rxLlikNbinom()
+  expect_true(is.na(llikNbinomInternal(5, 10, 0)$fx))
+  expect_true(is.na(llikNbinomInternal(5, 10, -0.5)$fx))
+  expect_true(is.na(llikNbinomInternal(5, 10, 1.5)$fx))
+})
+
+test_that("llikNbinom is degenerate at zero for prob == 1", {
+  # prob == 1 puts all the mass at zero, as stats::dnbinom() does; Stan cannot
+  # evaluate the point (mu == 0), so the value is filled in and the gradient,
+  # which does not exist here, is NA
+  res <- llikNbinomInternal(c(0, 1, 5), rep(10, 3), rep(1, 3))
+  expect_equal(res$fx, stats::dnbinom(c(0, 1, 5), size = 10, prob = 1, log = TRUE))
+  expect_equal(res$fx, c(0, -Inf, -Inf))
+  expect_true(all(is.na(res$dProb)))
+  # a continuous size is degenerate at zero too
+  expect_equal(llikNbinomInternal(c(0, 1), c(0.5, 0.5), c(1, 1))$fx,
+               stats::dnbinom(c(0, 1), size = 0.5, prob = 1, log = TRUE))
+  # and it is reachable from the R-level interface
+  expect_equal(llikNbinom(0L, 10, 1)$fx, 0)
+})
+
+test_that("llikNbinomInternal returns NA when size/prob overflows the mean", {
+  # size is no longer bounded by INT_MAX, so size*(1-prob)/prob can overflow
+  # even though both size and prob are finite and in range
+  expect_true(is.na(llikNbinomInternal(5, 1e300, 1e-10)$fx))
+})
+
+test_that("llikNbinomInternal returns NA when a denormal size underflows the mean", {
+  # size*(1-prob)/prob underflows to exactly 0 for a denormal size even with a
+  # prob well under 1.  That is not the degenerate prob == 1 point -- the true
+  # log-likelihood is finite (about -745 here) -- so it must return NA rather
+  # than the -Inf used at the point mass
+  expect_true(is.finite(stats::dnbinom(1, size = 5e-324, prob = 0.5, log = TRUE)))
+  expect_true(is.na(llikNbinomInternal(1, 5e-324, 0.5)$fx))
+  expect_true(is.na(llikNbinomInternal(0, 5e-324, 0.5)$fx))
 })
 
 test_that("llikNbinomMuInternal returns NA for x > INT_MAX", {
@@ -307,11 +376,36 @@ test_that("llikNbinomMuInternal returns NA for x > INT_MAX", {
   expect_true(is.na(res$dMu))
 })
 
-test_that("llikNbinomMuInternal returns NA for size > INT_MAX", {
+test_that("llikNbinomMuInternal accepts size > INT_MAX (size is a real dispersion)", {
   big <- 2^31
   res <- llikNbinomMuInternal(5, big, 40)
-  expect_true(is.na(res$fx))
-  expect_true(is.na(res$dMu))
+  expect_equal(res$fx, stats::dnbinom(5, size = big, mu = 40, log = TRUE))
+})
+
+test_that("llikNbinomMuInternal returns NA for size <= 0", {
+  expect_true(is.na(llikNbinomMuInternal(5, 0, 40)$fx))
+  expect_true(is.na(llikNbinomMuInternal(5, -1, 40)$fx))
+})
+
+test_that("llikNbinomMuInternal returns NA for mu < 0", {
+  # neg_binomial_2_lpmf() throws on a negative mean, which aborts the process
+  # when it escapes rxLlikNbinomMu()
+  expect_true(is.na(llikNbinomMuInternal(5, 10, -1)$fx))
+})
+
+test_that("llikNbinomMu is degenerate at zero for mu == 0", {
+  # mu == 0 puts all the mass at zero, as stats::dnbinom() does; Stan cannot
+  # evaluate the point, so the value is filled in and the gradient, which does
+  # not exist here, is NA
+  res <- llikNbinomMuInternal(c(0, 1, 5), rep(10, 3), rep(0, 3))
+  expect_equal(res$fx, stats::dnbinom(c(0, 1, 5), size = 10, mu = 0, log = TRUE))
+  expect_equal(res$fx, c(0, -Inf, -Inf))
+  expect_true(all(is.na(res$dMu)))
+  # a continuous size is degenerate at zero too
+  expect_equal(llikNbinomMuInternal(c(0, 1), c(0.5, 0.5), c(0, 0))$fx,
+               stats::dnbinom(c(0, 1), size = 0.5, mu = 0, log = TRUE))
+  # and it is reachable from the R-level interface
+  expect_equal(llikNbinomMu(0L, 10, 0)$fx, 0)
 })
 
 ## Large-vector test for R_xlen_t fix