Persistent explicit R interface

gsuite provides an explicit, experimental R interface to the gsolve, greml, and gbits native libraries. It is deliberately below a formula or data-frame layer: callers supply fixed-effect matrices, covariance components, response order, BED dimensions, and numerical controls directly.

Workflow

bed <- gs_bed("cohort.bed", n_samples = 100, n_markers = 1000)
model <- gs_model(
  fixed = cbind(intercept = 1, age = age),
  covariance = list(gs_genomic(bed), gs_residual()),
  method = "reml",
  route = "matrix_free"
)
fit <- gs_model_fit(
  model,
  response = y,
  start = c(genetic = 0.5, residual = 0.5),
  control = gs_control(threads = 4, multiple_rhs = "block")
)
summary(fit)

gs_bed() uses explicit physical dimensions and does not parse BIM or FAM files. Optional sample and marker selections are one-based, unique, ordered R indices. They are checked once and converted to the zero-based native contract. The selected order is preserved. Genotype dosage orientation and missing-call decoding are inherited from gbits.

The BED object owns immutable prepared packed storage, its selection, and its transformation. It supports raw, centered, scaled, and centered-scaled transforms and the public gbits missing policies. Centered-scaled mean imputation is the default. Marker weights are nonnegative, and genomic normalization is explicitly positive; by default it is the sum of weights. The matrix-free route uses gbits::gsolve_adapter and never forms a dense GRM.

Ownership and reuse

gs_bed, gs_model, and gs_fit are ordinary R lists containing registered, finalized native handles. Copying an R object shares its immutable native handle. Compiled models retain their BED operators, and fits retain their compiled model and response, so correctness does not depend on R garbage- collection order. Each fit or post-fit inference request creates an independent mutable greml::ModelWorkspace; stale workspace state is not reused.

Handles are process-local and are not serialized for reuse across R sessions. Save the explicit BED path, dimensions, selections, and model specification and reconstruct handles in a new process. Users never manage raw addresses.

Structured terms

The ordered-parent gs_model_pedigree() constructor below serves the persistent model workflow. The data-frame gs_pedigree() constructor serves the separate formula workflow; neither dispatches to the other. See qualification limits for the focused repair verification and remaining package limitations.

gs_dense_covariance() supplies a finite symmetric dense component. gs_genomic() supplies a prepared matrix-free genomic component. gs_residual() supplies identity residual covariance or a known positive diagonal. Residual weights are relative covariance diagonal values multiplied by the fitted residual variance; they are not precisions.

gs_group() takes one positive one-based group index per observation. It keeps declared level order, permits unobserved trailing levels through n_levels, and treats optional positive level_weights as relative level variances. The native implementation uses an independent-group operator or sparse incidence term and does not construct an R incidence matrix.

gs_model_pedigree() takes ordered one-based parent references. NA_integer_ and zero mean unknown; every known parent must precede its offspring, and gsuite does not reorder pedigrees. observation_animal supports animals without records and repeated records. Caller-supplied Mendelian variances must be positive. When absent, greml’s public pedigree helper derives them using zero supplied inbreeding coefficients. Matrix-free paths use pedigree recursion and sparse paths use Henderson contributions; neither dense A nor dense $A^{-1}$ is formed.

gs_sparse_random() copies the complete CSC slots of a Matrix::dgCMatrix. It represents $Z Q^{-1} Z^{\mathsf T}$, preserving coefficient-column order. precision is the positive diagonal of Q; alternatively, level_weights are positive relative coefficient variances and are inverted once to form diagonal Q. Identity and positive diagonal precision are supported by this constructor.

Models and routes

gs_model() inserts no intercept, drops no columns, and performs no rank repair. fixed = NULL means zero fixed-effect columns. Exact dense and matrix-free stochastic, sparse-MME, and hybrid routes are selected explicitly. Dense components are rejected by matrix-free and sparse routes. Sparse MME is currently REML-only and accepts structured terms plus exactly one residual. Hybrid models accept structured terms, one residual, and prepared genomic matrix-free terms. The native hybrid parameter order is structured terms, residual, then matrix-free terms; declared_components records caller order when this canonical order differs. Named starts avoid ambiguity.

Sparse-MME fixed-parameter evaluation returns exact coefficient-space solves, REML likelihood, and determinant ingredients. Its optimizer uses greml’s stochastic score machinery while retaining the exact sparse likelihood. A fixed-parameter hybrid call returns projection solves and predictions but no likelihood through the current facade; hybrid likelihood and optimization diagnostics are available through the stochastic optimizer path. Hybrid optimization is REML-only. No route is selected automatically.

