Corrected credible set using Z-scores and MAFs

Corrected credible set using Z-scores and MAFs

corrected_cs(z, f, N0, N1, Sigma, W = 0.2, lower = 0, upper = 1,
  desired.cov, acc = 0.005, max.iter = 20, pp0min = 0.001)

Arguments

z	Z-scores
f	Minor allele frequencies
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

Value

List of variants in credible set, required threshold, the corrected coverage and the size of the credible set

Examples

# \donttest{
 # this is a long running example

# In this example, the function is used to find a corrected 95% credible set
# using Z-scores and MAFs, 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(2)
nsnps = 200
N0 = 1000
N1 = 1000
z_scores <- rnorm(nsnps, 0, 1) # simulate a vector of Z-scores

## 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)

names(z_scores) <- seq(1,length(z_scores))

corrected_cs(z = z_scores, f = maf, N0, N1, Sigma = LD, desired.cov = 0.9, max.iter = 1)
#> [1] "thr:  0.5 , cov:  0.471664648318262"
#> $credset
#>  [1] "26"  "16"  "85"  "139" "168" "21"  "167" "52"  "106" "62"  "134" "108"
#> [13] "67"  "46"  "9"   "56"  "42"  "97"  "24"  "44"  "86"  "61"  "15"  "130"
#> [25] "129" "37"  "123" "82"  "92"  "60"  "180" "91"  "23"  "3"   "95"  "140"
#> [37] "80"  "84"  "114" "94"  "156" "193" "165" "77"  "155" "191" "54"  "182"
#> [49] "169" "121"
#> 
#> $req.thr
#> [1] 0.5
#> 
#> $corr.cov
#> [1] 0.4716646
#> 
#> $size
#> [1] 0.5004184
#> 
# max.iter set low for speed, should be set to at least
# the default to ensure convergence to desired coverage
# }