Add mr.mash-rss

DESCRIPTION

@@ -1,8 +1,8 @@
 Encoding: UTF-8
 Type: Package
 Package: mr.mash.alpha
-Version: 0.2.32
-Date: 2023-12-11
+Version: 0.3.22
+Date: 2024-04-30
 Title: Multiple Regression with Multivariate Adaptive Shrinkage
 Description: Provides an implementation of methods for multivariate 
     multiple regression with adaptive shrinkage priors.
@@ -11,6 +11,7 @@ Authors@R: c(person("Fabio","Morgante",role=c("cre","aut"),
                    email="[email protected]"),
              person("Jason","Willwerscheid",role="ctb"),
              person("Gao","Wang",role="ctb"),
+             person("Deborah","Kunkel",role="aut"),
              person("Peter","Carbonetto",role="aut"),
              person("Matthew","Stephens",role="aut"))
 License: MIT + file LICENSE
@@ -19,13 +20,13 @@ Imports:
     stats,
     Rcpp (>= 1.0.7),
     RcppParallel (>= 5.1.5),
-    MBSP (>= 3.0),
     mvtnorm,
     matrixStats,
     mashr (>= 0.2.73),
     ebnm,
-    flashier (>= 0.2.50),
-    parallel
+    flashier (>= 1.0.7),
+    parallel,
+    Rfast
 Suggests:
     testthat,
     varbvs,
@@ -35,5 +36,4 @@ LinkingTo:
     Rcpp,
     RcppArmadillo (>= 0.10.4.0.0),
     RcppParallel
-VignetteBuilder: knitr
-RoxygenNote: 7.2.2
+RoxygenNote: 7.2.3

LICENSE

@@ -1,3 +1,3 @@
-YEAR: 2020-2023
-COPYRIGHT HOLDER: Fabio Morgante, Jason Willwerscheid, Gao Wang, Peter Carbonetto, Matthew Stephens
-ORGANIZATION: Fabio Morgante, Jason Willwerscheid, Gao Wang, Peter Carbonetto, Matthew Stephens
+YEAR: 2020-2024
+COPYRIGHT HOLDER: Fabio Morgante, Jason Willwerscheid, Gao Wang, Deborah Kunkel, Peter Carbonetto, Matthew Stephens
+ORGANIZATION: Fabio Morgante, Jason Willwerscheid, Gao Wang, Deborah Kunkel, Peter Carbonetto, Matthew Stephens

LICENSE.md

@@ -1,5 +1,5 @@
-Copyright (c) 2020-2023, Fabio Morgante, Jason Willwerscheid, Gao Wang, 
-Peter Carbonetto, Matthew Stephens.
+Copyright (c) 2020-2024, Fabio Morgante, Jason Willwerscheid, Gao Wang, 
+Deborah Kunkel, Peter Carbonetto, Matthew Stephens.
 All rights reserved. 
 
 Redistribution and use in source and binary forms, with or without 
@@ -11,9 +11,9 @@ modification, are permitted provided that the following conditions are met:
    notice, this list of conditions and the following disclaimer in the 
    documentation and/or other materials provided with the distribution. 
  * Neither the name of Fabio Morgante, Jason Willwerscheid, Gao Wang, 
-   Peter Carbonetto, Matthew Stephens, nor the names of its contributors 
-   may be used to endorse or promote products derived from this software 
-   without specific prior written permission. 
+   Deborah Kunkel, Peter Carbonetto, Matthew Stephens, nor the names of 
+   its contributors may be used to endorse or promote products derived 
+   from this software without specific prior written permission. 
 
 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 

NAMESPACE

@@ -1,21 +1,26 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(coef,mr.mash)
+S3method(coef,mr.mash.rss)
 S3method(predict,mr.mash)
+S3method(predict,mr.mash.rss)
 export(autoselect.mixsd)
 export(coef.mr.mash)
+export(coef.mr.mash.rss)
 export(compute_canonical_covs)
 export(compute_data_driven_covs)
 export(compute_univariate_sumstats)
 export(expand_covs)
 export(mr.mash)
+export(mr.mash.rss)
 export(predict.mr.mash)
+export(predict.mr.mash.rss)
 export(simulate_mr_mash_data)
-importFrom(MBSP,matrix_normal)
 importFrom(Rcpp,evalCpp)
 importFrom(RcppParallel,RcppParallelLibs)
 importFrom(RcppParallel,defaultNumThreads)
 importFrom(RcppParallel,setThreadOptions)
+importFrom(Rfast,is.symmetric)
 importFrom(ebnm,ebnm_normal)
 importFrom(ebnm,ebnm_normal_scale_mixture)
 importFrom(flashier,flash)
