Methods

A suite-wide catalogue of implemented, experimental and planned capabilities.

Catalog authority. This page contains route-specific method descriptions and qualification notes. The current suite-wide LMM status is defined by the canonical model catalog, with evidence in model qualification. An exclusion below applies only to the workflow and version being described; it is not a global gsuite capability statement.

This catalogue maps development capabilities to their owning repository. Implemented means the capability is implemented in development code and exercised by relevant repository validation; it does not claim production-scale qualification. Experimental identifies a present but deliberately provisional integration surface. Planned identifies work with no current implementation. Availability through R is stated separately. Seven standalone genomic R packages are public source releases. Mixed-model and solver examples use the development integration; greml/gsolve R packaging is deferred, and a root gsuite R release is not currently planned.

The connected genomic example exercises representative routes across six modules using one 10K-individual, 50K-marker simulation. Its qualification record distinguishes execution checks and truth comparisons from untested capabilities.

The standalone command is gsuite run model.yml for both solving and REML. It calls the native gsuite_model_driver executable. The greml and gsolve library names identify their mixed-model and numerical-solver responsibilities; they are not separate user commands. Existing R interfaces are unchanged.

Capability overview

The native conditional ssSNPBLUP PCG route optionally replaces its full non-genotyped pedigree factor with an ancestral factor plus exact pedigree sweeps. This preserves the supplied covariance model and all individual coefficients. Supplied metafounder models additionally retain every correlated latent base in the small factor. Native dense-oracle checks and the ordinary saved 50K comparison are described in the factor-reduction study. The YAML frontend selects this opt-in route with execution.pedigree_factorization: ancestral_sweeps; existing R interfaces remain unchanged. The source bundle includes a complete GenomeDK 1M/16-thread model.

The experimental compact YAML frontend supports supplied-covariance multitrait pedigree, packed ssSNPBLUP and joint polynomial random-regression models with named input/effect mappings. It includes model inspection and labelled coefficient/curve predictions. Wide and long phenotype input preserves residual covariance within explicit measurement IDs, including missing traits and irregular visits. The long-input example (models/multitrait-single-step-long.yml) uses the same supplied-covariance model as the wide five-trait example. Preparation spools source rows and normalized occasions to temporary SQLite storage, and streams record prediction grids through prepared JSONL files. The optional execution.preparation.directory selects existing scratch space. Small storage tests check bounded retained Python row memory, cleanup and prediction-file integrity; pedigree/effect-level indexes and dense supplied matrices remain in memory. This is not industrial parser-memory qualification. Small independent dense-equation comparisons cover the adapters. Genetic channels can use separate incidence IDs within one shared relationship and joint covariance, including direct/maternal models. The complete example (models/direct-maternal-single-step.yml) reports direct and maternal coefficients separately and their mapped record contributions. This is supplied-covariance prediction; it adds no REML or PEV precision work. Selected ordinary packed genetic values additionally support exact PEV, prediction-error standard errors and reliability conditional on supplied covariances, with explicit solve-work caps and dense-inverse checks. Larger target sets can use variance-reduced stochastic PEV or prior-minus-posterior information probing with explicit probe/solve budgets, MCSE and separate prediction/uncertainty acceptance. A local 50K study measures overhead and precision; it does not establish industrial uncertainty throughput or calibration. The matched information-probing comparison uses the same posterior query budget and reports additional prior-product work. Stochastic PEV also supports trait-isolated residual probes at fixed query cost, with MCSE over complete independent groups. This removes cross-trait noise but reduces the independent sample count; its benefit is model-dependent. The fixed-budget 50K comparison records grouped MCSE, exact-target errors, runtime and memory. At 32 posterior solves, trait isolation worsened median MCSE and the eight-target exact error on this fixture; unstructured probing remains the default. Neither plan met the population-wide precision requirement at that budget. REML, broader uncertainty targets and industrial parser-memory qualification remain outside this frontend. The interface additionally supports supplied dense relationship covariance, fixed interactions, composite random-effect groups, Legendre/Chebyshev bases, group-specific trait residual covariance and univariate AR1 or continuous-time exponential residuals on the compiled route. Packed ssSNPBLUP also supports univariate AR1/exponential covariance through exact sparse Markov innovations; actual observed gaps, missing visits and singleton series retain their supplied covariance model. The AR1 (models/single-step-ar1.yml) and exponential (models/single-step-exponential.yml) examples include complete PLINK/YAML workflows. Grouped residual covariance also uses the existing ssSNPBLUP prepared-record path. Independent small dense references check the new compiled routes; this does not qualify their large-data memory use or new combinations of covariance structures. AR1/exponential residuals can also use a supplied trait covariance times a common time correlation. Compiled and packed prediction routes prepare exact conditional innovations over observed traits, preserving marginal covariance with intermittent missing measurements. The multivariate AR1 example (models/multitrait-single-step-ar1.yml) includes PLINK inputs and covariance matrices. Small dense references check wide/long data, irregular gaps and both packed pedigree factors/preconditioners. No dense time-series covariance is prepared; sparse filter design history may grow within incomplete series. Serial REML and industrial multivariate preparation throughput remain outside this qualification. Random-effect channel covariance also accepts supplied factor loadings and positive specific variances, giving K = Lambda Lambda' + diag(psi). Compiled FA1/FA2 models reuse native FA parameterization; packed ssSNPBLUP uses the same validated implied channel covariance. Environment-labelled long input represents genetic covariance across environments, with separate fixed coefficients per declared channel. Small independent dense references check coefficients and equivalent unstructured models, including partial genotyping and both pedigree factors/preconditioners. The multi-environment example (models/multienvironment-factor-analytic.yml) supplies all inputs. This adds prediction with supplied FA parameters; YAML FA estimation, uncertainty and industrial scaling remain unqualified. The simulator and user-data commands share this one YAML schema; numerical controls are specified in its execution block.

Supplied-genomic ssGBLUP is available through compiled PCG, with either a dense effective G covariance or explicit APY core/cross covariance blocks and marginal variances. Native greml retains exact marginal pedigree precision and applies implicit H; APY is explicitly a genomic covariance approximation. Small dense references and equivalent ssSNPBLUP models check predictions, including mixed genotyping and covariance-base settings. YAML also connects raw PLINK data to native streamed APY preparation, with ordered core IDs, explicit covariance scaling/ridge, bounded decoded tiles and file-backed cross products. Independent small dense equations check missing genotypes, ID mapping, APY predictions and the all-core exact boundary. Preparation avoids full G or decoded-panel storage, but retains quadratic core memory and O(noncore*core) disk. Compiled prediction defaults to omitting unused H covariance diagonals; exact precision diagonals/determinants and the preconditioner are unchanged. execution.relationship_metadata: complete restores exact covariance metadata. Industrial APY throughput remains unqualified; packed ssSNPBLUP retains its separate HPC evidence.

The YAML interface also exposes supplied metafounder Gamma with explicit unknown-parent assignments, or explicit independent random genetic group means. They support pedigree, supplied-genomic ssGBLUP/APY and packed ssSNPBLUP models. Metafounder base covariance and additive group-mean contributions remain distinct. Small independent dense equations check coefficients and predictions, including latent levels and individual a + F gamma reporting. Compiled fixed inherited group contrasts, explicit phantom group-parent absorption and combined metafounder/group models are also exposed. Packed Gamma pseudo-EM requires an explicit absolute marker base and returns conditional predictions at the estimated Gamma. The compiled group adapter materializes ancestry proportions; the packed adapter instead applies native F/F’ sweeps to cached original records and prepares exact group diagonals once with operator probes. Packed groups require matching genetic channels and a final group-mean effect, and support full or ancestral Qnn. Supplied metafounders also support ancestral reduction with their complete correlated base block retained; Gamma and the original covariance model are unchanged. The metafounder reduction study records exact references and the matched 50K comparison. Small dense references check both preconditioners and consistent individual reporting. Their industrial memory/runtime qualification remains separate. Existing R interfaces are unchanged.

The packed pedigree-base study records the focused dense/analytic checks and the 50K fixture measurements. Metafounder factors use an exact diagonal congruence for balance; native sparse verification uses compensated residuals and bounded retained-factor correction when needed. Original model residual thresholds and covariance parameters are unchanged. These workstation measurements do not establish 1M/HPC performance.

The optional standalone HPC ssSNPBLUP prediction frontend accepts supplied variance–covariance matrices, with either independent genomic residual covariance or a pedigree A22 blend. The source-only blend adapter has a focused dense-model oracle for both packed PCG routes; its representative five-trait workload has no industrial throughput or memory qualification yet. It performs fixed-covariance prediction and adds no REML or R selector. See the profile settings and limits.

Native greml has a separate opt-in two-component genomic REML route: generic PCG trace/AI identities and analytical profiling of the common variance scale for V = g K + e I. It supports a zero-eligible genetic component and positive residual, with explicit stationarity and numerical/stochastic precision gates. This development C++ interface has no new R selector; fixed absolute variances, general multi-component profiling and large-panel performance qualification remain outside its scope.

Scalar individual BED BayesC/BayesR in standalone gbayes expose opt-in scheduled marker updates, exact scheduler-state continuation, and realized genetic variance decomposition (Vg = VLE + VLD). Scheduling approximates full-sweep sampling; negative LD contributions are valid. Per-chain full-sweep variance and mixture traces preserve burn-in separately from post-burnin means. See the controls and interpretation. Their fit$posterior table includes VG, VLE, VLD and draw-wise heritability, with native uncertainty and convergence summaries.

