check_null_threshold

[jerome-f]

Hi,

what is the behavior of check_null_threshold=0 from the previous commits this was changed from check_null_tol=0.1, how does it affect the credible sets? The explanation in docs is cryptic. Is setting check_null_threshold=0.1 more strict? i.e. will lead to less false discoveries?

Best Jerome

[jerome-f]

reason I am bringing this up is cos I am running into a issue

Error in if (loglik(0, betahat, shat2, prior_weights) + check_null_threshold >=  : 
  missing value where TRUE/FALSE needed
Calls: susie ... single_effect_regression -> optimize_prior_variance
Execution halted

susie ran with individual level data (pheno and dosage) with the following setting:

susie( X, y, L = 30, \
       estimate_prior_method = "EM", \
       verbose = T, refine = T, max_iter = 1000) 

[pcarbo]

@jerome-f Just to check whether the problem is indeed due to estimating the prior variance, does the error go away if you set estimate_prior_variance = FALSE?

Setting check_null_threshold > 0 will have the effect of pruning the CSs slightly more aggressively.

[jerome-f]

Setting the gave estimate_prior_variance = FALSE me an error in ELBO

[1] "objective:-34159.5469366425"
[1] "objective:-34159.546639484"
Error in if ((elbo[i + 1] - elbo[i]) < tol) { : 
  missing value where TRUE/FALSE needed
Calls: susie -> susie
Execution halted 

[pcarbo]

@jerome-f So the issue does not seem to be the estimation of the prior. If you set L to a smaller value, such as L = 10, does it still produce an error?

[jerome-f]

Set the L=5 and still get the same errors. elbo[i+1] is going to NA, and likely checking for it somewhere in the eblo module or single_effect_regression module might help catching this.

[pcarbo]

Okay, thanks for troubleshooting. If you can share your data set here, or by email if you prefer, we would be happy to track down the exact source of the error and update susieR with a more useful error message.

[jerome-f]

@pcarbo attaching the data here, I have anonymized the sample and the variant ids. column_1 is sample_id, column_2 is Y (phenotype), column_3...column_N is X (per snp dosage)

[jerome-f]

Just granted you access. Let me know if you have trouble accessing

[pcarbo]

@jerome-f Thanks — can you please instead share this as an .RData file with the inputs needed to call susie (X, y)?

[jerome-f]

Also the same data was finemapped with FINEMAP resulting in 9 credible sets (8.9). Also this data is among the loci that takes forever to converge

[jerome-f]

@pcarbo shared file in .rds

the following should work ...

data <- readRDS("annonymized_data_for_test.rds")
fit <- susie(as.matrix( data[, -seq( 2 ), with = F ] ), data[[ 2 ]], 
               L = 10, estimate_prior_method = "EM", 
               verbose = T, refine = T, max_iter = 1000)

[pcarbo]

@jerome-f I ran this and it worked (no error):

library(susieR)
set.seed(1)
dat <- readRDS("annonymized_data_for_test.rds")
X <- as.matrix(dat[,-c(1,2)])
markers <- 7300:7400
X <- X[,markers]
y <- dat$PHENOTYPE
fit <- susie(X,y,L = 10, estimate_prior_method = "EM",
             verbose = TRUE,refine = FALSE,max_iter = 10)

However, when I set refine = TRUE then I got an error:

Error in if (loglik(0, betahat, shat2, prior_weights) + check_null_threshold >=  :
  missing value where TRUE/FALSE needed

Can you confirm?

Note I am using susieR 0.12.37.

[jerome-f]

@pcarbo I will try with refine = FALSE I ran with refine=TRUE so something is going wrong in the refinement step?

[pcarbo]

@jerome-f Yes, exactly, there must be a bug somewhere in the refinement step. If you don't mind, I will leave this Issue open as a reminder to fix the bug.

[jerome-f]

Thanks and yes please keep this open

[pcarbo]

@jerome-f Apologies for the delay. This should be fixed now. Can you please try again with the latest version of susieR on GitHub?

[jerome-f]

Will try it out. Thanks for fixing this.

[jerome-f]

Hi Peter,

The bug is fixed and I think we can close this.

Best Jerome

[pcarbo]

Thanks for checking.