Formula and variance-component compiled-MME REML

For the new data-frame shorthand, optional starts, fixed solves and genetic random regression, begin with Pedigree analysis. Both formula and formulas now use the same general adapter. This page documents the specialized compiled-MME route, explicitly requested with a gs_fit_control() object, complete explicit covariance starts, and task = "reml" (or greml() with those same controls). Its result layout remains specialized and does not support the gsolve(fit) handoff. Supplying starts alone does not select it. General fits accept named controls for the selected pedigree or component fitter. Compiled controls now default to CHOLMOD selected inversion.

The data-frame gs_pedigree() constructor described here is distinct from gs_model_pedigree(), the ordered-parent constructor used by the persistent-model API. The two signatures are not interchangeable and do not dispatch implicitly.

greml() is the estimation entry point; gfit() remains available with task = "reml". Model incidence comes from formulas, pedigree relationships come from gs_pedigree(), and explicit positive-interior covariance starts come from named gvc() components. The interface lowers one immutable model to greml’s public compiled-MME plan, creates one numerical session, optimizes, and finalizes BLUPs from the retained accepted response state.

PED <- gs_pedigree(pedigree, id = "id", sire = "sire", dam = "dam")

fit <- gfit(
  BW ~ sex + batch + (1 | id),
  records,
  list(
    animal = gvc("id", "BW", PED, start = 1),
    residual = gvc("residual", "BW", start = 1)
  ),
  control = gs_fit_control(direct_backend = "cholmod")
)

For multiple traits, formulas is a named list and each component start is an SPD matrix in the declared trait order. Each trait has a separate fixed-effect block. Scalar records are ordered by formula-list order and original data-row order. An NA response is absent; observed traits from the same data row form one complete or incomplete patterned residual unit. Missing fixed covariates, nonfinite observed responses, and missing pedigree identifiers are errors.

Runnable single-trait example

From the repository root, with gsuite installed, run:

source("tools/examples/run_single_trait_pedigree_example.R")

The complete example starts with two ordinary data frames: data has one observation per row (id, BW, factor sex, continuous covariate x); pedigree has one animal per row (id, sire, dam). It reuses the small deterministic BW formula-test illustration, not a random draw from the fitted pedigree model. Its artificial sex labels vary across records of an animal; real data should use consistent animal sex. These toy values illustrate the API, not biological inference or variance recovery.

The formula BW ~ sex + x + (1 | id) supplies fixed effects and animal incidence. A separate gs_pedigree() relationship and named animal and residual components supply covariance structure and positive starting variances. Under this model, animal covariance is the pedigree relationship times the animal variance; residuals are independent with a common variance. Repeated observations share the animal effect, with no additional permanent environment effect. Starts are not fixed variances.

IDs, not row positions, align the tables. Pedigree rows may be reordered; unknown parents use NA, while known parents must be listed. Observation IDs absent from the pedigree are errors. Animal a8 has no observations but remains in the pedigree and receives a conditional animal prediction from relatives.

summary(fit) displays covariance estimates and convergence, coef(fit) returns fixed effects, and the public fit$random_effects$animal table contains IDs and animal BLUPs (not full predicted body weights). fit$converged and fit$convergence_reason expose the stopping result. Prediction tables follow pedigree order, so compare them by ID. Prediction for new data and prediction uncertainty are not supported on this formula route.

For development verification only, select the existing isolated installation before sourcing the same example:

.libPaths(c(file.path("build", "r-library"), .libPaths()))
source("tools/examples/run_single_trait_pedigree_example.R")

Formula lowering

The only random grammar is exactly (1 | variable). A scalar formula has one or two distinct terms, with at most one pedigree component and one independent identity component. Multitrait formulas retain one shared pedigree term. Random slopes, more than two terms, duplicate terms, ||, nesting, crossing, functions, interactions, and expressions on the grouping side are rejected. After removing the random terms, base R model.frame() and model.matrix() build the fixed design and preserve coefficient names, contrasts, factor levels, and the formula environment.

Random incidence is never materialized as an R matrix. Each observed scalar record supplies a pedigree entity ID, trait/member ID, residual unit and coordinate, and a unit random-intercept weight directly to the compiled-MME specification. This first version may materialize a dense R fixed-effect model.matrix() and therefore targets moderate in-memory data frames.

Components, resources, and results

Scalar fits accept one identity component, one pedigree component, or one of each, plus one homogeneous residual component. Formula order determines the random coefficient blocks and covariance-parameter order. Pedigree levels use external pedigree order, including unphenotyped animals. Identity levels use the sorted observed character or integer domain. Every record contributes a unit weight directly to each present component, including the data-induced pedigree/identity cross block; no R incidence matrix is built.

For example, repeatability and a common-environment model use the same bounded composition rule:

ID <- gs_identity()
repeatability <- gfit(
  y ~ sex + (1 | animal_id) + (1 | permanent_id), records,
  list(
    animal = gvc("animal_id", "y", PED, start = 0.8),
    permanent = gvc("permanent_id", "y", ID, start = 0.4),
    residual = gvc("residual", "y", start = 0.3)
  ),
  control = gs_fit_control(direct_backend = "cholmod")
)

common_environment <- gfit(
  y ~ sex + (1 | animal_id) + (1 | environment), records,
  list(
    animal = gvc("animal_id", "y", PED, start = 0.8),
    environment = gvc("environment", "y", ID, start = 0.4),
    residual = gvc("residual", "y", start = 0.3)
  ),
  control = gs_fit_control(direct_backend = "cholmod")
)

Starts remain mandatory finite positive scalars or finite symmetric positive-definite matrices. Covariance boundaries and automatic rank repair are not attempted.

gs_fit_control() maps backend, direct/support threads, exact trace route, optimizer bounds, coefficient dimension, factor memory, selected-inverse workspace, exact trace systems, RHS workspace, and linear-solve controls directly to the public greml session. A failed preflight does not silently change backend, trace route, batch width, or solver. Diagnostics report the required and configured RHS bytes, peak RHS width, resident coefficient vectors, workspace category, limit result, overflow result, and zero numerical work after rejection.

PARDISO uses its own requested internal thread count and coordinate-exact traces. CHOLMOD uses serial support work and selected inversion. Generic support execution remains independent; an unsupported request is returned as not_supported before solver work. No GNU OpenMP layer is introduced.

The returned gs_fit owns covariance estimates, fixed effects, animal BLUPs in external pedigree order, caller-record mappings, optional fitted values, residuals and Py, score/standard-AI diagnostics, optimizer history, resource preflight, work counters, and retained-state diagnostics. print(), summary(), coef(), fitted(), residuals(), logLik(), and nobs() are provided. Record methods fail clearly when record values were omitted.

Genomic relationships, random slopes, correlated random components, more than one pedigree or identity component, multitrait repeatability, heterogeneous or serial residual covariance, covariance boundaries, PEV/reliability, prediction from new data, disk-backed data, formula dependencies, Bayesian tasks, JSON, Python, and R-side public plan/session objects remain deliberately outside this compiled-MME formula interface.