diff --git a/NAMESPACE b/NAMESPACE index 81962649a..a56557c2d 100644 --- a/NAMESPACE +++ b/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) diff --git a/NEWS.md b/NEWS.md index b9d933b9d..02d50a85f 100644 --- a/NEWS.md +++ b/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 diff --git a/R/RcppExports.R b/R/RcppExports.R index e776c74b3..30cc3e92f 100644 --- a/R/RcppExports.R +++ b/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) +} + diff --git a/R/confint.R b/R/confint.R index 9fb7ee507..8094f9d58 100644 --- a/R/confint.R +++ b/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) diff --git a/R/rxsolve.R b/R/rxsolve.R index 670ef939b..34a52723e 100644 --- a/R/rxsolve.R +++ b/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 diff --git a/R/utils.R b/R/utils.R index df510653b..b7694a96c 100644 --- a/R/utils.R +++ b/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 +} + diff --git a/inst/include/rxode2_RcppExports.h b/inst/include/rxode2_RcppExports.h index c50a082a7..cc821793b 100644 --- a/inst/include/rxode2_RcppExports.h +++ b/inst/include/rxode2_RcppExports.h @@ -1073,6 +1073,27 @@ namespace rxode2 { return Rcpp::as(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(Rcpp::wrap(x)), Shield(Rcpp::wrap(probs)), Shield(Rcpp::wrap(naRm)), Shield(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(rcpp_result_gen).c_str()); + return Rcpp::as(rcpp_result_gen); + } + } #endif // RCPP_rxode2_RCPPEXPORTS_H_GEN_ diff --git a/man/meanProbs.Rd b/man/meanProbs.Rd new file mode 100644 index 000000000..0dcec476a --- /dev/null +++ b/man/meanProbs.Rd @@ -0,0 +1,65 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils.R +\name{meanProbs} +\alias{meanProbs} +\alias{meanProbs.default} +\title{Calculate expected quantiles with normal or t sampling distribution} +\usage{ +meanProbs(x, ...) + +\method{meanProbs}{default}( + x, + probs = seq(0, 1, 0.25), + na.rm = FALSE, + names = TRUE, + useT = TRUE, + useProbs = FALSE +) +} +\arguments{ +\item{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’.} + +\item{probs}{numeric vector of probabilities with values in \link{0,1}.} + +\item{na.rm}{logical; if true, any NA and NaN's are removed from +\code{x} before the quantiles are computed.} + +\item{names}{logical; if true, the result has a names attribute.} + +\item{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.} + +\item{useProbs}{logical; if true, only return the probability +based confidence interval estimates, otherwise return} +} +\description{ +The generic function \code{meanProbs} produces expected quantiles under +either the t distribution or the normal sampling distribution. This +uses \code{qnorm()} or \code{qt()} with the mean and standard deviation. +} +\details{ +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 +} +\examples{ + +meanProbs(x <- rnorm(1001)) + +} +\author{ +Matthew L. Fidler +} diff --git a/man/rxode2.Rd b/man/rxode2.Rd index b16f51d22..2d8108f38 100644 --- a/man/rxode2.Rd +++ b/man/rxode2.Rd @@ -327,7 +327,7 @@ compilation model. \if{html}{\out{
}}\preformatted{## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’ -## rxode2 2.0.13.9000 model named rx_dbd9d9052eeef81ccd5e07cc7df253d2 model (ready). +## rxode2 2.0.13.9000 model named rx_a80a19966fa254e55fa958b934a783f9 model (ready). ## x$state: depot, center ## x$stateExtra: cp ## x$params: tka, tcl, tv, add.sd, eta.ka, eta.cl, eta.v, rxerr.cp @@ -340,7 +340,7 @@ mod$simulationIniModel \if{html}{\out{
}}\preformatted{## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’ -## rxode2 2.0.13.9000 model named rx_bd7124e0e4de9b93afb333b63054f58d model (ready). +## rxode2 2.0.13.9000 model named rx_efba4df1b26a5a9e54ae25cd8fe7247f model (ready). ## x$state: depot, center ## x$stateExtra: cp ## x$params: tka, tcl, tv, add.sd, eta.ka, eta.cl, eta.v, rxerr.cp diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index e4b5100bf..00a4ac6d5 100644 --- a/src/RcppExports.cpp +++ b/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)"); 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; } diff --git a/src/init.c b/src/init.c index 54cfbb141..5ecee8f32 100644 --- a/src/init.c +++ b/src/init.c @@ -282,6 +282,7 @@ SEXP _rxode2_trans(SEXP parse_file, SEXP prefix, SEXP model_md5, SEXP parseStr, SEXP isEscIn, SEXP inME, SEXP goodFuns, SEXP fullPrintIn); SEXP _rxode2_assignSeedInfo(void); SEXP _rxSetSeed(SEXP); +SEXP _rxode2_meanProbs_(SEXP x, SEXP probs, SEXP naRm, SEXP useT); typedef SEXP (*lotriMat_type) (SEXP, SEXP, SEXP); typedef SEXP (*asLotriMat_type) (SEXP, SEXP, SEXP); @@ -359,6 +360,7 @@ extern SEXP chin(SEXP x, SEXP table); SEXP _rxode2_RcppExport_registerCCallable(void); void R_init_rxode2(DllInfo *info){ R_CallMethodDef callMethods[] = { + {"_rxode2_meanProbs_", (DL_FUNC) &_rxode2_meanProbs_, 4}, {"_rxode2_getEtRxsolve", (DL_FUNC) &_rxode2_getEtRxsolve, 1}, {"_rxode2_assignSeedInfo", (DL_FUNC) &_rxode2_assignSeedInfo, 0}, {"_rxProgress", (DL_FUNC) &_rxProgress, 2}, diff --git a/src/utilcpp.cpp b/src/utilcpp.cpp index cf17e5119..af77f7f9c 100644 --- a/src/utilcpp.cpp +++ b/src/utilcpp.cpp @@ -127,3 +127,71 @@ NumericVector rxErf(NumericVector v) { } return ret; } + +#define min2( a , b ) ( (a) < (b) ? (a) : (b) ) +#define max2( a , b ) ( (a) > (b) ? (a) : (b) ) + +//[[Rcpp::export]] +NumericVector meanProbs_(NumericVector x, NumericVector probs, bool naRm, bool useT) { + double oldM = 0.0, newM = 0.0, + oldS = 0.0, newS = 0.0, mx=R_NegInf, mn=R_PosInf; + int n = 0; + for (int i = 0; i < x.size(); ++i) { + double cur = x[i]; + if (ISNA(cur)) { + if (naRm) { + continue; + } else { + NumericVector ret(6+probs.size()); + for (int j = 0; j < ret.size(); ++j) { + ret[j] = NA_REAL; + } + return ret; + } + } + n++; + mn = min2(cur, mn); + mx = max2(cur, mx); + if (n == 1) { + oldM = newM = cur; + oldS = 0.0; + } else { + newM = oldM + (cur - oldM)/n; + newS = oldS + (cur - oldM)*(cur - newM); + oldM = newM; + oldS = newS; + } + } + double var; + if (n > 1) { + var = newS/(n-1); + } else { + var = 0.0; + } + double sd = sqrt(var); + // mean, var, sd, min, max + NumericVector ret(6+probs.size()); + ret[0] = oldM; + ret[1] = var; + ret[2] = sd; + ret[3] = mn; + ret[4] = mx; + ret[5] = (double)n; + double c = sd/sqrt((double)(n)); + for (int i = 0; i < probs.size(); ++i) { + double p = probs[i]; + std::string str = std::to_string(p*100) + "%"; + if (p == 0) { + ret[i+6] = mn; + } else if (p == 1) { + ret[i+6] = mx; + } else if (p == 0.5) { + ret[i+6] = oldM; + } else if (useT) { + ret[i+6] = oldM + c * Rf_qt(p, (double)(n-1), 1, 0); + } else { + ret[i+6] = oldM + c * Rf_qnorm5(p, 0.0, 1.0, 1, 0); + } + } + return ret; +} diff --git a/tests/testthat/test-meanProb.R b/tests/testthat/test-meanProb.R new file mode 100644 index 000000000..322c539a3 --- /dev/null +++ b/tests/testthat/test-meanProb.R @@ -0,0 +1,25 @@ +test_that("test meanProb()", { + x <- rnorm(10) + m <- mean(x) + s <- sd(x) + v <- var(x) + mn <- min(x) + mx <- max(x) + t1 <- c("0%"=mn, "25%"=m+(s/sqrt(10))*qt(0.25, 9), "50%"=m, "75%"=m+(s/sqrt(10))*qt(0.75, 9), "100%"=mx) + expect_equal(meanProbs(x), t1) + t2 <- c(c("mean"=m, "var"=v, "sd"=s, "min"=mn, "max"=mx, "n"=10), t1) + expect_equal(meanProbs(x, onlyProbs=FALSE), t2) + x2 <- c(x, NA_real_) + setNames(rep(NA_real_, length(t1)),names(t1)) + + expect_equal(meanProbs(x2), setNames(rep(NA_real_, length(t1)),names(t1))) + + expect_equal(meanProbs(x2, onlyProbs=FALSE), + setNames(rep(NA_real_, length(t2)),names(t2))) + + expect_equal(meanProbs(x2,na.rm=TRUE), t1) + + expect_equal(meanProbs(x2, onlyProbs=FALSE, na.rm=TRUE), t2) + + +})