diff --git a/DESCRIPTION b/DESCRIPTION index 30cea6e..19ee42f 100644 --- a/DESCRIPTION +++ b/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="fabiom@clemson.edu"), 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 diff --git a/LICENSE b/LICENSE index ffad13d..ebf7e16 100644 --- a/LICENSE +++ b/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 diff --git a/LICENSE.md b/LICENSE.md index 1d117e6..867475e 100644 --- a/LICENSE.md +++ b/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 diff --git a/NAMESPACE b/NAMESPACE index 972d8f2..391bc1f 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/RcppExports.R b/R/RcppExports.R index e418778..bae2e1d 100644 --- a/R/RcppExports.R +++ b/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) +} + diff --git a/R/bayes_reg_mv_rss.R b/R/bayes_reg_mv_rss.R new file mode 100644 index 0000000..dab1eed --- /dev/null +++ b/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)) +} \ No newline at end of file diff --git a/R/elbo_rss.R b/R/elbo_rss.R new file mode 100644 index 0000000..d12e19a --- /dev/null +++ b/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)) +} + diff --git a/R/misc.R b/R/misc.R index e17de48..6247218 100644 --- a/R/misc.R +++ b/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)) +} + diff --git a/R/mr_mash.R b/R/mr_mash.R index 520bf30..7c503e6 100644 --- a/R/mr_mash.R +++ b/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) diff --git a/R/mr_mash_rss.R b/R/mr_mash_rss.R new file mode 100644 index 0000000..5600223 --- /dev/null +++ b/R/mr_mash_rss.R @@ -0,0 +1,606 @@ +#' @title Multiple Regression with Multivariate Adaptive Shrinkage +#' from summary data. +#' +#' @description Performs multivariate multiple regression with +#' mixture-of-normals prior. +#' +#' @param Bhat p x r matrix of regression coefficients from univariate +#' simple linear regression. +#' +#' @param Shat p x r matrix of standard errors of the regression coefficients +#' from univariate simple linear regression. +#' +#' @param Z p x r matrix of Z-scores from univariate +#' simple linear regression. +#' +#' @param R p x p correlation matrix among the variables. +#' +#' @param covY r x r covariance matrix across responses. +#' +#' @param n scalar indicating the sample size. +#' +#' @param V r x r residual covariance matrix. +#' +#' @param S0 List of length K containing the desired r x r prior +#' covariance matrices on the regression coefficients. +#' +#' @param w0 K-vector with prior mixture weights, each associated with +#' the respective covariance matrix in \code{S0}. +#' +#' @param mu1_init p x r matrix of initial estimates of the posterior +#' mean regression coefficients. These should be on the same scale as +#' the X provided. If \code{standardize=TRUE}, mu1_init will be scaled +#' appropriately after standardizing X. +#' +#' @param convergence_criterion Criterion to use for convergence check. +#' +#' @param tol Convergence tolerance. +#' +#' @param max_iter Maximum number of iterations for the optimization +#' algorithm. +#' +#' @param update_w0 If \code{TRUE}, prior weights are updated. +#' +#' @param update_w0_method Method to update prior weights. Only EM is +#' currently supported. +#' +#' @param w0_threshold Drop mixture components with weight less than this value. +#' Components are dropped at each iteration after 15 initial iterations. +#' This is done to prevent from dropping some poetentially important +#' components prematurely. +#' +#' @param update_V if \code{TRUE}, residual covariance is updated. +#' +#' @param update_V_method Method to update residual covariance. So far, +#' "full" and "diagonal" are supported. If \code{update_V=TRUE} and V +#' is not provided by the user, this option will determine how V is +#' computed (and fixed) internally from \code{mu1_init}. +#' +#' @param compute_ELBO If \code{TRUE}, ELBO is computed. +#' +#' @param standardize If \code{TRUE}, X is "standardized" using the +#' sample means and sample standard deviations. Standardizing X +#' allows a faster implementation, but the prior has a different +#' interpretation. Coefficients and covariances are returned on the +#' original scale. +#' +#' @param version Whether to use R or C++ code to perform the +#' coordinate ascent updates. +#' +#' @param verbose If \code{TRUE}, some information about the +#' algorithm's process is printed at each iteration. +#' +#' @param e A small number to add to the diagonal elements of the +#' prior matrices to improve numerical stability of the updates. +#' +#' @param ca_update_order The order with which coordinates are +#' updated. So far, "consecutive", "decreasing_logBF", +#' "increasing_logBF", "random" are supported. +#' +#' @param X_colmeans a p-vector of variable means. +#' +#' @param Y_colmeans a r-vector of response means. +#' +#' @param check_R If \code{TRUE}, R is checked to be positive semidefinite. +#' +#' @param R_tol tolerance to declare positive semi-definiteness of R. +#' +#' @param nthreads Number of RcppParallel threads to use for the +#' updates. When \code{nthreads} is \code{NA}, the default number of +#' threads is used; see +#' \code{\link[RcppParallel]{defaultNumThreads}}. This setting is +#' ignored when \code{version = "R"}. +#' +#' @return A mr.mash.rss fit, stored as a list with some or all of the +#' following elements: +#' +#' \item{mu1}{p x r matrix of posterior means for the regression +#' coeffcients.} +#' +#' \item{S1}{r x r x p array of posterior covariances for the +#' regression coeffcients.} +#' +#' \item{w1}{p x K matrix of posterior assignment probabilities to the +#' mixture components.} +#' +#' \item{V}{r x r residual covariance matrix} +#' +#' \item{w0}{K-vector with (updated, if \code{update_w0=TRUE}) prior mixture weights, each associated with +#' the respective covariance matrix in \code{S0}}. +#' +#' \item{S0}{r x r x K array of prior covariance matrices +#' on the regression coefficients}. +#' +#' \item{intercept}{r-vector containing posterior mean estimate of the +#' intercept, if \code{X_colmeans} and \code{Y_colmeans} are provided. +#' Otherwise, \code{NA} is output.} +#' +#' \item{ELBO}{Evidence Lower Bound (ELBO) at last iteration.} +#' +#' \item{progress}{A data frame including information regarding +#' convergence criteria at each iteration.} +#' +#' \item{converged}{\code{TRUE} or \code{FALSE}, indicating whether +#' the optimization algorithm converged to a solution within the chosen tolerance +#' level.} +#' +#' \item{Y}{n x r matrix of responses at last iteration (only relevant when missing values +#' are present in the input Y).} +#' +#' @examples +#' ###Set seed +#' set.seed(123) +#' +#' ###Simulate X and Y +#' ##Set parameters +#' n <- 1000 +#' p <- 100 +#' p_causal <- 20 +#' r <- 5 +#' +#' ###Simulate data +#' out <- simulate_mr_mash_data(n, p, p_causal, r, pve=0.5, B_cor=1, +#' B_scale=1, X_cor=0, X_scale=1, V_cor=0) +#' +#' ###Split the data in training and test sets +#' Ytrain <- out$Y[-c(1:200), ] +#' Xtrain <- out$X[-c(1:200), ] +#' Ytest <- out$Y[c(1:200), ] +#' Xtest <- out$X[c(1:200), ] +#' +#' ###Specify the covariance matrices for the mixture-of-normals prior. +#' univ_sumstats <- compute_univariate_sumstats(Xtrain, Ytrain, +#' standardize=TRUE, standardize.response=FALSE) +#' grid <- autoselect.mixsd(univ_sumstats, mult=sqrt(2))^2 +#' S0 <- compute_canonical_covs(ncol(Ytrain), singletons=TRUE, +#' hetgrid=c(0, 0.25, 0.5, 0.75, 1)) +#' S0 <- expand_covs(S0, grid, zeromat=TRUE) +#' +#' ###Fit mr.mash +#' covY <- cov(Ytrain) +#' corX <- cor(Xtrain) +#' n_train <- nrow(Ytrain) +#' fit <- mr.mash.rss(Bhat=univ_sumstats$Bhat, Shat=univ_sumstats$Shat, S0=S0, +#' covY=covY, R=corX, n=n_train, V=covY, update_V=TRUE, +#' X_colmeans=colMeans(Xtrain), Y_colmeans=colMeans(Ytrain)) +#' +#' # Predict the multivariate outcomes in the test set using the fitted model. +#' Ytest_est <- predict(fit,Xtest) +#' plot(Ytest_est,Ytest,pch = 20,col = "darkblue",xlab = "true", +#' ylab = "predicted") +#' abline(a = 0,b = 1,col = "magenta",lty = "dotted") +#' +#' @importFrom stats cov +#' @importFrom Rfast is.symmetric +#' @importFrom RcppParallel defaultNumThreads +#' @importFrom RcppParallel setThreadOptions +#' +#' @export +#' +mr.mash.rss <- function(Bhat, Shat, Z, R, covY, n, S0, w0=rep(1/(length(S0)), length(S0)), V, + mu1_init=NULL, tol=1e-4, convergence_criterion=c("mu1", "ELBO"), + max_iter=5000, update_w0=TRUE, update_w0_method="EM", + w0_threshold=0, compute_ELBO=TRUE, standardize=FALSE, 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", "random"), + X_colmeans=NULL, Y_colmeans=NULL, check_R=TRUE, R_tol=1e-08, + nthreads=as.integer(NA)) { + + if(verbose){ + tic <- Sys.time() + cat("Processing the inputs... ") + } + + # CHECK AND PROCESS INPUTS + # ------------------------ + ###Select method to check for convergence (if not specified by user, mu1 + ###will be used) + convergence_criterion <- match.arg(convergence_criterion) + + ###Select method to update the weights (if not specified by user, EM + ###will be used) + update_w0_method <- match.arg(update_w0_method) + + ###Select method to update the residual covariance (if not specified by user, full + ###will be used) + update_V_method <- match.arg(update_V_method) + + ###Select version of the inner loop (if not specified by user, Rcpp + ###will be used) + version <- match.arg(version) + + ###Select ordering of the coordinate ascent updates (if not specified by user, + ###consecutive will be used + ca_update_order <- match.arg(ca_update_order) + + ###Initialize the RcppParallel multithreading using a pre-specified number + ###of threads, or using the default number of threads when nthreads is NA. + if(version=="Rcpp"){ + if (is.na(nthreads)) { + setThreadOptions() + nthreads <- defaultNumThreads() + } else + setThreadOptions(numThreads = nthreads) + } + + ###Check that the inputs are in the correct format + if (sum(c(missing(Z), missing(Bhat) || missing(Shat))) != 1) + stop("Please provide either Z or (Bhat, Shat), but not both.") + + if(missing(Z)){ + if(!is.matrix(Bhat)) + stop("Bhat must be a matrix.") + if(!is.matrix(Shat)) + stop("Shat must be a matrix.") + if(any(is.na(Bhat)) || any(is.na(Shat))) + stop("Bhat, Shat must not contain missing values.") + if(any(Shat <= 0)) + stop("Shat cannot have zero or negative elements.") + } else { + if(!is.matrix(Z)) + stop("Z must be a matrix.") + if(any(is.na(Z))) + stop("Z must not contain missing values.") + } + if(!is.matrix(V) || !is.symmetric(V)) + stop("V must be a symmetric matrix.") + if(!missing(covY)){ + if(!is.matrix(covY) || !is.symmetric(covY)) + stop("covY must be a symmetric matrix.") + } + if(!is.matrix(R) || !is.symmetric(R)) + stop("R must be a symmetric matrix.") + if(!is.list(S0)) + stop("S0 must be a list.") + if(!is.vector(w0)) + stop("w0 must be a vector.") + if(length(w0)<2) + stop("At least 2 mixture components must be present.") + if(abs(sum(w0) - 1) > 1e-10) + stop("Elements of w0 must sum to 1.") + if(length(S0)!=length(w0)) + stop("S0 and w0 must have the same length.") + if(!is.null(mu1_init) && !is.matrix(mu1_init)) + stop("mu1_init must be a matrix.") + if(convergence_criterion=="ELBO" && !compute_ELBO) + stop("ELBO needs to be computed with convergence_criterion=\"ELBO\".") + if(isTRUE(standardize)) + stop("standardize=TRUE has not been implemented yet. Please use standardize=FALSE.") + + + # PRE-PROCESSING STEPS + # -------------------- + + ###Compute Z scores + if(missing(Z)){ + Z <- Bhat/Shat + } + + Z[is.na(Z)] <- 0 + + ###Compute pve-adjusted Z scores + adj <- (n-1)/(Z^2 + n - 2) + Z <- sqrt(adj) * Z + + ###Obtain dimensions and store dimensions names of the inputs + p <- nrow(Z) + r <- ncol(Z) + K <- length(S0) + Z_colnames <- colnames(Z) + Z_rownames <- rownames(Z) + + ###If covariance of Y and standard errors are provided, + ###the effects are on the *original scale*. + if(!missing(Shat) & !missing(covY)){ + XtXdiag <- rowMeans(matrix(diag(covY), nrow=p, ncol=r, byrow=TRUE) * adj/(Shat^2)) + XtX <- t(R * sqrt(XtXdiag)) * sqrt(XtXdiag) + XtX <- (XtX + t(XtX))/2 + XtY <- Z * sqrt(adj) * matrix(diag(covY), nrow=p, ncol=r, byrow=TRUE) / Shat + } else { + ###The effects are on the *standardized* X, y scale. + XtX <- R*(n-1) + XtY <- Z*sqrt(n-1) + covY <- cov2cor(V) + } + + YtY <- covY*(n-1) + + ###Check whether XtX is positive semidefinite + if(check_R){ + semi_pd <- check_semi_pd(XtX, R_tol) + if (!semi_pd$status) + stop("XtX is not a positive semidefinite matrix") + } + + ###Adjust XtX and XtY if X is standardized + if(standardize){ + dXtX <- diag(XtX) + sx <- sqrt(dXtX/(n-1)) + sx[sx == 0] <- 1 + XtX <- t((1/sx) * XtX) / sx + XtY <- XtY / sx + } + + ###Add number to diagonal elements of the prior matrices (improves + ###numerical stability) + S0 <- lapply(S0, makePD, e=e) + + ###Initialize regression coefficients to 0 if not provided + if(is.null(mu1_init)){ + mu1_init <- matrix(0, nrow=p, ncol=r) + } + + ###Scale mu1_init, if X is standardized + if(standardize) + mu1_init <- mu1_init*sx + + ###Initialize mu1, S1, w1, delta_mu1, delta_ELBO, delta_conv, ELBO, iterator, progress + mu1_t <- mu1_init + delta_mu1 <- matrix(Inf, nrow=p, ncol=r) + delta_ELBO <- Inf + if(convergence_criterion=="mu1") + delta_conv <- max(delta_mu1) + else if(convergence_criterion=="ELBO") + delta_conv <- delta_ELBO + ELBO <- -Inf + t <- 0 + progress <- as.data.frame(matrix(as.numeric(NA), nrow=max_iter, ncol=3)) + colnames(progress) <- c("iter", "timing", "mu1_max.diff") + if(compute_ELBO){ + progress$ELBO_diff <- as.numeric(NA) + progress$ELBO <- as.numeric(NA) + } + + ###Compute V, if not provided by the user + # if(is.null(V)){ + # # How to do so with sumstats?? + # + # if(update_V_method=="diagonal") + # V <- diag(diag(V)) + # } + + ###Set eps + eps <- .Machine$double.eps + + ###Precompute quantities + comps <- precompute_quants(n, V, S0, standardize, version) + if(!standardize){ + comps$xtx <- diag(XtX) + } + + if(compute_ELBO || !standardize) + ###Compute inverse of V (needed for the ELBO and unstandardized X) + Vinv <- chol2inv(comps$V_chol) + else { + if(version=="R") + Vinv <- NULL + else if(version=="Rcpp") + Vinv <- matrix(0, nrow=r, ncol=r) + } + + if(compute_ELBO) + ###Compute log determinant of V (needed for the ELBO) + ldetV <- chol2ldet(comps$V_chol) + else + ldetV <- NULL + + ###Set the ordering of the coordinate ascent updates + if(ca_update_order=="consecutive"){ + update_order <- 1:p + } else if(ca_update_order=="decreasing_logBF"){ + update_order <- compute_rank_variables_BFmix_rss(n, XtY, V, Vinv, w0, S0, comps, standardize, version, + decreasing=TRUE, eps, nthreads) + } else if(ca_update_order=="increasing_logBF"){ + update_order <- compute_rank_variables_BFmix_rss(n, XtY, 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(verbose) + cat("Done!\n") + + # MAIN LOOP + # --------- + if(verbose){ + if(version=="Rcpp" && nthreads>1) + cat(sprintf("Fitting the optimization algorithm using %d RcppParallel threads... \n", nthreads)) + else + cat("Fitting the optimization algorithm... \n") + cat(" iter mu1_max.diff") + if(compute_ELBO) + cat(" ELBO_diff ELBO\n") + else + cat("\n") + } + + ###Repeat the following until convergence, or until maximum number + ###of iterations is reached. + while(delta_conv>tol){ + + ##Start timing + time1 <- proc.time() + + ##Save current estimates. + mu1_old <- mu1_t + + ##Set last value of ELBO as ELBO_old + ELBO_old <- ELBO + + ##Update iterator + t <- t+1 + + ##Exit loop if maximum number of iterations is reached + if(t>max_iter){ + warning("Max number of iterations reached. Try increasing max_iter.") + break + } + + # M-STEP + # ------ + if(t > 1){ + ##Update V if requested + if(update_V){ + V <- update_V_rss_fun(n, RbartRbar, var_part_ERSS) + if(update_V_method=="diagonal") + V <- diag(diag(V)) + + #Recompute precomputed quantities after updating V + comps <- precompute_quants(n, V, S0, standardize, version) + if(!standardize) + comps$xtx <- diag(XtX) + if(compute_ELBO || !standardize) + Vinv <- chol2inv(comps$V_chol) + if(compute_ELBO) + ldetV <- chol2ldet(comps$V_chol) + } + + ##Update w0 if requested + if(update_w0){ + w0 <- update_weights_em(w1_t) + + #Drop components with mixture weight <= w0_threshold + if(t>15 && any(w0 < w0_threshold)){ + to_keep <- which(w0 >= w0_threshold) + w0 <- w0[to_keep] + w0 <- w0/sum(w0) + S0 <- S0[to_keep] + if(length(to_keep) > 1){ + comps <- filter_precomputed_quants(comps, to_keep, standardize, version) + } else if(length(to_keep) == 1 & all((S0[[to_keep]] - (diag(nrow(S0[[to_keep]]))*e)) < eps)){ #null component is the only one left + mu1_t <- matrix(0, nrow=p, ncol=r) + S1_t <- array(0, c(r, r, p)) + w1_t <- matrix(1, nrow=p, ncol=1) + warning("Only the null component is left. Estimated coefficients are set to 0.") + break + } else { #some other component is the only one left + stop("Only one component (different from the null) left. Consider lowering w0_threshold.") + } + } + } + } + + # E-STEP + # ------ + ###Update variational parameters + ups <- mr_mash_update_general_rss(n=n, XtX=XtX, XtY=XtY, YtY=YtY, + mu1_t=mu1_t, V=V, + Vinv=Vinv, ldetV=ldetV, w0=w0, S0=S0, + precomp_quants=comps, + compute_ELBO=compute_ELBO, + standardize=standardize, + update_V=update_V, version=version, + update_order=update_order, eps=eps, + nthreads=nthreads) + mu1_t <- ups$mu1_t + S1_t <- ups$S1_t + w1_t <- ups$w1_t + if(compute_ELBO) + ELBO <- ups$ELBO + if(update_V) + var_part_ERSS <- ups$var_part_ERSS + if(compute_ELBO || update_V) + RbartRbar <- ups$RbartRbar + + ##End timing + time2 <- proc.time() + + ##Compute difference in mu1 and ELBO between two successive iterations, + ##and assign the requested criterion to delta_conv + delta_mu1 <- abs(mu1_t - mu1_old) + delta_ELBO <- ELBO - ELBO_old + if(convergence_criterion=="mu1") + delta_conv <- max(delta_mu1) + else if(convergence_criterion=="ELBO") + delta_conv <- delta_ELBO + + ##Update progress data.frame + progress[t, c(1:3)] <- c(t, time2["elapsed"] - time1["elapsed"], max(delta_mu1)) + if(compute_ELBO) + progress[t, c(4, 5)] <- c(delta_ELBO, ELBO) + + if(verbose){ + ##Print out useful info + cat(sprintf("%4d %9.2e", t, max(delta_mu1))) + if(compute_ELBO) + cat(sprintf(" %9.2e %0.20e\n", delta_ELBO, ELBO)) + else + cat("\n") + } + } + + ###Record convergence status + if(t>max_iter) + converged <- FALSE + else + converged <- TRUE + + if(verbose){ + cat("Done!\n") + cat("Processing the outputs... ") + } + + # POST-PROCESSING STEPS + # -------------------- + if(standardize){ + ###Rescale posterior means and covariance of coefficients. In the + ###context of predicting Y, this rescaling is equivalent to + ###rescaling each column j of a given matrix, Xnew, by sx[j]. + post_rescaled <- rescale_post_mean_covar_fast(mu1_t, S1_t, sx) + mu1_t <- post_rescaled$mu1_orig + S1_t <- post_rescaled$S1_orig + } + + ###Compute posterior mean estimate of intercept. Note that when + ###columns of X are standardized, the intercept should be computed + ###with respect to the *rescaled* coefficients to recover the + ###correct fitted values. This is why this is done after rescaling + ###the coefficients above. + if(!is.null(X_colmeans) & !is.null(Y_colmeans)){ + intercept <- drop(Y_colmeans - X_colmeans %*% mu1_t) + names(intercept) <- Z_colnames + } + + ###Assign names to outputs dimensions + S0_names <- names(S0) + rownames(mu1_t) <- Z_rownames + colnames(mu1_t) <- Z_colnames + dimnames(S1_t)[[1]] <- Z_colnames + dimnames(S1_t)[[2]] <- Z_colnames + dimnames(S1_t)[[3]] <- Z_rownames + S0 <- lapply(S0, function(x){rownames(x) <- colnames(x) <- Z_colnames; return(x)}) + rownames(w1_t) <- Z_rownames + colnames(w1_t) <- S0_names + names(w0) <- S0_names + rownames(V) <- Z_colnames + colnames(V) <- Z_colnames + + ###Remove unused rows of progress + progress <- progress[rowSums(is.na(progress)) != ncol(progress), ] + + ###Return the posterior assignment probabilities (w1), the + ###posterior mean of the coefficients (mu1), and the posterior + ###covariance of the coefficients (S1), the residual covariance (V), + ###the prior weights (w0), the intercept (intercept), the fitted values (fitted), + ###and the progress data frame (progress), the prior covariance (S0), convergence + ###status, the covariance of the fitted values (G), the proportion of variance explained (pve), + ###the Evidence Lower Bound (ELBO; if computed) and imputed responses (Y; if + ###missing values were present). + out <- list(mu1=mu1_t, S1=S1_t, w1=w1_t, V=V, w0=w0, S0=simplify2array_custom(S0), + intercept=NA, progress=progress, converged=converged) + if(compute_ELBO) + ###Append ELBO to the output + out$ELBO <- ELBO + if(!is.null(X_colmeans) & !is.null(Y_colmeans)) + out$intercept <- intercept + + class(out) <- c("mr.mash.rss", "list") + + if(verbose){ + cat("Done!\n") + toc <- Sys.time() + cat("mr.mash.rss successfully executed in", difftime(toc, tic, units="mins"), + "minutes!\n") + } + + return(out) +} diff --git a/R/mr_mash_rss_predict.R b/R/mr_mash_rss_predict.R new file mode 100644 index 0000000..b5aa144 --- /dev/null +++ b/R/mr_mash_rss_predict.R @@ -0,0 +1,47 @@ +#' @title Predict future observations from mr.mash.rss fit. +#' +#' @param object a mr.mash.rss fit. +#' +#' @param newx a new value for X at which to do predictions. +#' +#' @param \dots Additional arguments (not used). +#' +#' @return Matrix of predicted values. +#' +#' @export +#' @export predict.mr.mash.rss +#' +predict.mr.mash.rss <- function(object, newx, ...){ + if(!is.matrix(newx)) + stop("X must be a matrix.") + if (any(is.na(newx))) + stop("X must not contain missing values.") + + if(any(!is.na(object$intercept))) + return(with(object,addtocols(newx %*% mu1,intercept))) + else + return(with(object,newx %*% mu1)) +} + +#' @title Extract coefficients from mr.mash.rss fit. +#' +#' @param object a mr.mash fit. +#' +#' @param \dots Other arguments (not used). +#' +#' @return p x r or (p+1) x r matrix of coefficients, +#' depending on whether an intercept was computed. +#' +#' @export +#' @export coef.mr.mash.rss +#' +coef.mr.mash.rss <- function(object, ...){ + if(any(!is.na(object$intercept))){ + coeffs <- rbind(object$intercept, object$mu1) + rownames(coeffs)[1] <- "(Intercept)" + } else { + coeffs <- object$mu1 + } + + return(coeffs) +} diff --git a/R/mr_mash_rss_updates.R b/R/mr_mash_rss_updates.R new file mode 100644 index 0000000..6282d33 --- /dev/null +++ b/R/mr_mash_rss_updates.R @@ -0,0 +1,184 @@ +###Update variational parameters, expected residuals, and ELBO components with or without scaling X +inner_loop_general_rss_R <- function(n, XtX, XtY, XtRbar, mu1, V, Vinv, w0, S0, ###note: V is only needed when not scaling X + precomp_quants, standardize, compute_ELBO, update_V, + update_order, eps){ + ###Create variables to store quantities + r <- ncol(mu1) + p <- nrow(mu1) + K <- length(S0) + S1 <- array(0, c(r, r, p)) + w1 <- matrix(0, nrow=p, ncol=K) + var_part_tr_wERSS <- 0 + neg_KL <- 0 + var_part_ERSS <- matrix(0, nrow=r, ncol=r) + + ##Loop through the variables + for(j in update_order){ + + if(standardize){ + xtx <- n-1 + } else { + xtx <- precomp_quants$xtx[j] + } + + #Remove j-th effect from expected residuals + XtRbar <- XtRbar + outer(XtX[, j], mu1[j, ]) + xtRbar_j <- XtRbar[j, ] + + #Run Bayesian SLR + if(standardize){ + bfit <- bayes_mvr_mix_standardized_X_rss(n, xtRbar_j, w0, S0, precomp_quants$S, precomp_quants$S1, + precomp_quants$SplusS0_chol, precomp_quants$S_chol, eps) + } else { + bfit <- bayes_mvr_mix_centered_X_rss(xtRbar_j, V, w0, S0, xtx, Vinv, + precomp_quants$V_chol, precomp_quants$d, + precomp_quants$QtimesV_chol, eps) + } + + #Update variational parameters + mu1[j, ] <- bfit$mu1 + S1[, , j] <- bfit$S1 + w1[j, ] <- bfit$w1 + + #Compute ELBO params + if(compute_ELBO){ + ELBO_parts <- compute_ELBO_rss_terms(var_part_tr_wERSS, neg_KL, xtRbar_j, bfit, xtx, Vinv) + var_part_tr_wERSS <- ELBO_parts$var_part_tr_wERSS + neg_KL <- ELBO_parts$neg_KL + } + + #Compute V params + if(update_V){ + var_part_ERSS <- compute_var_part_ERSS(var_part_ERSS, bfit, xtx) + } + + #Update expected residuals + XtRbar <- XtRbar - outer(XtX[, j], mu1[j, ]) + } + + ###Return output + if(compute_ELBO && update_V){ + return(list(mu1=mu1, S1=S1, w1=w1, var_part_tr_wERSS=var_part_tr_wERSS, neg_KL=neg_KL, var_part_ERSS=var_part_ERSS)) + } else if(compute_ELBO && !update_V){ + return(list(mu1=mu1, S1=S1, w1=w1, var_part_tr_wERSS=var_part_tr_wERSS, neg_KL=neg_KL)) + } else if(!compute_ELBO && update_V) { + return(list(mu1=mu1, S1=S1, w1=w1, var_part_ERSS=var_part_ERSS)) + } else { + return(list(mu1=mu1, S1=S1, w1=w1)) + } +} + +### Wrapper for the Rcpp function to update variational parameters, +### expected residuals, and ELBO components with or without scaling X. + +#' @importFrom Rcpp evalCpp +#' @importFrom RcppParallel RcppParallelLibs +#' @useDynLib mr.mash.alpha +#' +inner_loop_general_rss_Rcpp <- function(n, XtX, XtY, XtRbar, mu1, V, Vinv, w0, S0, precomp_quants, + standardize, compute_ELBO, update_V, update_order, + eps, nthreads){ + + out <- inner_loop_general_rss_rcpp(n, XtX, XtY, XtRbar, mu1, V, Vinv, w0, S0, precomp_quants, + standardize, compute_ELBO, update_V, update_order, + eps, nthreads) + + ###Return output + if(compute_ELBO && update_V){ + return(list(mu1=out$mu1, S1=out$S1, w1=out$w1, var_part_tr_wERSS=out$var_part_tr_wERSS, + neg_KL=out$neg_KL, var_part_ERSS=out$var_part_ERSS)) + } else if(compute_ELBO && !update_V){ + return(list(mu1=out$mu1, S1=out$S1, w1=out$w1, var_part_tr_wERSS=out$var_part_tr_wERSS, + neg_KL=out$neg_KL)) + } else if(!compute_ELBO && update_V) { + return(list(mu1=out$mu1, S1=out$S1, w1=out$w1, var_part_ERSS=out$var_part_ERSS)) + } else { + return(list(mu1=out$mu1, S1=out$S1, w1=out$w1)) + } +} + +###Wrapper of the inner loop with R or Rcpp +inner_loop_general_rss <- function(n, XtX, XtY, XtRbar, mu1, V, Vinv, w0, S0, precomp_quants, + standardize, compute_ELBO, update_V, version, + update_order, eps, nthreads){ + if(version=="R"){ + out <- inner_loop_general_rss_R(n, XtX, XtY, XtRbar, mu1, V, Vinv, w0, S0, precomp_quants, + standardize, compute_ELBO, update_V, update_order, eps) + } else if(version=="Rcpp"){ + update_order <- as.integer(update_order-1) + out <- inner_loop_general_rss_Rcpp(n, XtX, XtY, XtRbar, mu1, V, Vinv, w0, simplify2array_custom(S0), precomp_quants, + standardize, compute_ELBO, update_V, update_order, eps, nthreads) + } + + return(out) +} + + +###Perform one iteration of the outer loop with or without scaling X +mr_mash_update_general_rss <- function(n, XtX, XtY, YtY, mu1_t, V, Vinv, ldetV, w0, S0, + precomp_quants, compute_ELBO, standardize, + update_V, version, update_order, eps, + nthreads){ + + + + ##Compute expected residuals + XtRbar <- XtY - XtX %*% mu1_t + + ##Update variational parameters, expected residuals, and ELBO components + updates <- inner_loop_general_rss(n=n, XtX=XtX, XtY=XtY, XtRbar=XtRbar, mu1=mu1_t, V=V, Vinv=Vinv, w0=w0, S0=S0, + precomp_quants=precomp_quants, standardize=standardize, + compute_ELBO=compute_ELBO, update_V=update_V, version=version, + update_order=update_order, eps=eps, nthreads=nthreads) + mu1_t <- updates$mu1 + S1_t <- updates$S1 + w1_t <- updates$w1 + + out <- list(mu1_t=mu1_t, S1_t=S1_t, w1_t=w1_t) + + if(compute_ELBO || update_V){ + RbartRbar <- YtY - crossprod(mu1_t, XtY) - crossprod(XtY, mu1_t) + crossprod(mu1_t, XtX)%*%mu1_t + } + + if(compute_ELBO && update_V){ + ##Compute ELBO + var_part_tr_wERSS <- updates$var_part_tr_wERSS + neg_KL <- updates$neg_KL + out$ELBO <- compute_ELBO_rss_fun(n=n, RbartRbar=RbartRbar, Vinv=Vinv, ldetV=ldetV, + var_part_tr_wERSS=var_part_tr_wERSS, neg_KL=neg_KL) + + out$var_part_ERSS <- updates$var_part_ERSS + out$RbartRbar <- RbartRbar + + } else if(compute_ELBO && !update_V){ + ##Compute ELBO + var_part_tr_wERSS <- updates$var_part_tr_wERSS + neg_KL <- updates$neg_KL + out$ELBO <- compute_ELBO_rss_fun(n=n, RbartRbar=RbartRbar, Vinv=Vinv, ldetV=ldetV, + var_part_tr_wERSS=var_part_tr_wERSS, neg_KL=neg_KL) + + } else if(!compute_ELBO && update_V){ + out$var_part_ERSS <- updates$var_part_ERSS + out$RbartRbar <- RbartRbar + } + + return(out) +} + + +###Update V +update_V_rss_fun <- function(n, RbartRbar, var_part_ERSS){ + + ERSS <- RbartRbar + var_part_ERSS + V <- ERSS/n + + return(V) +} + + +###Update mixture weights +update_weights_em <- function(x){ + w <- colSums(x) + w <- w/sum(w) + return(w) +} diff --git a/R/simulate_demo_data.R b/R/simulate_demo_data.R index 66f02f0..78f4859 100644 --- a/R/simulate_demo_data.R +++ b/R/simulate_demo_data.R @@ -32,6 +32,7 @@ #' @param X_scale scalar indicating the diagonal value for Gamma. #' #' @param V_cor scalar indicating the positive correlation [0, 1] between residuals +#' #' #' @return A list with some or all of the #' following elements: @@ -59,13 +60,13 @@ #' mixture components each causal effect comes.} #' #' @importFrom mvtnorm rmvnorm -#' @importFrom MBSP matrix_normal #' @importFrom matrixStats colVars #' #' @export #' #' #' @examples +#' set.seed(1) #' dat <- simulate_mr_mash_data(n=50, p=40, p_causal=20, r=5, #' r_causal=list(1:2, 3:4), intercepts=rep(1, 5), #' pve=0.2, B_cor=c(0, 1), B_scale=c(0.5, 1), @@ -100,6 +101,7 @@ simulate_mr_mash_data <- function(n, p, p_causal, r, r_causal=list(1:r), interce Sigma[[i]] <- matrix(Sigma_offdiag, nrow=r_mix_length, ncol=r_mix_length) diag(Sigma[[i]]) <- B_scale[i] } + #Sample effects from a mixture of MVN distributions or a single MVN distribution B_causal <- matrix(0, nrow=p_causal, ncol=r) if(K>1){ @@ -119,10 +121,15 @@ simulate_mr_mash_data <- function(n, p, p_causal, r, r_causal=list(1:r), interce B[causal_variables, ] <- B_causal ##Simulate X from N_r(0, Gamma) where Gamma is a given covariance matrix across variables - Gamma_offdiag <- X_scale*X_cor - Gamma <- matrix(Gamma_offdiag, nrow=p, ncol=p) - diag(Gamma) <- X_scale - X <- rmvnorm(n=n, mean=rep(0, p), sigma=Gamma) + if(X_cor != 0){ + Gamma_offdiag <- X_scale*X_cor + Gamma <- matrix(Gamma_offdiag, nrow=p, ncol=p) + diag(Gamma) <- X_scale + X <- rmvnorm(n=n, mean=rep(0, p), sigma=Gamma) + } else { + X <- sapply(1:p, sample_norm, n=n, m=0, s2=X_scale) + Gamma <- paste0("I_", p) + } X <- scale_fast2(X, scale=FALSE)$M ##Compute G and its variance @@ -138,7 +145,7 @@ simulate_mr_mash_data <- function(n, p, p_causal, r, r_causal=list(1:r), interce V <- D %*% V_cor_mat %*% D ##Simulate Y from MN(XB, I_n, V) where I_n is an nxn identity matrix and V is the residual covariance - Y <- matrix_normal(G + matrix(intercepts, n, r, byrow=TRUE), diag(n), V) + Y <- matrix_normal_indep_rows(M=(G + matrix(intercepts, n, r, byrow=TRUE)), V=V) ##Compile output causal_responses <- r_causal @@ -165,8 +172,27 @@ simulate_mr_mash_data <- function(n, p, p_causal, r, r_causal=list(1:r), interce } +#' @importFrom stats rnorm +sample_norm <- function(i, n, m, s2){ + x <- rnorm(n=n, mean=m, sd=sqrt(s2)) + + return(x) +} - +#' @importFrom stats rnorm +matrix_normal_indep_rows = function(M, V){ + a <- nrow(M) + b <- ncol(M) + + # Draw Z from MN(O, I, I) + Z <- matrix(rnorm(a*b,0,1), a, b) + + # Cholesky decomposition of V + L2 <- chol(V) + + # Return draw from MN(M,I,V) + return(M + Z %*% L2) +} diff --git a/R/update_order_rss.R b/R/update_order_rss.R new file mode 100644 index 0000000..c2a6bbe --- /dev/null +++ b/R/update_order_rss.R @@ -0,0 +1,44 @@ +###Compute logbf from Bayesian multivariate simple regression with mixture prior +compute_logbf_rss_R <- function(n, XtY, V, Vinv, w0, S0, precomp_quants, standardize, eps, nthreads){ + p <- nrow(XtY) + logbf <- vector("numeric", p) + + for(j in 1:p){ + #Run Bayesian SLR + if(standardize){ + bfit <- bayes_mvr_mix_standardized_X_rss(n, XtY[j, ], w0, S0, precomp_quants$S, precomp_quants$S1, + precomp_quants$SplusS0_chol, precomp_quants$S_chol, eps) + } else { + bfit <- bayes_mvr_mix_centered_X_rss(XtY[j, ], V, w0, S0, precomp_quants$xtx[j], Vinv, + precomp_quants$V_chol, precomp_quants$d, + precomp_quants$QtimesV_chol, eps) + } + + logbf[j] <- bfit$logbf + } + + return(logbf) +} + +###Compute rank of logbf from Bayesian multivariate simple regression with mixture prior +compute_rank_variables_BFmix_rss <- function(n, XtY, V, Vinv, w0, S0, precomp_quants, standardize, version, + decreasing, eps, nthreads){ + if(version=="R"){ + logbfs <- compute_logbf_rss_R(n, XtY, V, Vinv, w0, S0, precomp_quants, standardize, eps) + } else if(version=="Rcpp"){ + logbfs <- compute_logbf_rss_rcpp(n, XtY, V, Vinv, w0, simplify2array_custom(S0), precomp_quants, + standardize, eps , nthreads) + logbfs <- drop(logbfs) + } + + if(decreasing){ + ##Negative sign is needed because rank() by default ranks from smallest to largest + ##while we want from largest to smallest + rank_variables_BFmix <- rank(-logbfs, ties.method="first", na.last="keep") + } else { + rank_variables_BFmix <- rank(logbfs, ties.method="first", na.last="keep") + } + + return(rank_variables_BFmix) +} + diff --git a/README.md b/README.md index e01f8ab..55938fc 100644 --- a/README.md +++ b/README.md @@ -35,8 +35,8 @@ repository useful for your work, please cite: ## License -Copyright (c) 2020-2023, Fabio Morgante, Peter Carbonetto and Matthew -Stephens. +Copyright (c) 2020-2024, Fabio Morgante, Jason Willwerscheid, Gao Wang, Deborah Kunkel, +Peter Carbonetto and Matthew Stephens. All source code and software in this repository are made available under the terms of the [MIT license][mit-license]. diff --git a/docs/404.html b/docs/404.html index f6c7607..0ff09a9 100644 --- a/docs/404.html +++ b/docs/404.html @@ -32,7 +32,7 @@ mr.mash.alpha - 0.2-27 + 0.3.21 @@ -83,12 +83,12 @@

