Numerical and mixed-model routes
Numerical and mixed-model routes
This guide explains route choice and scaling. The model catalog owns capability status, the primitive registry owns reusable mechanism mappings, and architecture owns package dependencies. A native primitive does not establish a qualified public R route.
Estimate or solve from R
These examples use the existing root development integration. They are not provided by the seven public source packages. Future greml/gsolve R packaging is deferred, and its interface and separate or combined distribution are not yet decided.
greml() estimates variance components; gsolve() holds them fixed and returns fixed effects and random-effect predictions. Both use the same formula/data/pedigree and explicit gvc() specifications as gfit().
fit <- greml(formula = BW ~ sex + (1 | id),
data = data, pedigree = pedigree)
summary(fit)
fit$variance_components
# The same model/data specification, now with fixed covariance values.
solved <- gsolve(formula = BW ~ sex + (1 | id), data = data,
pedigree = pedigree, vc = c(id = 1, residual = 1))
# Optional shortcut: reuse the prepared model and its fitted covariance values.
resolved <- gsolve(fit)
resolved_pcg <- gsolve(fit, control = list(linear_solver = "pcg"))Starting values are optional for greml(). Direct gsolve() calls require all covariance values. For gsolve(fit), the default is to use the fitted values, requiring successful status and REML convergence. A complete named vc list/vector can override them; names identify components, and matrix coordinates retain the fitted trait/basis order. Such an override solves the stated fixed model and does not turn a failed REML fit into a successful one.
The handoff reuses prepared records, sparse design and ID mappings. It supports saving/reloading the fit with saveRDS() / readRDS(). Numerical factors are rebuilt for the new solve. Fits therefore retain their analysis data; consider this when sharing saved objects. Supplying a new control replaces the route-specific controls; omitting it inherits the fit’s controls.
Single-trait, multitrait and longitudinal pedigree formulas support the existing PCG and sparse-direct routes. CHOLMOD is the default direct backend. The explicit single-response component route, including correlated channels and grouped/AR(1)/exponential residuals, uses its existing direct solvers. Backend availability depends on the installed build. General gs_model workflows and low-level numerical gsuite_solve() retain their own APIs.
gfit(task = "reml") and gfit(task = "solve") remain available. The compiled route selected by gs_fit_control() is REML-only: greml() can use it, but its fit cannot be handed to gsolve(). This interface change adds no new numerical algorithms or qualification claims. The pedigree tutorial shows the three main data-frame workflows. REML fits already include predictions, so a second solve is optional.
Numerical and statistical boundary
gsolve owns checked column-major matrix views, leading dimensions, caller-owned workspaces, square and rectangular operators, dense/sparse SPD factors, PCG, preconditioners, deterministic probes, SLQ, and numerical diagnostics. Its operator contracts declare symmetry, definiteness, concurrency, and available diagonals. Rectangular operators provide forward/transpose application and bounded column materialization.
greml owns covariance models, fixed-effect projection, likelihood, score, standard average information (AI), optimization, prediction, and uncertainty. gsuite validates explicit R inputs and translates installed public results. Packed genotype operations belong to gbits; resource identity and LD belong to gmat. The optional gbits adapter supplies numerical operators without a reverse dependency from gsolve or greml to genotype libraries.
Covariance and coefficient routes
For single-trait covariance components, \(V = \sum_c \theta_c K_c\) and
\[ P=V^{-1}-V^{-1}F(F^{\mathsf T}V^{-1}F)^{-1}F^{\mathsf T}V^{-1}. \]
| Route | Computation and storage | Practical boundary |
|---|---|---|
| Exact dense | Factors observation covariance V; quadratic storage | Small reference models; prepared genomic BED terms are rejected |
| Matrix-free covariance | PCG applies inverses/projections; common Hutchinson probes estimate traces and SLQ estimates log determinants | Operator and bounded \(O(nk)\) Krylov/probe storage; solver and stochastic precision must both be adequate |
| Sparse MME | Assembles coefficient-space statistics and uses a retained sparse factor or a qualified coefficient operator | Fill and resource limits matter; backend, trace, inference, and uncertainty support depend on the specific API |
| Hybrid covariance | Uses structured sparse terms as a state-matched preconditioner for separately declared matrix-free terms | Components enter once; a pedigree plus genomic sum is not automatic single-step GBLUP |
| Iterative large-p projection | Nested PCG avoids retaining \(V^{-1}F\) and the full fixed-normal matrix; fixed determinant uses SLQ | Selected contrasts may be available while full fixed-effect covariance is unavailable |
The low-level gs_model() sparse/hybrid interface is distinct from the specialized pedigree and component-MME interfaces. Its limitations do not apply to all sparse-direct APIs. CHOLMOD selected inversion and PARDISO selected_solve target uncertainty are documented in prediction error. PARDISO selected inversion is unsupported.
Unbalanced multivariate covariance uses the observed principal submatrix of a trait-by-unit covariance: scatter, base-operator application, trait mixing, and gather. Random regression uses supplied basis columns. Neither operation requires a dense production incidence or Kronecker matrix. Selected prediction uncertainty computes only requested target blocks.
Scaling and execution
Prepared packed BED storage trades preparation work and memory for repeated applications; genomic production routes form no dense GRM. Sparse MME cost depends on coefficient dimension and fill. Matrix-free memory is dominated by operators and bounded Krylov/probe buffers. Nested fixed-effect solves add their own tolerance and determinant-precision requirements.
Independent multiple-RHS PCG solves columns independently. Explicit true block PCG shares a block Krylov recurrence and reports rank compression, deflation, restarts, and each original RHS’s true residual. Its orchestration owns one outer thread. A two-level preconditioner uses a caller-supplied, copied, rank-compressed coarse basis tied to one covariance state; it does not discover a basis automatically. Outer RHS/probe parallelism and operator-local OpenMP must not be nested. Requested and resolved execution are reported separately.
Route selection is API-specific. Low-level model routes are explicit; pedigree direct controls may resolve automatic to an available backend according to their documented policy. An explicit unsupported request does not silently change the statistical route.
Interfaces and evidence
The persistent R interface specifies handles, lifetime, ordering, controls, and low-level examples. The formula interface is a separate bounded lowering API. The package qualification limits remain; that page also records the focused constructor naming repair.
Model qualification, qgdata evidence, and the DMU workflow retain scoped numerical evidence and negative outcomes. They do not establish general biobank readiness or universal performance. Build and packaging commands belong in development and binary distribution.