<jupyter_start><jupyter_text>Decreasing ELBO issuePossibly caused by comparing and setting prior variances to zero when estimating them with EM.  ContextIn mvSuSiE model with mixture prior, to estimate multivariate prior $b \sim \sum_k \pi_k MVN(0, \sigma_0^2U_k)$ we fix $U_k$ and $\pi_k$ and only estimate the scalar $\sigma_0^2$.Recent updates to `mvsusieR` package implemented EM updates for $\sigma_0^2$ in both MASH regression and Bayesian multivariate regression (hereafter BMR, essentially MASH with just one component but we have separate, simpler code for it written in R for prototyping purpose). I verified with some tests that for IBSS using MASH mixture with only one matrix the result agree with IBSS using BMR. And BMR, when reduced to just one response, agrees with `susieR::susie` result. The EM estimate for prior variance scalar are identical and ELBO are non-decreasing. ProblemHowever, as I moved on to more simulations I get warnings for decreasing ELBO with the EM estimates. It seems to happen only when there are many effect variables relative to conditions. For example my toy example below has 2 conditions each having 2 effect variables out of 50 variables. The problem is that EM method, as it is currently implemented, estimates $\sigma_0^2$ then at the end of an iteration, compares the update with zero and accept whichever has better likelihood for that single effect model. This resulted in decreasing ELBO. So far we haven't observed the behavior in SuSiE -- no matter how many effect variables are there, the ELBO is non-decreasing using EM updates. I think both Yuxin and Kaiqian has done some benchmark with SuSiE's EM for prior but don't see an issue.  Possible causeIf I switch off the comparison between current estimate and zero at the end of each EM update ([comment out this line](https://github.com/stephenslab/mvsusieR/blob/390eb23852fa5ed5c17b5062deb6f3e8494ba76e/R/single_effect_model.RL22)), then in this example the algorithm takes longer to converge but ELBO is non-decreasing, and is eventually larger than any other approaches (eg, `simple` and `optim` methods). The estimates make more sense compared to the "truth" we know for this simulated data.I think since EM is not directly maximizing that single effect loglik, we should not compare it with zero at each iteration? Otherwise comparing it and setting $\sigma_0^2 = 0$ will provide a very bad initialization for upcoming iterations.To see the practical impact: a major difference between EM update and maximizing the single effect loglik directly is that the former is a function of posterior quantities whereas the latter only involves the data, with the prior to be estimated ([here for instance](https://github.com/stephenslab/mvsusieR/blob/390eb23852fa5ed5c17b5062deb6f3e8494ba76e/R/bayesian_simple_regression.RL73)). As a result, in EM updates if in a previous iteration sets the estimate to zero then it will impact the posterior of the next iteration and the estimate will never come back from zero. This is not the case for when prior is directly estimated as in `optim`. IllustrationHere I use BMR not MASH regression to illustrate. **To reproduce you need to run on R 3.6.1+**.<jupyter_code>devtools::load_all('~/GIT/software/mvsusieR')
set.seed(1)
dat = mvsusie_sim1(n=50,p=50,r=2,s=2)
L = 10<jupyter_output>Loading mvsusieR

Loading required package: mashr

Loading required package: ashr<jupyter_text>There are 2 conditions, each having 2 effect variables.<jupyter_code>sum(rowSums(dat$b)>0)<jupyter_output><empty_output><jupyter_text>Fit without estimating priorWe use fixed prior here provided by simulated data object; everything looks fine:<jupyter_code>dat$V
fit0 = mvsusie(dat$X,dat$y,L=L,prior_variance=dat$V, 
              compute_objective=T, estimate_residual_variance=F, 
              estimate_prior_variance=F)
fit0$elbo
fit0$V
susieR::susie_plot(fit0, 'PIP', b = rowSums(dat$b))<jupyter_output><empty_output><jupyter_text>`simple` method Using `simple` method, only 1 out of 4 effects was captured.<jupyter_code>fit1 = mvsusie(dat$X,dat$y,L=L,prior_variance=dat$V, 
              compute_objective=T, estimate_residual_variance=F, 
              estimate_prior_variance=T, estimate_prior_method='simple',
              track_fit=TRUE)
fit1$elbo
fit1$V<jupyter_output><empty_output><jupyter_text>Looking at prior estimates in each iteration,<jupyter_code>fit1$prior_history<jupyter_output><empty_output><jupyter_text>`EM` methodNow we can see decreasing ELBO: from -183 to -185, compared to previously the ELBO was -182 at convergence,<jupyter_code>fit2 = mvsusie(dat$X,dat$y,L=L,prior_variance=dat$V, 
              compute_objective=T, estimate_residual_variance=F, 
              estimate_prior_variance=T, estimate_prior_method = 'EM',
              track_fit=TRUE)