Group, pedigree, and general sparse terms can also enter an ordinary matrix-free covariance sum. Hybrid construction pairs the structured sparse system with remaining matrix-free operators and prevents concurrent outer use of its shared sparse factorization. Use one outer thread and operator-local threads for prepared BED work.

A compiled model can be reused with compatible response vectors and starting parameters. Sparse symbolic analysis is performed once per independent fit workspace and reused across its numerical refactorizations; separate fits use separate mutable workspaces by design.

Controls and diagnostics

gs_control() exposes requested threads and seed; covariance and fixed-normal PCG tolerances; independent or true block multiple-RHS solves; automatic, identity, Jacobi, sparse-structured, or caller-basis two-level preconditioning; probe count, batching and Lanczos depth; adaptive precision caps; reference or safeguarded optimization; and bounded history, prediction, and uncertainty requests. Defaults are compact development settings, not universal production recommendations.

Invalid input and structural setup failures are R errors. Legitimate optimizer nonconvergence, iteration exhaustion, stochastic-precision failure, and linear-solve precision failure remain structured fit results (with an optional warning). diagnostics() keeps optimizer, Monte Carlo precision, linear-solve, preconditioner, execution, and linked-library information distinct. Matrix-free diagnostics include the preparation, response, and information PCG batches, per-column iteration and true-residual records, operator and preconditioner applications, processed columns, and block-PCG rank/activity counters rather than only aggregate wall time.

gs_ranef() returns named coefficient blocks in structured-term order. Sparse MME exposes group, pedigree, and general sparse-Z coefficient predictions. The current hybrid facade exposes coefficient predictions only for general sparse-Z terms; observation-space component predictions remain available for every hybrid component. fitted() excludes the declared residual component, and residuals() returns response minus those fitted values.

vcov(fit, type = "covariance") returns covariance-parameter uncertainty when available. type = "fixed" returns the full fixed-effect covariance only when formed; iterative large-p projection reports explicit unavailability instead of returning an empty matrix. gs_contrast() supports selected fixed-effect inference, including iterative contrast solves. gs_wald() supports fixed or covariance-parameter hypotheses and deliberately does not add p-values.

Multivariate and longitudinal construction

gs_multivariate_index() preserves caller row order and accepts each unit-trait pair at most once. Missing phenotypes are absent pairs. It records trait-local row positions, counts, declared dimensions, and balanced-complete status. gs_longitudinal_index() adds explicit one-based occasion indices and finite times; unit-trait pairs may repeat at distinct occasions, but a unit-trait-occasion triple is unique. Times are neither sorted nor transformed.

gs_trait_covariance() supports fixed full covariance, fixed positive diagonal covariance, a caller-supplied lower Cholesky factor, and optimized packed log-Cholesky coordinates. Packed order is column-major lower order: (1,1), (2,1), ..., (q,1), (2,2), .... Diagonal Cholesky coordinates are exponentiated, so optimized covariances are strictly SPD and are never silently projected to a PSD boundary.

gs_multivariate_model() combines unit-level dense, group, pedigree, sparse, or prepared-genomic operators with small trait covariances through greml’s indexed observed-principal-submatrix operator. It constructs neither a dense Kronecker matrix nor a global dense residual. Trait-specific fixed matrices follow index$trait_positions; global coefficient order is trait order and then caller column order. Global names use traitN: prefixes while local names remain in model metadata. Complete prediction order is unit + (trait - 1) * n_units in one-based R notation.

gs_random_regression() accepts a caller-supplied observation-by-basis matrix, a unit covariance, and a covariance over trait-by-basis axes. Latent order is unit + (((trait - 1) * n_basis + basis - 1) * n_units). The package does not generate polynomial or spline bases. gs_longitudinal_residual() exposes independent occasions, compound symmetry, discrete AR(1), continuous-time AR(1), and heterogeneous discrete AR(1). Its correlation argument is the native unconstrained transformed coordinate: logistic for compound symmetry, tanh for discrete AR(1), and log decay for continuous AR(1).

