Merge pull request #25 from weberse2/issue-release-1-8-2

release 1.8-2

Author: weberse2

.Rbuildignore

@@ -23,6 +23,7 @@ src/stan_files/.*\.hpp$
 ^.brms-cache
 ^.gitignore$
 ^.gitattributes$
+^vignettes/.gitignore$
 ^_pkgdown\.yml$
 ^docs$
 ^pkgdown$

.gitignore

@@ -7,3 +7,5 @@ src/*.so
 src/*.o
 src/stan_files/*.hpp
 src/stan_files/*.o
+*.*~
+inst/sbc/sbc*html

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.8-1
-Date: 2025-01-20
+Version: 1.8-2
+Date: 2025-04-25
 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")
@@ -39,7 +39,9 @@ Imports:
     utils,
     matrixStats,
     abind,
-    rlang
+    rlang,
+    jsonlite,
+    lifecycle
 LinkingTo:
     BH (>= 1.72.0),
     Rcpp (>= 0.12.0),
@@ -68,9 +70,11 @@ Suggests:
     parallel,
     brms,
     glue,
-    ragg
+    ragg,
+    withr
 VignetteBuilder: knitr
 SystemRequirements: GNU make, pandoc (>= 1.12.3), pngquant, C++17
 Encoding: UTF-8
 RoxygenNote: 7.3.2
 Config/testthat/edition: 3
+Roxygen: list(markdown = TRUE)

Makefile

@@ -5,7 +5,7 @@ TARGET = r-source
 OUTDIR = ./build
 
 # includes all src dirs excluding R/
-SRCDIR = ./demo ./inst/stan ./inst/stan/include ./man-roxygen
+SRCDIR = ./demo ./inst/stan ./inst/stan/include ./man-roxygen ./vignettes
 ##DIR_OBJ = ./obj
 
 OUTDIR_ABS=$(abspath $(OUTDIR))
@@ -86,7 +86,10 @@ man/package-doc: $(R_PKG_SRCS) $(BIN_OBJS)
 	"${R_HOME}/bin/Rscript" -e 'roxygen2::roxygenize()'
 	touch man/package-doc
 
-inst/sbc/sbc_report.html : inst/sbc/calibration.rds
+inst/sbc/sbc_report.html : inst/sbc/sbc_report.R inst/sbc/calibration.rds
+	cd $(@D); echo running $(RCMD) -e "rmarkdown::render('$(<F)', output_format=rmarkdown::html_vignette(self_contained=TRUE))"
+	cd $(@D); $(RCMD) -e "rmarkdown::render('$(<F)', output_format=rmarkdown::html_vignette(self_contained=TRUE))"
+
 
 inst/sbc/calibration.rds :
 	echo "Please run inst/sbc/make_reference_rankhist.R"
@@ -166,6 +169,12 @@ build/installed/$(RPKG)/DESCRIPTION : build/r-source-fast
 	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
 
+docs/index.html : $(SRCS)
+	NOT_CRAN=true $(RCMD) -e 'pkgdown::build_site()'
+
+PHONY += pkgdown
+pkgdown: docs/index.html
+
 PHONY += dev-install
 dev-install: build/installed/$(RPKG)/DESCRIPTION
 
@@ -185,19 +194,19 @@ PHONY += check-winbuilder-devel
 check-winbuilder-devel : r-source-release
 	cd build; $(RCMD) -e 'target <- tempdir()' \
 			  -e 'untar("$(RPKG)_$(PKG_VERSION).tar.gz", exdir=target)' \
-			  -e 'devtools::check_win_devel(pkg=file.path(target, "RBesT"))'
+			  -e 'devtools::check_win_devel(pkg=file.path(target, "$(RPKG)"))'
 
 PHONY += check-winbuilder-release
 check-winbuilder-release : r-source-release
 	cd build; $(RCMD) -e 'target <- tempdir()' \
 			  -e 'untar("$(RPKG)_$(PKG_VERSION).tar.gz", exdir=target)' \
-			  -e 'devtools::check_win_release(pkg=file.path(target, "RBesT"))'
+			  -e 'devtools::check_win_release(pkg=file.path(target, "$(RPKG)"))'
 
 PHONY += check-winbuilder-oldrelease
 check-winbuilder-oldrelease : r-source-release
 	cd build; $(RCMD) -e 'target <- tempdir()' \
 			  -e 'untar("$(RPKG)_$(PKG_VERSION).tar.gz", exdir=target)' \
-			  -e 'devtools::check_win_oldrelease(pkg=file.path(target, "RBesT"))'
+			  -e 'devtools::check_win_oldrelease(pkg=file.path(target, "$(RPKG)"))'
 
 PHONY += check-winbuilder
 check-winbuilder : check-winbuilder-devel check-winbuilder-release check-winbuilder-oldrelease

NAMESPACE

@@ -147,6 +147,7 @@ export(mn2gamma)
 export(mn2norm)
 export(ms2beta)
 export(ms2gamma)
+export(msr2mvnorm)
 export(oc1S)
 export(oc1Sdecision)
 export(oc2S)
@@ -159,10 +160,12 @@ export(postmix)
 export(preddist)
 export(qmix)
 export(qmixdiff)
+export(read_mix_json)
 export(rmix)
 export(rmixdiff)
 export(robustify)
 export(sigma)
+export(write_mix_json)
 import(Formula)
 import(Rcpp)
 import(abind)
@@ -175,6 +178,9 @@ import(rstantools)
 import(stats)
 importFrom(RcppParallel,CxxFlags)
 importFrom(RcppParallel,RcppParallelLibs)
+importFrom(jsonlite,fromJSON)
+importFrom(jsonlite,toJSON)
+importFrom(lifecycle,deprecated)
 importFrom(matrixStats,colLogSumExps)
 importFrom(matrixStats,colSums2)
 importFrom(matrixStats,logSumExp)

NEWS.md

@@ -1,3 +1,16 @@
+# RBesT 1.8-2 - April 25th, 2025
+
+## Enhancements
+
+* Add experimental `write_mix_json` and `read_mix_json` functions
+  which write and read mixture objects as JSON to/from files.
+* Add `"msr"` parametrization to `mixmvnorm` which allows to consturct
+  multi-variate normal mixtures in a parametrization using the mean
+  vector, standard deviations and the correlations. This functionality
+  is based on the new `msr2mvnorm` utility function.
+* Allow standard deviations of 0 for components of a multi-variate
+  normal mixture.
+
 # RBesT 1.8-1 - January 20th, 2025
 
 ## Bugfixes

R/AS.R

@@ -12,7 +12,7 @@
 #'   \item{n}{study size}
 #'   \item{r}{number of events}
 #' }
-#' @references Baeten D. et. al, \emph{The Lancet}, 2013, (382), 9906, p 1705
+#' @references Baeten D. et. al, *The Lancet*, 2013, (382), 9906, p 1705
 #'
 #' @template example-start
 #' @examples

R/BinaryExactCI.R

@@ -6,7 +6,7 @@
 #' @param r Number of success or responder
 #' @param n Sample size
 #' @param alpha  confidence level
-#' @param drop Determines if \code{\link{drop}} will be called on the result
+#' @param drop Determines if [drop()] will be called on the result
 #'
 #' @details
 #' Confidence intervals are obtained by a procedure first given in
@@ -33,7 +33,10 @@ BinaryExactCI <- function(r, n, alpha = 0.05, drop = TRUE) {
   pLow <- qbeta(Low, r, n - r + 1)
   pHigh <- qbeta(High, r + 1, n - r)
 
-  nms <- c(paste(round(100 * Low, 1), "%", sep = ""), paste(round(100 * High, 1), "%", sep = ""))
+  nms <- c(
+    paste(round(100 * Low, 1), "%", sep = ""),
+    paste(round(100 * High, 1), "%", sep = "")
+  )
 
   CI <- cbind(pLow, pHigh)
   colnames(CI) <- nms

R/EM_bmm_ab.R

@@ -1,6 +1,17 @@
 ## EM for Beta Mixture Models (BMM) with Nc components
 
-EM_bmm_ab <- function(x, Nc, mix_init, Ninit = 50, verbose = FALSE, Niter.max = 500, tol, Neps, eps = c(w = 0.005, a = 0.005, b = 0.005), constrain_gt1 = TRUE) {
+EM_bmm_ab <- function(
+  x,
+  Nc,
+  mix_init,
+  Ninit = 50,
+  verbose = FALSE,
+  Niter.max = 500,
+  tol,
+  Neps,
+  eps = c(w = 0.005, a = 0.005, b = 0.005),
+  constrain_gt1 = TRUE
+) {
   N <- length(x)
   assert_that(N + Nc >= Ninit)
 
@@ -11,12 +22,20 @@ EM_bmm_ab <- function(x, Nc, mix_init, Ninit = 50, verbose = FALSE, Niter.max =
   ## ensures proper handling during fit.
   x0 <- x == 0
   if (any(x0)) {
-    message("Detected ", sum(x0), " value(s) which are exactly 0.\nTo avoid numerical issues during EM such values are moved to smallest eps on machine.")
+    message(
+      "Detected ",
+      sum(x0),
+      " value(s) which are exactly 0.\nTo avoid numerical issues during EM such values are moved to smallest eps on machine."
+    )
     x[x0] <- .Machine$double.eps
   }
   x1 <- x == 1
   if (any(x1)) {
-    message("Detected ", sum(x1), " value(s) which are exactly 1.\nTo avoid numerical issues during EM such values are moved to one minus smallest eps on machine.")
+    message(
+      "Detected ",
+      sum(x1),
+      " value(s) which are exactly 1.\nTo avoid numerical issues during EM such values are moved to one minus smallest eps on machine."
+    )
     x[x1] <- 1 - .Machine$double.eps
   }
 
@@ -32,7 +51,10 @@ EM_bmm_ab <- function(x, Nc, mix_init, Ninit = 50, verbose = FALSE, Niter.max =
     ## abmEst[,1] <- 1/Nc
     ## assume that the sample is ordered randomly
     ind <- seq(1, N - Nc, length = Ninit)
-    knnInit <- list(mu = matrix(0, nrow = Nc, ncol = 1), p = rep(1 / Nc, times = Nc))
+    knnInit <- list(
+      mu = matrix(0, nrow = Nc, ncol = 1),
+      p = rep(1 / Nc, times = Nc)
+    )
     for (k in seq(Nc)) {
       knnInit$mu[k, 1] <- mean(x[ind + k - 1])
     }
@@ -52,7 +74,11 @@ EM_bmm_ab <- function(x, Nc, mix_init, Ninit = 50, verbose = FALSE, Niter.max =
     cmin <- which.min(KNN$p)
     muInit[cmin] <- sum(KNN$p * KNN$center)
     ## muInit[cmin] <- mean(x) ## could be considered here
-    nInit[cmin] <- pmax(muInit[cmin] * (1 - muInit[cmin]) / var(x) - 1, 1, na.rm = TRUE)
+    nInit[cmin] <- pmax(
+      muInit[cmin] * (1 - muInit[cmin]) / var(x) - 1,
+      1,
+      na.rm = TRUE
+    )
     ## Nmax <- max(2, max(nInit))
     ## ensure n is positive for each cluster; if this is not the
     ## case, sample uniformly from the range of n we have
@@ -126,7 +152,11 @@ EM_bmm_ab <- function(x, Nc, mix_init, Ninit = 50, verbose = FALSE, Niter.max =
   traceMix <- list()
   traceLli <- c()
   Dlli <- Inf
-  runMixPar <- array(-Inf, dim = c(Neps, 3, Nc), dimnames = list(NULL, rownames(mixEstPar), NULL))
+  runMixPar <- array(
+    -Inf,
+    dim = c(Neps, 3, Nc),
+    dimnames = list(NULL, rownames(mixEstPar), NULL)
+  )
   runOrder <- 0:(Neps - 1)
   Npar <- Nc + 2 * Nc
   if (Nc == 1) Npar <- Npar - 1
@@ -160,10 +190,14 @@ EM_bmm_ab <- function(x, Nc, mix_init, Ninit = 50, verbose = FALSE, Niter.max =
 
       trig_n <- trigamma(n)
 
-      grad1 <- 2 * sqTerm1 * (trigamma(ab[1]) * ab[1] - trig_n * ab[1]) -
+      grad1 <- 2 *
+        sqTerm1 *
+        (trigamma(ab[1]) * ab[1] - trig_n * ab[1]) -
         2 * sqTerm2 * (trig_n * ab[1])
 
-      grad2 <- -2 * sqTerm1 * (trig_n * ab[2]) +
+      grad2 <- -2 *
+        sqTerm1 *
+        (trig_n * ab[2]) +
         2 * sqTerm2 * (trigamma(ab[2]) * ab[2] - trig_n * ab[2])
 
       c(grad1, grad2)
@@ -182,7 +216,14 @@ EM_bmm_ab <- function(x, Nc, mix_init, Ninit = 50, verbose = FALSE, Niter.max =
     a <- mixEst[2, ]
     b <- mixEst[3, ]
     ## lli <- sweep( sweep(Lx, 2, a - 1, "*", check.margin=FALSE) + sweep(LxC, 2, b - 1, "*", check.margin=FALSE), 2, log(w) + lgamma(a + b) - lgamma(a) - lgamma(b), "+", check.margin=FALSE)
-    lli <- sweep(sweep(Lx, 2, a - 1, "*", check.margin = FALSE) + sweep(LxC, 2, b - 1, "*", check.margin = FALSE), 2, log(w) - lbeta(a, b), "+", check.margin = FALSE)
+    lli <- sweep(
+      sweep(Lx, 2, a - 1, "*", check.margin = FALSE) +
+        sweep(LxC, 2, b - 1, "*", check.margin = FALSE),
+      2,
+      log(w) - lbeta(a, b),
+      "+",
+      check.margin = FALSE
+    )
     ## lli <- t(matrix(log(mixEst[1,]) + dbeta(xRep, mixEst[2,], mixEst[3,], log=TRUE), nrow=Nc))
 
     ## ensure that the log-likelihood does not go out of numerical
@@ -202,15 +243,36 @@ EM_bmm_ab <- function(x, Nc, mix_init, Ninit = 50, verbose = FALSE, Niter.max =
       Dlli <- (traceLli[iter + 1] - traceLli[iter - 1]) / 2
     }
     if (Nc > 1) {
-      smean <- apply(runMixPar[order(runOrder), , , drop = FALSE], c(2, 3), function(x) mean(abs(diff(x))))
+      smean <- apply(
+        runMixPar[order(runOrder), , , drop = FALSE],
+        c(2, 3),
+        function(x) mean(abs(diff(x)))
+      )
       eps.converged <- sum(sweep(smean, 1, eps, "-") < 0)
     } else {
-      smean <- apply(runMixPar[order(runOrder), -1, , drop = FALSE], c(2, 3), function(x) mean(abs(diff(x))))
+      smean <- apply(
+        runMixPar[order(runOrder), -1, , drop = FALSE],
+        c(2, 3),
+        function(x) mean(abs(diff(x)))
+      )
       eps.converged <- sum(sweep(smean, 1, eps[-1], "-") < 0)
     }
     if (is.na(eps.converged)) eps.converged <- 0
     if (verbose) {
-      message("Iteration ", iter, ": log-likelihood = ", lliCur, "; Dlli = ", Dlli, "; converged = ", eps.converged, " / ", Npar, "\n", sep = "")
+      message(
+        "Iteration ",
+        iter,
+        ": log-likelihood = ",
+        lliCur,
+        "; Dlli = ",
+        Dlli,
+        "; converged = ",
+        eps.converged,
+        " / ",
+        Npar,
+        "\n",
+        sep = ""
+      )
     }
     if (checkTol & Dlli < tol) {
       break
@@ -250,7 +312,13 @@ EM_bmm_ab <- function(x, Nc, mix_init, Ninit = 50, verbose = FALSE, Niter.max =
       ## Default would be Nelder-Mead
       Lest <- optim(theta, bmm_ml(c1[i], c2[i]))
       if (Lest$convergence != 0 & Lest$value > 1E-4) {
-        warning("Warning: Component", i, "in iteration", iter, "had convergence problems!")
+        warning(
+          "Warning: Component",
+          i,
+          "in iteration",
+          iter,
+          "had convergence problems!"
+        )
       }
       if (constrain_gt1) {
         mixEst[2:3, i] <- 1 + pmax(exp(Lest$par), c(1E-8, 1E-8))
@@ -297,6 +365,11 @@ EM_bmm_ab <- function(x, Nc, mix_init, Ninit = 50, verbose = FALSE, Niter.max =
 
 #' @export
 print.EMbmm <- function(x, ...) {
-  cat("EM for Beta Mixture Model\nLog-Likelihood = ", logLik(x), "\n\n", sep = "")
+  cat(
+    "EM for Beta Mixture Model\nLog-Likelihood = ",
+    logLik(x),
+    "\n\n",
+    sep = ""
+  )
   NextMethod()
 }