Reviewed against the development interfaces on 30 September 2026. These are available workflows, not claims of full upstream-software parity or broad statistical calibration. Each guide records its supported inputs and limits. See the estimate/solve contract for the mixed-model entry points. | Analysis or resource | Native owner | Current R access and guide | | — | — | — | | Pedigree, multivariate and longitudinal mixed models | greml, gsolve | greml() for variance estimation, gsolve() for fixed-covariance solving, gfit() and persistent model interfaces; pedigree tutorial, LMM catalog | | Genotypes, marker QC, sparse LD, independent blocks, LD scores and regional eigenvalues/eigenpairs | gmat, gbits, gsolve | gbase::gprep() prepares and augments Glist resources; task = "eigen" prepares reusable full regional BED eigenvalues through the smaller sample/marker Gram matrix; validate_Glist() checks their lightweight descriptor and verify_ld() performs explicit full integrity verification; define_ld_blocks() verifies a block partition once; gbase::getChromosomeSets() extracts disjoint chromosome marker sets from Glist metadata; gbase::gfilter(), ld_scores(), ld_region() and ld_eigen() consume the same Glist; genotype preparation, LD resources | | Summary-statistics alignment, LD consistency checking and Z-score imputation | gbase integration over gbits, gmat and gsolve | checkStat(Glist, stat, ldcheck=TRUE, impute=TRUE) handles data frames, trait matrices and gact Z/N lists; impute(Glist, stat) predicts missing scores alone. Complete regional BED LD, ancestry sample groups, shared trait factorizations, approximate DENTIST-style diagnostics and quality thresholds keep accepted scores intact. Explicit mixtures assume independent standardized meta-analysis coefficients constant across markers; missing effect estimates and sample sizes are not inferred. Preparation and limits | | LDSC, GNOVA, SumHer, HESS/rho-HESS, SUPERGNOVA, LAVA and HDL-L | gcorr | gcorr(Glist, stat, method=...) fits ordinary LDSC, GNOVA or an illustrative equal-marker SumHer baseline from stored scores and gact Z/N matrices, or joint HESS from stored independent eigen-region definitions and reference samples with compact moments returned for reuse; ldsc(z=..., n=..., ldscores=...) accepts direct named arrays. ldsc(stat, Glist, ...) also provides the checked single- or multi-trait Glist workflow, including standard stratified and qgg-style grouped marker-set models with compact integer membership for disjoint groups; gcorr_set_scores() supplies the explicit compact input; LDSC, partitioned LDSC, GNOVA and SumHer pass numeric arrays and integer block membership to native code, retaining marker names in R and compact fit summaries; explicit focused functions remain selectable directly or through gcorr(method=...). Glist regional preparation reuses one sparse-LD reader and verifies one immutable spectrum per region for independent/joint HESS, SUPERGNOVA, HDL-L, and LAVA marginal, binary, simulation, partial-correlation and regression workflows. Task and uncertainty coverage. SUPERGNOVA exposes covariance only; HDL-L exposes two continuous traits per fit, with inference under evaluation | | SEM, partial-correlation networks and fixed-graph GGM | gcorr | gcorr_covariance(), gcorr_sem(), gcorr_ggm(), gcorr_ggm_network() and native test functions. Grouped interfaces extract GNOVA categories and SumHer components/categories with full cross-group jackknife covariance for joint SEM/GGM contrasts. See covariance models. | | GSMR, GSMR2, IVW, Egger and weighted-median MR | gcorr | gcorr_prepare_mr() plus focused method functions and thin gcorr(..., method=) dispatch consume aligned ordinary summary-statistics data frames, explicit effect units and signed LD. See the MR workflow. | | BayesC/BayesR from individual BED or aligned summaries, annotation-informed block-eigen regression, several traits or independent studies | gbayes | Learned scalar and multivariate summary/BED BayesC/BayesR, verified block-eigen scalar BayesC/BayesR through mcmc-eigen, plus single-trait sbayesrc over verified Glist LD blocks with retained eigensystems and component annotations. Shared-LD or independent sampling, study-specific references, covariance learning, posterior capture and continuation, genetic covariance, reporting partitions, selected-individual score uncertainty and four focused summary engines are available; Bayesian regression. | | Regional Bayesian mapping and posterior sets | gbayes, gmap | gmap::gmap() fits regional multi-effect, BayesC or BayesR models using full BED LD for dense MCMC or explicit sparse-LD eigen preparation; familiar iteration/prior arguments, regional post/conv/stat/cs outputs and standalone getCredibleSets() are available; gmap_bed() interprets a saved individual-BED posterior without refitting. gmap_sets() presents labelled set/member tables and gmap_bfdr() summarizes PIP-threshold Bayesian FDR. Mapping workflow. Marginal-PIP groups and component sets have different meanings | | ABF and component-pair colocalisation | gmap | gmap_coloc() provides the single-causal ABF route; gmap_coloc_components() compares saved multi-effect components without refitting. Colocalisation has method-specific assumptions | | Linear/logistic GWAS and quantitative-trait LOCO association | glma | glma::glma() provides ordinary, conditional, batched linear, logistic, infinitesimal LOCO and prior-selected LOCO scans over gbase Glist inputs. Binary logistic tests do not adjust for relatedness | | Genomic scoring and posterior prediction | gscore, model-fitting packages | gscore(weights, Glist) applies portable allele-labelled weights and can accumulate overlapping, cross-file marker-set scores from Glist-resolved marker names. clump_threshold(stat, Glist) constructs conventional C+T weights from reported association effects and prepared sparse LD. glma, gbayes and other model owners fit their statistical effects. Selected-individual score uncertainty is available through gbayes::gbayes(); guide | | Genotype, pedigree and phenotype simulation | gsim | Separate public R package; gsim documentation | | Marker-set association evidence | glma, gsea, gsolve | gsea::vegas() re-exports glma::vegas(). It consumes ordinary checked summary-statistics data frames and Glist sparse LD and provides LD-aware weighted quadratic tests. prepare_vegas() computes eigenvalue-only spectra in bounded parallel batches and returns a serializable object reusable across traits and sessions; gsea_input() validates the gene-labelled enrichment handoff. ACAT remains native-only. Default tails are approximate except for equal positive eigenvalues; the controlled positive-series evaluator reports accuracy and term-budget failures | | Gene-set, genomic-feature and gene-prioritization analysis | gsea | gsea() provides ORA, preranked GSEA, competitive GLS, bounded MAGMA and optional scalar/multitrait Bayesian pathway regression. Familiar magma() arguments adapt the existing competitive/Bayesian engines. pops() fits chromosome-held-out ridge, BayesC or BayesR regression and returns feature-based gene rankings with compact per-fold posteriors; marginal selection is optional for Bayesian fits. Explicit universe/set preparation, native aggregate/block resampling and optional BH/BY adjustment remain available. Pathway and gene-ranking workflow |

The native MAGMA comparison covers SNP-wise mean gene evidence and competitive pathway analysis with MAGMA 1.10 preprocessing. A small pinned comparison checked 40 genes and two pathways. This is not complete MAGMA equivalence: other gene models, rare-variant preprocessing, missing-data handling and singular-matrix pruning remain outside the route. Its controlled tail evaluator may return unavailable for difficult spectra or tails rather than use a silent fallback. Scalar and joint multitrait Bayesian pathway analysis is exposed through gsea(), reusing gbayes for sampling and parameter learning. PoPS-style gene ranking is exposed separately through gsea::pops(). It repeats marginal feature selection and ridge fitting, or native Bayesian regression with optional training-only selection, after excluding each target chromosome, accepts explicit blockwise gene covariance, and returns ranked predictions plus fold and feature diagnostics. It is a focused PoPS-style model rather than a claim of output equivalence to every upstream PoPS preprocessing and locus-reporting step.

Selected native operations remain outside the R bindings, including gscore pseudo-summary generation. The standalone gbayes::gbayes_single_effect() interface exposes fixed-variance single-effect inference from ordinary summary statistics. Checks cover independent small calculations, pinned ACAT comparison and a connected native example, not broad GWAS calibration. glma::vegas() implements native LD-aware quadratic gene evidence from aligned GWAS Z statistics. prepare_vegas() reuses one validated sparse-LD reader and prepares eigenvalue-only gene spectra in bounded parallel batches. Its ordinary R result can be saved and reused across traits, with explicit tail status during evaluation. It also consumes the full regional BED spectra stored by gprep(task="eigen") without repeating decomposition, for constant positive within-set weights. VEGAS selects checked statistic rows by marker ID; extra rows and another row order are accepted, while all prepared markers and their matching reference alleles remain required. The sparse-LD route retains nonuniform weighting. Existing tail defaults and accuracy limits are unchanged. The gsea package re-exports both interfaces and gsea_input() validates the gene-statistic handoff to gsea() without automatic harmonization or inferred cross-gene dependence. Methods that require gene covariance need it supplied separately. ACAT and other gene-evidence combinations remain native-only. Competitive GLS uses supplied gene-covariance blocks; its Student-t inference is a working approximation for non-Gaussian gene scores. It does not claim full MAGMA or Pascal equivalence. Python interfaces remain deferred. Multivariate sampling covariance can be learned only for the shared-LD route with an explicit response cross-product; arbitrary sampling overlap across different LD references remains outside the interface. See the current qualification limits before interpreting capability as validated use.

