Merge pull request #591 from nlmixr2/590-confint-for-binary-probability

590 confint for binary probability

NAMESPACE

@@ -41,6 +41,7 @@ S3method(as.rxUi,rxModelVars)
 S3method(as.rxUi,rxUi)
 S3method(as.rxUi,rxode2)
 S3method(as.rxUi,rxode2tos)
+S3method(binomProbs,default)
 S3method(coef,rxode2)
 S3method(confint,rxSolve)
 S3method(dimnames,rxSolve)
@@ -246,6 +247,7 @@ export(assertRxUiPrediction)
 export(assertRxUiRandomOnIdOnly)
 export(assertRxUiSingleEndpoint)
 export(assertRxUiTransformNormal)
+export(binomProbs)
 export(cvPost)
 export(erf)
 export(et)

NEWS.md

@@ -90,13 +90,22 @@ mu-referencing style to run the optimization.
   your own `ci=0.99` for instance
 
 - A new function was introduced `meanProbs()` which calculates the
-  mean and expected quantiles under either the normal or t
+  mean and expected confidence bands under either the normal or t
   distribution
 
+- A related new function was introduced that calculates the mean and
+  confidence bands under the binomial/Bernoulli distribution
+  (`binomProbs()`)
+
 - When calculating the intervals for `rxode2` simulated objects you
   can also use `mean=TRUE` to use the mean for the first level of
   confidence using `meanProbs()`
 
+- Also when calculating the intervals for `rxode2` simulated object
+  you can also use `mean="prob"` to use the binomial distributional
+  information (and ci) for the first level of confidence using
+  `binomProbs()`
+
 - When plotting the `confint` derived intervals from an `rxode2`
   simulation, you can now subset based on a simulated value like
   `plot(ci, Cc)` which will only plot the variable `Cc` that you

R/RcppExports.R

@@ -589,7 +589,11 @@ rxErf <- function(v) {
     .Call(`_rxode2_rxErf`, v)
 }
 
-meanProbs_ <- function(x, probs, naRm, useT) {
-    .Call(`_rxode2_meanProbs_`, x, probs, naRm, useT)
+binomProbs_ <- function(x, probs, naRm) {
+    .Call(`_rxode2_binomProbs_`, x, probs, naRm)
+}
+
+meanProbs_ <- function(x, probs, naRm, useT, useBinom) {
+    .Call(`_rxode2_meanProbs_`, x, probs, naRm, useT, useBinom)
 }
 

R/confint.R

@@ -32,9 +32,16 @@ confint.rxSolve <- function(object, parm = NULL, level = 0.95, ...) {
     }
   }
   .mean <- FALSE
