diff --git a/R/RcppExports.R b/R/RcppExports.R index d27334da0..7b91d96bd 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -589,8 +589,8 @@ rxErf <- function(v) { .Call(`_rxode2_rxErf`, v) } -binomProbs_ <- function(x, probs, naRm, pred, nIn, mIn) { - .Call(`_rxode2_binomProbs_`, x, probs, naRm, pred, nIn, mIn) +binomProbs_ <- function(x, probs, naRm, nIn, cont) { + .Call(`_rxode2_binomProbs_`, x, probs, naRm, nIn, cont) } meanProbs_ <- function(x, probs, naRm, useT, pred, nIn) { diff --git a/R/confint.R b/R/confint.R index fd5356151..078e7ad3e 100644 --- a/R/confint.R +++ b/R/confint.R @@ -63,6 +63,10 @@ confint.rxSolve <- function(object, parm = NULL, level = 0.95, ...) { if (any(names(.args) == "m")) { .m <- unique(.args$m) } + .method <- "wald" + if (any(names(.args) == "method")) { + .method <- .args$method + } .stk <- rxStack(object, parm, doSim=.doSim) if (!any(names(.stk) == "id") && any(names(.stk) == "sim.id")) { @@ -112,7 +116,7 @@ confint.rxSolve <- function(object, parm = NULL, level = 0.95, ...) { } else if (.binom) { .stk <- .stk[, list( p1 = .p, eff = rxode2::binomProbs(.SD$value, probs = .p, na.rm = TRUE, - n=.nC, m=.m, pred=.pred), + n=.nC, method=.method), Percentile = sprintf("%s%%", .p * 100) ), by = c("time", "trt", .by) @@ -154,7 +158,7 @@ confint.rxSolve <- function(object, parm = NULL, level = 0.95, ...) { print(.nC) .ret <- .ret[, list(p1 = .p, eff = rxode2::binomProbs(.SD$value, probs = .p, na.rm = TRUE, - n=.nC, m=.m, pred=.pred)), + n=.nC, method=.method)), by = c("id", "time", "trt", .by)] } else { .ret <- .ret[, list(p1 = .p, diff --git a/R/utils.R b/R/utils.R index 98c6f3784..8c99f8f48 100644 --- a/R/utils.R +++ b/R/utils.R @@ -879,14 +879,6 @@ meanProbs.default <- function(x, probs=seq(0, 1, 0.25), na.rm=FALSE, #' Calculate expected confidence bands with binomial sampling distribution #' -#' The generic function `binomProbs` produces expected confidence bands -#' using the Wilson score interval with continuity correction -#' (described in Newcombe 1998). -#' -#' The expected prediction intervals use the Frequentist Nelson's -#' predictional interval (described by Lu 2020). -#' -#' #' 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. @@ -896,33 +888,55 @@ meanProbs.default <- function(x, probs=seq(0, 1, 0.25), na.rm=FALSE, #' #' @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’. +#' 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) +#' 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/prediction interval estimates, -#' otherwise return extra statistics +#' otherwise return extra statistics. +#' #' @param n integer/integerish; this is the n used to calculate the #' prediction or confidence interval. When `n=0` (default) use the #' number of non-`NA` observations. When calculating the prediction #' interval, this represents the number of observations used in the #' input ("true") distribution. -#' @param m integer/integerish; this is the m used to calculate the -#' prediction interval. It represents the expected number of -#' samples in the new statistic. When `1=0` (default) use the -#' number of non-`NA` observations in the input `x`. +#' +#' @param method gives the method for calculating the confidence +#' interval. +#' +#' Can be: +#' +#' - "argestiCoull" -- Agresti-Coull method. For a 95\% confidence +#' interval, this method does not use the concept of "adding 2 +#' successes and 2 failures," but rather uses the formulas explicitly +#' described in the following link: +#' +#' https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval#Agresti-Coull_Interval. +#' +#' - "wilson" -- Wilson Method +#' +#' - "wilsonCorrect" -- Wilson method with continuity correction +#' +#' - "wald" -- Wald confidence interval or standard z approximation. +#' #' @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 @@ -962,19 +976,19 @@ binomProbs <- function(x, ...) { #' @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, pred=FALSE, n=0L, m=0L) { + names=TRUE, onlyProbs=TRUE, n=0L, + method=c("wilson", "wilsonCorrect", "agrestiCoull", "wald")) { + method <- match.arg(method) + method <- setNames(c("wilson"=1L, "wilsonCorrect"=0L, "agrestiCoull"=3, "wald"=2L)[method], NULL) checkmate::assertNumeric(x, min.len=1, lower=0.0, upper=1.0) x <- as.double(x) checkmate::assertIntegerish(n, min.len=1, lower=0, any.missing=FALSE) - checkmate::assertIntegerish(m, min.len=1, lower=0, any.missing=FALSE) n <- as.integer(n) - m <- as.integer(m) 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) - checkmate::assertLogical(pred, any.missing=FALSE, len=1) - .ret <- .Call(`_rxode2_binomProbs_`, x, probs, na.rm, pred, n, m) + .ret <- .Call(`_rxode2_binomProbs_`, x, probs, na.rm, n, method) .names <- NULL if (names) { .names <- paste0(probs*100, "%") diff --git a/inst/include/rxode2_RcppExports.h b/inst/include/rxode2_RcppExports.h index 44318dfb7..31115c902 100644 --- a/inst/include/rxode2_RcppExports.h +++ b/inst/include/rxode2_RcppExports.h @@ -1073,17 +1073,17 @@ namespace rxode2 { return Rcpp::as(rcpp_result_gen); } - inline NumericVector binomProbs_(NumericVector x, NumericVector probs, bool naRm, bool pred, int nIn, int mIn) { - typedef SEXP(*Ptr_binomProbs_)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP); + inline NumericVector binomProbs_(NumericVector x, NumericVector probs, bool naRm, int nIn, int cont) { + typedef SEXP(*Ptr_binomProbs_)(SEXP,SEXP,SEXP,SEXP,SEXP); static Ptr_binomProbs_ p_binomProbs_ = NULL; if (p_binomProbs_ == NULL) { - validateSignature("NumericVector(*binomProbs_)(NumericVector,NumericVector,bool,bool,int,int)"); + validateSignature("NumericVector(*binomProbs_)(NumericVector,NumericVector,bool,int,int)"); p_binomProbs_ = (Ptr_binomProbs_)R_GetCCallable("rxode2", "_rxode2_binomProbs_"); } RObject rcpp_result_gen; { RNGScope RCPP_rngScope_gen; - rcpp_result_gen = p_binomProbs_(Shield(Rcpp::wrap(x)), Shield(Rcpp::wrap(probs)), Shield(Rcpp::wrap(naRm)), Shield(Rcpp::wrap(pred)), Shield(Rcpp::wrap(nIn)), Shield(Rcpp::wrap(mIn))); + rcpp_result_gen = p_binomProbs_(Shield(Rcpp::wrap(x)), Shield(Rcpp::wrap(probs)), Shield(Rcpp::wrap(naRm)), Shield(Rcpp::wrap(nIn)), Shield(Rcpp::wrap(cont))); } if (rcpp_result_gen.inherits("interrupted-error")) throw Rcpp::internal::InterruptedException(); diff --git a/man/binomProbs.Rd b/man/binomProbs.Rd index 84e52ecfd..09e5b2bc8 100644 --- a/man/binomProbs.Rd +++ b/man/binomProbs.Rd @@ -13,15 +13,14 @@ binomProbs(x, ...) na.rm = FALSE, names = TRUE, onlyProbs = TRUE, - pred = FALSE, n = 0L, - m = 0L + method = c("wilson", "wilsonCorrect", "agrestiCoull", "wald") ) } \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’.} +vectors unless \code{na.rm} is \code{TRUE}.} \item{...}{Arguments passed to default method, allows many different methods to be applied.} @@ -29,7 +28,7 @@ different methods to be applied.} \item{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)} +represents the expected probability (mean).} \item{na.rm}{logical; if true, any NA and NaN's are removed from \code{x} before the quantiles are computed.} @@ -38,7 +37,7 @@ represents the expected probability (mean)} \item{onlyProbs}{logical; if true, only return the probability based confidence interval/prediction interval estimates, -otherwise return extra statistics} +otherwise return extra statistics.} \item{n}{integer/integerish; this is the n used to calculate the prediction or confidence interval. When \code{n=0} (default) use the @@ -46,10 +45,23 @@ number of non-\code{NA} observations. When calculating the prediction interval, this represents the number of observations used in the input ("true") distribution.} -\item{m}{integer/integerish; this is the m used to calculate the -prediction interval. It represents the expected number of -samples in the new statistic. When \code{1=0} (default) use the -number of non-\code{NA} observations in the input \code{x}.} +\item{method}{gives the method for calculating the confidence +interval. + +Can be: +\itemize{ +\item "argestiCoull" -- Agresti-Coull method. For a 95\\% confidence +interval, this method does not use the concept of "adding 2 +successes and 2 failures," but rather uses the formulas explicitly +described in the following link: +} + +https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval#Agresti-Coull_Interval. +\itemize{ +\item "wilson" -- Wilson Method +\item "wilsonCorrect" -- Wilson method with continuity correction +\item "wald" -- Wald confidence interval or standard z approximation. +}} } \value{ By default the return has the probabilities as names (if @@ -59,18 +71,11 @@ located given the sampling mean and standard deviation. If deviation, minimum, maximum and number of non-NA observations. } \description{ -The generic function \code{binomProbs} produces expected confidence bands -using the Wilson score interval with continuity correction -(described in Newcombe 1998). -} -\details{ -The expected prediction intervals use the Frequentist Nelson's -predictional interval (described by Lu 2020). - This is meant to perform in the same way as \code{quantile()} so it can be a drop in replacement for code using \code{quantile()} but using distributional assumptions. - +} +\details{ It is used for confidence intervals with rxode2 solved objects using \code{confint(mean="binom")} } diff --git a/man/meanProbs.Rd b/man/meanProbs.Rd index c250915ef..9e5bffb65 100644 --- a/man/meanProbs.Rd +++ b/man/meanProbs.Rd @@ -3,7 +3,7 @@ \name{meanProbs} \alias{meanProbs} \alias{meanProbs.default} -\title{Calculate expected confidence bands with normal or t sampling distribution} +\title{Calculate expected confidence bands or prediction intreval with normal or t sampling distribution} \usage{ meanProbs(x, ...) diff --git a/man/rxode2.Rd b/man/rxode2.Rd index da9339e0a..4c8582f5b 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_62be9c773d9bf813d22d7e58e66428af model (ready). +## rxode2 2.0.13.9000 model named rx_0b39f807dc809326dc25830c32ea3aae 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_6585156186f2c1e9b419c439d782b6d2 model (ready). +## rxode2 2.0.13.9000 model named rx_c0dcddf06898af4878966b21750d95b6 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 0f3483ea8..01889cea3 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -1825,25 +1825,24 @@ RcppExport SEXP _rxode2_rxErf(SEXP vSEXP) { return rcpp_result_gen; } // binomProbs_ -NumericVector binomProbs_(NumericVector x, NumericVector probs, bool naRm, bool pred, int nIn, int mIn); -static SEXP _rxode2_binomProbs__try(SEXP xSEXP, SEXP probsSEXP, SEXP naRmSEXP, SEXP predSEXP, SEXP nInSEXP, SEXP mInSEXP) { +NumericVector binomProbs_(NumericVector x, NumericVector probs, bool naRm, int nIn, int cont); +static SEXP _rxode2_binomProbs__try(SEXP xSEXP, SEXP probsSEXP, SEXP naRmSEXP, SEXP nInSEXP, SEXP contSEXP) { 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 pred(predSEXP); Rcpp::traits::input_parameter< int >::type nIn(nInSEXP); - Rcpp::traits::input_parameter< int >::type mIn(mInSEXP); - rcpp_result_gen = Rcpp::wrap(binomProbs_(x, probs, naRm, pred, nIn, mIn)); + Rcpp::traits::input_parameter< int >::type cont(contSEXP); + rcpp_result_gen = Rcpp::wrap(binomProbs_(x, probs, naRm, nIn, cont)); return rcpp_result_gen; END_RCPP_RETURN_ERROR } -RcppExport SEXP _rxode2_binomProbs_(SEXP xSEXP, SEXP probsSEXP, SEXP naRmSEXP, SEXP predSEXP, SEXP nInSEXP, SEXP mInSEXP) { +RcppExport SEXP _rxode2_binomProbs_(SEXP xSEXP, SEXP probsSEXP, SEXP naRmSEXP, SEXP nInSEXP, SEXP contSEXP) { SEXP rcpp_result_gen; { Rcpp::RNGScope rcpp_rngScope_gen; - rcpp_result_gen = PROTECT(_rxode2_binomProbs__try(xSEXP, probsSEXP, naRmSEXP, predSEXP, nInSEXP, mInSEXP)); + rcpp_result_gen = PROTECT(_rxode2_binomProbs__try(xSEXP, probsSEXP, naRmSEXP, nInSEXP, contSEXP)); } Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error"); if (rcpp_isInterrupt_gen) { @@ -1957,7 +1956,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(*binomProbs_)(NumericVector,NumericVector,bool,bool,int,int)"); + signatures.insert("NumericVector(*binomProbs_)(NumericVector,NumericVector,bool,int,int)"); signatures.insert("NumericVector(*meanProbs_)(NumericVector,NumericVector,bool,bool,bool,int)"); } return signatures.find(sig) != signatures.end(); diff --git a/src/init.c b/src/init.c index 49597c71f..8e1d0ec5f 100644 --- a/src/init.c +++ b/src/init.c @@ -283,7 +283,7 @@ SEXP _rxode2_trans(SEXP parse_file, SEXP prefix, SEXP model_md5, SEXP parseStr, SEXP _rxode2_assignSeedInfo(void); SEXP _rxSetSeed(SEXP); SEXP _rxode2_meanProbs_(SEXP x, SEXP probs, SEXP naRm, SEXP useT, SEXP pred, SEXP inN); -SEXP _rxode2_binomProbs_(SEXP x, SEXP probs, SEXP naRm, SEXP pred, SEXP inN, SEXP inM); +SEXP _rxode2_binomProbs_(SEXP x, SEXP probs, SEXP naRm, SEXP inN, SEXP cont); typedef SEXP (*lotriMat_type) (SEXP, SEXP, SEXP); typedef SEXP (*asLotriMat_type) (SEXP, SEXP, SEXP); @@ -361,7 +361,7 @@ extern SEXP chin(SEXP x, SEXP table); SEXP _rxode2_RcppExport_registerCCallable(void); void R_init_rxode2(DllInfo *info){ R_CallMethodDef callMethods[] = { - {"_rxode2_binomProbs_", (DL_FUNC) &_rxode2_binomProbs_, 6}, + {"_rxode2_binomProbs_", (DL_FUNC) &_rxode2_binomProbs_, 5}, {"_rxode2_meanProbs_", (DL_FUNC) &_rxode2_meanProbs_, 6}, {"_rxode2_getEtRxsolve", (DL_FUNC) &_rxode2_getEtRxsolve, 1}, {"_rxode2_assignSeedInfo", (DL_FUNC) &_rxode2_assignSeedInfo, 0}, diff --git a/src/utilcpp.cpp b/src/utilcpp.cpp index caa80c22f..a2dc45b9f 100644 --- a/src/utilcpp.cpp +++ b/src/utilcpp.cpp @@ -132,26 +132,10 @@ NumericVector rxErf(NumericVector v) { #define max2( a , b ) ( (a) > (b) ? (a) : (b) ) //[[Rcpp::export]] -NumericVector binomProbs_(NumericVector x, NumericVector probs, bool naRm, bool pred, - int nIn, int mIn) { +NumericVector binomProbs_(NumericVector x, NumericVector probs, bool naRm, int nIn, + int cont) { double oldM = 0.0, newM = 0.0; int n = 0; - // For CI use Wilson score interval with continuity correction. - // Described in 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. - // For PI use the Frequentist PI described in Hezhi Lu, Hua Jin, - // A new prediction interval for binomial random variable based on inferential models, - // Journal of Statistical Planning and Inference, - // Volume 205, - // 2020, - // Pages 156-174, - // ISSN 0378-3758, - // https://doi.org/10.1016/j.jspi.2019.07.001. - // in that case we are using the Nelson's prediction interval - // With the case where m=n=1 because this is the bernoulli case for (int i = 0; i < x.size(); ++i) { double cur = x[i]; if (ISNA(cur)) { @@ -179,12 +163,6 @@ NumericVector binomProbs_(NumericVector x, NumericVector probs, bool naRm, bool } else { nC = nIn; } - int m = 1; - if (mIn == 0) { - m = n; - } else { - m = mIn; - } double var=0; double sd=0; if (n > 1) { @@ -217,41 +195,61 @@ NumericVector binomProbs_(NumericVector x, NumericVector probs, bool naRm, bool } else { z = Rf_qnorm5(1.0-p, 0.0, 1.0, 1, 0); } - if (pred) { -#define y oldM - double z2 = z*z; - double A = m*nC*(2*y*z2*(nC+z2+m)+(2*y+z2)/((m+nC)*(m+nC))); - double B1 = m*nC*(m+nC)*z2*(m+nC+z2)*(m+nC+z2); - double B2 = 2*(nC-y)*(nC*nC*(2*y+z2)+4.0*m*nC*y+2.0*m*m*y); - double B3 = nC*z2*(nC*(2*y+z2)+3.0*m*nC+m*m); - double B = sqrt(B1*(B2+B3)); - double C1 = (nC+z2)*(m*m+nC*(nC+z2)); - double C2 = m*nC*(2*nC+3*z2); - double C = 2*nC*(C1+C2); - if (p < 0.5) { - ret[i+4] = max2(0.0, (A-B)/C); - } else { - ret[i+4] = min2(1.0, (A+B)/C); - } -#undef y - } else { - double z2 = z*z; - double coef1 = 2*nC*oldM + z2; - double coef2 = z*sqrt(z2 - 1.0/nC + 4*nC*var + (4*oldM-2))+1.0; - double coef3 = 2*(nC+z2); + double z2 = z*z; + double c1, c2, c3; +#define contWilsonContinuityCorrection 0 +#define contWilson 1 +#define contWald 2 +#define contAgrestiCoull 3 + switch (cont) { + case contWilsonContinuityCorrection: + c1 = 2*nC*oldM + z2; + c2 = z*sqrt(z2 - 1.0/nC + 4*nC*var + (4*oldM-2))+1.0; + c3 = 2*(nC+z2); if (p < 0.5) { if (p == 0.0) { ret[i+4] = 0.0; } else { - ret[i+4] = max2(0, (coef1-coef2)/coef3); + ret[i+4] = max2(0, (c1-c2)/c3); } } else { if (p == 1.0) { ret[i+4] = 1.0; } else { - ret[i+4] = min2(1, (coef1+coef2)/coef3); + ret[i+4] = min2(1, (c1+c2)/c3); } } + break; + case contWilson: + c1 = 1.0/(1.0+z2/nC)*(oldM+z2/(2.0*nC)); + c2 = z/(1.0+z2/nC)*sqrt(oldM*(1.0-oldM)/nC + z2/(4.0*nC*nC)); + if (p < 0.5) { + ret[i+4] = c1-c2; + } else { + ret[i+4] = c1+c2; + } + break; + case contWald: + c1 = oldM; + c2 = z*sqrt(oldM*(1-oldM)/nC); + if (p < 0.5) { + ret[i+4] = c1-c2; + } else { + ret[i+4] = c1+c2; + } + break; + case contAgrestiCoull: + // p = ns/n; p*n = ns + c1 = oldM*nC; + nC = nC + z2; + c1 = 1.0/nC*(c1+z2/2.0); + c2 = z*sqrt(c1/nC*(1-c1)); + if (p < 0.5) { + ret[i+4] = c1-c2; + } else { + ret[i+4] = c1+c2; + } + break; } } } diff --git a/tests/testthat/test-binomProb.R b/tests/testthat/test-binomProb.R index 8b5f10480..f389249d9 100644 --- a/tests/testthat/test-binomProb.R +++ b/tests/testthat/test-binomProb.R @@ -27,19 +27,20 @@ test_that("test binomProb()", { "50%"=m, "95%"=z5[2], "97.5%"=z2.5[2]) - expect_equal(binomProbs(x), t1) + expect_equal(binomProbs(x, method="wilsonCorrect"), t1) t2 <- c(c("mean"=m, "var"=v, "sd"=s, "n"=n), t1) - expect_equal(binomProbs(x, onlyProbs=FALSE), t2) + expect_equal(binomProbs(x, onlyProbs=FALSE, method="wilsonCorrect"), t2) x2 <- c(x, NA_real_) setNames(rep(NA_real_, length(t1)),names(t1)) - expect_equal(binomProbs(x2), setNames(rep(NA_real_, length(t1)),names(t1))) + expect_equal(binomProbs(x2, method="wilsonCorrect"), + setNames(rep(NA_real_, length(t1)),names(t1))) - expect_equal(binomProbs(x2, onlyProbs=FALSE), + expect_equal(binomProbs(x2, onlyProbs=FALSE, method="wilsonCorrect"), setNames(rep(NA_real_, length(t2)),names(t2))) - expect_equal(binomProbs(x2,na.rm=TRUE), t1) - expect_equal(binomProbs(x2, onlyProbs=FALSE, na.rm=TRUE), t2) + expect_equal(binomProbs(x2,na.rm=TRUE, method="wilsonCorrect"), t1) + expect_equal(binomProbs(x2, onlyProbs=FALSE, na.rm=TRUE, method="wilsonCorrect"), t2) n <- 50 @@ -65,8 +66,6 @@ test_that("test binomProb()", { "95%"=z5[2], "97.5%"=z2.5[2]) - expect_equal(binomProbs(x, n=50), t1) - - ## now get test binomial predictions + expect_equal(binomProbs(x, n=50, method="wilsonCorrect"), t1) })