stephenslab / mashr

An R package for multivariate adaptive shrinkage.
https://stephenslab.github.io/mashr
Other
87 stars 19 forks source link

obtaining both terms going into the local false sign rate calculation #64

Open hillarykoch opened 5 years ago

hillarykoch commented 5 years ago

This issue is more related to my own personal interest-

I have developed a statistical test that I am now trying to apply to some real data. I am trying it with a couple of different models, but one of those includes mashr. While it has been very easy to use mashr to fit a model, in order to apply my statistical test, I need both terms that go into computing the local false sign rate for any observation. That is, as far as I can tell, mashr returns the local false sign rate, but I need to know both P(β ≥ 0) and P(β ≤ 0). I can't figure out how I can get these terms easily. Would it be possible to make these more accessible to the user?

Thanks!

gaow commented 5 years ago

@hillarykoch thanks for the message. We do have the information you are looking for, it's just we do not have an interface to extract them. Take the minimal working example below:

> library(ashr)
> library(mashr)
> set.seed(1)
> simdata = simple_sims(500,5,1)
> data = mash_set_data(simdata$Bhat, simdata$Shat)
> U.c = cov_canonical(data)  
> m.c = mash(data, U.c)
 - Computing 2000 x 151 likelihood matrix.
 - Likelihood calculations took 0.03 seconds.
 - Fitting model with 151 mixture components.
 - Model fitting took 0.24 seconds.
 - Computing posterior matrices.
 - Computation allocated took 0.01 seconds.
> names(m.c)
[1] "result"            "loglik"            "vloglik"          
[4] "null_loglik"       "alt_loglik"        "fitted_g"         
[7] "posterior_weights" "alpha"            
> names(m.c$result)
[1] "PosteriorMean" "PosteriorSD"   "lfdr"          "NegativeProb" 
[5] "lfsr"

You see here lfdr is the probability of \beta = 0, NegativeProb is probability of \beta < 0. We can potentially add some get_* function to extract them formally, but for now you can work out the quantities you need from there. Does it help?

hillarykoch commented 5 years ago

Yes, thank you! I just did

get_prob_positive <- function(m) {
    prob0 <- m$result$lfdr
    probneg <- m$result$NegativeProb
    1 - prob0 - probneg
}

I used fdr and NegativeProb and tested against the lfsr output by mash. It all seems to agree.

Thanks again!