Merge pull request #578 from nlmixr2/577-allow-confint-to-apply-to-me…

…ans-as-well

577 allow confint to apply to means as well

NAMESPACE

@@ -52,6 +52,7 @@ S3method(getBaseSimModel,default)
 S3method(getBaseSymengineModel,default)
 S3method(ini,default)
 S3method(ini,rxUi)
+S3method(meanProbs,default)
 S3method(model,"function")
 S3method(model,default)
 S3method(model,rxModelVars)
@@ -302,6 +303,7 @@ export(logit)
 export(logitNormInfo)
 export(lotri)
 export(lowergamma)
+export(meanProbs)
 export(model)
 export(modelExtract)
 export(odeMethodToInt)

NEWS.md

@@ -83,11 +83,19 @@ mu-referencing style to run the optimization.
   example you can now stratify by gender and race by: `confint(sim,
   "sim", by=c("race", "gender"))`
 
-  - When calculating the intervals for `rxode2` simulated objects you
+- When calculating the intervals for `rxode2` simulated objects you
   can now use `ci=FALSE` so that it only calculates the default
   intervals without bands on each of the percentiles; You can also
   choose not to match the secondary bands limits with `levels` but use
   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
+  distribution
+
+- 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()`
 
 ## Internal new features
 

R/RcppExports.R

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

R/confint.R

@@ -23,14 +23,19 @@ confint.rxSolve <- function(object, parm = NULL, level = 0.95, ...) {
   if (any(names(.args) == "ci")) {
     .ci <- .args$ci
     if (inherits(.ci, "logical")) {
-      checkmate::assertLogical(.ci, len=1, any.missing=FALSE)
+      checkmate::assertLogical(.ci, len=1, any.missing=FALSE, .var.name="ci")
       if (!.ci) {
         .ci <- 0.0
       }
     } else {
-      checkmate::assertNumeric(.ci, lower=0, upper=1, finite=TRUE, any.missing=FALSE)
+      checkmate::assertNumeric(.ci, lower=0, upper=1, finite=TRUE, any.missing=FALSE, .var.name="ci")
     }
   }
+  .mean <- FALSE
+  if (any(names(.args) == "mean")) {
+    .mean <- .args$mean
+    checkmate::assertLogical(.mean, len=1, any.missing=FALSE, .var.name="mean")
+  }
   .stk <- rxStack(object, parm, doSim=.doSim)
   for(.v in .by) {
     .stk[[.v]] <- object[[.v]]
@@ -44,7 +49,8 @@ confint.rxSolve <- function(object, parm = NULL, level = 0.95, ...) {
     lvl = paste0("p", .p * 100),
     ci = paste0("p", .p2 * 100),
     parm = levels(.stk$trt),
-    by = .by
+    by = .by,
+    mean = .mean
   )
   class(.lst) <- "rxHidden"
   if (.ci ==0 || !any(names(.stk) == "sim.id")) {
@@ -61,14 +67,23 @@ confint.rxSolve <- function(object, parm = NULL, level = 0.95, ...) {
     if (.ci == 0 || .ntot < 2500) {
       if (.ci != 0.0) {
         .mwarn("in order to put confidence bands around the intervals, you need at least 2500 simulations")
-        message("summarizing data...", appendLF = FALSE)
       }
-      .stk <- .stk[, list(
-        p1 = .p, eff = stats::quantile(.SD$value, probs = .p, na.rm = TRUE),
-        Percentile = sprintf("%s%%", .p * 100)
-      ),
-      by = c("time", "trt", .by)
-      ]
+      message("summarizing data...", appendLF = FALSE)
+      if (.mean) {
+        .stk <- .stk[, list(
+          p1 = .p, eff = rxode2::meanProbs(.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),
+          Percentile = sprintf("%s%%", .p * 100)
+        ),
+        by = c("time", "trt", .by)
+        ]
+      }
       if (requireNamespace("tibble", quietly = TRUE)) {
         .stk <- tibble::as_tibble(.stk)
       }
@@ -88,7 +103,15 @@ confint.rxSolve <- function(object, parm = NULL, level = 0.95, ...) {
   }
   message("summarizing data...", appendLF = FALSE)
   .ret <- .stk[, id := sim.id %% .n]
-  .ret <- .ret[, list(p1 = .p, eff = stats::quantile(.SD$value, probs = .p, na.rm = TRUE)), by = c("id", "time", "trt", .by)]
+  if (.mean) {
+    .ret <- .ret[, list(p1 = .p,
+                        eff = rxode2::meanProbs(.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)),
                by = c("p1", "time", "trt", .by)

R/rxsolve.R

@@ -1,4 +1,4 @@
-#' Options, Solving & Simulation of an ODE/solved system
+ #' Options, Solving & Simulation of an ODE/solved system
 #'
 #' This uses rxode2 family of objects, file, or model specification to
 #' solve a ODE system.  There are many options for a solved rxode2

R/utils.R

@@ -772,3 +772,95 @@ is.latex <- function() {
     knitr::kable(table)
   }
 }
+
+#' Calculate expected quantiles 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.
+#'
+#' For a single probability, p, it uses either:
+#'
+#' mean + qt(p, df=n)*sd/sqrt(n)
+#'
+#' or
+#'
+#' mean + qnorm(p)*sd/sqrt(n)
+#'
+#' The smallest observation corresponds to a probability of 0 and the
+#' largest to a probability of 1 and the mean corresponds to 0.5.
+#'
+#' The mean and standard deviation of the sample is calculated based
+#' on Welford's method for a single pass.
+#'
+#' 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 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 useT logical; if true, use the t-distribution to calculate
+#'   the confidence-based estimates. If false use the normal
+#'   distribution to calculate the confidence based estimates.
+#' @param onlyProbs logical; if true, only return the probability based
+#'   confidence interval estimates, otherwise return
+#' @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
+#' @examples
+#'
+#' quantile(x<- rnorm(1001))
+#' meanProbs(x)
+#'
+#' # Can get some extra statistics if you request onlyProbs=FALSE
+#' meanProbs(x, onlyProbs=FALSE)
+#'
+#' x[2] <- NA_real_
+#'
+#' meanProbs(x, onlyProbs=FALSE)
+#'
+#' quantile(x<- rnorm(42))
+#'
+#' meanProbs(x)
+#' 
+#' meanProbs(x, useT=FALSE)
+#' 
+meanProbs <- function(x, ...) {
+  UseMethod("meanProbs")
+}
+
+#' @rdname meanProbs
+#' @export
+meanProbs.default <- function(x, probs=seq(0, 1, 0.25), na.rm=FALSE,
+                              names=TRUE, useT=TRUE, onlyProbs=TRUE) {
+  checkmate::assertNumeric(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(useT, any.missing=FALSE, len=1)
+  checkmate::assertLogical(onlyProbs, any.missing=FALSE, len=1)
+  .ret <- .Call(`_rxode2_meanProbs_`, x, probs, na.rm, useT)
+  .names <- NULL
+  if (names) {
+    .names <- paste0(probs*100, "%")
+  }
+  if (onlyProbs) {
+    .ret <- .ret[-1L:-6L]
+    if (names) {
+      names(.ret) <- .names
+    }
+  } else if (names) {
+    names(.ret) <- c("mean","var", "sd", "min", "max", "n", .names)
+  }
+  .ret
+}
+

inst/include/rxode2_RcppExports.h

@@ -1073,6 +1073,27 @@ 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);
+        static Ptr_meanProbs_ p_meanProbs_ = NULL;
+        if (p_meanProbs_ == NULL) {
+            validateSignature("NumericVector(*meanProbs_)(NumericVector,NumericVector,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)));
+        }
+        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);
+    }
+
 }
 
 #endif // RCPP_rxode2_RCPPEXPORTS_H_GEN_

src/RcppExports.cpp

@@ -1824,6 +1824,43 @@ RcppExport SEXP _rxode2_rxErf(SEXP vSEXP) {
     UNPROTECT(1);
     return rcpp_result_gen;
 }
+// meanProbs_
+NumericVector meanProbs_(NumericVector x, NumericVector probs, bool naRm, bool useT);
+static SEXP _rxode2_meanProbs__try(SEXP xSEXP, SEXP probsSEXP, SEXP naRmSEXP, SEXP useTSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP);
+    Rcpp::traits::input_parameter< NumericVector >::type probs(probsSEXP);
+    Rcpp::traits::input_parameter< bool >::type naRm(naRmSEXP);
+    Rcpp::traits::input_parameter< bool >::type useT(useTSEXP);
+    rcpp_result_gen = Rcpp::wrap(meanProbs_(x, probs, naRm, useT));
+    return rcpp_result_gen;
+END_RCPP_RETURN_ERROR
+}
+RcppExport SEXP _rxode2_meanProbs_(SEXP xSEXP, SEXP probsSEXP, SEXP naRmSEXP, SEXP useTSEXP) {
+    SEXP rcpp_result_gen;
+    {
+        Rcpp::RNGScope rcpp_rngScope_gen;
+        rcpp_result_gen = PROTECT(_rxode2_meanProbs__try(xSEXP, probsSEXP, naRmSEXP, useTSEXP));
+    }
+    Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error");
+    if (rcpp_isInterrupt_gen) {
+        UNPROTECT(1);
+        Rf_onintr();
+    }
+    bool rcpp_isLongjump_gen = Rcpp::internal::isLongjumpSentinel(rcpp_result_gen);
+    if (rcpp_isLongjump_gen) {
+        Rcpp::internal::resumeJump(rcpp_result_gen);
+    }
+    Rboolean rcpp_isError_gen = Rf_inherits(rcpp_result_gen, "try-error");
+    if (rcpp_isError_gen) {
+        SEXP rcpp_msgSEXP_gen = Rf_asChar(rcpp_result_gen);
+        UNPROTECT(1);
+        Rf_error(CHAR(rcpp_msgSEXP_gen));
+    }
+    UNPROTECT(1);
+    return rcpp_result_gen;
+}
 
 // validate (ensure exported C++ functions exist before calling them)
 static int _rxode2_RcppExport_validate(const char* sig) { 
@@ -1879,6 +1916,7 @@ static int _rxode2_RcppExport_validate(const char* sig) {
         signatures.insert("RObject(*rxSymInvCholEnvCalculate)(List,std::string,Nullable<NumericVector>)");
         signatures.insert("LogicalVector(*isNullZero)(RObject)");
         signatures.insert("NumericVector(*rxErf)(NumericVector)");
+        signatures.insert("NumericVector(*meanProbs_)(NumericVector,NumericVector,bool,bool)");
     }
     return signatures.find(sig) != signatures.end();
 }
@@ -1935,6 +1973,7 @@ RcppExport SEXP _rxode2_RcppExport_registerCCallable() {
     R_RegisterCCallable("rxode2", "_rxode2_rxSymInvCholEnvCalculate", (DL_FUNC)_rxode2_rxSymInvCholEnvCalculate_try);
     R_RegisterCCallable("rxode2", "_rxode2_isNullZero", (DL_FUNC)_rxode2_isNullZero_try);
     R_RegisterCCallable("rxode2", "_rxode2_rxErf", (DL_FUNC)_rxode2_rxErf_try);
+    R_RegisterCCallable("rxode2", "_rxode2_meanProbs_", (DL_FUNC)_rxode2_meanProbs__try);
     R_RegisterCCallable("rxode2", "_rxode2_RcppExport_validate", (DL_FUNC)_rxode2_RcppExport_validate);
     return R_NilValue;
 }