From d53e6f559bf672e4ec367afafaa4bf32aef7ddcf Mon Sep 17 00:00:00 2001 From: Sebastian Weber Date: Wed, 20 Nov 2024 15:12:30 +0100 Subject: [PATCH 1/9] 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 From 0902b470263e8081ebc8526b7d30f891f9c9fb16 Mon Sep 17 00:00:00 2001 From: Sebastian Weber Date: Wed, 20 Nov 2024 15:40:22 +0100 Subject: [PATCH 2/9] update --- DESCRIPTION | 8 ++++---- Makefile | 5 +++++ NEWS.md | 2 ++ R/mixstanvar.R | 8 ++++---- tests/testthat/helper-utils.R | 10 ++++++++++ tools/make-ds.R | 2 +- 6 files changed, 26 insertions(+), 9 deletions(-) create mode 100644 tests/testthat/helper-utils.R diff --git a/DESCRIPTION b/DESCRIPTION index b1eb6c6..34bc702 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -9,7 +9,7 @@ Description: Tool-set to support Bayesian evidence synthesis. This and Schmidli et al. (2014) explain details on the methodology. Version: 1.7-4 -Date: 2024-11-19 +Date: 2024-11-20 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") @@ -25,7 +25,7 @@ Imports: methods, Rcpp (>= 0.12.0), RcppParallel (>= 5.0.1), - rstan (>= 2.26.0), + rstan (>= 2.32.0), rstantools (>= 2.4.0), assertthat, mvtnorm, @@ -44,8 +44,8 @@ LinkingTo: Rcpp (>= 0.12.0), RcppEigen (>= 0.3.3.3.0), RcppParallel (>= 5.0.1), - rstan (>= 2.26.0), - StanHeaders (>= 2.26.0) + rstan (>= 2.32.0), + StanHeaders (>= 2.32.0) License: GPL (>=3) LazyData: true Biarch: true diff --git a/Makefile b/Makefile index 59dfa27..3692867 100644 --- a/Makefile +++ b/Makefile @@ -156,6 +156,11 @@ r-source-check : r-source cd build; tar xvzf $(RPKG)-source.tar.gz cd build; NOT_CRAN=true $(RCMD) CMD check $(RPKG) +PHONY += r-source-release-check +r-source-release-check : r-source-release + cd build; tar xvzf $(RPKG)_$(PKG_VERSION).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 diff --git a/NEWS.md b/NEWS.md index a2013eb..0d2baf7 100644 --- a/NEWS.md +++ b/NEWS.md @@ -13,6 +13,8 @@ * Avoid over and underflow of beta-binomial distribution. * Replace calls to deprecated `ggplot2::qplot` with respective calls to `ggplot2::ggplot`. +* Use new `array` notation for `mixstanvar` generated Stan code to + make generated Stan programs compatible with Stan >= 2.33 # RBesT 1.7-3 - January 2nd, 2024 diff --git a/R/mixstanvar.R b/R/mixstanvar.R index 4894a20..6da68e8 100644 --- a/R/mixstanvar.R +++ b/R/mixstanvar.R @@ -136,7 +136,7 @@ mixstanvar <- function(..., verbose=FALSE) { if(includes_density("mvnormMix")) { sv <- sv + brms::stanvar(name="mixmvnorm_lpdf", scode=" -real mixmvnorm_lpdf(vector y, vector w, vector[] m, matrix[] L) { +real mixmvnorm_lpdf(vector y, vector w, array[] vector m, array[] matrix L) { int Nc = rows(w); vector[Nc] lp_mix; for(i in 1:Nc) { @@ -218,10 +218,10 @@ mix2brms.mvnormMix <- function(mix, name, verbose=FALSE) { mvprior <- brms::stanvar(Nc, glue::glue("{prefix}Nc")) + brms::stanvar(p, glue::glue("{prefix}p")) + brms::stanvar(array(mix[1,,drop=TRUE], dim=Nc), glue::glue("{prefix}w"), scode=glue::glue("vector[{prefix}Nc] {prefix}w;")) + - brms::stanvar(t(mix[2:(p+1),,drop=FALSE]), glue::glue("{prefix}m"), scode=glue::glue("vector[{prefix}p] {prefix}m[{prefix}Nc];")) + - brms::stanvar(Sigma, glue::glue("{prefix}sigma"), scode=glue::glue("matrix[{prefix}p, {prefix}p] {prefix}sigma[{prefix}Nc];")) + + brms::stanvar(t(mix[2:(p+1),,drop=FALSE]), glue::glue("{prefix}m"), scode=glue::glue("array[{prefix}Nc] vector[{prefix}p] {prefix}m;")) + + brms::stanvar(Sigma, glue::glue("{prefix}sigma"), scode=glue::glue("array[{prefix}Nc] matrix[{prefix}p, {prefix}p] {prefix}sigma;")) + brms::stanvar(scode=glue::glue(" -matrix[{{prefix}}p, {{prefix}}p] {{prefix}}sigma_L[{{prefix}}Nc]; +array[{{prefix}}Nc] matrix[{{prefix}}p, {{prefix}}p] {{prefix}}sigma_L; for (i in 1:{{prefix}}Nc) { {{prefix}}sigma_L[i] = cholesky_decompose({{prefix}}sigma[i]); }", .open="{{", .close="}}"), block="tdata") diff --git a/tests/testthat/helper-utils.R b/tests/testthat/helper-utils.R new file mode 100644 index 0000000..95db7f2 --- /dev/null +++ b/tests/testthat/helper-utils.R @@ -0,0 +1,10 @@ + +source_example <- function(example, env=parent.frame(), disable_plots=TRUE) { + ex_source <- readLines(system.file("examples", example, package="RBesT", mustWork=TRUE)) + if(disable_plots) { + ex_source <- grep("plot\\(", ex_source, value=TRUE, invert=TRUE) + } + suppressMessages(ex <- source(textConnection(ex_source), + local=env, echo=FALSE, verbose=FALSE)) + invisible(ex) +} diff --git a/tools/make-ds.R b/tools/make-ds.R index 1ceff44..38f9787 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 <- "87eef77" + pkg_sha <- "1f6edf8" if (gsub("\\$", "", pkg_sha) == "Format:%h") { pkg_sha <- system("git rev-parse --short HEAD", intern=TRUE) From 5b8363212c825cc581ef20311b6d38d3d682cd3b Mon Sep 17 00:00:00 2001 From: Sebastian Weber Date: Wed, 20 Nov 2024 15:58:01 +0100 Subject: [PATCH 3/9] adding missing examples --- inst/examples/ess.R | 55 ++++++++++++++++++++++++++++++++++++ inst/examples/mixcombine.R | 6 ++++ inst/examples/predict_gMAP.R | 33 ++++++++++++++++++++++ inst/examples/robustify.R | 22 +++++++++++++++ 4 files changed, 116 insertions(+) create mode 100644 inst/examples/ess.R create mode 100644 inst/examples/mixcombine.R create mode 100644 inst/examples/predict_gMAP.R create mode 100644 inst/examples/robustify.R diff --git a/inst/examples/ess.R b/inst/examples/ess.R new file mode 100644 index 0000000..61e051d --- /dev/null +++ b/inst/examples/ess.R @@ -0,0 +1,55 @@ +# 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 + + diff --git a/inst/examples/mixcombine.R b/inst/examples/mixcombine.R new file mode 100644 index 0000000..f2ca36f --- /dev/null +++ b/inst/examples/mixcombine.R @@ -0,0 +1,6 @@ +# 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)) diff --git a/inst/examples/predict_gMAP.R b/inst/examples/predict_gMAP.R new file mode 100644 index 0000000..56dc59c --- /dev/null +++ b/inst/examples/predict_gMAP.R @@ -0,0 +1,33 @@ +# 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 diff --git a/inst/examples/robustify.R b/inst/examples/robustify.R new file mode 100644 index 0000000..3b4c40e --- /dev/null +++ b/inst/examples/robustify.R @@ -0,0 +1,22 @@ +bmix <- mixbeta(inf1=c(0.2, 8, 3), inf2=c(0.8, 10, 2)) +plot(bmix) + +rbmix <- robustify(bmix, weight=0.1, mean=0.5) +rbmix +plot(rbmix) + + +gmnMix <- mixgamma(inf1=c(0.2, 2, 3), inf2=c(0.8, 2, 5), param="mn") +plot(gmnMix) + +rgmnMix <- robustify(gmnMix, weight=0.1, mean=2) +rgmnMix +plot(rgmnMix) + + +nm <- mixnorm(inf1=c(0.2, 0.5, 0.7), inf2=c(0.8, 2, 1), sigma=2) +plot(nm) + +rnMix <- robustify(nm, weight=0.1, mean=0, sigma=2) +rnMix +plot(rnMix) From dee0638553577d5d120f50ea1039a9c0cb301638 Mon Sep 17 00:00:00 2001 From: Sebastian Weber Date: Thu, 21 Nov 2024 11:33:06 +0100 Subject: [PATCH 4/9] more updates to the release --- Makefile | 2 +- NEWS.md | 7 ++++++- R/plot_gMAP.R | 2 +- inst/examples/ess.R | 2 -- inst/sbc/calibration.md5 | 6 +++--- inst/sbc/calibration.rds | Bin 9861 -> 9826 bytes inst/sbc/make_reference_rankhist.R | 4 ++-- inst/sbc/sbc_tools.R | 8 ++++++-- inst/stan/gMAP.stan | 32 ++++++++++++++++++----------- tools/make-ds.R | 2 +- 10 files changed, 40 insertions(+), 25 deletions(-) diff --git a/Makefile b/Makefile index 3692867..ad4eb77 100644 --- a/Makefile +++ b/Makefile @@ -79,7 +79,7 @@ src/package-binary: R/stanmodels.R echo "importFrom(rstan, sampling)" >> NAMESPACE echo "useDynLib($(RPKG), .registration = TRUE)" >> NAMESPACE install -d src - "${R_HOME}/bin/Rscript" -e 'pkgbuild::compile_dll()' + "${R_HOME}/bin/Rscript" -e 'pkgbuild::compile_dll(debug=FALSE)' touch src/package-binary man/package-doc: $(R_PKG_SRCS) $(BIN_OBJS) diff --git a/NEWS.md b/NEWS.md index 0d2baf7..b5fc06c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,12 +2,17 @@ ## Enhancements -* Upgraded testhat edition to version 3 * Added for `mixstanvar` automatic generation of distribution functions allowing truncated priors in `brms` with mixtures +* Slight speed increase for Stan model by more efficient construction + of likelihood. Also reducing object size of `gMAP` objects by + avoiding to store redundant variables. +* Upgraded testhat edition to version 3. ## Bugfixes +* Fix issue #18 of ESS ELIR aborting whenever one mixture component + has zero weight. * 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. diff --git a/R/plot_gMAP.R b/R/plot_gMAP.R index 73d4c5a..a9a9d83 100644 --- a/R/plot_gMAP.R +++ b/R/plot_gMAP.R @@ -11,7 +11,7 @@ #' #' @template plot-help #' -#' @return The function returns a list of \code{\link{ggplot}} +#' @return The function returns a list of \code{\link[ggplot2:ggplot]{ggplot}} #' objects. #' #' @method plot gMAP diff --git a/inst/examples/ess.R b/inst/examples/ess.R index 61e051d..da70031 100644 --- a/inst/examples/ess.R +++ b/inst/examples/ess.R @@ -24,14 +24,12 @@ 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) diff --git a/inst/sbc/calibration.md5 b/inst/sbc/calibration.md5 index e5527cc..90b6094 100644 --- a/inst/sbc/calibration.md5 +++ b/inst/sbc/calibration.md5 @@ -1,3 +1,3 @@ -Created: 2024-11-20 13:57:46 UTC -git hash: c17af53badaeb321d3de398b84d998dff5cf685f -MD5: 95c9b0f2f899285fb53b68663ec1b6e3 +Created: 2024-11-21 10:30:15 UTC +git hash: 5f5afb101f98533e4229f9dc267b3540b172c573 +MD5: c8f566f591c59887653bac53ab03bb45 diff --git a/inst/sbc/calibration.rds b/inst/sbc/calibration.rds index 970c6e14b20faaafbd46e918feee4a8416803c8f..2fc287a292727e1e7279c1888f1aedadfe94753a 100644 GIT binary patch literal 9826 zcmZ8_X*?9%8+Ii;B|H&AW}cF4A!HvjNvKhT>|2>fk`QBOW{L{gT4thT56Qk{H)G4x z5JF}wV`LlqSjG%ydp-a6{q&ws=f2Om?hog8&ULQ)0w(=+=zoK6Eh{qFgBUplOce7s zNG%_|v35DJbQLO7>Y}FXt9|pOK-kNXheMJ9Lm#Bo-y}X(jjwy;C$FsDHMaW`{BbAq z&trY300wt}GQ|e5bO~b!g zPSn!^#yq^g+{ElzY00`WbZMo*_7#LS`KIt}zB}(#wx@o!W-@{?rNJv|%zpS26+7># zmiOVq9u>j%0d4bBks`Ofzw|}xZ3m3^9^lFwCXEITA5@E^-JiM@-9)f&+IbgEHMoK) z(V~(2T6{Gen^YQ*!G5!UQ*hq+b>o(>%e1mBdo00mPsMY-l)e80@;fHC(A#2r{l7Q{ zQ~BRwWhR@&SReZ}R8vfF7*lceE`5Y4uR*;S3!mivZ^QpZ$2`aVzsdh#os|D=v#&qE zxT>fklxC(gv3Y*DWf%8SztHlE-(Q;@i2nsEhK*R~|1akk>aIbYVUmg1czIz--=gdO zW-`h{S6)9flTjZkIhUn26!BVJa_Rc({iMq);j??TuP?l;bNlqx%y5>vV=JR(>-l;i z>GGQHz#jL>3TS>Sozwq$Vi7I}HCUJlRw^@{d0Gi)}%ujUS~Zpp}h zoW8r~T)a4|TC}DxKtWGL|DuT)0{!4uo5FfPKyiQ_jy2(!e`^rdLrZI1=Wmc=dW1Fc z&x2~?1)wm4waBKqO-2j(1~K#(42C-cH(o6tE~JVpxa=5xZKC*NAq&j|(I?9Ajs5Q) zXf11MT9c8qx(VmU9FrWUc=y7{Ga;FH3@D#607HZ}<2I*iSY6{;n_n${IHop>T2{cu z)@nVX$UV;KrI}#U?8+MrQZ%_3o`AC8XuQvMFN&<>foHdV@yybEOq%NBV>Fh1!-3<*xXQ-BjWG|1Y;acD7?NT{ z1&3!`SOe&_?!J^^6%uJXAMTH6FW`I-Ss`ST;o0@x1Xd^kl|!e%02Vkrf6lXFbCjd| zK#A?%7T7ExITC4(2C1#9u{cVye)qo75#-PjcK;jOQMJrpyRUCDl>v=mTmz!Yv%>eui2ZQr0iui^{XhMYXrQnlCF< zZV-Z)UZ_fU@}wsj$OTQftMeoZTQfH~(~$$J^Aowj5Hm!zqystZR#VTf0&`FUo##Ez zP5T-?7Mfe|Q6)vb2j%difHd4i*oi%JAhY$2V@L#e#-mJblIxpe0FF zJ{d$ib@=4`^t>uGD|oWbNA%{t`>NIh?U|K)ThW9jf4CgE5P3pr=bwMa%X;rJF%vn# zS)qmg@yS=hlfo-h{MAD?y$9r3MHM63l8A3bUEKByos82vPU=lMp12Omf=Ftu<~g5Pv6P9#mgM4mNheoP&Q%=@wo-- z&>%5ywaam62=DRF(<{h4IPR6Y;I`065J`yy`;NZBR2bbdAR6V>(O!T?gNCE+f52Ld z8S2A`e%gr$k=TVbXtw7&{!#$`PDu{QtP3leWz`H^{537TvHM}$pJ05&y#Q+LYhS|3 zwQ$>P`QExT1nO?gQTyqNO62mbPt;H!T>WAF>XvpPP6PcwdX~>m-`lXhWvvTpRP6!f zt0<@#Nn3$NH%^4s0)DMMgtGpHmgJPF?8x=R8b3px-2OWK?A9e3w_~`<_}3=st;t(} zg04+?%Ro!2{cuuPe-mJ%_%H-N3YuxkQpY6jt7Cd@zw=#1J(`y<+Pgs%M_fxkwznNS zw19shSg@$N!*9&^EN1b+W;}j#Jy3P$+UOzHGY)b76wSC`B&*8kPtdB-CG6)F9fc|W z0rRl-xK}A5B-@Myd`DEjw0|trgWOn?ch}J!qK}x2m&OWRT8*EsC zdW823s2ipQoZ%R+syekxD1`CV>-^E8f&?xl`+)ju*CA1+a0~(}w0Cw>E;QTu0LDwG*s+*@Ckx^xH zCGzW{^%q$QJFyqzG+@Ll;k)`P(sDi|e^@3%N%Pvm{qGzX`|(c}XpkhCk1UuR_w`tS z=Bk!I6zlm;?r!nw-ZlYYxMEoCA8ia~95&oZ7Sq0_k>#))Xh%RKobX!9+LG`f7SvIy zIaPeA={*khg3)L{dh1ZO8P8jN`W_w~A)T>wb;+=@PTBlp{!}`KE#`7%;$GoLk5Ctf z@%MvaPQGsod@7V8+u=~TXF8sCTIApvo$Kjj+`TY?` z8S(wu8-?op0fu+S>B>Pl%`16wtEJc#u0RVd%s5m*!#R{yTA??Ugg7#FW)?gbJ}94U zPcWFX8Z=&ru_-C;$7G-F`g3;hu$zb*`%SQWb$=0>J)KX9?jVXoS4+EgzGl9-o0=@q zC*qzGtXjwg-WmU>3T=9-X1h(Yz!-YeJ#S0wD^d2}F-xO8Oazl$NHupYzYaw>GZ$}? z;NMgE6(Fh;^>^=opEy>Cpp5(U!n@^OH`}-9nfI|t-kSc30Pz5N)3Lj zv~Iqt^ui=CLQ_{H%ArpJbchTY?Gu7Fy{xda+oF1u|IU&YVqDwEp9aQ5Lyk`5jJP0l zdiHf1be6cOnCx2G{f{2DUZzyC#$qztv}YEsPi;=knA$ZlM2FFpiF_A4?8Z^rMM3Vq zqC}AT25Lj(UlX*^14Au~#iJv{!jW5>uWQVkCj)!5mg>jwyFd*U4vh4Mv($l$fzpmkgKCN2T!1``6U99SL#^8B1f`MNV@ds zw4j{?xgUq2o_}Fa(;tOYc@&CF&4NeK`qpFdw(Az@E6zJHci;AWK-eU@cE)cyk6`pA z5u!4!MI+Cf(b$(AwwY_toL*nz^R{Hj$?>OYZC~i76-t^?m&8I!rJ468*Wv?y0t~3s zEG%W6Rf+C7_bn+T78YoTny)bXBwJ4_H^X7OI> zSH0^G?UjJ^HshkXG^)FG1662p#!P`Bz$GqYe^09K@>Bua_QRDG-O(iK_a1Ozz+d-k zzMyxeuWZ;dWMWcutjUu}hUzWZLN2|s9oSy%Gkfot<6!YiLK}rOB0)}x5+X`OS^cPF zAw=6ka9%jTajPyVdqA8JnkIV4=!EfKu8*OF-7#`@UYGC~q}HQ6 z?SPe6+rv=TkH&hRM)D4Bmkfh=vM`VzWfXkM<2jJeb^4Qd^KRs@m$Nh~RygHRmW|hE z%j{m~;Tp##Y&Dt0x9BSe%}%v?_R#{I5|zB+y?SIL0UV` z`s5_Qxqn}FljH}@G9JCLnifV(6qR;sBpx&JyHOJOQw1SHq8777!Xn7DqghpjR!t}> z8tMA@G6T3R0&Mzs{2nbI_RiZ_fAwtIVU*(K$|lZ-p$mAyaFeB$gq`^HFq;z5Y{&iC zUU*nH0{Lt99N_NFMA)pNfa>9-BZ3)@J#rYtK;PJBEy?4+!-fMl6Gtmm5hpi>!iLwL z4f*KgqB9o-7EQXKJmB^d?6>SSo1XbxAK~;V({yP&Le^*T_wju%_i(FG)zX5Nj#Wl; z(nxi>nAiE{NygzcS=_zSrNlGsIMk3g0SWNpe+X@naD~-ht}H7J@&*C6bQae5hmD3L z%D?)(`YKf?`L@%F-_P}S&zW~Jkf`+XD?rf=SJYY*C5gL!OC3Y(px6!?L5U9Lc z&pE!9(Bxs;v*2#OE*PO!x2&|Iy+G9&_P^na%VRYdWZCV%sM^aHRkiy5Y|*0ZXBRYE z_j4WI>n5U*;pw<8$c9f#Eiy@M#!$qi(Jn`H%`06mkt0O?ipk=;uSm!N<=MV{4A^|3 zt>~+!vL@E4VUZ9I<>XeT*Ns+K;DEQ6PyOObK;nzYP}$t%lrs)JwY9(wqfJ%%gYq+X zD>^o$ty@J=->^zB0X%xs?R&a2q-7?|({4Yi2)%zn_xd2en`o940M*G!3~^SLeN{^) zkbFYOHDpnSh^LuU=O{k?@l;*ANINx_14Y5yRts$k%f78EvGN}x>+{iOnpSfv+$<~a zkKdjRjpl^&29*6r7HPH6zS@S%BL6Cr4r9F+QvzFx=1fB8n_k!q&zm)T{gOO;g(x>P zp0QaI$o(W1?%^lOAySItuA`0jd{jrb)Affd&y< zZ8)A%UeIZuqQ?;R`5j$Li{Nx(hD#)<E$EI&s zU>HMhD!p9}c<=n;8}C*=PpcL&TBqVQ(3P|yi`Xp4I&MH2f%RSm6V5miMCs%Y6P~xn zyA}%;k8WHVwPJIcwR*mgIp$K%gD!}ZlV_TrO8rB8y>y_&rKSw>N4g}@_1|ukrwnwx z^RVxcecJF-btf`{AFNfSU)Awh&?B^nzJ5#rKzisBCs?LbeQd*NqfZsOI{gfvRGEW| z3#+zh@2T>C#ZNjqeR=nOg2_GMKSaFMdGHbAv-=wwZjdnkW_rXD&0AqzAV4_v_=Qfc z8LpjMs>~GAUG4?I(nkpz4DPiOWBR#E&G#CF7%y`EB#5mz1x5R>)4aV;bV(P2^*(hU za*l&l@VWcmDUr{Y6=@cZ7}8W=JSl?|B+quI9~pg$IA*SU-*AS{1Q4PedAGs0Dt`{q zjlA13t(UMlW0|m$KQ!ugrw23L(DCAWrwapt~G10M=-CtJ0|<$KPnG z{A*M1o2OlBuouu#d~6wMXUF}9UwkbXZj!knG^A>BS3(E;Qf7AXbJ_maP?a+ycmsSb15R=Gvhx+chMUK8O$guCBb;dVi znspo1>df!YcFE8)3Miz*gHjWcuKEjGWwzAQ_B|s{Ek|zAn#EkNaeFd4)xA7LJ~)^X zoH`Wd$4w|%cV=R|+vK8suO+@tQ=w+!iqTGz{@D)e&#emG#Sk7e$Y|K$ImQueWiM7L z_%9S<5WMR*h$|A%NZQi0c85*@cBWMm5mF${QTvvyNeTVarynnNc-oGKOsD)B-7eg3 zABM^08Qp96lZouu(WG=2UsMfs+<3RcRb{mgn!Lj_i-#fIQ`&%a-T2e+)9{xU6It$1 zjyMNQm}MnvRa?LAB^B;A6k)i}YeXO2kxoAqv)M~9NN!eF)r?s%%61Rod{<;^{le8} z=1!DvnK~ZINj!D==*rc!kB-aUp%A=aGCA#Myl0-pl2IB$sNcSandZ$l>} z`UKogQP!e({5eMz<#mP6oNh5<+}PO3NXosps$C_2DaWKR=Cg9>l?b7Ffh|6qEEuqB zGLpuXT1nH-YGxgNmq1b6IL+P>g1W|s&48O{zu=8>M{bO=M2(&_URwQ=(fg{{L?dQQ zQ_36YG59v(xsvSm-z(=C3Y(|ACq*CqGu}xQpO1D&^(k^kKMVmQr5P^6*7F* z+MVX=2;&aUE<@eVWS(g|n~mqXghI1ibbfJh!072%WyYF;A*J5c+OCcq*Dz zY?+t?D-d6K7N{T5S-17^9Q*#ryqA?%TldGflWw7ia%dkQ^ONfj*#$J|>L(vEn4A1| zpMvuhn?28c(8f#XIPoDr`lOJ(W6#0Q52J=Eg-i$!`h)~o;i*HP~5 zDqjo$z@K(jrNt!+ewVMTB3*c@0~uj8o819VXNiYfiR*TRDF<|?P!-xQ&>O0`)DxYj_(S7>IL2l zWfmC|G?<0x2PfkagF0wJqSAM=nI1MhI@dg76-O}pN2B<5Ztfy3VMTcEsXiS7$w3!Y z2CT}Nz=(%^{QJ6W_z&d;&AFZA9Hz3j7B-g@-W?*|AbmN&erCrcJY@H9VT{Q@J0<-3 z-Yv11lx{y*bnMALwGYExyZwmhhB#p7^nm`vX*>CBjML3a zJXHcb!x7`8ioM4c7z7iZ-%N-E6@yG)ctDDhm(CZ-6Q1>3ozP&%bNW6c=z0u}*9OL; zDX|yvurStM6z?Q2ao{|6-7$ObgrVS1Q|bbU5ASO7`@0x?u#&D&U@F!XZY3Rk6nyN3rkjlA zh)|^xt0mU+S=$loFZA}W6nM9F?0^7JhymIx-vmYYlm6;Lj+V1~tiJosfvC5BQC^F7 z`&g?~%RL%AJnxJF&olM=RdGAr(qD%OG^G(V_GvGHR9^!atU$ZZ>fLH%ay1xAyQo`{ zqcRs289+m0#zdW%>T&RR9!NI%XQ=Is=*L|o;{5iVw}je5cW(_7Nb_()ZFBebH%eK7 z(y)t4(izy#er23o*N2y?KDuA_2DW`owhtKhwH<*La8ypdr@jbXY%=UU&l?_|*Am18 z|1;%o(}T>{cs#1~X!j&hjwa}ts&hpBcLBTS2$GLJxH!84B?p>B3iv~n`8LjqhwB!z zi>4rm|+#dSM}k;@Cr0J7eCOBnZQu5 zBmikaM$>N*>ZYK-Ng6&9(|bOXw5|8aO53gY%@1N4si7XYWBX(kNC-DWXy&M?Mthp| z&D*P?tGtc9e8r=cA};a53AA@Q+p#NU$~cpR@~0Z_;TIPDR*45kSy+WVqJ5NUTw2FB zD7@fKBcgp4q%2`av`@5>va-B{7_&iJ@@)ijuuL7>Do(nDJDyMbNIW6es^(Sv4S^}|-#S>3ZfhX;kJ z=z5E?lvzUHA(-P6pck(FZGz%m;i?6%i;MZ)=l^?^hBY;3YBkG2(eY%-I+2OoisXSW zdw>H=xGvw#w3~k#=WDl&`KUH$T=m(>^nk{V4KMhd9+Vn9hRyFz7QkqoA2;e6Ms+1nbjQvcsVC#V6VWCtKimBEgTY_~GvL>>cs? z{)+tZr-C_uOHOUs?yF5;%{znm=9wrd4A+9*iL?)Ty4fZ27uW;bAOR0Gmv(K+uPcWg zG>&nyGc+-_Pn<3G}fn%VT5!?v#NtC-Va-(&N*Kn+l=TuV3 z!Z;nF-YrOqybIHboJ&%y2eawY(Evll9rAg&|f&W z@I+f#s6vPXY<(92VcXUcoT?Y4*a@e@F8PcKBmR86o--ukey2zE4AytKjfpy$`P^D^ ztSO(3nmU)(Cjgh9Qry}lx?(7*$3AB4z5x@H%0=JBBd79`B~rDdFq$8Qkg0)++zV5) z=+Pd25hYIN^x$G;cfPRinloeh<71@Vxd+nG7ydQCqH}!*KfCWC9zDL}RUbD+h0Ys< zKetYRn^R+e+3oM2WMozedriebMt)DPozyLl~87c;le zbW{;L_14CrkYKO%hgd1r=~e)!+CK5O7$r$va)qz$Oj)+mzpIe`#7>H-Rh${hW{*5A zPsp$z?d^rG>dy&A+-uq5wfs3_IWKzM=+LS<1z@)$I(%D%d#2g^M%sW0ViUG!5l}Gd~T0y60BJ+m3)*Jf(A9ST`2XKf+BB zdKR6oGlF5_o2bQpcPQjA64fslw}+RQIxZl-@hr|{SV`)v{reYVm@kDv_rH*;i0;ydLk~?X={foFX>}%K<%c+@#22xR>D0t4Ccen z{2EBmGpCt*bXCICo6(d!!C~WNzY$^H{ac%A8I4xTIPaUuS99JB+)N*znK(ObOA_cd zyh*vZ8SvDzkJ-588ex*YkyWRRTl+ia@|@rI+w8sWSg%>}v^Z3%lL_0X80sIi@W@2E z*W%IfuyIQQ&8SDFz2YA%D}mMNQm)LA5hB#`-BayT!YQ=3bzyN+sP2GxImtTtBB>bg z;`0IFGxZf;DSsbnICkQ*XNmk1oOpU%O2k!Z_=$H3{rIM#X+ce1DH2izrl0JWJDRh! z!7@hbO~r+01#vv>dp>mH6NFu9xJuQfPJ`nMko9k8;sgGsAJ4vR*Yo>HIqbkA#YxE4 z{4@Id{iN4k$Glm^JKb2VY%2E!a0#& z%qV34!eUq{yNYin-1T*U<7|(dS##}m*(5slTpR9sxWp2l*?z@2M3s#LP-Q4SDJ1+g z*X=bo<~7*%6HR%DDY>L^YN^0%-|02B@;a*eI>qx7P-#e1a%t9V-_oqP_&VgyC!qRJ z+NmXrBs$D&fAbvY+Us?-A#2H{;+ORFB>J$~{y{=@$WwWUV-}Tr{Q#PG&e-jf#Y_5+ zbH*;8@Gt2RuaQO9A+DcrE}w7?pK?@f?R2;g*}r#7HyvjydpuKcGN}I&D}L!rzdo91#qB+T(sfHh{YG-KuO1Z zGrRoEgCGqw60wz~S6Q&)&SHQNz`B%ga8}1mPQ%oGQAN&^a-&Swj6(s)e=ba86`2zA zNs|pPz4uR6eru^Lx3EL>rjnz>)*ejlWi8}91P>0?6htz^GfondEMqoWBzDQ4wkKs|5wD>qL!ozJr&d-5_V)E+5V*?q)gr};XivY&z%=U5oQD&v zZ6J?Zq$iri{svq_G%uxPgSEG>`9}`}futQv6fcl_br{I~j)ycax$#D}U{>S9(*k%KlMdTnDni>Uoc7~d@veR++yNCCC@4rlS^#C_JsgB}d0 z@%pl>X>Ck=um^=pUI-~*`xlPncBv2Njd4)l9T#i2^<;=@%rZF*`uGerhxKrFUz381 zR+H1!VeJJ@*4BgQl(4EGjsz9BxpUYPvbJKU&Qqq+ir?Rjke#xHQ7MPhq90PS`Jx*F zCtBBTFvwgrBkhK~@8kqIbMFQVZ-t8eQJ~&B#ax~eYMK`cYpOSYLuNEg@35Gwf$Gh* zvqBS0Ra?#qKVLfW6~;Vpk?r&*CV#WP{lsnU-T%_qV)Lzk zAwPDDBWh8d6C}tuGfvZv9ddHQw6b zxB?>R9pz~tm&aAY#O1v`J?`%n-hBB$KufB5Lv4TdZeJyGzt(k3TDH>L@){!ZNLjm) zS09%pulBQru|A2dBz^j6y*4#mgTg!1&-mjAK->l+mQsz{re?Z=t$SUW{()0CglpXn P?oh}zB?DWj!-xJ43P`=# 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 diff --git a/inst/sbc/make_reference_rankhist.R b/inst/sbc/make_reference_rankhist.R index b4c57ab..20fc4e3 100755 --- a/inst/sbc/make_reference_rankhist.R +++ b/inst/sbc/make_reference_rankhist.R @@ -23,7 +23,7 @@ library(data.table) library(knitr) sbc_tools <- new.env() source(here("inst", "sbc", "sbc_tools.R"), local=sbc_tools) -set.seed(453453) +set.seed(45348346) sbc_tools$rbest_lib_dir <- rbest_lib_dir sbc_tools$rbest_source_dir <- rbest_source_dir @@ -80,7 +80,7 @@ num_simulations <- nrow(scenarios) cat("Total number of jobs to dispatch:", num_simulations, "\n") RNGkind("L'Ecuyer-CMRG") -set.seed(56969) +set.seed(269698974) rng_seeds <- sbc_tools$setup_lecuyer_seeds(.Random.seed, num_simulations) sim_result <- Q_rows(scenarios, sbc_tools$run_sbc_case, const=list(base_scenarios=base_scenarios, seeds=rng_seeds), export=as.list(sbc_tools), n_jobs=n_jobs, pkgs=pkg) diff --git a/inst/sbc/sbc_tools.R b/inst/sbc/sbc_tools.R index c2b3cde..db65eb4 100644 --- a/inst/sbc/sbc_tools.R +++ b/inst/sbc/sbc_tools.R @@ -87,12 +87,16 @@ fit_rbest <- function(fake, draw, draw_theta, family, prior_mean_mu, prior_sd_mu gaussian=cbind(y, y_se) ~ 1 | group) options(RBesT.MC.warmup=2000, RBesT.MC.iter=4000, RBesT.MC.thin=1, RBesT.MC.init=0.1, - RBesT.MC.control=list(adapt_delta=0.999, stepsize=0.01, max_treedepth=10, + RBesT.MC.control=list(##adapt_delta=0.999, ## 2024-11-21 lowered to 0.95 for the sake of performance + adapt_delta=0.95, + stepsize=0.01, max_treedepth=10, adapt_init_buffer=100, adapt_term_buffer=300)) if(Ng > 5) { ## for the dense case we can be a bit less aggressive with the sampling tuning parameters - options(RBesT.MC.control=list(adapt_delta=0.95, stepsize=0.01, max_treedepth=10, + options(RBesT.MC.control=list(##adapt_delta=0.95, + adapt_delta=0.90, + stepsize=0.01, max_treedepth=10, adapt_init_buffer=100, adapt_term_buffer=300)) } fit <- gMAP(model, data=fake, family=family, diff --git a/inst/stan/gMAP.stan b/inst/stan/gMAP.stan index 597642f..b9bcc4d 100644 --- a/inst/stan/gMAP.stan +++ b/inst/stan/gMAP.stan @@ -134,10 +134,8 @@ parameters { } transformed parameters { vector[H] theta; - vector[n_groups] eta; vector[mX] beta; vector[n_tau_strata] tau; - vector[n_groups] tau_group; beta = beta_raw_guess[1] + beta_raw_guess[2] .* beta_raw; @@ -147,14 +145,24 @@ transformed parameters { else tau = exp(tau_raw_guess[1] + tau_raw_guess[2] * tau_raw); - tau_group = tau[tau_strata_gindex]; - - if (ncp) // NCP - eta = xi_eta .* tau_group; - else // CP places overall intercept into random effect - eta = beta_raw_guess[1,1] + beta_raw_guess[2,1] * xi_eta; - - theta = X_param * beta + eta[group_index]; + // expand random effect to groups in loop for performance reasons + if(ncp) { + if(n_tau_strata == 1) { + // most common case of just one stratum which simplifies things + // and in ncp mode + for(h in 1:H) { + theta[h] = X_param[h] * beta + xi_eta[ group_index[h] ] * tau[1]; + } + } else { + for(h in 1:H) { + theta[h] = X_param[h] * beta + xi_eta[ group_index[h] ] * tau[ tau_strata_gindex[ group_index[h] ] ]; + } + } + } else { + for(h in 1:H) { + theta[h] = X_param[h] * beta + beta_raw_guess[1,1] + beta_raw_guess[2,1] * xi_eta[ group_index[h] ]; + } + } } model { if (ncp) { @@ -163,8 +171,8 @@ model { if(re_dist == 1) xi_eta ~ student_t(re_dist_t_df, 0, 1); } else { // random effect distribution - if(re_dist == 0) xi_eta ~ normal( (beta[1] - beta_raw_guess[1,1])/beta_raw_guess[2,1], tau_group / beta_raw_guess[2,1]); - if(re_dist == 1) xi_eta ~ student_t(re_dist_t_df, (beta[1] - beta_raw_guess[1,1])/beta_raw_guess[2,1], tau_group / beta_raw_guess[2,1]); + if(re_dist == 0) xi_eta ~ normal( (beta[1] - beta_raw_guess[1,1])/beta_raw_guess[2,1], tau[tau_strata_gindex] / beta_raw_guess[2,1]); + if(re_dist == 1) xi_eta ~ student_t(re_dist_t_df, (beta[1] - beta_raw_guess[1,1])/beta_raw_guess[2,1], tau[tau_strata_gindex] / beta_raw_guess[2,1]); } // assign priors to coefficients diff --git a/tools/make-ds.R b/tools/make-ds.R index 38f9787..200cf30 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 <- "1f6edf8" + pkg_sha <- "a0c573a" if (gsub("\\$", "", pkg_sha) == "Format:%h") { pkg_sha <- system("git rev-parse --short HEAD", intern=TRUE) From 6d7eb25d7d8c86a209554648d451b4ac6b40e89c Mon Sep 17 00:00:00 2001 From: Sebastian Weber Date: Thu, 21 Nov 2024 12:15:12 +0100 Subject: [PATCH 5/9] regenerate doc --- DESCRIPTION | 2 +- NAMESPACE | 1 + NEWS.md | 4 ++-- R/RBesT-package.R | 1 + R/sysdata.rda | Bin 6091 -> 6069 bytes data/AS.rda | Bin 238 -> 241 bytes data/colitis.rda | Bin 233 -> 228 bytes data/crohn.rda | Bin 257 -> 257 bytes data/transplant.rda | Bin 267 -> 268 bytes man/RBesT-package.Rd | 1 + man/ess.Rd | 2 -- man/integrate_density_log.Rd | 2 +- man/mixbeta.Rd | 4 ++-- man/mixcombine.Rd | 5 ++--- man/mixgamma.Rd | 4 ++-- man/mixmvnorm.Rd | 4 ++-- man/mixnorm.Rd | 4 ++-- man/mixplot.Rd | 6 +++--- man/oc1S.Rd | 2 +- man/oc2S.Rd | 2 +- man/plot.gMAP.Rd | 2 +- man/pos1S.Rd | 2 +- man/pos2S.Rd | 2 +- man/predict.gMAP.Rd | 1 - tools/make-ds.R | 2 +- 25 files changed, 26 insertions(+), 27 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 34bc702..d4ead43 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -9,7 +9,7 @@ Description: Tool-set to support Bayesian evidence synthesis. This and Schmidli et al. (2014) explain details on the methodology. Version: 1.7-4 -Date: 2024-11-20 +Date: 2024-11-21 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") diff --git a/NAMESPACE b/NAMESPACE index 377f92c..103bfa0 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -76,6 +76,7 @@ S3method(print,betaBinomialMix) S3method(print,betaMix) S3method(print,decision1S) S3method(print,decision2S) +S3method(print,dlink) S3method(print,gMAP) S3method(print,gMAPpred) S3method(print,gMAPsummary) diff --git a/NEWS.md b/NEWS.md index b5fc06c..bcc1ea5 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,4 @@ -# RBesT 1.7-4 - November 20th, 2024 +# RBesT 1.7-4 - November 21st, 2024 ## Enhancements @@ -7,7 +7,7 @@ * Slight speed increase for Stan model by more efficient construction of likelihood. Also reducing object size of `gMAP` objects by avoiding to store redundant variables. -* Upgraded testhat edition to version 3. +* Upgraded `testthat` edition to version 3. ## Bugfixes diff --git a/R/RBesT-package.R b/R/RBesT-package.R index d128d85..14cc304 100644 --- a/R/RBesT-package.R +++ b/R/RBesT-package.R @@ -38,6 +38,7 @@ #' \tab \code{rel.tol=.Machine$double.eps^0.25,} \tab \cr #' \tab \code{abs.tol=.Machine$double.eps^0.25,} \tab \cr #' \tab \code{subdivisions=1E3)} \tab \cr +#' \code{RBesT.integrate_prob_eps} \tab \code{1E-6} \tab probability mass left out from tails if integration needs to be restricted in range \cr #' } #' #' @section Version History: diff --git a/R/sysdata.rda b/R/sysdata.rda index b219798bbefeac697afe55983e477f116dc179f9..2064e02ddf83cd5c2c67f24f3d0315fc846a5522 100644 GIT binary patch literal 6069 zcmV;m7fR?tT4*^jL0KkKS^fi;t^iLg|NsC0|NsC0|NrzN*h;_u-~a#r|NsC0|Ns6R zfFM8uKnvg$e^3iOQu8QmR;E>KlQU&)S+yu5QX11qRU0hXX{^@5LNmYs02HZ-wK5M) zG@1YbkN^WiL7)Rb000d%84V2?14fzyAoVn9h|mB400Szg{$y#OG-wkHgFqkzX)+3W zjXfrr13&-(0B8UJ00001qeCDEr~qYD|Jn?IOc^v}OieOigFS+B_(i&-_AZXJ@O*V<9&>#&OX`ll@4XLKodV|yqhMrSodZyI$(a=T8B~evW zRHCU=Qn}zKLP;P(1OW&NCnrK60uTjr5H&5+Av*O~>;s|&OI54YBU5>EdlnmJjhG>S z-{b!8FKy)6!_4QHv&tKJBp6B|J>dkEu5n?2}9ERs|D?PPmyw_XPJO z_eAT-?#b8_-V?DWy(hjWWKPVP-u#~cp9G%>pA4O$J|R9SFm|N)#O;ak$?^%?lej0z zCvs1ePUxIccVzDg^2y$l<`cdr7@he%0(d0wiP?6&r)_4>H@^O#0KN{6=9`*WijaIU z-;;Y!Z42HMkpiFAX~W_N!n5P58O!QVeyg$>MJ{h7AVg1i(2uxpO9l7bk z7IhI+s@OOdTg65!z6PP9#a6TS{1=?@qis`Y@iF0|J_?-!)w@P5M}qb)U=K!**MEZr zZK1f`b<^&{He902F4nn((fB1yz+yvBwwEv@55ckhKi9)MMdc@H@NWd|&1{8nx=1v2%Rm>$Uyv;a4 zLnJl0xza*nXbP%iS~2rgU<>?_6HnBfNNy-&fnKG$BwmxG>%w)Tm5^Jl=&s}PSqZ68 zh0EPm_=LOQqW%0CNZb_DT>0OYW*S7i&7d%9hT0W#N@7D=knqY38?QI=c!V^1-Pfgi zr{{7)d3p1f6!9X4Q6p}Ugf8}Bg_#XvZFCPqFa~&^%yYu}R)Lhn`E(7=DgaQRLZu?> z>vvw5-HTk7+^ku%1qtfX!jZztbRSHYF18k-S?G<9T+U8ygrulWx-=NfHsWBgg(x+^|*bZ5(fq`W;OpjAE~T1CI8Z zaTO_(YP?g|d^)x;-Vxh1DUO$o3Nq*3>qWMr)q7K92?KWdmug*#W5qNwq2h8(B`ulc zLnHzvFcWbaFlk&2;)_$n3_pG>L!BHg(*3=#;JX<`uH>A}>a4VmL}@9_mWyf=hKs2y z#OI+8=Y|GjQ!wE5<&@8r5J}M5Q*}=W)ZVs3xN1XOQ&xMBP65A=+O0+5*shyrkGI%H zdDfog0tA3}s>RR{|y*v@5GO^smwK-!Ni7YL`xALqgwRhbvL z5mCu-_R?4U4l*c0 zAQe;x3rpcMHJ;g=;@gi&T5YM|xw#6T1FxyB4XQ1ZO9Ei@Anc&Cl1e$wKj9Ozt#hH< zk9764vntU{V;ML|&px`^trm|66i}<)gzG_#`qb!T>1nx_wL?^abq*>+)bRn9Im|(p z@QUEckCxn~lF}148qONP(J!Z5v0EYSv0*{Hs>h8L2PAfSYtpGm zTddL79qF-Qv~H#T$4QNHb^!XAnP2FGH!~cDq!4U|VzR`A3t>xD>o5b%q@dzby6le@;rq|&=wV8@e*N!Y(Z6@J=q1HzM5&3OUAKz1{KfzI2qyy^`^(K|a}uz5t7u?rE3 z7UvMz&2V|L`n|0$0k(5lt9sKlak3lA+s&5f1Qcx}I- zW}+MjXv5*za;{@AFxlI;bizImHpdG~XCd&aJ0c33L!OIPdvV<4&<)s?E0oS%(>oy_ zWjiwx<%BxmJHQ%N@eHjvijX-9E$tBBHV&F$y08i1!rhZ|Ek^Bu0alDRgRniR>aitu zE{ijHmFJA(g!`;^W0cFW6wsvOG8uQ}7f^cJE2|MxB*CP##-jJ@dZ1<}XN$X8sTRGC zE~f@Gw~p?H2GtP}8xD&0$Nq0s=N@LW^E89LDpZ;~BUsT)NV{v5V7T?xu6>G)DYeWT zQinnT#WO{5RT?%vRGQP<1VnT&=CQZ&jYcDVr`Z}k5Xg+>zWQ1e^|!x}j?SXLQo|)m zJD*54(`$p|5cQX6RqQl$5Zs$fwLLJ{%RPFbZ&DJbcc|(mL8?*9ONTB?*DYa?{5M^i zP$d_2%yy-+s(QzheL5>9=i8-03ZNMgj1JNB~fwIdp z!PN6?%$kQK&`F6@DBIGOo)!3>h-zFWk*RY$qd4^DHA@W`I}vB=fJyc0C1TUQIvdQ6 z3skgag%l3JCufL(r2Aw9S?GLB)`S{5=%%BSJ&O_|*CMVn(l zyYHv`Qnpdgdb`@wh$U6qla2*yf;D4>p5YpLKuZ(BQQ5YQbj7*iFj!O^eGKact})>b z8EsL<87y#btgz7F(Av_5CK6N+k)Er4Cu9rouk?dDiA0k(S!qoq3g1%WgJx z-i5?k+%Fw&VyJ9gRAVH`^m99>7R(*P%JQhiXmVK7hh?LlzUD0_JJZEoYdp9i-t6zK ztSX2+RT+Rr9#ISIjZ7TEyA0-r5Ts4;wadcz-Gz!97Tm&ac)AAj>lJ1@?suP`q0MHl z=*W2MB>Spkd*g+4GdgH&isa4*G`rNcJsKw1CoNq_b}hQf6Yuz9@KxMJ(0tn4F5v7M zh~r+u>4h2NsBg(Zw)fY23+uv?>Qg``k|k&~^J-2zdr8+;i9i;I#Vm?`W&vChl`bv> z5|XVIjH{vQflBRTK;n~xt4ym*k~?5hs50rSAeJ6m9i6>Lk0ipIwSeM^@JfNj!MC-! zWJZ|zwLH3pOA@P^<1;QiFTM#)Mw(VRs&_cSj6}+=K-H@X89Ya1F=)L5*k`(i2P6x* zh@x%wh!w{LMrwV!Wu*Q!LFk<+y61_$=xE#^DvgJWpDJN&OR>7@w}LE$za=^@4?5M; zQ>()t33U*#S0DtzflLsv&a0JUWg>`UZoNEj*!tM`@@zp|HPfZf^h+)3vc%0>3Uo5E zpcK`%xEdUjLb)_J(9pxTTJrKg4auu`EY~{u6|H)zXf5EP+o=>=Jo8aas+a ziEXEOE1ycRIOf+h(l=mDh5L(1woL_6Um+20<$F@ltk;XW5kf8Wh{L=*BR*^j?|O0492q2aCv+Q9Gdf#MqEN9t1r_=)6lio z22j4m31oPCZ-vH{JS`qg=ci%Kg=x|wij_)CpyIjgc=xvhr2u#5Yu(wF#_rXL49XKK z?+jKYD`5dg5{}$QR_eLo$tQ7%^&Tp2S4$EQBvS>eo=Br!r;1=F zB?5OP9*j-i?pMFu>(j`FZha6LnJ-5O)zutA!QF=MMrmNpvB?^R)a-R#inUL1tqehB zs3w~+5VCF~d8yztre3N_xSBL>h>EDnN{b(!$bj-EIEHwL`to&HAdR%O09v~*#zqFU z{QYa_uBJ`aov=rWAXpcwSJR`yf-B;v>S{X_RsV4 zhRC6IyY25`ycw6?&dSl-rCbot6kO}i;)cr$G+UN3a3#W&01j3kG)UqCDvBswVw0x8 zQwJgWD^C4S*OFGS(R1O_=H6?@okzat2Jzm0rs8$_Tv?k5w{BoQwKw#%EYj6VC6u_m z9D@~B(8F)5vZ^C5CSlm!a@wr~Q;ES8ll4W5Q_V5c+ZEe&2u;U(550lm4vQAM41Awz z=ow>N53q@Fay#3dS-387_@L|00<8QCh> z+C`a60oWfWun;|*dV8VAmsPY@S_qEw%94F?&idOHYd$UQn4P>jJWCJS=2?ru39NBo zo#(N(+~X9Q-$|s5t0&!*C}D_rakFzUi!9IX%e#i!vKMk<6*CG*(365N_6J4m2d!(wt;mHAWv&*vDd=WwQqtf z+}D!bWJZg%oj1zPP_X(5FCxz5pp+3KHjdO?rx--uST9b&H|EsQz)az(1iU#H47>E) z6qRX=cD>f@^*FIP@yJa*hXx)zxhkfpR=JhIJuy}jP^R`+?LaieQn#;o+8Wv?xNl)> zY)ZF-1Gu4S=6*E{_Y@+2Rkm?R)`}njH06|CES`-N6~O?B_z%j*Mp*XaN0{69m z9k?jdP6?^S)Ey4|?=_=^s}xE)-r!vc&voKNIuA^#$5k)Ow@@Uo<);1awPeJ_Dj;{4 zOs@!;vC=@00v7=j6|6W|cJ?8FLD<|3G?1^s^*m5dmBB^`M??z+#w?>M>yTl@PcfKF z8a;N+CA9K}Kr;&bPL!mBS$e^Dj95376F?Kqc@1rMrOm;|K*HvXtr-qM2!{+5U4qnA zRij5lg53qjB2Ra|Xfa6L?OVXo^SZA7tu6QLkUIA+0hMkg?oQD5uBp(I0<)W*m$*ewxLkrxM~k?2%&tSVgXZT2|6poa1N93?CxsTpFtvzTH8sV7ZiK7fp{ta zv8ga(yat_TfpYqT15dSM1Hg&88=EEYoG{W#QfiFvvk>?>B685H6H!x6@`h>>=5rbR z5ae-4v-YRX&{GCFZ_{2DcgK5t14FG|6^Sn%#w9ECUh(@kvRofr5P+KOjjNvOrhSnU zU}tpYoJ_n#s(kpGODGdaZSL(&_*m~QCz{)nY`!v0sH%;=BaFLJT#4j|lUKe}X{QTr zVqK|Rs;fb~4U=Kytzs1jPdBzVRjr}rb^wBPs0djV1J_E@y0fexDzn-XVHyFgmc>*= zt;YINl)8h6D4|iY!ajzTS!$xM1JovVyNvzgw@pF%Gigz$bTo_{FexluWG~BB#ei9( z21QK`!vJ11j8BD;us|Vt`vJ-`A%bO;%J3Aio;+Pn@zftqgDg5ClQsj4LQ-s=&zU zM@e1VTu&#KbPmP$;{$62@D{~4&>$rmx+sNo@b2zju{iNpv~LA4ihM3LE;m!Yi)6P83D=_5<_c-Miq~|!I`;^%KGwGJia85VGFg> z#8-o#R0^h`cs!2Dba`1>sa^c4M=EZv5=2vo?v#oKrWUvWX$&0A;>@((VM6v%uNq}i zAPhQeVg(0tWCd3C*-H_1^6@BwNqLyc0rkMHt|_5mZJdT00~^5QAh;?mCa+BI4=uEU z=C#_XM{KjwC3}v|zoH^-%abJvDkmM9A}i6VsbF3e1w3H3bCO5gYO|9(is-GE*jxf0 zO?-M9e#SNTcppz%H{+&cqfcAm*vxcQB^4TZLy{1^lG1Htquxv9n_it-sPE~HN9?(W zre=h$$39!aGo959J4f`6M`L8%pzKW5++@|j;@v?UMRH6uMYK(GHTS{6sepxMUsNu9Kf~yR5Z~! zsWL;&SKlddc97~4^5{Vp*EKpUyas`Tofoy92)qYdUYxqFZxx>5>lFffR zxxR-a2oA0T~YUZT^ig6Fxz`waNvnG3Axu_ zk+YUAB_QnWZ@Xm+kQ>U9()P1#Nrn_eP`m{fEf0p&Of`kQmvbuAqFUM|&MW!6icVzV z-Xu6XDPQO5O@a^&wi}&rNn{h5(!!qlDP#xLx;cV`M_}UlANQ}k579sq_<1obC+2>; zBzPEdt%=D~mB;j|+R4t@xsgajh$^8ZCohqE&;Ge0{US*gXt>>-z$PZu_a*?r&q*Yn z6{#eW$DM8GNoHg%jHwuEAcX$%5G4|zsEHA!K_Ws~WQ!oON+lH}kV|SLAuORJf^%q} zUXdCV3#$60g|w0{+hP$KnjC78Xr4CiylY8sI;+L;?j~*jM>{LWzDnt_U1JXN`7PDI ze`@hO{J*uC_VMLgRf5kW2{LA8-@|yDYuN8zyiK*$b~}lG9bIM{h;!;8rzDa{+{dA{ zuxj%TcP-3yIChR^uG5!I>OE$zj<-#z((7w$dR1a&Ya$0gMAc0@hhGtpPRkpbQ$@IxGyE6ar@A ztqokn%|wzZ5=f-Kb?^n?65Q|%;0%rPErge2%2F#9-pNR&3h66wN&vK3QK?A}$x$6` vLL7T^C{kap(lO5yPn6G5VKX>WjG|oG3_t0n66_NVj0w literal 6091 zcmV;+7c}TXT4*^jL0KkKS?IXd%>Yj=fB*mg|NsC0|NrzN*h;_e-~a#r|NsC0|Ns8L z03<*pKnvg$e*in(c^1u4C0aD9DOxo^`(i*8>=`8zXwr(-4%t!;_y7O`s-K}e5cLe2 zVKl{0WXQmpX{m)89-tZk0000002%-Q00006fB*w3ssHQ{Gy$L*009PoniEY_{i=GW zsNSRv0000000D>q000000QCWmnrWsdq5|9T3l_bn32%N84PcIX9&9kzp`)d|Wu{l=5a7ukl&;8Y!Hv1xN{?%nDp3t7u zmV z@X6j2;uGSNy(fH6j8BeEK|Atz1o;H;N#PTcPm)iRPYj*k>+ZXb_rBwCd(J4N-=YKwT%(@{4Djxxe(DjRd2JkCDT7(t$ivy$pGb94V<%i8_h#m$iF(?iN zjRpvaNa>$A4N9M~88rLcz^DN7tiStzatFjR&R5E0YSv81^D@&W@Ev&9S$_6+g@d#4j&EJ>sK^aE5naGXWY_Kb1A6aHDH$4;?GHEGQ1YP^-;y*C{i?FssiEV zI4bhv(hDL@f!S>kV?TFf!gt{LRWMV?NCuiX;*S!gqOB%p!)h4G;%^xnj8qnVrKB+| z`Fe2v%jQYTt^XHDdy4ZT2E`IBmk(yjB(p&zGtWxia$CDsdtu|yXPU*=TQ#pK!+56S zIr^k9u`cFnFwWl1#f|s^(w2t+D?AS*j9CIrtEOjrB$g?|-pxBo(Kwk}FVL{d5Q~ce zOSnvI#7NfQB@NOHEVv~D!xOf?riae6%^D}gJepad`_yKaCal9LoISyD5{{PGG+>{8 zRTYs?8gpRU(Lkb$5i7Sw-qZseSNE5+lv1s!LqL1$F&>AzGj$*|FVfYtv_(@CFvQ|R zl>REr5+Vvq5?g2+w6i1QT|D?ebyd+F3V?JqP=_hBSBZ&EDMZ*%I_{nblgeZqjO7!=aX5@>CNWx)c*lI%DkO z_@_!x4FyXcw%%IqLKjxjVQuWTUAhxg=eU?u#sRT1E74IELRBxI86p{@eGbynwm2dcuqMXj~HMkah$*lZtEF08cJ4P(hm&gU}EvWpV&jv%AVYSJ?3RsLRMI#pLE zi_~0EJ6|Pc>G9os)bldVwuSX5@%Qvr6$ZM&hb*MJ5r|wI+GGP%}~5dK!S%y zvsS7ALA`o&2(W~Ty+XOu&50@Zw$d~+-EdiEmFLl}?{q^qU8Y%ZxWXvvdv(=cA#K?V z9Q$ePPBlSUx0rihXdFdHln(ZL+V*R&F;$UKD${f^Wg5y&0l7#!G3RpA9}}e@4BMU4 zMZwcdux@8d1UaKBl%wYziJ+Q}$Ix*+vl0@#nIPy*1H`lWHLQ<(iSz2f!b-I~y(4a3 z8wE-Y^;mrTn5l^|mvD_F*^E%C2TNICb|sl59GhzMhKv||&a{biPswJNL>V(u$q#?^ zYHKN{ zY~pOj{J?WpY)#pto>#3XdA?l>ps+a5-as;%8;O(PshC=Z{$i1e2Mc7PoKbKwqHmL~ zNtNgZRf3RLVOzmfPKnq*IP||V(!$YYcPAZ0#Ez_rgdWh*?3NJN?j?MD48-nWJuuOD)q_OFDeXbE^Jb5FN5?UDwGm3y&2b4J8uVmue9NCTGC2>Z79mmZ5$w-br=UEUt`jKFe>{#1m9CxDk@V0rMSj_}Qp6oy$(>^AQK$wxaka0KMPLkj4OUNXP3bcVU9pz#maV)O*>Xur z$2qaM-~+lHt_q~uY{sE-kSqZZ^x4l|I7iz^<>Rwq=Cm7z4fg11P>AIDX22OYX*u(% zaw*A-3FI964JTUJHkEDy(&?C54nb8)eCKRWs7mGWq-~?I1=<}1-?^kvf`2&<@`o!! zWiOSp!K6jr#z1LPsdLxN^>!`F$zl)v5z2aT0n#M4vqlnQZOB72_u{-}Y6e2h5Co|x zW15w_v?8pQ1cnrF1^p#Sgm?girxmqi_hHXCFeRXf>`Ir$C6IHXv#`|=jYKXf8Q7(I zabO}b4u<9QnJnlQR3{wFY(p}KpICy0hcie6Iaqc1hKowbtK$=63DY%D%retbBuh+N zvdAIVI1-OGc;8anC(YvDs*~)BuvO07G`Xat&gM)@I9`}Dy_aHci=e#>qG;->&c&Nh zVcno&UsqP7@p?5$WyN9G;eCGwS}r#OS4g@19phM?^4K+UI7aTCu3qVuv{3ER|0P)H z7-`2VplU`fIkDNxA^=A1CeJz2$+6O&M}@}Z2LS^&FxY4KP!(&56gN$IeLp5N1v#GUX%{W|u3NN+VbRd+XH-FqPi}F1yc1l9iei zELGU5v4Kqe9<}#dW_cDnVspFH)!q9NrDep`yseM}!16)>Gw-mD1R2^qTt> zizRZ9EHNwzk~?29xNoB+CuStrHLbvID-jyIHnwiF0@_vD<A`Z1gPWx@RX6$O^T@ zaVHJ=77kghlgOn=ir(ugI~rR=6&aw293)9{j=yf)){e%_&O$D%d@9>im4!*`B3GhO z5&?%{0MlWsKcebJVniS6Y~N zrFagmuoXbe0I$HgT_{Gj8(?Jwk>!eXzzxFQkv6hU!%1C5l2AR2fbdAbV~G0KgwiiA zT?;>e&ieynqGBzd{rIoaKNs;{Lz&#FfB=G(nPBxMIt2=eMxGkRF3-@jv8#e_p zd0cVlHtT3H$B(#2$++B9=({9*OgM*E09O)`D>#=U_aY9?T9|C!9SnQ1CJHd2TzE@? z@cO^P;1kF(2W#ZY)Rc%N{n(DxJa^FJiyaj?vXflDpEW@orJAT(Xl-v{*mF;vm zG%)9uB6u2sIpQ^eqlDN$IKpl5dL4p301f+LV|G-vONCWTT*ii6vs$k{Nh+KV(rCJXk(Gs3^@1pmd~g!#2r~K3qJLo80P~TnV2vJT2{$`^h}X6Dn0X5RIsn z=R7!aPLw*Y`RPg00j@dE#WGO@+hGvVu=OP33C>*&TQS5K_X>d(^|l!~crk3z<2uMT`#yr-N)&;*qq&rYKOcibFdkId zRTQ$SOD^EeRW+xNMcj?~d)Cv%#uYY~#}u+`OAN3=gQop2O=CM8_g|^l zk~pi3kUx119f{Xv2}Y(@-+pzNDc}lOnAD8WA@cPy^bQfMq`A}>Mi`^2X`g#?0iM7* zs!F_OO5bqo5vmn$Legf9aq4FdU4-x?S2GSd0kr_7=1_}ZrXIXo7BuHmY_isx#K%^f zgO=5>Ul;Dg3RPQW5ZT4W`g?MR(-3SpYc-QR?AL7!&hWTfkdNAH1UXen@h&;JLP0z= z3o2bV5qfm=fzYxZZ1oB&vfp#>#t9lhb

=QsfwsZISI^>P%M^=GYT{4a|D4Bp{AY zMh_}V&$~Q#D#G&fB7Ly%gnI02LQKAU2_Tdvkld-4eJm7t{vGl{jt1i~yQ@;Tl>^lz z6`e~C`glV50YZsRGy9Ej9a9ybGbC~wBDiRH*WXp^3htwck6gJxE0@S(_!8?CoB5?K zc~q_YU4qMe2{$Y0taZ=~d4f$B8Tpga?3Qr+6(0Q zal@!5?j|s>OvePu^0XyERZQWSr8f{@rz^HAwm%=VQVdrHQ!@lVVXtxuTRhh}uVk+T7fdRh# z5HkbtMRH@d`JL6ofQ~Ez38)_iacjT2+T6|sZ34JZ*^HKw%d*3-Zj&Y0E?JCM`lneP zoX(R_q@GaQGd2c5FBN#Q5n`RhWLzu^fyNfNc~y%tOqGc8=45Wxgu-=bUYFz0EjMOnF=tPg z!`lf5&I}-40yC>VZM%e)s(FC!%CW4WHy|5gmTRpCqtn((bvR)*-iGXYb3{w;Wi=Kv z4$hPrk>P_E8ICmay75e`5s6q>v%0V@WH+cYz@jB4W3M1;*s!m*yJ z8?@XvF6$hTQ%NCjx3SryB@U#+nUd|uY9blQC?2JSVI|eFap)T^0j4i(rR;eKo9Tat zW$QW3&{VULbe0qWRdWGbVNT8WKV!!)UV70J&uKh%kwPN&DH2%55EP681uv)zh@=Hc zYjJ1Ens7`RF$EryCY*#D3rJKorvR*_qR|g}OR+4RYshoeOHO#M7Wb7^$+s^*fY_ZB zUF%j|{R}rTyp+s}*B(kEAX;h&!D+Z6rbuT@^y*~jmnUUy=Nn$P)9)x3ST{EUx~VoR zQ%u6SDIr>ziIs3w8Uh}fjT8%>_FLC;i(YZT;3w8DliJSgX@;!iKLx7X-Uk{%OV}*$ zL;pbXo_>{ucdoHPvqCab?3J$tO>vy7KQNk|B)Fnx)9Zsv#*;8llukI4ika)pjlxTQ zso*sP?K7T%A7Zoc_2b;&b8YGXBG~eL0v0G=ug%@xs}fH|!U$XDhMdJgl7+lY4e`C- z(!~wuDuZh}E@Q!i1ciTem4j`9w#fKK!E>!VP}L(j&mFI|y46*iCrxpj2I1ap&!e|m z3vYeB?h^UBLv5v>wiK!Y?{z-ZOEYQ`M`7+!rJ#!3YIZld3p({R&|cv*t9!;fy^?AUo?r4K8qcaK2NjR9smF{&P)D;3WF;$qNEsd~Cns(lR7UXng zexQJlg7yHrbiDj`N<|;6scK3CF{cnx)nPffgri+)gGf}$(*TKbC>4qekOCvQ zay)Mf*xm-ly?_&A$%$b(k0!|<>lO=QY$^Nvu8mt?FSGG(WKkd_Js&SK{OD%I`b3f} z(RzG~fJ?DDHgaGLu98V4>DrP>EP6KHZfwkjv4agsAqnQlK$J>?q9jJ*f<%O}$rWT3 zB@(4cBof*~9Y{!sbSR&4gl#HC`oO6nZ6u5q*d&oDCbTgf!%mNgZT*O`ZYA)aRe`MbB`tY@Gex$+||C>rbrG)YIQ(Sl z+McrX&Apv@bc=2OpdD++#OtmAw%(u#m&404FR`)!Z8LIhy&w#XT030+&4;15_1LyL z97uHc{5IXEywJ?ap@)0Mvj?fq0W!Lpc1a|XHL_p3)HLkd=jVW!rtukAFK{$JXC0*g zXEvD;i~`UN>djyb`_3o_TdS(;fHUd8shNNR%B)s)GrKqhV~+G%DcY9x_Jl0_x^uYfNAmgj(H0A$~>u#)V& z3Pobt*(nt_6{1Q2#Vn}Y(Fd@hBca_E4PYf>UW>eM7nYA01)s0r8h766u_I#nCC RJ8l2s?ntK!5*-&B*`O4RRYd>* diff --git a/data/AS.rda b/data/AS.rda index 290f944f8bbfad2436a28498921933ea808bbf52..97965a40b0a7bef5cad4d6c354a277d899e395f9 100644 GIT binary patch delta 203 zcmV;+05t#Z0r3G1LRx4!F+o`-Q&~il?U4;Ie^kk)6HORHMus3VF$S72NhJZHqd)qAnSRn}^5;F?`W+y!p!^neNFrP}4|2`N5DH>^90wQE5%^I(( zZs_sv;i3Y-Ldh+F43ZfP3Cx5!nS=*w&L}WLBtS4q{zEGvDkL`AfzVKb5fNorM!y}v zH5jgoQlND0J-G z9xdDRyvY(ltO0~GAuuo}Fp%aD9bs`nf+;M3!7ThovV^FRwWt?L5JDm{s)*O;zyXf5 zGNw<-yJ-1sSbwM~N$9qbHTzAcmhLM9>e+z6fmP%RD^kZ~tM<`F{x0N-aG@a+4Waf( C9aiK3 diff --git a/data/colitis.rda b/data/colitis.rda index 7af25b08d11eb751da4c1784c5a8451b0d8ce635..937189f482cb03c823d01200efa39c7e00edacc1 100644 GIT binary patch delta 209 zcmV;?051RO0ptM=LRx4!F+o`-Q(2ce*wK*=9e=4qQ&M?BrcX)gc!)Fr84pI8K$O)2 zVKmSH2AT#y1JoD|0iXa4G#UY*4G*T3FJmy221+!ANCMh{LBmEsVC0ZQBdaNZ2rJ6O z08p46>nSu(8I%;-2+$CNA(%`iXOjT(5@utHp#$h7(9@PWs1%k-GbM1YsLWf3}#jLrF)Yba23;4T| LDZ+$=xz5Ij_Hk$s+7HK!eALN(iI>JY6S-k7{Lc5!4Qo( z0R@*CQd;g&rkb)z<0%m7VpznO389J&AcJirs?!CTuB*0kIJ)zC z7gc?r4rz(td4S4Ffdm54dlJMV=}e1>*L~)Dz<(*-l9UFvtr|fG-1F4anB0i*XbNVe mwVMS%0tvtj5oJnr`bVquqC-Ye4ZVPX?k?ntaG@a@4F96m5@|L7 delta 236 zcmV0GK8rlO~u1#L49_Pf3WW`lcq( zMgR$u6954ihD;<5wqUuZ$;V2S;Sl5rkAh zw5s`JeJh?smKG|jhwA`Pio&%`JA7jAjK$g@@DgGh5@5v!2!ly2zHnLEw=I0n=T?7W zT~%-gG{o>dN(xHA90D{4jU1v63EirYr&rTTED`P>F;&MkN zv&~~pfSS`~Oi#8|QUFFjGbE8)07@Cn{))P~`#W=L?(flwa4Oj)GC;(WLlT1lRb)&M zZC#QAiG*_^RFBTFO5DLosD=#RW};Fu5Fo*_KeiX%eFFgPTHIZBk)UjUr4j{#_k5$7 uSo&r$Ui>T9d)p-xxz)Z26atliqmqxhpW$^U>&Pli*3?1=E!0g<9@|WS)mQs9upcDOL(OC;h*gIwiJW*9 Date: Thu, 21 Nov 2024 12:32:52 +0100 Subject: [PATCH 6/9] silence expected warnings in tests --- tests/testthat/test-gMAP.R | 46 +++++++++++++++++++------------------ tests/testthat/test-pos2S.R | 22 ++++++++++-------- tools/make-ds.R | 2 +- 3 files changed, 37 insertions(+), 33 deletions(-) diff --git a/tests/testthat/test-gMAP.R b/tests/testthat/test-gMAP.R index 61f0aba..be4b222 100644 --- a/tests/testthat/test-gMAP.R +++ b/tests/testthat/test-gMAP.R @@ -105,14 +105,14 @@ rate <- round(-log(0.05)/2, 1) test_that("gMAP matches RStanArm binomial family", { skip("RStanArm has issues loading since 2024-01-02 on CI/CD systems.") skip_on_cran() - best_run <- gMAP(cbind(r, n-r) ~ 1 | study, + suppressWarnings( best_run <- gMAP(cbind(r, n-r) ~ 1 | study, data=AS, family=binomial, tau.dist="Exp", tau.prior=c(rate), beta.prior=cbind(0, 2) - ) - out <- capture.output(rstanarm_run <- make_rstanarm_ref( + ) ) + suppressWarnings( out <- capture.output(rstanarm_run <- make_rstanarm_ref( stan_glmer(cbind(r, n-r) ~ 1 + (1|study), data=AS, family=binomial, @@ -126,6 +126,7 @@ test_that("gMAP matches RStanArm binomial family", { prior_intercept=normal(0,2,autoscale=FALSE), prior_covariance=decov(1, 1, 1, 1/rate) ))) + ) cmp_reference(best_gmap=best_run, OB_ref=rstanarm_run) }) @@ -147,9 +148,10 @@ test_that("gMAP processes single trial case", { ) test_that("gMAP processes not continuously labeled studies", { - out <- capture.output(map1 <- gMAP(cbind(r, n-r) ~ 1 | study, data=AS[-1,], + suppressWarnings( out <- capture.output(map1 <- gMAP(cbind(r, n-r) ~ 1 | study, data=AS[-1,], family=binomial, tau.dist="HalfNormal", tau.prior=0.5, iter=100, warmup=50, chains=1, thin=1)) + ) expect_true(nrow(fitted(map1)) == nrow(AS) - 1) }) @@ -162,7 +164,7 @@ test_that("gMAP reports divergences", { tau.dist="Uniform", tau.prior=cbind(0, 1000), beta.prior=cbind(0,1E5), iter=1000, warmup=0, chains=1, thin=1, init=10))) - sp <- rstan::get_sampler_params(mcmc_div$fit)[[1]] + sp <- rstan::get_sampler_params(mcmc_div$fit, inc_warmup=FALSE)[[1]] expect_true(sum(sp[,"divergent__"]) > 0) }) @@ -172,31 +174,31 @@ do.call(options, std_sampling) test_that("gMAP handles extreme response rates", { n <- 5 data1 <- data.frame(n=c(n,n,n,n),r=c(5,5,5,5), study=1) - map1 <- gMAP(cbind(r, n-r) ~ 1 | study, family=binomial, - data=data1, tau.dist="HalfNormal", - tau.prior=2.0, beta.prior=2, - warmup=100, iter=200, chains=1, thin=1) + suppressWarnings( map1 <- gMAP(cbind(r, n-r) ~ 1 | study, family=binomial, + data=data1, tau.dist="HalfNormal", + tau.prior=2.0, beta.prior=2, + warmup=100, iter=200, chains=1, thin=1) ) expect_true(nrow(fitted(map1)) == 4) data2 <- data.frame(n=c(n,n,n,n),r=c(0,0,0,0), study=1) - map2 <- gMAP(cbind(r, n-r) ~ 1 | study, family=binomial, - data=data2, tau.dist="HalfNormal", - tau.prior=2.0, beta.prior=2, - warmup=100, iter=200, chains=1, thin=1) + suppressWarnings( map2 <- gMAP(cbind(r, n-r) ~ 1 | study, family=binomial, + data=data2, tau.dist="HalfNormal", + tau.prior=2.0, beta.prior=2, + warmup=100, iter=200, chains=1, thin=1) ) expect_true(nrow(fitted(map2)) == 4) data3 <- data.frame(n=c(n,n,n,n),r=c(5,5,5,5), study=c(1,1,2,2)) - map3 <- gMAP(cbind(r, n-r) ~ 1 | study, family=binomial, - data=data3, tau.dist="HalfNormal", - tau.prior=2.0, beta.prior=2, - warmup=100, iter=200, chains=1, thin=1) + suppressWarnings( map3 <- gMAP(cbind(r, n-r) ~ 1 | study, family=binomial, + data=data3, tau.dist="HalfNormal", + tau.prior=2.0, beta.prior=2, + warmup=100, iter=200, chains=1, thin=1) ) expect_true(nrow(fitted(map3)) == 4) }) test_that("gMAP handles fixed tau case", { - map1 <- gMAP(cbind(r, n-r) ~ 1 | study, family=binomial, - data=AS, tau.dist="Fixed", - tau.prior=0.5, beta.prior=2, - warmup=100, iter=200, chains=1, thin=1) - expect_true(map1$Rhat.max >= 1) + suppressWarnings( map1 <- gMAP(cbind(r, n-r) ~ 1 | study, family=binomial, + data=AS, tau.dist="Fixed", + tau.prior=0.5, beta.prior=2, + warmup=100, iter=200, chains=1, thin=1) ) + expect_true(map1$Rhat.max >= 1) }) test_that("gMAP labels data rows correctly when using covariates", { diff --git a/tests/testthat/test-pos2S.R b/tests/testthat/test-pos2S.R index 049086d..b6eef52 100644 --- a/tests/testthat/test-pos2S.R +++ b/tests/testthat/test-pos2S.R @@ -1,4 +1,3 @@ -context("pos2S: 2-Sample Probability of Success") ## test the analytical OC function via brute force simulation set.seed(12354) @@ -51,11 +50,12 @@ test_pos2S <- function(prior1, prior2, ia_dist1, ia_dist2, n1, n2, dec, decU) { ia1 <- postmix(prior1, m=0.2, se=s/sqrt(15)) ia2 <- postmix(prior2, m=0, se=s/sqrt(15)) -test_that("Normal PoS 2 sample function matches MC integration of CPO", - test_pos2S(prior1, prior2, - ia1, ia2, - N1, N2, - dec, decU)) +test_that("Normal PoS 2 sample function matches MC integration of CPO",{ + test_pos2S(prior1, prior2, + ia1, ia2, + N1, N2, + dec, decU) +}) ## also run a MC comparison pos2S_normal_MC <- function(prior1, prior2, N1, N2, dtheta1, dtheta2, pcrit=0.975, qcrit=0) { @@ -98,11 +98,12 @@ test_that("Normal PoS 2 sample function matches MC integration", beta_prior <- mixbeta(c(1, 1, 1)) beta_ia1 <- postmix(beta_prior, r=20, n=50) beta_ia2 <- postmix(beta_prior, r=30, n=50) -test_that("Binomial PoS 2 sample function matches MC integration of CPO", +test_that("Binomial PoS 2 sample function matches MC integration of CPO", { test_pos2S(beta_prior, beta_prior, beta_ia1, beta_ia2, N1, N2, - dec, decU)) + dec, decU) +}) gamma_prior <- mixgamma(c(1, 1, 1), param="mn") @@ -111,11 +112,12 @@ dec_count <- decision2S(1-alpha, 0, lower.tail=TRUE) dec_countU <- decision2S(1-alpha, 0, lower.tail=FALSE) gamma_ia1 <- postmix(gamma_prior, m=0.7, n=60) gamma_ia2 <- postmix(gamma_prior, m=1.2, n=60) -test_that("Poisson PoS 2 sample function matches MC integration of CPO", +test_that("Poisson PoS 2 sample function matches MC integration of CPO", { test_pos2S(gamma_prior, gamma_prior, gamma_ia1, gamma_ia2, N1, N2, - dec_count, dec_countU)) + dec_count, dec_countU) +}) test_that("Binomial PoS 2 with IA returns results", { ## reported by user diff --git a/tools/make-ds.R b/tools/make-ds.R index 66911bd..a06c447 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 <- "ff4fcab" + pkg_sha <- "bd767e7" if (gsub("\\$", "", pkg_sha) == "Format:%h") { pkg_sha <- system("git rev-parse --short HEAD", intern=TRUE) From 30521b80601d3842c199ca54e861e002881096f9 Mon Sep 17 00:00:00 2001 From: Sebastian Weber Date: Thu, 21 Nov 2024 13:06:24 +0100 Subject: [PATCH 7/9] do not run lengthy predicitve ess elir demonstration in non-interactive mode --- inst/examples/ess.R | 15 +++++++++------ tools/make-ds.R | 2 +- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/inst/examples/ess.R b/inst/examples/ess.R index da70031..264acba 100644 --- a/inst/examples/ess.R +++ b/inst/examples/ess.R @@ -24,12 +24,15 @@ round(sum(ab_matched)) ess(bmix, method="morita") # Predictive consistency of elir -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 +if(interactive()) { + # takes a moment to run + 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) diff --git a/tools/make-ds.R b/tools/make-ds.R index a06c447..0013846 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 <- "bd767e7" + pkg_sha <- "3377b7c" if (gsub("\\$", "", pkg_sha) == "Format:%h") { pkg_sha <- system("git rev-parse --short HEAD", intern=TRUE) From e64f7cf5c900ab53a8389190a65c7c9cdcc55221 Mon Sep 17 00:00:00 2001 From: Sebastian Weber Date: Thu, 21 Nov 2024 13:23:32 +0100 Subject: [PATCH 8/9] bump the Rd files as well --- R/sysdata.rda | Bin 6069 -> 6079 bytes inst/examples/ess.R | 17 ++++++++--------- man/ess.Rd | 6 ++++-- tools/make-ds.R | 2 +- 4 files changed, 13 insertions(+), 12 deletions(-) diff --git a/R/sysdata.rda b/R/sysdata.rda index 2064e02ddf83cd5c2c67f24f3d0315fc846a5522..f58c2da24d9098386cb203e905c75c71c9db3ae9 100644 GIT binary patch delta 5637 zcmV+g7W(P6FTXE;LRx4!F+o`-Q(02Og|q-qE`R_3|NsC0|NsB=BG^j5@8AFb|NsC0 z|NsB~Jb)lT0zeDk6n_9F^-Ij5v09l`vP{jDwPw_yjYw-vB~)y)Wu~)R2?);s004oW zqf^k2HBVE@WMXL027nD3G-v<-00uzQKr{^i00001p`ZYN13=I+s(N7?X^=fY000000000000xGEp!EO-RZsn($OOTYMnuHZCI~bpffU$`rjJn2 z0000000000000004FIaDfBFWDCXE^ZObL@fOcZ%FHBU+EXlMWc00000XaE2J00000 zXb2h_Gy@2K4F-XM00w~5AT$O60h3Jt01YrrG{FIrMnR(p05TXuOn|8Yu7@l|NMV8foa#85m6-rkY{~jXfqz(i&|hKo3x9lhk^H)Hb2LQ`Fw3Pt&zt6Fd3JJk9c?P5H|Ml_5?bGDNhFdt_sL|7+tOL;S!KOZ zI_wu05EFb#AWp265FKGD0B^}sOeE_`_x-UpvV6@wi(9HPadgD{*274Y;$i7}vvF@? zO&_^`Sqdj+PjpMI~|6U8TtPUE}oxHfh> zejo7upS!-U|ApsxzeBy{^D}FigZW?3@P;veVj(Kc1kCV@0k(Q-Y1W?oL^2gsTTwzxwN(bjiU%=Wm~%aw zV1%InPjFHICNHoK-1OlKI*6)OY#a+M;v*Jc15nXot6BR#1I~Wf+f>>-Ja}l2f~P=# z^={FNQQ*Ccm;=$H^xvw2xGBm&@hDYkFIbuv}7={DwkmRAF-WmA)Y^_$a@B22wW#G}k_N<(Y<&FEeNi8lkp@T+*13 z)}%Z#g9hu({9Yjq9=CPrUg`PVke*(A<;6Tmp_EA5q#+Bvm|NV@vn*QPgO*Cn?r7Hq*ndbF^laI&2T(K{m^*$J^{9JnK(#fdW7~)nf05Tk)J5R!|oZCW}5~BudmZAfrZ= z$T~|J7XwU!zPPgJ7t1Op*A6$Is`g)h2j&Y)THFFbPA2j!q+T%Q6n3?)bue};>}N8p zCdROT5N$`63xreTkMiMv4XVtG+=!^;xO-_U{s$Qu{@myS?D5q-j%;qUEYAg^bD)MI zkx2vUq}Zc&&=RWJfI#Mi$~R^&U62Z@1BIpVnVQdR&T(zWq^&m8@Lb%5Pl4Ce*9O%V z$)$lXdJuL{S;-|F=O6Eh+19zx?Z>)$+F6xorZJ41Bxj#pZB~naM}!I}RqsM|pvL`b zbTV|b+{@abszADj6(Q<)fXf`_Aj^0~aCSyR2|I&5b>PJ<4`sD03l_D&q+0H{W5ye?^r0bGF04vk04Zc|BV37ZXP4Pa=O z)2>*pkoMTHpx#w~W5$Yuk~=*$=~SdG)@bXF^w_Z4H&Xv&q{g{B0DVkMuk%5hnT|tJ z2sT48Sz|daYzil`H;Yq7zynx`K zI~l*g=WW?ubq1nnot>~)Jfcikg^0xqbBJu_xIEc?Ue=d?fZI8&)xBw&xY-TmanK}S zYlGB~p69b9)3p-*otUUahVLe@{6cF zZI#uCsghvQT4Par^}SFt6SKwMtkjEM#}`wB8e7MIcS8eeh=`4cMSEla7pn7*GgG*sR>g%)O8Y|)hOkq!Up+iO+%9CB*dx|ZRtx-3j9w*H7*lK z)VZEfoO*K_rG|_hh_m&;B>MFdv1#8O4dzFG1*%#y!iooA6SKrX(tWZ5Ec8AmYeEej zbW>5up2dj~YmrwO>5>Qa3~4-@VB02;hC!h&-S^Y|>02o0yvq8Xl-d_B4{(^Fv|e7x&(4@k>xZO**Ps%5tuI`2Z_Ep8W%w=q;UE~+t-WcoRs(~D*f z;pKT$VzfCdX~VM7&tG#Elbz|}uC<RTR-6KoTfuB1B_-DL^*JTZ7G?jq%E2b;YoEV zpcBawv>JJ}Cmp?{>#Ibd3q#_TML#nDt_ez)7Xk@MR*J?|(Dgv2cCnywNy1fsrd6g% z9k3}>8FbbVOAjrM&fcTPl3`8Sz;Q)*B|zff+uGc+BTRhSo?SyFiB-*UnU@|H-vp+k zO)DJLJDgy~B4t+~YSo1do+Gjtv|fSiGu=Z2k_FtvQ8xQT3gd#KH9p-k(tjGD^iGuB zbHv|tG;R0bg*q8oPzq|>Tn!FMp;l+*s7!}rTI+)-Ut)x^JUzF<<4T?uk0$fe zu;#+F=@G?Br6y2uT=qPB+kw)6JM*>f^vh#+YQ%H%j`uaP z5;K+&gj3TdFvNnW&ETmyhTBfS;BJ*o>LVVg$WsMyM$W`!5`rs#GXmMG<$*XLq`_&# zhjo+0b=P`!AEuH< zK8Os=m!pJg>W(2`?!$MZG_Ypa)bTBo?yh9I)k6HSp5JuWsfGu5@<0AuF{(iOeS5qeIPS_*G5G)JT zE9ueUK^5^-bu}G|D_T0<+=VZB^5?vrLwPkMg{?>*rkq z0{%_jyz=EBE0miLHo&FIpUkmXV|~1{FqN05OI^Exd`RGbF~ej~yIuD8u-**I?`LIb z?b5CYXNoR$=kY^jg_2q&2<4&XBbOU(rKT~l!{Vpucgxj|;A6lFGS{7+)r4q_qUXDSEt7u`j)!9`MmyL<1ET`U{{=-?Wuy{1jz|YUlKk?X zT3*uC?G=`JaO7Z`)KCnt^KMqosTW>inA@vX&k0rfhp-Jhy%*Mv8pUR=oA4Do>E~{J z3g#hyDw1MqND5Ilnedy7%Ezoh@4a({i7pH1f>)W?wGb4Rmn~QF(84eVls3<~c^-C7 zvfZU&_E$SdDv}BTB#@4}hP$J<8fIyF3fz0DUYir!G>9pmOjMp**w!r7jR=96_}R=B zXLdkC0P=-*o~EB95BF4Z!@_v`>bL+&TR^%DhY%;VwOH%nY1+5ClV}f)4NuC)Mp*X< zMXab~@7R+A5Zwa)V3Xnyj(_{1$l{S_?N6PcrVMo7ro1fh4)*v4hg!TV5?(!wN>}K; zTYp2!>;VMpP!O^z z2dx{kK5-4_|xJi6!5ye4F=?Ot+C>GF0 zyS5^qA+gL%e5%GMtBp(`Dp6F%XuKhZSp>Vn$PBA-SgdZjb&3IOtAAdcZ8cd?EQ0)N z**oer-~psCbALCBGShj53)x1zX_ZKTFzK;~6dlcw67B?mIUAh>5o@Oq3|7oOWo4uSTk+ zfp}CD@q*dTNgr{l&P?(vqPAaQa0q!d@#ty$7}wq4eSbY@-;SA*jXiINV=>WGlvHWu z4oE`wOG&kok9jYaZF+TTqrs**AF}2inVJ&29Qkhu&UaKa?JcKaC<|kv#Zb9en5-w% z$vhkA(3u#pe0Ou%S8Os{i`t>@qBQ0c1#)E)h%!yv?@8aMC>=er#RJwN{A#=M)K{l8-n;4D zcp~r}ZF+L*xx7|@O|78tKH{&uLqP#{VGG*3PuI?H7yR;s+oB&BJ0|1 zUAO`?IyYmEdQ#}t;5mld+Uth|NvKWEy8Mltv2iH}XKQ`iC|rQvRF;>un_^5bq9TRh zD1W(Vd^VzCtS#)jnO3C|*3m9;U(MuFb0-e*A;H;7|36A>5P)s4-0Ol%Ae_#Y6!+0f zAU>_p%oHL!SA)mla4v(9Koj$P_tUE*Yaj5c1gD(oSpB$7zo$2*^}Y4cC}h0S$1Ip;S|p7ZIvMy~He zQ`YW%W~Qd6q15W@S&HWSs=AX7z8^#bUG{RC_W&)o(*ZJiyt6TRg@^+)D@Q*bYn%a% zmyx5T%>RDW!XB5f-_0x+HGlgnHFCg9zMk&MB$7tsI^8vD=zy2zqB@52d4mtK z)7tFYxC1qxjTnF}0M=Qm0iSi*0Bh`QF#u@OTBaM|7GtmliUCO^l9q0|dzXML0XqVF zdvZWBfHQzCWs@?{6R7pj21fVX=rS0!=mgBkto_y$%tVqY5=f-Ia_|D^307`s1~3dw zW?KO+WnWS&Cql_crbTI1;FJMqvZHfSBRD5y f={sN(W`IkGj=rN&i`BVv|BJaIoG3^sVZzz~hp?Vc delta 5627 zcmV1Y*waKY?Cu(ZCSM_BT^dENmUyx*=elS!a_5^000!J ziM28hO*EPS0gwO#LqVVeKmY&@G#L#I83RU|10eM@X^7B&00003DyRNrX`nP{6AXhu zAOvYL3VMw_CYb|30001J0000000E;zAP1-bWmNy#41i1-G-OOoGGK#3Xc0|_%4qcs z0000000000000000MH7mr~ja6!f4T;1i+az1i?p>Q&jYxriOq30000027mwn00000 z13-bHqd+i!fY4|d5CCWmG6O(h5E(Sk007ej(@YQt8fl{-XwybbHi@RtAPpL6paVb+sixF=gVYR$o>OFcrquM& z&_&86QB_n_qN!9;x!@;4NgzT50SF2wCqf_s5CwC85H&5+Av*O~>;s|&OI54YBU5>E zdlnmJjhG>S-{b!8FKy)6!_4QHv&tKJBp6B|J>dkEu5n?2}9ERs|D? zPPmzWME3;uB=wU3G&I_ljak?Cm5aiJOX$m@QK-Wy{Bzv&o{pQpX2#= z72Nk;cZ2hLeoretxmbU;R^Po9uJlI~MaCB1Iu5mc(!I2K#QMl8Mtp`yiBv-bQKobjV=Q)uxq;i5hYodeaoMlDBw zg7z+84@QsIe}e^Wp}5_3)9%AIT%yb_*13ey_$5rhVna{1moOv`!Lj}26WUc~K5cIJ zw{%3VGp+g;uD*s&%m)ssUbn`EQImfBin0e>P~&~UPErSnLaRx7#L+xRz%jB06o2;S zP6M-5%q1+m%{W0rBsI9X(n4Zr3aVs(S~2rgU<>?_6HnBfNNy-&fnKG$BwmxG>%w)T zm5^Jl=&s}PSqZ68h0EPm_=LOQqW%0CNZb_DT>0OYW*S7i&7d%9hT0W#N@7D=knqY3 z8?QI=c!V^1-Pfgir{{7)d3p1f6!9X4Q6p}Ugf8}Bg_#XvZFCPqFa~&^%yYtj`c{FI z#QAg$&ME*+5%3nB9w9mfWmavjqw2(!!C#%5)!0moByzp;_pSj$F=8ZRD%7 zD;C(U8D&1)xe60bEKd8$yMoo*xji4p3*FA-$Juj6Pm^xbA4w7svm?j;QQWXq>}?!x zkop}>BaC9Negls7nsF5=lWM$wQ`dYtwlLli+chbUmyHTC=ici@wxZR0Q)CGPcKMfT zU5aDHG%}&$a!e&HndCzx0wpjLaT+jbTnyriQ^X8Eek?TwsVl_ip%3SV24hn&;PvH{&y^5K(AraVPYBfBwnMmoYC~L8R(p_6 z0l$#itwrJ3uA66%x7bE`)}G}81b}y{#orRQ<2W{~pe`Uy7JSG^m8fb#MvW_wbe1$O z2AKtYab?gimQ+oy9B)5W?7sdF%odilxCDZnP2^cfykX2J?Q2}>VC-4g&Sh9ljbQ#j z+K((32&c&(=fWFRnHRZ#5mCu-_R?4U4l*c0AQe;x3rpcMHJ;g=;@gi&T5YM|xw#6T1FxyB4XQ1ZO9Ei@ zAnc&Cl1e$wKj9Ozt#hHk1jT+#zd;p#+HeEHlUB`n0l`3aGk<~3 z+p@gs4Mfp9J7BPRM3}J)5sDV)5ZTRed9wPwtuFz#b6Kl@deb#=vKz|dph&>i2dN=F z&t^%dY9~BASXn5}usY?Et;3W1M%_jZUVO^cwZeRlk;^m=&NaO?G4IbldXi_@6ox5} zUvoR@R##BYIvEtrjg`}QZNH&rq8tcl!{OO-u46DT+1s~t!afi-#|ukmA@Hg@A_|*B zo{Ls{aops8&<)s?E0oS%(>oy_Wjiwx<%BxmJHQ%N@eHjvijX-9E$tBBHV&F$y08i1 z!rhZ|Ek^Bu0alDRgRniR>aituE{ijHmFJA(g!`;^W0cFW6wsvOG8uQ}7f^cJE2|Mx zB*CP##-jJ@dZ1<}XN$X8sTRGCE~f@Gw~p?H2GtRN5gQJQ_Q(EjRp%aNv-32AzA99j zJ0n=pOh~(Hm0-B_)~I+W zrgx}+>Lo#{QOiq*E=$)fVUYYcU7AoO7j?{brLwAe$CG_JD<^NrlXTR zixMK&BCa#jBoFEs(s?$)woM@ngF;=q@2C7ywo%V|yV}!;B~{y#jsz#kt}zSX3N+4C@B2G2spwZBfP7R(*P%JQhi zXmVK7hh?LlzUD0_JJZEoYdp9i-t6zKtSX2+RT+Rr9#ISIjZ7TEyA0-r5Ts4;wadcz z-Gz!97Tm&ac)AAj>lJ1@?suP`q0MH0uIR{k>m>WCV|(L;bTc|=Y>MR02Q<6Xwmljq z*e5MrNOmo{$`kMSV(?YmMbLcO+b-bj8i?aw!s&$>VZn_V?g4QgsV)eOp-f)U{a_u z>8v1@9$Ouqy+@BE!ke{#;)?J}fyKeMwYg+QnEAClx`s;f?5_vjtS`DCyZKruFpGvSe=GQdRH(*VFh5L(1woL_6 zUm+20<$F@ltk;XW5kf8Wh{L=*BR*^ zj?|O0492q2aCv+Q9Gdfgltx@aB&#peM$^!>*9K6&#R+71dvArtl{_sTP3Na!&4p>w zBZ`$uOrYYq?0EOL1Em0W=WE^Bmd5VYi44jUD(?(dB`aY8M-q{P?rUTuXDlNK zr>0C{i3L%c!BTS#ww-~%-71^ZMmnDio zuJr6fn-R$;af$UFDsERx5)dR)1*@J&qh6- zM+nu`974g}hVMpcV9l|~8imyCbzO?JPjRgbL1m~Wn=ufwZX3MJ0&bkH#{F}Xb z<;p-;DK;K$flHM?nPRZU`*~(zD=$u#yLSTkk-%ey$f0(ByY25`ycw6?&dSl-rCbot z6kO}i;)cr$G+UN3a3#W&01j3kG)UqCDvBswVw0x8QwJgWD^C4S*OFGS(R1O_=H6?@ zokzat2Jzm0rs8$_Tv?k5w{BoQwKw#%EYj6VC6u_m9D@~B(8F)5vZ^C5CSlm!a@wr~ zQ;ES8ll4V^i&M=p)7ur>bqGzzdk?*V;SP%yy9|7vYUmkbTo15`aB{Q*R=IVt+-oi6wp$4eckc;QU4-IOR{ zhbLNFRdCiip^a&;3{|1&fNMH%tBQp z#MF>~6ryc2;WrnRk644>dglugTo=&_TCe4ygkTLQZJ%=TJnWriyGq0C zu6B@BBoqQkAsuxMcSmqE%+mA~xc5}OHYc`e5K}&wsXVu_tXZoX5d$;vvzRQ-?0|*= zSaY(cFr_Rt*20CxkUKV%9dwc^!tzH$0FCNAvEA(FR`!}*& zA6yWCn(U3Mp6aH3krQBNbmg2(yhN&e_?k;76G(0D?M?Vt?=B~r+mmd*GEJzejlLs{ zyHZ?<2w4>a*Gkg5 zv#cO0v)U738Ud}A#Z*MC#`;o}x`T))p;58IK8BT9YND?L)FyVjjQ!)cO+orIX;G(i zG>jZDDJ)%NFUwZNfLWsm>jOC^P*Qv*kwi1PWo}z#97~sbBr9f8WlTw4rh#g&6|8D} zk4@Zmwe5O>f3cFXrz`j?Xbz01QK`!vJ11j8BD;us|Vt`vJ-`A%bO;%J3Aio;+Pn@zf ztqgDg5ClQsj4LQ-s=&zUM@e1VTu&#KbPmP$;{$62fAAK?H_#v@8oDTjbnx!(Ua>gw zSF{WpyyCqK#;Y-4K-)>%lw_r-9VXJlDMTLos}2(=WI^H+v1@H{c^Lu96%s>hhDH^S zyuq2dj>`J-S3JHXkYNk8)5KSUpHvE_pm;ov$#i*HS*cz8sz)krt`bC3i0+h%2BsFc z0BH;yf6d~|wBBJt_EE1IWl|svI&5ME2XkZvR`%IT5q0wMD1u3On92e5z^<+-pAh9hxF5(W}1UY%N~@9B<5?74@gW`wTCK3l>woz)FHOKI3j0@&!WR4!H~D+%>7PX78d zCPpkD9o+U++YFZC_NaTvjX8wDT$x0o43l>I(s${~2TyFVK=p{f8m{DNTkF3K)4XCO zef`6M`L8%pzKW5++@|j;@v?UMRH6uMYK( zGHTS{6sepxMUsNu9Kf~yR5Z~!sWL;&SKlddc97~4^5{Vp*EKpUyas`Tofoy92)qYd zUYxqFZxx>5>lFffRxxR-ae+Uk)1LJ~S^wpIQ!LH#+*yQ2|I`usJ3(Idt zx^Xvx0Qlj&=|+X73mttnS0LmVQ;b{H;PQz-rAT`rNWi|vrJ>sOQ!gdtU3*Qdw*W?m zM(lCVN?jV<2Qb@vU2x!uH3_-bUy-wxE+ruB?Qgqf3y>SilG66GY)OU`L{Pj1e-|we zhSW?og}s+^D%7G{+9l2_`MiowWZ~W$N$y1fb^s3s)&e^$hi4mni zB0^bYiy*Q}B^4x)OKKz`ETJTVb7-Glks1{Xs`{jbw308|Vi6je9BPqho;L2hYe{c9 ztHtr|CT;&mJ1fV&O6jp(V-E89E!Dq&YVkb$zqOh6@#R}pg3lufGG=Dqf5Uj2YuN8z zyiK*$b~}lG9bIM{h;!;8rzDa{+{dA{uxj%TcP-3yIChR^uG5!I>OE$zj<-#z((7w$ zdRqH+d3=^n-l_O;;ju_#LYyKDH2Gezjg2hQQ#8X@C@J#jq@#p zmt)FOD;D0#NTv$uD{x8xv{_N9Ne{_U9c@A*%A|?jRER*HolYG~Ds){JVghmK1i2Bz VsNFI8mtE?Mxgwk>NPhv#*8uz}iv0ip diff --git a/inst/examples/ess.R b/inst/examples/ess.R index 264acba..dfe433e 100644 --- a/inst/examples/ess.R +++ b/inst/examples/ess.R @@ -24,15 +24,14 @@ round(sum(ab_matched)) ess(bmix, method="morita") # Predictive consistency of elir -if(interactive()) { - # takes a moment to run - 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 -} +n_forward <- 1E1 +bmixPred <- preddist(bmix, n=n_forward) +pred_samp <- rmix(bmixPred, 1E2) +# use more samples here for greater accuracy +# 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) diff --git a/man/ess.Rd b/man/ess.Rd index 717c6ee..464388a 100644 --- a/man/ess.Rd +++ b/man/ess.Rd @@ -120,9 +120,11 @@ round(sum(ab_matched)) ess(bmix, method="morita") # Predictive consistency of elir -n_forward <- 1E2 +n_forward <- 1E1 bmixPred <- preddist(bmix, n=n_forward) -pred_samp <- rmix(bmixPred, 1E3) +pred_samp <- rmix(bmixPred, 1E2) +# use more samples here for greater accuracy +# 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 diff --git a/tools/make-ds.R b/tools/make-ds.R index 0013846..14cd272 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 <- "3377b7c" + pkg_sha <- "91e8f18" if (gsub("\\$", "", pkg_sha) == "Format:%h") { pkg_sha <- system("git rev-parse --short HEAD", intern=TRUE) From 138baff4b4d5721c26e7d894b8d4bc0b1fa2c4f3 Mon Sep 17 00:00:00 2001 From: Sebastian Weber Date: Thu, 21 Nov 2024 14:17:21 +0100 Subject: [PATCH 9/9] fix Makefile --- Makefile | 2 +- tools/make-ds.R | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Makefile b/Makefile index ad4eb77..ce32c5c 100644 --- a/Makefile +++ b/Makefile @@ -13,7 +13,7 @@ PROJROOT_ABS=$(abspath .) RPKG=$(patsubst ‘%’, %, $(word 2, $(shell grep ^Package: DESCRIPTION))) INCS = -R_PKG_SRCS = $(wildcard R/*.R) +R_PKG_SRCS = $(wildcard R/*.R inst/examples/*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) diff --git a/tools/make-ds.R b/tools/make-ds.R index 14cd272..9fd1a33 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 <- "91e8f18" + pkg_sha <- "896a402" if (gsub("\\$", "", pkg_sha) == "Format:%h") { pkg_sha <- system("git rev-parse --short HEAD", intern=TRUE)