Merge pull request #11 from weberse2/issue-rel-1-70

CRAN 1.7-0 accepted release

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.6-7
-Date: 2023-06-26
+Version: 1.7-0
+Date: 2023-07-19
 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")
@@ -36,7 +36,8 @@ Imports:
     dplyr,
     stats,
     utils,
-    matrixStats
+    matrixStats,
+    abind
 LinkingTo:
     BH (>= 1.72.0),
     Rcpp (>= 0.12.0),

NAMESPACE

@@ -18,6 +18,7 @@ S3method(dmix,betaMix)
 S3method(dmix,default)
 S3method(dmix,gammaMix)
 S3method(dmix,gammaPoissonMix)
+S3method(dmix,mvnormMix)
 S3method(dmix,normMix)
 S3method(ess,betaMix)
 S3method(ess,default)
@@ -41,11 +42,13 @@ S3method(oc2S,normMix)
 S3method(plot,EM)
 S3method(plot,gMAP)
 S3method(plot,mix)
+S3method(plot,mvnormMix)
 S3method(pmix,betaBinomialMix)
 S3method(pmix,betaMix)
 S3method(pmix,default)
 S3method(pmix,gammaMix)
 S3method(pmix,gammaPoissonMix)
+S3method(pmix,mvnormMix)
 S3method(pmix,normMix)
 S3method(pos1S,betaMix)
 S3method(pos1S,default)
@@ -62,10 +65,12 @@ S3method(postmix,normMix)
 S3method(preddist,betaMix)
 S3method(preddist,default)
 S3method(preddist,gammaMix)
+S3method(preddist,mvnormMix)
 S3method(preddist,normMix)
 S3method(predict,gMAP)
 S3method(print,EMbmm)
 S3method(print,EMgmm)
+S3method(print,EMmvnmm)
 S3method(print,EMnmm)
 S3method(print,betaBinomialMix)
 S3method(print,betaMix)
@@ -78,34 +83,40 @@ S3method(print,gammaExpMix)
 S3method(print,gammaMix)
 S3method(print,gammaPoissonMix)
 S3method(print,mix)
+S3method(print,mvnormMix)
 S3method(print,normMix)
 S3method(qmix,betaBinomialMix)
 S3method(qmix,betaMix)
 S3method(qmix,default)
 S3method(qmix,gammaMix)
 S3method(qmix,gammaPoissonMix)
+S3method(qmix,mvnormMix)
 S3method(qmix,normMix)
 S3method(rmix,betaBinomialMix)
 S3method(rmix,betaMix)
 S3method(rmix,default)
 S3method(rmix,gammaMix)
 S3method(rmix,gammaPoissonMix)
+S3method(rmix,mvnormMix)
 S3method(rmix,normMix)
 S3method(robustify,betaMix)
 S3method(robustify,default)
 S3method(robustify,gammaMix)
 S3method(robustify,normMix)
+S3method(sigma,mvnormMix)
 S3method(sigma,normMix)
 S3method(summary,betaBinomialMix)
 S3method(summary,betaMix)
 S3method(summary,gMAP)
 S3method(summary,gMAPpred)
 S3method(summary,gammaMix)
 S3method(summary,gammaPoissonMix)
+S3method(summary,mvnormMix)
 S3method(summary,normMix)
 S3method(support,betaMix)
 S3method(support,default)
 S3method(support,gammaMix)
+S3method(support,mvnormMix)
 S3method(support,normMix)
 export("likelihood<-")
 export("sigma<-")
@@ -127,6 +138,7 @@ export(mixbeta)
 export(mixcombine)
 export(mixfit)
 export(mixgamma)
+export(mixmvnorm)
 export(mixnorm)
 export(mn2beta)
 export(mn2gamma)
@@ -151,6 +163,7 @@ export(robustify)
 export(sigma)
 import(Formula)
 import(Rcpp)
+import(abind)
 import(assertthat)
 import(checkmate)
 import(ggplot2)

NEWS.md

@@ -1,3 +1,16 @@
+# RBesT 1.7-0 - July 19th, 2023
+
+## Enhancements
+
+* Implementation multivariate normal mixture support in a first
+  version. This includes density evaluation, basic summary functions
+  and multivariate normal EM fitting. Support is not yet as complete
+  as for other densites, but will be expanded in upcoming releases.
+* Change the default for the option `constrain_gt1` of the EM for beta
+  mixtures to `TRUE`. This will be default constrain the fitted
+  parameters of each beta mixture component to be greater than unity
+  as required for finite ESS elir calculations.
+
 # RBesT 1.6-7 - June 26th, 2023
 
 ## Bug fixes

R/EM_bmm_ab.R

@@ -1,16 +1,10 @@
 
 ## 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)
-{
+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)
 
