stephenslab / susieR

R package for "sum of single effects" regression.
https://stephenslab.github.io/susieR
Other
174 stars 44 forks source link

add functionality to sample from susie approximate posterior #96

Closed KangchengHou closed 3 years ago

KangchengHou commented 4 years ago

Hi,

I add a file R/sample_posterior.R which performs approximate posterior samping to fix issue #66. We could discuss in more details in this pull request.

Related to this, we should remove effects having estimated prior variance equals zero in the following two functions for consistency.

https://github.com/stephenslab/susieR/blob/3ceeb94c51c93dc9b6720381bad083273f215132/R/susie_utils.R#L491-L500

The follows is an informal unit test. I don't know what exact testing we should have here. Since samples are generated randomly, so the assertion would be probabilistic. I would like to know your ideas here.

library(susieR)

# simulation adopted from https://stephenslab.github.io/susie-paper/manuscript_results/motivating_example.html
set.seed(1234)
n = 100
p = 10
b = rep(0,p)
b[2] = 1
b[8] = 1
X = matrix(rnorm(n*p),nrow=n,ncol=p)
X[,2] = X[,4]
X[,6] = X[,8]
y = X %*% b + rnorm(n)

# fit susie
susie_fit = susieR::susie(X, y, L=5,
    estimate_residual_variance=TRUE, 
    scaled_prior_variance=0.2,
    tol=1e-3, track_fit=TRUE, min_abs_corr=0.1)

# get posterior samples
num_samples <- 100000
posterior_samples <- susie_get_posterior_samples(susie_fit, num_samples)
print(susie_get_pip(susie_fit))
print(rowMeans(posterior_samples$gamma))
print(susie_get_posterior_mean(susie_fit))
print(rowMeans(posterior_samples$b))

This is the first time i contributed to a open source repository so feel free to let me know if there is any code practice i should follow here :)