Longitudinal sparse-pedigree REML and BLUP
gs_longitudinal_pedigree_records() stores only observed scalar longitudinal records. Each row supplies an animal, phenotype, reporting time, fixed-design row, caller-supplied random-regression basis row, and residual-group label. Missing occasions are absent rows. Records may be irregular, repeated, arbitrarily ordered, and uneven across animals; pedigree-only animals remain in the fitted domain. Character residual groups use deterministic lexical order and integer groups use deterministic numeric order, so record permutation does not change covariance-parameter positions.
ped <- gs_pedigree_data(
individual = c("dam", "sire", "offspring", "relative"),
sire = c(NA, NA, "sire", "sire"),
dam = c(NA, NA, "dam", "dam")
)
animal <- c("dam", "dam", "sire", "offspring", "offspring")
time <- c(0, 1.4, 0.2, 0.7, 2.1)
y <- c(1.1, 1.8, 0.9, 1.4, 2.2)
records <- gs_longitudinal_pedigree_records(
ped, animal, y, time,
basis = cbind(intercept = 1, time = time),
fixed = cbind(intercept = 1, time = time),
residual_group = ifelse(time < 1, "early", "late")
)
start <- gs_longitudinal_pedigree_start(records)
control <- gs_longitudinal_pedigree_reml_control(
linear_solver = "sparse_direct", direct_backend = "pardiso",
requested_threads = 8,
exact_trace_method = "coordinate",
maximum_direct_coefficient_dimension = 5000,
max_exact_trace_systems = 5000
)
fit <- gs_longitudinal_pedigree_fit(
records, start, covariance = "optimized", control = control,
prediction = list(
animal = c("offspring", "relative"), time = c(1.5, 1.5),
basis = cbind(1, c(1.5, 1.5)),
fixed = cbind(1, c(1.5, 1.5))
)
)The deterministic non-RHE start fits the supplied fixed design, calculates \(s_r^2=\frac{\mathrm{RSS}}{N-\operatorname{rank}(X)}\) and mean basis magnitude \(m_b=\frac{1}{N}\sum_j\lVert b_j\rVert^2\), then sets \(Q_0=\frac{s_r^2}{2m_b}I\) and each residual-group variance to \(s_r^2/2\). User-supplied SPD Q and positive residual starts override it.
Both solvers use one canonical animal-major sparse MME. linear_solver = "pcg" is portable and is the large fixed-covariance BLUP route. "sparse_direct" is an optional moderate-problem route with one reused PARDISO or CHOLMOD symbolic analysis, exact retained-factor likelihood, and bounded complete coordinate traces. Automatic selection prefers PARDISO; CHOLMOD and support work are serial. Explicit direct selection never falls back: an unavailable backend returns not_supported and the gsolve message. PCG factors the bounded fixed-effect preconditioner block and one basis-rank block per pedigree animal, and verifies the true MME residual.
Coordinate traces remain the default. Explicit CHOLMOD supernodal selected inversion replaces only the exact coordinate trace systems. Deterministic positive-definite planning compiles derivative support once, then reuses one trace plan, derivative mapping, and selected-inverse workspace across covariance states. The result reports calls, selected and filled-factor entries, workspace/request/output and projected peak bytes, systems avoided, and serial selected-inverse support threads. AI derivative-response systems, exact determinant work, optimizer acceptance, and retained BLUP finalization are unchanged. CHOLMOD orchestration and support remain serial, while selected-inverse dense kernels use an independently requested oneMKL thread count. PARDISO selected inversion, including automatic when it resolves to PARDISO, is explicitly not_supported without fallback after symbolic capability resolution and before numerical factorization. PCG remains a separate fixed-covariance route with zero direct work.
The result reports starting/final covariance, fixed effects, animal-by-basis latent coefficients, observed genetic trajectories, fitted values, residuals, Py, optional prediction trajectories, covariance/optimizer diagnostics, and complete factorization/trace/AI work. Successful optimization reuses the accepted response coefficients for one BLUP finalization and performs no additional response solve. The shared typed serializer translates these diagnostics without numerical work. In the direct diagnostics, sparse_direct$requested is a logical direct-route flag, whereas requested_backend and resolved_backend are backend strings.
Resource preflight bounds coefficient dimension, coordinate systems, RHS batch width, workspace, and factor memory. Construction and PCG diagnostics expose the single record pass, zero optimizer record revisits, stored sparse-statistic values/nonzeros, fixed and animal preconditioner blocks, RHS width, and true residual. CHOLMOD selected-inverse trace and retained-factor prediction-error paths are available through this API. Exact PEV/reliability requires retained selected-inverse state.
PARDISO selected inversion, MUMPS, covariance boundaries, RHE, serial residual correlation, multiple relationship components, genomic models, marker effects, and ssGBLUP are not supported by this interface. The prediction-error contract supports bounded PARDISO selected_solve coefficient and trajectory uncertainty from retained factors. This does not add PARDISO selected inversion or PCG exact PEV.