-    if(missing(constrain_gt1)) {
-        message("Since version 1.6-0 of RBesT the EM for beta mixtures constrains a > 1 & b > 1 by default.\nUse constrain_gt1=FALSE for an unconstrained fit. To avoid this message use constrain_gt1=TRUE.\nIn a future version the new default will be used without this message.")
-        constrain_gt1 <- TRUE
-    }
-
     assert_logical(constrain_gt1, any.missing=FALSE, len=1)
 
     ## check data for 0 and 1 values which are problematic, but may be

R/EM_gmm.R

@@ -1,6 +1,5 @@
 ## EM for GMM with Nc components
-EM_gmm <- function(x, Nc, mix_init, Ninit=50, verbose=FALSE, Niter.max=500, tol, Neps, eps=c(weight=0.005,alpha=0.005,beta=0.005))
-{
+EM_gmm <- function(x, Nc, mix_init, Ninit=50, verbose=FALSE, Niter.max=500, tol, Neps, eps=c(weight=0.005,alpha=0.005,beta=0.005)) {
     N <- length(x)
     assert_that(N+Nc >= Ninit)
 

R/EM_mnmm.R

@@ -4,8 +4,7 @@
 ## cluster weights and covariance matrices (taken from the knn
 ## determined clusters)
 
-EM_mnmm <- function(X, Nc, mix_init, Ninit=50, verbose=FALSE, Niter.max=500, tol, Neps, eps=c(w=0.005,m=0.005,s=0.005))
-{
+EM_mnmm <- function(X, Nc, mix_init, Ninit=50, verbose=FALSE, Niter.max=500, tol, Neps, eps=c(w=0.005,m=0.005,s=0.005)) {
     ## in case X is no matrix, interpret it as uni-variate case
     if(!is.matrix(X))
         X <- matrix(X,ncol=1)
@@ -175,8 +174,10 @@ EM_mnmm <- function(X, Nc, mix_init, Ninit=50, verbose=FALSE, Niter.max=500, tol
     ##    tauEst <- sqrt(as.vector(covEst))
     ##}
 
-    mixEst <- list(p=pEst, mean=muEst, sigma=covEst)
+    ##mixEst <- list(p=pEst, mean=muEst, sigma=covEst)
 
+    mixEst <- do.call(mixmvnorm, lapply(1:Nc, function(i) c(pEst[i], muEst[i,,drop=FALSE], matrix(covEst[i,,], Nd, Nd))))
+
     ## give further details
     attr(mixEst, "df") <- df
     attr(mixEst, "nobs") <- N
@@ -187,25 +188,22 @@ EM_mnmm <- function(X, Nc, mix_init, Ninit=50, verbose=FALSE, Niter.max=500, tol
     attr(mixEst, "tol") <- tol
     attr(mixEst, "traceLli") <- traceLli
     attr(mixEst, "traceMix") <- traceMix
+    attr(mixEst, "x") <- X
+
+    class(mixEst) <- c("EM", "EMmvnmm", "mvnormMix", "mix")
 
-    invisible(mixEst)
+    mixEst
 }
 
 ## uni-variate case
-EM_nmm <- function(X, Nc, mix_init, verbose=FALSE, Niter.max=500, tol, Neps, eps=c(w=0.005,m=0.005,s=0.005))
-{
-    res <- EM_mnmm(X=X, Nc=Nc, mix_init=mix_init, verbose=verbose, Niter.max=Niter.max, tol=tol, Neps=Neps, eps=eps)
-    convert <- function(est) suppressWarnings(do.call(mixnorm, lapply(1:Nc, function(i) c(est$p[i], est$mean[i,1], sqrt(est$sigma[i,,])))))
-    mixEst <- convert(res)
-    attr(mixEst, "df") <- attr(res, "df")
-    attr(mixEst, "nobs") <- attr(res, "nobs")
-    attr(mixEst, "lli") <- attr(res, "lli")
-    attr(mixEst, "Nc") <- attr(res, "Nc")
-    attr(mixEst, "tol") <- attr(res, "tol")
-    attr(mixEst, "traceLli") <- attr(res, "traceLli")
-    attr(mixEst, "traceMix") <-  lapply(attr(res, "traceMix"), convert)
-    attr(mixEst, "x") <- X
+EM_nmm <- function(X, Nc, mix_init, verbose=FALSE, Niter.max=500, tol, Neps, eps=c(w=0.005,m=0.005,s=0.005)) {
+    if(is.matrix(X)) {
+        assert_matrix(X, any.missing=FALSE, ncols=1)
+    }
+    mixEst <- EM_mnmm(X=X, Nc=Nc, mix_init=mix_init, verbose=verbose, Niter.max=Niter.max, tol=tol, Neps=Neps, eps=eps)
+    rownames(mixEst) <- c("w", "m", "s")
     class(mixEst) <- c("EM", "EMnmm", "normMix", "mix")
+    likelihood(mixEst) <- "normal"
     mixEst
 }
 
@@ -215,6 +213,12 @@ print.EMnmm <- function(x, ...) {
     NextMethod()
 }
 
+#' @export
+print.EMmvnmm <- function(x, ...) {
+    cat("EM for Multivariate Normal Mixture Model\nLog-Likelihood = ", logLik(x), "\n\n",sep="")
+    NextMethod()
+}
+
 ## given a vector of means and a covariance matrix for a multi-variate
 ## normal (and optionally a vector of df), return a vector of
 ## coefficients with a deterministic mapping

R/EM_plot.R

@@ -52,6 +52,7 @@
 #' @export
 plot.EM <- function(x, size=1.25, link=c("identity", "logit", "log"), ...) {
     pl <- list()
+    if(inherits(x, "mvnormMix")) stop("Multivariate normal mixtures plotting not supported.")
     pl$mixdist <- plot.mix(x, size=size, ...)
     ## in verbose mode we output EM fit diagnostics
     if(getOption("RBesT.verbose", FALSE)) {

R/RBesT-package.R

@@ -64,4 +64,5 @@
 #' @import ggplot2
 #' @import Formula
 #' @import checkmate
+#' @import abind
 NULL

R/mix.R

@@ -105,11 +105,13 @@ rmix <- function(mix, n) UseMethod("rmix")
     cl <- grep("mix$", class(mix), ignore.case=TRUE, value=TRUE)
     dl <- dlink(mix)
     if(inherits(mix, "normMix")) s <- sigma(mix)
+    if(inherits(mix, "mvnormMix")) s <- sigma(mix)
     mix <- mix[,...,drop=FALSE]
     if(rescale) mix[1,] <- mix[1,] / sum(mix[1,])
     class(mix) <- cl
     dlink(mix) <- dl
     if(inherits(mix, "normMix")) sigma(mix) <- s
+    if(inherits(mix, "mvnormMix")) sigma(mix) <- s
     mix
 }
 #' @export
@@ -126,7 +128,7 @@ rmix <- function(mix, n) UseMethod("rmix")
 ## IMPLEMENTATION DETAILS
 
 #' @export
-dmix.default <- function(mix, x, log=FALSE) "Unknown mixture"
+dmix.default <- function(mix, x, log=FALSE) stop("Unknown mixture")
 
 ## default implementation which only needs the density function;
 ## assumption is that the first argument of the density corresponds to
@@ -165,9 +167,34 @@ dmix.betaBinomialMix <- function(mix, x, log=FALSE) dmix_impl(Curry(dBetaBinomia
 #' @export
 dmix.gammaPoissonMix <- function(mix, x, log=FALSE) dmix_impl(Curry(.dnbinomAB, n=attr(mix, "n")), mix, x, log)
 
+#' @export
+dmix.mvnormMix <- function(mix, x, log=FALSE) {
+    p <- mvnormdim(mix[-1,1])
+    Nc <- ncol(mix)
+    if(is.matrix(x)) {
+        Nx <- nrow(x)
+        assert_matrix(x, any.missing=FALSE, nrows=Nx, ncols=p)
+    } else if(is.vector(x)) {
+        Nx <- 1
+        assert_vector(x, any.missing=FALSE, len=p)
+    } else {
+        stop("x must a vector or a matrix.")
+    }
+    assert_that(is.mixidentity_link(mix))
+    comp_res <- matrix(NA_real_, nrow=Nx, ncol=Nc)
+    for(i in 1:Nc) {
+        S <- mvnormsigma(mix[-1,i])
+        comp_res[,i] <- log(mix[1,i]) + mvtnorm::dmvnorm(x, mix[2:(p+1), i], sigma=S, log=TRUE, checkSymmetry=FALSE)
+    }
+    res <- matrixStats::rowLogSumExps(comp_res)
+    if(!log)
+        res <- exp(res)
+    return(res)
+}
+
 ## DISTRIBUTION FUNCTIONS
 #' @export
-pmix.default <- function(mix, q, lower.tail = TRUE, log.p=FALSE) "Unknown mixture"
+pmix.default <- function(mix, q, lower.tail = TRUE, log.p=FALSE) stop("Unknown mixture")
 
 pmix_impl <- function(dist, mix, q, lower.tail = TRUE, log.p=FALSE) {
     Nc <- ncol(mix)
@@ -217,10 +244,14 @@ pmix.betaBinomialMix <- function(mix, q, lower.tail = TRUE, log.p=FALSE) {
 #' @export
 pmix.gammaPoissonMix <- function(mix, q, lower.tail = TRUE, log.p=FALSE) pmix_impl(Curry(.pnbinomAB, n=attr(mix, "n")), mix, q, lower.tail, log.p)
 
+#' @export
+pmix.mvnormMix <- function(mix, q, ...) stop("Multivariate normal mixture cumulative density not supported.")
+
+
 ## QUANTILE FUNCTION
 
 #' @export
-qmix.default <- function(mix, p, lower.tail = TRUE, log.p=FALSE) "Unknown mixture"
+qmix.default <- function(mix, p, lower.tail = TRUE, log.p=FALSE) stop("Unknown mixture")
 
 qmix_impl <- function(quant, mix, p, lower.tail = TRUE, log.p=FALSE) {
     Nc <- ncol(mix)
@@ -286,11 +317,11 @@ qmix.betaBinomialMix  <- function(mix, p, lower.tail = TRUE, log.p=FALSE) {
 ## internal redefinition of negative binomial
 ##.qnbinomAB <- function(p, a, b, lower.tail = TRUE, log.p = FALSE ) qnbinom(p, size=a, prob=b/(b+1), lower.tail = lower.tail, log.p = log.p )
 .qnbinomAB <- function(p, a, b, n=1, lower.tail = TRUE, log.p = FALSE ) qnbinom(p, size=a, prob=b/(b+n), lower.tail = lower.tail, log.p = log.p )
-#' @export
 ##qmix.gammaPoissonMix <- function(mix, p, lower.tail = TRUE, log.p=FALSE) qmix_impl(Curry(.qnbinomAB, n=attr(mix, "n")), mix, p, lower.tail, log.p, discrete=TRUE)
 
 ## switched to numeric implementation as discretization seems to cause
 ## some trouble in the above definitions
+#' @export
 qmix.gammaPoissonMix <- function(mix, p, lower.tail = TRUE, log.p=FALSE) {
     assert_that(is.dlink_identity(attr(mix, "link")))
     ## numerical evaulation
@@ -308,11 +339,13 @@ qmix.gammaPoissonMix <- function(mix, p, lower.tail = TRUE, log.p=FALSE) {
     ind
 }
 
+#' @export
+qmix.mvnormMix <- function(mix, p, ...) stop("Multivariate normal mixture quantiles not supported.")
 
 ### RANDOM NUMBER GENERATION
 
 #' @export
-rmix.default <- function(mix, n) "Unknown mixture"
+rmix.default <- function(mix, n) stop("Unknown mixture")
 
 rmix_impl <- function(rng, mix, n) {
     ind <-  sample.int(ncol(mix), n, replace = TRUE, prob = mix[1,])
@@ -343,6 +376,29 @@ rmix.betaBinomialMix  <- function(mix, n) {
 #' @export
 rmix.gammaPoissonMix  <- function(mix, n) rmix_impl(Curry(.rnbinomAB, n=attr(mix, "n")),  mix, n)
 
+#' @export
+rmix.mvnormMix <- function(mix, n) {
+    ## sample the mixture components
+    ind <-  sample.int(ncol(mix), n, replace = TRUE, prob = mix["w",])
+    ## sort these
+    sidx <- order(ind)
+    ## ensure we can sort into the original random order
+    oidx <- seq(1, n)[sidx]
+    sind <- ind[sidx]
+    ## now sind[oidx] == ind
+    ## count how many times we need to sample which component
+    r <- rle(sind)
+    p <- mvnormdim(mix[-1,1])
+    samp <- do.call(rbind,
+                    mapply(function(comp, cn) {
+                        m <- mix[2:(p+1),comp]
+                        S <- mvnormsigma(mix[-1,comp])
+                        rmvnorm(cn, m, S, checkSymmetry=FALSE)
+                    },
+                    r$values, r$lengths, SIMPLIFY=FALSE))[oidx,,drop=FALSE]
+    attr(samp, "ind") <- ind
+    samp
+}
 
 #' @export
 print.mix <- function(x, digits, ...) {

R/mixbeta.R

@@ -56,6 +56,7 @@ NULL
 #' @export
 mixbeta <- function(..., param=c("ab", "ms", "mn")) {
     mix <- mixdist3(...)
+    assert_matrix(mix, nrows=3, any.missing=FALSE)
     param <- match.arg(param)
     mix[c(2,3),] <- switch(param,
                            ab=mix[c(2,3),],

R/mixcombine.R

@@ -47,5 +47,6 @@ mixcombine <- function(..., weight, rescale=TRUE) {
     dlink(mix) <- dl
     likelihood(mix) <- lik
     if("normMix" %in% cl) sigma(mix) <- sigma(comp[[1]])
+    if("mvnormMix" %in% cl) sigma(mix) <- sigma(comp[[1]])
     mix
 }