EM update for prior variance when prior cannot be directly inverted¶

The EM approach previously shown for estimating prior variance scalar has a problem: its current form does not work when input prior matrices has non-invertable component. I used generalized inverse in the code but it does not work.

Illustration¶

devtools::load_all('~/GIT/software/mvsusieR')
set.seed(1)
dat = mvsusie_sim1(n=50,p=50,r=2,s=2)
L = 10

Loading mvsusieR

Loading required package: mashr

Loading required package: ashr

V = matrix(1,2,2)

fit = mvsusie(dat$X,dat$y,L=L,prior_variance=V, 
              compute_objective=T, estimate_residual_variance=F, 
              estimate_prior_variance=T, estimate_prior_method = 'EM',
              track_fit=TRUE)

Loading mvsusieR

Warning message in m$get_objective(dump = TRUE):
“Objective is not non-decreasing”

fit$elbo

fit$prior_history

It does not work at all -- in the end all esimates of V are very small and no signal is captured.

A workaround implemented here avoids inverse on input prior when it is not invertable. The disadvantage here is that some quantities have to be recomputed which could otherwise be avoided using original direct inverse implementation.

Verify the workaround agrees with direct method when prior is invertable¶

But before using this approach we need to establish that the result is the same as the original implementation when prior is invertable. To achieve this I use a diagonal prior matrix which is invertable. But to trigger the workaround code, I manually mess up this line to take a non-existing function for inverse, so I can run the workaround even when prior is invertable.

V = diag(2)

# This is the direct inverse version
devtools::load_all('~/GIT/software/mvsusieR')
fit1 = mvsusie(dat$X,dat$y,L=L,prior_variance=V, 
              compute_objective=T, estimate_residual_variance=F, 
              estimate_prior_variance=T, estimate_prior_method = 'EM',
              track_fit=TRUE)

fit1$V

# This is the indirect inverse version
# after manually messed up the version for direct inverse
devtools::load_all('~/GIT/software/mvsusieR')
fit2 = mvsusie(dat$X,dat$y,L=L,prior_variance=V, 
              compute_objective=T, estimate_residual_variance=F, 
              estimate_prior_variance=T, estimate_prior_method = 'EM',
              track_fit=TRUE)

fit2$V

Let's try one more with the original simulated prior,

# This is the direct inverse version
devtools::load_all('~/GIT/software/mvsusieR')
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 = 'EM',
              track_fit=TRUE)

fit1$V

# This is the indirect inverse version
# after manually messed up the version for direct inverse
devtools::load_all('~/GIT/software/mvsusieR')
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$V

It seems to work fine.

Comparing the new EM implementation with `optim`¶

V = matrix(1,2,2)
devtools::load_all('~/GIT/software/mvsusieR')
fit3 = mvsusie(dat$X,dat$y,L=L,prior_variance=V, 
              compute_objective=T, estimate_residual_variance=F, 
              estimate_prior_variance=T, estimate_prior_method = 'EM',
              track_fit=TRUE)

Loading mvsusieR

fit3$elbo

fit3$V

It seems to be working, and detects a signal with the 6th effect.

To compare I also run optim directly,

fit4 = mvsusie(dat$X,dat$y,L=L,prior_variance=V,
              compute_objective=T, estimate_residual_variance=F, 
              estimate_prior_variance=T, estimate_prior_method = 'optim',
              track_fit=TRUE)

fit4$elbo

fit4$V

It takes less iterations (though not much less computational time, if not more) and the ELBO achieved is slightly higher.

An additional comparison between EM and optim can be found here. optim seems better in that scenario.

Verify that MASH prior is also correctly implemented¶

When V is not invertable,

m_init = MashInitializer$new(list(V), 1)
fit5 = mvsusie(dat$X,dat$y,L=L,prior_variance=m_init, 
              compute_objective=T, estimate_residual_variance=F, precompute_covariances=T,
              estimate_prior_variance=T, estimate_prior_method = 'EM',
              track_fit=TRUE)

fit5$V

When V is invertable,

m_init = MashInitializer$new(list(diag(2)), 1)
fit6 = mvsusie(dat$X,dat$y,L=L,prior_variance=m_init, 
              compute_objective=T, estimate_residual_variance=F, precompute_covariances=T,
              estimate_prior_variance=T, estimate_prior_method = 'EM',
              track_fit=TRUE)

fit6$V

This result agrees with above from BMR. Now I comment out the line to precompute inverse so it will use the trick around inverse implemented for MASH,

devtools::load_all('~/GIT/software/mvsusieR')
m_init = MashInitializer$new(list(diag(2)), 1)
fit7 = mvsusie(dat$X,dat$y,L=L,prior_variance=m_init,
              compute_objective=T, estimate_residual_variance=F, precompute_covariances=T,
              estimate_prior_variance=T, estimate_prior_method = 'EM',
              track_fit=TRUE)

Loading mvsusieR

NULL

fit7$V

EM update for prior variance when prior cannot be directly inverted¶

Illustration¶

Verify the workaround agrees with direct method when prior is invertable¶

Comparing the new EM implementation with optim¶

Verify that MASH prior is also correctly implemented¶

Comparing the new EM implementation with `optim`¶