MTM: Fits a (Bayesian) Multivariate Gaussian Mixed Effects Model...
In QuantGen/MTM: MTM

Description Usage Arguments Details Value Examples

View source: R/MTM.R

Data equation: Y = 1mu' + XB + U1 + ... + Uq + E

where:

Y (numeric, nxp) is a matrix of phenotypes (individuals in rows, traits in columns, NAs accepted),
mu is a vector of (p) intercepts (included by default),
X (nxq) is an incidence matrix for q fixed effects,
B is a matrix of fixed effects,
U1, ..., Uq are (nxq) matrices of random effects with vec(Uj)~N(0, kronecker(Gj, Kj)), where Gj is a pxp (unknown) co-variance matrix, Kj is an nxn user-defined covariance matrix (kernel) which must by numeric, symmetric positive semi-definite,
E (nxp) is a matrix of model residuals, assumed to follow a MVN distribution vec(E)~N(0, kronecker(R0, I)).

1
2
3

MTM(Y, XF = NULL, K = NULL, resCov = list(type = "UN", df0 = 0, S0 =
  diag(0, ncol(as.matrix(Y)))), nIter = 110, burnIn = 10, thin = 2,
  saveAt = "", tolD = 1e-05)

`Y`	Phenotype matrix (nxp numeric, traits (p) in columns, individuals (n) in rows).
`XF`	A numeric design matrix (nxq) for fixed effects. For factors use XF=as.matrix(model.matrix(~x+y...))[, -1].
`K`	A 2-level list, 1st level defines random effects, inside each level a list is used to provide the kernel (K), a covariance structure (type, 'UN', 'DIAG', 'FA' supported) and hyper-parameters (degree of freedom, df0, and scale, S0).
`resCov`	A list used to define the co-variance matrix for model residuals (R0). Example: resCov=list(type='UN', df0=x, S0=V) specifies an un-structured covariance matrix, with an Inverse Whishart prior with degree of freedom df0 (scalar) and scale matrix (pxp) V.
`nIter`	The number of iterations (integer).
`burnIn`	The number of iterations to be discarded as burn-in (integer).
`thin`	Thinin interval (integer).
`saveAt`	A character path and a prefix used to define where to store samples (e.g., saveAt='c:/mtmFit/test_'. By default samples are saved in the current directory and filenames have no prefix.
`tolD`	A numeric parameter used to define the minimum eigenvalue to be maintained in the model. Eigenvectors of kernels smaller than tolD are removed. The default value is tolD=1e-6.

Intercepts are included by default, fixed and random effects are optional. The same fixed effects are applied to all traits. For each random effect the user must provide a kernel (K). By default the residual co-variance matrix (R0) and the co-variance matrices of random effects are un-structured; however users can specify other models (e.g., DIAG=diagonal, FA=factor analysis, and REC=recursive).

List containing estimated posterior means and estimated posterior standard deviations, including: $yHat.

## Not run: 
### Example: ##################################################################
# Model with an un-structured covariance matrix for the random effects of     #
# line and a diagonal residual co-variance matrix.                            #
###############################################################################
rm(list = ls())
library(BGLR)  # only needed to access wheat data. MTM does not depend on BGLR.
data(wheat)
X <- scale(wheat.X, center = TRUE, scale = TRUE)
Y <- wheat.Y
G <- tcrossprod(X)/ncol(X)
nTraits <- ncol(wheat.Y)
df0 <- 3
R2 <- 0.2
S0 <- var(Y)

# Defining the specification of the residual covariance matrix (diagonal).
R0 <- list(type = "DIAG", df0 = rep(df0, nTraits),
           S0 = diag(S0) * (1 - R2) * (df0 + 2))
R0 # Note, when type='DIAG', independent scaled-inv. chi-sq. prior are assigned
# to the diagonal elements, hence, df0 and S0 are vectors, each entry indexes
# one prior.

# Defining the specification of the covariance matrix of random effects
# (un-structured)
KList <- list(list(K = G, COV = list(type = "UN", S0 = diag((1 - R2) * (df0 + 2),
              nTraits), df0 = df0)))
KList[[1]]$COV
# Note for un-structured covariance matrix the prior is Inverse Wishart therefore,
# df0 is a scalar and S0 a pd matrix.

# Fitting the model
fm1 <- MTM(Y = Y, K = KList, resCov = R0, nIter = 1200, burnIn = 200, saveAt = "ex1_")

# Some estimates
fm1$K[[1]]$G  # Estimated genomic covariance matrix
fm1$mu  # Estimated intercepts
fm1$K[[1]]$U  # Predicted genomic values
fm1$resCov$R  # Estimated residual covariance matrix

# Samples and trace plots.  Genomic covariance matrix
G <- read.table("ex1_G_1.dat")
head(G)
plot(G[, 1], type = "o", cex = 0.5, col = 4)  # genomic variance 1st trait.
plot(G[, 2], type = "o", cex = 0.5, col = 4)  # genomic co-variance trait 1 and 2.
# ...
plot(G[, 5], type = "o", cex = 0.5, col = 4)  # genomic variance trait 2.
plot(G[, 6], type = "o", cex = 0.5, col = 4)  # genomic co-variance trait 2 and 3.
library(MCMCpack)
GHat <- xpnd(colMeans(G))

## Residual covariance matrix
R <- read.table("ex1_R.dat")
head(R) # Note, since type='DIAG' several entries (the covariances) are all equal
# to zero.

# Computing Samples for genomic heritabilities (example with trait 1)
H2_1 <- G[, 1]/(G[, 1] + R[, 1])
plot(density(H2_1))

## End(Not run)