From d53e6f559bf672e4ec367afafaa4bf32aef7ddcf Mon Sep 17 00:00:00 2001 From: Sebastian Weber Date: Wed, 20 Nov 2024 15:12:30 +0100 Subject: [PATCH] prep for 1.7.4 release --- DESCRIPTION | 9 +-- Makefile | 79 ++++++++++++++++++--- NEWS.md | 16 +++++ R/dlink.R | 1 + R/integrate_logit_log.R | 32 +++++---- R/mix.R | 40 +++++++---- R/mixcombine.R | 10 +-- R/mixess.R | 62 ++-------------- R/mixplot.R | 4 +- R/mixstanvar.R | 27 ++++++- R/predict_gMAP.R | 36 +--------- R/zzz.R | 1 + configure | 1 - demo/robustMAP.R | 24 +++---- inst/sbc/calibration.md5 | 6 +- inst/sbc/calibration.rds | Bin 9811 -> 9861 bytes inst/sbc/make_reference_rankhist.R | 13 +++- inst/sbc/sbc_tools.R | 67 ++--------------- tests/testthat/test-EM.R | 30 ++++---- tests/testthat/test-gMAP.R | 1 - tests/testthat/test-mixdiff.R | 29 ++++---- tests/testthat/test-mixdist.R | 65 ++++++++++------- tests/testthat/test-mixstanvar.R | 73 +++++++++++++------ tests/testthat/test-oc1S.R | 39 +++++----- tests/testthat/test-oc2S.R | 52 +++++++++----- tests/testthat/test-pos1S.R | 11 ++- tests/testthat/test-pos2S.R | 16 +++++ tests/testthat/test-postmix.R | 1 - tests/testthat/test-preddist.R | 13 ++-- tests/testthat/test-utils.R | 74 +++++++++++-------- tools/make-ds.R | 2 +- vignettes/articles/introduction_normal.Rmd | 32 +++++---- vignettes/introduction.Rmd | 16 +++-- 33 files changed, 470 insertions(+), 412 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 5315e06..b1eb6c6 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -8,8 +8,8 @@ Description: Tool-set to support Bayesian evidence synthesis. This for details on applying this package while Neuenschwander et al. (2010) and Schmidli et al. (2014) explain details on the methodology. -Version: 1.7-3 -Date: 2024-01-02 +Version: 1.7-4 +Date: 2024-11-19 Authors@R: c(person("Novartis", "Pharma AG", role = "cph") ,person("Sebastian", "Weber", email="sebastian.weber@novartis.com", role=c("aut", "cre")) ,person("Beat", "Neuenschwander", email="beat.neuenschwander@novartis.com", role="ctb") @@ -26,7 +26,7 @@ Imports: Rcpp (>= 0.12.0), RcppParallel (>= 5.0.1), rstan (>= 2.26.0), - rstantools (>= 2.3.1), + rstantools (>= 2.4.0), assertthat, mvtnorm, Formula, @@ -71,4 +71,5 @@ Suggests: VignetteBuilder: knitr SystemRequirements: GNU make, pandoc (>= 1.12.3), pngquant, C++17 Encoding: UTF-8 -RoxygenNote: 7.2.3 +RoxygenNote: 7.3.1 +Config/testthat/edition: 3 diff --git a/Makefile b/Makefile index 2c955c4..59dfa27 100644 --- a/Makefile +++ b/Makefile @@ -15,6 +15,9 @@ RPKG=$(patsubst ‘%’, %, $(word 2, $(shell grep ^Package: DESCRIPTION))) INCS = R_PKG_SRCS = $(wildcard R/*.R) R_SRCS = $(wildcard *.R $(foreach fd, $(SRCDIR), $(fd)/*.R)) +R_TEST_SRCS = $(wildcard tests/testthat/test*.R) +R_TEST_OBJS = $(R_TEST_SRCS:.R=.Rtest) +R_TESTFAST_OBJS = $(R_TEST_SRCS:.R=.Rtestfast) RMD_SRCS = $(wildcard *.Rmd $(foreach fd, $(SRCDIR), $(fd)/x*.Rmd)) STAN_SRCS = $(wildcard *.stan $(foreach fd, $(SRCDIR), $(fd)/*.stan)) SRCS = $(R_PKG_SRCS) $(R_SRCS) $(RMD_SRCS) $(STAN_SRCS) @@ -28,6 +31,7 @@ PKG_VERSION ?= $(patsubst ‘%’, %, $(word 2, $(shell grep ^Version DESCRIPTIO GIT_TAG ?= v$(PKG_VERSION) MD5 ?= md5sum +TMPDIR := $(realpath $(shell mktemp -d)) all : $(TARGET) @@ -45,31 +49,47 @@ all : $(TARGET) cd $(@D); echo running $(RCMD) -e "rmarkdown::render('$( $@ 2>&1 + @printf "Test summary for $( $@ 2>&1 + @printf "Test summary for $( NAMESPACE echo "import(Rcpp)" >> NAMESPACE echo "import(methods)" >> NAMESPACE echo "importFrom(rstan, sampling)" >> NAMESPACE echo "useDynLib($(RPKG), .registration = TRUE)" >> NAMESPACE + install -d src "${R_HOME}/bin/Rscript" -e "rstantools::rstan_config()" + touch R/stanmodels.R -src/package-binary: $(STAN_SRCS) R/stanmodels.R +src/package-binary: R/stanmodels.R ## ensure that NAMESPACE contains load directive echo "# Generated by roxygen2: do not edit by hand" > NAMESPACE echo "import(Rcpp)" >> NAMESPACE echo "import(methods)" >> NAMESPACE echo "importFrom(rstan, sampling)" >> NAMESPACE echo "useDynLib($(RPKG), .registration = TRUE)" >> NAMESPACE - "${R_HOME}/bin/R" --slave -e 'library(pkgbuild); pkgbuild::compile_dll()' + install -d src + "${R_HOME}/bin/Rscript" -e 'pkgbuild::compile_dll()' touch src/package-binary man/package-doc: $(R_PKG_SRCS) $(BIN_OBJS) - "${R_HOME}/bin/R" --slave -e 'library(roxygen2); roxygen2::roxygenize()' + "${R_HOME}/bin/Rscript" -e 'roxygen2::roxygenize()' touch man/package-doc +inst/sbc/sbc_report.html : inst/sbc/calibration.rds + inst/sbc/calibration.rds : - @echo Please run inst/sbc/make_reference_rankhist.R + echo "Please run inst/sbc/make_reference_rankhist.R" exit 1 R/sysdata.rda: inst/sbc/calibration.rds @@ -85,12 +105,21 @@ NAMESPACE: man/package-doc PHONY := $(TARGET) -$(TARGET): build/r-source +$(TARGET): build/r-source-fast -build/r-source : $(BIN_OBJS) $(DOC_OBJS) $(SRCS) +build/r-source-fast : $(BIN_OBJS) $(DOC_OBJS) $(SRCS) install -d build - cd build; $(RCMD) CMD build .. --no-build-vignettes --no-manual - cd build; mv $(RPKG)_$(PKG_VERSION).tar.gz $(RPKG)-source.tar.gz + git archive --format=tar.gz --prefix $(RPKG)-$(GIT_TAG)/ HEAD > build/$(RPKG)-$(GIT_TAG).tar.gz + rm -rf build/$(RPKG)-$(GIT_TAG) + cd build; tar x -C $(TMPDIR) -f $(RPKG)-$(GIT_TAG).tar.gz + rm -f build/$(RPKG)-$(GIT_TAG).tar.gz + cp -v NAMESPACE $(TMPDIR)/$(RPKG)-$(GIT_TAG) + install -d $(TMPDIR)/$(RPKG)-$(GIT_TAG)/man + cp -v man/*.Rd $(TMPDIR)/$(RPKG)-$(GIT_TAG)/man + cd $(TMPDIR)/$(RPKG)-$(GIT_TAG); "${R_HOME}/bin/R" --slave --file=tools/make-ds.R + cd $(TMPDIR); NOT_CRAN=false "${R_HOME}/bin/R" CMD build $(RPKG)-$(GIT_TAG) --no-build-vignettes --no-manual + rm -rf $(TMPDIR)/$(RPKG)-$(GIT_TAG) + mv $(TMPDIR)/$(RPKG)_$(PKG_VERSION).tar.gz build/$(RPKG)-source.tar.gz touch build/r-source-fast build/r-source-release : $(BIN_OBJS) $(DOC_OBJS) $(SRCS) inst/sbc/sbc_report.html @@ -124,8 +153,28 @@ derived : NAMESPACE $(BIN_OBJS) $(DOC_OBJS) PHONY += r-source-check r-source-check : r-source - cd build; tar xvzf RBesT-source.tar.gz - cd build; NOT_CRAN=true $(RCMD) CMD check RBesT + cd build; tar xvzf $(RPKG)-source.tar.gz + cd build; NOT_CRAN=true $(RCMD) CMD check $(RPKG) + +build/installed/$(RPKG)/DESCRIPTION : build/r-source-fast + rm -rf build/installed + install -d build/installed + cd build; $(RCMD) CMD INSTALL --library=./installed --no-docs --no-multiarch --no-test-load --no-clean-on-error $(RPKG)-source.tar.gz + +PHONY += dev-install +dev-install: build/installed/$(RPKG)/DESCRIPTION + +PHONY += test-all +test-all : $(R_TEST_OBJS) + +PHONY += testfast-all +testfast-all : $(R_TESTFAST_OBJS) + +PHONY += retestfast-all +retestfast-all : clean-test $(R_TESTFAST_OBJS) + +PHONY += retest-all +retest-all : clean-test $(R_TEST_OBJS) #$(DIR_OBJ)/%.o: %.c $(INCS) # mkdir -p $(@D) @@ -146,6 +195,14 @@ clean: rm -f vignettes/*.html rm -f vignettes/*.docx rm -rf .Rd2pdf* + rm -f $(R_TEST_OBJS) + rm -f $(R_TESTFAST_OBJS) + rm -rf src + rm -f R/stanmodels.R + +clean-test: + rm -f $(R_TEST_OBJS) + rm -f $(R_TESTFAST_OBJS) PHONY += doc doc: $(DOC_OBJS) diff --git a/NEWS.md b/NEWS.md index f08c827..a2013eb 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,19 @@ +# RBesT 1.7-4 - November 20th, 2024 + +## Enhancements + +* Upgraded testhat edition to version 3 +* Added for `mixstanvar` automatic generation of distribution + functions allowing truncated priors in `brms` with mixtures + +## Bugfixes + +* Fixed a rare issue when estimating ESS ELIR for beta mixtures. The + calculation aborted due to boundary issues which are now avoided. +* Avoid over and underflow of beta-binomial distribution. +* Replace calls to deprecated `ggplot2::qplot` with respective calls + to `ggplot2::ggplot`. + # RBesT 1.7-3 - January 2nd, 2024 ## Enhancements diff --git a/R/dlink.R b/R/dlink.R index 19d6653..7f1aac4 100644 --- a/R/dlink.R +++ b/R/dlink.R @@ -97,5 +97,6 @@ is.dlink <- function(x) is.dlink_identity <- function(x) is.dlink(x) & x$name == "identity" +#' @export print.dlink <- function(x, ...) cat("Link:", x$name, "\n") diff --git a/R/integrate_logit_log.R b/R/integrate_logit_log.R index 651be98..653ff4f 100644 --- a/R/integrate_logit_log.R +++ b/R/integrate_logit_log.R @@ -16,8 +16,8 @@ #' @param Lpupper logit of upper cumulative density #' #' @keywords internal -integrate_density_log <- function(log_integrand, mix, Lplower=-Inf, Lpupper=Inf, eps=1E-9) { - .integrand_comp <- function(mix_comp) { +integrate_density_log <- function(log_integrand, mix, Lplower=-Inf, Lpupper=Inf, eps=getOption("RBesT.integrate_prob_eps", 1E-6)) { + .integrand_comp_logit <- function(mix_comp) { function(l) { u <- inv_logit(l) lp <- log_inv_logit(l) @@ -37,10 +37,9 @@ integrate_density_log <- function(log_integrand, mix, Lplower=-Inf, Lpupper=Inf, upper <- inv_logit(Lpupper) return(sum(vapply(1:Nc, function(comp) { mix_comp <- mix[[comp, rescale=TRUE]] - mix_comp_identity <- mix_comp - dlink(mix_comp_identity) <- identity_dlink - if (all(dmix(mix_comp_identity, support(mix_comp_identity)) == 0 )) { - return(.integrate(.integrand_comp(mix_comp), Lplower, Lpupper)) + fn_integrand_comp_logit <- .integrand_comp_logit(mix_comp) + if (all(!is.na(fn_integrand_comp_logit(c(Lplower, Lpupper))))) { + return(.integrate(fn_integrand_comp_logit, Lplower, Lpupper)) } lower_comp <- ifelse(Lplower==-Inf, qmix(mix_comp, eps), qmix(mix_comp, lower)) upper_comp <- ifelse(Lpupper==Inf, qmix(mix_comp, 1-eps), qmix(mix_comp, upper)) @@ -48,8 +47,8 @@ integrate_density_log <- function(log_integrand, mix, Lplower=-Inf, Lpupper=Inf, }, c(0.1)) * mix[1,])) } -integrate_density <- function(integrand, mix, Lplower=-Inf, Lpupper=Inf, eps=1E-9) { - .integrand_comp <- function(mix_comp) { +integrate_density <- function(integrand, mix, Lplower=-Inf, Lpupper=Inf, eps=getOption("RBesT.integrate_prob_eps", 1E-6)) { + .integrand_comp_logit <- function(mix_comp) { function(l) { u <- inv_logit(l) lp <- log_inv_logit(l) @@ -57,18 +56,19 @@ integrate_density <- function(integrand, mix, Lplower=-Inf, Lpupper=Inf, eps=1E- exp(lp + lnp) * integrand(qmix(mix_comp, u)) } } - Nc <- ncol(mix) lower <- inv_logit(Lplower) upper <- inv_logit(Lpupper) + return(sum(vapply(1:Nc, function(comp) { mix_comp <- mix[[comp, rescale=TRUE]] - mix_comp_identity <- mix_comp - dlink(mix_comp_identity) <- identity_dlink - if (all(dmix(mix_comp_identity, support(mix_comp_identity)) == 0 )) { - return(.integrate(.integrand_comp(mix_comp), Lplower, Lpupper)) + ## ensure that the integrand is defined at the boundaries... + fn_integrand_comp_logit <- .integrand_comp_logit(mix_comp) + if (all(!is.na(fn_integrand_comp_logit(c(Lplower, Lpupper))))) { + return(.integrate(fn_integrand_comp_logit, Lplower, Lpupper)) } + ## ... otherwise we avoid the boundaries by eps prob density: lower_comp <- ifelse(Lplower==-Inf, qmix(mix_comp, eps), qmix(mix_comp, lower)) upper_comp <- ifelse(Lpupper==Inf, qmix(mix_comp, 1-eps), qmix(mix_comp, upper)) return(.integrate(function(x) integrand(x) * dmix(mix_comp, x), lower_comp, upper_comp)) @@ -80,7 +80,8 @@ integrate_density <- function(integrand, mix, Lplower=-Inf, Lpupper=Inf, eps=1E- args <- modifyList(list(lower=lower, upper=upper, rel.tol=.Machine$double.eps^0.25, abs.tol=.Machine$double.eps^0.25, - subdivisions=1000), + subdivisions=1000, + stop.on.error=TRUE), integrate_args_user) integrate(integrand, @@ -88,5 +89,6 @@ integrate_density <- function(integrand, mix, Lplower=-Inf, Lpupper=Inf, eps=1E- upper=args$upper, rel.tol=args$rel.tol, abs.tol=args$abs.tol, - subdivisions=args$subdivisions)$value + subdivisions=args$subdivisions, + stop.on.error=args$stop.on.error)$value } diff --git a/R/mix.R b/R/mix.R index 1ac1da1..5388492 100644 --- a/R/mix.R +++ b/R/mix.R @@ -222,20 +222,32 @@ pmix.normMix <- function(mix, q, lower.tail = TRUE, log.p=FALSE) pmix_impl(pnor ## pmix.betaBinomialMix <- function(mix, q, lower.tail = TRUE, log.p=FALSE) pmix_impl(Curry(pBetaBinomial, n=attr(mix, "n")), mix, q, lower.tail, log.p) pmix.betaBinomialMix <- function(mix, q, lower.tail = TRUE, log.p=FALSE) { assert_that(is.dlink_identity(attr(mix, "link"))) -## ## the analytic solution needs the generalized hypergeometric -## ## function which is only available in the hyper-geo package which -## ## is why a numeric integration is performed here -## ## treat out-of-bounds quantiles special - out_low <- q<0 - out_high <- q>attr(mix, "n") - q[out_low | out_high] <- 0 - dist <- cumsum(dmix.betaBinomialMix(mix, seq(0,max(q)))) - p <- dist[q+1] - p[out_low] <- 0 - p[out_high] <- 1 - if(!lower.tail) p <- 1-p - if(log.p) p <- log(p) - return(p) + ## ## the analytic solution needs the generalized hypergeometric + ## ## function which is only available in the hyper-geo package which + ## ## is why a numeric integration is performed here + ## ## treat out-of-bounds quantiles special +## if (requireNamespace("extraDistr", quietly = TRUE)) { +## Nq <- length(q) +## size <- attr(mix, "n") +## if(!log.p) +## return(matrixStats::colSums2(matrix(mix[1,] * extraDistr::pbbinom(rep(q, each=Nc), size, rep(mix[2,], times=Nq), rep(mix[3,], times=Nq), lower.tail=lower.tail, log.p=FALSE), nrow=Nc))) + ## log version is slower, but numerically more stable +## res <- matrixStats::colLogSumExps(matrix(log(mix[1,]) + dist(rep(q, each=Nc), rep(mix[2,], times=Nq), rep(mix[3,], times=Nq), lower.tail=lower.tail, log.p=TRUE), nrow=Nc)) +## return(res) +## } + + ## default implementation uses brute-force integration + out_low <- q<0 + out_high <- q>attr(mix, "n") + q[out_low | out_high] <- 0 + dist <- cumsum(dmix.betaBinomialMix(mix, seq(0,max(q)))) + dist <- pmin(pmax(dist, 0), 1) ## avoid over and underflows + p <- dist[q+1] + p[out_low] <- 0 + p[out_high] <- 1 + if(!lower.tail) p <- 1-p + if(log.p) p <- log(p) + return(p) } ## internal redefinition of negative binomial diff --git a/R/mixcombine.R b/R/mixcombine.R index 43033c9..c1e5310 100644 --- a/R/mixcombine.R +++ b/R/mixcombine.R @@ -19,14 +19,8 @@ #' @family mixdist #' @seealso \code{\link{robustify}} #' -#' @examples -#' # beta with two informative components -#' bm <- mixbeta(inf=c(0.5, 10, 100), inf2=c(0.5, 30, 80)) -#' -#' # robustified with mixcombine, i.e. a 10% uninformative part added -#' unif <- mixbeta(rob=c(1,1,1)) -#' mixcombine(bm, unif, weight=c(9, 1)) -#' +#' @example inst/examples/mixcombine.R +#' #' @export mixcombine <- function(..., weight, rescale=TRUE) { comp <- list(...) diff --git a/R/mixess.R b/R/mixess.R index 829d4b1..03ae16a 100644 --- a/R/mixess.R +++ b/R/mixess.R @@ -57,62 +57,8 @@ #' Predictively Consistent Prior Effective Sample Sizes. #' \emph{pre-print} 2019; arXiv:1907.04185 #' -#' @examples -#' # Conjugate Beta example -#' a <- 5 -#' b <- 15 -#' prior <- mixbeta(c(1, a, b)) -#' -#' ess(prior) -#' (a+b) -#' -#' # Beta mixture example -#' bmix <- mixbeta(rob=c(0.2, 1, 1), inf=c(0.8, 10, 2)) -#' -#' ess(bmix, "elir") -#' -#' ess(bmix, "moment") -#' # moments method is equivalent to -#' # first calculate moments -#' bmix_sum <- summary(bmix) -#' # then calculate a and b of a matching beta -#' ab_matched <- ms2beta(bmix_sum["mean"], bmix_sum["sd"]) -#' # finally take the sum of a and b which are equivalent -#' # to number of responders/non-responders respectivley -#' round(sum(ab_matched)) -#' -#' ess(bmix, method="morita") -#' -#' # Predictive consistency of elir -#' \donttest{ -#' n_forward <- 1E2 -#' bmixPred <- preddist(bmix, n=n_forward) -#' pred_samp <- rmix(bmixPred, 1E3) -#' pred_ess <- sapply(pred_samp, function(r) ess(postmix(bmix, r=r, n=n_forward), "elir") ) -#' ess(bmix, "elir") -#' mean(pred_ess) - n_forward -#' } -#' -#' # Normal mixture example -#' nmix <- mixnorm(rob=c(0.5, 0, 2), inf=c(0.5, 3, 4), sigma=10) -#' -#' ess(nmix, "elir") -#' -#' ess(nmix, "moment") -#' -#' ## the reference scale determines the ESS -#' sigma(nmix) <- 20 -#' ess(nmix) -#' -#' # Gamma mixture example -#' gmix <- mixgamma(rob=c(0.3, 20, 4), inf=c(0.7, 50, 10)) -#' -#' ess(gmix) ## interpreted as appropriate for a Poisson likelihood (default) -#' -#' likelihood(gmix) <- "exp" -#' ess(gmix) ## interpreted as appropriate for an exponential likelihood -#' -#' +#' @example inst/examples/ess.R +#' #' @export ess <- function(mix, method=c("elir", "moment", "morita"), ...) UseMethod("ess") #' @export @@ -156,12 +102,12 @@ mixInfo <- function(mix, x, dens, gradl, hessl) { dhl <- (hessl(x,a,b) + dgl^2) ## attempt numerically more stable log calculations if possible, ## i.e. if all sings are the same - if(all(dgl < 0) || all(dgl > 0)) { + if(all(!is.na(dgl)) && ( all(dgl < 0) || all(dgl > 0))) { gsum <- exp(2*matrixStats::logSumExp(lwdensComp + log(abs(dgl)))) } else { gsum <- (sum(exp(lwdensComp)*dgl))^2 } - if(all(dhl < 0) || all(dhl > 0)) { + if(all(!is.na(dhl)) && (all(dhl < 0) || all(dhl > 0))) { hsum <- sign(dhl[1]) * exp(matrixStats::logSumExp(lwdensComp + log(abs(dhl)))) } else { hsum <- (sum(exp(lwdensComp)*dhl)) diff --git a/R/mixplot.R b/R/mixplot.R index 974098e..972ef1a 100644 --- a/R/mixplot.R +++ b/R/mixplot.R @@ -12,7 +12,7 @@ #' which will display colour-coded each mixture component of the #' density in addition to the density. #' @param size controls the linesize in plots. -#' @param ... extra arguments passed on to the \code{\link[ggplot2]{qplot}} call. +#' @param ... extra arguments passed on to the plotted function. #' #' @details Plot function for mixture distribution objects. It shows #' the density/quantile/cumulative distribution (corresponds to @@ -85,7 +85,7 @@ plot.mix <- function(x, prob=0.99, fun=dmix, log=FALSE, comp=TRUE, size=1.25, .. plot_fun <- function(x, ...) fun(floor(x), ...) plot_geom <- "step" } else { - plot_fun <- fun + plot_fun <- function(x, ...) fun(x, ...) plot_geom <- "line" } n_fun <- 501 diff --git a/R/mixstanvar.R b/R/mixstanvar.R index db1756d..4894a20 100644 --- a/R/mixstanvar.R +++ b/R/mixstanvar.R @@ -154,7 +154,32 @@ real {{mixdens}}_lpdf(real y, vector w, vector {{arg1}}, vector {{arg2}}) { lp_mix[i] = {{standens}}_lpdf(y | {{arg1}}[i], {{arg2}}[i]); } return log_sum_exp(log(w) + lp_mix); -}", .open="{{", .close="}}"), block="functions") +} +real {{mixdens}}_lcdf(real y, vector w, vector {{arg1}}, vector {{arg2}}) { + int Nc = rows(w); + vector[Nc] lp_mix; + for(i in 1:Nc) { + lp_mix[i] = {{standens}}_lcdf(y | {{arg1}}[i], {{arg2}}[i]); + } + return log_sum_exp(log(w) + lp_mix); +} +real {{mixdens}}_lccdf(real y, vector w, vector {{arg1}}, vector {{arg2}}) { + int Nc = rows(w); + vector[Nc] lp_mix; + for(i in 1:Nc) { + lp_mix[i] = {{standens}}_lccdf(y | {{arg1}}[i], {{arg2}}[i]); + } + return log_sum_exp(log(w) + lp_mix); +} +real {{mixdens}}_cdf(real y, vector w, vector {{arg1}}, vector {{arg2}}) { + int Nc = rows(w); + vector[Nc] p_mix; + for(i in 1:Nc) { + p_mix[i] = {{standens}}_cdf(y | {{arg1}}[i], {{arg2}}[i]); + } + return sum(w + p_mix); +} +", .open="{{", .close="}}"), block="functions") } if(includes_density("normMix")) { sv <- sv + mixdensity_template("mixnorm", "normal", "m", "s") diff --git a/R/predict_gMAP.R b/R/predict_gMAP.R index c7cd732..b9f19c4 100644 --- a/R/predict_gMAP.R +++ b/R/predict_gMAP.R @@ -23,41 +23,7 @@ #' @seealso \code{\link{gMAP}}, \code{\link{predict.glm}} #' #' @template example-start -#' @examples -#' # create a fake data set with a covariate -#' trans_cov <- transform(transplant, country=cut(1:11, c(0,5,8,Inf), c("CH", "US", "DE"))) -#' set.seed(34246) -#' map <- gMAP(cbind(r, n-r) ~ 1 + country | study, -#' data=trans_cov, -#' tau.dist="HalfNormal", -#' tau.prior=1, -#' # Note on priors: we make the overall intercept weakly-informative -#' # and the regression coefficients must have tighter sd as these are -#' # deviations in the default contrast parametrization -#' beta.prior=rbind(c(0,2), c(0,1), c(0,1)), -#' family=binomial, -#' ## ensure fast example runtime -#' thin=1, chains=1) -#' -#' # posterior predictive distribution for each input data item (shrinkage estimates) -#' pred_cov <- predict(map) -#' pred_cov -#' -#' # extract sample as matrix -#' samp <- as.matrix(pred_cov) -#' -#' # predictive distribution for each input data item (if the input studies were new ones) -#' pred_cov_pred <- predict(map, trans_cov) -#' pred_cov_pred -#' -#' -#' # a summary function returns the results as matrix -#' summary(pred_cov) -#' -#' # obtain a prediction for new data with specific covariates -#' pred_new <- predict(map, data.frame(country="CH", study=12)) -#' pred_new -#' +#' @example inst/examples/predict_gMAP.R #' @template example-stop #' @rdname predict.gMAP #' @method predict gMAP diff --git a/R/zzz.R b/R/zzz.R index 1a480f6..083a9df 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -8,3 +8,4 @@ ver <- utils::packageVersion("RBesT") packageStartupMessage("This is RBesT version ", ver, " (released ", format(pkg_create_date, "%F"), ", git-sha ", pkg_sha, ")") } + diff --git a/configure b/configure index 7e23937..0304fc5 100755 --- a/configure +++ b/configure @@ -3,4 +3,3 @@ # Generated by rstantools. Do not edit by hand. "${R_HOME}/bin/Rscript" -e "rstantools::rstan_config()" - diff --git a/demo/robustMAP.R b/demo/robustMAP.R index 2f7e040..5b4e0da 100644 --- a/demo/robustMAP.R +++ b/demo/robustMAP.R @@ -170,15 +170,14 @@ kable(ocAdaptSamp, digits=1, caption="Sample size") #' ## Additional power Figure under varying pc #' -qplot(pt-pc, power, data=power, geom=c("line"), colour=prior) + +ggplot(power, aes(pt-pc, power, colour=prior)) + geom_line() + facet_wrap(~pc) + scale_y_continuous(breaks=seq(0,1,by=0.2)) + - scale_x_continuous(breaks=seq(-0.8,0.8,by=0.2)) + - coord_cartesian(xlim=c(-0.15,0.5)) + - geom_hline(yintercept=0.025, linetype=2) + - geom_hline(yintercept=0.8, linetype=2) + - ggtitle("Prob. for alternative for different pc") - + scale_x_continuous(breaks=seq(-0.8,0.8,by=0.2)) + + coord_cartesian(xlim=c(-0.15,0.5)) + + geom_hline(yintercept=0.025, linetype=2) + + geom_hline(yintercept=0.8, linetype=2) + + ggtitle("Prob. for alternative for different pc") #' #' ## Bias and rMSE, Figure 1 @@ -265,10 +264,8 @@ est <- foreach(case=names(map), .combine=rbind) %do% { } - -qplot(p, 100*bias, data=est, geom="line", colour=prior, main="Bias") -qplot(p, 100*rMSE, data=est, geom="line", colour=prior, main="rMSE") - +ggplot(est, aes(p, 100*bias, colour=prior)) + geom_line() + ggtitle("Bias") +ggplot(est, aes(p, 100*rMSE, colour=prior)) + geom_line() + ggtitle("rMSE") #' @@ -330,9 +327,8 @@ post <- foreach(prior=names(mapCol), .combine=rbind) %do% { res } -qplot(r, mean, data=post, colour=prior, shape=prior) + geom_abline(slope=1/20) -qplot(r, sd, data=post, colour=prior, shape=prior) + coord_cartesian(ylim=c(0,0.17)) - +ggplot(post, aes(r, mean, colour=prior, shape=prior)) + geom_point() + geom_abline(slope=1/20) +ggplot(post, aes(r, sd, colour=prior, shape=prior)) + geom_point() + coord_cartesian(ylim=c(0,0.17)) sessionInfo() diff --git a/inst/sbc/calibration.md5 b/inst/sbc/calibration.md5 index fd88f2f..e5527cc 100644 --- a/inst/sbc/calibration.md5 +++ b/inst/sbc/calibration.md5 @@ -1,3 +1,3 @@ -Created: 2023-08-03 13:56:36 UTC -git hash: 47e501da2dc04c25171cd0a1b14da6cf5115b0ed -MD5: 067f0aaa153b8c6ee2d100b34fc7dd99 +Created: 2024-11-20 13:57:46 UTC +git hash: c17af53badaeb321d3de398b84d998dff5cf685f +MD5: 95c9b0f2f899285fb53b68663ec1b6e3 diff --git a/inst/sbc/calibration.rds b/inst/sbc/calibration.rds index 8447e65197b97173848b903fc5a25e4a7ffb1a2c..970c6e14b20faaafbd46e918feee4a8416803c8f 100644 GIT binary patch literal 9861 zcmY*-t^Sdtc|D``qVq{y5jUKIh!$zGah7@cei2ZRYv}f+j=r zz=>zn3==DR9xN35dz4tB3>xH|K3vuI{{(thT9P|a>ybWw!}I2ucWkIqW7y$%*UG!M}rpf%-Yf?n0oh} zb3WdY8~QkfCjq?j;Z88hx5>!KX+dyQ=4Zrg@~_fp#?B@m0?= zi=ss0|DHud&5JkHTBjeC-Ol)OhGgjuA>=K$+m!+BN_7?Xmmx7K)(XUq#Ky5lRW^3d zi_pX3$qh(Mp_BV6i-$X)%LN8SQRJs^s&a13!4Wsz_5v|eXAsJO;b}c}M)ztThLS|`0 z!^E=dfA=b`|8q#ZwhG;tnw_3#Xa(B$sW>T=KgpncR1$|^nWBEhOw_wHUEEvs1y_gC+$=I+%Yw~pwZPQ zw|D+exY%;iHJf;dS@;xvp6)1ikQ7*P!q&`jrxkbT-&W@sKPBg))@4mj+%eDf!Y9ywDF4rut=99A?T2a?KVn& zSiT^;l|=d3d?+rB+ueO2_%kmI^WaMmLIy$v&{ro>c@U1b#o&+3j9pki+%>FKo}|Cj zBB4GC?}4XYZ$)kGw6}EM8ZTRhV(;y_1m%i}=yBVEf0{>4z_WMS70Br=SJ)3S7rNCm zJ0{sfi8h$0NQgj-IcvD7+Ir&4x^PQnbKYxl*8t?KQyFY56cwWXmiv1B(KSDuDZCf>@fG)dUa-!zrc<;!L`7wSnS;z0!q{rd^I~%SZdXov#`M<> zquJuV8d7z$rEOKw79y3I+j%BcJ>q`RpxyZM-yx`Lmc=s|8FHmzUQqP2g#z4sMB_RV z_xE(i@^#n@7<+vg+s5J({TO_Zy*j>7yA~EPny#K14eV>xXx4#p&uRJ9D& ziC7cgu?#i6YLGCMet0|LvUa$Ve6k*a+SO9w!O!u&1h7qlh@LkfQ3cE}CcMFE8T>LbSNK3KV zWfh?#y^v8Vd7;g-;N8RVA>H{fdb4Y%U{#hmQr5?FY!xvn{pa<-R;biRNZl6Y%_?a% zc<$S_;W9zBzdhh5rBKkQeDQ6&zkKkan2B?svjdym7OWL<`tnd0>+k?zJCt!;O4p=cw?8_#uf4{~*vXWN4P<#z=S)eR`G z-%ZlKyoI`k5S+&A`r$Rh9zzTi6}AgxDvVVS9+(&|G+b5k{9uci*o-W!s6*G9iC~@{ zVGSH(pD$ly=Q{R~0H4Pu2lc;eLr(0vB8T9X(VXvtJ^7_=*fi%QkD~NV>`P_Bh zow=nm(g%l8L{E923MK$TtC3t8hGoBetjff2dMj;wKsLWC8HszVZ{k`;zFoR?ym9x4G{#Y(6BL6I z8F;gN2O|Bzna226_*dCP~{}Q z+Y>j#+^Z%mCr_eLod*WXbg<1<-?m89S`$SUDk&i%eQD*w?&7nnI>xOAh*3d2exHa$1MvF7|;fvz&_Vgg)BId$9)%@#)RP2ePx7_V8R zN|*5!`CfYn5ViT5pu4G=4VBr;y^L0ufE#%~%H5}&MF?IbXW7W6M;ao(b{DKCLWfj6 z3yz~*QmgS2d>uh zPqc9s8IX?_$5w3q-AyVi1jHmXRZ*QzNnk{fX$~qAiS8l01IT3`B ztzU&mITp3c)R|?5c6mtRLFv1{G7MJQaQvv6!Q`||Y}ohEubc$r!CW@M$#(7gTT$b= zwq5zm$&9QYW@z>J4efaYk<5d@yUo#cC57ek_b&4ZM1!ui=LNxS`wzp14o!!X8#W!P zD4@oO--oOA{X#{+i_eK6=JcB+mxnvr8@#>jjZ7inX%F+nd-4`Po2yq=#d*N+YUM#xD_K^j!mC-?z`_dR0AKj$bfE*RERkpGx6=ewQEn!%8Jo z?HcIX?Bk^lp(`=!!DhV% z@4ne#g*G<4Hx@2yB=OPb_fpZoSjVkQ2|ts?Z4Xg{l0{czAnWAUiKbOO`!&9=w~f52 z!CIU74$}*7*75S*>tDuKj?E6GSFX3MXWOl%r~;c$H?L&9bI83qE329zyXp1w2ZZ)d z;3-<^!j^L*0eUml=dVE?kahtwj)S%lza8!ASEL%-n6AVetU054F$g@5pmL)tf9;Zu zAIze4(v`@7$cd))a%FPx=X~2OZYWxMC5d(7Lpo|FThh;Jd`LWWUuj>i#%1~94i8Jr zBkX44diw_3qxdqs+)&ohh_P`?yT4cnm}ifbtErO-X%%7SCIhdEfzn;3KlNPbj|wFE zN!aduQpVx$Vv#!|ZTXrBR^_7zChDkMm}lUz1J6;z^8JVCExJR7k3@tzR)ui^b|Z`^ z56lBp_qcUt)xb5O{fs5$$=~D-mutdS*-!EBN1O#R@S!7*OCO$(#FY_l6xQ2WFk`g< zMKh$gyK=*-3Xf7Vynh%pmF+2BxNdv}xY?;QM*jLD4gcURp-`iftKKp*4wU{@+|yFK zX{)=dLlj0WwAqNt1c8f&9cDVnNIjJcJOKW7^!B>#mNpfuuJ9_x#bncj9Vh6m5GPo9^~-^p^db9YOT zzrYiItReXe-D#<9MSB)xJ+CZ?j5L6aMJRwyQ<5V)hvr=Ib}$#*moa}Q2&csc(E>SEM9iZ#31(WQ~59AR$xEa4{f1u4cgFXWK==Q0{_vMS7_bl0ji zeYiDsAkD%>AWF(OhA99@usWq<@?_!*q_ zXs+stqeDMH%l~Yc;Y!V1gY?6yvwI$> zJ>^L?1O2v-L5}wBmn)nXu|iwxV(6A9w<8=9`mOZKtituDGUW z4?zQN=l9NnT=In6M#YgAD1~DP{t2U!p6TzsR)$zW*YwlVWF+Da>;B|dV*jIKt$d1K zl$XqX&wBu8$nj&?6QCROv`dDMbrX$-TRQy5Bh6Bq6b-!cSqW;HK`)EDw}c#TStc4! zWB_caZ}`*`+JoP4gMqzaI;D9DY&EEJZNw~?Tm4j=;%Z; zkn~z+DXN@_K%O<4f453|v7Tplup_MEcdi$O-x;+OaGBsmgIZQPgj+huY{TgnemuT` z_Nmv!A(z%zmL{y0G0w5A_D7P*CpOS23c-(&!v-r4t8Ot!*653xe{|fF3&3;v$y-(u z$&$+Z3iI++f)Sw)OIbf#?$b81-$CU~Q(K}!R-6=w!;vn5tWk=3EBq)zc&8=yv4q== zxk0m|VBe+QRkzUIJH6E`B&4O+aCZoM_x9k){g& z#t5LY*JY_#be`rcHGkjiHZe;x1PPaGu8Y)3Avj;X+M&sL2 z`sXNMl_js!iW!b7w4@01%(;Epu-Q2wyffbGd#$~zzy$%=d4n4*e0tw2Y|K};ca@i1 zPE98OLhkh8Y6i?duJBeIK2yLLxE@G*;6PYn;Aq72=STjUDU77l>(Y%i~}W5`bix+U+Qfy+PlB(^_^o@G=mf{ zO@+Ve;_RtiXMy^NhoyYqq4NNiGNiP4@wYWmtkL7AXYNCHy%C?ArQ-6f*ro#THNnCy zkw(q_lT5~&G2$_u#h1Mg@CnbXToWfRImc|`>XvK(6NcfE!6zGThUb7;r&y|JhXSpU7%KOq2hw+k=!LG^# zc}P`;H^!=YB3#h$X|)Yb#(t;wnRdV3XtjU(R7Jj7yC4!Er|G|Om#=vJ%^Bz8ah+RI zONHy7NSf(2hUQYW#%|=K$f}ko>__b+wexKm?acUfK9~pR8rjh*bjg6n#|(ol7GEk5 ztwW4tzA#^DpF8Gxm2czH&QDePm9^DS3{}Dn`~s?`TstO%jvhOX)&NU%yNl{=*6LON zDj(1j+j|_~8e~}0I)Q9Z+IwHEvUGcy6Qocr327J42tHjen2!@l$Q5iiAl@RU)~-D4 z7wnyLyNEY1p9SpMJaO1dK`4m45wrG|H=7eJ2S&U!{_}Isd{O+t;ysApN zn(iV&=X@L!7yd{wL;F5**fLHLuJ4-Qp^2|LWHek0H_kfM6q5N)T zzDT-R`nzhG5ZD|X`Map+#`emsYQ8$7SD5=&YwMlmH^0Hq&y==$D z*2&~G?V`=}4ULQy!!_FN{5RlbB@fp&sMzY=tRjC|GCLo)X z0B77X2NSH9ns3##p9NObSR|^tH%jQ_zByZHNqsLf2?4E|OqKL#e1UvlIJ-;tTtX!? zkUJY2Z)`wmJ0pB}hd-Fz=X=U#9qH{w;}O*cDC_CZ17VlL4QPUkW0%n#C!!zO43%h( zzoF0V_TIoleEN59-Rsvn!KcpWtiEH!zA_!Zd5pDaj(I^0MW-9>_)wiTR@^{tP)}oF z*QDS%sHQzNFfg?aD9g_QNw!~pA5!<$|~kgY>r+;OdcI2X3t z<?S9W3vCBkviC*8mWs$^Lx)9_5eiiKoNU zrzND>?v~1`gQ}jl0xvyDQxbPM=UQuI(wC4f9xlw`nN^1rfrpib|Ni0f{b`J=JOY<} z@}fj5Ki>bd575xEqP;+WM+&dRkZ7n6noN-8jmpt{WPNaTF1?|GjUBUAh`xC~QCr&J z?%+!Za7Gp;2Ec2h+c%lcg)^H$wZs~0kCI#5e;WUBkhQP6>w8dZ56kEn*0iq`h+FO) zjvK}w@&w*{q;m~uqC=a%v%rnzq^q4WF?tfbGIX(w*;eCcWXpoZ*wr-)8JlGM z67_8-S&6|%g9dW7mzv`>Ghsyf+~wZAZBY7Hg16gopH?h{zcPW&Y8rTM@G?mJve@d~@Kpfh3YJ^5a~fc9 zJj;JDpH~6gDf;YA#8XLG)5Tw%HW8x@OlJk6U*HePwAwGqOx2DEH?!KO?v@#~zgxy= zvLPwj7eap6OVi!2YN*oVFHp>urjGX~1xrkN^h2olvMH~Fj0S-yrR}~-ri8)qj{9>Q z2B~wh#_}kH2$}lVLHQw3=K2854-Aa?iRP=`GQkt|HnJ|iTl~}~E_ul>BT4H$A~C09 z0ErdpBw>z$XFlIGZ}#zp*_JWVnzyoO)5Gd9uhWJAgs!y}x$+D}r|AWv)*in*z_`oM zDR`f2GrujKq2tOflS3YoTMwo3J=wb71Gl-!U;DFog1xGDoOSVhRM+mz`mG)Jv^hv! z8PcZ}c_3&c5*_U8D0SFXJVZ49d@()S#++PnzCoILj_tTDdh>yYC{9p&=g(v#Sf98O zGx3ukzw~;4-E^mryGs{N z30qTaWD0%#G*Y6?*!vVw$|o_%BKm4R>3jNLh8moV zadla^@I~%x(n&;Ghb(M7TXOix^{2fFjQAH8xqFBH==_EFflV6(`L(_-{2sE^Se<51 zjimW`U)<$!`YZZ0;|51+C6SS_r;Pdd8pX-pnYo>!yDlXWTlYlMtyc~G2I!oT8 z>Zs~fT#V#Hq-p`%WplI+RveXO+Cq2QH#IgIY7tc7?*UiFYCpWl{QAnHs;6mQ{EUz4 z9Lk%bG+@gja^nOQm!z$+xmm`gseYdKE}*b=9E!?*4$dxe{~z4WUtNIXW!0 z1$77nyG6z7mB(rH{SSk`0CeL1bh_Q5K4j>7!R7($THb>Z0#xuu0P4pQd`K-w!0lGV zmiBx}1{wZc0&h{H^&?6M{+EntJO5r^i%njNL2b}&g5G*aSHuy2dz@Q;S@-wcvNqaZ zJnGVlgolFz0w-pS_;*$O%=&K1&c*q=OHiD<|Boif#A(CJ*Ka0_si(l-wm&IDHWnMP zV0V>aDXI}RJsYwm8QdZPOf7FFW>jO3f+u^V;GOkNJ6pPzMQKvQae1BcB6?OSybW2q zK<8@!+=kF1pL~}<%rL1;($8o^oN?|eZsA_atnQN&fzP9CY0YqBAa~JBA+L5+BPzIt zGB4D{v|qfre9Q|R_kr(u&s%v>Ko&|ZSzDSt7bD4q{*yiyUGL1C{ih~!@kh&;2(L`v zRKWs}6hGc?%xmO*d$sjGtVXD_TIb-UhH%t*mLUSASv*If{yT(SXw@w=lI?Y>xe!1* z{Ujvh^;sZEo)(uSswu-zj%8vg@*HW_f4-i)!0#A{6pNbIM_U?dMR?Mw59oM^UBx*7NSfW3sTJt0rZLslm zMZQDor&%)LG#6#!QH9mP>M_DbPRXRS**Un!)jkzfuBAJ+CIeyM#Z)Gwk8WCAp^UMT zT@Yh0DOawut*!96HOGSmcwL}o5f8+;V$0%Q!i8)R6P$&9e^+CZ=`+RG0HbO)VtmP2 ze%3GO$D^d}9`$5O59VJ^fX2A>Wnf zPaS|A_ln-W%xCD9cgbN+Jnd@Y{o#*iYCa#ez^LOltJ7;!Yf=AlxaWp6V{W<7>h1wq z*E>XlL`~jvoG;J#1^<{%)E2w$kaE{;(@pc%aDm;3Aw$Irptm4S_t={Aic?dsX~^DT zrD>?_*Brf$96xM#E**btvgqqAQSwC@V-iEsN4)Ot=ZeL$HP_ikvqlBRxeHz9n7Pb~ zZZBjuH<-#aTI|2qE$`t2D#NaFewvRg5+qlz4Yg*^gHf zRNlYqo0&njnxYH3>gath*opqCrV?~(NHM?FbFKYQV~;0LY={>-eX65B9Lwl=WJ+9i z5pRm3{#MZD1f(nRp=2>|wn&!M#?sth+L>n;&O!cidasGCI#GowWaphRl8G=Owio^q z^g@y>+rIRaD+auM)!2RR#`lKgisK-e`gF~#-MVvV5YEZcw*^HvC^IK;n|4ksUXTxe zFqoPxi2*b1mW&E6QgbhO(Neq#Wnwh-Gb8oYC5{M!YYM8VXtZpag}v9FZjemD^H)Jo&<36I)9b_&E=gj^ewNm(V7@ zTWlmV7uI`z@Um~T$8$F66-sr^`M4Gkk#428Lrx9nLCD}a+b6pO13&E*^qY_S6P$O_FtZsiH z7{Ade%;z_x^T?g2zv_yAoqH=Za`xFk#1TcU7#z`xnlOe3`(YQeF6QX`4mF~$z|~;= zoo`d*Z(Yg##W&wlbOQSP5bJ8z)OsQ(nVKdyV;c`W|HsvgktY{LIcSZ9iV2=;AYDgp0{s?K(2}S4>8| zOEH;xM@&Yo3#{5@CH}|OjInOUkTUZwI{@!v)}nuN%eZ#ja*o4c^^qdlW1m}PSWV|v zt;5@34!c7ibf9Y>61!vb=jKo6>4t9`&xwQRvSS#&-YmX=Z7F-^(oW@G2SXl+}l$0kEDOfA;KBOP?L!9*7TXZt)dyY zm8r%SndS=1fgY|2hrr3iS#XEwGUO+y?&4 zY04_5?7Otl529}TZfw|^Zm;_mfrIA!|34Ze&jV>#SYyMtO^(R>RCp5%Hc-{@AxA-{ zaq}6GrTCfSOMC?0CK2|8T`}1pD7yqj<5zAPOrRVNbw#pdU{BiL-~jE_pJs4dL2RbO z#KN%t;YRyn5~Pje2@i!|FbU-oXtgHGeCfj9&YgPuP3;ws=yOwK5$b z&F((ZEjjNRnj%{!X=Lw2vYku(9VTCd@HNeSdGx=8;uH`>nO-)tr1Brn$^N>vlwywZJaJzx?1uPpBn6`gBGf7pXTI*2N@F756E=TV8v5f+ zbij%pIgPI^x46qDXtg=+P>(Jb6R`?C@NSpNkpQ7Tn@c99?#=2hTD*`?{o7*tjO)Uc zwvxkg*dx>`_GD|~)Gpd@6Gmv*@b>c|ea&$?sM^oj7Jqi>`pQ8SeB9mt@ds|T)Kg|L zm+OR}!(ZMU3-f5C5IUfAnQ=klE>5IVYAbwW($J}AeRcSa?;8h%JAp6X7s37{*r*c_ z;6HURhm*x~=4c-TW#LcO<@52Jm7;XtX!dqOgE5CHFKP^`XGZEr@6tIxC)J^cqeD{O zA&R>jT3jk5pEIvxk6m8xJ7Y~M;L)O1t<*(h8Y?5ndD3dIWV2~@Pu1*E+z5_I(5N)o m==Ic-r0I^OzU2nR(g*_a!!qr|ru8J3_lEp4sEs%;&;J8JRpnLy literal 9811 zcmZ9yc{mhY8$X^DZ3q<-Q+YhL&=9hXBufepvSp2sY{|aOSSnOHpAFvhFQP8@9+2h=X?IUKj*&Bxvq1a^SSTqb8bNTsbl{eyjw-F0VQkU z1XQE4tYv5jX?xnoX>urEqs2=b`IocSbFr{fXPW2oFDbM>s*vqpz5XCyUbn?|0Ty0i zmTmBso<0ary>4vt*=*D2e(MSv*lD2YF`E$1d*+tBrbjO=B{@i3Eo>Jr12zDu+J!rp zm9_TK`SoIsx+wdg+|c4Ry|p-Ov|Icx?|-(ZTnpYpw8S9K_v_++iBXstVFTU9J`X~) z;Ph`5x)q49$P4px&4`LGm`5KK_miKtHj@Jv9f)6Da$-$iQo2r5tJ$q~7es5-5)~bY znvb`Oa*x?nV!9&LDU~fQ#IK7cL*4r#rGJ9yL94jt4>c}13(q z`{l!j)qdRnY(@9~P7-b&(YJ7uixZ7i82dRzXYZPToZ3c6)?DNw=l>72#GUNU#gvcNmD?BV_pRyI=#hA$>0+*B zoBA5(CyuMW?SNaU$-xwbA*u=i5l8}Yl@*B~P(l}W$z22b|9q9RsTrbO@EfMM$WK&E z+%|Rn8;@F!)t3%FU^}l`z?{QllU3B8Ge+~lEAX13Kx&6UIk;+DN0%MB)6wfy5O$Os zU$^>FsnU2Ve%xEA!1Ap0%qJVB^?>#Mvje`zw2wj#Z2UW>hV}BvgLKY~4!9x#NvvG+&11`b2 z1L`ZZg*ewom_!*cWFqh$%B#C5q(l%=?0_@Dtvt4YSGet;ivGb;qC>uiXGM{VwrAdL z{QQNS!f`Za5KpB|^w(UGGx3+HExp)EHW)Re+`bK>hue3&fit?RNsT2I1%CsxR?Isz zNfDcz*nhO!9*@+iOJg$zW>$cLd6h<(4xCq~R5&*Vx28@1*by@tFG1B+H^*BcES}ES z+i(O5MqS~$jF7Rw9_D0aWJcSFcW%+&E?C&Y6>H3Ur;0_17lb#9C_mJ~M(mbK=!jR^ zZf!gAG!Z|!(Wrsa`byoWe%^R#*|C}ei(O1T(%70>i&2fz^9qvLjB)cP)ed>-eRFB+ z#&nhP#keg5FSZ-BOBpzgJ07)Suh|Ft~d;hUf0p|L! zBWv%aPye}nxX~o$QQ3kopg3Hz>LEf@dZpAaS(?WkvAEW`5wW;MN1%GD9JNvwDZsXn z!LW@8$UDM79ORKg@`n+jX<9SW(E(l{4D}uJs{Kf_1%t@Q2BGGU$SPb_DbZ#_>L4YA z?!Ty~O=)#lY?Gd0YvfXtGaDA-1&_{?G^hLa;iL?v2*quuNtn(Qpes#0FlT3l>wGQL z^>sDSSyZ40ntZA2I}|F&s%ilTUW53sS?wD)H#^rJ3a#F04*jPCoQ+#A*hkXx`q$ln zR%3l`cAKfWI4g33FU!rf!_umN->}@i_TRVZtHOR5YT=F}R*`6PpuT#M89ZF`-dtbl z*^U4a^s*mf)rq&OVBK34?8Jj688EI7NTqX8yJ9p1?M+9Johr5(OAE2rk|$qc9By#q}(x@CW6v|>5QB;k5ASB!l{g@r(Y=S%G0;{f=F8W zm3klYL<7y))mHyT9!xflXc3VS2qE%@VU?(Pe6TTzK7H1*4uM;P={W3c8P6b=LQ2!= zO5WRFk5XCd4jA6)qieM@`V3d_gkCDA+vPGr0wuP$mN)TcTd~3p)mFC1e5YW$4f#en zvO(}-^e?)ubWQL{+HV%sAR&G6&mHDPH7YC;kyy@oqRt{^*9M{V7&U;E-vTo1Lfyx8 z>%HVyK~4rbuFIrXHYs;#0Q*hu?u}J+6CBz)E}NV4(sz+<%J=X{+CP_5#y#lR2IivH zSziZ|wx#rRmYLn7p^)#CY1m$VGm>VK{f2>t<83fQ3Aej zPSPcv>sBK^=uPLk`23@=-~%Hacgm^_yW5dVRiF~Zav}!5m8iD6K8u_1v-{#3t`OdB zS<#gjp=(X+5Ye~Kp>2J(aWm?boh-#QKewHaoPAxVu=a;YwaHVy^O-SSz2+T&z(mwN z$=TmapT<}JZ1-^E^N?yjUJ27(R&#HSca_*Up^`yT0#m*bi~z@y3nSgjaOv=LYOHyO z6_9eN1J-g;dJzBcOC9o&BCkMkB;uQB znCSFQs!e}sZ2#=LqV+PP3J_E=X-JpoT`3sFT7Kt~@>&@XNP)1>1i!haRtJGoxIv z9|1<@MG<%X9I~`gE@jFuVclhuSpG%9hPx_|Jhv2fm&KPVDRlupSCx50 zsmst{Z=RgH!~9yf9C=5@jPl9nNKz|Kf>;Rv!CU*H6{Ma_; z>X7O0DCyfWd#wZEqY7&};LOPios>8Md}%@PKg-x9p+RG^%Cz?Fvk?zJ+sR(i-=sO} zdIVVyn18`vBjS&9E=J#c7cnW#dZQmaeQ;;ni*)_LN_Qwmztxts`$Y}#>vZRwq5lz<1fGN<_vU)>CjTl0{mD2web5|?K zZW9)0?!>Rz13C?gxZ6>ejPL}2YNoP#gd`^45_G@1qzv;;|KcIW_;OE_qfI_4vnCrWe4fyx-&RJ=xB)=bVuFw#gX^S!YWd24YQ}SX+TDXDFoWH zARf1dDv5tcT4?pO=0CK#Rt0*AaRubLy#k^~)+C4dPYpqz9C=B<7$_^v17Gw7xvGrQ zx@-onC}qBno#j3J%X1j35c1WsVoRJ-eimAwlJUgNuXk8JWL$Hm&C{(6)mjkb0NXuJ z_8K%j$jY~&pGE@fJ_Qa+vpV(7zj^ksy)0&6;7=g>YhPQ8etbVwpdE79&kxhljo@#^ zzGL=Z_K{kVF(ja&BP;ONUs@H$ke3j?U{lVD9XjP1sIbj#?2)+C)RuTG|Savo-mOwH3UjlxfIGM zv=P&VFEYnaJ~`k)ApO%S5ctND2nwCuf`v2WnMo7d{$5AxWu~Zk;d3_yp;$scSZ}%H1rqW=vSYFiXshN88-tqo6h(w^RXOCHfx!;%Y*) zmi6*n(ga-hSBS|+dla*ALXh{pp~-oNad<}b$6EqB`jMT|FG6Tv*Rgi8qXYSwOw$~- zvkS5J)Be@BEn7%hFtxl|m;DQiSiDvN&bs1dA)w^7W8#9Pfz%cx-ys)gh3gn5d#&Ox``(<9bFTK?~#lXtREQ^_ru>?7b?R zl;aNCew_(NCO_oHDyHIlG~N+x%>e)~K7YcEv_k+&@ATzQXcq3r4hqsv)U%u>r9Ylg z(W-*IUvERVK2)Of)&J+qMEm#Zv6=6N-%m1gqIeiau3z}egmN0Be@pC1fOB_SB4mrt za3UY2C62|9C6cP-)&GvUJ~q29g*nRHmfrH&ajzDkDR0QG?uZviJiAe zVEFA{rc#f#x0VRJ=FIw4#v7Kj<5sH8w6wz}Q#JQ+00jC4qUZEW#sz7(&&6inN0HQ0 z9vmoA?84~<9pVs_E4~1%lpW*++7-n9?EguwVzq+AmeLa@(|#W~K0BVcPYE0r@;_@pl#chs1A1 z*-pi~bzqaJV|;lpp)UbEUY8?c2k(2HxfIX+td{jbXC9a03sc*DJE5w*tRy~eOq{E| z?_Ky?mWZ6Nj+u61NM8Bm{-xxXRd;E|ziB}IIz1Wnnkn3@)`smI|Dw~fD|DLqdDiHS zd7!F-w;J+6{i+Ho)JY+B3>cHch+(k>Py@v$^uJ$N7mBz)SC!cBjy%V7-EsLVVrlSEck3*jUdnw`HuZYZws6^9@OVC{gs8&Se>*Z)2tDXI6; zGWyL|_q3dhoZ^)M6Q~E5JjDpj`4&n(HlJwfA>(R}vf81gPK;g#ZdIDEe0O+lk;;Ku zoOyGB+6YZZPfR%{4>W#Cd=>pKUi&i@xpTUroAr%rvYN=B(1S)ez+(S7fRuB`g2&Rv zW4}Wh!gs?P_)Bv=Y6dQQD$^9>&kUL4Q1HTn-{Y}k;XB!xid&$R(wv{Rj|xp`Eawixb8MH#KJFr`_wD^ zyeFZVvDcJ(u6I5lCQpA_OFnDt6d?OsD^@~B1N_(OrBMDjy?J!xNx?J$$Nmt{?Uv_( z@X~bA;7Q_|0}R=1i4toa(U7fd9f$15tUnxnv?KIQ|1cp~Vt!ffj#eh!GhW72Ui9$P z?TcWU0E^$T78`9Kx94r!G`^gee;fB6D(eV?j*+_LjN z+pyZU-ZVcnxSM-#+cPCWKcS(I@QK1lF^Y(6mm+6W%aMCjh#Wv58l1tt>|LK}pA&Wc z@XW5=tp2(TTB6X}ljN}(HVQpM`rIpm20?SSg~^Z}S72x%3!TLqg9u<(O=m6!RKu007gE{;06lXkzNu$r8p zrwT1mT0a?a|BCN7b=Hq>^hScXu4GZRcyl9`CobS&v~;6|1n(;YZQ5!%rpY?c93spi zdshSrho`96*cJoQj^5yp@eOe9k%dCT%DFoaW^xFU*WsS_X@(7#b~Bflnf_$DA^m}Q z2Cqv-)?)^3ff=Wdizr9K;`G6=$*2 z>bgP7dU>q|uE=gal=;45CGRDY>d+OD-i(j}51B(Q&tBOzYQpS4ot<7!APcYmGJ-2j zOY1$i>p$~^`Gk-Xlw#=+iW|b@9seD%b7T7X@Z=MeWGD~ezP1~$glMSZYS`iWocoSn zVzl&LtR!2Ebp6oP0W5c_BxIN;t_r=f?3K|H=e7|#07_@oPw$_|xKc9GtcGvDLK!bd z3n!z7%%!B~q_7MakiQmlJ26l4TVw!s$NGv5SU0PE=vmzacS7(~3_o@V_}A>5;r`ND z@Of5zPhDwv8^PD&WQl!_e`dV2RMM>2!MUfFB_jN((GkWQ6ke^XAmvqFV9eAck&~pX zoQqlE9)Iu`{y4&?oSQVTJ{2w-aY;|cO&R_P%0*5;v%297Q(bVJ^+ZSg67sOR_O}3+ zW~T6q+uNTGYJT305=9&@aj-aI13m4>PiDI|2wkxq-LjS)?TAu(g|Gm7DE-k!u1?=@ z{PR|(f08iqV(#*@mScd2@=t-BSJf%n7bq^a7!k7eZ_H_8hX{G(-NMcd5sv>3UY4j` zE+6S~BqdKvE4Gp7H+7l&=v;}9Qj`ExeIWkx#TlP)X7W->%`GyGbB^I(4Q*TI;lZs@V>H)6M<7ab_MI4<2+Y zI`sEJnoYMyQiE!}me-n>DY{8QmSaYi**#&smb5o4)$qdlO}ljcHBYnenWiwpU-)&D zP!3Y;!mYPakHbhy8~%^Ewfl}Fpu=;_p}+y6y^)q=V8}Ek9#IZ;55>E#@uiEC##EhV z)^Z4m;&x!V*dh=UXCYz4uW3j==If9;A_Kny;RCu4_hJMGo5c4qEyJB?A=EcLRZErx zOXlwz&n}bR-gck`jXrpF?i84hn`G+I8Hfq)s8giBmzRuLaq?!p*txyxWZ~Zqa z$brf;!ngkX05J9A5@+>f+RVEamQso@u~ddwOSncgS|gmb-5+XBaM@8tgi`>7EZ(1-JER**SkejHJFdQYf$B3r*k62DZbC>|Y`){TkR$Daj#g6?E+gpCA7uf+Cw(+i?*dUoHY9Fx8f z%MZkOm<*FLIlL6!CLw7RdF+-aLUjB{Kz`EJz~$GJ zA;x=#p#)?Ad9HJXuKh{RH{5ZfSK`rp8M;D5OvKOb$Re{kMZXG;K3rGdg=^GuK^K6; z@!+sQjG&pj5+`odny0AdOQiUf4i{Fa80+IN~kwZxn0w&RDU+tJlyzCT4kT- z;@~WU`#Z{>TIZQVPKozrHSb>$23Co|xLo1e_ANrcKbF*EJ*@{3tn}%}rVsSg@r!aQV zpijO1;jOO+Jbm4cB_7EkeejjxYpmsmfx}S8_aRv+>L;^(-INuvoyGQ>f>Ic539r)` zJOzBNEtLm;UDXGoVK_A^o@wBDBb7(??^8u8&!i6CQjgo8Gf3hNzF1CA6Ba@$%X5W( z2zwVEmF^XCgmZ-of}GNL>GIxpr-Zvel^FthrOo<7yJg+m{sj}rt4P`Tw4*yN43@I> zHtp~U|0~NY%v-(+Y4mjJy^I%wKRc(HZ?%O?7Q=U_)q4d-LzsPikA4ke+iqFy5~;~t z;<#oS_lh(mG9`<7Rd8V-Ub0olDuuE&zEcc~F(;FjSukNJDM2big zM*Ggq4SR2LOxqq4w0qa|MGZaEf#Wa zwrKX-F6$4;YB4Tivj`%@AShrrbdcgGYUm%<>9}&u_UHre?)RnF=rF!`ZQIAT%dKPy zZ!+L$ynO24+cSuNVx;%qG$k*MqZ2}=U$`A9>mFDv+a1zQMJR#fe~)Y`+fur@y@JQm zy{t?l&+1=&7VGk5=%hYW8?b+UnpRCP9i4`Oepd?>2(qeekG8HJrmA+p^+eB8Ura;W z$ClhvMBW$M%$Iub+hW8o(lidH)iGSiDar6Z`J6AD*mlh^_FBQv8Y`9Hi!wpg_nQ)XfJ{6*{2r0kCu?%$9DsXwP~6k>i_%k`5WV$`ZcInS z*RT@pL8pji3opE|zd6g-RYLVL_xwDz4w%t^TKKC)X)?Bw}Uk3dqwD?_jk)!?UvnhU1nA^s{U#-WS zj23%MxP38l#Y%^2cIHDEb>1GWgFE&4^}p)5siNt=g}a%No8~uAl?ONz(9wzUSI+Gg@?h5h~#}E4+{DqO93^bU7f*eP{mr64LLsh8fcJ8_>n6vew}dPswNOG z%h8%Q`FkoS^tAhhC~sGZwdQfiEYI%3qtCj{5i9WPfhrRpu6Ng&QB)^z#*91J;A%5~ zORNJxWy17rm{m5o-ORw|b8PxSxiRCcF(ciW@xYidWXw=9W;9%9M`eQ^e?~q2486Ef zcX31R?ZNJK_S0?e&{UDtTDQrGXEiMkBY;5=vj zJgnJu_yeFC`XMt?GXszrDliC-JA7)J&L0j920!y%;q-4Q>*hyjZTsU_3Y?LPhfnw6 zzJo2~MI>rhVGWP+!rCI}1nyr)h?za7rC#Ljy-7|#q>=hS3&C8npnw^bk$}z(TuN{Z zJ-YxOe^41|5ITt{A#R&aXD;&u5nAd(=ng@pBpj0#M+HBVe7}}8y!Of*#q=0e}Vf$j?SW{c}|82aWCa(g5v+7GZ;!Hf|wz893 z%vhCW#JH3mZr)tGRQtC|kX5jAupd2>yGDx2uD#V&k&J?Jn$1{x#IM|Y)sLGA^TpG1MR+1rinY+Rb zV7T>e1V^wZtys_*$*efCLouF1$-SI)h5cI!z*_wI=KhTOBr{!^A55VWq_K(o)ANu` z_d3@S0wN;*2it#-nqF7y3YlMFeyeleJ(M0}kQRQ^qVTMb?uS&i)TBT#X9CLUVW)H9 zzecR|>CkN!)j<-!NZYUB?s;T$eIa^H{(H}R3sUS2_KOTSRsxWNX(jGv)dwxDxV23P zp50KzwdgH?YWz_{pd6E4Ks*T%8O-Sh_#KW~f$4}69KLhw{XS%8x44O5hU9!fS!f-_ zPA-Ja5<|9GX)UA`f0x|eqLU9#Mb88W~KlBa3hHhaN zi;TB@mTxK|WE7R||7_zh=h+sw_KdV}%|Xbz+J%2L*=mhoZZSL4I5NK^Xr3N5xy?qh z$0dV;gZ;oSW5=7Mm<6^|`g2{ev9{Nylr>8NSX*1%!!(nx?Y^~dT{DeO5}AOsnTV}S ztu_PN6ufTnH;vQ3>Z(NsIK$sJ^f{x8q`4mq%9WS6D(E&A%6o&(d=@C~>``Z9D+BE} z0q+YC$Cn1n|IKbQX#ZnU7m3&x7Vr@j^`py#87%?#+qv zr6dH+qNd}u69*NPM~vifv3nzQu`~CFGJBDbHKcy5MTvzaA@s~#z25?9z!B|#gnjQ? zRTsDX$+STFkMJB|-*6M=&kjl_JnU@&Zm&c}eD`#45<9;^R zIeyNBQt*?FsA&ENA(M4`1BBSvA6gUw=bqP)XmgjByX^{Ry&6e}9Fy~pylErGbL{^B D{IS{@ diff --git a/inst/sbc/make_reference_rankhist.R b/inst/sbc/make_reference_rankhist.R index 91c46d1..b4c57ab 100755 --- a/inst/sbc/make_reference_rankhist.R +++ b/inst/sbc/make_reference_rankhist.R @@ -5,14 +5,19 @@ start_time <- Sys.time() here::i_am("inst/sbc/make_reference_rankhist.R") library(here) +rbest_lib_dir <- here("build", "installed") +rbest_source_dir <- here() + +## ensure that the current dev version of OncoBayes2 is build and +## loaded. +## code to use if one wants the more conservative load_OB2_dev_install routine to be used setwd(here()) -system("make binary") +system("make dev-install", ignore.stdout=TRUE, ignore.stderr=TRUE) setwd(here("inst", "sbc")) pkg <- c("assertthat", "rstan", "mvtnorm", "checkmate", "Formula", "abind", "dplyr", "tidyr", "here", "bayesplot") sapply(pkg, require, character.only=TRUE) - library(clustermq) library(data.table) library(knitr) @@ -20,6 +25,10 @@ sbc_tools <- new.env() source(here("inst", "sbc", "sbc_tools.R"), local=sbc_tools) set.seed(453453) +sbc_tools$rbest_lib_dir <- rbest_lib_dir +sbc_tools$rbest_source_dir <- rbest_source_dir +sbc_tools$load_rbest_dev() + scheduler <- getOption("clustermq.scheduler") if(is.null(scheduler)) { diff --git a/inst/sbc/sbc_tools.R b/inst/sbc/sbc_tools.R index 453f6e6..c2b3cde 100644 --- a/inst/sbc/sbc_tools.R +++ b/inst/sbc/sbc_tools.R @@ -2,13 +2,13 @@ #' Utilities for SBC validation #' load_rbest_dev <- function() { - here::i_am("inst/sbc/sbc_tools.R") + if(rbest_lib_dir != .libPaths()[1]) { + cat("Unloading RBesT and setting libPaths to dev RBesT install.\n") + unloadNamespace("RBesT") + .libPaths(c(rbest_lib_dir, .libPaths())) + } if(!("RBesT" %in% .packages())) { - cat("RBesT not yet loaded, bringing up local dev version.\n") - library(devtools) - devtools::load_all(here()) - } else { - cat("RBesT is already loaded.\n") + require(RBesT, lib.loc=rbest_lib_dir) } } @@ -159,58 +159,3 @@ run_sbc_case <- function(job.id, repl, data_scenario, family, mean_mu, sd_mu, sd c(list(job.id=job.id, time.running=runtime["elapsed"]), fit) } -## Submits to batchtools cluster with fault tolerance, i.e. -## resubmitting failed jobs max_num_tries times -auto_submit <- function(jobs, registry, resources=list(), max_num_tries = 10) { - all_unfinished_jobs <- jobs - - num_unfinished_jobs <- nrow(all_unfinished_jobs) - num_all_jobs <- num_unfinished_jobs - remaining_tries <- max_num_tries - all_jobs_finished <- FALSE - while (remaining_tries > 0 && !all_jobs_finished) { - remaining_tries <- remaining_tries - 1 - - message("Submitting jobs at ", Sys.time()) - # Once things run fine let's submit this work to the cluster. - submitJobs(all_unfinished_jobs, resources=resources) - # Wait for results. - waitForJobs() - message("Finished waiting for jobs at ", Sys.time()) - - # Check status: - print(getStatus()) - - # Ensure that all jobs are done - if (nrow(findNotDone()) != 0) { - not_done_jobs <- findNotDone() - print(getErrorMessages(not_done_jobs)) - ##browser() - ##invisible(readline(prompt="Press [enter] to continue")) - - message("Some jobs did not complete. Please check the batchtools registry ", registry$file.dir) - all_unfinished_jobs <- inner_join(not_done_jobs, all_unfinished_jobs) - - if (num_unfinished_jobs == nrow(all_unfinished_jobs) && nrow(all_unfinished_jobs) > 0.25 * num_all_jobs) - { - # Unfinished job count did not change -> retrying will probably not help. Abort! - warning("Error: unfinished job count is not decreasing. Aborting job retries.") - remaining_tries <- 0 - } - - if (num_unfinished_jobs == nrow(jobs)) - { - # All jobs errored -> retrying will probably not help. Abort! - warning("Error: all jobs errored. Aborting job retries.") - remaining_tries <- 0 - } - - num_unfinished_jobs <- nrow(all_unfinished_jobs) - message("Trying to resubmit jobs. Remaining tries: ", remaining_tries, " / ", max_num_tries) - } else { - all_jobs_finished <- TRUE - } - } - - invisible(all_jobs_finished) -} diff --git a/tests/testthat/test-EM.R b/tests/testthat/test-EM.R index 3fd0705..c1eae92 100644 --- a/tests/testthat/test-EM.R +++ b/tests/testthat/test-EM.R @@ -1,4 +1,4 @@ -context("EM: Expectation-Maximization") +#context("EM: Expectation-Maximization") ## test that the EM algorithms recover reliably test distributions; ## test criterium is a "sufficiently" small KL divergence @@ -145,23 +145,23 @@ EM_mvn_test <- function(mixTest, seed, Nsim=1e4, verbose=FALSE, ...) { expect_true(all(EMmix1 == EMmix2), info="Result of EM is independent of random seed.") } -test_that("Normal EM fits single component", EM_test(ref$norm_single, 3453563, Nsim, verbose)) -test_that("Normal EM fits heavy-tailed mixture", EM_test(ref$norm_heavy, 9275624, Nsim, verbose)) -test_that("Normal EM fits bi-modal mixture", EM_test(ref$norm_bi, 9345726, Nsim, verbose)) +test_that("Normal EM fits single component", { EM_test(ref$norm_single, 3453563, Nsim, verbose) }) +test_that("Normal EM fits heavy-tailed mixture", { EM_test(ref$norm_heavy, 9275624, Nsim, verbose) }) +test_that("Normal EM fits bi-modal mixture", { EM_test(ref$norm_bi, 9345726, Nsim, verbose) }) -test_that("Multivariate Normal EM fits single component", EM_mvn_test(ref$mvnorm_single, 3453563, Nsim, verbose)) -test_that("Multivariate Normal EM fits heavy-tailed mixture", EM_mvn_test(ref$mvnorm_heavy, 9275624, Nsim, verbose)) -test_that("Multivariate Normal EM fits bi-modal mixture", EM_mvn_test(ref$mvnorm_bi, 9345726, Nsim, verbose)) -test_that("Multivariate Normal EM fits bi-modal mixture 1D", EM_mvn_test(ref$mvnorm_bi_1D, 9345726, Nsim, verbose)) +test_that("Multivariate Normal EM fits single component", { EM_mvn_test(ref$mvnorm_single, 3453563, max(1E4, Nsim), verbose) }) +test_that("Multivariate Normal EM fits heavy-tailed mixture", { EM_mvn_test(ref$mvnorm_heavy, 9275624, max(1E4, Nsim), verbose) }) +test_that("Multivariate Normal EM fits bi-modal mixture", { EM_mvn_test(ref$mvnorm_bi, 9345726, max(1E4, Nsim), verbose) }) +test_that("Multivariate Normal EM fits bi-modal mixture 1D", { EM_mvn_test(ref$mvnorm_bi_1D, 9345726, max(1E4, Nsim), verbose) }) -test_that("Gamma EM fits single component", EM_test(ref$gamma_single, 9345835, Nsim, verbose)) -test_that("Gamma EM fits heavy-tailed mixture", EM_test(ref$gamma_heavy, 5629389, Nsim, verbose)) -test_that("Gamma EM fits bi-modal mixture", EM_test(ref$gamma_bi, 9373515, Nsim, verbose)) +test_that("Gamma EM fits single component", { EM_test(ref$gamma_single, 9345835, Nsim, verbose) }) +test_that("Gamma EM fits heavy-tailed mixture", { EM_test(ref$gamma_heavy, 5629389, Nsim, verbose) }) +test_that("Gamma EM fits bi-modal mixture", { EM_test(ref$gamma_bi, 9373515, Nsim, verbose) }) -test_that("Beta EM fits single component", EM_test(ref$beta_single, 7265355, Nsim, verbose)) -test_that("Beta EM fits single component with mass at boundary", EM_test(ref$beta_single_alt, 7265355, Nsim, verbose, constrain_gt1=FALSE)) -test_that("Beta EM fits heavy-tailed mixture", EM_test(ref$beta_heavy, 2946562, Nsim, verbose)) -test_that("Beta EM fits bi-modal mixture", EM_test(ref$beta_bi, 9460370, Nsim, verbose)) +test_that("Beta EM fits single component", { EM_test(ref$beta_single, 7265355, Nsim, verbose) }) +test_that("Beta EM fits single component with mass at boundary", { EM_test(ref$beta_single_alt, 7265355, Nsim, verbose, constrain_gt1=FALSE) }) +test_that("Beta EM fits heavy-tailed mixture", { EM_test(ref$beta_heavy, 2946562, Nsim, verbose) }) +test_that("Beta EM fits bi-modal mixture", { EM_test(ref$beta_bi, 9460370, Nsim, verbose) }) test_that("Constrained Beta EM respects a>1 & b>1", { unconstrained <- mixbeta(c(0.6, 2.8, 64), c(0.25, 0.5, 0.92), c(0.15, 3, 15)) diff --git a/tests/testthat/test-gMAP.R b/tests/testthat/test-gMAP.R index 512ee92..61f0aba 100644 --- a/tests/testthat/test-gMAP.R +++ b/tests/testthat/test-gMAP.R @@ -1,4 +1,3 @@ -context("gMAP: Generalized Meta Analytic Predictive") ## test gMAP results using SBC and with matching rstanarm models diff --git a/tests/testthat/test-mixdiff.R b/tests/testthat/test-mixdiff.R index fbe7fae..154d05c 100644 --- a/tests/testthat/test-mixdiff.R +++ b/tests/testthat/test-mixdiff.R @@ -1,4 +1,3 @@ -context("mixdiff: Mixture Difference Distribution") ## test that calculations for the cumulative distribution function of ## differences in mixture distributions are correct by comparison with @@ -64,18 +63,18 @@ mixdiff_cmp_norm <- function(case, mixref) { expect_true(all(res_quants < eps)) } -test_that("Difference in beta variates evaluates correctly", mixdiff_cmp(beta)) -test_that("Difference in beta mixture variates evaluates correctly", mixdiff_cmp(betaMix)) +test_that("Difference in beta variates evaluates correctly", { mixdiff_cmp(beta) }) +test_that("Difference in beta mixture variates evaluates correctly", { mixdiff_cmp(betaMix) }) -test_that("Difference in gamma variates evaluates correctly", mixdiff_cmp(gamma)) -test_that("Difference in gamma mixture variates evaluates correctly", mixdiff_cmp(gammaMix)) +test_that("Difference in gamma variates evaluates correctly", { mixdiff_cmp(gamma) }) +test_that("Difference in gamma mixture variates evaluates correctly", { mixdiff_cmp(gammaMix) }) -test_that("Difference in normal variates evaluates correctly", mixdiff_cmp(norm)) -test_that("Difference in normal mixture variates evaluates correctly", mixdiff_cmp(normMix)) +test_that("Difference in normal variates evaluates correctly", { mixdiff_cmp(norm) }) +test_that("Difference in normal mixture variates evaluates correctly", { mixdiff_cmp(normMix) }) ## for the normal we can use exact analytical results -test_that("Difference in normal variates evaluates (analytically) correctly", mixdiff_cmp_norm(norm,norm_ref)) -test_that("Difference in normal mixture variates evaluates (analytically) correctly", mixdiff_cmp_norm(normMix,normMix_ref)) +test_that("Difference in normal variates evaluates (analytically) correctly", { mixdiff_cmp_norm(norm,norm_ref) }) +test_that("Difference in normal mixture variates evaluates (analytically) correctly", { mixdiff_cmp_norm(normMix,normMix_ref) }) ## now test difference distributions on the link-transformed scales ## (the cannonical cases) @@ -86,14 +85,14 @@ apply_link <- function(dists, link) { } ## log-odds -test_that("Difference in beta variates with logit link evaluates correctly", mixdiff_cmp(apply_link(beta, RBesT:::logit_dlink))) -test_that("Difference in beta mixture variates with logit link evaluates correctly", mixdiff_cmp(apply_link(betaMix, RBesT:::logit_dlink))) +test_that("Difference in beta variates with logit link evaluates correctly", { mixdiff_cmp(apply_link(beta, RBesT:::logit_dlink)) }) +test_that("Difference in beta mixture variates with logit link evaluates correctly", { mixdiff_cmp(apply_link(betaMix, RBesT:::logit_dlink)) }) ## relative risk -test_that("Difference in beta variates with log link evaluates correctly", mixdiff_cmp(apply_link(beta, RBesT:::log_dlink))) -test_that("Difference in beta mixture variates with log link evaluates correctly", mixdiff_cmp(apply_link(betaMix, RBesT:::log_dlink))) +test_that("Difference in beta variates with log link evaluates correctly", { mixdiff_cmp(apply_link(beta, RBesT:::log_dlink)) }) +test_that("Difference in beta mixture variates with log link evaluates correctly", { mixdiff_cmp(apply_link(betaMix, RBesT:::log_dlink)) }) ## relative counts -test_that("Difference in gamma variates with log link evaluates correctly", mixdiff_cmp(apply_link(gamma, RBesT:::log_dlink))) -test_that("Difference in gamma mixture variates with log link evaluates correctly", mixdiff_cmp(apply_link(gammaMix, RBesT:::log_dlink))) +test_that("Difference in gamma variates with log link evaluates correctly", { mixdiff_cmp(apply_link(gamma, RBesT:::log_dlink)) }) +test_that("Difference in gamma mixture variates with log link evaluates correctly", { mixdiff_cmp(apply_link(gammaMix, RBesT:::log_dlink)) }) diff --git a/tests/testthat/test-mixdist.R b/tests/testthat/test-mixdist.R index 202e618..84a9cdf 100644 --- a/tests/testthat/test-mixdist.R +++ b/tests/testthat/test-mixdist.R @@ -1,4 +1,3 @@ -context("mixdist: Mixture Distribution") ## various tests around mixture distributions set.seed(234534) @@ -41,14 +40,14 @@ pmix_lower_tail_test <- function(mix, N=Nsamp_quant) { do_test(preddist(mix, n=100)) } -test_that("Cumulative beta distribution function evaluates lower.tail correctly", pmix_lower_tail_test(beta)) -test_that("Cumulative beta mixture distribution function evaluates lower.tail correctly", pmix_lower_tail_test(betaMix)) +test_that("Cumulative beta distribution function evaluates lower.tail correctly", { pmix_lower_tail_test(beta) }) +test_that("Cumulative beta mixture distribution function evaluates lower.tail correctly", { pmix_lower_tail_test(betaMix) }) -test_that("Cumulative normal distribution function evaluates lower.tail correctly", pmix_lower_tail_test(norm)) -test_that("Cumulative normal mixture distribution function evaluates lower.tail correctly", pmix_lower_tail_test(normMix)) +test_that("Cumulative normal distribution function evaluates lower.tail correctly", { pmix_lower_tail_test(norm) }) +test_that("Cumulative normal mixture distribution function evaluates lower.tail correctly", { pmix_lower_tail_test(normMix) }) -test_that("Cumulative gamma distribution function evaluates lower.tail correctly", pmix_lower_tail_test(gamma)) -test_that("Cumulative gamma mixture distribution function evaluates lower.tail correctly", pmix_lower_tail_test(gammaMix)) +test_that("Cumulative gamma distribution function evaluates lower.tail correctly", { pmix_lower_tail_test(gamma) }) +test_that("Cumulative gamma mixture distribution function evaluates lower.tail correctly", { pmix_lower_tail_test(gammaMix) }) ## tests the quantile and distribution function against simulated samples mix_simul_test <- function(mix, eps, qtest, ptest = seq(0.1, 0.9, by=0.1), S=Nsamp_equant) { @@ -63,40 +62,52 @@ mix_simul_test <- function(mix, eps, qtest, ptest = seq(0.1, 0.9, by=0.1), S=Nsa expect_true(all(res_probs < eps)) } -test_that("Beta quantile function is correct", mix_simul_test(beta, eps, c(0.1, 0.9))) -test_that("Beta mixture quantile function is correct", mix_simul_test(betaMix, eps, c(0.1, 0.9))) +test_that("Beta quantile function is correct", { mix_simul_test(beta, eps, c(0.1, 0.9)) }) +test_that("Beta mixture quantile function is correct", { mix_simul_test(betaMix, eps, c(0.1, 0.9)) }) -test_that("Normal quantile function is correct", mix_simul_test(norm, eps, c(-1, 0))) -test_that("Normal mixture quantile function is correct", mix_simul_test(normMix, eps, c(4, 1))) -test_that("Normal mixture with very weak component quantile function is correct", mix_simul_test(normMixWeak, eps, c(4, 1))) +test_that("Normal quantile function is correct", { mix_simul_test(norm, eps, c(-1, 0)) }) +test_that("Normal mixture quantile function is correct", { mix_simul_test(normMix, eps, c(4, 1)) }) +test_that("Normal mixture with very weak component quantile function is correct", { mix_simul_test(normMixWeak, eps, c(4, 1)) }) -test_that("Gamma quantile function is correct", mix_simul_test(gamma, eps, c(2, 7))) -test_that("Gamma mixture quantile function is correct", mix_simul_test(gammaMix, eps, c(2, 7), ptest = seq(0.2, 0.8, by=0.1))) +test_that("Gamma quantile function is correct", { mix_simul_test(gamma, eps, c(2, 7)) }) +test_that("Gamma mixture quantile function is correct", { mix_simul_test(gammaMix, eps, c(2, 7), ptest = seq(0.2, 0.8, by=0.1)) }) ## problematic gamma (triggers internally a fallback to root finding) gammaMix2 <- mixgamma(c(8.949227e-01, 7.051570e-01, 6.125121e-02), c(1.049106e-01, 3.009986e-01, 5.169626e-04), c(1.666667e-04, 1.836051e+04, 1.044005e-02)) -test_that("Singular gamma mixture quantile function is correct", mix_simul_test(gammaMix2, 10*eps, c(1, 1E3), ptest = seq(0.2, 0.8, by=0.1))) +test_that("Singular gamma mixture quantile function is correct", { mix_simul_test(gammaMix2, 10*eps, c(1, 1E3), ptest = seq(0.2, 0.8, by=0.1)) }) consistent_cdf <- function(mix, values) { dens <- dmix(mix, values) cdf <- pmix(mix, values) + lcdf <- pmix(mix, values, log.p=TRUE) expect_true(all(diff(cdf) >= 0)) - expect_numeric(dens, any.missing=FALSE) - expect_numeric(cdf, any.missing=FALSE) + expect_numeric(dens, lower=0, finite=TRUE, any.missing=FALSE) + expect_numeric(cdf, lower=0, upper=1, finite=TRUE, any.missing=FALSE) + expect_numeric(lcdf, upper=0, finite=FALSE, any.missing=FALSE) } -test_that("Beta CDF function is consistent", consistent_cdf(beta, seq(0.1, 0.9, by=0.1))) -test_that("Beta mixture CDF function is consistent", consistent_cdf(betaMix, seq(0.1, 0.9, by=0.1))) +consistent_ccdf <- function(mix, values) { + dens <- dmix(mix, values) + ccdf <- pmix(mix, values, FALSE) + lccdf <- pmix(mix, values, FALSE, TRUE) + expect_true(all(diff(ccdf) <= 0)) + expect_numeric(dens, lower=0, finite=TRUE, any.missing=FALSE) + expect_numeric(ccdf, lower=0, upper=1, finite=TRUE, any.missing=FALSE) + expect_numeric(lccdf, upper=0, finite=FALSE, any.missing=FALSE) +} + +test_that("Beta CDF function is consistent", { consistent_cdf(beta, seq(0.1, 0.9, by=0.1)) }) +test_that("Beta mixture CDF function is consistent", { consistent_cdf(betaMix, seq(0.1, 0.9, by=0.1)) }) -test_that("Normal CDF is consistent", consistent_cdf(norm, seq(-2, 2, by=0.1))) -test_that("Normal mixture CDF is consistent", consistent_cdf(norm, seq(-2, 2, by=0.1))) +test_that("Normal CDF is consistent", { consistent_cdf(norm, seq(-2, 2, by=0.1)) }) +test_that("Normal mixture CDF is consistent", { consistent_cdf(norm, seq(-2, 2, by=0.1)) }) -test_that("Gamma CDF function is consistent", consistent_cdf(gamma, seq(2, 7, by=0.1))) -test_that("Gamma mixture CDF function is consistent", consistent_cdf(gammaMix, seq(2, 7, by=0.1))) +test_that("Gamma CDF function is consistent", { consistent_cdf(gamma, seq(2, 7, by=0.1)) }) +test_that("Gamma mixture CDF function is consistent", { consistent_cdf(gammaMix, seq(2, 7, by=0.1)) }) ## problematic beta which triggers that the cumulative of the @@ -104,8 +115,12 @@ test_that("Gamma mixture CDF function is consistent", consistent_cdf(gammaMix, s ## again once 2.18 is out) ## problematic Beta density -bm <- mixbeta(c(1.0, 298.30333970, 146.75306521)) -test_that("Problematic (1) BetaBinomial CDF function is consistent", consistent_cdf(preddist(bm, n=50), 0:50)) +bm1 <- mixbeta(c(1.0, 298.30333970, 146.75306521)) +test_that("Problematic (1) BetaBinomial CDF function is consistent", { consistent_cdf(preddist(bm1, n=50), 0:50) }) +test_that("Problematic (1) BetaBinomial CCDF function is consistent", { consistent_ccdf(preddist(bm1, n=50), 0:50) }) +bm2 <- mixbeta(c(1.0, 3 + 1/3, 47 + 1/3)) +test_that("Problematic (2) BetaBinomial CDF function is consistent", { consistent_cdf(preddist(bm2, n=50), 0:50) }) +test_that("Problematic (2) BetaBinomial CCDF function is consistent", { consistent_ccdf(preddist(bm2, n=50), 0:50) }) ## tests for the multivariate normal mixture density p <- 4 diff --git a/tests/testthat/test-mixstanvar.R b/tests/testthat/test-mixstanvar.R index b5b463c..813e814 100644 --- a/tests/testthat/test-mixstanvar.R +++ b/tests/testthat/test-mixstanvar.R @@ -1,4 +1,4 @@ -context("mixstanvar: Mixture Distribution brms Adapter") + skip_if_not_installed("brms") @@ -56,9 +56,17 @@ mixstanvar_simul_test <- function(mix, mixstanvar_test <- function(mix, brms_args) { brms_prior_empty <- do.call(brms::brm, c(brms_args, list(seed=1423545, refresh=0, sample_prior="only", stanvars=mixstanvar(prior=mix), empty=TRUE))) - stan_dist <- paste0("mix", gsub("Mix$", "", class(mix)[1]), "_lpdf") + mix_class <- gsub("Mix$", "", class(mix)[1]) + stan_dist_lpdf <- paste0("mix", mix_class, "_lpdf") + stan_dist_lcdf <- paste0("mix", mix_class, "_lcdf") + stan_dist_lccdf <- paste0("mix", mix_class, "_lccdf") + stan_dist_cdf <- paste0("mix", mix_class, "_cdf") ## look for the declared density in Stan - expect_true(grep(stan_dist, brms::stancode(brms_prior_empty)) == 1, info="Looking for declared Stan mixture density in generated brms Stan code.") + stan_code <- brms::stancode(brms_prior_empty) + expect_true(grep(stan_dist_lpdf, stan_code) == 1, info="Looking for declared Stan mixture density pdf in generated brms Stan code.") + expect_true(grep(stan_dist_lcdf, stan_code) == 1, info="Looking for declared Stan mixture density cdf in generated brms Stan code.") + expect_true(grep(stan_dist_lccdf, stan_code) == 1, info="Looking for declared Stan mixture density ccdf in generated brms Stan code.") + expect_true(grep(stan_dist_cdf, stan_code) == 1, info="Looking for declared Stan mixture density natural scale cdf in generated brms Stan code.") ## now check for the mixture being passed to Stan as data stan_data <- brms::standata(brms_prior_empty) for(i in 1:3) { @@ -71,27 +79,46 @@ brms_beta_args <- list(formula=brms::bf(r | trials(n) ~ 1, family=brms::brmsfami data=data.frame(r=0, n=0), prior=brms::prior(mixbeta(prior_w, prior_a, prior_b), coef=Intercept)) -test_that("Beta quantiles are correct for brms sampled prior", mixstanvar_simul_test(beta, brms_beta_args, eps, c(0.1, 0.9))) -test_that("Beta prior is declared correctly in brms generated model and data", mixstanvar_test(beta, brms_beta_args)) -test_that("Beta mixture quantiles are correct for brms sampled prior", mixstanvar_simul_test(betaMix, brms_beta_args, eps, c(0.1, 0.9))) -test_that("Beta mixture prior is declared correctly in brms generated model and data", mixstanvar_test(betaMix, brms_beta_args)) +test_that("Beta quantiles are correct for brms sampled prior", { mixstanvar_simul_test(beta, brms_beta_args, eps, c(0.1, 0.9)) }) +test_that("Beta prior is declared correctly in brms generated model and data", { mixstanvar_test(beta, brms_beta_args) }) +test_that("Beta mixture quantiles are correct for brms sampled prior", { mixstanvar_simul_test(betaMix, brms_beta_args, eps, c(0.1, 0.9)) }) +test_that("Beta mixture prior is declared correctly in brms generated model and data", { mixstanvar_test(betaMix, brms_beta_args) }) + +brms_beta_trunc_args <- list(formula=brms::bf(r | trials(n) ~ 1, family=brms::brmsfamily("binomial", link="identity"), center=FALSE), + data=data.frame(r=0, n=0), + prior=brms::prior(mixbeta(prior_w, prior_a, prior_b), class=b, lb=0.1, ub=0.9)) +test_that("Beta (truncated) prior is declared correctly in brms generated model and data", { mixstanvar_test(beta, brms_beta_trunc_args) }) +test_that("Beta mixture (truncated) prior is declared correctly in brms generated model and data", { mixstanvar_test(betaMix, brms_beta_trunc_args) }) brms_normal_args <- list(formula=brms::bf(y ~ 1, family=brms::brmsfamily("gaussian", link="identity"), center=FALSE), data=data.frame(y=0), prior=brms::prior(mixnorm(prior_w, prior_m, prior_s), coef=Intercept) + brms::prior(constant(1), class=sigma)) -test_that("Normal quantiles are correct for brms sampled prior", mixstanvar_simul_test(norm, brms_normal_args, eps, c(-1, 0))) -test_that("Normal prior is declared correctly in brms generated model and data", mixstanvar_test(norm, brms_normal_args)) -test_that("Normal mixture quantiles are correct for brms sampled prior", mixstanvar_simul_test(normMix, brms_normal_args, eps, c(2, 1), ptest=c(0.3, 0.5, 0.7))) -test_that("Normal mixture prior is declared correctly in brms generated model and data", mixstanvar_test(normMix, brms_normal_args)) +test_that("Normal quantiles are correct for brms sampled prior", { mixstanvar_simul_test(norm, brms_normal_args, eps, c(-1, 0)) }) +test_that("Normal prior is declared correctly in brms generated model and data", { mixstanvar_test(norm, brms_normal_args) }) +test_that("Normal mixture quantiles are correct for brms sampled prior", { mixstanvar_simul_test(normMix, brms_normal_args, eps, c(2, 1), ptest=c(0.3, 0.5, 0.7)) }) +test_that("Normal mixture prior is declared correctly in brms generated model and data", { mixstanvar_test(normMix, brms_normal_args) }) + +brms_normal_trunc_args <- list(formula=brms::bf(y ~ 1, family=brms::brmsfamily("gaussian", link="identity"), center=FALSE), + data=data.frame(y=0), + prior=brms::prior(mixnorm(prior_w, prior_m, prior_s), class=b, lb=-5, ub=5) + brms::prior(constant(1), class=sigma)) +test_that("Normal (truncated) prior is declared correctly in brms generated model and data", { mixstanvar_test(norm, brms_normal_trunc_args) }) +test_that("Normal mixture (truncated) prior is declared correctly in brms generated model and data", { mixstanvar_test(normMix, brms_normal_trunc_args) }) brms_gamma_args <- list(formula=brms::bf(y ~ 1, family=brms::brmsfamily("gaussian", link="identity"), center=FALSE), data=data.frame(y=1), prior=brms::prior(mixgamma(prior_w, prior_a, prior_b), coef=Intercept) + brms::prior(constant(1), class=sigma)) -test_that("Gamma quantiles are correct for brms sampled prior", mixstanvar_simul_test(gamma, brms_gamma_args, eps, c(2, 7))) -test_that("Gamma prior is declared correctly in brms generated model and data", mixstanvar_test(gamma, brms_gamma_args)) -test_that("Gamma mixture quantile function is correct for brms sampled prior", mixstanvar_simul_test(gammaMix, brms_gamma_args, eps, c(2, 7), ptest = seq(0.2, 0.8, by=0.2))) -test_that("Gamma mixture prior is declared correctly in brms generated model and data", mixstanvar_test(gammaMix, brms_gamma_args)) +test_that("Gamma quantiles are correct for brms sampled prior", { mixstanvar_simul_test(gamma, brms_gamma_args, eps, c(2, 7)) }) +test_that("Gamma prior is declared correctly in brms generated model and data", { mixstanvar_test(gamma, brms_gamma_args) }) +test_that("Gamma mixture quantile function is correct for brms sampled prior", { mixstanvar_simul_test(gammaMix, brms_gamma_args, eps, c(2, 7), ptest = seq(0.2, 0.8, by=0.2)) }) +test_that("Gamma mixture prior is declared correctly in brms generated model and data", { mixstanvar_test(gammaMix, brms_gamma_args) }) + +brms_gamma_trunc_args <- list(formula=brms::bf(y ~ 1, family=brms::brmsfamily("gaussian", link="identity"), center=FALSE), + data=data.frame(y=1), + prior=brms::prior(mixgamma(prior_w, prior_a, prior_b), class=b, lb=0.1, ub=10) + brms::prior(constant(1), class=sigma)) + +test_that("Gamma (truncated) prior is declared correctly in brms generated model and data", { mixstanvar_test(gamma, brms_gamma_trunc_args) }) +test_that("Gamma mixture (truncated) prior is declared correctly in brms generated model and data", { mixstanvar_test(gammaMix, brms_gamma_trunc_args) }) # Here we approximate the samples using a multi-variante normal via a # moment based approxmation and compare this to the respective @@ -170,10 +197,10 @@ brms_mvn_4_args <- list(formula=brms::bf(y ~ 1 + l1 + l2 + l3, family=brms::brms data=data.frame(y=1, l1=0, l2=0, l3=0), prior=brms::prior(mixmvnorm(prior_w, prior_m, prior_sigma_L), class=b) + brms::prior(constant(1), class=sigma)) -test_that("Multivariate normal (4D) is correct for brms sampled prior", mixstanvar_simul_mv_test(mvnorm_single_4, brms_mvn_4_args, eps)) -test_that("Multivariate normal (4D) prior is declared correctly in brms generated model and data", mixstanvar_test_mvnormMix(mvnorm_single_4, brms_mvn_4_args)) -test_that("Multivariate normal with heavy (4D) tails is correct for brms sampled prior", mixstanvar_simul_mv_test(mvnorm_heavy_4, brms_mvn_4_args, eps)) -test_that("Multivariate normal with heavy (4D) tails is declared correctly in brms generated model and data", mixstanvar_test_mvnormMix(mvnorm_heavy_4, brms_mvn_4_args)) +test_that("Multivariate normal (4D) is correct for brms sampled prior", { mixstanvar_simul_mv_test(mvnorm_single_4, brms_mvn_4_args, eps) }) +test_that("Multivariate normal (4D) prior is declared correctly in brms generated model and data", { mixstanvar_test_mvnormMix(mvnorm_single_4, brms_mvn_4_args) }) +test_that("Multivariate normal with heavy (4D) tails is correct for brms sampled prior", { mixstanvar_simul_mv_test(mvnorm_heavy_4, brms_mvn_4_args, eps) }) +test_that("Multivariate normal with heavy (4D) tails is declared correctly in brms generated model and data", { mixstanvar_test_mvnormMix(mvnorm_heavy_4, brms_mvn_4_args) }) mvnorm_single_2 <- mixmvnorm(c(1, m1[1:2], 5), param="mn", sigma=S[1:2,1:2]) @@ -187,7 +214,7 @@ brms_mvn_2_args <- list(formula=brms::bf(y ~ 1 + l1, family=brms::brmsfamily("ga data=data.frame(y=1, l1=0), prior=brms::prior(mixmvnorm(prior_w, prior_m, prior_sigma_L), class=b) + brms::prior(constant(1), class=sigma)) -test_that("Multivariate normal (2D) is correct for brms sampled prior", mixstanvar_simul_mv_test(mvnorm_single_2, brms_mvn_2_args, eps)) -test_that("Multivariate normal (2D) is declared correctly in brms generated model and data", mixstanvar_test_mvnormMix(mvnorm_single_2, brms_mvn_2_args)) -test_that("Multivariate normal with heavy (2D) tails is correct for brms sampled prior", mixstanvar_simul_mv_test(mvnorm_heavy_2, brms_mvn_2_args, eps)) -test_that("Multivariate normal with heavy (2D) is declared correctly in brms generated model and data", mixstanvar_test_mvnormMix(mvnorm_heavy_2, brms_mvn_2_args)) +test_that("Multivariate normal (2D) is correct for brms sampled prior", { mixstanvar_simul_mv_test(mvnorm_single_2, brms_mvn_2_args, eps) }) +test_that("Multivariate normal (2D) is declared correctly in brms generated model and data", { mixstanvar_test_mvnormMix(mvnorm_single_2, brms_mvn_2_args) }) +test_that("Multivariate normal with heavy (2D) tails is correct for brms sampled prior", { mixstanvar_simul_mv_test(mvnorm_heavy_2, brms_mvn_2_args, eps) }) +test_that("Multivariate normal with heavy (2D) is declared correctly in brms generated model and data", { mixstanvar_test_mvnormMix(mvnorm_heavy_2, brms_mvn_2_args) }) diff --git a/tests/testthat/test-oc1S.R b/tests/testthat/test-oc1S.R index 0d16cc5..910470c 100644 --- a/tests/testthat/test-oc1S.R +++ b/tests/testthat/test-oc1S.R @@ -1,4 +1,3 @@ -context("oc1S: 1-Sample Operating Characteristics") ## test the analytical OC function via brute force simulation set.seed(12354) @@ -34,13 +33,13 @@ test_scenario <- function(oc_res, ref) { expect_true(all(abs(resA) < eps)) } -test_that("Classical NI design critical value", expect_true( abs(decision1S_boundary(prior, 155, decA) - c1) < eps)) +test_that("Classical NI design critical value", { expect_true( abs(decision1S_boundary(prior, 155, decA) - c1) < eps) }) ## n set to give power 80% to detect 0 and type I error 5% for no ## better than theta_ni -test_that("Classical NI design at target sample size for OCs", test_scenario(oc1S(prior, 155, decA)(thetaA), c(1-beta, alpha)) ) -test_that("Classical NI design at increased sample size for OCs", test_scenario(oc1S(prior, 233, decA)(thetaA), c(1-0.08, alpha)) ) -test_that("Classical NI design at decreased sample size for OCs", test_scenario(oc1S(prior, 77, decA)(thetaA), c(1-0.45, alpha)) ) +test_that("Classical NI design at target sample size for OCs", { test_scenario(oc1S(prior, 155, decA)(thetaA), c(1-beta, alpha)) }) +test_that("Classical NI design at increased sample size for OCs", { test_scenario(oc1S(prior, 233, decA)(thetaA), c(1-0.08, alpha)) }) +test_that("Classical NI design at decreased sample size for OCs", { test_scenario(oc1S(prior, 77, decA)(thetaA), c(1-0.45, alpha)) }) ## now double criterion with indecision point (mean estimate must be ## lower than this) @@ -57,17 +56,17 @@ thetaD <- c(theta_c, theta_ni) ## since theta_c == c1, both decision criteria are the same for n = ## 155 -test_that("Double criterion NI design at target sample size for OCs, combined ", test_scenario(oc1S(prior, 155, decComb)(thetaD), c(0.50, alpha)) ) -test_that("Double criterion NI design at target sample size for OCs, stat criterion", test_scenario(oc1S(prior, 155, dec1)(thetaD), c(0.50, alpha)) ) -test_that("Double criterion NI design at target sample size for OCs, mean criterion", test_scenario(oc1S(prior, 155, dec2)(thetaD), c(0.50, alpha)) ) +test_that("Double criterion NI design at target sample size for OCs, combined ", { test_scenario(oc1S(prior, 155, decComb)(thetaD), c(0.50, alpha)) }) +test_that("Double criterion NI design at target sample size for OCs, stat criterion", { test_scenario(oc1S(prior, 155, dec1)(thetaD), c(0.50, alpha)) }) +test_that("Double criterion NI design at target sample size for OCs, mean criterion", { test_scenario(oc1S(prior, 155, dec2)(thetaD), c(0.50, alpha)) }) ## at an increased sample size only the mean criterion is active -test_that("Double criterion NI design at increased sample size for OCs, combined ", test_scenario(oc1S(prior, 233, decComb)(thetaD), c(0.50, 0.02)) ) -test_that("Double criterion NI design at increased sample size for OCs, mean criterion", test_scenario(oc1S(prior, 233, dec2)(thetaD), c(0.50, 0.02)) ) +test_that("Double criterion NI design at increased sample size for OCs, combined ", { test_scenario(oc1S(prior, 233, decComb)(thetaD), c(0.50, 0.02)) }) +test_that("Double criterion NI design at increased sample size for OCs, mean criterion", { test_scenario(oc1S(prior, 233, dec2)(thetaD), c(0.50, 0.02)) }) ## at a decreased sample size only the stat criterion is active -test_that("Double criterion NI design at decreased sample size for OCs, combined ", test_scenario(oc1S(prior, 78, decComb)(thetaD), c(1-0.68, alpha)) ) -test_that("Double criterion NI design at decreased sample size for OCs, stat criterion", test_scenario(oc1S(prior, 78, dec1)(thetaD), c(1-0.68, alpha)) ) +test_that("Double criterion NI design at decreased sample size for OCs, combined ", { test_scenario(oc1S(prior, 78, decComb)(thetaD), c(1-0.68, alpha)) }) +test_that("Double criterion NI design at decreased sample size for OCs, stat criterion", { test_scenario(oc1S(prior, 78, dec1)(thetaD), c(1-0.68, alpha)) }) ## test type 1 error and correctness of critical values wrt to ## lower.tail=TRUE/FALSE @@ -95,12 +94,12 @@ design_binaryB <- oc1S(beta_prior, 1000, dec1b) crit1 <- decision1S_boundary(beta_prior, 1000, dec1) crit1B <- decision1S_boundary(beta_prior, 1000, dec1b) posterior_binary <- function(r) postmix(beta_prior, r=r, n=1000) -test_that("Binary type I error rate", test_scenario(design_binary(theta_ni), alpha)) -test_that("Binary crticial value, lower.tail=TRUE", test_critical_discrete(crit1, dec1, posterior_binary)) -test_that("Binary crticial value, lower.tail=FALSE", test_critical_discrete(crit1B, dec1b, posterior_binary)) +test_that("Binary type I error rate", { test_scenario(design_binary(theta_ni), alpha) }) +test_that("Binary crticial value, lower.tail=TRUE", { test_critical_discrete(crit1, dec1, posterior_binary) }) +test_that("Binary crticial value, lower.tail=FALSE", { test_critical_discrete(crit1B, dec1b, posterior_binary) }) -test_that("Binary boundary case, lower.tail=TRUE", expect_numeric(design_binary( 1), lower=0, upper=1, finite=TRUE, any.missing=FALSE)) -test_that("Binary boundary case, lower.tail=FALSE", expect_numeric(design_binaryB(0), lower=0, upper=1, finite=TRUE, any.missing=FALSE)) +test_that("Binary boundary case, lower.tail=TRUE", { expect_numeric(design_binary( 1), lower=0, upper=1, finite=TRUE, any.missing=FALSE) }) +test_that("Binary boundary case, lower.tail=FALSE", { expect_numeric(design_binaryB(0), lower=0, upper=1, finite=TRUE, any.missing=FALSE) }) ## poisson case gamma_prior <- mixgamma(c(1, 2, 2)) @@ -111,6 +110,6 @@ design_poissonB <- oc1S(gamma_prior, 1000, dec_countB) pcrit1 <- decision1S_boundary(gamma_prior, 1000, dec_count) pcrit1B <- decision1S_boundary(gamma_prior, 1000, dec_countB) posterior_poisson <- function(m) postmix(gamma_prior, m=m/1000, n=1000) -test_that("Poisson type I error rate", test_scenario(design_poisson(1), alpha) ) -test_that("Poisson critical value, lower.tail=TRUE", test_critical_discrete(pcrit1, dec_count, posterior_poisson)) -test_that("Poisson critical value, lower.tail=FALSE", test_critical_discrete(pcrit1B, dec_countB, posterior_poisson)) +test_that("Poisson type I error rate", { test_scenario(design_poisson(1), alpha) } ) +test_that("Poisson critical value, lower.tail=TRUE", { test_critical_discrete(pcrit1, dec_count, posterior_poisson) }) +test_that("Poisson critical value, lower.tail=FALSE", { test_critical_discrete(pcrit1B, dec_countB, posterior_poisson) }) diff --git a/tests/testthat/test-oc2S.R b/tests/testthat/test-oc2S.R index 4255f1c..1978bef 100644 --- a/tests/testthat/test-oc2S.R +++ b/tests/testthat/test-oc2S.R @@ -1,4 +1,3 @@ -context("oc2S: 2-Sample Operating Characteristics") ## test the analytical OC function via brute force simulation set.seed(12354) @@ -20,8 +19,7 @@ theta2 <- 0.5 Nsim <- 1e4 -run_on_cran <- function() -{ +run_on_cran <- function() { if (identical(Sys.getenv("NOT_CRAN"), "true")) { return(FALSE) } @@ -87,7 +85,7 @@ test_that("Type I error is matching between MC and analytical computations in th test_that("Power is matching between MC and analytical computations in the normal mixture case", { skip_on_cran() - power <- oc2S(prior1, prior2, N1, N2, decision2S(pcrit, qcrit))(theta1, theta2) + power <- oc2S(prior1, prior2, N1, N2, decision2S(pcrit, qcrit), sigma1=sigma(prior1), sigma2=sigma(prior2))(theta1, theta2) powerMC <- oc2S_normal_MC(prior1, prior2, N1, N2, theta1, theta2, pcrit, qcrit) res <- 100 * abs(power - powerMC) expect_equal(sum(res > 2) , 0) @@ -125,10 +123,10 @@ test_that("Gsponer et al. results match (normal end-point)", { ## spline function ## Table 1, probability for interim for success - oc$success <- oc2S(priorP, priorT, nP1, nT1, successCrit, Ngrid=1)(-49, -49-oc$delta) + oc$success <- oc2S(priorP, priorT, nP1, nT1, successCrit, Ngrid=1, sigma1=sigmaFixed, sigma2=sigmaFixed)(-49, -49-oc$delta) ## Table 1, probability for interim for futility - oc$futile <- oc2S(priorP, priorT, nP1, nT1, futilityCrit, Ngrid=1)(-49, -49-oc$delta) + oc$futile <- oc2S(priorP, priorT, nP1, nT1, futilityCrit, Ngrid=1, sigma1=sigmaFixed, sigma2=sigmaFixed)(-49, -49-oc$delta) ## Table 1, first three columns, page 74 oc[-1] <- lapply(100*oc[-1], round, 1) @@ -236,9 +234,9 @@ expect_equal_each <- function(test, expected) { ## design object, decision function, posterior function must return ## posterior after updatding the prior with the given value; we assume ## that the priors are the same for sample 1 and 2 -test_critical_discrete <- function(design, decision, posterior, y2) { +test_critical_discrete <- function(boundary_design, decision, posterior, y2) { lower.tail <- attr(decision, "lower.tail") - crit <- design(y2=y2) + crit <- boundary_design(y2) post2 <- posterior(y2) if(lower.tail) { expect_equal(decision(posterior(crit-1), post2), 1) @@ -263,13 +261,20 @@ decB <- decision2S(1-alpha, 0, lower.tail=FALSE) beta_prior <- mixbeta(c(1, 1, 1)) if(!run_on_cran()) { design_binary <- oc2S(beta_prior, beta_prior, 100, 100, dec) + boundary_design_binary <- decision2S_boundary(beta_prior, beta_prior, 100, 100, dec) design_binaryB <- oc2S(beta_prior, beta_prior, 100, 100, decB) + boundary_design_binaryB <- decision2S_boundary(beta_prior, beta_prior, 100, 100, decB) +} else { + design_binary <- function(...) { return(0.1) } + design_binaryB <- function(...) { return(0.1) } + boundary_design_binary <- function(...) { return(0.1) } + boundary_design_binaryB <- function(...) { return(0.1) } } posterior_binary <- function(r) postmix(beta_prior, r=r, n=100) p_test <- 1:9 / 10 test_that("Binary type I error rate", { skip_on_cran(); test_scenario(design_binary(p_test, p_test), alpha) }) -test_that("Binary crticial value, lower.tail=TRUE", { skip_on_cran(); test_critical_discrete(design_binary, dec, posterior_binary, 30) }) -test_that("Binary crticial value, lower.tail=FALSE", { skip_on_cran(); test_critical_discrete(design_binaryB, decB, posterior_binary, 30) }) +test_that("Binary crticial value, lower.tail=TRUE", { skip_on_cran(); test_critical_discrete(boundary_design_binary, dec, posterior_binary, 30) }) +test_that("Binary crticial value, lower.tail=FALSE", { skip_on_cran(); test_critical_discrete(boundary_design_binaryB, decB, posterior_binary, 30) }) test_that("Binary boundary case, lower.tail=TRUE", { skip_on_cran(); expect_numeric(design_binary( 1, 1), lower=0, upper=1, finite=TRUE, any.missing=FALSE) }) test_that("Binary boundary case, lower.tail=FALSE", { skip_on_cran(); expect_numeric(design_binaryB(0, 0), lower=0, upper=1, finite=TRUE, any.missing=FALSE) }) @@ -280,9 +285,11 @@ beta_prior1 <- mixbeta(c(1, 0.9, 1000), param="mn") beta_prior2 <- mixbeta(c(1, 0.1, 1000), param="mn") design_lower <- oc2S(beta_prior1, beta_prior2, 20, 20, dec) ## always 0 design_upper <- oc2S(beta_prior1, beta_prior2, 20, 20, decB) ## always 1 +boundary_design_lower <- decision2S_boundary(beta_prior1, beta_prior2, 20, 20, dec) ## always 0 +boundary_design_upper <- decision2S_boundary(beta_prior1, beta_prior2, 20, 20, decB) ## always 1 -test_that("Binary case, no decision change, lower.tail=TRUE, critical value", { skip_on_cran(); expect_equal_each(design_lower(y2=0:20), -1) }) -test_that("Binary case, no decision change, lower.tail=FALSE, critical value", { skip_on_cran(); expect_equal_each(design_upper(y2=0:20), 21) }) +test_that("Binary case, no decision change, lower.tail=TRUE, critical value", { skip_on_cran(); expect_equal_each(boundary_design_lower(0:20), -1) }) +test_that("Binary case, no decision change, lower.tail=FALSE, critical value", { skip_on_cran(); expect_equal_each(boundary_design_upper(0:20), 21) }) test_that("Binary case, no decision change, lower.tail=TRUE, frequency=0", { skip_on_cran(); expect_equal_each(design_lower(p_test, p_test), 0.0) }) test_that("Binary case, no decision change, lower.tail=FALSE, frequency=1", { skip_on_cran(); expect_equal_each(design_upper(p_test, p_test), 1.0) }) @@ -290,10 +297,17 @@ test_that("Binary case, no decision change, lower.tail=FALSE, frequency=1", { s if(!run_on_cran()) { design_lower_rev <- oc2S(beta_prior2, beta_prior1, 20, 20, dec) ## always 1 design_upper_rev <- oc2S(beta_prior2, beta_prior1, 20, 20, decB) ## always 0 + boundary_design_lower_rev <- decision2S_boundary(beta_prior2, beta_prior1, 20, 20, dec) ## always 1 + boundary_design_upper_rev <- decision2S_boundary(beta_prior2, beta_prior1, 20, 20, decB) ## always 0 +} else { + design_lower_rev <- function(...) return(1) + design_upper_rev <- function(...) return(0) + boundary_design_lower_rev <- function(...) return(1) + boundary_design_upper_rev <- function(...) return(0) } -test_that("Binary case, no decision change (reversed), lower.tail=TRUE, critical value", { skip_on_cran(); expect_equal_each(design_lower_rev(y2=0:20), 20) }) -test_that("Binary case, no decision change (reversed), lower.tail=FALSE, critical value", { skip_on_cran(); expect_equal_each(design_upper_rev(y2=0:20), -1) }) +test_that("Binary case, no decision change (reversed), lower.tail=TRUE, critical value", { skip_on_cran(); expect_equal_each(boundary_design_lower_rev(0:20), 20) }) +test_that("Binary case, no decision change (reversed), lower.tail=FALSE, critical value", { skip_on_cran(); expect_equal_each(boundary_design_upper_rev(0:20), -1) }) test_that("Binary case, no decision change (reversed), lower.tail=TRUE, frequency=0", { skip_on_cran(); expect_equal_each(design_lower_rev(p_test, p_test), 1.0) }) test_that("Binary case, no decision change (reversed), lower.tail=FALSE, frequency=1", { skip_on_cran(); expect_equal_each(design_upper_rev(p_test, p_test), 0.0) }) test_that("Binary case, log-link", { @@ -346,11 +360,13 @@ gamma_prior <- mixgamma(c(1, 2, 2)) design_poisson <- oc2S(gamma_prior, gamma_prior, 100, 100, dec) design_poissonB <- oc2S(gamma_prior, gamma_prior, 100, 100, decB) +boundary_design_poisson <- decision2S_boundary(gamma_prior, gamma_prior, 100, 100, dec) +boundary_design_poissonB <- decision2S_boundary(gamma_prior, gamma_prior, 100, 100, decB) posterior_poisson <- function(m) postmix(gamma_prior, m=m/100, n=100) lambda_test <- seq(0.5, 1.3, by=0.1) test_that("Poisson type I error rate", { skip_on_cran(); test_scenario(design_poisson(lambda_test, lambda_test), alpha) }) -test_that("Poisson crticial value, lower.tail=TRUE", { skip_on_cran(); test_critical_discrete(design_poisson, dec, posterior_poisson, 90) }) -test_that("Poisson crticial value, lower.tail=FALSE", { skip_on_cran(); test_critical_discrete(design_poissonB, decB, posterior_poisson, 90) }) +test_that("Poisson crticial value, lower.tail=TRUE", { skip_on_cran(); test_critical_discrete(boundary_design_poisson, dec, posterior_poisson, 90) }) +test_that("Poisson crticial value, lower.tail=FALSE", { skip_on_cran(); test_critical_discrete(boundary_design_poissonB, decB, posterior_poisson, 90) }) ## 22 Nov 2017: disabled test as we trigger always calculation of the ## boundaries as of now. ##test_that("Poisson results cache expands", { @@ -413,8 +429,8 @@ test_that("Normal OC 2-sample avoids undefined behavior, example 1", { dec <- decision2S(1-alpha, 0, lower.tail=FALSE) n <- 58 k <- 2 - design_map <- oc2S(prior_flat, map_ref, n, n/k, dec) - design_map_2 <- oc2S(prior_flat, map_ref, n, n/k, dec) + design_map <- oc2S(prior_flat, map_ref, n, n/k, dec, sigma1=sigma_ref, sigma2=sigma_ref) + design_map_2 <- oc2S(prior_flat, map_ref, n, n/k, dec, sigma1=sigma_ref, sigma2=sigma_ref) x <- seq(-2.6, -1.6, by=0.1) expect_numeric(design_map(x, x), lower=0, upper=1, any.missing=FALSE) diff --git a/tests/testthat/test-pos1S.R b/tests/testthat/test-pos1S.R index 45c620b..151c469 100644 --- a/tests/testthat/test-pos1S.R +++ b/tests/testthat/test-pos1S.R @@ -1,5 +1,3 @@ -context("pos1S: 1-Sample Probability of Success") - ## test the analytical OC function via brute force simulation set.seed(12354) @@ -50,19 +48,18 @@ test_pos1S <- function(prior, ia_dist, n, dec, decU) { } -test_that("Normal PoS 1 sample function matches MC integration of CPO", test_pos1S(prior, ia_dist, n1, decA, decAU)) +test_that("Normal PoS 1 sample function matches MC integration of CPO", { test_pos1S(prior, ia_dist, n1, decA, decAU) }) beta_prior <- mixbeta(c(1, 1, 1)) beta_ia <- postmix(beta_prior, r=20, n=50) -test_that("Binomial PoS 1 sample function matches MC integration of CPO", test_pos1S(beta_prior, beta_ia, n1, decA, decAU)) +test_that("Binomial PoS 1 sample function matches MC integration of CPO", { test_pos1S(beta_prior, beta_ia, n1, decA, decAU) }) gamma_prior <- mixgamma(c(1, 1, 1), param="mn") dec_count <- decision1S(1-alpha, 1, lower.tail=TRUE) dec_countU <- decision1S(1-alpha, 1, lower.tail=FALSE) gamma_ia <- postmix(gamma_prior, m=0.9, n=40) -test_that("Poisson PoS 1 sample function matches MC integration of CPO", test_pos1S(gamma_prior, gamma_ia, n1, dec_count, dec_countU)) +test_that("Poisson PoS 1 sample function matches MC integration of CPO", { test_pos1S(gamma_prior, gamma_ia, n1, dec_count, dec_countU) }) prior_unit_inf <- mixnorm(c(1, 0, 1), sigma=s, param="mn") post_ia_unit_inf <- postmix(prior_unit_inf, m=-1, n=162) -test_pos1S(prior_unit_inf, post_ia_unit_inf, 459-162, decA, decAU) -test_that("Normal PoS 1 sample function matches MC integration of CPO (more extreme case)", test_pos1S(prior_unit_inf, post_ia_unit_inf, 459-162, decA, decAU)) +test_that("Normal PoS 1 sample function matches MC integration of CPO (more extreme case)", { test_pos1S(prior_unit_inf, post_ia_unit_inf, 459-162, decA, decAU) }) diff --git a/tests/testthat/test-pos2S.R b/tests/testthat/test-pos2S.R index 2f52a97..049086d 100644 --- a/tests/testthat/test-pos2S.R +++ b/tests/testthat/test-pos2S.R @@ -117,3 +117,19 @@ test_that("Poisson PoS 2 sample function matches MC integration of CPO", N1, N2, dec_count, dec_countU)) +test_that("Binomial PoS 2 with IA returns results", { + ## reported by user + successCrit <- decision2S(c(0.9), c(0), lower.tail=FALSE) + n0 <- 50 + n <- 100 + n_alt <- 140 + neutr_prior <- mixbeta(c(1,1/3,1/3)) + post_placeboIA <- postmix(neutr_prior, r=13, n=n0) + post_treatIA <- postmix(neutr_prior, r=3, n=n0) + # Criterion for PPoS at IA + pos_final <- pos2S(post_treatIA, post_placeboIA, n-n0, n-n0, successCrit) + pos_final_alt <- pos2S(post_treatIA, post_placeboIA, n_alt-n0, n_alt-n0, successCrit) + #Predictive proba of success at the end + expect_number(pos_final_alt(post_treatIA, post_placeboIA), na.ok=FALSE, lower=0, upper=1, finite=TRUE, null.ok=FALSE) + expect_number(pos_final(post_treatIA, post_placeboIA), na.ok=FALSE, lower=0, upper=1, finite=TRUE, null.ok=FALSE) +}) diff --git a/tests/testthat/test-postmix.R b/tests/testthat/test-postmix.R index c2a8ad8..2492979 100644 --- a/tests/testthat/test-postmix.R +++ b/tests/testthat/test-postmix.R @@ -1,4 +1,3 @@ -context("postmix: Posterior Mixture Distribution") norm <- mixnorm(c(1, 0, 0.5), sigma=1) diff --git a/tests/testthat/test-preddist.R b/tests/testthat/test-preddist.R index 635eebc..0da4b37 100644 --- a/tests/testthat/test-preddist.R +++ b/tests/testthat/test-preddist.R @@ -1,4 +1,3 @@ -context("preddist: Mixture Predictive Distribution") ## check that predictive distributions hold what they promise, ## i.e. that they describe the sum of n new data points. @@ -71,11 +70,11 @@ preddist_cmp <- function(mix, n, n_rng, N=Nsamp, qntls=p_quants, stat=c("sum", } -test_that("Predictive for a beta evaluates correctly (binary)", preddist_cmp(beta, n, Curry(rbinom, n=n, size=1))) -test_that("Predictive for a beta mixture evaluates correctly (binary)", preddist_cmp(betaMix, n, Curry(rbinom, n=n, size=1))) +test_that("Predictive for a beta evaluates correctly (binary)", { preddist_cmp(beta, n, Curry(rbinom, n=n, size=1)) }) +test_that("Predictive for a beta mixture evaluates correctly (binary)", { preddist_cmp(betaMix, n, Curry(rbinom, n=n, size=1)) }) -test_that("Predictive for a gamma evaluates correctly (poisson)", preddist_cmp(gamma, n, Curry(rpois, n=n))) -test_that("Predictive for a gamma mixture evaluates correctly (poisson)", preddist_cmp(gammaMix, n, Curry(rpois, n=n), Teps=1E-1)) +test_that("Predictive for a gamma evaluates correctly (poisson)", { preddist_cmp(gamma, n, Curry(rpois, n=n)) }) +test_that("Predictive for a gamma mixture evaluates correctly (poisson)", { preddist_cmp(gammaMix, n, Curry(rpois, n=n), Teps=1E-1) }) -test_that("Predictive for a normal evaluates correctly (normal)", preddist_cmp(norm, n, Curry(rnorm, n=n, sd=sigma(norm)), stat="mean")) -test_that("Predictive for a normal mixture evaluates correctly (normal)", preddist_cmp(normMix, n, Curry(rnorm, n=n, sd=sigma(normMix)), stat="mean", Teps=1E-1)) +test_that("Predictive for a normal evaluates correctly (normal)", { preddist_cmp(norm, n, Curry(rnorm, n=n, sd=sigma(norm)), stat="mean") }) +test_that("Predictive for a normal mixture evaluates correctly (normal)", { preddist_cmp(normMix, n, Curry(rnorm, n=n, sd=sigma(normMix)), stat="mean", Teps=1E-1) }) diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R index 2d0e391..81ea59e 100644 --- a/tests/testthat/test-utils.R +++ b/tests/testthat/test-utils.R @@ -1,10 +1,9 @@ -context("Utilities: predict") ## run the example from predict.gMAP -suppressMessages(example("predict.gMAP", package="RBesT", echo=FALSE, ask=FALSE, verbose=FALSE)) +source_example("predict_gMAP.R") ## check that we got for each input data item a prediction -test_that("correct # of predictions are generated", expect_equal(nrow(map$data), ncol(samp))) +test_that("correct # of predictions are generated", { expect_equal(nrow(map$data), ncol(samp)) }) ## check that the predictive distribution has a variance which is ## larger in accordance to the betwee-trial heterogeniety (needs to be @@ -28,7 +27,7 @@ test_that("variances have correct ordering", { ## new prediction was done for a single data item -test_that("correct # of new predictions are generated", expect_equal(ncol(pred_new), 1)) +test_that("correct # of new predictions are generated", { expect_equal(ncol(pred_new), 1) } ) ## must have larger sd than between-trial alone (on link scale) test_that("predictive variances have correct ordering",{ @@ -46,8 +45,6 @@ test_that("predictive distributions for the same study & covariate must match ex expect_equal(post_trans[,1], post_trans[,2]) }) -context("Utilities: (auto)mixfit") - test_that("automixfit attempts K=4 different models and returns best fitting", { auto_map <- automixfit(map, Nc=1:4, k=6) models <- attr(auto_map, "models") @@ -75,16 +72,14 @@ test_that("mixfit for prediction handles response and link scale", { }) -context("Utilities: mixcombine") - -example("mixcombine", package="RBesT", echo=FALSE, ask=FALSE, verbose=FALSE) +source_example("mixcombine.R") test_that("combination of mixtures", { m1 <- mixcombine(bm, unif, weight=c(9, 1)) m2 <- mixcombine(bm, unif, unif, weight=c(8, 1, 1)) - expect_equivalent(m1[1,], c(bm[1,] - 0.1/2, 0.1)) - expect_equivalent(m1[2:3,1:2], bm[2:3,1:2]) - expect_equivalent(m2[2:3,1:2], bm[2:3,1:2]) + expect_equal(m1[1,], c(bm[1,] - 0.1/2, 0.1), ignore_attr=TRUE) + expect_equal(m1[2:3,1:2], bm[2:3,1:2], ignore_attr=TRUE) + expect_equal(m2[2:3,1:2], bm[2:3,1:2], ignore_attr=TRUE) }) test_that("throws an error if more weights than mixtures given", { @@ -98,50 +93,47 @@ test_that("combination of normal mixtures without default sigma works", { expect_true(ncol(norm_ui_mix) == 2) }) -context("Utilities: robustify") - -example("robustify", package="RBesT", echo=FALSE, ask=FALSE, verbose=FALSE) +source_example("robustify.R") test_that("beta mixture is robustified with Beta(1,1)", { expect_equal(ncol(bmix)+1, ncol(rbmix)) - expect_equivalent(rbmix[,ncol(rbmix)], c(0.1, 1, 1)) + expect_equal(rbmix[,ncol(rbmix)], c(0.1, 1, 1), ignore_attr=TRUE) }) test_that("beta mixture is robustified with Beta(0.5,0.5)", { - rbmix2 <- robustify(bmix, w=0.1, n=0) + rbmix2 <- robustify(bmix, w=0.1, n=0, mean=0.5) expect_equal(ncol(bmix)+1, ncol(rbmix2)) - expect_equivalent(rbmix2[,ncol(rbmix2)], c(0.1, 0.5, 0.5)) + expect_equal(rbmix2[,ncol(rbmix2)], c(0.1, 0.5, 0.5), ignore_attr=TRUE) }) test_that("gamma mixture is robustified with n=1 equivalent prior", { m <- summary(gmnMix)["mean"] nr <- ncol(rgmnMix) - expect_equivalent(rgmnMix[[nr, rescale=TRUE]], mixgamma(c(1, m, 1), param="mn")) + expect_equal(rgmnMix[[nr, rescale=TRUE]], mixgamma(c(1, m, 1), param="mn"), ignore_attr=TRUE) expect_equal(rgmnMix[1,nr], 0.1) }) test_that("gamma mixture is robustified with n=5 equivalent prior", { m <- summary(gmnMix)["mean"] - rgmnMix2 <- robustify(gmnMix, w=0.1, n=5) + rgmnMix2 <- robustify(gmnMix, w=0.1, n=5, mean=2) nr <- ncol(rgmnMix2) - expect_equivalent(rgmnMix2[[nr, rescale=TRUE]], mixgamma(c(1, m, 5), param="mn")) + expect_equal(rgmnMix2[[nr, rescale=TRUE]], mixgamma(c(1, m, 5), param="mn"), ignore_attr=TRUE) expect_equal(rgmnMix2[1,nr], 0.1) }) test_that("normal mixture is robustified with n=1 equivalent prior", { nr <- ncol(rnMix) - expect_equivalent(rnMix[[nr, rescale=TRUE]], mixnorm(c(1, 0, 1), param="mn", sigma=sigma(nm))) + expect_equal(rnMix[[nr, rescale=TRUE]], mixnorm(c(1, 0, 1), param="mn", sigma=sigma(nm)), ignore_attr=TRUE) expect_equal(rnMix[1,nr], 0.1) }) test_that("normal mixture is robustified with n=5 equivalent prior", { rnMix2 <- robustify(nm, w=0.1, mean=0, n=5, sigma=sigma(nm)) nr <- ncol(rnMix2) - expect_equivalent(rnMix2[[nr, rescale=TRUE]], mixnorm(c(1, 0, 5), param="mn", sigma=sigma(nm))) + expect_equal(rnMix2[[nr, rescale=TRUE]], mixnorm(c(1, 0, 5), param="mn", sigma=sigma(nm)), ignore_attr=TRUE) expect_equal(rnMix2[1,nr], 0.1) }) -context("Utilities: Plotting of Mixtures") test_that("plotting of normal mixtures without default sigma works", { norm_ui <- mixnorm(c(1, 0, 2)) norm_mix_ui <- mixcombine(norm_ui, norm_ui, weight=c(0.5,0.5)) @@ -149,9 +141,7 @@ test_that("plotting of normal mixtures without default sigma works", { expect_true(inherits(pl, "ggplot")) }) -context("Utilities: Mixture Effective Sample Size") - -example("ess", package="RBesT", echo=FALSE, ask=FALSE, verbose=FALSE) +source_example("ess.R") test_that("conjugate beta case matches canonical formula", { expect_equal(a+b, ess(prior, "moment")) @@ -173,7 +163,7 @@ test_that("ess elir for beta mixtures gives a warning for a<1 & b<1 densities", test_that("ess elir for normal mixtures returns correct values", { mix <- mixnorm( inf1=c(0.5026,-191.1869,127.4207),inf2=c(0.2647,-187.5895,31.6130),inf3=c(0.2326,-184.7445,345.3849), sigma=270.4877) - expect_gt(ess(mix), 0) + expect_gt(ess(mix, sigma=270.4877), 0) }) test_that("moment matching for beta mixtures is correct", { @@ -297,3 +287,31 @@ test_that("ESS elir is predictively consistent for gamma mixtures (Poisson likel gmixP <- mixgamma(rob=c(0.3, 20, 4), inf=c(0.7, 50, 10), likelihood="poisson") elir_predictive_consistent(gmixP, m=1E2, Nsim=1E3, seed=355435, stat="m", n=1E2) }) + +test_that("ess elir for problematic beta mixtures gives correct result 1", { + ## by user reported beta mixture density which triggers this erros + ## with RBesT 1.7.2 & 1.7.3 (others not tested): + ## Error in if (all(dgl < 0) || all(dgl > 0)) { : + ## missing value where TRUE/FALSE needed + + mixmat <- matrix(c(0.06429517, 0.03301215, 0.00269268, 0.90000000, + 437.32302999, 64.04211307, 5.92543558, 1.00000000, + 10.71709277, 2.14157953, 1.00000001, 1.00000000), byrow=TRUE, ncol=4) + + mixb <- do.call(mixbeta, apply(mixmat,2,c,simplify=FALSE)) + + expect_double(ess(mixb), lower=0, finite=TRUE, any.missing=FALSE, len=1) +}) + +test_that("ess elir for problematic beta mixtures gives correct result 2", { + mixmat <- matrix(c(0.7237396, 0.1665037, 0.1097567, + 53.3721902, 44.3894573, 9.8097062, + 1.4301638, 4.3842200, 1.8492197 + ), byrow=TRUE, ncol=3) + + mixb <- do.call(mixbeta, apply(mixmat,2,c,simplify=FALSE)) + + expect_double(ess(robustify(mixb, 0.05, 0.5)), lower=0, finite=TRUE, any.missing=FALSE, len=1) + expect_double(ess(robustify(mixb, 0.95, 0.5)), lower=0, finite=TRUE, any.missing=FALSE, len=1) +}) + diff --git a/tools/make-ds.R b/tools/make-ds.R index a7274c9..1ceff44 100644 --- a/tools/make-ds.R +++ b/tools/make-ds.R @@ -43,7 +43,7 @@ make_internal_ds <- function() { calibration_meta["MD5"] <- vals["MD5"] pkg_create_date <- Sys.time() - pkg_sha <- "511a0f1" + pkg_sha <- "87eef77" if (gsub("\\$", "", pkg_sha) == "Format:%h") { pkg_sha <- system("git rev-parse --short HEAD", intern=TRUE) diff --git a/vignettes/articles/introduction_normal.Rmd b/vignettes/articles/introduction_normal.Rmd index bf49a9e..b102d97 100644 --- a/vignettes/articles/introduction_normal.Rmd +++ b/vignettes/articles/introduction_normal.Rmd @@ -237,11 +237,13 @@ ocI <- rbind(data.frame(cfb_truth=cfb_truth, typeI=typeI1, data.frame(cfb_truth=cfb_truth, typeI=typeI4, design="40:20 with robust prior for placebo") ) -qplot(cfb_truth, typeI, data=ocI, colour=design, geom="line", main="Type I Error") + +ggplot(ocI, aes(cfb_truth, typeI, colour=design)) + + geom_line() + + ggtitle("Type I Error") + xlab(expression(paste('True value of change from baseline ', mu[act] == mu[pbo]))) + - ylab('Type I error') + - coord_cartesian(ylim=c(0,0.2)) + - theme(legend.justification=c(1,1),legend.position=c(0.95,0.85)) + ylab('Type I error') + + coord_cartesian(ylim=c(0,0.2)) + + theme(legend.justification=c(1,1),legend.position=c(0.95,0.85)) ``` @@ -274,17 +276,17 @@ ocP <- rbind(data.frame(cfb_truth1=cfb_truth1, cfb_truth2=cfb_truth2, design="40:20 with non-robust prior for placebo"), data.frame(cfb_truth1=cfb_truth1, cfb_truth2=cfb_truth2, delta=delta, power=power4, - design="40:20 with robust prior for placebo") -) - -qplot(delta, power, data=ocP, colour=design, geom="line", main="Power") + - xlab('True value of difference (act - pbo)')+ ylab('Power') + - scale_y_continuous(breaks=c(seq(0,1,0.2),0.9)) + - scale_x_continuous(breaks=c(seq(-80,0,20),-70)) + - geom_hline(yintercept=0.9,linetype=2) + - geom_vline(xintercept=-70,linetype=2) + - theme(legend.justification=c(1,1),legend.position=c(0.95,0.85)) - + design="40:20 with robust prior for placebo")) + +ggplot(ocP, aes(delta, power, colour=design)) + + geom_line() + + ggtitle("Power") + + xlab('True value of difference (act - pbo)')+ ylab('Power') + + scale_y_continuous(breaks=c(seq(0,1,0.2),0.9)) + + scale_x_continuous(breaks=c(seq(-80,0,20),-70)) + + geom_hline(yintercept=0.9,linetype=2) + + geom_vline(xintercept=-70,linetype=2) + + theme(legend.justification=c(1,1),legend.position=c(0.95,0.85)) ``` diff --git a/vignettes/introduction.Rmd b/vignettes/introduction.Rmd index 2d811c3..5e6cbf9 100644 --- a/vignettes/introduction.Rmd +++ b/vignettes/introduction.Rmd @@ -214,8 +214,8 @@ ess_weight <- rbind(ess_weight, data.frame(weight=c(0, 1), ess=c(ess(map), ess(mixbeta(c(1,1,1)))))) -qplot(weight, ess, data=ess_weight, geom=c("point", "line"), - main="ESS of robust MAP for varying weight of robust component") + +ggplot(ess_weight, aes(weight, ess)) + geom_point() + geom_line() + + ggtitle("ESS of robust MAP for varying weight of robust component") + scale_x_continuous(breaks=seq(0, 1, by=0.1)) + scale_y_continuous(breaks=seq(0, 40, by=5)) ``` @@ -290,7 +290,7 @@ ocI <- rbind(data.frame(theta=theta, typeI=typeI_robust, prior="robust"), data.frame(theta=theta, typeI=typeI_classic, prior="uniform 24:24") ) -qplot(theta, typeI, data=ocI, colour=prior, geom="line", main="Type I Error") +ggplot(ocI, aes(theta, typeI, colour=prior)) + geom_line() + ggtitle("Type I Error") ``` Note that observing response rates greater that 50% is highly @@ -304,8 +304,9 @@ Hence, it is resonable to restrict the response rates $\theta$ for which we evaluate the type I error to a a range of plausible values: ```{r} -qplot(theta, typeI, data=subset(ocI, theta < 0.5), colour=prior, geom="line", - main="Type I Error - response rate restricted to plausible range") +ggplot(ocI, aes(theta, typeI, colour=prior)) + geom_line() + + ggtitle("Type I Error - response rate restricted to plausible range") + + coord_cartesian(xlim=c(0, 0.5)) ``` ### Power @@ -331,7 +332,8 @@ ocP <- rbind(data.frame(theta_active, theta_control, delta=delta, power=power_ro data.frame(theta_active, theta_control, delta=delta, power=power_classic, prior="uniform 24:24") ) -qplot(delta, power, data=ocP, colour=prior, geom="line", main="Power") +ggplot(ocP, aes(delta, power, colour=prior)) + geom_line() + + ggtitle("Power") ``` @@ -378,7 +380,7 @@ ocC <- rbind(data.frame(y2=treat_y2, y1_crit=crit_robust(treat_y2), prior="ro data.frame(y2=treat_y2, y1_crit=crit_uniform(treat_y2), prior="uniform") ) -qplot(y2, y1_crit, data=ocC, colour=prior, geom="step", main="Critical values y1(y2)") +ggplot(ocC, aes(y2, y1_crit, colour=prior)) + geom_step() + ggtitle("Critical values y1(y2)") ``` The graph shows that the decision will always be negative if there are