rrm: Reduced-rank covariance structure

In sommer: Solving Mixed Model Equations in R

Reduced-rank covariance structure

Description

Reduced-rank covariance structure. The returned dimensionless covariance shape is intended for use inside vsm, which supplies the single overall variance scale.

Usage

rrm(x, k = 1L, loadings = NULL, fixed = NULL)

Arguments

x	Variable defining q covariance levels.
k	Reduced rank, satisfying 1 <= k < q.
loadings	Optional finite q by k starting loading matrix. The leading k by k block is made lower triangular and must have positive diagonal entries. Loadings are divided by their first diagonal entry to remove the internal scale.
fixed	Logical vector of length qk-k(k-1)/2, controlling the triangular loading entries in column-major order. The first entry is always fixed: the reference loading is one, and vsm() estimates the overall variance scale. Other fixed entries constrain relative loadings.

Details

The reduced-rank structure uses

K=\Lambda\Lambda^{\mathsf T},\qquad \Lambda_{11}=1,

and vsm() supplies the overall variance scale:

G=\sigma^2 K.

All specific variances are fixed at exactly zero; no diagonal remainder or numerical nugget is added. This is the same covariance family as ASReml's default rr(x, k) model. Individual loadings can differ by sign or rotation, so compare reconstructed covariance matrices when comparing fits.

The leading k\times k loading block is lower triangular for rotational identification and its diagonal entries are positive. There are qk-k(k-1)/2-1 free shape parameters plus the overall scale. Because k<q, the random-effect covariance is positive semidefinite and singular by design.

The Henderson engine fits this structure through latent factors rather than inverting the singular environment covariance. The direct-inversion formulation (control$formulation="direct") also supports it, provided the total observation covariance is positive definite. rrm() is intended for random effects, not residual covariance structures.

Native latent-factor Henderson fits support control$solve$method="ldlt" and control$solve$method="cholmod" (including automatic selection). They use loading-aware likelihood scores and average-information updates in the original loading/variance coordinates. The loading derivatives account for the parameter-dependent latent design; the singular environment covariance is never inverted. control$proposal$stepWeight damps covariance-parameter proposals, while control$update$blendWeight blends EM/GEM or scalar PX-GEM with average-information update directions, selected through control$update$recovery, as in ordinary mmes() fitting. Parameter standard errors use the unblended average-information matrix at the final parameter values. Prediction errors are mapped back from the latent coefficient equations to the original random effects.

The likelihood trace terms are computed exactly by sharing random-design and residual-block contractions across covariance parameters. Bounded panels of Henderson solves avoid constructing a dense record-by-record covariance or projection matrix. Random trace solves use the loading span (rank k), rather than every environment direction. Reducing the coefficient dimension does not guarantee a faster fit for every dataset.

Section-specific spatial residuals supplied by dsumm() are evaluated on their observed residual blocks, including grids with missing cells and unlisted sections with independent residuals. Residual covariance derivatives and trace contractions never expand the full section-by-grid covariance. engineDiagnostics$residualEvaluation reports "observed-blocks" and engineDiagnostics$largestResidualBlock records the largest observed block dimension used by the native RR AI evaluator.

Likelihood and AI calculations share the selected residual precision operator and the live coefficient factorization, including grid residuals with missing-cell absorption and eligible block-Schur solves. No separate coefficient factorization or dense observed-block precision is built for AI. Random trace panels use the latent Henderson identity to cancel relationship covariance actions. Grid residual traces use block-local projection panels rather than repeated full-observation precision applications. engineDiagnostics$residualEngine identifies the shared residual engine, and sharedConditionalOperators confirms its use. aiElapsedSeconds gives elapsed timings for the final AI evaluation's setup, random traces, residual traces, and information matrix.

The positive-specific-variance Woodbury Schur formulation used for fam() is unnecessary for exact RR: the specific-effect blocks are absent, not inverted at zero. The smaller latent coefficient system still uses the ordinary LDLT/CHOLMOD and eligible block-Schur dispatch. With control$diagnostics$verbose=TRUE, the native fit reports its RR path, selected coefficient solver, and then the progress table. Each accepted iteration reports likelihood, wall-clock time, cumulative CPU seconds, the scheduled EMweight, accepted-step trustAlpha and lineAlpha, and score norm. The step scales describe the update that reached the current row (both are one at initialization); lineAlpha includes control$proposal$stepWeight and any backtracking. Boundary-restraint and notPD columns are omitted because the latent optimizer does not use those ordinary covariance-update diagnostics. Numeric weights, step scales and CPU times are retained in engineDiagnostics$iterationProgress.

For native Henderson fits with weights W, let U=\mathrm{chol}(W) in observation order. The effective residual precision is U R_0^{-1} U^{\mathsf T}, where R_0 is the unweighted residual covariance. The score and AI use derivatives of the corresponding covariance U^{-\mathsf T}R_0U^{-1}.

The native path does not currently support PCG, vcc constraints, covariance-update acceleration, control$task$mode="solve", or control$representation$factorScores augmentation. These combinations produce explicit errors. The stored coefficient matrix C describes the latent equations; coefficientMap maps latent coefficients to the original effect space. Use predict(), postPEV(), and the fitted uPevList for prediction uncertainty rather than directly indexing C by original effect indices.

Earlier versions used \Lambda\Lambda^{\mathsf T}+I_q, estimating a nonzero common specific variance on the fitted scale. To fit that former covariance family, use fam(x, k, variance="homogeneous").

Value

A list containing the incidence/design matrix in Z and a compiled CovarianceFactor v2 descriptor in covFactor, for use inside vsm.

See Also

fam, ownm, vsm.

Examples

## Not run: 
vsm(rrm(environment, k=2), ism(genotype))
# Former rrm() model: an estimated common specific variance
vsm(fam(environment, k=2, variance="homogeneous"), ism(genotype))

## End(Not run)

sommer documentation built on Oct. 24, 2026, 1:08 a.m.