The gbayes::gbayes() choices include point_normal, continuous_shrinkage, variational and mixture_variational. These retain method-specific controls/results, rather than inheriting all BayesC/BayesR features. gbase::checkStat() provides the shared ordinary summary-statistics preparation step. See the summary interface guide.

Pedigree formula workflows

The pedigree tutorial demonstrates the data-frame interface on 10,000 simulated animals. It reuses existing numerical engines:

Formula model	Scope
Pedigree animal model	One or multiple traits; missing responses omitted trait by trait
Genetic random regression	One response, one pedigree effect and homogeneous residuals
Independent components	Multiple pedigree/identity intercepts, including permanent environment and common groups
Correlated channels	One response with named channels under a shared kernel, such as direct/maternal effects
Structured residuals	Grouped variances, integer-gap AR(1), or continuous-time exponential correlation

Short and explicit forms share normalization. formula and formulas are aliases; starts do not select the engine. REML initializes omitted starts, while supplied values remain authoritative. Fixed solves require complete covariance values. The tutorial gives the route-specific initialization rules.

Distinct descriptors are independent. Correlated channels must be declared explicitly and have known IDs for observed responses. Correlated residual blocks require unique positions within one variance group. Multi-component multitrait/random regression, multi-response residual descriptors, nuggets and unknown-time imputation remain outside this adapter. Covariance boundaries and unsupported combinations are not silently repaired. These examples do not extend workload or platform qualification.

Association mapping

Association mapping belongs in native glma and the Glist-based glma() R interface; fine-mapping and colocalisation remain in gmap. The current quantitative-trait methods are exact linear regression, Gaussian/ridge-only infinitesimal LOCO offset testing, and prior-selected LOCO offset testing with automatic structure-dependent scaling or a manual override. The LOCO methods can fit on a selected subset of Glist markers while testing every marker in that Glist; the subset must cover every tested chromosome. Association scans also accept a phenotype matrix with named traits and return one association result per trait. Exact linear scans share genotype reads and dosage decoding. LOCO traits share immutable genotype preparation but each retains its own variance fit, prior selection and calibration. Both LOCO methods use the joint prediction engine. Infinitesimal LOCO fixes one Gaussian/ridge prior and equal standardized-marker variances; prior-selected LOCO compares ridge/elastic priors. Blocked RHE initializes a shared three-point Monte Carlo heritability grid. The global variance pair remains fixed for chromosome-excluded refits. Diagnostics identify this approximate estimator, its grid fallback and fitting phases. The slow per-chromosome variance search and its controls are retired from the standalone interface. Portable matrix products support operation-local OpenMP threads. See the LOCO variance procedure. Exact linear scans can include selected marker genotypes as fixed covariates through condition_on; this is a complete-case conditional test and does not reuse the LOCO null model. Binary traits can use the distinct logistic method for ordinary covariate-adjusted per-marker log-odds Wald tests. Sparse case/control counts, separation and non-convergence have explicit statuses. This baseline does not adjust for relatedness or claim LOCO binary calibration. glma provides linear, logistic, infinitesimal LOCO and prior-selected LOCO functions plus the glma(..., method=) dispatcher. vegas() consumes validated gbase summaries and Glist sparse LD. These methods are separate from gmap_stat() summary-statistic preparation.

Marker sets can also define separate genetic background components for either LOCO method. Named disjoint groups plus an automatic unassigned background are fitted jointly using greml’s multi-component RHE variance profile and the joint engine’s scalar Monte Carlo scale refinement. Prepared storage is required. Component estimates and chromosome exclusions are returned in fit$components and fit$loco_components. See the component model and its approximation limits.

Computational building blocks

Library	Method or capability	Status	Local documentation
glma	Glist-based exact linear and binary logistic association, infinitesimal LOCO, and prior-selected LOCO with conditional ridge-reference scaling	Broad mixed-model null-tail validation remains open	Association scope; installed help("glma", package="glma")
gbits	Strict PLINK BED access, decoding, ordered selections, marker statistics, and affine missing-aware transformations	Implemented	gbits methods
gbits	Direct and prepared \(Xv\), \(X^{\mathsf T}v\), \(XB\), \(X^{\mathsf T}Y\); source-independent binary/ternary packed matrices	Implemented	gbits methods
gbits	Weighted genotype covariance, normal-equation, diagonal, and trace operators	Implemented	gbits methods
gbits	gbits::gsolve_adapter genomic covariance and packed rectangular design operators	Experimental	gbits methods
gmat	BedMatrix, ordered BedPanel, in-memory and disk-backed sparse-LD CSR resources and algebra	Implemented	gmat methods
gmat	BED/panel sparse-LD construction; ordinary, finite-sample-corrected, and annotation-specific LD scores	Implemented	gmat methods
gsolve	Dense Cholesky and sparse LDL’ SPD factorization and solves	Implemented	gsolve methods
gsolve	Matrix-free PCG, independent and true block multiple-RHS orchestration, SolvePlan, preconditioners, and execution diagnostics	Implemented	gsolve methods
gsolve	Square/rectangular structured operators, reproducible Rademacher probes, and stochastic Lanczos log determinants	Implemented	gsolve methods
greml	Exact dense single-trait ML/REML evaluation and reference or safeguarded constrained optimization	Implemented	greml methods
greml	Matrix-free stochastic ML/REML, adaptive precision, sparse MME, and hybrid covariance routes	Implemented	greml methods
gsuite	Eligible affine spectral reuse, purpose-specific staged stochastic REML, and opt-in determinant-deferred packed AI-REML	Experimental; explicit controls	persistent R interface
greml/gsuite	Packed coefficient/sufficient-statistic reuse for single-trait and complete-record multivariate stochastic score/AI, with bounded prefix-aware probe cache and fused shared-source reconstruction	Experimental execution optimization	persistent R interface
greml/gsuite	Balanced structured longitudinal RHE to determinant-deferred AI-REML to latent-coefficient and trajectory BLUP	Experimental complete-record workflow	persistent R interface
greml/gsuite	Sparse-pedigree determinant-deferred AI-REML and retained-state all-animal BLUP, with coordinate or CHOLMOD selected-inverse exact direct traces	Experimental positive-interior workflow	Sparse-direct pedigree methods
greml/gsuite	Positive-interior multitrait pedigree sparse-direct REML with missing-trait residual blocks, bounded coordinate or CHOLMOD selected-inverse exact traces, exact PARDISO/CHOLMOD determinants, and retained BLUP	Experimental moderate-problem workflow	multitrait pedigree workflow
greml/gsuite	Irregular scalar longitudinal pedigree PCG/PARDISO/CHOLMOD MME, positive-interior random-regression REML, grouped diagonal residuals, selected-inverse exact traces through CHOLMOD, and retained coefficient/trajectory BLUP	Experimental; direct optimization is bounded moderate-problem work	longitudinal pedigree workflow
greml	Unbalanced multivariate and longitudinal/random-regression covariance inference and optimization	Implemented	greml methods
greml	Unified ModelSpec → CompiledModel → ModelFit facade, prediction, contrasts, and uncertainty	Implemented	greml methods
gcorr	Ordinary, standard stratified and grouped single-trait LD-score regression, enrichment, jackknife, and diagnostics	Implemented	gcorr methods
gcorr	Bivariate and serial multivariate LDSC with genetic covariance/correlation and paired jackknife uncertainty	Implemented	gcorr methods
gcorr	Partitioned bivariate LDSC with annotation covariance coefficients, total/category genetic covariance and correlation, and paired jackknife; annotation-specific two-update weights	Implemented	gcorr design
gcorr	Native SNP summary-statistic preparation: explicit allele orientation, missingness, variable sample sizes and a common all-trait marker mask	Implemented	gcorr native workflow
gcorr	Ordinary multivariate labelled joint sampling covariance from paired block deletes	Implemented	gcorr native workflow
gcorr	Partitioned multivariate LDSC with an all-trait marker mask, annotation and total/category matrices, and labelled joint sampling covariance from coordinated block deletions	Implemented	gcorr design
gcorr	GNOVA, SumHer, HESS/rho-HESS, SUPERGNOVA, LAVA, SEM/GGM and MR	Focused interfaces and gcorr(method=...) dispatch	Capability overview and linked method contracts
gscore	Tiled multi-score BED application with standardized or raw mean-imputed coding and operation-local OpenMP	Implemented	gscore methods
gscore	Portable C+T weight construction and allele-labelled multi-score application to gbase resources	clump_threshold(statistics, Glist) uses validated sparse LD; gscore(weights, Glist, sets=) supports exact/swapped allele alignment, selected-sample order, sparse marker-set accumulation and standardized or raw mean-imputed coding. Statistical model fitting remains outside this R package	scoring workflow
gsuite	Linked version/capability reporting and gsolve::SolvePlan SPD solving from R	Experimental	R package overview
gsuite	Exact dense and prepared matrix-free BED ML/REML bindings through the greml model facade	Experimental	R package overview
gsuite	Persistent prepared-BED, compiled single-trait model, fit, contrast, Wald, uncertainty, and optimizer-diagnostic interface	Experimental	persistent R interface
gsuite	Independent-group, ordered-pedigree, and general sparse-Z terms through matrix-free covariance, sparse MME, and structured-plus-genomic hybrid routes	Experimental	persistent R interface
gsuite	Explicit unbalanced multivariate indexed covariance, longitudinal random regression, temporal residual structures, complete predictions, and selected-target uncertainty	Experimental	persistent R interface
gsuite	Canonical gsolve/greml route, scaling, R-interface, and bounded qualification overview	Experimental	gsolve/greml capabilities
gsuite	Serializable BED/BIM/FAM resources and ordered panel selection through gmat::BedMatrix and gmat::BedPanel	Experimental	gmat R integration
gbase	Serializable qgg-familiar genotype preparation with exact packed BED summaries, QC aliases, marker-set helpers and stored LD scores; qualified gmat/gbits build/load guards	Experimental	genotype preparation
gbase	Summary-only marker QC with exact HWE, INFO identity alignment, allele, missingness, and assembly-aware MHC rules	Task-neutral interface	genotype filtering
gbase	gbase::gprep(task = "sparseld") orchestration returning an augmented Glist with selected physical BED mappings, transactional gmat CSR, aligned AF, ordinary LD scores, lightweight descriptor validation and explicit full verify_ld() integrity checks	Experimental	sparse-LD preparation
gmap	Regional BayesC/BayesR and marginal/component posterior sets via unified native gmap built with Bayesian support	Broader native controls remain outside the focused R interface	Regional posterior interpretation
gbayes	Single-/multi-effect fitting with optional prior-variance learning	Implemented	Single-/multi-effect inference; Bayesian workflow
gbayes	Scalar BayesC/BayesR plus point-normal, continuous-shrinkage and variational summaries; multivariate fits, independent studies, annotation and parameter learning, posterior genetic summaries	Summary/BED fitting, independent study references, annotations, covariance learning, posterior capture and continuation, genetic-covariance reporting, partitions and selected-individual uncertainty are implemented.	Bayesian workflow