Page not found (404)

diff --git a/docs/LICENSE-text.html b/docs/LICENSE-text.html index 1f2defd..19c2cb2 100644 --- a/docs/LICENSE-text.html +++ b/docs/LICENSE-text.html @@ -17,7 +17,7 @@ mr.mash.alpha - 0.2-27 + 0.3.21 @@ -46,8 +46,9 @@

License

-
YEAR: 2020
-COPYRIGHT HOLDER: Fabio Morgante, 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
 

 
   
@@ -61,11 +62,11 @@ 
License

 
 
       
diff --git a/docs/LICENSE.html b/docs/LICENSE.html
new file mode 100644
index 0000000..320ffa6
--- /dev/null
+++ b/docs/LICENSE.html
@@ -0,0 +1,84 @@
+
+NA • mr.mash.alpha
+    
+
+    

+      

+
+      
+
+      

+  

+    

+      
NA

+    

+
+
+
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 modification, are permitted provided that the following conditions are met:

+
  • Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
    • 
+
  • Redistributions in binary form must reproduce the above copyright 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, 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 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

+
+
+  

+
+  
+
+

+
+
+
+      

+  
Developed by Fabio Morgante, Deborah Kunkel, Peter Carbonetto, Matthew Stephens.

+

+
+

+  
Site built with pkgdown 2.0.7.

+

+
+      

+
+  
+
+
+  
+
+  
+
diff --git a/docs/articles/index.html b/docs/articles/index.html
index ad70b11..526e98d 100644
--- a/docs/articles/index.html
+++ b/docs/articles/index.html
@@ -17,7 +17,7 @@
       
       
         mr.mash.alpha
-        0.2-27
+        0.3.21
       
     
 
@@ -58,11 +58,11 @@ 
All vignettes

 
 
       
diff --git a/docs/articles/mr_mash_intro.html b/docs/articles/mr_mash_intro.html
index 97a5aeb..49e634d 100644
--- a/docs/articles/mr_mash_intro.html
+++ b/docs/articles/mr_mash_intro.html
@@ -33,7 +33,7 @@
       
       
         mr.mash.alpha
-        0.2-27
+        0.3.21
       
     
 
@@ -70,7 +70,7 @@ 
Introduction to mr.mash

                         
Peter
 Carbonetto & Fabio Morgante

             
-            
2023-05-22

+            
2024-04-28

       
       Source: vignettes/mr_mash_intro.Rmd
       
@@ -86,7 +86,7 @@ 
2023-05-22

 and we load the “mr.mash.alpha” package.

 
 library(mr.mash.alpha)
-set.seed(123)

+set.seed(12)

We illustrate the application of mr.mash to a data set simulated from a multivariate, multiple linear regression with 5 responses in which the coefficients are the same for all responses. In the target application @@ -206,29 +206,29 @@

Fit a mr.mash model to the data 
 head(sort(fit$w0,decreasing = TRUE),n = 10)
-#             null   shared1_grid13   shared1_grid12   shared1_grid11 
-#     9.925505e-01     2.156501e-03     1.279544e-03     2.432624e-04 
-#   shared1_grid14   shared1_grid10 singleton1_grid1 singleton3_grid1 
-#     1.673357e-04     1.171963e-04     6.237657e-05     6.236753e-05 
-# singleton2_grid1 singleton1_grid2 
-#     6.211754e-05     6.208252e-05
+# null shared1_grid12 singleton1_grid10 singleton5_grid8 +# 9.822578e-01 5.431954e-03 4.477378e-03 3.180946e-03 +# singleton5_grid7 singleton5_grid6 singleton1_grid11 singleton5_grid5 +# 1.531493e-03 3.655315e-04 1.818372e-04 1.496145e-04 +# singleton5_grid4 singleton5_grid3 +# 9.529390e-05 7.656325e-05

