Improve handling of singular LD matrices (#22)

Found empirically that the previous regularization (eps = 1e-4) applied to singular LD matrices was insufficient. Increasing to eps = 1 resolves the issue, even in the extreme case that an entire row and column of the LD matrix is 0.

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",

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,

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.

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))
 
-})
+})