+  .prob <- FALSE
   if (any(names(.args) == "mean")) {
     .mean <- .args$mean
-    checkmate::assertLogical(.mean, len=1, any.missing=FALSE, .var.name="mean")
+    if (inherits(.mean, "character") &&
+          any(.mean == c("binom", "prob"))) {
+      .prob <- TRUE
+      .mean <- FALSE
+    } else {
+      checkmate::assertLogical(.mean, len=1, any.missing=FALSE, .var.name="mean")
+    }
   }
   .stk <- rxStack(object, parm, doSim=.doSim)
   for(.v in .by) {
@@ -50,7 +57,8 @@ confint.rxSolve <- function(object, parm = NULL, level = 0.95, ...) {
     ci = paste0("p", .p2 * 100),
     parm = levels(.stk$trt),
     by = .by,
-    mean = .mean
+    mean = .mean,
+    prob=.prob
   )
   class(.lst) <- "rxHidden"
   if (.ci ==0 || !any(names(.stk) == "sim.id")) {
@@ -75,7 +83,14 @@ confint.rxSolve <- function(object, parm = NULL, level = 0.95, ...) {
           Percentile = sprintf("%s%%", .p * 100)
         ),
         by = c("time", "trt", .by)
-        ]        
+        ]
+      } else if (.prob) {
+        .stk <- .stk[, list(
+          p1 = .p, eff = rxode2::binomProbs(.SD$value, probs = .p, na.rm = TRUE),
+          Percentile = sprintf("%s%%", .p * 100)
+        ),
+        by = c("time", "trt", .by)
+        ]
       } else {
         .stk <- .stk[, list(
           p1 = .p, eff = stats::quantile(.SD$value, probs = .p, na.rm = TRUE),
@@ -107,10 +122,14 @@ confint.rxSolve <- function(object, parm = NULL, level = 0.95, ...) {
     .ret <- .ret[, list(p1 = .p,
                         eff = rxode2::meanProbs(.SD$value, probs = .p, na.rm = TRUE)),
                  by = c("id", "time", "trt", .by)]
+  } else if (.prob) {
+    .ret <- .ret[, list(p1 = .p,
+                        eff = rxode2::binomProbs(.SD$value, probs = .p, na.rm = TRUE)),
+                 by = c("id", "time", "trt", .by)]
   } else {
     .ret <- .ret[, list(p1 = .p,
                         eff = stats::quantile(.SD$value, probs = .p, na.rm = TRUE)), by = c("id", "time", "trt", .by)]
-    
+
   }
   .ret <- .ret[, setNames(as.list(stats::quantile(.SD$eff, probs = .p2, na.rm = TRUE)),
                           sprintf("p%s", .p2 * 100)),

R/utils.R

@@ -747,15 +747,15 @@ is.latex <- function() {
   ret
 }
 #' Print out a table in the documentation
-#'  
+#'
 #' @param table data frame
 #' @param caption a character vector representing the caption for the latex table
 #' @return based on the `knitr` context:
 #' - output a `kableExtra::kbl` for `latex` output
 #' - output a `DT::datatable` for html output
 #' - otherwise output a `knitr::kable`
 #' @keywords internal
-#' @export 
+#' @export
 #' @author Matthew L. Fidler
 #' @examples
 #' .rxDocTable(rxReservedKeywords)
@@ -773,11 +773,12 @@ is.latex <- function() {
   }
 }
 
-#' Calculate expected quantiles with normal or t sampling distribution
+#' Calculate expected confidence bands with normal or t sampling distribution
 #'
-#' The generic function `meanProbs` produces expected quantiles under
-#' either the t distribution or the normal sampling distribution. This
-#' uses `qnorm()` or `qt()` with the mean and standard deviation.
+#' The generic function `meanProbs` produces expected confidence bands
+#' under either the t distribution or the normal sampling
+#' distribution. This uses `qnorm()` or `qt()` with the mean and
+#' standard deviation.
 #'
 #' For a single probability, p, it uses either:
 #'
@@ -796,11 +797,11 @@ is.latex <- function() {
 #' This is meant to perform in the same way as `quantile()` so it can
 #' be a drop in replacement for code using `quantile()` but using
 #' distributional assumptions.
-#'  
+#'
 #' @param x numeric vector whose mean and probability based confidence
 #'   values are wanted, NA and NaN values are not allowed in numeric
 #'   vectors unless ‘na.rm’ is ‘TRUE’.
-#' @param probs numeric vector of probabilities with values in [0,1].
+#' @param probs numeric vector of probabilities with values in the interval from 0 to 1 .
 #' @param na.rm logical; if true, any NA and NaN's are removed from
 #'   `x` before the quantiles are computed.
 #' @param names logical; if true, the result has a names attribute.
@@ -809,6 +810,8 @@ is.latex <- function() {
 #'   distribution to calculate the confidence based estimates.
 #' @param onlyProbs logical; if true, only return the probability based
 #'   confidence interval estimates, otherwise return
+#' @param ... Arguments passed to default method, allows many
+#'   different methods to be applied.
 #' @return By default the return has the probabilities as names (if
 #'   named) with the points where the expected distribution are
 #'   located given the sampling mean and standard deviation. If
@@ -831,9 +834,9 @@ is.latex <- function() {
 #' quantile(x<- rnorm(42))
 #'
 #' meanProbs(x)
-#' 
+#'
 #' meanProbs(x, useT=FALSE)
-#' 
+#'
 meanProbs <- function(x, ...) {
   UseMethod("meanProbs")
 }
@@ -864,3 +867,84 @@ meanProbs.default <- function(x, probs=seq(0, 1, 0.25), na.rm=FALSE,
   .ret
 }
 
+#' Calculate expected confidence bands with binomial sampling distribution
+#'
+#' The generic function `binomProbs` produces expected confidence bands
+#' using the
+#'
+#' This is meant to perform in the same way as `quantile()` so it can
+#' be a drop in replacement for code using `quantile()` but using
+#' distributional assumptions.
+#'
+#' It is used for confidence intervals with rxode2 solved objects using
+#' `confint(mean="prob")` or `confint(mean="binom")
+#'
+#' @param x numeric vector whose mean and probability based confidence
+#'   values are wanted, NA and NaN values are not allowed in numeric
+#'   vectors unless ‘na.rm’ is ‘TRUE’.
+#' @param probs numeric vector of probabilities with values in the
+#'   interval 0 to 1, inclusive. When 0, it represents the maximum
+#'   observed, when 1, it represents the maximum observed. When 0.5 it
+#'   represents the expected probability (mean)
+#' @param na.rm logical; if true, any NA and NaN's are removed from
+#'   `x` before the quantiles are computed.
+#' @param names logical; if true, the result has a names attribute.
+#' @param onlyProbs logical; if true, only return the probability
+#'   based confidence interval estimates, otherwise return
+#' @param ... Arguments passed to default method, allows many
+#'   different methods to be applied.
+#' @return By default the return has the probabilities as names (if
+#'   named) with the points where the expected distribution are
+#'   located given the sampling mean and standard deviation. If
+#'   `onlyProbs=FALSE` then it would prepend mean, variance, standard
+#'   deviation, minimum, maximum and number of non-NA observations.
+#' @export
+#' @author Matthew L. Fidler
+#' @references - Newcombe,
+#'   R. G. (1998). "Two-sided confidence intervals for the single proportion: comparison of seven methods". Statistics
+#'   in Medicine. 17 (8):
+#'   857–872. doi:10.1002/(SICI)1097-0258(19980430)17:8<857::AID-SIM777>3.0.CO;2-E. PMID
+#'   9595616.
+#' @examples
+#'
+#' quantile(x<- rbinom(7001, p=0.375, size=1))
+#' binomProbs(x)
+#'
+#' # Can get some extra statistics if you request onlyProbs=FALSE
+#' binomProbs(x, onlyProbs=FALSE)
+#'
+#' x[2] <- NA_real_
+#'
+#' binomProbs(x, onlyProbs=FALSE)
+#'
+#' binomProbs(x, na.rm=TRUE)
+#'
+binomProbs <- function(x, ...) {
+  UseMethod("binomProbs")
+}
+
+#' @rdname binomProbs
+#' @export
+binomProbs.default <- function(x, probs=c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm=FALSE,
+                               names=TRUE, onlyProbs=TRUE) {
+  checkmate::assertIntegerish(x, min.len=1, lower=0.0, upper=1.0)
+  x <- as.double(x)
+  checkmate::assertNumeric(probs, min.len=1, any.missing = FALSE, lower=0.0, upper=1.0)
+  checkmate::assertLogical(na.rm, any.missing=FALSE, len=1)
+  checkmate::assertLogical(names, any.missing=FALSE, len=1)
+  checkmate::assertLogical(onlyProbs, any.missing=FALSE, len=1)
+  .ret <- .Call(`_rxode2_binomProbs_`, x, probs, na.rm)
+  .names <- NULL
+  if (names) {
+    .names <- paste0(probs*100, "%")
+  }
+  if (onlyProbs) {
+    .ret <- .ret[-1L:-4L]
+    if (names) {
+      names(.ret) <- .names
+    }
+  } else if (names) {
+    names(.ret) <- c("mean","var", "sd", "n", .names)
+  }
+  .ret
+}

_pkgdown.yml

@@ -65,7 +65,9 @@ reference:
   - rxPkg
 - title: Specialized Simulation functions
   contents:
+  - binomProbs
   - cvPost
+  - meanProbs
   - rinvchisq
   - rxGetSeed
   - rxPp

inst/include/rxode2_RcppExports.h

@@ -1073,17 +1073,38 @@ namespace rxode2 {
         return Rcpp::as<NumericVector >(rcpp_result_gen);
     }
 
-    inline NumericVector meanProbs_(NumericVector x, NumericVector probs, bool naRm, bool useT) {
-        typedef SEXP(*Ptr_meanProbs_)(SEXP,SEXP,SEXP,SEXP);
+    inline NumericVector binomProbs_(NumericVector x, NumericVector probs, bool naRm) {
+        typedef SEXP(*Ptr_binomProbs_)(SEXP,SEXP,SEXP);
+        static Ptr_binomProbs_ p_binomProbs_ = NULL;
+        if (p_binomProbs_ == NULL) {
+            validateSignature("NumericVector(*binomProbs_)(NumericVector,NumericVector,bool)");
+            p_binomProbs_ = (Ptr_binomProbs_)R_GetCCallable("rxode2", "_rxode2_binomProbs_");
+        }
+        RObject rcpp_result_gen;
+        {
+            RNGScope RCPP_rngScope_gen;
+            rcpp_result_gen = p_binomProbs_(Shield<SEXP>(Rcpp::wrap(x)), Shield<SEXP>(Rcpp::wrap(probs)), Shield<SEXP>(Rcpp::wrap(naRm)));
+        }
+        if (rcpp_result_gen.inherits("interrupted-error"))
+            throw Rcpp::internal::InterruptedException();
+        if (Rcpp::internal::isLongjumpSentinel(rcpp_result_gen))
+            throw Rcpp::LongjumpException(rcpp_result_gen);
+        if (rcpp_result_gen.inherits("try-error"))
+            throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
+        return Rcpp::as<NumericVector >(rcpp_result_gen);
+    }
+
+    inline NumericVector meanProbs_(NumericVector x, NumericVector probs, bool naRm, bool useT, bool useBinom) {
+        typedef SEXP(*Ptr_meanProbs_)(SEXP,SEXP,SEXP,SEXP,SEXP);
         static Ptr_meanProbs_ p_meanProbs_ = NULL;
         if (p_meanProbs_ == NULL) {
-            validateSignature("NumericVector(*meanProbs_)(NumericVector,NumericVector,bool,bool)");
+            validateSignature("NumericVector(*meanProbs_)(NumericVector,NumericVector,bool,bool,bool)");
             p_meanProbs_ = (Ptr_meanProbs_)R_GetCCallable("rxode2", "_rxode2_meanProbs_");
         }
         RObject rcpp_result_gen;
         {
             RNGScope RCPP_rngScope_gen;
-            rcpp_result_gen = p_meanProbs_(Shield<SEXP>(Rcpp::wrap(x)), Shield<SEXP>(Rcpp::wrap(probs)), Shield<SEXP>(Rcpp::wrap(naRm)), Shield<SEXP>(Rcpp::wrap(useT)));
+            rcpp_result_gen = p_meanProbs_(Shield<SEXP>(Rcpp::wrap(x)), Shield<SEXP>(Rcpp::wrap(probs)), Shield<SEXP>(Rcpp::wrap(naRm)), Shield<SEXP>(Rcpp::wrap(useT)), Shield<SEXP>(Rcpp::wrap(useBinom)));
         }
         if (rcpp_result_gen.inherits("interrupted-error"))
             throw Rcpp::internal::InterruptedException();