fit2$elbo
fit2$V<jupyter_output><empty_output><jupyter_text>And no effect was captured! It seems in the first iteration both prior are estimated non-zero. But the 2nd iteration they both became zero because zero maximizes the loglik in that iteration. However apparently in the next iteration setting prior to zero decreased the ELBO:<jupyter_code>fit2$prior_history<jupyter_output><empty_output><jupyter_text>However, when I set `L = 5` it can capture 1 effect,<jupyter_code>fit21 = mvsusie(dat$X,dat$y,L=4,prior_variance=dat$V, 
              compute_objective=T, estimate_residual_variance=F, 
              estimate_prior_variance=T, estimate_prior_method = 'EM',
              track_fit=T)
fit21$elbo
fit21$V
fit21$prior_history
susieR::susie_plot(fit21, 'PIP', b = rowSums(dat$b))
sum(fit21$pip)<jupyter_output><empty_output><jupyter_text>`optim` methodIn `optim` method we still have that check at the end of each update comparing current estimate with zero in terms of likelihood. Thus many `V` are set to zero. However, two effects were captured this time,<jupyter_code>fit3 = mvsusie(dat$X,dat$y,L=L,prior_variance=dat$V, 
              compute_objective=T, estimate_residual_variance=F, 
              estimate_prior_variance=T, estimate_prior_method='optim',
              track_fit=T)
fit3$elbo
fit3$V<jupyter_output><empty_output><jupyter_text>It seems unlike with the EM update, the estimates never hit zero:<jupyter_code>fit3$prior_history
susieR::susie_plot(fit3, 'PIP', b = rowSums(dat$b))<jupyter_output><empty_output><jupyter_text>`EM` update without checking with zeroI commented out [this line](https://github.com/stephenslab/mvsusieR/blob/390eb23852fa5ed5c17b5062deb6f3e8494ba76e/R/single_effect_model.RL22)) in my code, and run EM with `L=10`:<jupyter_code>devtools::load_all('~/GIT/software/mvsusieR')
fit4 = mvsusie(dat$X,dat$y,L=L,prior_variance=dat$V, 
              compute_objective=T, estimate_residual_variance=F, 
              estimate_prior_variance=T, estimate_prior_method='EM', max_iter = 200,
              track_fit=T)
fit4$elbo
fit4$niter<jupyter_output><empty_output><jupyter_text>It took 105 iterations to converge. The ELBO -181 is comparable to that obtained by `optim`.The prior estimates are:<jupyter_code>fit4$V
fit4$prior_history<jupyter_output><empty_output><jupyter_text>Here it captures 2 effects, as also reflected below that two effects have large PIP,<jupyter_code>susieR::susie_plot(fit4, 'PIP', b = rowSums(dat$b))<jupyter_output><empty_output><jupyter_text>Final solutionFollowing from [Matthew's suggestion](https://github.com/stephenslab/mvsusieR/issues/26issuecomment-612947198), we check with zero and set all posteriors to zero, only after 10 iterations. This is the current default implementation of `EM` method,<jupyter_code>fit5 = mvsusie(dat$X,dat$y,L=L,prior_variance=dat$V, 
              compute_objective=T, estimate_residual_variance=F, 
              estimate_prior_variance=T, estimate_prior_method='EM',
              track_fit=T)<jupyter_output><empty_output><jupyter_text>The convergence is now faster,<jupyter_code>fit5$niter
fit5$elbo<jupyter_output><empty_output><jupyter_text>This now achieves the same ELBO as `optim` method.<jupyter_code>fit5$V<jupyter_output><empty_output><jupyter_text>Here, two effects are captured, with very similar result compared to `optim`; **but they appear in effects 1 and 9!**<jupyter_code>susieR::susie_plot(fit5, 'PIP', b = rowSums(dat$b))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08008684115305753, 'fraction_numerical': 0.023157640815341936, 'mean_word_length': 5.607171314741036, 'pattern_counts': {'":': 0, '<': 68, '<?xml version=': 0, '>': 69, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '116975', 'license': 'MIT', 'n_tokens_mistral': 2794, 'n_tokens_neox': 2581, 'n_words': 840, 'path': 'inst/prototypes/elbo_prior_est.ipynb', 'repo_name': 'stephenslab/mmbr', 'size': '8291'}