Also, reassuringly, the estimated residual variance-covariance matrix is close to the matrix used to simulate the data:

 dat$V
-#           [,1]      [,2]      [,3]      [,4]      [,5]
-# [1,] 3.5145374 0.8786343 0.8786343 0.8786343 0.8786343
-# [2,] 0.8786343 3.5145373 0.8786343 0.8786343 0.8786343
-# [3,] 0.8786343 0.8786343 3.5145373 0.8786343 0.8786343
-# [4,] 0.8786343 0.8786343 0.8786343 3.5145373 0.8786343
-# [5,] 0.8786343 0.8786343 0.8786343 0.8786343 3.5145373
+#          [,1]     [,2]     [,3]     [,4]     [,5]
+# [1,] 6.622079 1.655520 1.655520 1.655520 1.655520
+# [2,] 1.655520 6.622079 1.655520 1.655520 1.655520
+# [3,] 1.655520 1.655520 6.622079 1.655520 1.655520
+# [4,] 1.655520 1.655520 1.655520 6.622079 1.655520
+# [5,] 1.655520 1.655520 1.655520 1.655520 6.622079
 fit$V
-#          [,1]      [,2]      [,3]      [,4]      [,5]
-# [1,] 3.461066 0.9454000 0.9492950 1.4169111 1.0618637
-# [2,] 0.945400 3.7718288 1.3080240 0.7038973 0.1283624
-# [3,] 0.949295 1.3080240 2.7266827 0.7350486 1.0846300
-# [4,] 1.416911 0.7038973 0.7350486 2.9549883 0.8032564
-# [5,] 1.061864 0.1283624 1.0846300 0.8032564 3.1519688
+# [,1] [,2] [,3] [,4] [,5] +# [1,] 5.7356423 1.585484 1.3494127 0.8463714 0.9642602 +# [2,] 1.5854844 6.578508 2.4138518 1.5273860 1.6646585 +# [3,] 1.3494127 2.413852 6.0258742 0.4755238 1.5601590 +# [4,] 0.8463714 1.527386 0.4755238 4.7464139 1.3982834 +# [5,] 0.9642602 1.664659 1.5601590 1.3982834 5.0244255

Use the fitted mr.mash model to make predictions @@ -268,12 +268,12 @@

Use the fitted mr.mash

-

Developed by Fabio Morgante, Peter Carbonetto, Matthew Stephens.

+

Developed by Fabio Morgante, Deborah Kunkel, Peter Carbonetto, Matthew Stephens.

-

Site built with pkgdown 2.0.6.

+

Site built with pkgdown 2.0.7.

diff --git a/docs/articles/mr_mash_intro_files/figure-html/plot-coefs-1.png b/docs/articles/mr_mash_intro_files/figure-html/plot-coefs-1.png index 8ec5aca..0972596 100644 Binary files a/docs/articles/mr_mash_intro_files/figure-html/plot-coefs-1.png and b/docs/articles/mr_mash_intro_files/figure-html/plot-coefs-1.png differ diff --git a/docs/articles/mr_mash_intro_files/figure-html/plot-pred-test-1.png b/docs/articles/mr_mash_intro_files/figure-html/plot-pred-test-1.png index 92909df..548333e 100644 Binary files a/docs/articles/mr_mash_intro_files/figure-html/plot-pred-test-1.png and b/docs/articles/mr_mash_intro_files/figure-html/plot-pred-test-1.png differ diff --git a/docs/authors.html b/docs/authors.html index effe8ee..c501766 100644 --- a/docs/authors.html +++ b/docs/authors.html @@ -17,7 +17,7 @@ mr.mash.alpha - 0.2-27 + 0.3.21

@@ -52,6 +52,18 @@

Authors

Fabio Morgante. Maintainer, author.