@@ -38,4 +43,5 @@ importFrom(stats,cov)
 importFrom(stats,cov2cor)
 importFrom(stats,lm)
 importFrom(stats,predict)
+importFrom(stats,rnorm)
 useDynLib(mr.mash.alpha)

R/RcppExports.R

@@ -9,6 +9,10 @@ impute_missing_Y_rcpp <- function(Y, mu, Vinv, miss, non_miss) {
     .Call('_mr_mash_alpha_impute_missing_Y_rcpp', PACKAGE = 'mr.mash.alpha', Y, mu, Vinv, miss, non_miss)
 }
 
+inner_loop_general_rss_rcpp <- function(n, XtX, XtY, XtRbar, mu1, V, Vinv, w0, S0, precomp_quants_list, standardize, compute_ELBO, update_V, update_order, eps, nthreads) {
+    .Call('_mr_mash_alpha_inner_loop_general_rss_rcpp', PACKAGE = 'mr.mash.alpha', n, XtX, XtY, XtRbar, mu1, V, Vinv, w0, S0, precomp_quants_list, standardize, compute_ELBO, update_V, update_order, eps, nthreads)
+}
+
 scale_rcpp <- function(M, a, b) {
     .Call('_mr_mash_alpha_scale_rcpp', PACKAGE = 'mr.mash.alpha', M, a, b)
 }
@@ -25,3 +29,7 @@ compute_logbf_rcpp <- function(X, Y, V, Vinv, w0, S0, precomp_quants_list, stand
     .Call('_mr_mash_alpha_compute_logbf_rcpp', PACKAGE = 'mr.mash.alpha', X, Y, V, Vinv, w0, S0, precomp_quants_list, standardize, eps, nthreads)
 }
 
+compute_logbf_rss_rcpp <- function(n, XtY, V, Vinv, w0, S0, precomp_quants_list, standardize, eps, nthreads) {
+    .Call('_mr_mash_alpha_compute_logbf_rss_rcpp', PACKAGE = 'mr.mash.alpha', n, XtY, V, Vinv, w0, S0, precomp_quants_list, standardize, eps, nthreads)
+}
+

R/bayes_reg_mv_rss.R

