This vignette is written like a short book. Each chapter builds on the previous one, and all chapters share the same notation so formulas can be compared directly across chapters:
vsm().Table of contents
ai_mme_sp2()\newpage
mmes(fixed = y ~ 1, random = ~ vsm(...), rcov = ~ vsm(...), data = ...)
fits
$$ y = X\beta + \sum_{i=1}^{r} Z_i u_i + e, \qquad u_i \sim N(0, G_i), \quad e \sim N(0, R) $$
random=~... entry.sommer never treats \(G_i\) or \(R\) as free unconstrained matrices. Every one is built from a small number of variance components \(\theta\) (Chapter 2), and REML estimates \(\theta\), not the matrix entries directly.
Rather than maximizing the marginal likelihood of \(y\) directly (which needs \(V = ZGZ' + R\), generally dense and huge), Henderson’s equations solve for \(\beta,u\) jointly using only \(R^{-1}\) and \(G_i^{-1}\):
\begin{pmatrix}X’R^{-1}y\ Z’R^{-1}y\end{pmatrix}, \qquad C = \begin{pmatrix} X’R^{-1}X & X’R^{-1}Z\ Z’R^{-1}X & Z’R^{-1}Z + G^{-1} \end{pmatrix} $$
with \(G^{-1} = \mathrm{blockdiag}(G_1^{-1}, \dots, G_r^{-1})\).
Why this matters for parameterization: every quantity sommer actually
needs — \(C\), its log-determinant, and the REML score/AI-matrix — is built
from inverses (\(R^{-1}\), \(G_i^{-1}\)) and their derivatives with
respect to \(\theta\), never from \(G_i\) or \(R\) themselves except to invert
them once per structure per iteration. This is why sommer’s internal name
for \(G_i^{-1}\) within a structure is Lambda (\(\Lambda = K^{-1}\)).
For a random term built with vsm(dsm(Env), usm(Trait), ism(Line)), sommer
constructs
$$ G = \sigma^2 \cdot K_{dsm} \otimes K_{usm} \otimes K_{ism} $$
vsm() call (its sigma2
argument).dsm, usm, ism, …) supplies a dimensionless
shape \(K\) with \(K_{11}=1\).This separation of scale (one parameter, always variance-like) from
shape (any number of parameters, some variance-like, some
correlation-like) is the central design decision explored in the rest of
this vignette, and it directly determines how derivatives are computed in
ai_mme_sp2().
| Variance-type parameter | Correlation-type parameter | |
|---|---|---|
| Example | sigma2, dsm() ratios, fam() specific variances |
csm()/ar1m() rho, usm()/corgm() off-diagonal terms |
| Effect on \(K\) | rescales one or few diagonal blocks | reshapes the entire matrix jointly |
| \(B_k=\partial K/\partial \eta_k\) | sparse / diagonal | dense |
| Average vs Expected Information | typically agree closely | can diverge — needs EM/AI blending |
| Practical consequence | fast trace/derivative shortcuts available | full dense chain-rule required |
\newpage
Imagine you are tuning two dials for a model:
If you let an optimization algorithm (the computer’s “hill-climbing” procedure that tries to improve the fit step by step) move these dials by any amount, in any direction, it will sooner or later try an illegal value: a negative variance, or a correlation bigger than 1. At that point the model is broken (the covariance matrix stops being a valid one — technically it is no longer positive definite, which is the mathematical way of saying “this could actually occur as a covariance of real data”).
sommer’s answer to this problem is always the same trick, used for every single covariance structure in the package:
Never let the optimizer touch the real dial (\(\theta\)). Instead, give it a fake dial (\(\eta\)) that can be turned all the way in either direction, forever, without ever breaking anything. Then use a fixed recipe (\(g\)) to translate the fake dial’s position into a legal value of the real dial.
Formally: \(\theta = g(\eta)\), where \(\eta\) ranges over all real numbers (\(-\infty\) to \(\infty\)) and \(g\) is built so that \(g(\eta)\) always lands in the legal range for \(\theta\), no matter what \(\eta\) is.
This is exactly the same idea used elsewhere in statistics: for example, when you model a probability (which must stay in \([0,1]\)) with a logistic regression, you don’t optimize the probability directly — you optimize an unconstrained “logit” score and convert it to a probability at the very end. sommer does the same thing, systematically, for every variance and correlation in every covariance structure it supports.
In sommer’s code, the collection of fake dials for one covariance structure
is called par (short for “working parameters”), and it always lives in
the CovarianceFactor descriptor returned by every constructor
(csm(), ar1m(), usm(), …). Three related pieces travel together:
par — the current numeric values of \(\eta\).free — a logical flag per entry: is this \(\eta_k\) estimated by REML, or
is the user holding it fixed (fixed=TRUE)?par_names — a human-readable label per entry, e.g. "rho" or
"variance_ratio[B]", used when you call summary() on a fitted model.REML (Chapter 4) works entirely with par/$\eta$. It never sees \(\theta\)
directly; \(\theta = g(\eta)\) is only computed when sommer needs to actually
build the covariance matrix \(K\), or when it reports results back to you in
a human-readable form.
The rule: \(\theta > 0\).
The recipe: \(\theta = \exp(\eta) = e^{\eta}\).
Why this works: no matter what real number \(\eta\) is — very negative, zero, very positive — \(e^\eta\) is always a positive number. As \(\eta\to-\infty\), \(\theta\to 0\) (but never reaches it); as \(\eta\to+\infty\), \(\theta\to\infty\). Every positive number is reachable by exactly one \(\eta\) (the inverse recipe is \(\eta=\log(\theta)\)).
A nice bonus: the “speed” at which \(\theta\) changes as you nudge \(\eta\) (its derivative, \(\partial\theta/\partial\eta = e^\eta = \theta\)) is proportional to \(\theta\) itself. This means equal steps in \(\eta\)-space are equal-percentage steps in \(\theta\)-space (multiplicative, like moving up a musical scale in octaves) rather than equal-amount steps — exactly the natural way variances tend to move (a variance going from 1 to 2 “feels” like the same size change as going from 10 to 20).
sommer almost never applies this to a raw variance, though. Instead it applies it to a ratio relative to a reference level — this detail matters because of scale confounding (see 2.6), so let’s define it now since you will see it in almost every structure in Chapter 3:
$$ \eta_j = \log!\left(\frac{\theta_j}{\theta_1}\right), \quad j=2,\dots,q $$
Level 1 is always the reference: \(\theta_1\) is fixed at 1 by definition, so only \(q-1\) working parameters are needed to describe \(q\) relative variances.
The rule: \(\theta \in (lo, hi)\), a bounded open interval — e.g. a correlation \(\rho\in(-1,1)\), or a compound-symmetry \(\rho\in(-1/(q-1), 1)\).
The recipe (general bounded interval):
$$ \theta = lo + (hi-lo)\cdot s(\eta), \qquad s(\eta) = \frac{1}{1+e^{-\eta}} \ \ (\text{the logistic/“sigmoid” function}) $$
The sigmoid function \(s(\eta)\) is the same S-shaped curve used in logistic regression: it always outputs a number strictly between 0 and 1, approaching 0 as \(\eta\to-\infty\) and approaching 1 as \(\eta\to+\infty\), and it is smooth and increasing everywhere. Stretching and shifting its output (\(lo + (hi-lo)\cdot s(\eta)\)) maps it onto any open interval \((lo,hi)\) instead of just \((0,1)\).
The recipe (symmetric interval \((-1,1)\)): when the bounds are exactly \(-1\) and \(1\) (as for a plain correlation), sommer often uses the equivalent but more compact hyperbolic tangent:
$$ \theta = \tanh(\eta) = \frac{e^{\eta}-e^{-\eta}}{e^{\eta}+e^{-\eta}} $$
which is just the logit recipe above rescaled — same S-shape, same guarantee (\(\tanh(\eta)\) is always strictly between \(-1\) and \(1\)), just centered at 0 instead of at \(½\).
Why bother with a whole interval instead of just “positive,” like variances? Because a correlation has two walls (an upper and a lower limit), while a variance only has one (a floor at 0, no ceiling). The sigmoid/tanh family is the natural tool whenever there are two walls to respect simultaneously.
A subtlety that matters in Chapter 4: near the middle of the interval, \(s(\eta)\) changes quickly as \(\eta\) moves — a small step in \(\eta\) makes a decent-sized change in \(\theta\). But near the two walls (as \(\eta\) gets very large in either direction), \(s(\eta)\) becomes almost flat: huge changes in \(\eta\) barely move \(\theta\) at all. This is deliberate and useful — it means that if REML tries to push a correlation toward its boundary, the working scale automatically resists further movement (the “gas pedal” stops doing much once you’re already flooring it), rather than needing an artificial rule bolted on afterward.
Families 1 and 2 work great for a single variance or a single correlation, one dial at a time. But some structures have many correlations that all have to be jointly consistent with each other — for example, an unstructured \(3\times 3\) correlation matrix has 3 pairwise correlations, but not every combination of 3 numbers in \((-1,1)\) is a valid correlation matrix (try \(\rho_{12}=0.9,\ \rho_{13}=0.9,\ \rho_{23}=-0.9\) — no data set can ever produce that).
So instead of transforming each correlation separately, sommer reparameterizes the whole matrix at once through a factorization — a recipe that builds the matrix out of simpler unconstrained ingredients in such a way that the result is automatically valid, no matter what those ingredients are. The two building blocks used throughout the package are:
usm() (Chapter 3.7) and, restricted to unit-diagonal
\(L\), by corgm() (Chapter 3.8).fam()/rrcm() (Chapter 3.9) and, in a
triangular/banded variant, by antem() (Chapter 3.10).The common thread across all three families: every recipe in this chapter is a one-way-safe translation from “any real number(s) at all” to “a legal \(\theta\)/$K$.” This is why sommer’s optimizer (Chapter 4) can take a Newton step of any size, in any direction, in \(\eta\)-space, and never produce an invalid model — the only thing that can still go wrong is that the fit gets worse, which is a completely different (and much easier to handle) problem than breaking the math.
Recall from Chapter 1 that a vsm() term is a Kronecker product
\(G = \sigma^2\cdot K_1\otimes K_2\otimes\cdots\). If every \(K_i\) were allowed
to carry its own free overall scale (instead of being pinned to \(K_{i,11}=1\)
via a reference level), the model would have infinitely many
\((\sigma^2, K_1,K_2,\dots)\) combinations that produce the exact same \(G\)
— for instance doubling \(\sigma^2\) while halving every entry of \(K_1\)
changes nothing observable. This is called non-identifiability: REML
would have no way to choose among the infinitely many equally-good answers,
and the optimizer would wander aimlessly or fail to converge.
sommer avoids this by a strict convention, applied without exception: every
covariance-shaping factor is normalized so its first diagonal entry is
exactly 1 (\(K_{11}=1\)), and the one and only free overall scale lives in
vsm()’s sigma2 argument. This is why, throughout Chapter 3, you will see
“$\theta_1$ is the reference, fixed at 1” repeated for structure after
structure — it is not a coincidence, it is the mechanism that keeps the
whole model identifiable.
\newpage
For each structure: the shape formula \(K(\eta)\) (normalized so \(K_{11}=1\)), the working parameters \(\eta\) and their transform (from Chapter 2), and whether it is fundamentally a variance-type or correlation/shape-type parameterization (this label is used again in Chapter 4 and in the matching worked example of Chapter 5).
ism(x) — identity$$K = I_q, \qquad \text{no free parameters.}$$
Simplest possible shape: independent, homoscedastic levels. Used as the
“main effect” incidence in most vsm() calls, or as a Kronecker factor when
no structure is desired on one dimension.
dsm(x, values, fixed, theta) / atm(x, levs, values, fixed) — heterogeneous diagonal variances$$ K = \mathrm{diag}(1, \theta_2/\theta_1, \dots, \theta_q/\theta_1), \qquad \eta_j = \log(\theta_j/\theta_1),\ j=2,\dots,q $$
Family-1 (exponential/log-ratio) transform, applied to every level except
the reference. Pure variance-type structure: \(q-1\) free log-variance
ratios, no correlation at all (off-diagonal is exactly 0). atm() is the
same construction restricted to a user-selected subset of levels (levs);
unselected levels get zero columns in the design matrix and are excluded
from \(\eta\).
csm(x, rho, fixed, variance, values) — compound symmetry$$ K_{ii}=1,\quad K_{ij}=\rho\ (i\neq j), \qquad lo=\frac{-1}{q-1} < \rho < 1 $$
Family-2 (bounded logit) transform: $\eta_\rho = \mathrm{logit}!\left( \frac{\rho-lo}{1-lo}\right)$.
variance="homogeneous" (default): only \(\eta_\rho\) is free — a single
correlation-type parameter.variance="heterogeneous": additionally appends \(q-1\) family-1
log-variance-ratios \(\eta_j=\log(\theta_j/\theta_1)\), so
\(K=\mathrm{diag}(\sqrt\theta)\, R\, \mathrm{diag}(\sqrt\theta)\) with \(R\)
the homogeneous compound-symmetric correlation. A mixed
parameterization: 1 correlation-type + \((q-1)\) variance-type parameters.ar1m/ar2m/ar3m/toeplitzm — autoregressive & general Toeplitz correlationAll parameterize an ordered correlation matrix through partial autocorrelations (PACF) \(\kappa\in(-1,1)^{order}\), family-2 (tanh) transform \(\eta=\mathrm{atanh}(\kappa)\), reconstructed to full correlations \(\rho_h\) (lag \(h\)) via the Durbin–Levinson recursion:
$$ K_{ij} = \rho_{|i-j|} $$
ar1m(): order 1, \(\rho_h=\rho^h\) (the classical AR(1) correlation), 1
correlation-type parameter.ar2m()/ar3m(): order ⅔, PACF vector of length ⅔ mapped through
Durbin–Levinson to \(\rho_h\) for all lags up to \(q-1\).toeplitzm(): order \(q-1\) (fully general stationary Toeplitz
correlation, no autoregressive-order restriction) — the most flexible
member of this family, still guaranteed valid for any PACF in
\((-1,1)^{q-1}\) (this is exactly the joint-validity guarantee discussed in
Chapter 2.5, achieved here through PACF rather than through a matrix
factorization).variance = "heterogeneous" (AR family only;
toeplitzm() is shape-only) the same way as csm(): append \(q-1\)
family-1 log-variance-ratios after the PACF working parameters.mam/ma1m/ma2m — moving-average covariance$$ K \propto \mathrm{Cov}!\left(e_t + \theta_1 e_{t-1} + \theta_2 e_{t-2}\right) $$
Working parameters are the raw MA coefficients themselves (\(\eta=\theta\), no transform at all) — unlike AR, any finite MA coefficients produce a valid covariance; invertibility of the MA polynomial is not required to define the covariance matrix, so no bounding transform is needed. This is the one “ordered” structure whose working parameter equals its natural parameter directly. Purely a correlation/shape-type structure (no variance option).
corgm(x, theta, fixed) — general correlation matrix$$ K = AA’, \qquad A \text{ unit-diagonal lower-triangular} \ (A_{ii}=1) $$
Family-3 (factorization) transform: the strictly-lower-triangular entries
of \(A\) are working parameters directly (no further per-entry transform is
needed — Chapter 2.5 explains why any real lower-triangular \(A\) with unit
diagonal automatically gives a valid \(K\)). \(q(q-1)/2\) purely
correlation-type parameters, no variance component at all (contrast
with usm() below, which also carries diagonal/variance information).
usm(x, theta, fixed) — general unstructured covariance$$ K = LL’, \qquad L \text{ lower-triangular}, L_{11}=1 $$
Family-3 (Cholesky factorization) transform: for each free entry of \(L\) below the diagonal — the off-diagonal entries are working parameters directly (\(\eta=L_{ij}\)), the diagonal entries go through the family-1 log transform (\(\eta=\log L_{ii}\), keeping \(L_{ii}>0\) hence \(K\) nonsingular). Total \(q(q+1)/2 - 1\) free parameters (the \(-1\) is \(L_{11}=1\), fixed by the Chapter 2.6 reference-level convention).
This is the maximally flexible valid shape for \(q\) levels — every possible covariance shape is reachable, at the cost of \(O(q^2)\) parameters. It mixes variance-type (log-diagonal) and correlation-type (off-diagonal) working parameters in a single structure.
fam(x, k, loadings, specific, fixed) / rrcm(x, k, loadings, fixed) — factor-analytic covariance$$ M = \Lambda\Lambda’ + \Psi, \qquad K = M / M_{11} $$
Family-3 (loadings + diagonal) transform:
A rank-$k$ correlation-type structure (loadings) plus a genuinely
variance-type component (specific variances) — a compact alternative to
usm() when \(q\) is large and a low-rank + diagonal approximation is
statistically reasonable (the standard “FA” model for many traits or many
environments). rrcm() is the same loading parameterization with
\(\Psi=I\) fixed (a pure reduced-rank approximation, no free specific
variances) — purely correlation/shape-type.
antem(x, order, beta, innovations, fixed) — antedependence$$ Ty = e,\quad \mathrm{Cov}(e)=D \text{ diagonal},\qquad K = T^{-1} D T^{-T},\quad D_{11}=1 $$
\(T\) is unit lower-triangular with regression coefficients
\(\beta_{ij}\) only in the requested order subdiagonals (working parameters
= \(\beta\) directly, unconstrained — a banded relative of the family-3
factorization idea). \(D\)’s remaining diagonal entries are \(q-1\) family-1
log-innovation-ratios (variance-type). A mixed structure like usm() but
banded/order-limited rather than fully dense, useful for
longitudinal/repeated-measures data where dependence should decay with
temporal order but need not follow a strict AR form.
maternm(x, range, nu, fixed, distance) — Matérn spatial correlation$$ K_{ij} = \frac{2^{1-\nu}}{\Gamma(\nu)} \left(\frac{\sqrt{2\nu},d_{ij}}{r}\right)^{!\nu} K_\nu!\left(\frac{\sqrt{2\nu},d_{ij}}{r}\right) $$
Family-1 (log) transform on both parameters: \(\eta = (\log r, \log \nu)\) (both must be positive: range and smoothness). Distance \(d_{ij}\) comes from supplied coordinates or a user-supplied distance matrix. Two correlation-type (in the “shape” sense — they reshape the whole matrix, not just rescale one entry) parameters governing decay rate and smoothness of spatial dependence.
sar(x, W, rho, fixed) — simultaneous autoregressive spatial covariance$$ B = I - \rho W, \qquad M = B^{-1}B^{-T}, \qquad K = M/M_{11} $$
Family-2 (bounded logit) transform: \(\rho\) bounded by the inverse spectral radius of \(W\) (guarantees \(B\) nonsingular). One correlation-type (shape) parameter that reshapes the entire spatial covariance through matrix inversion.
car(x, W, rho, fixed) — proper conditional autoregressive covariance$$ Q = D-\rho W,\quad D=\mathrm{diag}(\text{rowSums}(W)),\qquad M=Q^{-1},\qquad K=M/M_{11} $$
Family-2 (bounded logit) transform: \(\rho\) bounded by the reciprocal extreme
eigenvalues of the normalized adjacency \(S=D^{-½}WD^{-½}\) (the exact
condition for \(Q\) to stay valid). Same flavor as sar(): one
correlation/shape-type parameter, whole-matrix reshaping via inversion.
ownm(x, K, fun, par, fixed, dfun) — user-defined shapeEither a fixed known valid matrix K (zero free parameters, sommer only
normalizes it by its \([1,1]\) entry), or a user function fun(par) returning
a valid matrix, normalized the same way. Here the user is responsible for
choosing their own transform family from Chapter 2 (or inventing a new one);
sommer only imposes the final \(K_{11}=1\) normalization. If dfun(par,k)
(analytic \(B_k=\partial K/\partial\eta_k\)) is not supplied, ai_mme_sp2()
falls back to a numerical derivative. This is the escape hatch for any
structure not already covered — it plugs into the exact same
CovarianceFactor interface as maternm()/toeplitzm()/sar()/car().
covm(ran1, ran2, thetaC, theta, ...) — legacy two-random-effect helperPredates the general vsm() Kronecker interface; builds a
CovarianceFactor directly from two random-effect terms and an explicit
correlation-structure matrix thetaC/theta. Internally normalized and
compiled through the same descriptor contract as every structure above, so
it participates in REML identically once compiled — retained mainly for
backward compatibility with older model specifications.
| Constructor | Free working parameters \(\eta\) | Transform family (Ch. 2) | Type |
|---|---|---|---|
ism |
0 | — | — |
dsm/atm |
\(q-1\) | 1 (log-ratio) | variance |
csm (homogeneous) |
1 (rho) |
2 (bounded logit) | correlation |
csm (heterogeneous) |
\(1+(q-1)\) | 2 + 1 | mixed |
ar1m |
1 (rho) |
2 (tanh) | correlation |
ar2m/ar3m |
2 / 3 (PACF), +$(q-1)$ if heterogeneous | 2 (+1) | correlation / mixed |
toeplitzm |
\(q-1\) (PACF) | 2 (tanh) | correlation |
mam/ma1m/ma2m |
order (1 or 2) | none (identity) | correlation |
usm |
\(q(q+1)/2-1\) | 3 (Cholesky) + 1 (diag) | mixed (dense) |
corgm |
\(q(q-1)/2\) | 3 (unit-diag factor) | correlation (dense) |
fam |
loadings (\(k\)-rank) + \((q-1)\) specific | 3 (loadings) + 1 | mixed |
rrcm |
loadings (\(k\)-rank) | 3 (loadings) | correlation |
antem |
banded regression coeffs + \((q-1)\) innovations | 3 (banded) + 1 | mixed |
maternm |
2 (range, nu) | 1 (log) x2 | correlation (shape) |
sar/car |
1 (rho) |
2 (bounded logit) | correlation (shape) |
ownm |
user-defined | user-defined | user-defined |
\newpage
ai_mme_sp2()For the Henderson system of Chapter 1, the (restricted) log-likelihood
ai_mme_sp2() evaluates every iteration is
$$ \ell(\eta) = -\tfrac12\Big( \log|C| ;+; \log|R| ;+; \sum_i n_i\log|A_i| ;+; y’Py \Big) $$
mapped directly onto the code:
logDetC — \(\log|C|\), from the sparse Cholesky/LDLT factorization of
Henderson’s \(C\) matrix (or CHOLMOD/PCG-approximated for large sparse \(C\)).logDetR — \(\log|R|\), computed per residual-structure branch (diagonal
fast path, Kronecker-residual branch, or dense fallback).logDetA(i) — \(-\log|A_i^{-1}|\) per random effect (relationship/precision
matrices), via sparse LDLT.yPy — the quadratic form \(y'Py\) from absorbing \(y\) into the MME solve.REML maximizes \(\ell(\eta)\) over the working parameters \(\eta\) from Chapter 2 (not the natural parameters \(\theta\)), using a Newton-type Average-Information (AI) REML algorithm.
$$ \eta^{(t+1)} = \eta^{(t)} + \mathrm{InfMat}^{-1}, s(\eta^{(t)}) $$
InfMat: a convex blend of the Average Information matrix (avInf)
and the Expected Information matrix (emInf):InfMat = (weightAiInfMat * avInf) + (weightEmInfMat * emInf);
with per-iteration weights (emweight/weightEmInf, user/schedule
controlled) trading off AI-REML’s fast quadratic convergence against
classical EM-REML’s guaranteed monotone ascent and better global behavior
when the AI matrix is momentarily ill-conditioned. This blending is exactly
the mechanism that absorbs the numerical differences between variance-type
and correlation-type parameters described below — when the AI block for a
dense correlation structure is poorly conditioned, increasing the EM weight
locally stabilizes the step without touching the well-behaved variance
blocks.
For a covariance shape \(K(\eta)\), sommer works with its inverse (precision)
\(\Lambda = K^{-1}\) (lambdaDense in ai_mme_sp2()). All REML derivatives
route through the standard matrix-inverse identity:
$$ \frac{\partial \Lambda}{\partial \eta_k} = -\Lambda B_k \Lambda, \qquad B_k = \frac{\partial K}{\partial \eta_k} $$
implemented essentially verbatim:
arma::mat dLambda = -lambdaDense * Bk * lambdaDense; // Bk = dK/d(eta_k)
dLambda = 0.5 * (dLambda + dLambda.t()); // symmetrize (roundoff)
Bk (called cachedCovarianceD1()/covarianceD1() in the code) is where
every structure-specific formula from Chapter 3 actually enters the
optimizer — for analytic structures (diag, AR, CSM, US, …) it’s a
closed-form matrix; for generic-interface structures (maternm,
toeplitzm, sar, car, ownm without dfun) it’s a numerical
derivative of the R-level evaluator function. Either way, from this point
on ai_mme_sp2() treats Bk identically regardless of structure — the
only thing that changes the downstream numerics is what Bk looks like,
which is precisely the variance-vs-correlation distinction below.
For a pure variance-ratio working parameter (e.g. dsm()’s
\(\eta_j=\log(\theta_j/\theta_1)\), or vsm()’s own outer sigma2):
$$ K = \mathrm{diag}(1,\theta_2/\theta_1,\dots),\qquad B_j = \frac{\partial K}{\partial \eta_j} = \theta_j/\theta_1 \cdot e_j e_j’ $$
\(B_k\) is rank-one and diagonal — a single nonzero entry. Consequences:
avInf/emInf
(generically \(\mathrm{tr}(P\,\mathrm{d}C_k)\) and
\(\mathrm{tr}(P\,\mathrm{d}C_i\,P\,\mathrm{d}C_j)\)) reduce to a small
number of scalar row/column contractions instead of full dense matrix
products — this is exactly the “diagonal residual fast path”
specialization used for pure-variance residual structures, which only
exists because the derivative itself is diagonal, not because the
covariance is diagonal per se.sigma2 in vsm(), $\partial G/\partial\sigma^2 =
K$ exactly (the whole shape, undifferentiated) — the classical textbook
variance-component REML derivative. This is the parameter for which
Average Information and Expected Information coincide most closely in
well-behaved models, which is part of why classical variance-components
REML (pre-AI, EM-only) worked well for decades on models built purely
from ism()/dsm()-type structures.For a correlation/shape working parameter (e.g. csm()’s
\(\eta_\rho=\mathrm{logit}(\cdot)\), ar1m()’s \(\eta=\mathrm{atanh}(\rho)\),
usm()/corgm()’s Cholesky-factor entries):
$$ \frac{\partial K}{\partial \eta_\rho} = \frac{\partial K}{\partial \rho}\cdot\frac{\partial \rho}{\partial \eta_\rho} $$
Both factors are dense in general: changing a single correlation
parameter perturbs every off-diagonal entry of \(K\) simultaneously (e.g.
for ar1m(), \(K_{ij}=\rho^{|i-j|}\) so $\partial K_{ij}/\partial\rho =
|i-j|,\rho^{|i-j|-1}$ is nonzero for every \(i\neq j\)). Consequences:
avInf block for correlation parameters is fully dense and
more prone to near-singularity (e.g. when a correlation parameter is
weakly identified, or two correlation parameters trade off against each
other, as in usm() with many levels and modest sample size) — this is
precisely the situation the EM/AI blending in Section 4.2 exists to
stabilize: emInf (expected information) is guaranteed valid and
typically better conditioned than avInf in these dense/boundary-adjacent
regimes, so increasing its weight recovers a well-defined,
descent-guaranteeing step even when avInf alone would not.trust_cap: a second REML-stability mechanism tied to parameter typeIndependent of AI/EM blending, ai_mme_sp2() also enforces a per-parameter
trust-region cap (trust_cap in each CovarianceFactor) that bounds
how far a single working parameter may move in one iteration:
arma::vec caps = Rcpp::as<arma::vec>(f["trust_cap"]);
sigma2-like parameters can safely take large
multiplicative steps early in REML with little risk of overshoot causing
invalidity — only the likelihood itself, not validity, is at stake.maternm()‘s range/nu use
c(1.0, 0.75), sar()/car() use 1.0) are given tighter caps,
because a large step on the working scale, pushed through a saturating
transform (logit/tanh), can still correspond to a huge change in the
reconstructed matrix’s conditioning even though the value itself stays
valid — the cap prevents a technically-valid but wildly-overshooting
proposal from destabilizing the next likelihood evaluation.This is a purely numerical safeguard layered on top of the always-valid-by-construction guarantee from Chapter 2 — validity is never at risk, but good conditioning and monotone likelihood ascent are, and they are protected differently for the two parameter families.
flowchart TD
A["Working parameter eta (unconstrained, R)"] -->|transform g| B["Natural parameter theta (variance / correlation)"]
B -->|reconstruction formula| C["Covariance shape K(eta)"]
C -->|invert| D["Precision Lambda = K^-1"]
D -->|analytic identity dLambda = -Lambda Bk Lambda| E["dLambda/d(eta_k)"]
E --> F["Score s(eta), avInf, emInf"]
F -->|blend + trust_cap| G["AI-REML Newton step"]
G --> A
| Stage | Variance-type parameter | Correlation-type parameter |
|---|---|---|
| Transform (Ch. 2) | family 1: exponential (unbounded above) | family ⅔: logit/tanh/factorization (saturating, bounded) |
| \(B_k = \partial K/\partial \eta_k\) | sparse/diagonal, rank-1 | dense, full matrix |
| Cost of \(\mathrm{d}\Lambda_k\) | cheap (structured update) | \(O(q^3)\) dense product |
| avInf conditioning | usually well-behaved | can be near-singular, especially near boundary or weak identifiability |
| EM/AI blend role | rarely needed | often needed for stability |
trust_cap |
looser | tighter |
This is the concrete, code-level answer to “how does the REML process differ between a variance parameter and a correlation parameter”: both are optimized by the same AI-REML machinery over unconstrained working parameters (Chapter 2), but the structure of their derivative (\(B_k\), hence \(\mathrm{d}\Lambda_k\), hence the AI/EM information blocks) is fundamentally sparser and better-conditioned for variance parameters than for correlation parameters, which is why sommer layers EM/AI blending and per-parameter trust caps on top of a uniform derivative framework rather than hard-coding structure-specific optimizers.
\newpage
Every example below reuses the exact notation from Chapter 2 (\(\theta\), \(\eta\), \(g\), \(K\), \(\sigma^2\), \(\Lambda\), \(B_k\)) so you can trace each structure’s entry in Chapter 3 all the way to a concrete number.
ism()\(q=2\) levels, no working parameters at all: $K=I_2=\begin{pmatrix}1&0\0&1
\end{pmatrix}$, \(\Lambda = K^{-1} = I_2\). There is no \(\eta\) to iterate on
for this factor — it contributes only through the Kronecker product with
the other factors in the same vsm() term.
dsm()\(q=3\) levels, natural variances \(\theta=(1, 2.25, 0.81)\) (already showing \(\theta_1=1\), the reference from Chapter 2.6). Family-1 transform:
$$ \eta_2=\log(2.25/1)=0.8109,\qquad \eta_3=\log(0.81/1)=-0.2107 $$
$$ K=\mathrm{diag}(1,\ 2.25,\ 0.81),\qquad \Lambda = K^{-1} = \mathrm{diag}(1,\ 0.4444,\ 1.2346) $$
Structural derivative for \(\eta_2\): $B_2 = \partial K/\partial \eta_2 = 2.25\cdot e_2e_2’ = \mathrm{diag}(0,2.25,0)$ — a single nonzero entry, the sparse/diagonal pattern of Chapter 4.4.
csm(), homogeneous\(q=3\), \(lo=-1/(q-1)=-0.5\). Suppose \(\rho=0.4\). Family-2 transform:
$$ p=\frac{\rho-lo}{1-lo}=\frac{0.4+0.5}{1.5}=0.6,\qquad \eta_\rho=\mathrm{logit}(0.6)=\log\frac{0.6}{0.4}=0.4055 $$
$$ K=\begin{pmatrix}1&0.4&0.4\0.4&1&0.4\0.4&0.4&1\end{pmatrix} $$
\(B_1=\partial K/\partial\eta_\rho\) is dense (every off-diagonal entry moves together) — matching Chapter 4.5’s correlation-type pattern.
csm(), heterogeneousSame \(\rho=0.4\) (\(\eta_\rho=0.4055\)), plus variances \(\theta=(1,\ 1.5,\ 0.6)\) giving \(\eta_2=\log(1.5)=0.4055\), \(\eta_3=\log(0.6)=-0.5108\) (family-1, exactly as in 5.2). Full working vector: \(\eta=(0.4055,\ 0.4055,\ -0.5108)\) — 1 correlation-type entry followed by 2 variance-type entries, the “mixed” row of the Chapter 3.15 table.
ar1m()\(q=4\), \(\rho=0.5\). Family-2 (tanh) transform: $\eta=\mathrm{atanh}(0.5)= 0.5493$.
$$ K=\begin{pmatrix} 1 & 0.5 & 0.25 & 0.125\ 0.5 & 1 & 0.5 & 0.25\ 0.25 & 0.5 & 1 & 0.5\ 0.125 & 0.25 & 0.5 & 1 \end{pmatrix} $$
\(K_{ij}=\rho^{|i-j|}\); the single working parameter reshapes every off-diagonal band simultaneously (dense \(B_1\), Chapter 4.5).
ar2m() / ar3m() / toeplitzm()\(q=4\), order-2 PACF \(\kappa=(0.3,-0.2)\). Family-2 transform:
\(\eta=(\mathrm{atanh}(0.3), \mathrm{atanh}(-0.2)) = (0.3095,\ -0.2027)\). The
Durbin–Levinson recursion (Chapter 3.4) turns \(\kappa\) into full lag
correlations \(\rho_1,\rho_2,\rho_3\), which then fill \(K_{ij}=\rho_{|i-j|}\)
exactly like ar1m() above but allowing a richer decay pattern.
toeplitzm() with \(q=4\) would use 3 PACF values (one per lag up to \(q-1\))
instead of being capped at a fixed AR order.
mam() / ma1m() / ma2m()\(q=3\), order 1, \(\theta_1=0.3\) (the raw MA coefficient — recall from Chapter 3.5 there is no transform, so \(\eta=\theta=0.3\) directly). The resulting \(K\) is built from \(\mathrm{Cov}(e_t+\theta_1 e_{t-1})\), giving a banded correlation matrix with only lag-1 correlation nonzero. Because \(g\) is the identity map here, \(B_1=\partial K/\partial\eta_1 = \partial K/\partial\theta_1\) directly — no chain rule factor to track, unlike every bounded/positive structure above.
usm()\(q=2\) (the case from the original single-structure example, reproduced here in the shared notation). Working parameters \(\eta=(\eta_{21},\ \eta_{22})=(0.30,\ \log(0.95))\), i.e. \(L_{21}=0.30\) (family-3, off-diagonal, no transform) and \(L_{22}=\exp(\eta_{22})=0.95\) (family-3 diagonal, family-1 log transform):
$$ L=\begin{pmatrix}1&0\0.30&0.95\end{pmatrix},\qquad K=LL’=\begin{pmatrix}1.0000&0.3000\0.3000&0.9925\end{pmatrix} $$
$$ \Lambda=K^{-1}=\frac{1}{|K|}\begin{pmatrix}K_{22}&-K_{12}\-K_{12}&K_{11} \end{pmatrix},\qquad |K|=0.9025 $$
$$ \Lambda\approx\begin{pmatrix}1.0997&-0.3325\-0.3325&1.1080\end{pmatrix} $$
Structural derivatives (\(B_1=\partial K/\partial\eta_{21}\), \(B_2=\partial K/\partial\eta_{22}\), chain rule \(\partial K/\partial\eta_{22} = \partial K/\partial L_{22}\cdot L_{22}\) since \(L_{22}=\exp(\eta_{22})\)):
$$ B_1=\begin{pmatrix}0&1\1&0.60\end{pmatrix},\qquad B_2=\begin{pmatrix}0&0\0&1.805\end{pmatrix} $$
\(B_1\) (off-diagonal, correlation-type) is dense; \(B_2\) (log-diagonal, variance-type in flavor) is sparse/rank-one — Chapter 4.4/4.5’s pattern, both appearing inside the same structure.
corgm()\(q=3\), unit-diagonal lower factor entries \(A_{21}=0.2\), \(A_{31}=0.1\), \(A_{32}=0.15\) (family-3, no per-entry transform — any real numbers here give a valid unit-diagonal \(K\)):
$$ A=\begin{pmatrix}1&0&0\0.2&\sqrt{1-0.2^2}&0\0.1&\ (\dots)&(\dots) \end{pmatrix} $$
(rows are renormalized to unit length internally so \(K_{ii}=1\) exactly);
working parameters are simply \(\eta=(0.2,\ 0.1,\ 0.15)\), i.e. par equals
the lower-triangular entries directly, with no chain-rule factor to apply
when differentiating \(K\) with respect to \(\eta\).
fam() / rrcm()\(q=3\), rank \(k=1\). Loadings \(\Lambda=(0.6,\ 0.4,\ 0.2)'\) (with the leading entry \(\Lambda_{11}=0.6>0\) log-parameterized: \(\eta_{load,1}=\log(0.6)=-0.5108\); the remaining loadings unconstrained: \(\eta_{load,2}=0.4\), \(\eta_{load,3}=0.2\)), specific variances \(\psi=(1,\ 0.5,\ 0.3)\) (family-1 log-ratio on levels 2-3: \(\eta_{spec,2}=\log(0.5)=-0.6931\), \(\eta_{spec,3}=\log(0.3)=-1.2040\)):
$$ M=\Lambda\Lambda’+\Psi= \begin{pmatrix}0.36+1&0.24&0.12\0.24&0.16+0.5&0.08\0.12&0.08&0.04+0.3 \end{pmatrix}, \qquad K=M/M_{11} $$
rrcm() would drop the \(\psi\) working parameters entirely (fixing
\(\Psi=I\)) and keep only the loading parameters \(\eta_{load}\).
antem()\(q=3\), order 1. Regression coefficients on the first subdiagonal only,
\(\beta_{21}=0.4\), \(\beta_{32}=0.3\) (unconstrained, \(\eta=\beta\) directly —
same “no per-entry transform” logic as corgm()’s off-diagonal factor
entries), innovations \(D=\mathrm{diag}(1,\ 0.7,\ 0.5)\) (family-1 log-ratio:
\(\eta_{D,2}=\log(0.7)=-0.3567\), \(\eta_{D,3}=\log(0.5)=-0.6931\)):
$$ T=\begin{pmatrix}1&0&0\-0.4&1&0\0&-0.3&1\end{pmatrix},\qquad K=T^{-1}DT^{-T}\ /\ (T^{-1}DT^{-T})_{11} $$
maternm()\(q=3\) spatial points, distances \(d_{12}=1,\ d_{13}=2,\ d_{23}=1\). Suppose \(r=1.5\), \(\nu=0.5\) (Matérn with \(\nu=0.5\) reduces to the exponential correlation \(K_{ij}=\exp(-d_{ij}/r)\)). Family-1 transform on both: \(\eta=(\log 1.5,\ \log 0.5)=(0.4055,\ -0.6931)\).
$$ K_{12}=e^{-1/1.5}=0.5134,\quad K_{13}=e^{-2/1.5}=0.2636,\quad K_{23}=e^{-1/1.5}=0.5134 $$
Both working parameters are correlation/shape-type: changing either \(r\) or \(\nu\) reshapes every pairwise entry of \(K\) at once (dense \(B_k\), as in Chapter 4.5).
sar()\(q=3\) areal units, adjacency \(W\) with spectral radius (largest eigenvalue magnitude) \(=2\), so the valid range is \(\rho\in(-0.5,0.5)\). Suppose \(\rho=0.2\). Family-2 (symmetric bounded logit) transform:
$$ p=\frac{\rho+0.5}{1.0}=0.7,\qquad \eta=\mathrm{logit}(0.7)=0.8473 $$
$$ K = M/M_{11},\qquad M=(I-0.2W)^{-1}(I-0.2W)^{-T} $$
car()\(q=3\) areal units, symmetric non-negative adjacency \(W\), row sums \(D=\mathrm{diag}(3,2,2)\). Suppose the valid range works out to \(\rho\in(-1.2,\ 1.0)\) and \(\rho=0.3\). Family-2 transform:
$$ p=\frac{0.3+1.2}{2.2}=0.6818,\qquad \eta=\mathrm{logit}(0.6818)=0.7677 $$
$$ K = M/M_{11},\qquad M = (D-0.3W)^{-1} $$
ownm()A user supplies fun <- function(par) matrix(c(1, par[1], par[1], 1), 2, 2)
with par=c(0.25) — here the user has chosen to reuse the family-2
(bounded) idea manually by ensuring par[1] always stays in \((-1,1)\) some
other way (sommer does not enforce this for ownm()); \(\eta=\theta=0.25\)
directly, sommer only performs the \(K_{11}=1\) normalization from Chapter
2.6 and, since no dfun was supplied, computes \(B_1\) by a numerical
(finite-difference) derivative of fun.
\newpage
Chapters 2-5 showed how a working parameter \(\eta\) turns into a shape \(K\)
and a precision \(\Lambda\). This chapter shows where in the computation
those numbers actually get used, by carrying one tiny, fully-numeric example
all the way from \(\eta\) to the REML log-likelihood \(\ell(\eta)\) itself —
reusing the exact usm() numbers from Section 5.8.
Two lines (\(L_1,L_2\)), two traits (\(A,B\)), one observation per
line-trait combination, model \(y=\mathbf{1}\beta + u + e\) (a single fixed
intercept \(\beta\), one random effect per line-trait cell, and residual
noise). Take the random-effect design as vsm(usm(Trait), ism(Line)),
which makes each line’s 2-vector of trait effects
\(u_{l}=(u_{l,A},u_{l,B})\) independent across lines with common covariance
\(\sigma^2 K\):
$$ y = \begin{pmatrix}10\12\8\9\end{pmatrix}\ \ (\text{order: } L_1A,\ L_1B,\ L_2A,\ L_2B), \qquad X=\begin{pmatrix}1\1\1\1\end{pmatrix},\qquad Z=I_4,\qquad R=I_4 $$
To isolate exactly what the shape working parameters do, fix the scale at
\(\sigma^2=1\) and the residual variance at \(1\) (in a real mmes() fit these
would be extra variance-type working parameters of their own, following
Section 4.4).
Exactly the Section 5.8 usm() example: $\eta=(\eta_{21},\eta_{22}) =
(0.30,\ \log 0.95)$, so \(L_{21}=0.30\), \(L_{22}=0.95\), and
$$ K=\begin{pmatrix}1.0000&0.3000\0.3000&0.9925\end{pmatrix}, \qquad \Lambda=K^{-1}=\begin{pmatrix}1.0997&-0.3324\-0.3324&1.1080\end{pmatrix} $$
Because lines are independent (ism(Line)), the random-effect precision
\(G^{-1}=\sigma^{-2}\,(I_2\otimes\Lambda)\) is block-diagonal, one copy of
\(\Lambda\) per line — this is literally the only place \(\eta\) enters the
whole computation; every quantity from here on is arithmetic built on top of
these four numbers. Henderson’s \(C\) (Section 1.2), with \(R=I_4\) so
\(Z'R^{-1}Z=I_4\):
$$ C=\begin{pmatrix} 4 & 1&1&1&1\ 1 & 2.0997&-0.3324&0&0\ 1 & -0.3324&2.1080&0&0\ 1 & 0&0&2.0997&-0.3324\ 1 & 0&0&-0.3324&2.1080 \end{pmatrix} $$
rows/columns ordered \((\beta,\,u_{L_1A},\,u_{L_1B},\,u_{L_2A},\,u_{L_2B})\). The \(\Lambda\) blocks sit exactly where \(Z'R^{-1}Z\) would otherwise be left alone — \(\eta\)’s only job, computationally, is to determine these four numbers added into \(C\).
Solving \(C\,(\hat\beta,\hat u)' = (X'R^{-1}y,\ Z'R^{-1}y)' = (39,\,10,\,12,\,8,\,9)'\) gives
$$ \hat\beta \approx 9.752,\quad \hat u_{L_1A}\approx 0.294,\ \hat u_{L_1B}\approx 1.113,\quad \hat u_{L_2A}\approx -0.914,\ \hat u_{L_2B}\approx -0.501 $$
If \(\eta\) had been different, \(\Lambda\) would differ, the four numbers inserted in Section 6.3 would differ, and every one of these five solutions would come out different — this is the concrete meaning of “the BLUPs/BLUEs depend on the covariance parameters.”
Using the identity \(y'Py = y'R^{-1}y - (X'R^{-1}y,\,Z'R^{-1}y)\cdot(\hat\beta,\hat u)\):
$$ y’R^{-1}y = 10^2+12^2+8^2+9^2 = 389, \qquad y’Py \approx 389 - 384.81 \approx 4.19 $$
and, from the same \(C\) matrix of Section 6.3 (a block-determinant calculation, done in practice by sparse Cholesky, not shown here step by step),
$$ \log|C| \approx 3.48 $$
With ism(Line) contributing \(\log|A^{-1}|=0\) and \(R=I_4\) contributing
\(\log|R|=0\), Chapter 4.1’s REML log-likelihood becomes a single concrete
number:
$$ \ell(\eta) = -\tfrac12\big(\log|C| + \log|R| + \textstyle\sum_i n_i\log|A_i| + y’Py\big) \approx -\tfrac12(3.48+0+0+4.19) \approx -3.83 $$
This is the number AI-REML is trying to make as large as possible, and every term in it was produced starting from the two numbers in \(\eta\).
Section 4.3 gives the analytic derivative \(\partial\Lambda/\partial\eta_k = -\Lambda B_k\Lambda\). We can check this numerically instead of trusting the formula blindly, using the same finite- difference idea a high-school calculus class would use for any derivative: nudge \(\eta_{21}\) from \(0.30\) to \(0.31\) (holding \(\eta_{22}\) fixed), rebuild \(K\) and \(\Lambda\) exactly as in Section 6.2, and see how much \(\Lambda\) moved.
Nudged value: \(L_{21}=0.31,\ L_{22}=0.95\), giving \(K_{22}=0.31^2+0.95^2=0.9986\) (note \(\det K = (L_{11}L_{22})^2=0.9025\) is unchanged, since only the off-diagonal \(\eta_{21}\) moved) and
$$ \Lambda_{\text{nudged}} \approx \begin{pmatrix}1.1064&-0.3435\-0.3435&1.1080\end{pmatrix} $$
Finite-difference estimate of the derivative:
$$ \frac{\Lambda_{\text{nudged}}-\Lambda}{0.01} \approx \begin{pmatrix}0.668&-1.108\-1.108&0.000\end{pmatrix} $$
Compare this to the analytic formula from Section 4.⅗.8, \(B_1=\partial K/\partial\eta_{21}=\begin{pmatrix}0&1\\1&0.60\end{pmatrix}\), plugged into \(-\Lambda B_1 \Lambda\):
$$ -\Lambda B_1 \Lambda \approx \begin{pmatrix}0.665&-1.108\-1.108&0.000\end{pmatrix} $$
The two match (up to the expected finite-step-size rounding). This is
exactly the object ai_mme_sp2() forms internally as dLambda (Section
4.3), and it is what turns into the score \(s(\eta)\) and the avInf/emInf
blocks of Section 4.2 that decide the next value of \(\eta\) — closing the
loop: \(\eta \to K \to \Lambda \to C \to (\hat\beta,\hat u) \to \ell(\eta)\)
for the forward pass (Sections 6.2-6.5), and
$\eta \to B_k \to \mathrm{d}\Lambda_k \to s(\eta),,\mathrm{InfMat}
\to \eta^{(t+1)}$ for the update that drives REML to convergence (Chapter
4.2).
\newpage
?vsm and the
constructor’s own help page (e.g. ?csm, ?ar1m, ?usm).CovarianceFactor,
free, par_names, trust_cap), see using-covariance-structures.Rmd.ai_mme_sp2() referenced in Chapter 4 (sparse LDLT, CHOLMOD, PCG,
selected-inverse traces), see the package’s C++ source comments in
src/MNR.cpp.