Covariance structure spanning several random terms

Combines two or more random terms built with vsm into one random structure with an estimated covariance between the terms. Typical uses are direct-maternal animal models, correlated permanent-environment effects and indirect genetic effects.

Usage

strm(..., cov = usm, Gu = NULL, labels = NULL, sigma2 = NULL,
     fixedSigma2 = FALSE)

Arguments

...: Two or more random vsm() terms, optionally named (names become the term labels), e.g. dir = vsm(ism(id)), mat = vsm(ism(dam)).
cov: A covariance constructor applied to the term index, e.g. usm (default), dsm, csm, corgm, or function(x) fam(x, k = 1).
Gu: Precision matrix shared by all terms (attr(Gu, "inverse") = TRUE). If NULL, a Gu given inside the terms is used (all terms must then use the same one) or an identity over the union of the terms' levels.
labels: Term labels; default the argument names or t1, t2, ....
sigma2: Starting overall variance scale; NULL lets mmes choose.
fixedSigma2: Fix the overall variance scale.

Details

For terms $u_1,\dots,u_T$ with incidence matrices $Z_1,\dots,Z_T$ the structure is $$\mathrm{Var}(u_1,\dots,u_T)=\sigma^2\,K_{terms}\otimes K_{inner}\otimes A,$$ where $K_{terms}$ is given by cov, $K_{inner}$ is the Kronecker product of the covariance factors used inside the terms (e.g. dsm(env)), and $A$ is the relationship matrix. The requirements are those of a direct product: all terms share the coefficient levels (their levels are aligned on the union, or on the levels of Gu), the same Gu, and the same inner covariance constructors with the same levels; the parameters of the inner factors are shared. The sigma2 of individual terms is ignored. BLUPs are returned in one uList element with a column per term (and inner level). covm(ran1, ran2) is the two-term special case with cov = usm.

The terms do not need to be adjacent in the model formula and their order is checked by level names rather than assumed.

Value

A random structure for the random formula of mmes.

Examples

set.seed(3)
ids <- paste0("i", 1:30)
A <- diag(30) + diag(0.2, 30); dimnames(A) <- list(ids, ids)
Ai <- solve(A); attr(Ai, "inverse") <- TRUE
d <- data.frame(id = factor(sample(ids, 200, TRUE), levels = ids),
                dam = factor(sample(ids[1:15], 200, TRUE), levels = ids))
d$y <- rnorm(200)
fit <- mmes(y ~ 1, random = ~ strm(dir = vsm(ism(id)), mat = vsm(ism(dam)), Gu = Ai),
            data = d, verbose = FALSE)
covparams_mmes(fit, 1)