@@ -0,0 +1,111 @@
+# Bayesian multivariate regression with mixture-of-normals prior
+# (mixture weights w0 and covariance matrices S0) with standardized X.
+#
+# The outputs are: the log-Bayes factor (logbf), the posterior assignment probabilities
+# (w1), the posterior mean of the coefficients given that all the
+# coefficients are not nonzero (mu1), and the posterior covariance of
+# the coefficients given that all the coefficients are not zero (S1).
+bayes_mvr_mix_standardized_X_rss <- function (n, xtY, w0, S0, S, S1, SplusS0_chol, S_chol, eps) {
+
+
+  # Get the number of conditions (r), the number of mixture
+  # components (K).
+  r <- length(xtY)
+  K <- length(S0)
+
+  # Compute the least-squares estimate.
+  b <- xtY/(n-1)
+
+  # Compute the Bayes factors and posterior statistics separately for
+  # each mixture component.
+  bayes_mvr_ridge_lapply <- function(i){
+    bayes_mvr_ridge_standardized_X(b, S0[[i]], S, S1[[i]], SplusS0_chol[[i]], S_chol)
+  }
+  out <- lapply(1:K, bayes_mvr_ridge_lapply)
+
+  # Compute the posterior assignment probabilities for the latent
+  # indicator variable.
+  logbf <- sapply(out,function (x) x$logbf)
+  w1    <- softmax(logbf + log(w0 + eps))
+
+  # Compute the posterior mean (mu1_mix) and covariance (S1_mix) of the
+  # regression coefficients.
+  A   <- matrix(0,r,r)
+  mu1_mix <- rep(0,r)
+  for (k in 1:K) {
+    wk  <- w1[k]
+    muk <- out[[k]]$mu1
+    Sk  <- S1[[k]]
+    mu1_mix <- mu1_mix + wk*muk
+    A   <- A   + wk*(Sk + tcrossprod(muk))
+  }
+  S1_mix <- A - tcrossprod(mu1_mix)
+
+  # Compute the log-Bayes factor for the mixture as a linear combination of the
+  # individual BFs foreach mixture component.
+  u     <- max(logbf)
+  logbf_mix <- u + log(sum(w0 * exp(logbf - u)))
+
+  # Return the log-Bayes factor for the mixture (logbf), the posterior assignment probabilities (w1), the
+  # posterior mean of the coefficients (mu1), and the posterior
+  # covariance of the coefficients (S1).
+  return(list(logbf = logbf_mix,w1 = w1,mu1 = mu1_mix,S1 = S1_mix))
+}
+
+
+# Bayesian multivariate regression with mixture-of-normals prior
+# (mixture weights w0 and covariance matrices S0) and centered X
+#
+# The outputs are: the log-Bayes factor (logbf), the posterior assignment probabilities
+# (w1), the posterior mean of the coefficients given that all the
+# coefficients are not nonzero (mu1), and the posterior covariance of
+# the coefficients given that all the coefficients are not zero (S1).
+bayes_mvr_mix_centered_X_rss <- function (xtY, V, w0, S0, xtx, Vinv, V_chol, d, QtimesV_chol, eps) {
+
+  # Get the number of variables (n) and the number of mixture
+  # components (k).
+  r <- length(xtY)
+  K <- length(S0)
+
+  # Compute the least-squares estimate and covariance.
+  b <- xtY/xtx
+  S <- V/xtx
+
+  # Compute quantities needed for bayes_mvr_ridge_centered_X()
+  S_chol <- V_chol/sqrt(xtx)
+
+  # Compute the Bayes factors and posterior statistics separately for
+  # each mixture component.
+  bayes_mvr_ridge_lapply <- function(i){
+    bayes_mvr_ridge_centered_X(V, b, S, S0[[i]], xtx, Vinv, V_chol, S_chol, d[[i]], QtimesV_chol[[i]])
+  }
+  out <- lapply(1:K, bayes_mvr_ridge_lapply)
+
+  # Compute the posterior assignment probabilities for the latent
+  # indicator variable.
+  logbf <- sapply(out,function (x) x$logbf)
+  w1    <- softmax(logbf + log(w0 + eps))
+
+  # Compute the posterior mean (mu1_mix) and covariance (S1_mix) of the
+  # regression coefficients.
+  A   <- matrix(0,r,r)
+  mu1_mix <- rep(0,r)
+  for (k in 1:K) {
+    wk  <- w1[k]
+    muk <- out[[k]]$mu1
+    Sk  <- out[[k]]$S1
+    mu1_mix <- mu1_mix + wk*muk
+    A   <- A   + wk*(Sk + tcrossprod(muk))
+  }
+  S1_mix <- A - tcrossprod(mu1_mix)
+
+  # Compute the log-Bayes factor for the mixture as a linear combination of the
+  # individual BFs foreach mixture component.
+  u     <- max(logbf)
+  logbf_mix <- u + log(sum(w0 * exp(logbf - u)))
+
+  # Return the log-Bayes factor for the mixture (logbf), the posterior assignment probabilities (w1), the
+  # posterior mean of the coefficients (mu1), and the posterior
+  # covariance of the coefficients (S1).
+  return(list(logbf = logbf_mix,w1 = w1,mu1 = mu1_mix,S1 = S1_mix))
+}

R/elbo_rss.R

@@ -0,0 +1,26 @@
+###Compute ELBO from intermediate components
+compute_ELBO_rss_fun <- function(n, RbartRbar, Vinv, ldetV, var_part_tr_wERSS, neg_KL){
+  r <- ncol(RbartRbar)
+
+  tr_wERSS <- sum(Vinv*RbartRbar) + var_part_tr_wERSS
+
+  ELBO <- -log(n)/2 - (n*r)/2*log(2*pi) - n/2 * ldetV - 0.5*(tr_wERSS) + neg_KL
+
+  return(ELBO)
+}
+
+###Compute intermediate components of the ELBO
+compute_ELBO_rss_terms <- function(var_part_tr_wERSS, neg_KL, xtRbar_j, bfit, xtx, Vinv){
+  mu1_mat <- matrix(bfit$mu1, ncol=1)
+
+  var_part_tr_wERSS <- var_part_tr_wERSS + (sum(Vinv*bfit$S1)*xtx)
+
+  mu1xtRbar_j <- mu1_mat%*%xtRbar_j
+
+  Cm <- - mu1xtRbar_j - t(mu1xtRbar_j) + tcrossprod(mu1_mat)*xtx + bfit$S1*xtx
+
+  neg_KL <- neg_KL + (bfit$logbf + 0.5*(sum(Vinv*Cm)))
+
+  return(list(var_part_tr_wERSS=var_part_tr_wERSS, neg_KL=neg_KL))
+}
+

R/misc.R

@@ -37,9 +37,6 @@ removefromcols <- function (A, b)
   t(t(A) - b)
 
 ###Function to simulate from MN distribution
