Merge pull request #20 from weberse2/issue-prep-v174

1.7.4 release submit 1

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)
     <doi:10.1177/1740774509356002> and Schmidli et al. (2014)
     <doi:10.1111/biom.12242> explain details on the methodology.
-Version: 1.7-3
-Date: 2024-01-02
+Version: 1.7-4
+Date: 2024-11-21
 Authors@R: c(person("Novartis", "Pharma AG", role = "cph")
             ,person("Sebastian", "Weber", email="[email protected]", role=c("aut", "cre"))
             ,person("Beat", "Neuenschwander", email="[email protected]", role="ctb")
@@ -25,8 +25,8 @@ Imports:
     methods,
     Rcpp (>= 0.12.0),
     RcppParallel (>= 5.0.1),
-    rstan (>= 2.26.0),
-    rstantools (>= 2.3.1),
+    rstan (>= 2.32.0),
+    rstantools (>= 2.4.0),
     assertthat,
     mvtnorm,
     Formula,
@@ -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
@@ -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

Makefile

@@ -13,8 +13,11 @@ 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)
+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('$(<F)', output_format=rmarkdown::html_document(self_contained=TRUE))"
 	cd $(@D); $(RCMD) -e "rmarkdown::render('$(<F)', output_format=rmarkdown::html_document(self_contained=TRUE))"
 
-R/stanmodels.R:
+tests/%.Rtest : tests/%.R
+	NOT_CRAN=true $(RCMD) -e "devtools::load_all()" -e "test_file('$<')" > $@ 2>&1
+	@printf "Test summary for $(<F): "
+	@grep '^\[' $@ | tail -n 1
+
+tests/%.Rtestfast : tests/%.R
+	NOT_CRAN=false $(RCMD) -e "devtools::load_all()" -e "test_file('$<')" > $@ 2>&1
+	@printf "Test summary for $(<F): "
+	@grep '^\[' $@ | tail -n 1
+
+
+R/stanmodels.R: $(STAN_SRCS)
 	## 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
+	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(debug=FALSE)'
 	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,33 @@ 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)
+
+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
+	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 +200,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)

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)

NEWS.md

@@ -1,3 +1,26 @@
+# RBesT 1.7-4 - November 21st, 2024
+
+## Enhancements
+
+* 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 `testthat` 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.
+* 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
 
 ## Enhancements

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:

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")

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,38 +37,38 @@ 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))
         return(.integrate(function(x) exp(log_integrand(x) + dmix(mix_comp, x, log=TRUE)), lower_comp, upper_comp))
     }, 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)
             lnp <- log_inv_logit(-l)
             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,13 +80,15 @@ 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,
               lower=args$lower,
               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
 }

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

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(...)