Open kkdey opened 1 year ago
@kkdey I think it works if you set R to a "dummy value", say, the identity matrix. Here's an example:
set.seed(1)
n = 200
p = 1000
beta = rep(0,p+1)
beta[1] = 1
X = matrix(rnorm(n*p),nrow = n,ncol = p)
X = cbind(X,X[,1] + 0.1*rnorm(n))
X = scale(X,center = TRUE,scale = TRUE)
y = drop(X %*% beta + rnorm(n))
input_ss = compute_suff_stat(X,y,standardize = TRUE)
ss = univariate_regression(X,y)
R = with(input_ss,cov2cor(XtX))
zhat = with(ss,betahat/sebetahat)
res1 = susie_rss(zhat,R,n = n,L = 1,min_abs_corr = 0)
res2 = susie_rss(zhat,R = diag(p+1),n = n,L = 1,min_abs_corr = 0)
print(res1$sets)
print(res2$sets)
Both res1
and res2
should give you the same CS.
Note that you need to remove the "purity" filter for this to work correctly because of course the second call to susie_rss
does not have the right correlations.
Thanks Peter for the demo code and explanation. Will try it like you mention. Will just keep this thread open for a bit in case I face any issues in running it this way.
Hi, @pcarbo I'm wondering will susie_rss select the variant with the lowest p-value or highest z-score when we set L=1. I'm curious because since there is no LD reference in this case. Thanks!
@PhoebeGuo97 In most cases, the smallest p-value should correspond to the largest PIP in this case (so long as the association is "strong enough," i.e., significant).
Hi all,
Is it possible to have a susie_rss() model with 1 causal variant assumption, that can be applied on the GWAS summary statistics without needing to input an LD reference panel?