fix NA handling and website

R/AllGenerics.R

@@ -655,7 +655,9 @@ setGeneric("getRegion", function(x, ...) standardGeneric("getRegion"))
 #' @param x A \code{QtlDataset} object.
 #' @param ... Selection arguments: \code{traitId}, \code{region},
 #'   \code{cisWindow}, \code{phenotypeCovariatesToRemove},
-#'   \code{genotypeCovariatesToRemove}.
+#'   \code{genotypeCovariatesToRemove}, and \code{covariateNaAction}
+#'   (\code{"impute"}, the default, mean-imputes missing covariate cells;
+#'   \code{"drop"} removes samples with any missing covariate).
 #' @return A numeric matrix (samples x variants).
 #' @export
 setGeneric("getResidualizedGenotypes",
@@ -669,7 +671,9 @@ setGeneric("getResidualizedGenotypes",
 #' @param ... Selection arguments: \code{contexts} (required),
 #'   \code{traitId}, \code{region},
 #'   \code{phenotypeCovariatesToRemove},
-#'   \code{genotypeCovariatesToRemove}.
+#'   \code{genotypeCovariatesToRemove}, and \code{covariateNaAction}
+#'   (\code{"impute"}, the default, mean-imputes missing covariate cells;
+#'   \code{"drop"} removes samples with any missing covariate).
 #' @return A named list of numeric matrices keyed by context.
 #' @export
 setGeneric("getResidualizedPhenotypes",

R/QtlDataset.R

@@ -844,6 +844,37 @@ setMethod("getPhenotypeCovariates", "QtlDataset",
   do.call(cbind, blocks)
 }
 
+# Internal: resolve missing values in the residualization covariate matrix
+# `C` (samples x covariates) before it reaches `stats::lm.fit`, which does
+# not tolerate NA in the design and would otherwise error. Two strategies:
+#   - "impute" (default): replace each covariate column's NA cells with that
+#     column's observed (non-NA) mean. A wholly-missing column has an
+#     undefined mean and is filled with 0, so it contributes nothing to the
+#     fit rather than poisoning every row.
+#   - "drop": complete-case -- remove any sample (row) carrying an NA in any
+#     covariate. The caller's downstream sample intersection then narrows the
+#     response matrix (G or Y) to the retained samples.
+# Returns the cleaned matrix (fewer rows possible under "drop"), or `C`
+# unchanged when it is NULL or already NA-free.
+.qtlHandleCovariateNa <- function(C, action = c("impute", "drop")) {
+  action <- match.arg(action)
+  if (is.null(C) || !anyNA(C)) return(C)
+  if (action == "drop") {
+    keep <- rowSums(is.na(C)) == 0L
+    return(C[keep, , drop = FALSE])
+  }
+  for (j in seq_len(ncol(C))) {
+    col <- C[, j]
+    na <- is.na(col)
+    if (any(na)) {
+      mu <- mean(col[!na])
+      col[na] <- if (is.finite(mu)) mu else 0
+      C[, j] <- col
+    }
+  }
+  C
+}
+
 # Internal: resolve a (convenience, precise) flag pair to a single boolean.
 # `missing*` arguments are passed as the result of `missing()` evaluated in
 # the calling method to detect whether the user explicitly set the value.
@@ -879,7 +910,9 @@ setMethod("getResidualizedGenotypes", "QtlDataset",
            residualizeGenotypeCovariates  = TRUE,
            residualizePhenotypeCovariatesFromGenotypes = NULL,
            residualizeGenotypeCovariatesFromGenotypes  = NULL,
+           covariateNaAction = c("impute", "drop"),
            ...) {
+    covariateNaAction <- match.arg(covariateNaAction)
     if (missing(contexts) || is.null(contexts) || length(contexts) == 0L) {
       stop("`contexts` is required for getResidualizedGenotypes(QtlDataset). ",
            "Use getContexts(x) to list the available contexts. ",
@@ -924,6 +957,7 @@ setMethod("getResidualizedGenotypes", "QtlDataset",
       genoSelection  = genoSel,
       includePheno   = includePheno,
       includeGeno    = includeGeno)
+    C <- .qtlHandleCovariateNa(C, covariateNaAction)
     if (!is.null(C)) {
       common <- intersect(rownames(G), rownames(C))
       if (length(common) == 0L) {
@@ -948,10 +982,12 @@ setMethod("getResidualizedPhenotypes", "QtlDataset",
            residualizePhenotypeCovariatesFromPhenotypes = NULL,
            residualizeGenotypeCovariatesFromPhenotypes  = NULL,
            naAction = c("keep", "drop", "impute"),
+           covariateNaAction = c("impute", "drop"),
            outlierAction = c("keep", "drop"),
            outlierPvalThreshold = 1e-3,
            ...) {
     naAction <- match.arg(naAction)
+    covariateNaAction <- match.arg(covariateNaAction)
     outlierAction <- match.arg(outlierAction)
     if (missing(contexts) || is.null(contexts) || length(contexts) == 0L) {
       stop("`contexts` is required for getResidualizedPhenotypes().")
@@ -994,6 +1030,7 @@ setMethod("getResidualizedPhenotypes", "QtlDataset",
       genoSelection  = genoSel,
       includePheno   = includePheno,
       includeGeno    = includeGeno)
+    C <- .qtlHandleCovariateNa(C, covariateNaAction)
 
     out <- lapply(contexts, function(ctx) {
       se <- Yraw[[ctx]]

_pkgdown.yml

@@ -115,6 +115,8 @@ reference:
       reference.
     contents:
       - AnnotationMatrix
+      - CtwasResult
+      - CtwasResultEntry
       - FineMappingEntry
       - GenotypeHandle
       - GwasFineMappingResult
@@ -125,6 +127,7 @@ reference:
       - QtlDataset
       - QtlFineMappingResult
       - QtlSumStats
+      - qtlSumStatsFromZMatrix
       - SldscData
       - TwasWeights
       - TwasWeightsEntry
@@ -167,14 +170,22 @@ reference:
       - getSusieFit
       - getTopLoci
       - getVariantIds
+      - resolveWeights
 
   - subtitle: "FineMappingResultBase"
     contents:
       - getLdSketch
       - getMethodNames
+      - getRegion
       - getStudy
       - writeSumstatsVcf
 
+  - subtitle: "CtwasResult / CtwasResultEntry"
+    contents:
+      - getCtwasParam
+      - getFinemap
+      - getSusieAlpha
+
   - subtitle: "GenotypeHandle"
     contents:
       - getFormat
@@ -333,6 +344,7 @@ reference:
       - estCtwasParam
       - screenCtwasRegions
       - finemapCtwasRegions
+      - asCtwasResult
 
   - title: "Quality control"
     desc: >
@@ -391,6 +403,7 @@ reference:
       - regionToDf
       - regionsOverlap
       - findOverlappingRegions
+      - asGranges
 
   - title: "Fine-mapping"
     desc: >
@@ -422,6 +435,7 @@ reference:
       - formatFinemappingOutput
       - lbfToAlpha
       - postprocessFinemappingFits
+      - combineFineMappingResults
 
   - title: "TWAS"
     desc: >
@@ -475,6 +489,7 @@ reference:
     contents:
       - learnTwasWeights
       - ensembleWeights
+      - combineTwasWeights
       - twasWeightsCv
       - twasPredict
       - twasZ
@@ -515,6 +530,7 @@ reference:
       - h2EstimateToSldscTrait
       - isBinarySldscAnnot
       - metaSldscRandom
+      - sldscSubsetMeta
       - qtlEnrichment
 
   - title: "Bundled example datasets"

tests/testthat/test_QtlDataset.R

@@ -795,6 +795,49 @@ test_that(".qtlResolveResidualizationFlag: both set + conflicting errors", {
   )
 })
 
+# ===========================================================================
+# .qtlHandleCovariateNa (covariate missing-value strategy)
+# ===========================================================================
+
+test_that(".qtlHandleCovariateNa: NULL and NA-free inputs pass through unchanged", {
+  expect_null(pecotmr:::.qtlHandleCovariateNa(NULL))
+  C <- matrix(c(1, 2, 3, 4), nrow = 2)
+  expect_identical(pecotmr:::.qtlHandleCovariateNa(C), C)
+  expect_identical(pecotmr:::.qtlHandleCovariateNa(C, "drop"), C)
+})
+
+test_that(".qtlHandleCovariateNa: impute (default) fills NA with column mean", {
+  C <- matrix(c(1, NA, 3,    # col1 observed mean = (1 + 3)/2 = 2
+                10, 20, NA), # col2 observed mean = (10 + 20)/2 = 15
+              nrow = 3, ncol = 2,
+              dimnames = list(c("s1", "s2", "s3"), c("a", "b")))
+  out <- pecotmr:::.qtlHandleCovariateNa(C)
+  expect_false(anyNA(out))
+  expect_equal(dim(out), dim(C))
+  expect_equal(out["s2", "a"], 2)
+  expect_equal(out["s3", "b"], 15)
+  # observed cells are untouched
+  expect_equal(out["s1", "a"], 1)
+})
+
+test_that(".qtlHandleCovariateNa: impute fills a wholly-missing column with 0", {
+  C <- matrix(c(NA_real_, NA_real_, 5, 7), nrow = 2, ncol = 2,
+              dimnames = list(c("s1", "s2"), c("a", "b")))
+  out <- pecotmr:::.qtlHandleCovariateNa(C, "impute")
+  expect_false(anyNA(out))
+  expect_equal(out[, "a"], c(s1 = 0, s2 = 0))
+})
+
+test_that(".qtlHandleCovariateNa: drop removes any sample with a missing covariate", {
+  C <- matrix(c(1, NA, 3,
+                4, 5, 6),
+              nrow = 3, ncol = 2,
+              dimnames = list(c("s1", "s2", "s3"), c("a", "b")))
+  out <- pecotmr:::.qtlHandleCovariateNa(C, "drop")
+  expect_equal(rownames(out), c("s1", "s3"))
+  expect_false(anyNA(out))
+})
+
 # ===========================================================================
 # getResidualizedGenotypes (QtlDataset)
 # ===========================================================================
@@ -914,6 +957,34 @@ test_that("getResidualizedGenotypes: includes genotype covariates when supplied"
   }
 })
 
+test_that("getResidualizedGenotypes: mean-imputes missing covariates by default", {
+  gc <- matrix(rnorm(12 * 2), nrow = 12, ncol = 2,
+               dimnames = list(paste0("s", 1:12), c("pc1", "pc2")))
+  gc[3, "pc1"] <- NA  # missing covariate cell
+  qd <- .qr_makeDataset(contexts = "brain", geno_cov = gc)
+  local_mocked_bindings(extractBlockGenotypes = .qr_mockExtractor(),
+                        .package = "pecotmr")
+  # Default covariateNaAction = "impute": no error, all 12 samples retained.
+  G <- getResidualizedGenotypes(qd, contexts = "brain",
+                                genotypeCovariatesToResidualize = c("pc1", "pc2"))
+  expect_equal(nrow(G), 12L)
+  expect_false(anyNA(G))
+})
+
+test_that("getResidualizedGenotypes: covariateNaAction='drop' removes samples with missing covariates", {
+  gc <- matrix(rnorm(12 * 2), nrow = 12, ncol = 2,
+               dimnames = list(paste0("s", 1:12), c("pc1", "pc2")))
+  gc[3, "pc1"] <- NA
+  qd <- .qr_makeDataset(contexts = "brain", geno_cov = gc)
+  local_mocked_bindings(extractBlockGenotypes = .qr_mockExtractor(),
+                        .package = "pecotmr")
+  G <- getResidualizedGenotypes(qd, contexts = "brain",
+                                genotypeCovariatesToResidualize = c("pc1", "pc2"),
+                                covariateNaAction = "drop")
+  expect_equal(nrow(G), 11L)
+  expect_false("s3" %in% rownames(G))
+})
+
 # ===========================================================================
 # getResidualizedPhenotypes (QtlDataset)
 # ===========================================================================