This is prototype to a unit test to verify implementation of mash computation is correct, by comparing it to univariate case when $Y$ has one column and prior covariance matrices is fixed.
set.seed(1)
n = 1000
p = 1000
beta = rep(0,p)
beta[1:4] = 1
X = matrix(rnorm(n*p),nrow=n,ncol=p)
y = X %*% beta + rnorm(n)
#' res =susie(X,y,L=10)
library(mvsusieR)
prior_var = 0.2 * as.numeric(var(y))
residual_var = as.numeric(var(y))
data = mvsusieR:::DenseData$new(X,y)
residual_var
prior_var
m1 = mvsusieR:::BayesianSimpleRegression$new(ncol(X), residual_var, prior_var)
m1$fit(data, save_summary_stats = T)
head(m1$posterior_b1)
m1
# Assuming 1 out of $J$ are causal, we place a null weight $1-1/J$ a priori.
# This will lead to some shrinkage
# null_weight = 1 - 1 / ncol(X)
null_weight = 0
prior_covar = mvsusieR:::MashInitializer$new(list(0.2 * cov(y)), 1, prior_weights=1 - null_weight, null_weight=null_weight, alpha = 0)
residual_covar = cov(y)
prior_covar$prior_variance
residual_covar
m2 = mvsusieR:::MashRegression$new(ncol(X), residual_covar, prior_covar)
m2$fit(data, save_summary_stats = T)
head(m2$posterior_b1)
m2
All quantities seem to agree now.
ASH works well.
library(ashr)
a.out = ash(as.vector(m1$bhat), as.vector(m1$sbhat), mixcompdist = 'normal')
head(get_pm(a.out))
prior_covar$mash_prior
library(mashr)
data = mash_set_data(m2$bhat, m2$sbhat)
m.c = mash(data, g = prior_covar$mash_prior, fixg = TRUE, algorithm ='Rcpp')
head(get_pm(m.c))
Very similiar results to what I got with fixed g earlier.
U.c = cov_canonical(data)
print(names(U.c))
m.c = mash(data, U.c, , algorithm ='Rcpp')
head(get_pm(m.c))
m.c$fitted_g