The connected gcorr validation checks genotype-derived LD scores and fits against dense references. A shared unit floor on marginal working means keeps all prescribed draws finite on the fixture. The floor changes variance weights without clamping raw estimates, and activation is reported in native diagnostics. See the weighting contract. The subsequent statistical calibration study completes all fits and passes parameter-recovery screens, but flags undercoverage for 15 primary quantities. Uncertainty qualification remains incomplete. An independent trait-1 QR diagnostic agrees with the native fixed-weight procedure; re-estimating weights during block deletion improves coverage but does not fully resolve it. These capabilities have not been broadly qualified on real GWAS. The coordinated reliability study finds better coverage with known weights and analytical SEs, supporting an uncertainty-procedure limitation in the tested setting. Matched full LD gives stable estimates in the fixture; short windows and mismatched references cause bias. Current normal/jackknife intervals remain experimental; broad confidence-interval qualification is still open.

Detailed equations, layouts, ownership, and limitations remain in each repository’s local catalogue and linked design documents. The actual local dataset campaign is summarized separately in qgdata qualification.

Marker-space and RHE methods

With marker weights D, normalization s, genomic variance V_g, and residual variance V_e, gsuite derives \(\Lambda = (V_e s/V_g)D^{-1}\). Explicit marker_gs constructs gbits’ existing fused-coordinate PackedGsPreconditioner and uses it inside every covariance PCG solve in the REML evaluation and optimizer. It is therefore reported as marker_gs_preconditioned_pcg, not as a direct marker inverse solve. The fused coordinate cache is prepared once per fit; a changed parameter state builds a bounded penalty/preconditioner state, while identical common-precision reevaluation may reuse the exact immutable pairing.

marker_preconditioner_sweeps is the exact symmetric sweep count per preconditioner application. Diagnostics separate cache construction, parameter-state preparations, PCG applications, processed RHS columns, sweeps, marker visits, packed decode passes, factor lifecycle, and memory. Fused coordinate application uses no block-Cholesky factorization. The surrounding covariance route continues to report response/preparation/information PCG iterations and true residuals together with operator and preconditioner work. Multiple RHS, including block PCG, are supported by gbits; concurrent outer applications are rejected before fitting because the packed preconditioner is nonconcurrent.

The total genomic prediction is recovered as \(y-F\beta-V_e Py\), so no generic \(XB\) pass or final marker solve is performed for that result. A precise final marker solve is reserved for explicitly requested component-specific marker effects or validation.

initializer="rhe" asks greml to build projected response moments and common- probe trace moments, solve the unregularized small system, retain raw estimates, and make every zero/positive boundary adjustment explicit. RHE is a start, not an ML/REML estimate; the existing safeguarded AI-REML route remains responsible for the reported likelihood and convergence. RHE diagnostics cover initialization only. Evaluation at the RHE starting point, subsequent optimizer work, and optional final marker prediction are separate stages; whole-fit elapsed time is not labeled as RHE initialization.

genomic_solver="marker_batch" is distinct from the preconditioner. For every parameter state it derives one penalty vector per ordered genomic component and calls gbits 0.16 once with response plus the complete deterministic probe prefix. AI component right-hand sides are submitted in one additional batch. Each system owns independent fixed effects, marker coefficients, residual, convergence and warm state; all components in that system share the same observation residual. One decoded block is read by disjoint system-column owners, and exact forward/backward marker order is retained per system.

Response/common-probe and component-Py batches retain independent warm states. Existing ordered columns remain warm through compatible penalty changes and staged precision tightening; deterministic prefix extensions initialize only new columns at zero. An active component-set change invalidates the incompatible marker state. When all genetic variances are zero, gsuite handles the residual-plus-fixed projection analytically and never constructs an invalid marker penalty.

For precisely one genomic component plus identity residual under REML, the explicitly selected experimental option marker_score_ai="identity_optimized" uses the identities

\[ \begin{aligned} V_g\operatorname{tr}(PK)+V_e\operatorname{tr}(P)&=n-p,\\ V_g KPy+V_e Py&=y-X\widehat{\beta},\\ V_g P(KPy)+V_e P(Py)&=Py. \end{aligned} \]

to derive the genomic trace and score vector from the projected common probes and response, then projects only Py in the information batch and derives \(P(K Py)\). The mathematical score and full 2-by-2 AI definition do not change. marker_score_ai="direct", the public default, retains the component covariance applications and two-column information batch. The optimized route requires \(V_g>64\sqrt{\epsilon_{\mathrm{double}}}\max(\lvert V_g\rvert,\lvert V_e\rvert)\); near or at zero, for ML, or for any model other than exactly genomic-plus-identity, it falls back to the direct route. This is independent of the opt-in determinant- deferred acceptance policy.

With the direct statistical policy, the independent execution control marker_coefficient_reuse defaults to "automatic"; "off" forces the original component applications and "required" fails if a safely eligible packed component cannot supply all coefficient-derived quantities. This does not select another estimator. The packed MME equations use \(\Lambda_{cj} = V_e s_c/(V_c w_{cj})\) and therefore return \(a_c(b) = V_c D_c W_s^{\mathsf T}Pb/s_c\). gsuite caches \(W_s^{\mathsf T}Z\) once per prepared-source and probe identity, then evaluates the original direct probe realization \(a_c(z)^{\mathsf T}W_s^{\mathsf T}z/V_c\) and response quadratic \(s_c a_c(y)^{\mathsf T}D_c^{-1}a_c(y)/V_c^2\).

The cache is double precision and bounded by marker_coefficient_cache_memory_limit (512 MiB by default). Its identity includes the prepared packed storage, ordered sample/marker selection, transformation and missing policy, common-probe stream, dimensions, and actual probe columns. A compatible prefix hit reuses old columns and transposes only the extension; independent verification probes and incompatible sources form separate entries. Exceeding the budget falls back in automatic mode and fails in required mode. At zero or below the scale-aware guard \(64\sqrt{\epsilon}\,\texttt{parameter\_scale}\left(1+\frac{\texttt{achieved\_normal\_residual}}{\texttt{requested\_tolerance}}\right)\) the direct boundary component path is retained, so no division by a small variance occurs.

For one packed component plus identity, \(KPy\) uses the exact response residual identity. Multiple weight-defined components sharing a source concatenate their response marker coefficients and use one multi-RHS packed xmat, so the source is decoded once. Generic components retain direct operator application. The cache transpose and fused reconstruction use the existing operation-local packed thread controls and no nested OpenMP. Diagnostics report requested and resolved policy/components, fallback reasons, cache hits/misses/columns/bytes, reused/appended columns, reduction and reconstruction time, avoided and direct applications, information width, and reconciled marker visits/decode passes. Likelihood and determinant computation, including determinant-deferred’s one final likelihood-only verification, are unchanged. Longitudinal AI-REML coefficient reuse remains future work.

Fixed-parameter evaluation can omit likelihood and SLQ explicitly. The default optimizer remains likelihood-monotone. Opt-in acceptance_policy="determinant_deferred" uses curvature-scaled projected score and Newton decrement for backtracking, with fixed and active-boundary components removed from the free AI solve. Provisional convergence additionally requires projected KKT score, parameter change, solver residual, score SE, stable active set, a larger final probe prefix, and independent verification probes; a complete deterministic trace basis may bypass fresh probes. Iterative accepted/rejected/verification states do no determinant work. One generic covariance-PCG/SLQ likelihood-only verification runs after provisional convergence and is merged without repeating packed score/AI batches.

