No problem with EE model¶

I discovered while working on the example of "two SNP dispute" that for Exchangable Effect (EE) model, that is, setting MASH alpha=0, there seems convergence issues with the algorithm, for the specific simulation below. I cannot reproduce the problem with alpha=1 or with the 2nd block having positive sign of effect size ...

This notebook shows that the issue cannot be reproduced when the 2nd block has positive sign of effects

%cd /home/gaow/tmp/13-May-2019

/home/gaow/Documents/TempDir/13-May-2019

Find negatively correlation SNPs¶

genotype = readRDS('Multi_Tissues.ENSG00000145214.RDS')$X

Simulate multivariate phenotypes¶

P = ncol(genotype)
R = 6
eff_factor = 1.5
snp1 = 184
snp2 = 354
cor(genotype[,snp1], genotype[,snp2])

The two SNPs are negatively correlated.

set.seed(1)
b = matrix(0, P, R)
sharing = matrix(0.75, (R/2), (R/2))
diag(sharing) = 1
print(sharing)

     [,1] [,2] [,3]
[1,] 1.00 0.75 0.75
[2,] 0.75 1.00 0.75
[3,] 0.75 0.75 1.00

Here above I have a matrix block. I'm going to create a blocky structure next such that the two halves of traits are drawn from this sharing pattern, but there is no sharing in between them.

b[snp1, 1:(R/2)] = abs(MASS::mvrnorm(1, rep(0,(R/2)), sharing)) / eff_factor
b[snp2, (R/2+1):R] = abs(MASS::mvrnorm(1, rep(0,(R/2)), sharing)) / eff_factor
print(b[snp1, 1:(R/2)])
print(b[snp2, (R/2+1):R])

[1] 0.1951948 0.3530102 0.5955379
[1] 1.136650 1.039006 0.736915

g = genotype %*% b

e = MASS::mvrnorm(1, rep(0,R), diag(R) * 1)

Y = g + e

Prepare MASH mixture prior¶

U1 = matrix(0,R,R)
U1[1:(R/2),1:(R/2)] = sharing
U2 = matrix(0,R,R)
U2[(R/2+1):R,(R/2+1):R] = sharing
Ulist = list(U1=U1, U2=U2)
scaling = c(0.5,1) / eff_factor
mash_init = mvsusieR:::MashInitializer$new(Ulist, scaling, alpha = 0)

Notice here I'm using the EE model by setting alpha = 0.

Fit MV-SuSiE¶

Using EE model it takes quite a while to converge,

start_time <- Sys.time()
res = mvsusieR::susie(genotype,Y,
                  L=10,V=mash_init,
                  precompute_covariances=FALSE,
                  compute_objective=FALSE)
end_time <- Sys.time()
print(end_time - start_time)

Time difference of 2.75596 mins

Number of iterations,

res$niter

snps = sub("_[A-Z]*_[A-Z]*_b38", "", colnames(genotype))

rownames(res$coef) = c('intercept', snps)

Visualize results¶

pdf('susie_plot_toy.pdf', width=8, height=4)
susieR::susie_plot(res,y='PIP', main = 'Cross-condition Posterior Inclusion Probability', xlab = 'SNP positions', add_legend = F)
dev.off()

p = mvsusieR::mvsusie_plot(res)

Suggested PDF canvas width: 1 height: 5.4

pdf('mvsusie_plot_toy.pdf', width = 8, height = 4)
print(p$plot)
dev.off()

%preview susie_plot_toy.pdf -s png --dpi 150

%preview mvsusie_plot_toy.pdf -s png --dpi 150

colnames(genotype)[snp1]

colnames(genotype)[snp2]