corrcov_nvar_bhat {corrcoverage}

Corrected coverage estimate using estimated effect sizes and their standard errors (fixing nvar)

Description

Obtain corrected coverage estimate using estimated effect sizes and their standard errors (limiting simulations used for estimation to those with correct nvar)

Usage

corrcov_nvar_bhat(
  bhat,
  V,
  N0,
  N1,
  Sigma,
  nvar,
  thr,
  W = 0.2,
  nrep = 10000,
  pp0min = 0.001
)

Arguments

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

Details

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.

Value

Corrected coverage estimate

Author(s)

Anna Hutchinson

Examples

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