The policy is supported for single-trait, complete-record multivariate, balanced structured longitudinal, and general observed-record longitudinal matrix-free workflows with the safeguarded nonadaptive optimizer. Sparse, hybrid, legacy indexed/sparse longitudinal, adaptive, and other facade routes reject it. Work diagnostics separate all score/AI, likelihood, verification, SLQ, batch, thread, sweep, warm-system, decode, residual, and timing categories. Score/AI trials can cost more than likelihood-only rejected trials, so speedup is workload-dependent. Balanced-longitudinal BLUP requires REML convergence; an unavailable BLUP reports that failure explicitly and does not establish a converged boundary solution. Structured Q defaults to "positive" (unchanged log-Cholesky SPD behavior). Balanced "may_be_zero" blocks use symmetric covariance coordinates with PSD projection and AI quadratic proposals, accepted by likelihood backtracking. Convergence requires the AI-scaled projected-score mapping to satisfy the existing score tolerance; projection or failed backtracking alone never suffices. At zero, a positive score eigen-direction can be released. Complete explicit probes and full-dimensional Lanczos depth bound this initial deterministic route. RHE keeps its raw estimates and maps PSD starts without a positive eigenvalue floor. Residuals stay strictly positive; fixed and observed-record PSD blocks remain unsupported. Likelihood and BLUP apply Q without inverting it.

Packed-call diagnostics label response/common-probe, component-information, identity-information, and likelihood-inverse work and report RHS width, requested/resolved threads, marker visits, decode passes, decode time, compute time, and wall solve time. The compatibility packed_decode_passes field keeps its prior marker-batch meaning. batched_total_packed_decode_passes includes likelihood-inverse batches, covariance_packed_decode_passes includes packed covariance products used by score components, PCG, and SLQ, and the explicitly named total_packed_decode_passes is their sum; marker-visit and timing totals are exposed alongside them.

Multivariate RHE-to-BLUP workflow

gs_multivariate_rhe_blup() is a single bounded workflow: it estimates named component trait covariances and the residual trait covariance once with matrix-free randomized HE, then passes the explicitly PSD-adjusted estimates directly to one multivariate BLUP solve. The return value keeps raw RHE estimates, adjusted prediction parameters, component-specific values, and the combined value distinct.

Packed genomic terms created from the same gs_bed() resource are grouped by prepared-source identity. All marker-weight components share one decode for the phenotype moments, one decode per bounded common-probe batch, each batched multivariate covariance application, the diagonal preparation, and the final component-value application. Generic dense, pedigree, group, and sparse relationship terms retain the ordinary matrix-free operator route.

Controls expose probe count/batch/seed, packed threads, moment rank and conditioning limits, PSD floors, multiple-RHS strategy, PCG tolerances and iteration cap, block numerical-rank tolerance/restart cap, true-residual tolerance, and the explicit unavailable-diagonal fallback. Packed execution is the only OpenMP layer. Diagnostics attach a purpose to every traversal and report marker visits, decode passes, RHS width, component count, thread resolution, and decode/compute time.

The supported scope is complete traits on common individuals, one common full-rank fixed design, and arbitrary named PSD relationship components. Multivariate AI-REML, determinants, missing traits, marker effects, PEV, reliability, single-step GBLUP, and longitudinal models are not part of this workflow.

Complete-record multivariate AI-REML packed reuse

The existing gs_multivariate_model() plus gs_model_fit() stochastic AI-REML route retains its finite log-Cholesky parameters, derivative matrices, optimizer, likelihood-monotone acceptance, and determinant calculation. The execution-only marker_coefficient_reuse control now also applies to compatible packed terms in complete-record multivariate fits.

For \(K_c = W_c D_c W_c^{\mathsf T}/s_c\), gsuite obtains \(W_c^{\mathsf T}Z\), \(W_c^{\mathsf T}PZ\), and \(W_c^{\mathsf T}PY\) with the existing packed multi-RHS transpose. Components sharing a prepared BED source share the unweighted transposes; their weights and normalizations are applied only in marker-space reductions. The resulting small trait matrices are

\[ \begin{aligned} (PZ)^{\mathsf T}K_c Z&=(W_c^{\mathsf T}PZ)^{\mathsf T}D_c(W_c^{\mathsf T}Z)/s_c,\\ (PY)^{\mathsf T}K_c PY&=(W_c^{\mathsf T}PY)^{\mathsf T}D_c(W_c^{\mathsf T}PY)/s_c. \end{aligned} \]

greml contracts these matrices with the unchanged log-Cholesky derivative \(H_j\). gsuite reconstructs all \(K_c PY\) for a shared source with one fused packed output call; greml then forms \(K_c PY H_j\) and retains its existing batched information projection and AI definition. Native multivariate marker coefficients are not required or exposed.

The double-precision \(W^{\mathsf T}Z\) cache is bounded by marker_coefficient_cache_memory_limit, keyed by prepared-source and explicit probe identity, supports compatible prefix extension, and isolates independent probe streams. Automatic mode falls back term-by-term to direct covariance application; required mode fails for an unresolved packed term; off is the original direct route. Dense/operator-defined relationship terms remain direct.

PLINK mode validation remains in gbits and rejects non-SNP-major BED files. Packed transposes process all trait/probe columns together with disjoint marker rows, and shared-source reconstruction processes every component/trait column in one output call. Diagnostics report cache lifecycle, \(W^{\mathsf T}Q\) source, packed purposes, contiguous SNP runs, decoded records, fused width, components, requested/resolved threads, avoided/direct applications, and reconciled totals. Thread count does not multiply decode passes.

Experimental multivariate determinant-deferred acceptance composes with this reuse path without altering its sufficient statistics. It evaluates rejected trials with score only in the current accepted AI metric, completes AI once for an accepted state, and runs the unchanged likelihood exactly once after numerical convergence. Optimizer convergence and fixed-probe Monte Carlo adequacy are separate diagnostics; an iteration-limit exit has no finalized likelihood. Exact singular covariance boundaries, missing or unbalanced traits, unbalanced/indexed longitudinal structures, PEV/reliability, and marker-effect output remain outside this route.

Formula and variance-component pedigree REML

The persistent covariance constructor gs_model_pedigree(sire, dam, observation_animal, name, ...) supplies a gs_covariance term to gs_model() and gs_model_fit(). It preserves the ordered-parent and optional Mendelian variance contract in the persistent interface. The formula kernel below uses the distinct gs_pedigree() constructor for pedigree data frames.

gs_pedigree(), gvc(), gs_fit_control(), and gfit() provide the bounded in-memory formula interface documented in formula-reml.md. One single-trait formula or named multitrait formula list is lowered to greml’s installed public compiled-MME interface when compiled controls are explicitly requested. Those single-trait models accept one or two (1 | variable) terms: at most one pedigree and one independent identity component, plus a homogeneous residual. Multitrait models retain one shared pedigree intercept per trait. This route requires complete positive-interior starts. General formula calls accept either argument spelling and have optional starts and a broader independent-component scope, as described in the pedigree tutorial.

Fixed effects use base R model machinery. Random incidence is compiled from record IDs, trait/member identities, and unit weights; gsuite does not construct an R random-effect incidence matrix. Missing trait responses are absent scalar records in complete or incomplete patterned residual units. The immutable plan performs one record pass, a direct session retains one symbolic structure, and successful optimization reuses the accepted likelihood and response for one BLUP finalization.

The control maps direct and support thread domains, exact trace selection, and all public resource limits, including required/configured peak dense RHS workspace. Rejected preflight and unsupported support execution retain greml’s status and perform zero solver work. PARDISO coordinate-exact traces and CHOLMOD selected inversion are reported without capability substitution.

Balanced structured longitudinal workflow

gs_balanced_longitudinal_model() is one shared specification for the complete pipeline. It accepts a finite n by T phenotype matrix, common full-rank sample design X, full-rank time design F, and named components constructed by gs_longitudinal_component(). Component c supplies a relationship base, an explicit T by r_c basis B_c, and an interior r_c by r_c coefficient covariance Q_c. The model is

\(\operatorname{E}(Y) = X C F^{\mathsf T}\),

\[ V(A) = \sum_c K_c A B_c Q_c B_c^{\mathsf T} + \operatorname{residual}(A). \]

The residual is homogeneous independent by default or complete-record unstructured SPD. No dense production Kronecker design, GRM, or nT covariance is constructed.

gs_longitudinal_fit() runs RHE once, maps its separately preserved adjusted prediction covariances to explicitly reported finite log-Cholesky REML initial coordinates, runs AI-REML once, and optionally runs BLUP once. RHE returns raw and adjusted Q_c, observed-grid covariances, residual adjustment records, moment rank/condition/Monte Carlo diagnostics, fixed-design diagnostics, and packed traversals. REML retains the existing derivatives \(B_c \frac{\partial Q_c}{\partial\eta_j} B_c^{\mathsf T}\) and uses packed marker sufficient statistics without one source traversal per Q_c parameter.

The new workflow defaults to experimental determinant-deferred acceptance. Score stationarity is the configured raw-score rule or the joint curvature-scaled-score/Newton-decrement rule at a valid accepted state; it does not require another accepted micro-step. Line-search exhaustion is never a convergence condition. Numerical convergence and fixed-probe stochastic precision are separate. Successful convergence has zero iterative likelihood, determinant, and SLQ work and one final likelihood-only evaluation; final score and AI are retained.

