Fit mash model and estimate residual correlations using EM algorithm

Estimates a residual correlation matrix from data using an ad hoc EM algorithm.

mash_estimate_corr_em(
  data,
  Ulist,
  init,
  max_iter = 30,
  tol = 1,
  est_cor = TRUE,
  track_fit = FALSE,
  prior = c("nullbiased", "uniform"),
  details = TRUE,
  ...
)

Arguments

data: a mash data object, eg as created by mash_set_data
Ulist: a list of covariance matrices to use
init: the initial value for the residual correlation. If it is not given, we use result from estimate_null_correlation_simple
max_iter: maximum number of iterations to perform
tol: convergence tolerance
est_cor: whether to estimate correlation matrix (TRUE) or the covariance matrix (FALSE)
track_fit: add an attribute trace to output that saves current values of all iterations
prior: indicates what penalty to use on the likelihood, if any
details: whether to return details of the model, if it is TRUE, the mash model, the number of iterations and the value of objective functions will be returned
...: other parameters pass to mash

Value

the estimated correlation matrix and the fitted mash model

V: estimated residual correlation matrix
mash.model: fitted mash model

Details

Returns the estimated residual correlation matrix among conditions. We estimate the residual correlation matrix using an ad hoc em algorithm. The update in the ad hoc M step is not guaranteed to increase the likelihood, therefore, the EM algorithm is stopped before the likelihood drops. The residual correlation matrix V is estimated using the posterior second moment of the noise.

Warning: This method could take some time. The estimate_null_correlation_simple gives a quick approximation for the null correlation matrix.

Examples

simdata = simple_sims(100,5,1)
m.1by1 = mash_1by1(mash_set_data(simdata$Bhat,simdata$Shat))
strong.subset = get_significant_results(m.1by1,0.05)
random.subset = sample(1:nrow(simdata$Bhat),20)
data.strong = mash_set_data(simdata$Bhat[strong.subset,], simdata$Shat[strong.subset,])
data.tmp = mash_set_data(simdata$Bhat[random.subset,], simdata$Shat[random.subset,])
U_pca = cov_pca(data.strong, 3)
U_ed = cov_ed(data.strong, U_pca)
Vhat = mash_estimate_corr_em(data.tmp, U_ed)