Covariance and correlation matrix of a fitted term with standard errors

Returns the full covariance matrix of a random or residual term of a fitted mmes model over its covariance coordinates (for example the genetic trait-by-trait covariance matrix of a multi-trait model), together with the implied correlation matrix and delta-method standard errors of every element.

Usage

covmatrix_mmes(object, term, se = TRUE, max.dim = 500L, rel_step = 1e-5)

Arguments

object: A fitted mmes model.
term: Name (as in names(object$covStruct)) or position of the term. The residual term is the last one.
se: If TRUE, add standard errors.
max.dim: Refuse terms whose covariance coordinate space is larger than this (the matrices are dense).
rel_step: Relative step of the central differences used for the Jacobian.

Details

The matrix is $\sigma^2 K_1 \otimes \cdots \otimes K_m$ built from all the covariance factors of the vsm term (the relationship matrix Gu is not included), so for vsm(usm(trait), ism(id)) it is the trait covariance matrix, and for vsm(dsm(env), usm(trait), ism(id)) the environment-by-trait matrix. Row and column names are the coordinate levels.

Standard errors use the delta method, $\mathrm{Var}(g(\hat\theta)) \approx J V J^\top$, where $V$ is object$theta_se (the covariance of the reported parameters, as used by covparams_mmes_se) and $J$ the Jacobian of each element with respect to the reported parameters. Fixed parameters contribute no variance. Diagonal standard errors of the correlation matrix are zero. Terms split by dsumm sections are not supported.

Value

A list with covariance and correlation and, if se=TRUE, covariance.se and correlation.se.

Examples

set.seed(17)
n <- 60
DT <- data.frame(id = factor(rep(seq_len(n), each = 3)))
u <- matrix(rnorm(n * 2), n) 
DT$trait_a <- u[DT$id, 1] + rnorm(nrow(DT), sd = 0.7)
DT$trait_b <- u[DT$id, 2] + rnorm(nrow(DT), sd = 1)
DTL <- stackTraits(DT, traits = c("trait_a", "trait_b"), keep = "id")
fit <- mmes(value ~ trait,
            random = ~ vsm(usm(trait), ism(id)),
            rcov = ~ vsm(usm(trait), ism(record)),
            data = DTL, verbose = FALSE)
G <- covmatrix_mmes(fit, 1)       # genetic
G$covariance
G$correlation
G$correlation.se
R <- covmatrix_mmes(fit, 2)       # residual
R$correlation