-#
-#' @importFrom MBSP matrix_normal
-#'
 sim_mvr <- function (X, B, V) {
 
   # Get the number of samples (n) and conditions (m).
@@ -48,8 +45,7 @@ sim_mvr <- function (X, B, V) {
 
   # Simulate the responses, Y.
   M <- X%*%B
-  U <- diag(n)
-  Y <- matrix_normal(M, U, V)
+  Y <- matrix_normal_indep_rows(M, V)
 
   # Output the simulated responses.
   return(Y)
@@ -94,9 +90,8 @@ makePD <- function(S0, e){
 }
 
 ###Precompute quantities in any case
-precompute_quants <- function(X, V, S0, standardize, version){
+precompute_quants <- function(n, V, S0, standardize, version){
   if(standardize){
-    n <- nrow(X)
     xtx <- n-1
 
     ###Quantities that don't depend on S0
@@ -287,8 +282,8 @@ compute_cov_flash <- function(Y, error_cache = NULL){
       covar <- diag(fl$residuals_sd^2) + crossprod(t(fl$F_pm) * fsd)
     }
     if (nrow(covar) == 0) {
-    covar <- diag(ncol(Y))
-    stop("Computed covariance matrix has zero rows")
+      covar <- diag(ncol(Y))
+      stop("Computed covariance matrix has zero rows")
     }
   }, error = function(e) {
     if (!is.null(error_cache)) {
@@ -339,3 +334,13 @@ extract_missing_Y_pattern <- function(Y){
   return(list(miss=miss, non_miss=non_miss))
 }
 
+###Check whether a matrix is PSD
+check_semi_pd <- function (A, tol) {
+  attr(A,"eigen") <- eigen(A,symmetric = TRUE)
+  v <- attr(A,"eigen")$values
+  v[abs(v) < tol] = 0
+  return(list(matrix      = A,
+              status      = !any(v < 0),
+              eigenvalues = v))
+}
+

R/mr_mash.R

@@ -63,7 +63,7 @@
 #'   
 #' @param ca_update_order The order with which coordinates are
 #'   updated.  So far, "consecutive", "decreasing_logBF",
-#'   "increasing_logBF" are supported.
+#'   "increasing_logBF", "random" are supported.
 #'   
 #' @param nthreads Number of RcppParallel threads to use for the
 #'   updates. When \code{nthreads} is \code{NA}, the default number of
@@ -166,7 +166,7 @@ mr.mash <- function(X, Y, S0, w0=rep(1/(length(S0)), length(S0)), V=NULL,
                     max_iter=5000, update_w0=TRUE, update_w0_method="EM", 
                     w0_threshold=0, compute_ELBO=TRUE, standardize=TRUE, verbose=TRUE,
                     update_V=FALSE, update_V_method=c("full", "diagonal"), version=c("Rcpp", "R"), e=1e-8,
-                    ca_update_order=c("consecutive", "decreasing_logBF", "increasing_logBF"),
+                    ca_update_order=c("consecutive", "decreasing_logBF", "increasing_logBF", "random"),
                     nthreads=as.integer(NA)) {
 
   if(verbose){
@@ -325,7 +325,7 @@ mr.mash <- function(X, Y, S0, w0=rep(1/(length(S0)), length(S0)), V=NULL,
   eps <- .Machine$double.eps
 
   ###Precompute quantities
-  comps <- precompute_quants(X, V, S0, standardize, version)
+  comps <- precompute_quants(n, V, S0, standardize, version)
   if(!standardize){
     xtx <- colSums(X^2)
     comps$xtx <- xtx
@@ -356,7 +356,8 @@ mr.mash <- function(X, Y, S0, w0=rep(1/(length(S0)), length(S0)), V=NULL,
   } else if(ca_update_order=="increasing_logBF"){
     update_order <- compute_rank_variables_BFmix(X, Y, V, Vinv, w0, S0, comps, standardize, version, 
                                                  decreasing=FALSE, eps, nthreads)
-  }
+  } else if(ca_update_order=="random")
+    update_order <- sample(x=1:p, size=p)
 
   if(!Y_has_missing){
     Y_cov <- matrix(0, nrow=r, ncol=r)
@@ -418,7 +419,7 @@ mr.mash <- function(X, Y, S0, w0=rep(1/(length(S0)), length(S0)), V=NULL,
           V <- diag(diag(V))
 
         #Recompute precomputed quantities after updating V
-        comps <- precompute_quants(X, V, S0, standardize, version)
+        comps <- precompute_quants(n, V, S0, standardize, version)
         if(!standardize)
           comps$xtx <- xtx
         if(compute_ELBO || !standardize)