jagsAM | R Documentation |
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.
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
)
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 |
inits |
list of initial values (possible to use 1 sub-list per chain, see |
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 |
nb.adapt |
number of iterations for adaptation (see |
burnin |
number of initial iterations to discard (see the update function of the rjags package) |
nb.iters |
number of iterations to monitor (see |
thin |
thinning interval for monitored iterations (see |
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) |
mcmc.list
object
Timothee Flutre
lmerAM
, inlaAM
, stanAM
## 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.