Both families use the existing persistent gs_model/gs_fit handles and fresh workspaces. Exact-dense and matrix-free stochastic routes are explicit. Prepared genomic bases require matrix-free inference and form no dense GRM. predict() returns observed contributions, complete_predictions retains complete multivariate targets, gs_ranef() returns longitudinal latent coefficients, and gs_selected_prediction() performs bounded one-based target inference. Covariance summaries, genetic correlations, delta-method uncertainty when requested, Monte Carlo precision, and solver diagnostics stay separate.

Current limitations

This interface does not provide formulas, automatic intercepts or contrasts, data-frame alignment, BIM/FAM parsing, allele flipping, genome-build conversion, arbitrary sparse coefficient precision, automatic route selection, selected-inverse uncertainty, or downloadable binary distribution. A pedigree and genomic term in one hybrid model are separate covariance components and are not an automatic single-step H or $H^{-1}$ construction. Windows local installation is a maintainer/development workflow requiring access to the sibling native repositories. This persistent-model interface does not expose the specialized sparse pedigree MME fitters described in the following sections; those have separate public R constructors and fit functions. Complete-grid prediction is provided for cross-sectional multivariate terms; longitudinal prediction targets are explicit latent basis coefficients rather than an implicit time grid.

For current build requirements use package metadata. Package loading queries the installed libraries and requires compatible ABIs and adapter linkage. These checks do not promise a stable native C++ ABI.

Multivariate sparse-pedigree boundary

gs_multitrait_pedigree_records() represents only observed scalar traits and groups them into residual-independent observation units. gs_multitrait_pedigree_fit() supports fixed-SPD sparse-PCG MME/BLUP and fixed-parameter exact or stochastic REML. Covariance optimization supports complete deterministic coefficient traces under an explicit dimension limit or the qualified fixed-common-probe stochastic PCG route. The result distinguishes covariance-optimizer status, optional likelihood availability, BLUP finalization, logical PCG systems and recovered retries; an imprecise final SLQ estimate remains diagnostic and does not invalidate converged covariance or retained BLUP. See the dedicated interface note.

Longitudinal sparse-pedigree boundary

gs_longitudinal_pedigree_records() represents one observed scalar record per row with a pedigree animal, supplied fixed and random-regression basis rows, time for reporting, and a residual-group label. Missing occasions are absent. gs_longitudinal_pedigree_fit() supports portable fixed-covariance PCG BLUP or bounded exact-trace positive-interior PARDISO REML, returning all-animal latent coefficients plus observed and supplied prediction trajectories. The direct route has explicit coefficient, trace, batch and factor-memory limits and no silent PCG fallback. See the dedicated interface note.

All three pedigree REML controls expose exact_trace_method = "coordinate" | "selected_inverse" and maximum_selected_inverse_workspace_bytes. Coordinate traces remain the default. Selected inversion requires an explicitly resolved serial CHOLMOD factor, performs one derivative-support inverse request per full accepted-state trace evaluation, and leaves AI derivative-response solves unchanged. PARDISO, including automatic resolving to PARDISO, returns not_supported; neither the backend nor trace method falls back. The common sparse_direct$selected_inverse_trace diagnostics report requested/resolved method, capability and message, derivative/request/factor entry counts, workspace and projected peak bytes, calls/successes/failures, coordinate systems avoided, support threads, numerical validity, and retained-factor validity.

For a route-by-route numerical and statistical overview, see gsolve and greml capabilities. The bounded read-only dataset campaign and its explicit nonconvergence outcomes are recorded in qgdata qualification.

Packed genomic and initializer controls

gs_control() accepts:

initializer = "provided" (default) or "rhe";
rhe_probe_count, rhe_seed, and rhe_probe_batch_size;
genomic_solver = "covariance_pcg" (default), "marker_batch", "marker_gs", or "auto_iteration";
marker_threads, marker_warm_start, and final_likelihood;
marker_score_ai="direct" (default), plus the independent execution control marker_coefficient_reuse ("automatic", "off", or "required") and marker_coefficient_cache_memory_limit;
marker_block_size, gram_cache_policy, and gram_cache_memory_limit;
factor_cache_policy and factor_cache_memory_limit for the distinct per-penalty-state Cholesky cache;
marker_maximum_sweeps, marker relative/absolute tolerances, and marker_preconditioner_sweeps.

It also accepts acceptance_policy = "likelihood_monotone" (the unchanged default) or experimental "determinant_deferred", plus maximum curvature-scaled score, maximum Newton decrement, stable-active-set iterations, and independent verification probe count, seed, SE multiplier, and absolute tolerance.

