From fbe1e03b2f01f6f355f30725aaffa7d9dc1ce8b3 Mon Sep 17 00:00:00 2001 From: Daniel Nachun Date: Fri, 10 Jul 2026 11:22:29 -0700 Subject: [PATCH] fix NA handling and website --- R/AllGenerics.R | 8 +++- R/QtlDataset.R | 37 +++++++++++++++++ _pkgdown.yml | 16 +++++++ man/getResidualizedGenotypes.Rd | 5 ++- man/getResidualizedPhenotypes.Rd | 5 ++- tests/testthat/test_QtlDataset.R | 71 ++++++++++++++++++++++++++++++++ 6 files changed, 138 insertions(+), 4 deletions(-) diff --git a/R/AllGenerics.R b/R/AllGenerics.R index 55c6091f..1e92dab8 100644 --- a/R/AllGenerics.R +++ b/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", diff --git a/R/QtlDataset.R b/R/QtlDataset.R index 4d660845..27870385 100644 --- a/R/QtlDataset.R +++ b/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]] diff --git a/_pkgdown.yml b/_pkgdown.yml index 383dd0c3..89aac534 100644 --- a/_pkgdown.yml +++ b/_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" diff --git a/man/getResidualizedGenotypes.Rd b/man/getResidualizedGenotypes.Rd index a6e46b57..d95e3492 100644 --- a/man/getResidualizedGenotypes.Rd +++ b/man/getResidualizedGenotypes.Rd @@ -20,6 +20,7 @@ getResidualizedGenotypes(x, ...) residualizeGenotypeCovariates = TRUE, residualizePhenotypeCovariatesFromGenotypes = NULL, residualizeGenotypeCovariatesFromGenotypes = NULL, + covariateNaAction = c("impute", "drop"), ... ) } @@ -28,7 +29,9 @@ getResidualizedGenotypes(x, ...) \item{...}{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).} } \value{ A numeric matrix (samples x variants). diff --git a/man/getResidualizedPhenotypes.Rd b/man/getResidualizedPhenotypes.Rd index 3cdf022d..bc9a86f8 100644 --- a/man/getResidualizedPhenotypes.Rd +++ b/man/getResidualizedPhenotypes.Rd @@ -19,6 +19,7 @@ getResidualizedPhenotypes(x, ...) residualizePhenotypeCovariatesFromPhenotypes = NULL, residualizeGenotypeCovariatesFromPhenotypes = NULL, naAction = c("keep", "drop", "impute"), + covariateNaAction = c("impute", "drop"), outlierAction = c("keep", "drop"), outlierPvalThreshold = 0.001, ... @@ -30,7 +31,9 @@ getResidualizedPhenotypes(x, ...) \item{...}{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).} } \value{ A named list of numeric matrices keyed by context. diff --git a/tests/testthat/test_QtlDataset.R b/tests/testthat/test_QtlDataset.R index 8837f9b9..04e9b765 100644 --- a/tests/testthat/test_QtlDataset.R +++ b/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) # ===========================================================================