corrected_cs_bhat.RdCorrected credible set using estimated effect sizes and their standard errors
corrected_cs_bhat(bhat, V, N0, N1, Sigma, W = 0.2, lower = 0, upper = 1, desired.cov, acc = 0.005, max.iter = 20, pp0min = 0.001)
| bhat | Estimated effect sizes |
|---|---|
| V | Prior variance of estimated effect sizes |
| N0 | Number of controls |
| N1 | Number of cases |
| Sigma | Correlation matrix of SNPs |
| W | Prior for the standard deviation of the effect size parameter, beta (default 0.2) |
| lower | Lower threshold (default = 0) |
| upper | Upper threshold (default = 1) |
| desired.cov | The desired coverage of the causal variant in the credible set |
| acc | Accuracy of corrected coverage to desired coverage (default = 0.005) |
| max.iter | Maximum iterations (default = 20) |
| pp0min | Only average over SNPs with pp0 > pp0min |
List of variants in credible set, required threshold, the corrected coverage and the size of the credible set
# \donttest{ # this is a long running example # In this example, the function is used to find a corrected 95% credible set # using bhats and their standard errors, that is the smallest set of variants # required such that the resultant credible set has coverage close to (/within # some accuracy of) the "desired coverage" (here set to 0.95). Max.iter parameter # defines the maximum number of iterations to try in the root bisection algorithm, # this should be increased to ensure convergence to the desired coverage, but is set # to 1 here for speed (and thus the resultant credible set will not be accurate). set.seed(18) nsnps <- 100 N0 <- 500 # number of controls N1 <- 500 # number of cases # simulate fake haplotypes to obtain MAFs and LD matrix ## generate example LD matrix library(mvtnorm) nsamples = 1000 simx <- function(nsnps, nsamples, S, maf=0.1) { mu <- rep(0,nsnps) rawvars <- rmvnorm(n=nsamples, mean=mu, sigma=S) pvars <- pnorm(rawvars) x <- qbinom(1-pvars, 1, maf) } S <- (1 - (abs(outer(1:nsnps,1:nsnps,`-`))/nsnps))^4 X <- simx(nsnps,nsamples,S) LD <- cor2(X) maf <- colMeans(X) varbeta <- Var.data.cc(f = maf, N = N0 + N1, s = N1/(N0+N1)) bhats = rnorm(nsnps,0,0.2) # log OR names(bhats) <- seq(1,length(bhats)) corrected_cs_bhat(bhat = bhats, V = varbeta, N0, N1, Sigma = LD, desired.cov = 0.9, max.iter = 1)#> [1] "thr: 0.5 , cov: 0.877073980331893"#> $credset #> [1] "38" "90" "92" "10" "62" #> #> $req.thr #> [1] 0.5 #> #> $corr.cov #> [1] 0.877074 #> #> $size #> [1] 0.5140384 #># max.iter set low for speed, should be set to at # least the default to ensure convergence to desired coverage # }