Repository navigation
Conversation
- Add detect_family() for automatic phenotype type detection
* Detects {0,1}, {1,2}, {-1,1} as binomial; else gaussian
* Handles NA values correctly; respects verbose flag
- Add recode_phenotype() for standardising phenotype coding
* Maps {1,2} -> {0,1} via subtraction
* Maps {-1,1} -> {0,1} via (y+1)/2
* Preserves NAs; stores original_coding in attributes
- Add validate_binary_phenotype() for pre-analysis validation
* Hard checks: no NAs, exactly 2 unique values
* Soft checks: warns if too few cases/controls (default >=10)
- Export all three functions via exports.R
- Add roxygen2 documentation for all three functions
- Add 27 unit tests in tests/testthat/test_phenotype_utilities.R
Why: Foundation for logistic regression support; enables auto-detection
of binary vs. continuous phenotypes and proper recoding to {0,1} standard.
Validation:
- All 27 phenotype utility tests pass (FAIL 0 | WARN 0 | SKIP 0 | PASS 27)
- Convert 1x1 matrix beta to numeric scalar in .fastlmm_core - Resolves test failure: 'FaSTLMM LL estimation is numerically stable' - Ensures consistency with test expectations and user API
Logistic regression (Phase 2): - Implement logistic_c<T>() template function in src/assoc.cpp - Newton-Raphson IRLS with configurable max_iter=25 and tol=1e-8 - Firth penalised-likelihood correction (enabled by default) for complete/quasi-separation stability - Wald test (chi-squared 1 df) for significance: p = 2*Phi(-|z|) - SVD pseudo-inverse fallback when X'WX is singular - Monomorphic SNP detection (var < 1e-10) skips marker gracefully - Type dispatch for char/short/int/double genotype BigMatrix types - OpenMP parallelisation over markers within each batch - Returns m x 3 matrix: [Effect, SE, p-value] with NA for failed markers - Auto-generated RcppExports.cpp and R/RcppExports.R updated Bug fix (cherry-pick from master): - Fix FaSTLMM beta return: as.numeric(beta) instead of raw 1x1 matrix - Resolves test_regression_gwas.R:131 FaSTLMM LL stability failure Why: Core computational engine for binary phenotype GWAS. Handles numerical stability (separation) and convergence monitoring. Validation: - devtools::load_all() compiles without errors (DONE rMVP) - logistic_c is.function: TRUE - All 4 genotype types dispatched correctly - All 45 regression_gwas tests pass (FAIL 0 | WARN 0 | SKIP 1)
- Implement MVP.Logistic() high-level R wrapper in R/MVP.Logistic.r
- Auto-detect binary phenotype via detect_family() when family='auto'
- Reject gaussian phenotypes with informative error directing to MVP.GLM()
- Validate binary phenotype via validate_binary_phenotype() before scanning
- Recode phenotype to standard {0,1} via recode_phenotype()
- Build null design matrix: intercept + optional covariates (CV)
* Drop constant CV columns automatically
* Validate no NAs in covariates
- Call C++ logistic_c() backend with full parameter passthrough:
* ind_idx, mrk_idx for subsetting individuals/markers
* firth, max_iter, tol for algorithm control
* cpu/threads for parallelisation
- Return m x 3 matrix [Effect, SE, p-value] consistent with MVP.GLM()
- Attach result attributes: family, method, n_cases, n_controls, firth
- Full roxygen2 documentation with parameter descriptions and example
Why: User-facing API for logistic regression GWAS. Abstracts all
C++ complexity and handles phenotype standardisation transparently.
Validation:
- Smoke test on extdata: 15x3 output, p-values in [0,1], SE >= 0
- family attribute = 'binomial', n_cases/n_controls correct
- All 45 regression_gwas tests still pass (FAIL 0 | WARN 0 | SKIP 1)
- Add family parameter to MVP() signature (default: 'auto') * 'auto' -> detect_family() auto-detects binary vs. continuous * 'binomial' -> forces MVP.Logistic() regardless of phenotype values * 'gaussian' -> forces MVP.GLM() regardless of phenotype values * Invalid value raises informative error - Add @PARAM family to roxygen2 documentation - Insert phenotype auto-detection block (after model-indicator setup) * Calls detect_family() only when GLM is in method list * Logs 'Binary phenotype detected' when binomial is resolved - Replace hard-coded MVP.GLM() call with family-aware dispatch: * 'binomial' -> MVP.Logistic() with full parameter passthrough * 'gaussian' -> MVP.GLM() (original behaviour, unchanged) * Shared post-processing (colnames, lambda, file output) unchanged Why: Seamless user experience; binary phenotypes are automatically analysed with logistic regression while maintaining full backward compatibility for continuous phenotypes. Validation: - Test 1 (continuous auto): routed to MVP.GLM(), 3-col result - Test 2 (binary auto): family=binomial, method=Logistic - Test 3 (explicit binomial): family=binomial confirmed - Test 4 (explicit gaussian): family attr NULL (MVP.GLM() used) - All 45 regression_gwas tests pass (FAIL 0 | WARN 0 | SKIP 1)
56 tests across 6 sections — FAIL 0 | WARN 0 | SKIP 0 | PASS 56:
Section 1 — Phenotype utility unit tests (8 tests):
- detect_family: binomial for {0,1}, gaussian for continuous
- recode_phenotype: {1,2}->{0,1} and {-1,1}->{0,1} with coding attr
- validate_binary_phenotype: accepts balanced, rejects continuous,
warns on too few cases/controls
Section 2 — MVP.Logistic() output structure (11 tests):
- Returns numeric matrix with 3 columns [Effect, SE, p-value]
- Correct row count matching marker count
- p-values in [0,1], SEs >= 0
- Works with {1,2} and {-1,1} coded phenotypes
- Works with covariates (CV argument)
- Rejects continuous phenotype with informative error
- Rejects family='gaussian' with informative error
Section 3 — MVP() integration dispatch (5 tests):
- Auto-detects binary -> binomial -> logistic
- Continuous -> gaussian -> GLM
- Explicit family='binomial' respected
- Explicit family='gaussian' overrides auto-detection
- 3-column output confirmed for logistic path
Section 4 — Numerical stability / edge cases (4 tests):
- Monomorphic markers return NA without crash
- Extreme imbalance (90:10) handled without error, valid structure
- firth=FALSE produces valid structure and valid p-values
- Near-perfect separation does not throw error (Firth protects)
Section 5 — Validation against R glm() (3 tests):
- Effect estimates correlate r > 0.90 with glm()
- SE estimates correlate r > 0.85 with glm()
- p-values correlate r > 0.85 with glm() (log-scale)
Section 6 — Output attributes (6 tests):
- family='binomial', method='Logistic' attributes
- n_cases and n_controls non-null numeric attributes
- n_cases + n_controls == n
- firth attribute reflects input parameter
- cbind compatibility with map data frame
No regressions: all 45 regression_gwas tests still pass.
…_regression_logistic.R The pattern '*.log' in R's dir() is treated as regex, where '.' matches any character - accidentally matching 'test_regression_logistic.R'. Changed to the proper anchored regex '\\.log$' to only match files ending in .log.
…(Phase 6)
- Create vignettes/logistic_regression_gwas.Rmd with a comprehensive tutorial
on using the new binary trait functionality (phenotype codings, Firth
correction, output interpretation, etc.)
- Add missing roxygen2 @PARAM tags in R/MVP.Utility.r
- Regenerate man pages and NAMESPACE via devtools::document() to expose
MVP.Logistic, detect_family, recode_phenotype, and validate_binary_phenotype
- Fix validate_binary_phenotype to respect verbose flag on error logging - Fix detect_family test to correctly capture stdout output - Suppress bigmemory downcast warning globally in test_regression_logistic.R Why: Keeps devtools::test() output clean and prevents confusing warnings/logs from polluting test results.
|
Thank you for your excellent work; I apologize for taking so long to notice this PR. A member of our team, Chen @cgl17686 , is conducting research related to logistic regression, so I have invited him to help handle this PR. |
|
Hi team, No worries at all about the timeline, I completely understand and appreciate you taking the time to look into this! Welcome @cgl17686! I am looking forward to your review and collaborating with you on this. Please feel free to ping me if you have any questions, the Firth correction implementation, or the warm-start optimizations. I am happy to discuss your research and make any necessary adjustments to ensure this aligns perfectly with the direction you want to take rMVP. Also, as a quick housekeeping item: could you please close the existing repository issues that are officially resolved by other PR (or let me know if I should link them to auto-close upon merge)? Thanks again, and looking forward to the feedback! |
|
Hi @zankrut20, after discussing this with @cgl17686, we really appreciate the work you’ve put into this, but we don’t think this PR fits rMVP’s main codebase in its current form. rMVP is centered on mixed-model association for large-scale GWAS, whereas this PR implements a per-marker fixed-effect Firth logistic regression rather than a mixed-model solution. Statistically, combining Firth estimates with unpenalized Wald tests may miscalibrate p-values, especially for rare variants, and the current implementation does not explicitly address highly unbalanced phenotype/case-control designs. Computationally, fitting each marker fully rather than using a score test limits scalability on large datasets. This makes it valuable for small-scale exact/penalized analyses, but not well aligned with rMVP’s scope and performance goals. So we would prefer not to merge it into the main rMVP codebase. That said, we’d be happy to discuss whether it could live as a standalone package/tool or be adapted for a more targeted use case. Thank you again for your contribution and understanding. |
|
Hi, |
2 participants
Overview
This Pull Request introduces full, seamless support for binary phenotypes in
rMVPvia Logistic Regression. It automatically detects whether a phenotype is binary or continuous and dispatches the computation to a highly optimized, numerically stable C++ backend.Key Features & Enhancements
1. Automatic Phenotype Detection
detect_family()utility handles auto-detection, andrecode_phenotype()automatically forces common binary schemas (e.g.,{1, 2}or{-1, 1}) into the standard{0, 1}format required by logistic regression.MVP.GLMas before.2. Firth-Corrected Newton-Raphson Optimization (C++)
logistic_cbackend insrc/assoc.cpp.NAand emitting anRcpp::warning()in R.3. "Warm-Start" Performance Optimizations
bigmemorybackends.4. Seamless API Integration
MVP.Logistic()as a dedicated wrapper.MVP()workflow. Users simply provide a binary trait, and the pipeline gracefully outputsEffect(Log Odds) and standard errors, just like FarmCPU or MLMM.Validation & Testing
validate_binary_phenotype,detect_family, andMVP.Logistic()pipelines.rMVPtest suite (152 tests) runs perfectly green. Existing OLS, FarmCPU, and MLMM logic remains entirely uncompromised..ped/.mapcohorts.Documentation
roxygen2(including newly exported utilities likerecode_phenotype).R CMD check(metadata is strictly preserved and non-ASCII characters have been thoroughly purged).Note to Reviewers:
All R-hub CI errors related to
.onLoad failed in loadNamespace() for 'RhpcBLASctl'andlibR.so: cannot open shared object fileare known R-hub infrastructure bugs with Clang/R-devel linux containers installing upstream dependencies, and are entirely unrelated to the structural integrity of this PR.Request: Please close those issues which are already resolved.