BLUP returns latent \(U_c = K_c \alpha B_c Q_c\), component trajectories \(U_c B_c^{\mathsf T}\), their total, and fixed/fitted trajectories. Explicit compatible B_new and F_new matrices provide common-new-grid prediction. The balanced route remains the optimized common-grid convenience path.

Missing and irregular observed-record longitudinal workflow

gs_longitudinal_records_model() represents only observed scalar records. Its y, subject vector, optional reporting time, and rows of the general fixed design have common length Nobs. Each component supplies an Nobs by r_c basis and components may have different ranks. gs_record_longitudinal_residual() maps every record to a homogeneous or grouped diagonal variance and requires every group to occur. Absent records are never imputed or silently dropped.

The native operator scatters record RHS values through each component basis to an n by r_c subject matrix, applies the relationship to all coefficient/RHS columns together, multiplies by Q_c, and gathers to records. RHE, REML, likelihood and BLUP use this same operator. Raw RHE estimates, adjusted prediction estimates, and finite-interior REML initial values remain distinct. Record REML defaults to determinant-deferred acceptance with fixed common probes, zero iterative likelihood/SLQ calls, and one final likelihood after verified convergence.

Packed components sharing a prepared BED source share SNP-major \(W^{\mathsf T}Z\), \(W^{\mathsf T}PZ\), \(W^{\mathsf T}PY\), and fused reconstruction. Diagnostics report cache activity, purpose-specific traversals, fused RHS width, component count, marker visits, decode passes, and requested/resolved threads. gs_longitudinal_fit() returns latent subject coefficients, observed values and residuals, and optional values for explicit prediction-record subject/fixed/basis rows.

Deliberate exclusions are pedigree and sparse MME, multiple response channels per record, AR(1)/continuous-time residual correlation, exact singular boundaries, PEV/reliability, marker effects, and single-step models.

Sparse-direct pedigree methods

The three pedigree APIs use the same canonical lower-CSC MME for their supported direct backends. Retained factors supply exact determinants and bounded coordinate traces or CHOLMOD selected-inverse traces. Rejected likelihood-only trials omit trace/AI work; accepted response coefficients finalize BLUP without another response solve. Backend factor-memory limits are preflight bounds, not strict allocator quotas.

Use the model catalog for backend/inference status, multitrait workflow and longitudinal workflow for record contracts, and prediction error for complete-MME PEV, CHOLMOD selected inversion, and bounded PARDISO retained-factor target solves. PARDISO selected inversion remains unsupported. Direct and support thread domains and requested/resolved methods remain distinct diagnostics.

Independent LD references

Glist is the single public genotype and LD descriptor. Augment it directly with gmat’s signed, disk-backed sparse LD:

Glist <- gbase::gprep(Glist = Glist, task = "sparseld",
                      out_prefix = "reference_ld", r2 = 0.01,
                      max_distance_bp = 1000000)
saveRDS(Glist, "Glist-EUR.rds")
Glist <- readRDS("Glist-EUR.rds")
gbase::verify_ld(Glist) # once after transfer or download
scores <- gbase::ld_scores(Glist)

Glist$sparseLD retains authoritative marker order, alleles, chromosome, position, reference sample count, file paths, checksums and construction settings. ldfiles, idsLD, and rsidsLD identify the same persisted resource, reference samples, and matrix-ordered markers. RDS saves the descriptor, not the external LD files. validate_Glist() checks this identity and alignment without rescanning the CSR; verify_ld() performs the full native structure, checksum, and recorded-score verification after a transfer or when corruption is suspected. Neither check establishes population suitability. Scores are computed from that same resource:

scores <- gbase::ld_scores(Glist, annotations = annotation,
                           corrected = TRUE)

Annotation rows are aligned by exact marker identity. Corrected scores use the selected reference sample count and retain negative corrected contributions. All scores reflect retained sparse pairs; omitted correlations are not reconstructed. Summary harmonization remains separate.

Regional matrices and reusable eigensystems

Select a region using unique marker IDs already in reference order:

R_region <- gbase::ld_region(Glist, markers = region_markers)
eig <- gbase::ld_eigen(Glist, markers = region_markers)

gmat opens the checked resource once per high-level call and reuses one indexed reader across regions; no full CSR is loaded and only the active dense matrix is materialized. gsolve computes the complete symmetric eigensystem using its existing Eigen backend and verifies reconstruction, residual and orthogonality. The result contains descending values, column vectors with marker row names, and numerical diagnostics. Singular or indefinite matrices are not silently repaired: zero and negative eigenvalues remain available to method-specific checks. Thresholded or windowed sparse LD is not necessarily positive semidefinite.

The eigensystem is computed from the checked regional matrix on demand and can be reused across traits or fits for that region. The full decomposition takes quadratic storage and cubic computation; increasing the dimension limit is an explicit choice. Omitted sparse correlations remain zero, the diagonal is one, and no harmonization, clipping or PSD repair occurs.

For a BED Glist with independent regions saved by gbase::gprep(task="eigen"), gcorr(Glist, stat, method="hess", nthreads=4) reuses its region definitions and reference sample selection and returns estimates plus compact moments in fit$projection. Named Z/N arrays must already be oriented to Glist’s a1 alleles. Stored eigenvalues alone cannot project statistics, so the required vectors are recomputed per region and discarded. Both compact-fit overlap inputs default to zero and can be overridden.

For explicit disk-limited joint HESS, gcorr_prepare_hess() streams BED regions, uses validated physical BED/FAM rows without rebuilding chromosome-wide metadata or name lookups, and releases R matching buffers before decoding. It decomposes the smaller sample- or marker-space Gram and projects all supplied traits together. It discards large matrices/eigenvectors after each region; gcorr_fit_hess() reuses the serializable compact moments in the existing joint correction and conditional uncertainty calculations. This needs no persistent LD files. Spectral controls and trait statistics are fixed at preparation, and overlap inputs can be supplied at fit time. Automatic compact preparation limits the conservative total rank to half the smallest supplied N. Explicit fit-time genomic-control factors and provenance scale unchanged moments without changing N or repeating decomposition. Conditioning and region-specific inference blockers are reported. Compact fits default to uncertainty="hess", using the original software’s full-rank h2/covariance variance systems and genome rg region jackknife. Point estimates retain the documented row-rank equations. The software variance convention has no regional rg SEs or joint sampling covariance; uncertainty="moment" selects the previous signed conditional fourth-moment calculation, and "strict" retains its parameter checks. Prepared-LD interfaces keep their strict default. Invalid jackknife replicates remain missing; out-of-range rg receives only a descriptive SE, without Wald inference. Approximate uncertainty does not establish statistic/reference-scale validity or broad empirical calibration. See the R workflow.

The gcorr R interface also prepares all named regions from one Glist-backed native reader and retains one PreparedRegionalLd per region. HESS, rho-HESS, SUPERGNOVA and LAVA have single-region overloads that reuse a complete checked eigensystem; LDSC and GNOVA already accept prepared gmat score inputs. gsolve can verify imported eigenvalues and vectors against a regional matrix without rerunning the eigensolver. This dense verification is performed once on native import, followed by reuse across fits. The gcorr R interface connects selected methods to these resources, including explicit joint HESS/rho-HESS correction with regional and total uncertainty. It adds no cross-ancestry estimator and also exposes shared LAVA matrices and native conditional inference.

Across these methods, result tables use explicit reference, region, annotation-category and architecture-component labels. Named h2/cov/rg matrices share one layout; method-specific uncertainty, native diagnostics and inconsistent pairwise estimates remain visible. See result conventions.

The interface also exposes SumHer architecture-component and reference-total h2/cov/rg. Multiple traits use the native coordinated marginal and pairwise fits, with joint fixed-weight jackknife covariance when blocks are supplied. Tagging weights are explicit and separate from architecture-weighted scores; no cross-ancestry or full upstream-software parity is implied.

Regional LAVA also accepts explicit binary-trait metadata, reusing Glist allele frequencies for native logistic-Wald reconstruction. Observed and liability h2 remain distinct; binary/mixed pairs use the existing screening and optional bivariate simulation inference. LAVA shares each regional fit across traits and exposes partial correlations, multiple regression, explained variance and supported conditional/omnibus inference.

Native regional execution now supports multiple named references and regions. Each request explicitly assigns a reference to every trait; all traits in that request must agree. Independent requests may use different references. Results retain analysis, region and reference labels, along with method diagnostics or a per-request failure. Execution is sequential and writes no result files. Joint HESS and rho-HESS can reuse prepared regions within one selected reference while retaining their existing disjoint-region and joint uncertainty assumptions. Reference labels are caller-declared identities, not an ancestry validation step.

For continuous-trait regional LAVA correlations, gcorr_lava() accepts an explicit inference list. The native bivariate procedure supplies screened p-values, 95% simulation intervals and Monte Carlo diagnostics while retaining raw point estimates. See the regional inference contract. This adds no new numerical algorithm or statistical calibration claim.

Native colocalisation

