stephenslab / susieR

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

susie sometimes fail to compute univariate z score when refine = T #175

Open hsun3163 opened 1 year ago

hsun3163 commented 1 year ago

when running

test_output = susieR::susie(susie_list[[1]]$input_data$X_resid,susie_list[[1]]$input_data$Y_resid,
                           L=10,
                           max_iter=1000,
                           estimate_residual_variance=TRUE,
                           estimate_prior_variance=TRUE,
                           refine=TRUE,compute_univariate_zscore = TRUE, coverage = 0.95 )
test_output$z

return NULL

however, when running

test_output = susieR::susie(susie_list[[1]]$input_data$X_resid,susie_list[[1]]$input_data$Y_resid,
                           L=10,
                           max_iter=1000,
                           estimate_residual_variance=TRUE,
                           estimate_prior_variance=TRUE,
                           refine=F,compute_univariate_zscore = TRUE, coverage = 0.95 )
test_output$z

The z score can be successfully produced.

I believed this is due to the compute_univariate_zscore = FALSE option in the following code block which overwrite the calculated Z :

m = list()
            for (cs in 1:length(s$sets$cs)) {
                pw_cs = pw_s
                pw_cs[s$sets$cs[[cs]]] = 0
                if (all(pw_cs == 0)) {
                  break
                }
                s2 = susie(X, y, L = L, scaled_prior_variance = scaled_prior_variance, 
                  residual_variance = residual_variance, prior_weights = pw_cs, 
                  s_init = NULL, null_weight = null_weight, standardize = standardize, 
                  intercept = intercept, estimate_residual_variance = estimate_residual_variance, 
                  estimate_prior_variance = estimate_prior_variance, 
                  estimate_prior_method = estimate_prior_method, 
                  check_null_threshold = check_null_threshold, 
                  prior_tol = prior_tol, coverage = coverage, 
                  residual_variance_upperbound = residual_variance_upperbound, 
                  min_abs_corr = min_abs_corr, compute_univariate_zscore = FALSE, 
                  na.rm = na.rm, max_iter = max_iter, tol = tol, 
                  verbose = FALSE, track_fit = FALSE, residual_variance_lowerbound = var(drop(y))/10000, 
                  refine = FALSE)
                sinit2 = s2[c("alpha", "mu", "mu2")]
                class(sinit2) = "susie"
                s3 = susie(X, y, L = L, scaled_prior_variance = scaled_prior_variance, 
                  residual_variance = residual_variance, prior_weights = pw_s, 
                  s_init = sinit2, null_weight = null_weight, 
                  standardize = standardize, intercept = intercept, 
                  estimate_residual_variance = estimate_residual_variance, 
                  estimate_prior_variance = estimate_prior_variance, 
                  estimate_prior_method = estimate_prior_method, 
                  check_null_threshold = check_null_threshold, 
                  prior_tol = prior_tol, coverage = coverage, 
                  residual_variance_upperbound = residual_variance_upperbound, 
                  min_abs_corr = min_abs_corr, compute_univariate_zscore = FALSE, 
                  na.rm = na.rm, max_iter = max_iter, tol = tol, 
                  verbose = FALSE, track_fit = FALSE, residual_variance_lowerbound = var(drop(y))/10000, 
                  refine = FALSE)
                m = c(m, list(s3))
            }
pcarbo commented 1 year ago

@zouyuxin Is there a reason why you set compute_univariate_zscore = FALSE with refine = TRUE? Are you worried it will slow down the refinement procedure?

zouyuxin commented 1 year ago

Oh, sorry, it's a bug. The compute_univariate_zscore = FALSE for s3 should be compute_univariate_zscore = compute_univariate_zscore.

pcarbo commented 1 year ago

@hsun3163 Does the latest version fix the problem?