affine_spectral_reuse explicitly selects disabled, automatic, or required reuse for an eligible one-kernel plus identity-residual covariance. precision = "staged" uses the existing probe/PCG/Lanczos values as the initial state and accepts named working_precision and final_precision overrides. Precision tightens monotonically, common candidate comparisons use the same state, and final precision is required before convergence. Rejected candidates are likelihood-only and omit score, AI, prediction, and uncertainty work. Diagnostics report eligibility and fallback, spectral setup/reuse (including cumulative setup time), purpose-specific evaluation counts, avoided score/AI work, escalation reasons, and separate RHE, REML, and final-evaluation time.

The returned initializer$rhe object separates raw and constrained estimates, boundary adjustments, moment rank/conditioning, probe uncertainty, seed, batching, and covariance applications. These are initialization diagnostics, not whole-gs_model_fit() timing. Evaluation at the RHE start, optimization, and an optional final marker prediction remain separate stages.

marker_solver reports requested and resolved route, fallback reason, parameter-state preparation/construction counts, fixed symmetric sweep count, marker visits, decode passes, applications, processed RHS columns, factor lifecycle, and packed/cache/workspace memory. marker_gs resolves to marker_gs_preconditioned_pcg: covariance PCG with the live fused-coordinate gbits preconditioner, not a direct marker-space inverse. It rejects ineligible models and outer threads above one. covariance_pcg is unchanged; auto_iteration conservatively reports a covariance-PCG fallback because no size-based heuristic is validated. Total genomic prediction is y-F beta-V_e Py; component marker effects require a separate explicit result request. Trait-coupled multivariate marker updates and MCEM remain unsupported.

For marker_batch, the response and complete deterministic trace-probe prefix are one gbits call and the component-Py information systems are a second call. Complete system columns belong to disjoint workers; decoded packed blocks are shared read-only, while fixed/marker coefficients, residuals, convergence, and warm state remain independent. Multiple genomic components retain their own packed cache, selection, transform, weights, normalization, penalty, and coefficient block but share one observation residual in each system.

marker_solver additionally reports requested/resolved marker threads, response/probe/information system counts, warm-started systems, sweeps, maximum normal residual, shared decode/category/residual passes, timing, likelihood availability, and memory. Staged covariance tolerances also tighten the marker-MME solve, deterministic prefix growth retains old warm columns, and new probes start at zero.

For the default direct score/AI estimator, automatic marker-coefficient reuse uses the MME scaling $a_c = \theta_c D_c W^{\mathsf T}Pb/s_c$ to reproduce every original probe realization and response quadratic without applying K_c. $W^{\mathsf T}Z$ is a double-precision, prepared-source/probe-identity cache; compatible prefix growth appends only new columns, while independent verification probes use a separate entry. The configured memory bound causes a diagnosed direct fallback in automatic mode or an explicit failure in required mode. Zero and guarded near-zero component variances always retain direct boundary handling. Shared-source component values are reconstructed by one multi-RHS packed application. marker_solver$coefficient_reuse exposes policy resolution, per-component reasons, cache lifecycle and bytes, reduction/reconstruction time, avoided/direct applications, and the cache/reconstruction traversals included in overall totals. This control is not combined with the experimental identity_optimized trace estimator and does not change likelihood work.

final_likelihood=FALSE is supported for a fixed-parameter score/AI call and skips determinant work. Optimization still requires a final likelihood. With acceptance_policy="determinant_deferred", iterative accepted, rejected, and convergence-verification states compute projected score/AI without SLQ. Backtracking reduces the maximum of curvature-scaled projected score and the free-system Newton decrement. Fixed and active-boundary components do not enter the free AI solve.

Convergence also requires KKT score, parameter change, linear/MME residuals, score SE, active-set stability, a larger configured final probe prefix, and a distinct-seed independent verification set. complete_probe_basis=TRUE may bypass fresh probes. Exactly one final likelihood-only generic covariance-PCG/SLQ evaluation then fills the missing determinant fields without repeating packed Py, trace, score, AI, or probe batches. Failure of this final evaluation fails the fit.