gmap provides approximate-Bayes-factor colocalisation for one region and two continuous traits, assuming at most one causal variant per trait and non-overlapping samples. It returns H0-H4 probabilities and explicit prior-sensitivity summaries. The implementation has independent numerical checks and pinned coloc comparisons; these establish numerical agreement, not broad statistical calibration. The gmap::gmap_coloc() route accepts two aligned regional beta/SE tables, explicit phenotype SDs and priors without an LD object. It returns regional probabilities and conditional-H4 marker support, distinct from component-pair colocalisation. See the worked example. GWAS association is implemented separately in glma; Bayesian fine-mapping and posterior-set interpretation are described below.

Native Bayesian inference

The gbayes::gbayes() R workflow exposes single-trait and joint multivariate regression from a prepared individual-level BED Glist or from checked summaries and Glist sparse LD. The BED route uses matched pre-adjusted phenotype columns and native packed-genotype preparation; its joint route can learn residual covariance from observed trait cross-products. BayesC/BayesR support fixed parameters or optional effect-covariance and mixture-weight learning, exact chain seeds, chain-level threading and native posterior summaries. Multiple traits require explicit independent sampling or the common-N shared-LD covariance model. Variable sample sizes are an explicit approximation and, for multivariate fits, require independent sampling and study-size metadata. Shared-LD multivariate sampling covariance can be learned when the response cross-product and likelihood dimension are supplied. Scalar BayesC/BayesR expose fixed, exact in-sample and reference-LD adjusted residual policies. Exact learning uses an explicit response sum of squares, likelihood dimension, scaled inverse-chi-square prior and complete same-sample LD. The reference-LD route uses standardized summary normalization and records its base residual variance, genetic variance and adjusted working variance separately. Latent base-effect covariance is labelled separately from realised genetic covariance. The annotation argument exposes logistic inclusion, continuation-probit pattern/component probabilities and log-linear variance models, with fixed or learned coefficients and native diagnostics. The same interface now binds grouped latent covariance priors and optional native posterior genetic covariance, partition activity/covariance and selected-individual score uncertainty. Heritability normalization requires explicit phenotype variances. Group summaries retain cross-group LD; score uncertainty excludes residual noise. Named study-specific references additionally expose native independent-provider preparation with differing LD, sample sizes and marker subsets. Multiple studies may inform shared effects of one trait within a model of at least two traits. Study identity, trait identity and reference provenance remain distinct. This is not an overlap correction or an automatic ancestry-specific effect model.

gbayes provides single-trait BayesC and BayesR with fixed parameters or explicit learning of residual variance, slab variance and mixture weights. Prepared scores and cross-products can use dense storage, gmat LD operations with explicit marker scaling, or prepared factors with reduced residuals. Factorized execution avoids a dense cross-product matrix. Truncated inputs use explicitly projected scores; residual learning then requires retained-space normalization (dimension and response sum of squares). Results include posterior means, empirical posterior variances, marginal inclusion probabilities and compact per-chain summaries. BayesR also returns marker component probabilities; both methods summarize learned parameters. Residual-variance learning requires compatible phenotype sums of squares and an explicit Gaussian likelihood dimension, not just LD and scores. The gmap() R binding exposes BayesC/BayesR separately across named traits and regions, including optional variance and mixture-weight learning. MCMC learning requires explicit inverse-chi-square or Dirichlet hyperpriors. Residual learning also requires complete in-sample LD, phenotype sum of squares and likelihood dimension; these are never inferred from the reference panel.

The scalar reference-LD recipe has five-repetition qualification for matched continuous-trait BayesC and BayesR analyses, including heritability, effect and score recovery and four-chain diagnostics. Independent scalar calculations, dense/gmat agreement and matched sblr checks support its native implementation. See the qualification record for the exact LD construction, priors and chain lengths. Supplied cross-products must define a valid positive-semidefinite operator; structural checks do not repair or certify arbitrary thresholded LD. A separate five-repetition design qualifies disjoint-block residual learning, the BayesS MAF exponent, shared-LD latent effect and sampling-covariance learning, and strong-LD relocation under their stated contracts. Scalar and multivariate saved states also reproduce uninterrupted chains exactly. Broader annotation models and arbitrary overlap structures remain outside that evidence.

Native regional Bayesian analysis

The Bayesian-enabled gmap::gmap library fits or summarizes one-region BayesC/BayesR models through installed gbayes. It preserves marker alignment, region coordinates, trait/effect units and the complete native posterior. Results include marker rankings and the posterior expected active-marker count, calculated as the sum of marginal PIPs. Existing gbayes priors and learning controls pass through.

These posterior summaries describe marker support. Separate set helpers below interpret marginal PIPs or component assignments. These summaries do not themselves perform colocalisation or estimate regional heritability. The separate gmap_coloc_components() interface consumes multi-effect evidence. Inputs are already aligned. The gmap() R binding returns regional posterior estimates, native set controls and chain diagnostics for named regions and traits, with regional preparation shared across traits. One regional eigensystem supports an ordered retained-positive-mass schedule; the first threshold accepted for every trait is selected and every attempt remains visible. The colocalisation target acquires no Bayesian dependency. Focused checks cover posterior reuse, metadata/ranking semantics and installed interfaces.

Native single-effect inference

gbayes now provides a deterministic single-effect normal model for one prepared marker set. It reuses the normal conditional calculation and accepts the existing dense, LD-backed and factorized inputs. Outputs include marker assignment probabilities, conditional and averaged effect moments, and model evidence against no effect. Residual and effect-prior variances are fixed; marker prior weights can be specified explicitly.

The model assumes exactly one active marker, so assignments sum to one even when evidence is weak. They are not unconditional association probabilities. This is a reusable building block; the separate multi-effect entry point below adds iteration and convergence handling. gmap now consumes assignments for credible sets; multi-signal colocalisation remains separate. The standalone gbayes::gbayes_single_effect() R interface accepts an ordinary data frame with unique rsids and either b/seb or score/information, plus fixed positive residual/effect variances and optional marker prior weights. It returns per-marker assignment probabilities and moments. It does not expose every native input or parameter-learning option. Focused independent numerical integration and installed-interface checks cover the native calculation.

Native multi-effect inference

gbayes also provides single-trait sum-of-single-effects fitting through sequential IBSS updates. It reuses the existing single-effect calculation and dense, LD or factorized operators. Outputs include component assignments, combined effect means, variational variances/PIPs, and an objective trace with explicit convergence or sweep-limit status.

The objective is an evidence lower bound relative to the null likelihood, not an exact model Bayes factor. The approximation factorizes across components; convergence does not guarantee a global optimum or calibrated uncertainty. The maximum component count is supplied, and components may select the same marker. The default fixed-variance route leaves every component active. Optional EM learning estimates a variance for each component and compares it with zero; inactive components contribute no effects or PIPs. Optional residual-variance learning uses expected residual sums of squares with compatible normalization. Inner learning limits and outer convergence are reported explicitly. gmap supplies credible sets only for active components while retaining inactive indices and all component variances. The gmap() R binding exposes this route separately across regions and traits, with optional residual- and effect-variance learning. Residual learning requires complete in-sample LD and explicit phenotype normalization. Small analytic/evidence and operator checks cover the native calculation; full susieR equivalence is not claimed.

Native posterior sets

The optional gmap Bayesian component supports two explicitly labelled routes. Marginal-PIP sets reuse sblr’s LD-neighborhood and approximate multi-signal selection, including incomplete-set handling and alternative lead/removal rules. Their mass is an expected active-marker count, not a probability that at least one causal variant lies in the set. This route consumes BayesC/BayesR PIPs.

Component-based credible sets consume single-/multi-effect assignment rows. They retain each component, report achieved probability coverage and flag exact duplicate memberships. Shared diagnostics include lead-based and full pairwise LD summaries. Low purity is flagged without trimming away coverage. Sparse LD omissions are explicitly recorded; purity is evaluated on the supplied matrix.

These helpers validate probability domains and exact fit/LD marker order, preserve available convergence/likelihood status, and do not refit or align data. Learned zero-variance components are recorded without sets; the raw assignment route still treats every row as a component. Weak fixed-variance components can produce diffuse sets; a complete credible set does not establish a signal. The R binding exposes these controls through control$sets, with explicit set meanings and component eligibility diagnostics.

gmap_sets() converts the nested native result to set-level and marker-level tables without changing those meanings. gmap_bfdr() operates only on fitted marker PIPs: at a supplied threshold it reports the selected count, the expected false-discovery count sum(1 - PIP) and their ratio. A requested BFDR target selects the largest PIP-threshold group below that bound. It does not require LD, retained posterior draws or another model fit, and it is not a frequentist FDR guarantee.

Prior-variance learning is empirical Bayes, not posterior uncertainty over model size or a test proving signal existence. The EM search is local; inactive components are reconsidered from a positive starting variance each sweep. The reusable learned single-effect entry is available independently. Existing fixed-variance entry points retain their previous numerical behavior.

Regional Bayesian R workflow

For an individual-level BayesC/BayesR fit, gmap::gmap_bed(fit, Glist, regions) reuses the saved global posterior. The Glist must carry in-sample sparse LD prepared from exactly the fitted people and markers. Sample order, physical BED offsets, alleles and frequencies are checked. Within-region LD must be complete by default. Native gmap groups the fitted marginal PIPs; the result is neither a regional refit nor a component credible set. This route also accepts a joint-trait BED fit and retains trait labels. See the BED-to-gmap workflow.

