diff --git a/DESCRIPTION b/DESCRIPTION index dd53c04..8ff98c2 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: AllelicSeries Title: Allelic Series Test -Version: 0.1.0.0 +Version: 0.1.0.1 Authors@R: c(person(given = "Zachary", family = "McCaw", diff --git a/R/allelic_series_sumstats.R b/R/allelic_series_sumstats.R index c46ac1d..b507c44 100644 --- a/R/allelic_series_sumstats.R +++ b/R/allelic_series_sumstats.R @@ -14,7 +14,7 @@ DEFAULT_WEIGHTS <- c(1, 2, 3) #' @param se (snps x 1) vector of standard errors for the effect sizes. #' @param check Run input checks? Default: TRUE. #' @param eps Epsilon added to the diagonal of the LD matrix if not positive -#' definite. Note, epsilon should increase as the sample size decreases. +#' definite. Note, smaller values increase the chances of a false positive. #' @param lambda Optional genomic inflation factor. Defaults to 1, which #' results in no rescaling. #' @param ld (snps x snps) matrix of correlations among the genetic variants. @@ -45,7 +45,7 @@ ASBTSS <- function( beta, se, check = TRUE, - eps = 1e-4, + eps = 1, lambda = 1, ld = NULL, maf = NULL, @@ -69,7 +69,6 @@ ASBTSS <- function( is_pd <- TRUE } - n_snps <- length(anno) if (is.null(ld)) {ld <- diag(n_snps)} if (is.null(maf)) {maf <- rep(0, n_snps)} @@ -141,7 +140,7 @@ ASBTSS <- function( #' @param se (snps x 1) vector of standard errors for the effect sizes. #' @param check Run input checks? Default: TRUE. #' @param eps Epsilon added to the diagonal of the LD matrix if not positive -#' definite. Note, epsilon should increase as the sample size decreases. +#' definite. Note, smaller values increase the chances of a false positive. #' @param lambda Optional genomic inflation factor. Defaults to 1, which #' results in no rescaling. #' @param maf (snps x 1) vector of minor allele frequencies. Although ideally @@ -170,7 +169,7 @@ ASKATSS <- function( beta, se, check = TRUE, - eps = 1e-4, + eps = 1, lambda = 1, ld = NULL, maf = NULL, @@ -296,7 +295,7 @@ COASTSS <- function( beta, se, check = TRUE, - eps = 1e-4, + eps = 1, lambda = c(1, 1, 1), maf = NULL, ld = NULL, diff --git a/man/ASBTSS.Rd b/man/ASBTSS.Rd index 714c928..6fa81a0 100644 --- a/man/ASBTSS.Rd +++ b/man/ASBTSS.Rd @@ -9,7 +9,7 @@ ASBTSS( beta, se, check = TRUE, - eps = 1e-04, + eps = 1, lambda = 1, ld = NULL, maf = NULL, @@ -28,7 +28,7 @@ within a gene.} \item{check}{Run input checks? Default: TRUE.} \item{eps}{Epsilon added to the diagonal of the LD matrix if not positive -definite. Note, epsilon should increase as the sample size decreases.} +definite. Note, smaller values increase the chances of a false positive.} \item{lambda}{Optional genomic inflation factor. Defaults to 1, which results in no rescaling.} diff --git a/man/ASKATSS.Rd b/man/ASKATSS.Rd index c5b2863..0a33c19 100644 --- a/man/ASKATSS.Rd +++ b/man/ASKATSS.Rd @@ -9,7 +9,7 @@ ASKATSS( beta, se, check = TRUE, - eps = 1e-04, + eps = 1, lambda = 1, ld = NULL, maf = NULL, @@ -27,7 +27,7 @@ within a gene.} \item{check}{Run input checks? Default: TRUE.} \item{eps}{Epsilon added to the diagonal of the LD matrix if not positive -definite. Note, epsilon should increase as the sample size decreases.} +definite. Note, smaller values increase the chances of a false positive.} \item{lambda}{Optional genomic inflation factor. Defaults to 1, which results in no rescaling.} diff --git a/man/COASTSS.Rd b/man/COASTSS.Rd index d79df00..eacbbaf 100644 --- a/man/COASTSS.Rd +++ b/man/COASTSS.Rd @@ -9,7 +9,7 @@ COASTSS( beta, se, check = TRUE, - eps = 1e-04, + eps = 1, lambda = c(1, 1, 1), maf = NULL, ld = NULL, diff --git a/src/sumstats.cpp b/src/sumstats.cpp index 971d3c5..af25d13 100644 --- a/src/sumstats.cpp +++ b/src/sumstats.cpp @@ -91,7 +91,9 @@ SEXP IVWCpp( arma::colvec se_anno = se.elem(key); arma::mat d_anno = arma::diagmat(se_anno); arma::mat v_anno = d_anno * ld_anno * d_anno; - arma::mat v_anno_inv = arma::inv_sympd(v_anno); + + // Use of pseudo-inverse here protects against the case where `v_anno` is singular. + arma::mat v_anno_inv = arma::pinv(v_anno); double weight_sum = arma::accu(v_anno_inv); // Meta-analyze. diff --git a/tests/testthat/test-allelic_series_sumstats.R b/tests/testthat/test-allelic_series_sumstats.R index 0827335..58d5fa0 100644 --- a/tests/testthat/test-allelic_series_sumstats.R +++ b/tests/testthat/test-allelic_series_sumstats.R @@ -111,11 +111,15 @@ test_that("Check case of rank-deficient LD.", { withr::local_seed(101) data <- DGP(n = 1e2, prop_causal = 0) sumstats <- CalcSumstats(data = data) + ld <- sumstats$ld + + # Ensure the matrix is singular. + ld[1, 0] <- ld[0, 1] <- 0 # The LD matrix is singular. expect_false(isPD(sumstats$ld)) - # Without epsilon, matrix inversion will fail. + # Test COAST runs even *without* epsilon. expect_error( COASTSS( anno = sumstats$anno, @@ -125,7 +129,8 @@ test_that("Check case of rank-deficient LD.", { eps = 0, ld = sumstats$ld, maf = sumstats$maf - ) + ), + NA ) # With epsilon, the test runs. @@ -134,14 +139,13 @@ test_that("Check case of rank-deficient LD.", { anno = sumstats$anno, beta = sumstats$sumstats$beta, se = sumstats$sumstats$se, - eps = 1e-4, + eps = 1, ld = sumstats$ld, maf = sumstats$maf ) ) expect_true(all(results@Pvals$pval > 0.05)) - }) @@ -201,4 +205,4 @@ test_that("Check cases where not all annotaiton categories are provided.", { result12 <- RunCOAST(sumstats12) expect_true(all(result12$pval >= 0.05)) -}) \ No newline at end of file +})