Description Usage Arguments Details Value Author(s) Examples
View source: R/corrected_cov_functions.R
Obtain corrected coverage estimate using estimated effect sizes and their standard errors (limiting simulations used for estimation to those with correct nvar)
1 2 3 4 5 6 7 8 9 10 11 12 | corrcov_nvar_bhat(
bhat,
V,
N0,
N1,
Sigma,
nvar,
thr,
W = 0.2,
nrep = 10000,
pp0min = 0.001
)
|
bhat |
Estimated effect sizes from single-SNP logistic regressions |
V |
Variance of estimated effect sizes |
N0 |
Number of controls |
N1 |
Number of cases |
Sigma |
SNP correlation matrix |
nvar |
The number of variants that simulated credible sets used for estimation should contain |
thr |
Minimum threshold for fine-mapping experiment |
W |
Prior for the standard deviation of the effect size parameter, beta (default 0.2) |
nrep |
The number of simulated posterior probability systems to consider for the corrected coverage estimate (nrep = 10000 default due to trimming) |
pp0min |
Only average over SNPs with pp0 > pp0min |
This function requires the marginal summary statistics from GWAS and an nvar value. It should only be used when nvar is very low ($<3$) and there is some evidence to suggest that only simulated credible sets with this nvar value should be used to derive the corrected coverage estimate.
Corrected coverage estimate
Anna Hutchinson
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 | set.seed(1)
nsnps <- 100
N0 <- 5000 # number of controls
N1 <- 5000 # number of cases
## 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
corrcov_nvar_bhat(bhat = bhats, V = varbeta, N0, N1, Sigma = LD, thr = 0.95, nvar = 1, nrep = 1000)
# note that nrep should be at least the default value (nrep = 10000) but is
# lower here for speed of computation
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.