stanAM: Animal model (univariate)

View source: R/quantgen.R

stanAMR Documentation

Animal model (univariate)

Description

Given I genotypes, Q covariates and N=I*Q phenotypes for the trait, fit an "animal model" with the rstan 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. Missing phenotypes are jointly imputed with the other unknown variables, and the errors can follow a Student's t distribution to handle outliers.

Usage

stanAM(
  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 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)

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

lmerAM, inlaAM, jagsAM

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 rstan
modelA$data$geno.add <- modelA$data$geno; modelA$data$geno <- NULL
fitA <- stanAM(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=1, Q=3, A=A, V.G.A=15, V.E=5,
                            D=D, V.G.D=3, seed=1859)

## infer with rstan
modelAD$data$geno.add <- modelAD$data$geno; modelAD$data$geno <- NULL
fitAD <- stanAM(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. 7, 2024, 8:17 a.m.