stanAMmv: Animal model (multivariate)

In timflutre/rutilstimflutre: Timothee Flutre's personal R code

Animal model (multivariate)

Description

Given T traits, I genotypes, Q covariates and N=I*Q phenotypes per trait, fit an "animal model" with the rstan package via the following likelihood: Y = W C + Z G_A + Z G_D + E, where Y is NxT; W is NxQ; Z is NxI; G_A ~ Normal_IxT(0, A, V_G_A) with A the known matrix of additive genetic relationships; G_D ~ Normal_IxT(0, D, V_G_D) with D the known matrix of dominant genetic relationships; E ~ Normal_NxT(0, Id_N, V_E); Missing phenotypes are jointly imputed with the other unknown variables, and the errors can follow a Student's t distribution to handle outliers.

Usage

stanAMmv(
  data,
  relmat,
  errors.Student = FALSE,
  nb.chains = 1,
  nb.iters = 10^3,
  burnin = 10^2,
  thin = 10,
  out.dir = getwd(),
  task.id = format(Sys.time(), "%Y-%m-%d-%H-%M-%S"),
  compile.only = FALSE,
  rm.stan.file = FALSE,
  rm.sm.file = FALSE,
  verbose = 0
)

Arguments

data	data.frame containing the data corresponding to relmat; should have columns grep-able for "response" as well as a column "geno.add" used with matrix A; if a column "geno.dom" exists, it will be used with matrix D; any other column will be interpreted as corresponding to "fixed effects"
relmat	list containing the matrices of genetic relationships (see estimGenRel); additive relationships (matrix A) are compulsory, with name "geno.add"; dominant relationships (matrix D) are optional, with name "geno.dom"; can be in the "matrix" class (base) or the "dsCMatrix" class (Matrix package)
errors.Student	use a Student's t distribution for the errors (useful to handle outliers)
nb.chains	number of independent chains to run
nb.iters	number of iterations to monitor
burnin	number of initial iterations to discard
thin	thinning interval for monitored iterations
out.dir	working directory in which the stan and sm files will be written
task.id	identifier of the task (used in temporary file names)
compile.only	only compile the model, don't run it
rm.stan.file	remove the file specifying the model written in the STAN language
rm.sm.file	remove the file corresponding to the compiled model
verbose	verbosity level (0/1)

Value

path to compiled file if compile.only=TRUE, object of class stanfit-class otherwise

Author(s)

Timothee Flutre

See Also

stanAM

Examples

## Not run: ## simulate genotypes
set.seed(1859)
X <- simulGenosDose(nb.genos=200, nb.snps=2000)

## simulate phenotypes with only additive part of genotypic values
A <- estimGenRel(X, relationships="additive", method="vanraden1", verbose=0)
modelA <- simulAnimalModel(T=2, Q=3, A=A, nu.G.A=5, nu.E=3, seed=1859)

## infer with rstan
modelA$data$geno.add <- modelA$data$geno; modelA$data$geno <- NULL
fitA <- stanAMmv(data=modelA$data, relmat=list(geno.add=A))
fitA <- rstan::As.mcmc.list(fitA)
plotMcmcChain(fitA[[1]], "sigma_g_A", 1:4, sqrt(modelA$V.G.A))
cbind(truth=c(c(modelA$C), sqrt(modelA$V.G.A), sqrt(modelA$V.E)),
      summaryMcmcChain(fitA[[1]], c("c[1]", "c[2]", "c[3]", "sigma_g_A", "sigma_E")))

## simulate phenotypes with additive and dominant parts of genotypic values
D <- estimGenRel(X, relationships="dominant", method="vitezica", verbose=0)
modelAD <- simulAnimalModel(T=2, Q=3, A=A, nu.G.A=15, nu.E=5,
                            D=D, nu.G.D=3, seed=1859)

## infer with rstan
modelAD$data$geno.add <- modelAD$data$geno; modelAD$data$geno <- NULL
fitAD <- stanAMmv(data=modelAD$data, relmat=list(geno.add=A, geno.dom=D))
fitAD <- rstan::As.mcmc.list(fitAD)
plotMcmcChain(fitAD[[1]], "sigma_g_A", 1:4, sqrt(modelAD$V.G.A))
plotMcmcChain(fitAD[[1]], "sigma_g_D", 1:4, sqrt(modelAD$V.G.D))
cbind(truth=c(c(modelAD$C), sqrt(modelAD$V.G.A), sqrt(modelAD$V.G.D), sqrt(modelAD$V.E)),
      summaryMcmcChain(fitAD[[1]], c("c[1]", "c[2]", "c[3]", "sigma_g_A", "sigma_g_D", "sigma_E")))
plot(as.vector(fitAD[[1]][,"sigma_g_A"]), as.vector(fitAD[[1]][,"sigma_g_D"]))

## End(Not run)

timflutre/rutilstimflutre documentation built on Feb. 12, 2025, 11:35 p.m.