jagsAM: Animal model

View source: R/quantgen.R

jagsAMR Documentation

Animal model

Description

Given I genotypes, Q covariates and N=I*Q phenotypes for the trait, fit an "animal model" with the rjags package via the following likelihood: y = W c + Z g_A + Z g_D + epsilon, where y is Nx1; W is NxQ; Z is NxI; g_A ~ Normal_I(0, sigma_A^2 A) with A the known matrix of additive genetic relationships; g_D ~ Normal_I(0, sigma_D^2 D) with D the known matrix of dominant genetic relationships; epsilon ~ Normal_N(0, sigma^2 Id_N); Cov(g_A,g_D)=0; Cov(g_A,e)=0; Cov(g_D,e)=0.

Usage

jagsAM(
  data,
  relmat,
  inits = NULL,
  priors = list(fix = list(dist = "c", par = 5), vc = list(dist = "hc", par = 5)),
  nb.chains = 1,
  nb.adapt = 10^3,
  burnin = 10^2,
  nb.iters = 10^3,
  thin = 10,
  progress.bar = NULL,
  rm.jags.file = TRUE,
  verbose = 0
)

Arguments

data

data.frame containing the data corresponding to relmat; should have a column 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)

inits

list of initial values (possible to use 1 sub-list per chain, see jags.model); if NULL, JAGS will choose typical values (usually mean, median, or mode of the prior)

priors

list of specifying priors; each component is itself a list for which "dist" specifies the distribution and "par" the parameter; a component named "fix" is used for fixed effects, so that ("c",5) means c[q] ~ Cauchy(0,5) and ("n",10^6) means c[q] ~ Normal(0,10^6), but note that the global intercept always is Normal(mean(y),10^6); a component named "vc" is used for variance components, so that ("hc",5) means stdev ~ half-Cauchy(0,5) and ("ig",0.001) means var ~ InvGamma(shape=0.001,rate=0.001)

nb.chains

number of independent chains to run (see jags.model)

nb.adapt

number of iterations for adaptation (see jags.model)

burnin

number of initial iterations to discard (see the update function of the rjags package)

nb.iters

number of iterations to monitor (see coda.samples)

thin

thinning interval for monitored iterations (see coda.samples)

progress.bar

type of progress bar (text/gui/none or NULL)

rm.jags.file

remove the file specifying the model written in the JAGS-dialect of the BUGS language

verbose

verbosity level (0/1)

Value

mcmc.list object

Author(s)

Timothee Flutre

See Also

lmerAM, inlaAM, 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=1, Q=3, A=A, V.G.A=15, V.E=5, seed=1859)

## infer with rjags
modelA$data$geno.add <- modelA$data$geno; modelA$data$geno <- NULL
fitA <- jagsAM(data=modelA$data, relmat=list(geno.add=A))
plotMcmcChain(fitA[[1]], "V.G.A", 1:4, modelA$V.G.A)
cbind(truth=c(c(modelA$C), modelA$V.G.A, modelA$V.E),
      summaryMcmcChain(fitA[[1]], c("c[1]", "c[2]", "c[3]", "sigma.A2", "V.E")))

## End(Not run)

timflutre/rutilstimflutre documentation built on Feb. 7, 2024, 8:17 a.m.