+
  • +

    Jason Willwerscheid. Contributor. +

    +
    • +
  • +

    Gao Wang. Contributor. +

    +
    • +
  • +

    Deborah Kunkel. Author. +

    +

  • Peter Carbonetto. Author.

    @@ -69,15 +81,15 @@

    Citation

    -

    Morgante F, Carbonetto P, Stephens M (2023). +

    Morgante F, Kunkel D, Carbonetto P, Stephens M (2024). mr.mash.alpha: Multiple Regression with Multivariate Adaptive Shrinkage. -R package version 0.2-27, https://github.com/stephenslab/mr.mash.alpha. +R package version 0.3.21, https://github.com/stephenslab/mr.mash.alpha. 

    @Manual{,
   title = {mr.mash.alpha: Multiple Regression with Multivariate Adaptive Shrinkage},
-  author = {Fabio Morgante and Peter Carbonetto and Matthew Stephens},
-  year = {2023},
-  note = {R package version 0.2-27},
+  author = {Fabio Morgante and Deborah Kunkel and Peter Carbonetto and Matthew Stephens},
+  year = {2024},
+  note = {R package version 0.3.21},
   url = {https://github.com/stephenslab/mr.mash.alpha},
 }
    @@ -88,11 +100,11 @@

    Citation

    diff --git a/docs/index.html b/docs/index.html index a9640f8..3c7be0d 100644 --- a/docs/index.html +++ b/docs/index.html @@ -34,7 +34,7 @@ mr.mash.alpha - 0.2-27 + 0.3.21 @@ -69,7 +69,7 @@

    mr.mash.alpha

    - +

    R-CMD-check CircleCI build status

    R package implementing methods for multivariate multiple regression with adaptive shrinkage priors (mr.mash).

    Quick Start @@ -84,12 +84,12 @@

    Quick Start

    Citing this work

    -

    If you find the mr.mash.alpha package or any of the source code in this repository useful for your work, please cite: > Morgante, F., Carbonetto, P., Wang, G., Zou, Y., Sarkar, A. & > Stephens, M. (2022). A flexible empirical Bayes approach to > multivariate multiple regression, and its improved accuracy > in predicting multi-tissue gene expression from genotypes. > bioRxiv https://doi.org/10.1101/2022.11.22.517471

    +

    If you find the mr.mash.alpha package or any of the source code in this repository useful for your work, please cite: > Morgante, F., Carbonetto, P., Wang, G., Zou, Y., Sarkar, A. & > Stephens, M. (2023). A flexible empirical Bayes approach to > multivariate multiple regression, and its improved accuracy > in predicting multi-tissue gene expression from genotypes. > PLoS Genetics 19(7): e1010539. https://doi.org/10.1371/journal.pgen.1010539

    License

    -

    Copyright (c) 2020-2022, Fabio Morgante, Peter Carbonetto and Matthew Stephens.

    +

    Copyright (c) 2020-2024, Fabio Morgante, Jason Willwerscheid, Gao Wang, Deborah Kunkel, Peter Carbonetto and Matthew Stephens.

    All source code and software in this repository are made available under the terms of the MIT license.

    @@ -106,9 +106,8 @@

    Links

    License

    @@ -124,18 +123,14 @@

    Citation

    Developers

    -
    -

    Dev status

    - -
    + @@ -143,12 +138,12 @@

    Dev status

    diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml index 90f63a0..26aa5b2 100644 --- a/docs/pkgdown.yml +++ b/docs/pkgdown.yml @@ -1,7 +1,7 @@ -pandoc: 3.1.2 -pkgdown: 2.0.6 +pandoc: 3.1.1 +pkgdown: 2.0.7 pkgdown_sha: ~ articles: mr_mash_intro: mr_mash_intro.html -last_built: 2023-05-22T14:18Z +last_built: 2024-04-28T21:37Z diff --git a/docs/reference/Rplot002.png b/docs/reference/Rplot002.png index 67e94ad..cb645bc 100644 Binary files a/docs/reference/Rplot002.png and b/docs/reference/Rplot002.png differ diff --git a/docs/reference/autoselect.mixsd.html b/docs/reference/autoselect.mixsd.html index 08fedac..cce3122 100644 --- a/docs/reference/autoselect.mixsd.html +++ b/docs/reference/autoselect.mixsd.html @@ -18,7 +18,7 @@ mr.mash.alpha - 0.2-27 + 0.3.21 @@ -84,11 +84,11 @@

    Value

    diff --git a/docs/reference/coef.mr.mash.html b/docs/reference/coef.mr.mash.html index 34f438e..f0fcf7e 100644 --- a/docs/reference/coef.mr.mash.html +++ b/docs/reference/coef.mr.mash.html @@ -17,7 +17,7 @@ mr.mash.alpha - 0.2-27 + 0.3.21 @@ -82,11 +82,11 @@

    Value

    diff --git a/docs/reference/coef.mr.mash.rss.html b/docs/reference/coef.mr.mash.rss.html new file mode 100644 index 0000000..486646c --- /dev/null +++ b/docs/reference/coef.mr.mash.rss.html @@ -0,0 +1,101 @@ + +Extract coefficients from mr.mash.rss fit. — coef.mr.mash.rss • mr.mash.alpha + + +
    +
    + + + +
    +
    +
    +

    Extract coefficients from mr.mash.rss fit.

    + Source: R/mr_mash_rss_predict.R + +
    + +
    +

    Extract coefficients from mr.mash.rss fit.

    +
    + +
    + 
    # S3 method for mr.mash.rss
+coef(object, ...)
    +
    + +
    +

    Arguments

    +
    object
    +

    a mr.mash fit.

    + + +
    ...
    +

    Other arguments (not used).

    + +
    +
    +

    Value

    + + +

    p x r or (p+1) x r matrix of coefficients, + depending on whether an intercept was computed.

    +
    + +
    + +
    + + +
    +

    Developed by Fabio Morgante, Deborah Kunkel, Peter Carbonetto, Matthew Stephens.

    +
    + +
    +

    Site built with pkgdown 2.0.7.

    +
    + +
    + + + + + + + + diff --git a/docs/reference/compute_canonical_covs.html b/docs/reference/compute_canonical_covs.html index 2b530ec..a1205e6 100644 --- a/docs/reference/compute_canonical_covs.html +++ b/docs/reference/compute_canonical_covs.html @@ -17,7 +17,7 @@ mr.mash.alpha - 0.2-27 + 0.3.21 @@ -91,11 +91,11 @@

    Value

    diff --git a/docs/reference/compute_data_driven_covs.html b/docs/reference/compute_data_driven_covs.html index ef495a3..b78a441 100644 --- a/docs/reference/compute_data_driven_covs.html +++ b/docs/reference/compute_data_driven_covs.html @@ -18,7 +18,7 @@ mr.mash.alpha - 0.2-27 + 0.3.21 @@ -110,11 +110,11 @@

    Value

    diff --git a/docs/reference/compute_univariate_sumstats.html b/docs/reference/compute_univariate_sumstats.html index adea763..77de0e8 100644 --- a/docs/reference/compute_univariate_sumstats.html +++ b/docs/reference/compute_univariate_sumstats.html @@ -18,7 +18,7 @@ mr.mash.alpha - 0.2-27 + 0.3.21 @@ -110,11 +110,11 @@

    Value

    diff --git a/docs/reference/expand_covs.html b/docs/reference/expand_covs.html index a580d60..9c94ac8 100644 --- a/docs/reference/expand_covs.html +++ b/docs/reference/expand_covs.html @@ -17,7 +17,7 @@ mr.mash.alpha - 0.2-27 + 0.3.21 @@ -85,11 +85,11 @@

    Value

    diff --git a/docs/reference/index.html b/docs/reference/index.html index e2ed427..eaec2cb 100644 --- a/docs/reference/index.html +++ b/docs/reference/index.html @@ -17,7 +17,7 @@ mr.mash.alpha - 0.2-27 + 0.3.21 @@ -58,6 +58,10 @@

    All functions coef(<mr.mash>)

    Extract coefficients from mr.mash fit.

    + +

    coef(<mr.mash.rss>)

    + +

    Extract coefficients from mr.mash.rss fit.

    compute_canonical_covs()

    @@ -78,10 +82,19 @@

    All functions mr.mash()

    Multiple Regression with Multivariate Adaptive Shrinkage.

    + +

    mr.mash.rss()

    + +

    Multiple Regression with Multivariate Adaptive Shrinkage + from summary data.

    predict(<mr.mash>)

    Predict future observations from mr.mash fit.

    + +

    predict(<mr.mash.rss>)

    + +

    Predict future observations from mr.mash.rss fit.

    simulate_mr_mash_data()

    @@ -95,11 +108,11 @@

    All functions
    -

    Developed by Fabio Morgante, Peter Carbonetto, Matthew Stephens.

    +

    Developed by Fabio Morgante, Deborah Kunkel, Peter Carbonetto, Matthew Stephens.

    -

    Site built with pkgdown 2.0.6.

    +

    Site built with pkgdown 2.0.7.

    diff --git a/docs/reference/mr.mash-1.png b/docs/reference/mr.mash-1.png index 3c5fea3..342e032 100644 Binary files a/docs/reference/mr.mash-1.png and b/docs/reference/mr.mash-1.png differ diff --git a/docs/reference/mr.mash-2.png b/docs/reference/mr.mash-2.png index 072ccdb..28cfe1a 100644 Binary files a/docs/reference/mr.mash-2.png and b/docs/reference/mr.mash-2.png differ diff --git a/docs/reference/mr.mash.html b/docs/reference/mr.mash.html index 5ab9963..67307a5 100644 --- a/docs/reference/mr.mash.html +++ b/docs/reference/mr.mash.html @@ -18,7 +18,7 @@ mr.mash.alpha - 0.2-27 + 0.3.21 @@ -75,7 +75,7 @@

    Multiple Regression with Multivariate Adaptive Shrinkage.

    update_V_method = c("full", "diagonal"), version = c("Rcpp", "R"), e = 1e-08, - ca_update_order = c("consecutive", "decreasing_logBF", "increasing_logBF"), + ca_update_order = c("consecutive", "decreasing_logBF", "increasing_logBF", "random"), nthreads = as.integer(NA) ) @@ -181,7 +181,7 @@

    Arguments

    ca_update_order

    The order with which coordinates are updated. So far, "consecutive", "decreasing_logBF", -"increasing_logBF" are supported.

    +"increasing_logBF", "random" are supported.

    nthreads
    @@ -296,189 +296,189 @@

    Examples

    ###Fit mr.mash fit <- mr.mash(Xtrain, Ytrain, S0, update_V=TRUE) #> Processing the inputs... Done! -#> Fitting the optimization algorithm using 20 RcppParallel threads... +#> Fitting the optimization algorithm using 8 RcppParallel threads... #> iter mu1_max.diff ELBO_diff ELBO -#> 1 1.71e+00 Inf -1.17642528220990225236e+04 -#> 2 3.16e-01 3.86e+02 -1.13785636762104713853e+04 -#> 3 7.32e-02 2.38e+01 -1.13547399116407923430e+04 -#> 4 2.64e-02 2.75e+00 -1.13519874874270280998e+04 -#> 5 1.63e-02 1.28e+00 -1.13507113822810115380e+04 -#> 6 1.23e-02 7.76e-01 -1.13499353149131493410e+04 -#> 7 9.55e-03 5.23e-01 -1.13494120998063317529e+04 -#> 8 7.63e-03 3.77e-01 -1.13490351417299498280e+04 -#> 9 6.23e-03 2.85e-01 -1.13487496678169645747e+04 -#> 10 5.18e-03 2.25e-01 -1.13485247734950517042e+04 -#> 11 4.38e-03 1.83e-01 -1.13483418476869410370e+04 -#> 12 3.76e-03 1.53e-01 -1.13481890954040336510e+04 -#> 13 3.27e-03 1.30e-01 -1.13480587207075168408e+04 -#> 14 2.87e-03 1.13e-01 -1.13479453863624312362e+04 -#> 15 2.54e-03 1.00e-01 -1.13478453286182575539e+04 -#> 16 2.27e-03 8.95e-02 -1.13477558262726124667e+04 -#> 17 2.05e-03 8.10e-02 -1.13476748699693507660e+04 -#> 18 1.86e-03 7.39e-02 -1.13476009491374661593e+04 -#> 19 1.70e-03 6.80e-02 -1.13475329105017372058e+04 -#> 20 1.56e-03 6.30e-02 -1.13474698615490688098e+04 -#> 21 1.44e-03 5.88e-02 -1.13474111030748608755e+04 -#> 22 1.34e-03 5.50e-02 -1.13473560810572562332e+04 -#> 23 1.24e-03 5.17e-02 -1.13473043517017777049e+04 -#> 24 1.16e-03 4.88e-02 -1.13472555556662555318e+04 -#> 25 1.15e-03 4.62e-02 -1.13472093988159122091e+04 -#> 26 1.14e-03 4.38e-02 -1.13471656377071776660e+04 -#> 27 1.13e-03 4.16e-02 -1.13471240685484717687e+04 -#> 28 1.11e-03 3.95e-02 -1.13470845187506092770e+04 -#> 29 1.09e-03 3.77e-02 -1.13470468404263883713e+04 -#> 30 1.06e-03 3.59e-02 -1.13470109053693431633e+04 -#> 31 1.03e-03 3.43e-02 -1.13469766011625943065e+04 -#> 32 9.99e-04 3.28e-02 -1.13469438281548773375e+04 -#> 33 9.66e-04 3.13e-02 -1.13469124971045639541e+04 -#> 34 9.32e-04 3.00e-02 -1.13468825273393049429e+04 -#> 35 8.97e-04 2.87e-02 -1.13468538453142664366e+04 -#> 36 8.62e-04 2.75e-02 -1.13468263834788267559e+04 -#> 37 8.45e-04 2.63e-02 -1.13468000793820228864e+04 -#> 38 8.27e-04 2.52e-02 -1.13467748749628535734e+04 -#> 39 8.08e-04 2.42e-02 -1.13467507159838860389e+04 -#> 40 7.88e-04 2.32e-02 -1.13467275515757955873e+04 -#> 41 7.69e-04 2.22e-02 -1.13467053338681453170e+04 -#> 42 7.48e-04 2.13e-02 -1.13466840176871337462e+04 -#> 43 7.28e-04 2.05e-02 -1.13466635603056383843e+04 -#> 44 7.08e-04 1.96e-02 -1.13466439212341210805e+04 -#> 45 6.87e-04 1.89e-02 -1.13466250620438004262e+04 -#> 46 6.67e-04 1.81e-02 -1.13466069462152881897e+04 -#> 47 6.46e-04 1.74e-02 -1.13465895390076802869e+04 -#> 48 6.26e-04 1.67e-02 -1.13465728073440150183e+04 -#> 49 6.06e-04 1.61e-02 -1.13465567197101463535e+04 -#> 50 5.86e-04 1.55e-02 -1.13465412460646311956e+04 -#> 51 5.67e-04 1.49e-02 -1.13465263577577770775e+04 -#> 52 5.48e-04 1.43e-02 -1.13465120274583514401e+04 -#> 53 5.29e-04 1.38e-02 -1.13464982290868829296e+04 -#> 54 5.10e-04 1.33e-02 -1.13464849377545342577e+04 -#> 55 4.92e-04 1.28e-02 -1.13464721297067571868e+04 -#> 56 4.74e-04 1.23e-02 -1.13464597822712312336e+04 -#> 57 4.57e-04 1.19e-02 -1.13464478738093584980e+04 -#> 58 4.40e-04 1.15e-02 -1.13464363836710490432e+04 -#> 59 4.23e-04 1.11e-02 -1.13464252921522620454e+04 -#> 60 4.07e-04 1.07e-02 -1.13464145804550535104e+04 -#> 61 3.91e-04 1.03e-02 -1.13464042306498577091e+04 -#> 62 3.76e-04 1.00e-02 -1.13463942256396585435e+04 -#> 63 3.61e-04 9.68e-03 -1.13463845491259435221e+04 -#> 64 3.47e-04 9.36e-03 -1.13463751855761784100e+04 -#> 65 3.32e-04 9.07e-03 -1.13463661201927152433e+04 -#> 66 3.19e-04 8.78e-03 -1.13463573388829008763e+04 -#> 67 3.06e-04 8.51e-03 -1.13463488282303751475e+04 -#> 68 2.93e-04 8.25e-03 -1.13463405754673603951e+04 -#> 69 2.80e-04 8.01e-03 -1.13463325684479768825e+04 -#> 70 2.68e-04 7.77e-03 -1.13463247956223931396e+04 -#> 71 2.57e-04 7.55e-03 -1.13463172460118184972e+04 -#> 72 2.45e-04 7.34e-03 -1.13463099091843414499e+04 -#> 73 2.34e-04 7.13e-03 -1.13463027752314555983e+04 -#> 74 2.24e-04 6.94e-03 -1.13462958347453786700e+04 -#> 75 2.14e-04 6.76e-03 -1.13462890787970027304e+04 -#> 76 2.04e-04 6.58e-03 -1.13462824989146156440e+04 -#> 77 1.95e-04 6.41e-03 -1.13462760870632246224e+04 -#> 78 1.86e-04 6.25e-03 -1.13462698356246055482e+04 -#> 79 1.80e-04 6.10e-03 -1.13462637373779434711e+04 -#> 80 1.74e-04 5.95e-03 -1.13462577854811770521e+04 -#> 81 1.69e-04 5.81e-03 -1.13462519734529723792e+04 -#> 82 1.65e-04 5.68e-03 -1.13462462951552915911e+04 -#> 83 1.60e-04 5.55e-03 -1.13462407447766618134e+04 -#> 84 1.55e-04 5.43e-03 -1.13462353168159825145e+04 -#> 85 1.51e-04 5.31e-03 -1.13462300060670459061e+04 -#> 86 1.47e-04 5.20e-03 -1.13462248076035757549e+04 -#> 87 1.43e-04 5.09e-03 -1.13462197167649374023e+04 -#> 88 1.38e-04 4.99e-03 -1.13462147291423807474e+04 -#> 89 1.35e-04 4.89e-03 -1.13462098405659035052e+04 -#> 90 1.31e-04 4.79e-03 -1.13462050470916674385e+04 -#> 91 1.27e-04 4.70e-03 -1.13462003449899839325e+04 -#> 92 1.23e-04 4.61e-03 -1.13461957307338107057e+04 -#> 93 1.20e-04 4.53e-03 -1.13461912009878160461e+04 -#> 94 1.17e-04 4.45e-03 -1.13461867525979341735e+04 -#> 95 1.13e-04 4.37e-03 -1.13461823825814262818e+04 -#> 96 1.10e-04 4.29e-03 -1.13461780881174272508e+04 -#> 97 1.13e-04 4.22e-03 -1.13461738665379689337e+04 -#> 98 1.17e-04 4.15e-03 -1.13461697153194254497e+04 -#> 99 1.20e-04 4.08e-03 -1.13461656320744168625e+04 -#> 100 1.23e-04 4.02e-03 -1.13461616145441475965e+04 -#> 101 1.25e-04 3.95e-03 -1.13461576605910940998e+04 -#> 102 1.28e-04 3.89e-03 -1.13461537681921363401e+04 -#> 103 1.30e-04 3.83e-03 -1.13461499354320203565e+04 -#> 104 1.32e-04 3.77e-03 -1.13461461604971591441e+04 -#> 105 1.34e-04 3.72e-03 -1.13461424416698482673e+04 -#> 106 1.36e-04 3.66e-03 -1.13461387773226979334e+04 -#> 107 1.38e-04 3.61e-03 -1.13461351659134415968e+04 -#> 108 1.40e-04 3.56e-03 -1.13461316059800101357e+04 -#> 109 1.41e-04 3.51e-03 -1.13461280961358843342e+04 -#> 110 1.43e-04 3.46e-03 -1.13461246350657238509e+04 -#> 111 1.44e-04 3.41e-03 -1.13461212215212108276e+04 -#> 112 1.45e-04 3.37e-03 -1.13461178543171663478e+04 -#> 113 1.46e-04 3.32e-03 -1.13461145323278578871e+04 -#> 114 1.47e-04 3.28e-03 -1.13461112544835632434e+04 -#> 115 1.48e-04 3.23e-03 -1.13461080197672672512e+04 -#> 116 1.48e-04 3.19e-03 -1.13461048272115967848e+04 -#> 117 1.49e-04 3.15e-03 -1.13461016758959467552e+04 -#> 118 1.50e-04 3.11e-03 -1.13460985649437061511e+04 -#> 119 1.50e-04 3.07e-03 -1.13460954935197078157e+04 -#> 120 1.50e-04 3.03e-03 -1.13460924608278055530e+04 -#> 121 1.51e-04 2.99e-03 -1.13460894661085876578e+04 -#> 122 1.51e-04 2.96e-03 -1.13460865086372141377e+04 -#> 123 1.51e-04 2.92e-03 -1.13460835877214140055e+04 -#> 124 1.51e-04 2.89e-03 -1.13460807026995498745e+04 -#> 125 1.51e-04 2.85e-03 -1.13460778529388408060e+04 -#> 126 1.51e-04 2.82e-03 -1.13460750378336724680e+04 -#> 127 1.50e-04 2.78e-03 -1.13460722568039691396e+04 -#> 128 1.50e-04 2.75e-03 -1.13460695092937439767e+04 -#> 129 1.50e-04 2.71e-03 -1.13460667947696310875e+04 -#> 130 1.49e-04 2.68e-03 -1.13460641127195758600e+04 -#> 131 1.49e-04 2.65e-03 -1.13460614626515998680e+04 -#> 132 1.49e-04 2.62e-03 -1.13460588440925548639e+04 -#> 133 1.48e-04 2.59e-03 -1.13460562565870823164e+04 -#> 134 1.47e-04 2.56e-03 -1.13460536996965038270e+04 -#> 135 1.47e-04 2.53e-03 -1.13460511729978443327e+04 -#> 136 1.46e-04 2.50e-03 -1.13460486760829226114e+04 -#> 137 1.45e-04 2.47e-03 -1.13460462085574417870e+04 -#> 138 1.45e-04 2.44e-03 -1.13460437700401580514e+04 -#> 139 1.44e-04 2.41e-03 -1.13460413601621221460e+04 -#> 140 1.43e-04 2.38e-03 -1.13460389785659099289e+04 -#> 141 1.42e-04 2.35e-03 -1.13460366249049602629e+04 -#> 142 1.41e-04 2.33e-03 -1.13460342988429001707e+04 -#> 143 1.40e-04 2.30e-03 -1.13460320000529063691e+04 -#> 144 1.39e-04 2.27e-03 -1.13460297282171777624e+04 -#> 145 1.38e-04 2.25e-03 -1.13460274830263297190e+04 -#> 146 1.37e-04 2.22e-03 -1.13460252641788993060e+04 -#> 147 1.36e-04 2.19e-03 -1.13460230713808687142e+04 -#> 148 1.35e-04 2.17e-03 -1.13460209043451959587e+04 -#> 149 1.34e-04 2.14e-03 -1.13460187627913674078e+04 -#> 150 1.33e-04 2.12e-03 -1.13460166464450176136e+04 -#> 151 1.32e-04 2.09e-03 -1.13460145550375273160e+04 -#> 152 1.31e-04 2.07e-03 -1.13460124883056687395e+04 -#> 153 1.30e-04 2.04e-03 -1.13460104459912581660e+04 -#> 154 1.29e-04 2.02e-03 -1.13460084278408448881e+04 -#> 155 1.28e-04 1.99e-03 -1.13460064336054219893e+04 -#> 156 1.27e-04 1.97e-03 -1.13460044630401098402e+04 -#> 157 1.25e-04 1.95e-03 -1.13460025159039360005e+04 -#> 158 1.24e-04 1.92e-03 -1.13460005919595387240e+04 -#> 159 1.23e-04 1.90e-03 -1.13459986909729868785e+04 -#> 160 1.22e-04 1.88e-03 -1.13459968127135070972e+04 -#> 161 1.21e-04 1.86e-03 -1.13459949569533218892e+04 -#> 162 1.19e-04 1.83e-03 -1.13459931234674386360e+04 -#> 163 1.18e-04 1.81e-03 -1.13459913120334767882e+04 -#> 164 1.17e-04 1.79e-03 -1.13459895224315005180e+04 -#> 165 1.16e-04 1.77e-03 -1.13459877544438604673e+04 -#> 166 1.14e-04 1.75e-03 -1.13459860078550482285e+04 -#> 167 1.13e-04 1.73e-03 -1.13459842824515726534e+04 -#> 168 1.12e-04 1.70e-03 -1.13459825780218197906e+04 -#> 169 1.11e-04 1.68e-03 -1.13459808943559455656e+04 -#> 170 1.09e-04 1.66e-03 -1.13459792312457830121e+04 -#> 171 1.08e-04 1.64e-03 -1.13459775884847022098e+04 -#> 172 1.07e-04 1.62e-03 -1.13459759658675611718e+04 -#> 173 1.06e-04 1.60e-03 -1.13459743631906167138e+04 -#> 174 1.05e-04 1.58e-03 -1.13459727802514189534e+04 -#> 175 1.03e-04 1.56e-03 -1.13459712168487694726e+04 -#> 176 1.02e-04 1.54e-03 -1.13459696727826503775e+04 -#> 177 1.01e-04 1.52e-03 -1.13459681478541679098e+04 -#> 178 9.98e-05 1.51e-03 -1.13459666418654851441e+04 +#> 1 1.94e+00 Inf -1.17846672186051764584e+04 +#> 2 4.42e-01 4.05e+02 -1.13796024155407103535e+04 +#> 3 7.87e-02 2.57e+01 -1.13539336660288772691e+04 +#> 4 2.58e-02 2.82e+00 -1.13511145940660353517e+04 +#> 5 1.32e-02 1.31e+00 -1.13498072323880332988e+04 +#> 6 8.97e-03 7.79e-01 -1.13490277541913474124e+04 +#> 7 7.04e-03 5.18e-01 -1.13485098888134161825e+04 +#> 8 5.70e-03 3.70e-01 -1.13481397456728846009e+04 +#> 9 4.72e-03 2.79e-01 -1.13478605360439596552e+04 +#> 10 3.99e-03 2.19e-01 -1.13476411019760180352e+04 +#> 11 3.42e-03 1.78e-01 -1.13474630696575877664e+04 +#> 12 2.97e-03 1.48e-01 -1.13473149527450605092e+04 +#> 13 2.61e-03 1.26e-01 -1.13471892152355630969e+04 +#> 14 2.32e-03 1.09e-01 -1.13470807097645356407e+04 +#> 15 2.07e-03 9.49e-02 -1.13469857995247657527e+04 +#> 16 1.87e-03 8.40e-02 -1.13469018398800035357e+04 +#> 17 1.69e-03 7.50e-02 -1.13468268590668303659e+04 +#> 18 1.57e-03 6.75e-02 -1.13467593542255581269e+04 +#> 19 1.48e-03 6.12e-02 -1.13466981570112020563e+04 +#> 20 1.40e-03 5.58e-02 -1.13466423427090321638e+04 +#> 21 1.36e-03 5.12e-02 -1.13465911674013368611e+04 +#> 22 1.36e-03 4.71e-02 -1.13465440236996964813e+04 +#> 23 1.35e-03 4.36e-02 -1.13465004090367383469e+04 +#> 24 1.32e-03 4.05e-02 -1.13464599026111500280e+04 +#> 25 1.30e-03 3.78e-02 -1.13464221483868714131e+04 +#> 26 1.26e-03 3.53e-02 -1.13463868423840212927e+04 +#> 27 1.22e-03 3.31e-02 -1.13463537230474576063e+04 +#> 28 1.18e-03 3.12e-02 -1.13463225638458152389e+04 +#> 29 1.14e-03 2.94e-02 -1.13462931675035342778e+04 +#> 30 1.09e-03 2.78e-02 -1.13462653614406572160e+04 +#> 31 1.05e-03 2.64e-02 -1.13462389941154578992e+04 +#> 32 1.00e-03 2.51e-02 -1.13462139320497208246e+04 +#> 33 9.57e-04 2.39e-02 -1.13461900573764960427e+04 +#> 34 9.13e-04 2.28e-02 -1.13461672657931958383e+04 +#> 35 8.69e-04 2.18e-02 -1.13461454648335638922e+04 +#> 36 8.27e-04 2.09e-02 -1.13461245723942793120e+04 +#> 37 7.87e-04 2.01e-02 -1.13461045154679814004e+04 +#> 38 7.48e-04 1.93e-02 -1.13460852290460861695e+04 +#> 39 7.10e-04 1.86e-02 -1.13460666551632675692e+04 +#> 40 6.75e-04 1.79e-02 -1.13460487420616555028e+04 +#> 41 6.41e-04 1.73e-02 -1.13460314434572956088e+04 +#> 42 6.09e-04 1.67e-02 -1.13460147178949864610e+04 +#> 43 5.78e-04 1.62e-02 -1.13459985281798235519e+04 +#> 44 5.50e-04 1.57e-02 -1.13459828408760185994e+04 +#> 45 5.23e-04 1.52e-02 -1.13459676258648305520e+04 +#> 46 4.97e-04 1.48e-02 -1.13459528559546215547e+04 +#> 47 4.73e-04 1.43e-02 -1.13459385065370843222e+04 +#> 48 4.51e-04 1.40e-02 -1.13459245552843021869e+04 +#> 49 4.30e-04 1.36e-02 -1.13459109818820033979e+04 +#> 50 4.11e-04 1.32e-02 -1.13458977677949151257e+04 +#> 51 3.98e-04 1.29e-02 -1.13458848960604354943e+04 +#> 52 3.85e-04 1.25e-02 -1.13458723511075259012e+04 +#> 53 3.72e-04 1.22e-02 -1.13458601185977404384e+04 +#> 54 3.61e-04 1.19e-02 -1.13458481852858749335e+04 +#> 55 3.49e-04 1.16e-02 -1.13458365388979145791e+04 +#> 56 3.38e-04 1.14e-02 -1.13458251680241410213e+04 +#> 57 3.28e-04 1.11e-02 -1.13458140620256563125e+04 +#> 58 3.18e-04 1.09e-02 -1.13458032109525811393e+04 +#> 59 3.08e-04 1.06e-02 -1.13457926054725921858e+04 +#> 60 2.99e-04 1.04e-02 -1.13457822368083907350e+04 +#> 61 2.90e-04 1.01e-02 -1.13457720966830856923e+04 +#> 62 2.82e-04 9.92e-03 -1.13457621772723778122e+04 +#> 63 2.73e-04 9.71e-03 -1.13457524711627811485e+04 +#> 64 2.66e-04 9.50e-03 -1.13457429713149813324e+04 +#> 65 2.58e-04 9.30e-03 -1.13457336710317831603e+04 +#> 66 2.51e-04 9.11e-03 -1.13457245639299435425e+04 +#> 67 2.44e-04 8.92e-03 -1.13457156439154150576e+04 +#> 68 2.37e-04 8.74e-03 -1.13457069051615471835e+04 +#> 69 2.30e-04 8.56e-03 -1.13456983420898977784e+04 +#> 70 2.24e-04 8.39e-03 -1.13456899493531527696e+04 +#> 71 2.17e-04 8.23e-03 -1.13456817218200667412e+04 +#> 72 2.11e-04 8.07e-03 -1.13456736545620060497e+04 +#> 73 2.06e-04 7.91e-03 -1.13456657428409325803e+04 +#> 74 2.00e-04 7.76e-03 -1.13456579820986335108e+04 +#> 75 1.95e-04 7.61e-03 -1.13456503679470260977e+04 +#> 76 1.89e-04 7.47e-03 -1.13456428961593574058e+04 +#> 77 1.84e-04 7.33e-03 -1.13456355626622407726e+04 +#> 78 1.79e-04 7.20e-03 -1.13456283635283780313e+04 +#> 79 1.74e-04 7.07e-03 -1.13456212949698547163e+04 +#> 80 1.70e-04 6.94e-03 -1.13456143533320100687e+04 +#> 81 1.66e-04 6.82e-03 -1.13456075350877035817e+04 +#> 82 1.63e-04 6.70e-03 -1.13456008368320217414e+04 +#> 83 1.60e-04 6.58e-03 -1.13455942552773249190e+04 +#> 84 1.56e-04 6.47e-03 -1.13455877872485798434e+04 +#> 85 1.57e-04 6.36e-03 -1.13455814296789831133e+04 +#> 86 1.59e-04 6.25e-03 -1.13455751796058266336e+04 +#> 87 1.60e-04 6.15e-03 -1.13455690341665576852e+04 +#> 88 1.61e-04 6.04e-03 -1.13455629905950609100e+04 +#> 89 1.62e-04 5.94e-03 -1.13455570462180694449e+04 +#> 90 1.62e-04 5.85e-03 -1.13455511984517888777e+04 +#> 91 1.63e-04 5.75e-03 -1.13455454447986412561e+04 +#> 92 1.64e-04 5.66e-03 -1.13455397828441164165e+04 +#> 93 1.64e-04 5.57e-03 -1.13455342102538015752e+04 +#> 94 1.65e-04 5.49e-03 -1.13455287247704745823e+04 +#> 95 1.65e-04 5.40e-03 -1.13455233242113645247e+04 +#> 96 1.65e-04 5.32e-03 -1.13455180064654359740e+04 +#> 97 1.65e-04 5.24e-03 -1.13455127694908624107e+04 +#> 98 1.66e-04 5.16e-03 -1.13455076113125214761e+04 +#> 99 1.66e-04 5.08e-03 -1.13455025300196248281e+04 +#> 100 1.66e-04 5.01e-03 -1.13454975237633934739e+04 +#> 101 1.66e-04 4.93e-03 -1.13454925907548640680e+04 +#> 102 1.66e-04 4.86e-03 -1.13454877292627461429e+04 +#> 103 1.66e-04 4.79e-03 -1.13454829376113357284e+04 +#> 104 1.66e-04 4.72e-03 -1.13454782141785617569e+04 +#> 105 1.66e-04 4.66e-03 -1.13454735573940706672e+04 +#> 106 1.65e-04 4.59e-03 -1.13454689657373492082e+04 +#> 107 1.65e-04 4.53e-03 -1.13454644377359963983e+04 +#> 108 1.65e-04 4.47e-03 -1.13454599719639627438e+04 +#> 109 1.64e-04 4.40e-03 -1.13454555670399531664e+04 +#> 110 1.64e-04 4.35e-03 -1.13454512216258208355e+04 +#> 111 1.64e-04 4.29e-03 -1.13454469344250337599e+04 +#> 112 1.63e-04 4.23e-03 -1.13454427041812286916e+04 +#> 113 1.62e-04 4.17e-03 -1.13454385296767795808e+04 +#> 114 1.62e-04 4.12e-03 -1.13454344097314697137e+04 +#> 115 1.61e-04 4.07e-03 -1.13454303432011420227e+04 +#> 116 1.61e-04 4.01e-03 -1.13454263289764967340e+04 +#> 117 1.60e-04 3.96e-03 -1.13454223659818417218e+04 +#> 118 1.59e-04 3.91e-03 -1.13454184531739501836e+04 +#> 119 1.58e-04 3.86e-03 -1.13454145895409474178e+04 +#> 120 1.58e-04 3.82e-03 -1.13454107741012267070e+04 +#> 121 1.57e-04 3.77e-03 -1.13454070059024234070e+04 +#> 122 1.56e-04 3.72e-03 -1.13454032840204163222e+04 +#> 123 1.55e-04 3.68e-03 -1.13453996075583818310e+04 +#> 124 1.54e-04 3.63e-03 -1.13453959756458607444e+04 +#> 125 1.53e-04 3.59e-03 -1.13453923874378815526e+04 +#> 126 1.52e-04 3.55e-03 -1.13453888421141036815e+04 +#> 127 1.52e-04 3.50e-03 -1.13453853388780153182e+04 +#> 128 1.51e-04 3.46e-03 -1.13453818769560984947e+04 +#> 129 1.50e-04 3.42e-03 -1.13453784555971269583e+04 +#> 130 1.49e-04 3.38e-03 -1.13453750740713730920e+04 +#> 131 1.48e-04 3.34e-03 -1.13453717316699530784e+04 +#> 132 1.47e-04 3.30e-03 -1.13453684277041065798e+04 +#> 133 1.46e-04 3.27e-03 -1.13453651615045473591e+04 +#> 134 1.45e-04 3.23e-03 -1.13453619324208448234e+04 +#> 135 1.44e-04 3.19e-03 -1.13453587398207855585e+04 +#> 136 1.42e-04 3.16e-03 -1.13453555830897948908e+04 +#> 137 1.41e-04 3.12e-03 -1.13453524616303784569e+04 +#> 138 1.40e-04 3.09e-03 -1.13453493748615182994e+04 +#> 139 1.39e-04 3.05e-03 -1.13453463222182090249e+04 +#> 140 1.38e-04 3.02e-03 -1.13453433031508866407e+04 +#> 141 1.37e-04 2.99e-03 -1.13453403171249356092e+04 +#> 142 1.36e-04 2.95e-03 -1.13453373636202322814e+04 +#> 143 1.35e-04 2.92e-03 -1.13453344421306501317e+04 +#> 144 1.34e-04 2.89e-03 -1.13453315521636050107e+04 +#> 145 1.33e-04 2.86e-03 -1.13453286932396367774e+04 +#> 146 1.32e-04 2.83e-03 -1.13453258648919636471e+04 +#> 147 1.31e-04 2.80e-03 -1.13453230666660747374e+04 +#> 148 1.30e-04 2.77e-03 -1.13453202981193389860e+04 +#> 149 1.29e-04 2.74e-03 -1.13453175588206013344e+04 +#> 150 1.27e-04 2.71e-03 -1.13453148483498025598e+04 +#> 151 1.26e-04 2.68e-03 -1.13453121662976591324e+04 +#> 152 1.25e-04 2.65e-03 -1.13453095122652430291e+04 +#> 153 1.24e-04 2.63e-03 -1.13453068858636943332e+04 +#> 154 1.23e-04 2.60e-03 -1.13453042867138574366e+04 +#> 155 1.22e-04 2.57e-03 -1.13453017144459627161e+04 +#> 156 1.21e-04 2.55e-03 -1.13452991686993227631e+04 +#> 157 1.20e-04 2.52e-03 -1.13452966491220304306e+04 +#> 158 1.19e-04 2.49e-03 -1.13452941553706350533e+04 +#> 159 1.18e-04 2.47e-03 -1.13452916871099005220e+04 +#> 160 1.17e-04 2.44e-03 -1.13452892440124960558e+04 +#> 161 1.16e-04 2.42e-03 -1.13452868257587324479e+04 +#> 162 1.15e-04 2.39e-03 -1.13452844320363055886e+04 +#> 163 1.14e-04 2.37e-03 -1.13452820625400509016e+04 +#> 164 1.13e-04 2.35e-03 -1.13452797169716723147e+04 +#> 165 1.12e-04 2.32e-03 -1.13452773950395385327e+04 +#> 166 1.11e-04 2.30e-03 -1.13452750964584192843e+04 +#> 167 1.10e-04 2.28e-03 -1.13452728209492961469e+04 +#> 168 1.09e-04 2.25e-03 -1.13452705682390951551e+04 +#> 169 1.08e-04 2.23e-03 -1.13452683380605449202e+04 +#> 170 1.07e-04 2.21e-03 -1.13452661301519019617e+04 +#> 171 1.06e-04 2.19e-03 -1.13452639442568179220e+04 +#> 172 1.05e-04 2.16e-03 -1.13452617801240958215e+04 +#> 173 1.04e-04 2.14e-03 -1.13452596375075281685e+04 +#> 174 1.03e-04 2.12e-03 -1.13452575161657441640e+04 +#> 175 1.03e-04 2.10e-03 -1.13452554158619823284e+04 +#> 176 1.02e-04 2.08e-03 -1.13452533363639649906e+04 +#> 177 1.01e-04 2.06e-03 -1.13452512774437182088e+04 +#> 178 9.98e-05 2.04e-03 -1.13452492388774335268e+04 #> Done! #> Processing the outputs... Done! -#> mr.mash successfully executed in 0.08278321 minutes! +#> mr.mash successfully executed in 0.1288379 minutes! # Compare the "fitted" values of Y against the true Y in the training set. plot(fit$fitted,Ytrain,pch = 20,col = "darkblue",xlab = "true", @@ -503,11 +503,11 @@

    Examples

    diff --git a/docs/reference/mr.mash.rss-1.png b/docs/reference/mr.mash.rss-1.png new file mode 100644 index 0000000..664a6f3 Binary files /dev/null and b/docs/reference/mr.mash.rss-1.png differ diff --git a/docs/reference/mr.mash.rss.html b/docs/reference/mr.mash.rss.html new file mode 100644 index 0000000..0428f47 --- /dev/null +++ b/docs/reference/mr.mash.rss.html @@ -0,0 +1,555 @@ + +Multiple Regression with Multivariate Adaptive Shrinkage + from summary data. — mr.mash.rss • mr.mash.alpha + + +
    +
    + + + +
    +
    +
    +

    Multiple Regression with Multivariate Adaptive Shrinkage + from summary data.

    + Source: R/mr_mash_rss.R + +
    + +
    +

    Performs multivariate multiple regression with + mixture-of-normals prior.

    +
    + +
    + 
    mr.mash.rss(
+  Bhat,
+  Shat,
+  Z,
+  R,
+  covY,
+  n,
+  S0,
+  w0 = rep(1/(length(S0)), length(S0)),
+  V,
+  mu1_init = NULL,
+  tol = 1e-04,
+  convergence_criterion = c("mu1", "ELBO"),
+  max_iter = 5000,
+  update_w0 = TRUE,
+  update_w0_method = "EM",
+  w0_threshold = 0,
+  compute_ELBO = TRUE,
+  standardize = FALSE,
+  verbose = TRUE,
+  update_V = FALSE,
+  update_V_method = c("full", "diagonal"),
+  version = c("Rcpp", "R"),
+  e = 1e-08,
+  ca_update_order = c("consecutive", "decreasing_logBF", "increasing_logBF", "random"),
+  X_colmeans = NULL,
+  Y_colmeans = NULL,
+  check_R = TRUE,
+  R_tol = 1e-08,
+  nthreads = as.integer(NA)
+)
    +
    + +
    +

    Arguments

    +
    Bhat
    +

    p x r matrix of regression coefficients from univariate +simple linear regression.

    + + +
    Shat
    +

    p x r matrix of standard errors of the regression coefficients +from univariate simple linear regression.

    + + +
    Z
    +

    p x r matrix of Z-scores from univariate +simple linear regression.

    + + +
    R
    +

    p x p correlation matrix among the variables.

    + + +
    covY
    +

    r x r covariance matrix across responses.

    + + +
    n
    +

    scalar indicating the sample size.

    + + +
    S0
    +

    List of length K containing the desired r x r prior +covariance matrices on the regression coefficients.

    + + +
    w0
    +

    K-vector with prior mixture weights, each associated with +the respective covariance matrix in S0.

    + + +
    V
    +

    r x r residual covariance matrix.

    + + +
    mu1_init
    +

    p x r matrix of initial estimates of the posterior +mean regression coefficients. These should be on the same scale as +the X provided. If standardize=TRUE, mu1_init will be scaled +appropriately after standardizing X.

    + + +
    tol
    +

    Convergence tolerance.

    + + +
    convergence_criterion
    +

    Criterion to use for convergence check.

    + + +
    max_iter
    +

    Maximum number of iterations for the optimization +algorithm.

    + + +
    update_w0
    +

    If TRUE, prior weights are updated.

    + + +
    update_w0_method
    +

    Method to update prior weights. Only EM is +currently supported.

    + + +
    w0_threshold
    +

    Drop mixture components with weight less than this value. +Components are dropped at each iteration after 15 initial iterations. +This is done to prevent from dropping some poetentially important +components prematurely.

    + + +
    compute_ELBO
    +

    If TRUE, ELBO is computed.

    + + +
    standardize
    +

    If TRUE, X is "standardized" using the +sample means and sample standard deviations. Standardizing X +allows a faster implementation, but the prior has a different +interpretation. Coefficients and covariances are returned on the +original scale.

    + + +
    verbose
    +

    If TRUE, some information about the +algorithm's process is printed at each iteration.

    + + +
    update_V
    +

    if TRUE, residual covariance is updated.

    + + +
    update_V_method
    +

    Method to update residual covariance. So far, +"full" and "diagonal" are supported. If update_V=TRUE and V +is not provided by the user, this option will determine how V is +computed (and fixed) internally from mu1_init.

    + + +
    version
    +

    Whether to use R or C++ code to perform the +coordinate ascent updates.

    + + +
    e
    +

    A small number to add to the diagonal elements of the +prior matrices to improve numerical stability of the updates.

    + + +
    ca_update_order
    +

    The order with which coordinates are +updated. So far, "consecutive", "decreasing_logBF", +"increasing_logBF", "random" are supported.

    + + +
    X_colmeans
    +

    a p-vector of variable means.

    + + +
    Y_colmeans
    +

    a r-vector of response means.

    + + +
    check_R
    +

    If TRUE, R is checked to be positive semidefinite.

    + + +
    R_tol
    +

    tolerance to declare positive semi-definiteness of R.

    + + +
    nthreads
    +

    Number of RcppParallel threads to use for the +updates. When nthreads is NA, the default number of +threads is used; see +defaultNumThreads. This setting is +ignored when version = "R".

    + +
    +
    +

    Value

    + + +

    A mr.mash.rss fit, stored as a list with some or all of the +following elements:

    +
    mu1
    +

    p x r matrix of posterior means for the regression + coeffcients.

    + + +
    S1
    +

    r x r x p array of posterior covariances for the + regression coeffcients.

    + + +
    w1
    +

    p x K matrix of posterior assignment probabilities to the + mixture components.

    + + +
    V
    +

    r x r residual covariance matrix

    + + +
    w0
    +

    K-vector with (updated, if update_w0=TRUE) prior mixture weights, each associated with + the respective covariance matrix in S0

    +

    .

    +
    S0
    +

    r x r x K array of prior covariance matrices + on the regression coefficients

    +

    .

    +
    intercept
    +

    r-vector containing posterior mean estimate of the + intercept, if X_colmeans and Y_colmeans are provided. + Otherwise, NA is output.

    + + +
    ELBO
    +

    Evidence Lower Bound (ELBO) at last iteration.

    + + +
    progress
    +

    A data frame including information regarding + convergence criteria at each iteration.

    + + +
    converged
    +

    TRUE or FALSE, indicating whether + the optimization algorithm converged to a solution within the chosen tolerance + level.

    + + +
    Y
    +

    n x r matrix of responses at last iteration (only relevant when missing values + are present in the input Y).

    + +
    + +
    +

    Examples

    + 
    ###Set seed