gmap::gmap() accepts an aligned Glist and exact sufficient-statistic columns marker, allele1, allele2, score (\(X^T y\)) and marker_ss (\(\operatorname{diag}(X^T X)\)). It fits each named region/trait separately using BayesC, BayesR or multi-effect inference. It does not convert GWAS z scores or beta/SE automatically. The supplied LD and marker scales define the working cross-products; a reference approximation does not become exact cohort LD.

fit <- gmap::gmap(list(height = stat_height, weight = stat_weight), Glist,
            regions = list(locus1 = markers_locus1, locus2 = markers_locus2),
            method = "multi_effect",
            prior = list(residual_variance = 1, effect_variance = .02),
            control = list(effects = 4, estimate_effect_variance = TRUE,
              eigen = list(retained_positive_mass = 1)))
summary(fit)
plot(fit, region = "locus1", trait = "height", quantity = "pip")

Prior values above are illustrative and must match the supplied data scale. Controls and seeds are shared across fits; use separate calls for different priors or independent random streams. Each verified region is decomposed once across traits. The default retained-mass schedule is c(.995, .99, .95). diagnostic_summary records the selected requested/achieved mass, rank and projection change, while eigen_attempts records each attempted threshold. Set control$eigen$retained_positive_mass=1 for the exact full-rank route. Residual-variance learning is restricted to that full-rank route. Component colocalisation also requires a converged full-rank multi-effect fit; use the same explicit setting when the fit will be passed to gmap_coloc_components(). control$sets exposes method-specific native set controls. BayesC/BayesR retain pooled and per-chain PIPs, with native between-chain spread summaries; these are not convergence certificates or Monte Carlo standard errors. Scalar summary BayesC/BayesR in gmap() and gbayes() also expose optional control$ld_swap relocation moves on their dense/sparse and eigen likelihoods. The native gbayes sampler owns neighbour selection and Metropolis-Hastings acceptance. Probability defaults to zero; fit$ld_swap reports proposal counts, not convergence. The existing timing studies used swaps disabled.

For traits analysed in independent samples, the saved fits can be compared:

shared <- gmap::gmap_coloc_components(
  fit, fit, Glist, trait1 = "height", trait2 = "weight", regions = "locus1",
  control = list(complete_region_coverage = TRUE,
                 non_overlapping_samples = TRUE))
summary(shared)
plot(shared)

The two assertions require complete regional coverage and independent trait samples; they must describe the actual analysis. This route requires converged full-likelihood multi-effect fits, uniform marker priors and the same complete LD reference. It consumes stored component evidence without fitting or eigen computation. LD files must remain available for checked purity calculations. H0–H4 probabilities belong to individual component pairs, not the entire region. Inactive, low-purity, incomplete and duplicate components remain labelled; no eligible pairs is not evidence for H0. Summaries and plots of completed results work from saved R objects without accessing LD files.

OLS beta/SE/n preparation is available through gmap_stat(); see below and the worked simulation tutorial. No cross-reference or overlapping-sample colocalisation is exposed.

Linear-regression GWAS input for mapping

gmap::gmap_stat() performs a separate OLS summary conversion in the R layer: for sample phenotype variance \(s_y^2\), it reconstructs \(d_j=(n-1)s_y^2/[(n-2)SE_j^2+\hat\beta_j^2]\) and \(X_j^\mathsf{T}y=d_j\hat\beta_j\). The common-sample, intercept-only, continuous-outcome contract is explicit. LD supplies the off-diagonal cross-products; external LD is approximate. Numerical fitting and sets remain native. See the worked example.

MR through the R interface

gcorr exposes GSMR, GSMR2, IVW, Egger and weighted-median components through gcorr_prepare_mr() and method-specific fitters. Preparation consumes two aligned gbase effect/SE summaries and an explicit signed-LD matrix; it does not accept an LDlist, align alleles or infer ancestry/sample overlap. The serializable result records every selection and LD decision. Focused fitters and gcorr(prepared, method=...) return HEIDI, estimate-unit and method-specific diagnostics. See the MR workflow.

Covariance models through the R interface

gcorr_covariance() retains labelled continuous genetic covariance and full joint sampling covariance. Focused SEM and GGM functions expose recursive SEM, descriptive partial-correlation networks and fixed-graph GGM. Parameter tables, trait units and estimation metadata remain explicit. Fitted SEM paths and free GGM edges support native joint Wald tests. Coordinated multivariate GNOVA annotation categories and SumHer architecture components can reuse the same models, with cross-group contrasts based on the retained joint jackknife covariance. See the covariance-model workflow.

Native matrix-free execution experiments

The generic greml PCG route now supports same-state adaptive trace-prefix extension. Native-only opt-in controls add bounded raw inverse-probe warm starts, coarse component-action reuse, solution-subspace recycling and independent-basis residual trace estimation. The latter two require experimental evaluation; determinant precision and the packed coefficient path remain separate. See native contracts. No new R controls or general scaling qualification are implied.

Standalone gbayes::gbayes(data=..., Glist=..., method="bayesc") and BayesR accept automatic phenotype-scaled BED priors with parameter learning and retained diagnostics. Scalar BED responses are centred on their original units; see the defaults and controls for initialization, explicit-prior compatibility and convergence limits.

Fixed-covariance compiled prediction also omits unused global REML first/second covariance-derivative matrices and RHS arrays. It retains the same numerical MME, residual precision caches and preconditioner; full estimation keeps its existing derivative assembly. This removes quadratic-in-parameter-count storage from prediction without imposing a memory ceiling. The refreshed native capability is prediction_value_only_assembly.

Extended standalone estimation and model inputs

The YAML adapters add deterministic sparse-direct REML for unstructured/FA genetic covariance and scalar or separable multivariate AR1/exponential residual parameters, with natural starts, whole-block or named-coordinate fixing, positive-definite parameterization, likelihood acceptance and separate convergence reporting. Small balanced univariate and bivariate references recover analytic REML estimates on PARDISO and CHOLMOD; fixed blocks remain exactly unchanged. Packed EM-REML has its narrower scope below; industrial REML precision and runtime remain unqualified. Serial REML uses exact observed joint blocks, including missing traits/visits; it does not freeze a whitening transform while changing covariance parameters.

Weighted multiple-individual incidence works on compiled and packed prediction routes, with explicit literal or data-column weights and independently checked mixture predictions. Sparse external covariance files feed retained native sparse inverse operators for additional compiled PCG effects. CHOLMOD uses selected inversion for precision diagonals; PARDISO currently uses bounded unit-solve batches, which can dominate setup for large relationships. Input sparsity does not bound factor fill. Supplied joint trait-by-visit residual blocks support time-varying trait covariance and general finite nonseparable serial covariance, including missing measurements, on compiled models.

Fixed/combined group models and absorbed phantom group-parent slots match independently assembled equations. Packed learned Gamma is checked against an analytic absolute-base covariance and a refit supplying the prepared estimate. The added tests are small deterministic references; no new industrial performance claim follows. The complete examples live in the model example directory. Existing R interfaces are unchanged. Packed fixed/combined group models prepare ordinary inheritance columns once. PEV precision work has resumed: an optional independent pilot selects the final residual-probe count within the explicit solve cap. Its samples are discarded; the final estimator and MCSE use separate random streams. More work can reduce Monte Carlo error but does not guarantee meeting the target, improve a fixed-budget estimator, or establish coverage.

The standalone packed ssSNPBLUP REML route estimates unstructured genetic and independent residual covariance through native EM posterior moments. It retains the pedigree factor and packed genotype provider, supports missing traits and fixed additional identity-effect covariances, and reports exact or stochastic scores with explicit solve/MCSE controls. Optional start-to-final restricted-likelihood changes use exact references or paired SLQ, with separate sampling MCSE and Lanczos-depth sensitivity. Omission performs no global determinant work. There is no likelihood-based EM line search. Small independent references establish numerical correctness; industrial runtime, stochastic precision calibration and model extensions require further qualification.

The standalone application is organized by model route and task; see the greml/apps/model_driver/ organization. Source organization does not change the model-specific estimation or uncertainty qualifications. The standalone gsuite run model.yml command dispatches the task declared in that same YAML for every supported model route. It manages temporary preparation and native work, collects results and failure diagnostics, and records elapsed times, native peak RSS and separately qualified sampled process-tree RSS. gsuite validate model.yml and gsuite describe model.yml check or inspect inputs without fitting. Detailed scratch and reporting instructions accompany the source bundle. An optional YAML evaluation block compares a named individual-level prediction output with supplied long-format truth after an accepted fit. It reports Pearson/Spearman correlations, RMSE/MAE, mean prediction-minus-truth bias and truth-on-prediction calibration, optionally stratified by mapped genotype IDs. Held-out individuals have no observed training responses in any trait; the block does not create the holdout split. No truth-based centering or tuning is performed. Label/join acceptance is separate from numerical residual gates and does not establish accuracy or replicated precision. The updated frontend collects these outputs through the same gsuite run model.yml command. An optional private GenomeDK binary producer/installer preserves the same model frontend and native equations, with a minimal Conda runtime and fresh-prefix exact/backend checks before creating an archive. Packaging contracts tested on Windows are separate from actual Linux numerical qualification. This adds no R controls, covariance family, PEV precision claim or industrial throughput claim. See the standalone user guide for installation and a complete model, and the GenomeDK validation record for the accepted artifact, source revisions and scope of the operator-reported checks.