This policy requires a safeguarded nonadaptive matrix-free route. It is available for single-trait, complete-record multivariate, balanced structured, and general observed-record longitudinal models. Sparse MME, hybrid, legacy indexed/sparse longitudinal, and incompatible facade models reject it explicitly. diagnostics(fit)$evaluation_work separates accepted/rejected score/AI, convergence and independent verification, likelihood-only finalization, SLQ applications, time, requested/resolved packed threads, processed and warm systems, symmetric sweeps, and packed decode passes. Score/AI rejected trials can be more expensive than likelihood-only trials, so performance gains depend on the workload.

Multivariate RHE and genomic values

gs_multivariate_rhe_blup(Y, fixed, covariance, control) accepts a complete numeric phenotype matrix, a common full-rank fixed-effect matrix, and a named list of non-residual gs_covariance components. It rejects missing/non-finite phenotypes and incompatible individual dimensions rather than imputing or dropping rows. Trait, sample, fixed-effect, and component names are retained.

gs_multivariate_rhe_control() configures the common probes, packed probe-batch size and threads, moment rank/condition thresholds, explicit PSD adjustment, block or independent multiple-RHS PCG, block numerical-rank controls, true residual verification, and the documented unit-diagonal fallback. The default block rank tolerance is 1e-12 and the bounded restart cap is 64; these govern numerical deflation, not the configured residual convergence tolerance.

The returned rhe section reports status, moment rank/condition and collinearity, common-probe Monte Carlo diagnostics, raw component/residual matrices, their eigenvalues and PSD status, adjusted prediction matrices, and every adjustment. The prediction section reports B_hat, alpha, every named U_c, U_total, PCG/operator/preconditioner work, final true residuals, block regularization/fallbacks, and purpose-labelled packed traversals.

The RHE-to-BLUP helper remains a one-pass initializer/prediction workflow. The separate complete-record gs_multivariate_model() and gs_model_fit() route provides multivariate REML and likelihood evaluation; this does not add REML semantics to gs_multivariate_rhe_blup().

Multivariate packed sufficient-statistic reuse

For a matrix-free complete-record multivariate model, the existing marker_coefficient_reuse control selects only execution. "automatic" uses packed reuse term-by-term, "off" forces direct covariance applications, and "required" errors when a packed relationship term cannot be resolved. The log-Cholesky covariance parameters, derivative matrices, stochastic probes, score, AI, and determinant route are identical between reuse policies.

Packed terms sharing a gs_bed() source share $W^{\mathsf T}Z$, $W^{\mathsf T}PZ$, and $W^{\mathsf T}PY$ transposes. Marker weights are reduced in marker space and all component K_c PY matrices are reconstructed by one fused output traversal per source. The fit diagnostics expose coefficient_reuse cache hits, misses, extensions, invalidations, bytes, columns, $W^{\mathsf T}Q$ source, reduction/reconstruction time, fallback reasons, traversals by purpose, contiguous SNP runs, decoded records, fused RHS widths, and requested/resolved packed threads.

Reuse requires complete finite trait records and a compatible packed prepared source. Generic components remain on the direct route. Experimental acceptance_policy="determinant_deferred" is available for these complete-record multivariate matrix-free fits. It uses common probes, score-only rejected trials, retained-state AI completion, separate numerical/stochastic convergence fields, and one final likelihood-only evaluation after verified convergence. The likelihood-monotone default is unchanged. Singular covariance boundaries, missing patterns, longitudinal incidence, PEV/reliability, and marker-effect output remain unsupported.

Balanced structured longitudinal pipeline

gs_longitudinal_component(..., constraint = "positive") explicitly constrains the coefficient covariance Q, not the relationship operator K. The base’s scalar constraint belongs to gs_model() and is not a constraint on Q here. The default remains SPD/log-Cholesky. Balanced random blocks may explicitly use constraint = "may_be_zero", accepting PSD Q including zero and partial rank. Use gs_longitudinal_control(acceptance = "likelihood_monotone", ...) with an explicit complete REML probe matrix, complete_probe_basis = TRUE, and lanczos_depth >= length(Y). These probes must actually satisfy $ZZ^{\mathsf T}/p=I$; the flag declares that contract. Ordinary stochastic PSD optimization is outside this first version. RHE raw estimates are preserved; their PSD start mapping uses no positive eigenvalue floor. Residuals remain "positive", and "fixed" blocks and observed-record PSD optimization are rejected. Returned PSD REML parameters are lower-column-packed Q entries, with off-diagonals scaled by sqrt(2); SPD parameters retain their existing log-Cholesky meaning.

The balanced interface uses one model object for initialization, refinement, and prediction:

time <- c(0, 1, 2, 3)
X <- cbind(intercept = 1, exposure = exposure)
F <- cbind(intercept = 1, time = time)
B <- cbind(intercept = 1, time = time)

component <- gs_longitudinal_component(
  gs_genomic(bed, "genomic", marker_weights, sum(marker_weights)),
  basis = B,
  Q = matrix(c(0.45, 0.08, 0.08, 0.24), 2, 2),
  name = "genomic"
)
model <- gs_balanced_longitudinal_model(
  Y, X, F, list(component),
  residual = gs_balanced_longitudinal_residual("homogeneous", 0.32)
)
control <- gs_longitudinal_control(
  marker_coefficient_reuse = "automatic",
  acceptance = "determinant_deferred",
  maximum_iterations = 15L,
  seed = 20260826
)
fit <- gs_longitudinal_fit(
  model, control, predict = TRUE,
  B_new = list(cbind(1, c(-0.5, 1.5, 2.5))),
  F_new = cbind(1, c(-0.5, 1.5, 2.5))
)

Y must be a complete finite n by T matrix. X, F, and every supplied B_c must have full column rank; bases are never generated internally. The model is $\operatorname{E}(Y) = X C F^{\mathsf T}$ with covariance action

\[ V(A) = \sum_c K_c A B_c Q_c B_c^{\mathsf T} + \operatorname{residual}(A). \] $B_c=I_T$ and $F=I_T$ are the exact complete-record multivariate reduction.

fit$rhe distinguishes raw and adjusted Q_c, raw and adjusted residuals, adjustment records, the mapped REML initial parameters, moment-system and fixed-design diagnostics, and packed traversals. fit$reml returns final coefficient/grid covariances, residual, C, PY, score/AI, score/Newton stationarity diagnostics, likelihood-work counters, and coefficient-reuse cache/traversal diagnostics. fit$blup is computed once when requested and returns latent coefficients, observed-grid component trajectories, totals, fixed/fitted trajectories, true residuals, and packed traversals.

The new control defaults to experimental determinant-deferred acceptance. Numerical score stationarity does not require an additional accepted micro-step; line-search exhaustion alone cannot converge. Stochastic precision remains a separate flag. Successful convergence performs no iterative likelihood or SLQ work and exactly one final likelihood-only evaluation. Likelihood-monotone is available explicitly for reference validation.

The balanced route supports homogeneous independent and complete-record unstructured residuals. General missing/irregular records use the interface below.

Missing and irregular longitudinal records

Each row is one observed scalar response. Missing occasions are omitted rather than represented by NA:

subject <- c("a", "a", "b", "b", "b", "c", "c")
time <- c(0, 1.4, 0.2, 0.9, 2.1, 0, 1.8)
y <- c(1.2, 1.7, 0.8, 1.1, 1.6, 1.0, 1.5)
fixed <- cbind(intercept = 1, time = time,
               treated = subject %in% c("b", "c"))
B <- cbind(intercept = 1, time = time)

component <- gs_longitudinal_component(
  gs_genomic(bed, "genomic", marker_weights, sum(marker_weights)),
  basis = B,
  Q = matrix(c(0.45, 0.08, 0.08, 0.24), 2, 2),
  name = "genomic"
)
residual <- gs_record_longitudinal_residual(
  group = ifelse(time < 1, "early", "late"),
  variance = c(early = 0.25, late = 0.38)
)
model <- gs_longitudinal_records_model(
  y, subject, fixed, list(component), time, residual
)
fit <- gs_longitudinal_fit(
  model,
  gs_longitudinal_control(
    acceptance = "determinant_deferred",
    marker_coefficient_reuse = "automatic",
    maximum_iterations = 15L,
    seed = 20260827
  ),
  predict = TRUE,
  prediction_records = list(
    subject = c("a", "b"),
    fixed = cbind(1, c(2.2, 2.6), c(0, 1)),
    basis = list(cbind(1, c(2.2, 2.6)))
  )
)

Subject and record order are preserved. Component bases can have different column counts. Residual groups are diagonal and each estimated group must have at least one record. The result separates raw, adjusted and REML-initial RHE values; BLUP returns n by r_c latent coefficients, observed-record component/total/fixed/fitted values and residuals, and requested prediction records. Pedigree/sparse MME, multiple responses per record, correlated AR(1)/continuous-time residuals, exact singular boundaries, PEV/reliability, marker effects, and single-step models remain unsupported.