+set.seed(123)
+
+###Simulate X and Y
+##Set parameters
+n  <- 1000
+p <- 100
+p_causal <- 20
+r <- 5
+
+###Simulate data
+out <- simulate_mr_mash_data(n, p, p_causal, r, pve=0.5, B_cor=1,
+                             B_scale=1, X_cor=0, X_scale=1, V_cor=0)
+
+###Split the data in training and test sets
+Ytrain <- out$Y[-c(1:200), ]
+Xtrain <- out$X[-c(1:200), ]
+Ytest <- out$Y[c(1:200), ]
+Xtest <- out$X[c(1:200), ]
+
+###Specify the covariance matrices for the mixture-of-normals prior.
+univ_sumstats <- compute_univariate_sumstats(Xtrain, Ytrain,
+                   standardize=TRUE, standardize.response=FALSE)
+grid <- autoselect.mixsd(univ_sumstats, mult=sqrt(2))^2
+S0 <- compute_canonical_covs(ncol(Ytrain), singletons=TRUE,
+                             hetgrid=c(0, 0.25, 0.5, 0.75, 1))
+S0 <- expand_covs(S0, grid, zeromat=TRUE)
+
+###Fit mr.mash
+covY <- cov(Ytrain)
+corX <- cor(Xtrain)
+n_train <- nrow(Ytrain)
+fit <- mr.mash.rss(Bhat=univ_sumstats$Bhat, Shat=univ_sumstats$Shat, S0=S0, 
+                   covY=covY, R=corX, n=n_train, V=covY, update_V=TRUE,
+                   X_colmeans=colMeans(Xtrain), Y_colmeans=colMeans(Ytrain))
+#> Processing the inputs... Done!
+#> Fitting the optimization algorithm using 8 RcppParallel threads... 
+#>  iter    mu1_max.diff     ELBO_diff               ELBO
+#>    1       1.94e+00            Inf      -1.17846672186051746394e+04
+#>    2       4.42e-01       4.05e+02      -1.13796024155407103535e+04
+#>    3       7.87e-02       2.57e+01      -1.13539336660288754501e+04
+#>    4       2.58e-02       2.82e+00      -1.13511145940660335327e+04
+#>    5       1.32e-02       1.31e+00      -1.13498072323880351178e+04
+#>    6       8.97e-03       7.79e-01      -1.13490277541913437744e+04
+#>    7       7.04e-03       5.18e-01      -1.13485098888134161825e+04
+#>    8       5.70e-03       3.70e-01      -1.13481397456728827819e+04
+#>    9       4.72e-03       2.79e-01      -1.13478605360439578362e+04
+#>   10       3.99e-03       2.19e-01      -1.13476411019760143972e+04
+#>   11       3.42e-03       1.78e-01      -1.13474630696575877664e+04
+#>   12       2.97e-03       1.48e-01      -1.13473149527450568712e+04
+#>   13       2.61e-03       1.26e-01      -1.13471892152355649159e+04
+#>   14       2.32e-03       1.09e-01      -1.13470807097645320027e+04
+#>   15       2.07e-03       9.49e-02      -1.13469857995247657527e+04
+#>   16       1.87e-03       8.40e-02      -1.13469018398800035357e+04
+#>   17       1.69e-03       7.50e-02      -1.13468268590668321849e+04
+#>   18       1.57e-03       6.75e-02      -1.13467593542255581269e+04
+#>   19       1.48e-03       6.12e-02      -1.13466981570112002373e+04
+#>   20       1.40e-03       5.58e-02      -1.13466423427090303448e+04
+#>   21       1.36e-03       5.12e-02      -1.13465911674013350421e+04
+#>   22       1.36e-03       4.71e-02      -1.13465440236996946624e+04
+#>   23       1.35e-03       4.36e-02      -1.13465004090367347089e+04
+#>   24       1.32e-03       4.05e-02      -1.13464599026111482090e+04
+#>   25       1.30e-03       3.78e-02      -1.13464221483868677751e+04
+#>   26       1.26e-03       3.53e-02      -1.13463868423840194737e+04
+#>   27       1.22e-03       3.31e-02      -1.13463537230474557873e+04
+#>   28       1.18e-03       3.12e-02      -1.13463225638458152389e+04
+#>   29       1.14e-03       2.94e-02      -1.13462931675035342778e+04
+#>   30       1.09e-03       2.78e-02      -1.13462653614406572160e+04
+#>   31       1.05e-03       2.64e-02      -1.13462389941154560802e+04
+#>   32       1.00e-03       2.51e-02      -1.13462139320497190056e+04
+#>   33       9.57e-04       2.39e-02      -1.13461900573764942237e+04
+#>   34       9.13e-04       2.28e-02      -1.13461672657931958383e+04
+#>   35       8.69e-04       2.18e-02      -1.13461454648335620732e+04
+#>   36       8.27e-04       2.09e-02      -1.13461245723942793120e+04
+#>   37       7.87e-04       2.01e-02      -1.13461045154679777625e+04
+#>   38       7.48e-04       1.93e-02      -1.13460852290460861695e+04
+#>   39       7.10e-04       1.86e-02      -1.13460666551632675692e+04
+#>   40       6.75e-04       1.79e-02      -1.13460487420616536838e+04
+#>   41       6.41e-04       1.73e-02      -1.13460314434572974278e+04
+#>   42       6.09e-04       1.67e-02      -1.13460147178949864610e+04
+#>   43       5.78e-04       1.62e-02      -1.13459985281798235519e+04
+#>   44       5.50e-04       1.57e-02      -1.13459828408760222374e+04
+#>   45       5.23e-04       1.52e-02      -1.13459676258648323710e+04
+#>   46       4.97e-04       1.48e-02      -1.13459528559546197357e+04
+#>   47       4.73e-04       1.43e-02      -1.13459385065370825032e+04
+#>   48       4.51e-04       1.40e-02      -1.13459245552842985489e+04
+#>   49       4.30e-04       1.36e-02      -1.13459109818820052169e+04
+#>   50       4.11e-04       1.32e-02      -1.13458977677949133067e+04
+#>   51       3.98e-04       1.29e-02      -1.13458848960604336753e+04
+#>   52       3.85e-04       1.25e-02      -1.13458723511075259012e+04
+#>   53       3.72e-04       1.22e-02      -1.13458601185977368004e+04
+#>   54       3.61e-04       1.19e-02      -1.13458481852858749335e+04
+#>   55       3.49e-04       1.16e-02      -1.13458365388979145791e+04
+#>   56       3.38e-04       1.14e-02      -1.13458251680241428403e+04
+#>   57       3.28e-04       1.11e-02      -1.13458140620256526745e+04
+#>   58       3.18e-04       1.09e-02      -1.13458032109525793203e+04
+#>   59       3.08e-04       1.06e-02      -1.13457926054725903668e+04
+#>   60       2.99e-04       1.04e-02      -1.13457822368083907350e+04
+#>   61       2.90e-04       1.01e-02      -1.13457720966830856923e+04
+#>   62       2.82e-04       9.92e-03      -1.13457621772723759932e+04
+#>   63       2.73e-04       9.71e-03      -1.13457524711627756915e+04
+#>   64       2.66e-04       9.50e-03      -1.13457429713149795134e+04
+#>   65       2.58e-04       9.30e-03      -1.13457336710317813413e+04
+#>   66       2.51e-04       9.11e-03      -1.13457245639299435425e+04
+#>   67       2.44e-04       8.92e-03      -1.13457156439154114196e+04
+#>   68       2.37e-04       8.74e-03      -1.13457069051615508215e+04
+#>   69       2.30e-04       8.56e-03      -1.13456983420898977784e+04
+#>   70       2.24e-04       8.39e-03      -1.13456899493531509506e+04
+#>   71       2.17e-04       8.23e-03      -1.13456817218200649222e+04
+#>   72       2.11e-04       8.07e-03      -1.13456736545620042307e+04
+#>   73       2.06e-04       7.91e-03      -1.13456657428409307613e+04
+#>   74       2.00e-04       7.76e-03      -1.13456579820986353297e+04
+#>   75       1.95e-04       7.61e-03      -1.13456503679470242787e+04
+#>   76       1.89e-04       7.47e-03      -1.13456428961593537679e+04
+#>   77       1.84e-04       7.33e-03      -1.13456355626622389536e+04
+#>   78       1.79e-04       7.20e-03      -1.13456283635283762123e+04
+#>   79       1.74e-04       7.07e-03      -1.13456212949698528973e+04
+#>   80       1.70e-04       6.94e-03      -1.13456143533320100687e+04
+#>   81       1.66e-04       6.82e-03      -1.13456075350877017627e+04
+#>   82       1.63e-04       6.70e-03      -1.13456008368320217414e+04
+#>   83       1.60e-04       6.58e-03      -1.13455942552773249190e+04
+#>   84       1.56e-04       6.47e-03      -1.13455877872485780244e+04
+#>   85       1.57e-04       6.36e-03      -1.13455814296789812943e+04
+#>   86       1.59e-04       6.25e-03      -1.13455751796058266336e+04
+#>   87       1.60e-04       6.15e-03      -1.13455690341665595042e+04
+#>   88       1.61e-04       6.04e-03      -1.13455629905950590910e+04
+#>   89       1.62e-04       5.94e-03      -1.13455570462180639879e+04
+#>   90       1.62e-04       5.85e-03      -1.13455511984517888777e+04
+#>   91       1.63e-04       5.75e-03      -1.13455454447986412561e+04
+#>   92       1.64e-04       5.66e-03      -1.13455397828441164165e+04
+#>   93       1.64e-04       5.57e-03      -1.13455342102537997562e+04
+#>   94       1.65e-04       5.49e-03      -1.13455287247704764013e+04
+#>   95       1.65e-04       5.40e-03      -1.13455233242113608867e+04
+#>   96       1.65e-04       5.32e-03      -1.13455180064654359740e+04
+#>   97       1.65e-04       5.24e-03      -1.13455127694908624107e+04
+#>   98       1.66e-04       5.16e-03      -1.13455076113125232951e+04
+#>   99       1.66e-04       5.08e-03      -1.13455025300196248281e+04
+#>  100       1.66e-04       5.01e-03      -1.13454975237633934739e+04
+#>  101       1.66e-04       4.93e-03      -1.13454925907548658870e+04
+#>  102       1.66e-04       4.86e-03      -1.13454877292627425049e+04
+#>  103       1.66e-04       4.79e-03      -1.13454829376113320905e+04
+#>  104       1.66e-04       4.72e-03      -1.13454782141785635758e+04
+#>  105       1.66e-04       4.66e-03      -1.13454735573940706672e+04
+#>  106       1.65e-04       4.59e-03      -1.13454689657373510272e+04
+#>  107       1.65e-04       4.53e-03      -1.13454644377359927603e+04
+#>  108       1.65e-04       4.47e-03      -1.13454599719639609248e+04
+#>  109       1.64e-04       4.40e-03      -1.13454555670399549854e+04
+#>  110       1.64e-04       4.35e-03      -1.13454512216258190165e+04
+#>  111       1.64e-04       4.29e-03      -1.13454469344250355789e+04
+#>  112       1.63e-04       4.23e-03      -1.13454427041812286916e+04
+#>  113       1.62e-04       4.17e-03      -1.13454385296767795808e+04
+#>  114       1.62e-04       4.12e-03      -1.13454344097314642568e+04
+#>  115       1.61e-04       4.07e-03      -1.13454303432011420227e+04
+#>  116       1.61e-04       4.01e-03      -1.13454263289764949150e+04
+#>  117       1.60e-04       3.96e-03      -1.13454223659818380838e+04
+#>  118       1.59e-04       3.91e-03      -1.13454184531739501836e+04
+#>  119       1.58e-04       3.86e-03      -1.13454145895409455989e+04
+#>  120       1.58e-04       3.82e-03      -1.13454107741012248880e+04
+#>  121       1.57e-04       3.77e-03      -1.13454070059024215880e+04
+#>  122       1.56e-04       3.72e-03      -1.13454032840204181412e+04
+#>  123       1.55e-04       3.68e-03      -1.13453996075583800121e+04
+#>  124       1.54e-04       3.63e-03      -1.13453959756458571064e+04
+#>  125       1.53e-04       3.59e-03      -1.13453923874378760956e+04
+#>  126       1.52e-04       3.55e-03      -1.13453888421141036815e+04
+#>  127       1.52e-04       3.50e-03      -1.13453853388780134992e+04
+#>  128       1.51e-04       3.46e-03      -1.13453818769560984947e+04
+#>  129       1.50e-04       3.42e-03      -1.13453784555971269583e+04
+#>  130       1.49e-04       3.38e-03      -1.13453750740713730920e+04
+#>  131       1.48e-04       3.34e-03      -1.13453717316699530784e+04
+#>  132       1.47e-04       3.30e-03      -1.13453684277041047608e+04
+#>  133       1.46e-04       3.27e-03      -1.13453651615045455401e+04
+#>  134       1.45e-04       3.23e-03      -1.13453619324208448234e+04
+#>  135       1.44e-04       3.19e-03      -1.13453587398207819206e+04
+#>  136       1.42e-04       3.16e-03      -1.13453555830897948908e+04
+#>  137       1.41e-04       3.12e-03      -1.13453524616303748189e+04
+#>  138       1.40e-04       3.09e-03      -1.13453493748615182994e+04
+#>  139       1.39e-04       3.05e-03      -1.13453463222182090249e+04
+#>  140       1.38e-04       3.02e-03      -1.13453433031508866407e+04
+#>  141       1.37e-04       2.99e-03      -1.13453403171249356092e+04
+#>  142       1.36e-04       2.95e-03      -1.13453373636202322814e+04
+#>  143       1.35e-04       2.92e-03      -1.13453344421306464938e+04
+#>  144       1.34e-04       2.89e-03      -1.13453315521636031917e+04
+#>  145       1.33e-04       2.86e-03      -1.13453286932396331395e+04
+#>  146       1.32e-04       2.83e-03      -1.13453258648919636471e+04
+#>  147       1.31e-04       2.80e-03      -1.13453230666660765564e+04
+#>  148       1.30e-04       2.77e-03      -1.13453202981193371670e+04
+#>  149       1.29e-04       2.74e-03      -1.13453175588205958775e+04
+#>  150       1.27e-04       2.71e-03      -1.13453148483498025598e+04
+#>  151       1.26e-04       2.68e-03      -1.13453121662976573134e+04
+#>  152       1.25e-04       2.65e-03      -1.13453095122652466671e+04
+#>  153       1.24e-04       2.63e-03      -1.13453068858636925142e+04
+#>  154       1.23e-04       2.60e-03      -1.13453042867138537986e+04
+#>  155       1.22e-04       2.57e-03      -1.13453017144459627161e+04
+#>  156       1.21e-04       2.55e-03      -1.13452991686993227631e+04
+#>  157       1.20e-04       2.52e-03      -1.13452966491220286116e+04
+#>  158       1.19e-04       2.49e-03      -1.13452941553706368722e+04
+#>  159       1.18e-04       2.47e-03      -1.13452916871099023410e+04
+#>  160       1.17e-04       2.44e-03      -1.13452892440124960558e+04
+#>  161       1.16e-04       2.42e-03      -1.13452868257587288099e+04
+#>  162       1.15e-04       2.39e-03      -1.13452844320363037696e+04
+#>  163       1.14e-04       2.37e-03      -1.13452820625400454446e+04
+#>  164       1.13e-04       2.35e-03      -1.13452797169716723147e+04
+#>  165       1.12e-04       2.32e-03      -1.13452773950395367137e+04
+#>  166       1.11e-04       2.30e-03      -1.13452750964584192843e+04
+#>  167       1.10e-04       2.28e-03      -1.13452728209492906899e+04
+#>  168       1.09e-04       2.25e-03      -1.13452705682390969741e+04
+#>  169       1.08e-04       2.23e-03      -1.13452683380605431012e+04
+#>  170       1.07e-04       2.21e-03      -1.13452661301519037806e+04
+#>  171       1.06e-04       2.19e-03      -1.13452639442568161030e+04
+#>  172       1.05e-04       2.16e-03      -1.13452617801240921835e+04
+#>  173       1.04e-04       2.14e-03      -1.13452596375075318065e+04
+#>  174       1.03e-04       2.12e-03      -1.13452575161657441640e+04
+#>  175       1.03e-04       2.10e-03      -1.13452554158619823284e+04
+#>  176       1.02e-04       2.08e-03      -1.13452533363639613526e+04
+#>  177       1.01e-04       2.06e-03      -1.13452512774437182088e+04
+#>  178       9.98e-05       2.04e-03      -1.13452492388774317078e+04
+#> Done!
+#> Processing the outputs... Done!
+#> mr.mash.rss successfully executed in 0.1818886 minutes!
+
+# Predict the multivariate outcomes in the test set using the fitted model.
+Ytest_est <- predict(fit,Xtest)
+plot(Ytest_est,Ytest,pch = 20,col = "darkblue",xlab = "true",
+     ylab = "predicted")
+abline(a = 0,b = 1,col = "magenta",lty = "dotted")
+
+
+
    +
    +
    + +
    + + +
    +

    Developed by Fabio Morgante, Deborah Kunkel, Peter Carbonetto, Matthew Stephens.

    +
    + +
    +

    Site built with pkgdown 2.0.7.

    +
    + +
    + + + + + + + + diff --git a/docs/reference/predict.mr.mash.html b/docs/reference/predict.mr.mash.html index ae9182b..a6450c1 100644 --- a/docs/reference/predict.mr.mash.html +++ b/docs/reference/predict.mr.mash.html @@ -17,7 +17,7 @@ mr.mash.alpha - 0.2-27 + 0.3.21 @@ -86,11 +86,11 @@

    Value

    diff --git a/docs/reference/predict.mr.mash.rss.html b/docs/reference/predict.mr.mash.rss.html new file mode 100644 index 0000000..0cbf7ce --- /dev/null +++ b/docs/reference/predict.mr.mash.rss.html @@ -0,0 +1,104 @@ + +Predict future observations from mr.mash.rss fit. — predict.mr.mash.rss • mr.mash.alpha + + +
    +
    + + + +
    +
    +
    +

    Predict future observations from mr.mash.rss fit.

    + Source: R/mr_mash_rss_predict.R + +
    + +
    +

    Predict future observations from mr.mash.rss fit.

    +
    + +
    + 
    # S3 method for mr.mash.rss
+predict(object, newx, ...)
    +
    + +
    +

    Arguments

    +
    object
    +

    a mr.mash.rss fit.

    + + +
    newx
    +

    a new value for X at which to do predictions.

    + + +
    ...
    +

    Additional arguments (not used).

    + +
    +
    +

    Value

    + + +

    Matrix of predicted values.

    +
    + +
    + +
    + + +
    +

    Developed by Fabio Morgante, Deborah Kunkel, Peter Carbonetto, Matthew Stephens.

    +
    + +
    +

    Site built with pkgdown 2.0.7.

    +
    + +
    + + + + + + + + diff --git a/docs/reference/simulate_mr_mash_data.html b/docs/reference/simulate_mr_mash_data.html index bd98f04..dabdf2e 100644 --- a/docs/reference/simulate_mr_mash_data.html +++ b/docs/reference/simulate_mr_mash_data.html @@ -18,7 +18,7 @@ mr.mash.alpha - 0.2-27 + 0.3.21 @@ -181,7 +181,8 @@

    Value

    Examples

    - 
    dat <- simulate_mr_mash_data(n=50, p=40, p_causal=20, r=5,
+    
    set.seed(1)
+dat <- simulate_mr_mash_data(n=50, p=40, p_causal=20, r=5,
                              r_causal=list(1:2, 3:4), intercepts=rep(1, 5),
                              pve=0.2, B_cor=c(0, 1), B_scale=c(0.5, 1),
                              w=c(0.5, 0.5), X_cor=0.5, X_scale=1, 
@@ -198,11 +199,11 @@ 
    Examples
    
 
 
       
    
-  
    Developed by Fabio Morgante, Peter Carbonetto, Matthew Stephens.
    
+  
    Developed by Fabio Morgante, Deborah Kunkel, Peter Carbonetto, Matthew Stephens.
    
 
    
 
 
    
-  
    Site built with pkgdown 2.0.6.
    
+  
    Site built with pkgdown 2.0.7.
    
 
    
 
       
    
diff --git a/docs/sitemap.xml b/docs/sitemap.xml
index dd31b6a..8a03990 100644
--- a/docs/sitemap.xml
+++ b/docs/sitemap.xml
@@ -6,6 +6,9 @@
   
     /LICENSE-text.html
   
+  
+    /LICENSE.html
+  
   
     /articles/index.html
   
@@ -24,6 +27,9 @@
   
     /reference/coef.mr.mash.html
   
+  
+    /reference/coef.mr.mash.rss.html
+  
   
     /reference/compute_canonical_covs.html
   
@@ -42,9 +48,15 @@
   
     /reference/mr.mash.html
   
+  
+    /reference/mr.mash.rss.html
+  
   
     /reference/predict.mr.mash.html
   
+  
+    /reference/predict.mr.mash.rss.html
+  
   
     /reference/simulate_mr_mash_data.html
   
diff --git a/man/coef.mr.mash.rss.Rd b/man/coef.mr.mash.rss.Rd
new file mode 100644
index 0000000..42f9fee
--- /dev/null
+++ b/man/coef.mr.mash.rss.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mr_mash_rss_predict.R
+\name{coef.mr.mash.rss}
+\alias{coef.mr.mash.rss}
+\title{Extract coefficients from mr.mash.rss fit.}
+\usage{
+\method{coef}{mr.mash.rss}(object, ...)
+}
+\arguments{
+\item{object}{a mr.mash fit.}
+
+\item{\dots}{Other arguments (not used).}
+}
+\value{
+p x r or (p+1) x r matrix of coefficients,
+ depending on whether an intercept was computed.
+}
+\description{
+Extract coefficients from mr.mash.rss fit.
+}
diff --git a/man/mr.mash.Rd b/man/mr.mash.Rd
index 9f0c75f..8d6cb59 100644
--- a/man/mr.mash.Rd
+++ b/man/mr.mash.Rd
@@ -11,7 +11,7 @@ mr.mash(
   w0 = rep(1/(length(S0)), length(S0)),
   V = NULL,
   mu1_init = matrix(0, nrow = ncol(X), ncol = ncol(Y)),
-  tol = 0.0001,
+  tol = 1e-04,
   convergence_criterion = c("mu1", "ELBO"),
   max_iter = 5000,
   update_w0 = TRUE,
@@ -24,7 +24,7 @@ mr.mash(
   update_V_method = c("full", "diagonal"),
   version = c("Rcpp", "R"),
   e = 1e-08,
-  ca_update_order = c("consecutive", "decreasing_logBF", "increasing_logBF"),
+  ca_update_order = c("consecutive", "decreasing_logBF", "increasing_logBF", "random"),
   nthreads = as.integer(NA)
 )
 }
@@ -89,7 +89,7 @@ prior matrices to improve numerical stability of the updates.}
 
 \item{ca_update_order}{The order with which coordinates are
 updated.  So far, "consecutive", "decreasing_logBF",
-"increasing_logBF" are supported.}
+"increasing_logBF", "random" are supported.}
 
 \item{nthreads}{Number of RcppParallel threads to use for the
 updates. When \code{nthreads} is \code{NA}, the default number of
diff --git a/man/mr.mash.rss.Rd b/man/mr.mash.rss.Rd
new file mode 100644
index 0000000..ed6934b
--- /dev/null
+++ b/man/mr.mash.rss.Rd
@@ -0,0 +1,212 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mr_mash_rss.R
+\name{mr.mash.rss}
+\alias{mr.mash.rss}
+\title{Multiple Regression with Multivariate Adaptive Shrinkage
+ from summary data.}
+\usage{
+mr.mash.rss(
+  Bhat,
+  Shat,
+  Z,
+  R,
+  covY,
+  n,
+  S0,
+  w0 = rep(1/(length(S0)), length(S0)),
+  V,
+  mu1_init = NULL,
+  tol = 1e-04,
+  convergence_criterion = c("mu1", "ELBO"),
+  max_iter = 5000,
+  update_w0 = TRUE,
+  update_w0_method = "EM",
+  w0_threshold = 0,
+  compute_ELBO = TRUE,
+  standardize = FALSE,
+  verbose = TRUE,
+  update_V = FALSE,
+  update_V_method = c("full", "diagonal"),
+  version = c("Rcpp", "R"),
+  e = 1e-08,
+  ca_update_order = c("consecutive", "decreasing_logBF", "increasing_logBF", "random"),
+  X_colmeans = NULL,
+  Y_colmeans = NULL,
+  check_R = TRUE,
+  R_tol = 1e-08,
+  nthreads = as.integer(NA)
+)
+}
+\arguments{
+\item{Bhat}{p x r matrix of regression coefficients from univariate
+simple linear regression.}
+
+\item{Shat}{p x r matrix of standard errors of the regression coefficients 
+from univariate simple linear regression.}
+
+\item{Z}{p x r matrix of Z-scores from univariate
+simple linear regression.}
+
+\item{R}{p x p correlation matrix among the variables.}
+
+\item{covY}{r x r covariance matrix across responses.}
+
+\item{n}{scalar indicating the sample size.}
+
+\item{S0}{List of length K containing the desired r x r prior
+covariance matrices on the regression coefficients.}
+
+\item{w0}{K-vector with prior mixture weights, each associated with
+the respective covariance matrix in \code{S0}.}
+
+\item{V}{r x r residual covariance matrix.}
+
+\item{mu1_init}{p x r matrix of initial estimates of the posterior
+mean regression coefficients. These should be on the same scale as
+the X provided. If \code{standardize=TRUE}, mu1_init will be scaled
+appropriately after standardizing X.}
+
+\item{tol}{Convergence tolerance.}
+
+\item{convergence_criterion}{Criterion to use for convergence check.}
+
+\item{max_iter}{Maximum number of iterations for the optimization
+algorithm.}
+
+\item{update_w0}{If \code{TRUE}, prior weights are updated.}
+
+\item{update_w0_method}{Method to update prior weights. Only EM is
+currently supported.}
+
+\item{w0_threshold}{Drop mixture components with weight less than this value.
+Components are dropped at each iteration after 15 initial iterations.
+This is done to prevent from dropping some poetentially important 
+components prematurely.}
+
+\item{compute_ELBO}{If \code{TRUE}, ELBO is computed.}
+
+\item{standardize}{If \code{TRUE}, X is "standardized" using the
+sample means and sample standard deviations. Standardizing X
+allows a faster implementation, but the prior has a different
+interpretation. Coefficients and covariances are returned on the
+original scale.}
+
+\item{verbose}{If \code{TRUE}, some information about the
+algorithm's process is printed at each iteration.}
+
+\item{update_V}{if \code{TRUE}, residual covariance is updated.}
+
+\item{update_V_method}{Method to update residual covariance. So far,
+"full" and "diagonal" are supported. If \code{update_V=TRUE} and V 
+is not provided by the user, this option will determine how V is 
+computed (and fixed) internally from \code{mu1_init}.}
+
+\item{version}{Whether to use R or C++ code to perform the
+coordinate ascent updates.}
+
+\item{e}{A small number to add to the diagonal elements of the
+prior matrices to improve numerical stability of the updates.}
+
+\item{ca_update_order}{The order with which coordinates are
+updated.  So far, "consecutive", "decreasing_logBF",
+"increasing_logBF", "random" are supported.}
+
+\item{X_colmeans}{a p-vector of variable means.}
+
+\item{Y_colmeans}{a r-vector of response means.}
+
+\item{check_R}{If \code{TRUE}, R is checked to be positive semidefinite.}
+
+\item{R_tol}{tolerance to declare positive semi-definiteness of R.}
+
+\item{nthreads}{Number of RcppParallel threads to use for the
+updates. When \code{nthreads} is \code{NA}, the default number of
+threads is used; see
+\code{\link[RcppParallel]{defaultNumThreads}}. This setting is
+ignored when \code{version = "R"}.}
+}
+\value{
+A mr.mash.rss fit, stored as a list with some or all of the
+following elements:
+
+\item{mu1}{p x r matrix of posterior means for the regression
+  coeffcients.}
+
+\item{S1}{r x r x p array of posterior covariances for the
+  regression coeffcients.}
+
+\item{w1}{p x K matrix of posterior assignment probabilities to the
+  mixture components.}
+  
+\item{V}{r x r residual covariance matrix}
+
+\item{w0}{K-vector with (updated, if \code{update_w0=TRUE}) prior mixture weights, each associated with
+  the respective covariance matrix in \code{S0}}.
+  
+\item{S0}{r x r x K array of prior covariance matrices
+  on the regression coefficients}.
+
+\item{intercept}{r-vector containing posterior mean estimate of the
+  intercept, if \code{X_colmeans} and \code{Y_colmeans} are provided.
+  Otherwise, \code{NA} is output.}
+
+\item{ELBO}{Evidence Lower Bound (ELBO) at last iteration.}
+
+\item{progress}{A data frame including information regarding
+  convergence criteria at each iteration.}
+  
+\item{converged}{\code{TRUE} or \code{FALSE}, indicating whether
+  the optimization algorithm converged to a solution within the chosen tolerance
+  level.}
+  
+\item{Y}{n x r matrix of responses at last iteration (only relevant when missing values
+  are present in the input Y).}
+}
+\description{
+Performs multivariate multiple regression with
+  mixture-of-normals prior.
+}
+\examples{
+###Set seed
+set.seed(123)
+
+###Simulate X and Y
+##Set parameters
+n  <- 1000
+p <- 100
+p_causal <- 20
+r <- 5
+
+###Simulate data
+out <- simulate_mr_mash_data(n, p, p_causal, r, pve=0.5, B_cor=1,
+                             B_scale=1, X_cor=0, X_scale=1, V_cor=0)
+
+###Split the data in training and test sets
+Ytrain <- out$Y[-c(1:200), ]
+Xtrain <- out$X[-c(1:200), ]
+Ytest <- out$Y[c(1:200), ]
+Xtest <- out$X[c(1:200), ]
+
+###Specify the covariance matrices for the mixture-of-normals prior.
+univ_sumstats <- compute_univariate_sumstats(Xtrain, Ytrain,
+                   standardize=TRUE, standardize.response=FALSE)
+grid <- autoselect.mixsd(univ_sumstats, mult=sqrt(2))^2
+S0 <- compute_canonical_covs(ncol(Ytrain), singletons=TRUE,
+                             hetgrid=c(0, 0.25, 0.5, 0.75, 1))
+S0 <- expand_covs(S0, grid, zeromat=TRUE)
+
+###Fit mr.mash
+covY <- cov(Ytrain)
+corX <- cor(Xtrain)
+n_train <- nrow(Ytrain)
+fit <- mr.mash.rss(Bhat=univ_sumstats$Bhat, Shat=univ_sumstats$Shat, S0=S0, 
+                   covY=covY, R=corX, n=n_train, V=covY, update_V=TRUE,
+                   X_colmeans=colMeans(Xtrain), Y_colmeans=colMeans(Ytrain))
+
+# Predict the multivariate outcomes in the test set using the fitted model.
+Ytest_est <- predict(fit,Xtest)
+plot(Ytest_est,Ytest,pch = 20,col = "darkblue",xlab = "true",
+     ylab = "predicted")
+abline(a = 0,b = 1,col = "magenta",lty = "dotted")
+
+}
diff --git a/man/predict.mr.mash.rss.Rd b/man/predict.mr.mash.rss.Rd
new file mode 100644
index 0000000..193cbc1
--- /dev/null
+++ b/man/predict.mr.mash.rss.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mr_mash_rss_predict.R
+\name{predict.mr.mash.rss}
+\alias{predict.mr.mash.rss}
+\title{Predict future observations from mr.mash.rss fit.}
+\usage{
+\method{predict}{mr.mash.rss}(object, newx, ...)
+}
+\arguments{
+\item{object}{a mr.mash.rss fit.}
+
+\item{newx}{a new value for X at which to do predictions.}
+
+\item{\dots}{Additional arguments (not used).}
+}
+\value{
+Matrix of predicted values.
+}
+\description{
+Predict future observations from mr.mash.rss fit.
+}
diff --git a/man/simulate_mr_mash_data.Rd b/man/simulate_mr_mash_data.Rd
index 3e8b8e2..eee76b8 100644
--- a/man/simulate_mr_mash_data.Rd
+++ b/man/simulate_mr_mash_data.Rd
@@ -82,6 +82,7 @@ Function to simulate data from \eqn{MN_{nxr}(XB, I, V)}, where \eqn{X \sim N_p(0
   \eqn{B \sim \sum_k w_k N_r(0, Sigma_k)}, with \eqn{Gamma}, \eqn{w_k}, \eqn{Sigma_k}, and \eqn{V} defined by the user.
 }
 \examples{
+set.seed(1)
 dat <- simulate_mr_mash_data(n=50, p=40, p_causal=20, r=5,
                              r_causal=list(1:2, 3:4), intercepts=rep(1, 5),
                              pve=0.2, B_cor=c(0, 1), B_scale=c(0.5, 1),
diff --git a/src/Makevars b/src/Makevars
index 47251da..4bb1874 100644
--- a/src/Makevars
+++ b/src/Makevars
@@ -1,5 +1,6 @@
 CXX_STD = CXX11
-PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) -DARMA_WARN_LEVEL=1 \
-               -DARMA_NO_DEBUG -DARMA_DONT_USE_OPENMP -DARMA_USE_TBB_ALLOC
+PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) -DARMA_WARN_LEVEL=1 -DARMA_64BIT_WORD \
+               -DARMA_DONT_USE_OPENMP -DARMA_USE_TBB_ALLOC -DARMA_NO_DEBUG \
+               -DARMA_USE_BLAS -DARMA_USE_LAPACK -DRCPP_PARALLEL_USE_TBB=1
 PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) \
            $(shell ${R_HOME}/bin/Rscript -e "RcppParallel::RcppParallelLibs()")
diff --git a/src/Makevars.win b/src/Makevars.win
index fb8175a..77b5b0c 100644
--- a/src/Makevars.win
+++ b/src/Makevars.win
@@ -1,6 +1,7 @@
 CXX_STD = CXX11
 PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) -DARMA_WARN_LEVEL=1 \
-               -DARMA_NO_DEBUG -DARMA_DONT_USE_OPENMP \
-	       -DARMA_USE_TBB_ALLOC -DRCPP_PARALLEL_USE_TBB=1
+               -DARMA_64BIT_WORD -DARMA_DONT_USE_OPENMP \
+	       -DARMA_USE_TBB_ALLOC -DRCPP_PARALLEL_USE_TBB=1 \
+               -DARMA_NO_DEBUG -DARMA_USE_BLAS -DARMA_USE_LAPACK
 PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) \
            $(shell "${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" -e "RcppParallel::RcppParallelLibs()")
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index a8c4abd..6b5fae5 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -50,6 +50,32 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// inner_loop_general_rss_rcpp
+List inner_loop_general_rss_rcpp(unsigned int n, const arma::mat& XtX, const arma::mat& XtY, arma::mat& XtRbar, arma::mat& mu1, const arma::mat& V, const arma::mat& Vinv, const arma::vec& w0, const arma::cube& S0, const List& precomp_quants_list, bool standardize, bool compute_ELBO, bool update_V, const arma::vec& update_order, double eps, unsigned int nthreads);
+RcppExport SEXP _mr_mash_alpha_inner_loop_general_rss_rcpp(SEXP nSEXP, SEXP XtXSEXP, SEXP XtYSEXP, SEXP XtRbarSEXP, SEXP mu1SEXP, SEXP VSEXP, SEXP VinvSEXP, SEXP w0SEXP, SEXP S0SEXP, SEXP precomp_quants_listSEXP, SEXP standardizeSEXP, SEXP compute_ELBOSEXP, SEXP update_VSEXP, SEXP update_orderSEXP, SEXP epsSEXP, SEXP nthreadsSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< unsigned int >::type n(nSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type XtX(XtXSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type XtY(XtYSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type XtRbar(XtRbarSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type mu1(mu1SEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type V(VSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type Vinv(VinvSEXP);
+    Rcpp::traits::input_parameter< const arma::vec& >::type w0(w0SEXP);
+    Rcpp::traits::input_parameter< const arma::cube& >::type S0(S0SEXP);
+    Rcpp::traits::input_parameter< const List& >::type precomp_quants_list(precomp_quants_listSEXP);
+    Rcpp::traits::input_parameter< bool >::type standardize(standardizeSEXP);
+    Rcpp::traits::input_parameter< bool >::type compute_ELBO(compute_ELBOSEXP);
+    Rcpp::traits::input_parameter< bool >::type update_V(update_VSEXP);
+    Rcpp::traits::input_parameter< const arma::vec& >::type update_order(update_orderSEXP);
+    Rcpp::traits::input_parameter< double >::type eps(epsSEXP);
+    Rcpp::traits::input_parameter< unsigned int >::type nthreads(nthreadsSEXP);
+    rcpp_result_gen = Rcpp::wrap(inner_loop_general_rss_rcpp(n, XtX, XtY, XtRbar, mu1, V, Vinv, w0, S0, precomp_quants_list, standardize, compute_ELBO, update_V, update_order, eps, nthreads));
+    return rcpp_result_gen;
+END_RCPP
+}
 // scale_rcpp
 arma::mat scale_rcpp(const arma::mat& M, const arma::vec& a, const arma::vec& b);
 RcppExport SEXP _mr_mash_alpha_scale_rcpp(SEXP MSEXP, SEXP aSEXP, SEXP bSEXP) {
@@ -109,14 +135,36 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// compute_logbf_rss_rcpp
+arma::vec compute_logbf_rss_rcpp(unsigned int n, const arma::mat& XtY, const arma::mat& V, const arma::mat& Vinv, const arma::vec& w0, const arma::cube& S0, const List& precomp_quants_list, bool standardize, double eps, unsigned int nthreads);
+RcppExport SEXP _mr_mash_alpha_compute_logbf_rss_rcpp(SEXP nSEXP, SEXP XtYSEXP, SEXP VSEXP, SEXP VinvSEXP, SEXP w0SEXP, SEXP S0SEXP, SEXP precomp_quants_listSEXP, SEXP standardizeSEXP, SEXP epsSEXP, SEXP nthreadsSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< unsigned int >::type n(nSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type XtY(XtYSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type V(VSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type Vinv(VinvSEXP);
+    Rcpp::traits::input_parameter< const arma::vec& >::type w0(w0SEXP);
+    Rcpp::traits::input_parameter< const arma::cube& >::type S0(S0SEXP);
+    Rcpp::traits::input_parameter< const List& >::type precomp_quants_list(precomp_quants_listSEXP);
+    Rcpp::traits::input_parameter< bool >::type standardize(standardizeSEXP);
+    Rcpp::traits::input_parameter< double >::type eps(epsSEXP);
+    Rcpp::traits::input_parameter< unsigned int >::type nthreads(nthreadsSEXP);
+    rcpp_result_gen = Rcpp::wrap(compute_logbf_rss_rcpp(n, XtY, V, Vinv, w0, S0, precomp_quants_list, standardize, eps, nthreads));
+    return rcpp_result_gen;
+END_RCPP
+}
 
 static const R_CallMethodDef CallEntries[] = {
     {"_mr_mash_alpha_inner_loop_general_rcpp", (DL_FUNC) &_mr_mash_alpha_inner_loop_general_rcpp, 14},
     {"_mr_mash_alpha_impute_missing_Y_rcpp", (DL_FUNC) &_mr_mash_alpha_impute_missing_Y_rcpp, 5},
+    {"_mr_mash_alpha_inner_loop_general_rss_rcpp", (DL_FUNC) &_mr_mash_alpha_inner_loop_general_rss_rcpp, 16},
     {"_mr_mash_alpha_scale_rcpp", (DL_FUNC) &_mr_mash_alpha_scale_rcpp, 3},
     {"_mr_mash_alpha_scale2_rcpp", (DL_FUNC) &_mr_mash_alpha_scale2_rcpp, 3},
     {"_mr_mash_alpha_rescale_post_mean_covar_rcpp", (DL_FUNC) &_mr_mash_alpha_rescale_post_mean_covar_rcpp, 3},
     {"_mr_mash_alpha_compute_logbf_rcpp", (DL_FUNC) &_mr_mash_alpha_compute_logbf_rcpp, 10},
+    {"_mr_mash_alpha_compute_logbf_rss_rcpp", (DL_FUNC) &_mr_mash_alpha_compute_logbf_rss_rcpp, 10},
     {NULL, NULL, 0}
 };
 
diff --git a/src/bayes_reg_mv.cpp b/src/bayes_reg_mv.cpp
index 9081781..c660229 100644
--- a/src/bayes_reg_mv.cpp
+++ b/src/bayes_reg_mv.cpp
@@ -319,3 +319,122 @@ double bayes_mvr_mix_centered_X (const vec& x, const mat& Y, const mat& V,
   double u = max(logbfmix);
   return u + log(sum(exp(logbfmix - u)));
 }
+
+
+// Perform Bayesian multivariate simple regression with summary data
+// and mixture-of-normals prior with standardized x.
+double bayes_mvr_mix_standardized_X_rss (unsigned int n, const vec& xtY, const vec& w0,
+                                         const cube& S0, const mat& S,
+                                         const cube& S1, const cube& SplusS0_chol,
+                                         const mat& S_chol, double eps,
+                                         unsigned int nthreads, vec& mu1_mix,
+                                         mat& S1_mix, vec& w1) {
+  unsigned int k = w0.n_elem;
+  unsigned int r = xtY.n_elem;
+  
+  mat mu1mix(r,k);
+  vec logbfmix(k);
+  vec mu1(r);
+  
+  // Compute the least-squares estimate.
+  vec b = xtY/(n-1);
+  
+  // Compute the quantities separately for each mixture component.
+  if (nthreads > 1) {
+    bayes_mvr_mix_standardized_X_worker worker(b,S,S_chol,S0,S1,SplusS0_chol,
+                                                   logbfmix,mu1mix);
+    parallelFor(0,k,worker);
+  } else {
+    for (unsigned int i = 0; i < k; i++) {
+      logbfmix(i) = bayes_mvr_ridge_standardized_X(b,S0.slice(i),S,S1.slice(i),
+               SplusS0_chol.slice(i),
+               S_chol,mu1);
+      mu1mix.col(i) = mu1;
+    }
+  }
+  
+  // Compute the posterior assignment probabilities for the latent
+  // indicator variable.
+  logbfmix += log(w0 + eps);
+  w1 = logbfmix;
+  softmax(w1);
+  
+  // Compute the posterior mean (mu1) and covariance (S1_mix) of the
+  // regression coefficients.
+  S1_mix.fill(0);
+  mu1_mix.fill(0);
+  for (unsigned int i = 0; i < k; i++) {
+    b    = mu1mix.col(i);
+    mu1_mix += w1(i) * b;
+    S1_mix  += w1(i) * (S1.slice(i) + b * trans(b));
+  }
+  S1_mix -= mu1_mix * trans(mu1_mix);
+  
+  // Compute the log-Bayes factor as a linear combination of the
+  // individual Bayes factors for each mixture component.
+  double u = max(logbfmix);
+  return u + log(sum(exp(logbfmix - u)));
+}
+
+// Perform Bayesian multivariate simple regression with summary data
+// and mixture-of-normals prior with centered x.
+double bayes_mvr_mix_centered_X_rss (const vec& xtY, const mat& V,
+                                     const vec& w0, const cube& S0, double xtx, 
+                                     const mat& Vinv, const mat& V_chol,
+                                     const mat& d, const cube& QtimesV_chol,
+                                     double eps, unsigned int nthreads,
+                                     vec& mu1_mix, mat& S1_mix, vec& w1) {
+  unsigned int k = w0.n_elem;
+  unsigned int r = xtY.n_elem;
+  
+  mat  mu1mix(r,k);
+  cube S1mix(r,r,k);
+  vec  logbfmix(k);
+  vec  mu1(r);
+  mat  S1(r,r);
+  
+  // Compute the least-squares estimate.
+  vec b = xtY/xtx;
+  mat S = V/xtx;
+  
+  // Compute quantities needed for bayes_mvr_ridge_centered_X()
+  mat S_chol = V_chol/sqrt(xtx);
+  
+  // Compute the quantities separately for each mixture component.
+  if (nthreads > 1) {
+    bayes_mvr_mix_centered_X_worker worker(b,xtx,V,Vinv,V_chol,S,S_chol,d,
+                                               S0,QtimesV_chol,logbfmix,mu1mix,
+                                               S1mix);
+    parallelFor(0,k,worker);
+  } else {
+    for (unsigned int i = 0; i < k; i++) {
+      logbfmix(i) = bayes_mvr_ridge_centered_X(V,b,S,S0.slice(i),xtx,Vinv,  
+               V_chol,S_chol,d.col(i),
+               QtimesV_chol.slice(i),mu1,S1);
+      mu1mix.col(i)  = mu1;
+      S1mix.slice(i) = S1;
+    }
+  }
+  
+  // Compute the posterior assignment probabilities for the latent
+  // indicator variable.
+  logbfmix += log(w0 + eps);
+  w1 = logbfmix;
+  softmax(w1);
+  
+  // Compute the posterior mean (mu1) and covariance (S1_mix) of the
+  // regression coefficients.
+  S1_mix.fill(0);
+  mu1_mix.fill(0);
+  for (unsigned int i = 0; i < k; i++) {
+    b    = mu1mix.col(i);
+    mu1_mix += w1(i) * b;
+    S1_mix  += w1(i) * (S1mix.slice(i) + b * trans(b));
+  }
+  S1_mix -= mu1_mix * trans(mu1_mix);
+  
+  // Compute the log-Bayes factor as a linear combination of the
+  // individual Bayes factors for each mixture component.
+  double u = max(logbfmix);
+  return u + log(sum(exp(logbfmix - u)));
+}
diff --git a/src/bayes_reg_mv.h b/src/bayes_reg_mv.h
index 41c2c27..0b9f319 100644
--- a/src/bayes_reg_mv.h
+++ b/src/bayes_reg_mv.h
@@ -17,33 +17,62 @@ double bayes_mvr_ridge_centered_X (const arma::mat& V, const arma::vec& b, const
                                    arma::vec& mu1, arma::mat& S1);
 
 double bayes_mvr_mix_standardized_X (const arma::vec& x,
-				     const arma::mat& Y,
-				     const arma::vec& w0,
-				     const arma::cube& S0,
-				     const arma::mat& S,
-				     const arma::cube& S1,
-				     const arma::cube& SplusS0_chol,
-				     const arma::mat& S_chol,
-				     double eps,
-				     unsigned int nthreads,
-				     arma::vec& mu1_mix,
-				     arma::mat& S1_mix,
-				     arma::vec& w1);
+				                             const arma::mat& Y,
+				                             const arma::vec& w0,
+				                             const arma::cube& S0,
+				                             const arma::mat& S,
+				                             const arma::cube& S1,
+				                             const arma::cube& SplusS0_chol,
+				                             const arma::mat& S_chol,
+				                             double eps,
+				                             unsigned int nthreads,
+				                             arma::vec& mu1_mix,
+				                             arma::mat& S1_mix,
+				                             arma::vec& w1);
 
 double bayes_mvr_mix_centered_X (const arma::vec& x,
-				 const arma::mat& Y,
-				 const arma::mat& V,
+				                         const arma::mat& Y,
+				                         const arma::mat& V,
                                  const arma::vec& w0,
-				 const arma::cube& S0,
-				 double xtx, 
+				                         const arma::cube& S0,
+				                         double xtx, 
                                  const arma::mat& Vinv,
-				 const arma::mat& V_chol,
+				                         const arma::mat& V_chol,
                                  const arma::mat& d,
-				 const arma::cube& QtimesV_chol,
-				 double eps,
-				 unsigned int nthreads,
+				                         const arma::cube& QtimesV_chol,
+				                         double eps,
+				                         unsigned int nthreads,
                                  arma::vec& mu1_mix,
-				 arma::mat& S1_mix,
-				 arma::vec& w1);
+				                         arma::mat& S1_mix,
+				                         arma::vec& w1);
+
+double bayes_mvr_mix_standardized_X_rss (unsigned int n,
+                                         const arma::vec& xtY,
+                                         const arma::vec& w0,
+                                         const arma::cube& S0,
+                                         const arma::mat& S,
+                                         const arma::cube& S1,
+                                         const arma::cube& SplusS0_chol,
+                                         const arma::mat& S_chol,
+                                         double eps,
+                                         unsigned int nthreads,
+                                         arma::vec& mu1_mix,
+                                         arma::mat& S1_mix,
+                                         arma::vec& w1);
+
+double bayes_mvr_mix_centered_X_rss (const arma::vec& xtY,
+                                     const arma::mat& V,
+                                     const arma::vec& w0,
+                                     const arma::cube& S0,
+                                     double xtx, 
+                                     const arma::mat& Vinv,
+                                     const arma::mat& V_chol,
+                                     const arma::mat& d,
+                                     const arma::cube& QtimesV_chol,
+                                     double eps,
+                                     unsigned int nthreads,
+                                     arma::vec& mu1_mix,
+                                     arma::mat& S1_mix,
+                                     arma::vec& w1);
 
 #endif
diff --git a/src/misc.cpp b/src/misc.cpp
index c324190..310078f 100644
--- a/src/misc.cpp
+++ b/src/misc.cpp
@@ -87,3 +87,31 @@ void compute_var_part_ERSS (mat& var_part_ERSS, const mat& S1, double xtx){
 //   
 //   return var_part_ERSS;
 // }
+
+// Function to compute some terms of the ELBO with summary data
+void compute_ELBO_rss_terms (double& var_part_tr_wERSS, double& neg_KL, const mat& XtRbar_j,
+                             double logbf, const mat& mu1, const mat& S1, 
+                             double xtx, const mat& Vinv){
+  
+  var_part_tr_wERSS += (as_scalar(accu(Vinv % S1))*xtx);
+  
+  mat mu1XtRbar_j = mu1 * trans(XtRbar_j);
+  
+  neg_KL += (logbf + as_scalar(0.5*accu(Vinv % (- mu1XtRbar_j - trans(mu1XtRbar_j) + mu1 * trans(mu1) * xtx + S1*xtx))));
+}
+
+// // [[Rcpp::plugins("cpp11")]]
+// // [[Rcpp::depends(RcppArmadillo)]]
+// // [[Rcpp::export]]
+// Rcpp::List compute_ELBO_rss_terms_rcpp (double var_part_tr_wERSS_init, double neg_KL_init,
+//                               const arma::mat& xtRbar_j, double logbf, const arma::mat& mu1, const arma::mat& S1,
+//                               double xtx, const arma::mat& Vinv) {
+// 
+//   double var_part_tr_wERSS = var_part_tr_wERSS_init;
+//   double neg_KL = neg_KL_init;
+// 
+//   compute_ELBO_rss_terms(var_part_tr_wERSS, neg_KL, xtRbar_j, logbf, mu1, S1, xtx, Vinv);
+// 
+//   return Rcpp::List::create(Rcpp::Named("var_part_tr_wERSS") = var_part_tr_wERSS,
+//                             Rcpp::Named("neg_KL")            = neg_KL);
+// }
diff --git a/src/misc.h b/src/misc.h
index 2857c77..953df43 100644
--- a/src/misc.h
+++ b/src/misc.h
@@ -31,6 +31,10 @@ void compute_ELBO_terms (double& var_part_tr_wERSS, double& neg_KL, const arma::
                          const arma::mat& rbar_j, double logbf, const arma::mat& mu1, const arma::mat& S1, 
                          double xtx, const arma::mat& Vinv);
 
+void compute_ELBO_rss_terms (double& var_part_tr_wERSS, double& neg_KL,
+                             const arma::mat& XtRbar_j, double logbf, const arma::mat& mu1, const arma::mat& S1, 
+                             double xtx, const arma::mat& Vinv);
+
 void compute_var_part_ERSS (arma::mat& var_part_ERSS, const arma::mat& S1, double xtx);
 
 
diff --git a/src/mr_mash_updates.cpp b/src/mr_mash_updates.cpp
index 711f4fe..f359327 100644
--- a/src/mr_mash_updates.cpp
+++ b/src/mr_mash_updates.cpp
@@ -43,12 +43,20 @@ void inner_loop_general (const mat& X, mat& Rbar, mat& mu1, const mat& V,
                          cube& S1, mat& w1, double& var_part_tr_wERSS, 
                          double& neg_KL, mat& var_part_ERSS);
 
-
 // Missing Y imputation
 void impute_missing_Y (mat& Y, const mat& mu, const mat& Vinv, 
                        const mat& miss, const mat& non_miss,
                        mat& Y_cov, double& sum_neg_ent_Y_miss);
 
+// Inner loop rss
+void inner_loop_general_rss (unsigned int n, const mat& XtX, const mat& XtY, mat& XtRbar, 
+                             mat& mu1, const mat& V, const mat& Vinv, const vec& w0, 
+                             const cube& S0, const mr_mash_precomputed_quantities& precomp_quants,
+                             bool standardize, bool compute_ELBO, bool update_V,
+                             const vec& update_order, double eps, unsigned int nthreads,
+                             cube& S1, mat& w1, double& var_part_tr_wERSS, 
+                             double& neg_KL, mat& var_part_ERSS);
+
 
 // FUNCTION DEFINITIONS
 // --------------------
@@ -127,7 +135,7 @@ void inner_loop_general (const mat& X, mat& Rbar, mat& mu1, const mat& V,
     x = X.col(j);
     mu1_j = trans(mu1.row(j));
     
-    // Disregard the ith predictor in the expected residuals.
+    // Disregard the jth predictor in the expected residuals.
     Rbar_j = Rbar + (x * trans(mu1_j));
     
     // Update the posterior quantities for the jth
@@ -175,11 +183,6 @@ void inner_loop_general (const mat& X, mat& Rbar, mat& mu1, const mat& V,
   }
 }
 
-
-
-
-
-
 // Impute missing Y
 //
 // [[Rcpp::depends(RcppArmadillo)]]
@@ -241,3 +244,116 @@ void impute_missing_Y (mat& Y, const mat& mu, const mat& Vinv,
     }
   }
 }
+
+// Inner loop rss
+//
+// [[Rcpp::depends(RcppArmadillo)]]
+// [[Rcpp::depends(RcppParallel)]]
+// [[Rcpp::export]]
+List inner_loop_general_rss_rcpp (unsigned int n, const arma::mat& XtX, const arma::mat& XtY, 
+                                  arma::mat& XtRbar, arma::mat& mu1, const arma::mat& V, 
+                                  const arma::mat& Vinv, const arma::vec& w0,
+                                  const arma::cube& S0, const List& precomp_quants_list,
+                                  bool standardize, bool compute_ELBO, bool update_V,
+                                  const arma::vec& update_order, double eps, unsigned int nthreads) {
+  unsigned int r = mu1.n_cols;
+  unsigned int p = mu1.n_rows;
+  unsigned int k = w0.n_elem;
+  cube S1(r,r,p);
+  mat  w1(p,k);
+  mat  mu1_new  = mu1;
+  double var_part_tr_wERSS;
+  double neg_KL;
+  mat var_part_ERSS(r,r);
+  mr_mash_precomputed_quantities precomp_quants
+  (as(precomp_quants_list["S"]),
+   as(precomp_quants_list["S_chol"]),
+   as(precomp_quants_list["S1"]),
+   as(precomp_quants_list["SplusS0_chol"]),
+   as(precomp_quants_list["V_chol"]),
+   as(precomp_quants_list["d"]),
+   as(precomp_quants_list["QtimesV_chol"]),
+   as(precomp_quants_list["xtx"]));
+  inner_loop_general_rss(n, XtX, XtY, XtRbar, mu1_new, V, Vinv, w0, S0, precomp_quants,
+                         standardize, compute_ELBO, update_V, update_order, eps,
+                         nthreads, S1, w1, var_part_tr_wERSS, neg_KL, var_part_ERSS);
+  return List::create(Named("mu1")                = mu1_new,
+                      Named("S1")                 = S1,
+                      Named("w1")                 = w1,
+                      Named("var_part_tr_wERSS")  = var_part_tr_wERSS,
+                      Named("neg_KL")             = neg_KL,
+                      Named("var_part_ERSS")      = var_part_ERSS);
+}
+
+// Perform the inner loop rss
+void inner_loop_general_rss (unsigned int n, const mat& XtX, const mat& XtY, mat& XtRbar, 
+                             mat& mu1, const mat& V, const mat& Vinv, const vec& w0, 
+                             const cube& S0, const mr_mash_precomputed_quantities& precomp_quants,
+                             bool standardize, bool compute_ELBO, bool update_V, 
+                             const vec& update_order, double eps, unsigned int nthreads,
+                             cube& S1, mat& w1, double& var_part_tr_wERSS, 
+                             double& neg_KL, mat& var_part_ERSS) {
+  unsigned int p = mu1.n_rows;
+  unsigned int r = mu1.n_cols;
+  unsigned int k = w0.n_elem;
+  vec Xtx(p);
+  vec XtRbar_j(r);
+  vec mu1_j(r);
+  vec mu1_mix(r);
+  mat S1_mix(r,r);
+  vec w1_mix(k);
+  double logbf_mix;
+  double xtx_j;
+  
+  // Initialize ELBO parameters
+  var_part_tr_wERSS = 0;
+  neg_KL = 0;
+  
+  // Initialize V parameters
+  var_part_ERSS.zeros(r,r);
+  
+  // Repeat for each predictor.
+  for (unsigned int j : update_order) {
+    
+    if (standardize)
+      xtx_j = n-1;
+    else
+      xtx_j = precomp_quants.xtx(j);
+    
+    Xtx = XtX.col(j);
+    mu1_j = trans(mu1.row(j));
+    
+    // Disregard the jth predictor in the expected residuals.
+    XtRbar += (Xtx * trans(mu1_j));
+    XtRbar_j = trans(XtRbar.row(j));
+    
+    // Update the posterior quantities for the jth
+    // predictor.
+    if (standardize)
+      logbf_mix = bayes_mvr_mix_standardized_X_rss(n, XtRbar_j, w0, S0, precomp_quants.S,
+                                                   precomp_quants.S1,
+                                                   precomp_quants.SplusS0_chol,
+                                                   precomp_quants.S_chol, eps, nthreads,
+                                                   mu1_mix, S1_mix, w1_mix);
+    else
+      logbf_mix = bayes_mvr_mix_centered_X_rss(XtRbar_j, V, w0, S0, xtx_j, Vinv,
+                                               precomp_quants.V_chol, precomp_quants.d, 
+                                               precomp_quants.QtimesV_chol, eps, nthreads,
+                                               mu1_mix, S1_mix, w1_mix);
+    
+    mu1.row(j)  = trans(mu1_mix);
+    S1.slice(j) = S1_mix;
+    w1.row(j)   = trans(w1_mix);
+    
+    // Compute ELBO parameters
+    if (compute_ELBO)
+      compute_ELBO_rss_terms(var_part_tr_wERSS, neg_KL, XtRbar_j, logbf_mix, mu1_mix, S1_mix, xtx_j, Vinv);
+    
+    // Compute V parameters
+    if (update_V)
+      compute_var_part_ERSS(var_part_ERSS, S1_mix, xtx_j);
+    
+    // Update the expected residuals.
+    XtRbar -= (Xtx * trans(mu1_mix));
+  }
+}
diff --git a/src/update_order.cpp b/src/update_order.cpp
index d01b6e4..41bca8d 100644
--- a/src/update_order.cpp
+++ b/src/update_order.cpp
@@ -37,6 +37,12 @@ vec compute_logbf (const mat& X, const mat& Y, const mat& V,
                    const mr_mash_precomputed_quantities& precomp_quants,
                    bool standardize, double eps, unsigned int nthreads);
 
+vec compute_logbf_rss (unsigned int n, const mat& XtY, const mat& V,
+                       const mat& Vinv, const vec& w0, const cube& S0,
+                       const mr_mash_precomputed_quantities& precomp_quants,
+                       bool standardize, double eps, unsigned int nthreads);
+
+
 // FUNCTION DEFINITIONS
 // --------------------
 // Compute logbf for each variable
@@ -61,6 +67,7 @@ arma::vec compute_logbf_rcpp (const arma::mat& X, const arma::mat& Y, const arma
   return compute_logbf(X, Y, V, Vinv, w0, S0, precomp_quants, standardize, eps, nthreads);
 }
 
+
 vec compute_logbf (const mat& X, const mat& Y, const mat& V,
                    const mat& Vinv, const vec& w0, const cube& S0,
                    const mr_mash_precomputed_quantities& precomp_quants,
@@ -101,3 +108,69 @@ vec compute_logbf (const mat& X, const mat& Y, const mat& V,
   
   return logbf_mix_all;
 }
+
+
+// Compute logbf for each variable from summary data only
+// [[Rcpp::depends(RcppArmadillo)]]
+// [[Rcpp::export]]
+arma::vec compute_logbf_rss_rcpp (unsigned int n, const arma::mat& XtY, const arma::mat& V,
+                                  const arma::mat& Vinv, const arma::vec& w0, const arma::cube& S0,
+                                  const List& precomp_quants_list,
+                                  bool standardize, double eps, unsigned int nthreads) {
+  
+  mr_mash_precomputed_quantities precomp_quants
+  (as(precomp_quants_list["S"]),
+   as(precomp_quants_list["S_chol"]),
+   as(precomp_quants_list["S1"]),
+   as(precomp_quants_list["SplusS0_chol"]),
+   as(precomp_quants_list["V_chol"]),
+   as(precomp_quants_list["d"]),
+   as(precomp_quants_list["QtimesV_chol"]),
+   as(precomp_quants_list["xtx"]));
+  
+  
+  return compute_logbf_rss(n, XtY, V, Vinv, w0, S0, precomp_quants, standardize, eps, nthreads);
+}
+
+
+vec compute_logbf_rss (unsigned int n, const mat& XtY, const mat& V,
+                       const mat& Vinv, const vec& w0, const cube& S0,
+                       const mr_mash_precomputed_quantities& precomp_quants,
+                       bool standardize, double eps, unsigned int nthreads){
+  
+  unsigned int p = XtY.n_rows;
+  unsigned int r = XtY.n_cols;
+  unsigned int k = S0.n_slices;
+  vec XtY_j(r);
+  vec mu1_mix(r);
+  mat S1_mix(r,r);
+  vec w1_mix(k);
+  double logbf_mix;
+  vec logbf_mix_all(p);
+  
+  // Repeat for each predictor.
+  for (unsigned int j = 0; j < p; j++) {
+    
+    XtY_j = trans(XtY.row(j));
+
+    // Update the posterior quantities for the jth
+    // predictor.
+    if (standardize)
+      logbf_mix = bayes_mvr_mix_standardized_X_rss(n, XtY_j, w0, S0, precomp_quants.S,
+                                               precomp_quants.S1,
+                                               precomp_quants.SplusS0_chol,
+                                               precomp_quants.S_chol, eps, nthreads,
+                                               mu1_mix, S1_mix, w1_mix);
+    else {
+      double xtx_j = precomp_quants.xtx(j);
+      logbf_mix = bayes_mvr_mix_centered_X_rss(XtY_j, V, w0, S0, xtx_j, Vinv,
+                                               precomp_quants.V_chol, precomp_quants.d, 
+                                               precomp_quants.QtimesV_chol, eps, nthreads,
+                                               mu1_mix, S1_mix, w1_mix);
+    }
+    
+    logbf_mix_all(j) = logbf_mix;
+  }
+  
+  return logbf_mix_all;
+}
diff --git a/tests/testthat/test_compute_rank_variables_BFmix_R_vs_Rcpp.R b/tests/testthat/test_compute_rank_variables_BFmix_R_vs_Rcpp.R
index e34ecc7..1b75f18 100644
--- a/tests/testthat/test_compute_rank_variables_BFmix_R_vs_Rcpp.R
+++ b/tests/testthat/test_compute_rank_variables_BFmix_R_vs_Rcpp.R
@@ -34,26 +34,44 @@ test_that("R and Rcpp version of compute_rank_variables_BFmix return the same re
   ###Compute the inverse of V
   Vinv <- solve(V)
   
+  ###Compute quantities needed for rss version
+  XtY <- crossprod(X, scale(Y, center=TRUE, scale=FALSE))
+  
   ###Set eps
   eps <- .Machine$double.eps
   
   ###Compute logbf with standardize=TRUE
-  comps_r <- precompute_quants(X, V, S0mix, standardize=TRUE, version="R")
+  comps_r <- precompute_quants(n, V, S0mix, standardize=TRUE, version="R")
   ranks_r <- compute_rank_variables_BFmix(X, Y, V, Vinv, w0, S0mix, comps_r, standardize=TRUE, version="R", decreasing=TRUE, eps)
   
-  comps_rcpp <- precompute_quants(X, V, S0mix, standardize=TRUE, version="Rcpp")
+  comps_rcpp <- precompute_quants(n, V, S0mix, standardize=TRUE, version="Rcpp")
   ranks_rcpp <- compute_rank_variables_BFmix(X, Y, V, Vinv, w0, S0mix, comps_rcpp, standardize=TRUE, version="Rcpp", decreasing=TRUE, eps, nthreads=1)
+
+  ranks_rss_r <- compute_rank_variables_BFmix_rss(n, XtY, V, Vinv, w0, S0mix, comps_r, standardize=TRUE, version="R", decreasing=TRUE, eps)
+  
+  ranks_rss_rcpp <- compute_rank_variables_BFmix_rss(n, XtY, V, Vinv, w0, S0mix, comps_rcpp, standardize=TRUE, version="Rcpp", decreasing=TRUE, eps, nthreads=1)
   
+    
   ###Compute logbf with standardize=FALSE
-  comps1_r <- precompute_quants(X, V, S0mix, standardize=FALSE, version="R")
+  comps1_r <- precompute_quants(n, V, S0mix, standardize=FALSE, version="R")
   comps1_r$xtx <- colSums(X^2)
   ranks1_r <- compute_rank_variables_BFmix(X, Y, V, Vinv, w0, S0mix, comps1_r, standardize=FALSE, version="R", decreasing=TRUE, eps)
   
-  comps1_rcpp <- precompute_quants(X, V, S0mix, standardize=FALSE, version="Rcpp")
+  comps1_rcpp <- precompute_quants(n, V, S0mix, standardize=FALSE, version="Rcpp")
   comps1_rcpp$xtx <- colSums(X^2)
   ranks1_rcpp <- compute_rank_variables_BFmix(X, Y, V, Vinv, w0, S0mix, comps1_rcpp, standardize=FALSE, version="Rcpp", decreasing=TRUE, eps, nthreads=1)
+ 
+  ranks1_rss_r <- compute_rank_variables_BFmix_rss(n, XtY, V, Vinv, w0, S0mix, comps1_r, standardize=FALSE, version="R", decreasing=TRUE, eps)
   
+  ranks1_rss_rcpp <- compute_rank_variables_BFmix_rss(n, XtY, V, Vinv, w0, S0mix, comps1_rcpp, standardize=FALSE, version="Rcpp", decreasing=TRUE, eps, nthreads=1)
+  
+   
   ###Tests
   expect_equal(ranks_r, ranks_rcpp, tolerance = 1e-10, scale = 1)
   expect_equal(ranks1_r, ranks1_rcpp, tolerance = 1e-10, scale = 1)
+  expect_equal(ranks_rss_r, ranks_rss_rcpp, tolerance = 1e-10, scale = 1)
+  expect_equal(ranks1_rss_r, ranks1_rss_rcpp, tolerance = 1e-10, scale = 1)
+  expect_equal(ranks_rss_rcpp, ranks_rcpp, tolerance = 1e-10, scale = 1)
+  expect_equal(ranks1_rss_rcpp, ranks1_rcpp, tolerance = 1e-10, scale = 1)
+  
 })
\ No newline at end of file
diff --git a/tests/testthat/test_mr.mash_multithreading.R b/tests/testthat/test_mr.mash_multithreading.R
index 7f927b0..4a2bbfa 100644
--- a/tests/testthat/test_mr.mash_multithreading.R
+++ b/tests/testthat/test_mr.mash_multithreading.R
@@ -42,6 +42,12 @@ test_that("mr.mash with 1 or 2 thread(s) return the same results", {
   ###Estimate residual covariance
   V_est <- cov(Y)
   
+  ###Compute quantities needed for mr.mash.rss
+  out <- compute_univariate_sumstats(X=X, Y=Y, standardize=FALSE, standardize.response=FALSE, mc.cores=1)
+  R <- cor(X)
+  X_colMeans <- colMeans(X)
+  Y_colMeans <- colMeans(Y)
+  
   ###Fit with current implementation (1 thread)
   capture.output(
     fit_1 <- mr.mash(X, Y, S0mix, w0, V_est, update_w0=TRUE,
@@ -56,6 +62,14 @@ test_that("mr.mash with 1 or 2 thread(s) return the same results", {
   fit_1_miss$progress <- fit_1_miss$progress[, -2] ##This line is needed to remove the timing column -->  
   ##hopefully faster when using multiple threads
   
+  fit_1_rss <- mr.mash.rss(Bhat=out$Bhat, Shat=out$Shat, covY=V_est, R=R, n=n, S0=S0mix, 
+                           w0=w0, V=V_est, update_w0=TRUE, compute_ELBO=TRUE, standardize=TRUE,
+                           verbose=FALSE, update_V=TRUE, X_colmeans=X_colMeans, Y_colmeans=Y_colMeans, 
+                           nthreads=1)
+  fit_1_rss$progress <- fit_1_rss$progress[, -2] ##This line is needed to remove the timing column -->  
+  ##hopefully faster when using multiple threads
+
+  
   
   ###Fit with current implementation (2 threads)
   capture.output(
@@ -70,6 +84,14 @@ test_that("mr.mash with 1 or 2 thread(s) return the same results", {
                           verbose=FALSE, update_V=TRUE, nthreads=2))
   fit_2_miss$progress <- fit_2_miss$progress[, -2]
   
+  fit_2_rss <- mr.mash.rss(Bhat=out$Bhat, Shat=out$Shat, covY=V_est, R=R, n=n, S0=S0mix, 
+                           w0=w0, V=V_est, update_w0=TRUE, compute_ELBO=TRUE, standardize=TRUE,
+                           verbose=FALSE, update_V=TRUE, X_colmeans=X_colMeans, Y_colmeans=Y_colMeans, 
+                           nthreads=2)
+  fit_2_rss$progress <- fit_2_rss$progress[, -2] ##This line is needed to remove the timing column -->  
+  ##hopefully faster when using multiple threads
+  
+  
   
   ###Test
   expect_equal(fit_1, fit_2, tolerance=1e-10, scale=1)
diff --git a/tests/testthat/test_mr.mash_vs_mr.mash.rss.R b/tests/testthat/test_mr.mash_vs_mr.mash.rss.R
new file mode 100644
index 0000000..30aa9c1
--- /dev/null
+++ b/tests/testthat/test_mr.mash_vs_mr.mash.rss.R
@@ -0,0 +1,128 @@
+context("mr.mash and mr.mash.rss versions return same result")
+
+test_that("mr.mash and mr.mash.rss return the same results", {
+  ###Set seed
+  set.seed(123)
+  
+  ###Simulate X and Y
+  n  <- 100
+  p <- 10
+  
+  ###Set residual covariance
+  V  <- rbind(c(1.0,0.2),
+              c(0.2,0.4))
+  
+  ###Set true effects
+  B  <- matrix(c(-2, -2,
+                 5, 5,
+                 rep(0, (p-2)*2)), byrow=TRUE, ncol=2)
+  
+  ###Simulate X
+  X <- matrix(rnorm(n*p), nrow=n, ncol=p)
+  X <- scale(X, center=TRUE, scale=FALSE)
+  
+  ###Simulate Y from MN(XB, I_n, V) where I_n is an nxn identity
+  ###matrix and V is the residual covariance
+  Y <- sim_mvr(X, B, V)
+  
+  ###Specify the mixture weights and covariance matrices for the
+  ###mixture-of-normals prior
+  grid <- seq(1, 5)
+  S0mix <- compute_cov_canonical(ncol(Y), singletons=TRUE,
+                                 hetgrid=c(0, 0.25, 0.5, 0.75, 0.99), grid,
+                                 zeromat=TRUE)
+  
+  w0    <- rep(1/(length(S0mix)), length(S0mix))
+  
+  ###Estimate residual covariance
+  V_est <- cov(Y)
+  
+  ###Fit mr.mash
+  fit <- mr.mash(X, Y, S0mix, w0, V_est, update_w0=TRUE,
+                   update_w0_method="EM", compute_ELBO=TRUE, standardize=FALSE,
+                   verbose=FALSE, update_V=FALSE)
+  fit$progress <- fit$fitted <- fit$pve <- fit$G <- NULL
+  
+  fit_scaled <- mr.mash(X, Y, S0mix, w0, V_est, update_w0=TRUE,
+                          update_w0_method="EM", compute_ELBO=TRUE, 
+                          standardize=TRUE, verbose=FALSE, update_V=FALSE)
+  fit_scaled$progress <- fit_scaled$fitted <- fit_scaled$pve <- fit_scaled$G <- NULL
+  
+  fit_V <- mr.mash(X, Y, S0mix, w0, V_est, update_w0=TRUE,
+                     update_w0_method="EM", compute_ELBO=TRUE, 
+                     standardize=FALSE, verbose=FALSE, update_V=TRUE)
+  fit_V$progress <- fit_V$fitted <- fit_V$pve <- fit_V$G <- NULL
+  
+  fit_V_diag <- mr.mash(X, Y, S0mix, w0, V_est, update_w0=TRUE,
+                   update_w0_method="EM", compute_ELBO=TRUE, 
+                   standardize=FALSE, verbose=FALSE, update_V=TRUE,
+                   update_V_method="diagonal")
+  fit_V_diag$progress <- fit_V_diag$fitted <- fit_V_diag$pve <- fit_V_diag$G <- NULL
+  
+  fit_scaled_V <- mr.mash(X, Y, S0mix, w0, V_est, update_w0=TRUE,
+                            update_w0_method="EM", compute_ELBO=TRUE, 
+                            standardize=TRUE, verbose=FALSE, update_V=TRUE)
+  fit_scaled_V$progress <- fit_scaled_V$fitted <- fit_scaled_V$pve <- fit_scaled_V$G <- NULL
+  
+  fit_scaled_V_declogBF <- mr.mash(X, Y, S0mix, w0, V_est, update_w0=TRUE,
+                            update_w0_method="EM", compute_ELBO=TRUE, 
+                            standardize=TRUE, verbose=FALSE, update_V=TRUE,
+                            ca_update_order="decreasing_logBF")
+  fit_scaled_V_declogBF$progress <- fit_scaled_V_declogBF$fitted <- fit_scaled_V_declogBF$pve <- fit_scaled_V_declogBF$G <- NULL
+
+
+  ###Fit mr.mash.rss
+  out <- compute_univariate_sumstats(X=X, Y=Y, standardize=FALSE, standardize.response=FALSE, mc.cores=1)
+  R <- cor(X)
+  X_colMeans <- colMeans(X)
+  Y_colMeans <- colMeans(Y)
+  
+  fit_rss <- mr.mash.rss(Bhat=out$Bhat, Shat=out$Shat, covY=V_est, R=R, n=n, S0=S0mix, 
+                          w0=w0, V=V_est, update_w0=TRUE, compute_ELBO=TRUE, standardize=FALSE,
+                          verbose=FALSE, update_V=FALSE, X_colmeans=X_colMeans, Y_colmeans=Y_colMeans)
+  fit_rss$progress <- NULL
+  
+  fit_scaled_rss <- mr.mash.rss(Bhat=out$Bhat, Shat=out$Shat, covY=V_est, R=R, n=n, S0=S0mix, 
+                                w0=w0, V=V_est, update_w0=TRUE, compute_ELBO=TRUE, 
+                                standardize=TRUE, verbose=FALSE, update_V=FALSE,
+                                X_colmeans=X_colMeans, Y_colmeans=Y_colMeans)
+  fit_scaled_rss$progress <- NULL
+  
+  fit_V_rss <- mr.mash.rss(Bhat=out$Bhat, Shat=out$Shat, covY=V_est, R=R, n=n, S0=S0mix, 
+                          w0=w0, V=V_est, update_w0=TRUE, compute_ELBO=TRUE, 
+                          standardize=FALSE, verbose=FALSE, update_V=TRUE,
+                          X_colmeans=X_colMeans, Y_colmeans=Y_colMeans)
+  fit_V_rss$progress <- NULL
+  
+  fit_V_diag_rss <- mr.mash.rss(Bhat=out$Bhat, Shat=out$Shat, covY=V_est, R=R, n=n, S0=S0mix, 
+                           w0=w0, V=V_est, update_w0=TRUE, compute_ELBO=TRUE, 
+                           standardize=FALSE, verbose=FALSE, update_V=TRUE,
+                           X_colmeans=X_colMeans, Y_colmeans=Y_colMeans,
+                           update_V_method="diagonal")
+  fit_V_diag_rss$progress <- NULL
+  
+  fit_scaled_V_rss <- mr.mash.rss(Bhat=out$Bhat, Shat=out$Shat, covY=V_est, R=R, n=n, S0=S0mix, 
+                                  w0=w0, V=V_est, update_w0=TRUE, compute_ELBO=TRUE, 
+                                  standardize=TRUE, verbose=FALSE, update_V=TRUE,
+                                  X_colmeans=X_colMeans, Y_colmeans=Y_colMeans)
+  fit_scaled_V_rss$progress <- NULL
+  
+  
+  fit_scaled_V_declogBF_rss <- mr.mash.rss(Bhat=out$Bhat, Shat=out$Shat, covY=V_est, R=R, n=n, S0=S0mix, 
+                                          w0=w0, V=V_est, update_w0=TRUE, compute_ELBO=TRUE, 
+                                          standardize=TRUE, verbose=FALSE, update_V=TRUE,
+                                          ca_update_order="decreasing_logBF",
+                                          X_colmeans=X_colMeans, Y_colmeans=Y_colMeans)
+  fit_scaled_V_declogBF_rss$progress <- NULL
+  
+  
+  
+  
+  ###Tests
+  expect_equal(unclass(fit), unclass(fit_rss), tolerance=1e-10, scale=1)
+  expect_equal(unclass(fit_scaled), unclass(fit_scaled_rss), tolerance=1e-10, scale=1)
+  expect_equal(unclass(fit_V), unclass(fit_V_rss), tolerance=1e-10, scale=1)
+  expect_equal(unclass(fit_V_diag), unclass(fit_V_diag_rss), tolerance=1e-10, scale=1)
+  expect_equal(unclass(fit_scaled_V), unclass(fit_scaled_V_rss), tolerance=1e-10, scale=1)
+  expect_equal(unclass(fit_scaled_V_declogBF), unclass(fit_scaled_V_declogBF_rss), tolerance=1e-10, scale=1)
+})
diff --git a/vignettes/mr_mash_intro.Rmd b/vignettes/mr_mash_intro.Rmd
index f91abaa..761bfdf 100644
--- a/vignettes/mr_mash_intro.Rmd
+++ b/vignettes/mr_mash_intro.Rmd
@@ -28,7 +28,7 @@ and we load the "mr.mash.alpha" package.
 
 ```{r load-pkgs}
 library(mr.mash.alpha)
-set.seed(123)
+set.seed(12)
 ```
 
 We illustrate the application of mr.mash to a data set simulated from