simulAnimalModel: Animal model

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

Animal model

Description

Given T traits, I genotypes, Q covariates and N=I*Q phenotypes per trait, simulate phenotypes via the following "animal model": Y = W C + Z G_A + Z G_D + E, where Y is N x T; W is N x Q; Z is N x I; G_A ~ Normal_IxT(0, A, V_G_A) with A being IxI and V_G_A being TxT; G_D ~ Normal_IxT(0, D, V_G_D) with D being IxI and V_G_D being TxT; E ~ Normal_NxT(0, Id_N, V_E) with Id_N being NxN and V_E being TxT. Note that using the matrix Normal (MN) is equivalent to using the multivariate Normal (MVN) via the vec operator and Kronecker product: Y ~ MN(M, U, V) <=> vec(Y) ~ MVN(vec(M), V \otimes U). Spatial heterogeneity can be add (see Dutkowski et al (2002)).

Usage

simulAnimalModel(
  T = 1,
  Q = 3,
  mu = rep(50, T),
  mean.C = 5,
  sd.C = 2,
  A,
  V.G.A = NULL,
  scale.hC.G.A = NULL,
  nu.G.A = T,
  D = NULL,
  V.G.D = NULL,
  scale.hC.G.D = NULL,
  nu.G.D = T,
  V.E = NULL,
  scale.hC.E = NULL,
  nu.E = T,
  err.df = Inf,
  perc.NA = 0,
  seed = NULL,
  nb.rows = NULL,
  nb.cols = NULL,
  sigma2.ksi = NULL,
  rho.rows = NULL,
  rho.cols = NULL
)

Arguments

T	number of traits
Q	number of covariates (as "fixed effects", e.g. replicates)
mu	T-vector of overall means (one per trait), i.e. C[1,1:T]
mean.C	mean of the univariate Normal prior on C[2:Q,1:T] (ignored if Q=1)
sd.C	std dev of the univariate Normal prior on C[2:Q,1:T] (ignored if Q=1)
A	IxI matrix of additive genetic relationships between genotypes (see estimGenRel with VanRaden's estimator); if A is singular (i.e. has a large condition number), the Choleski decomposition will be performed with pivoting
V.G.A	scalar (if T=1) or TxT matrix of additive genetic variance-covariance between traits (e.g. 15 when T=1)
scale.hC.G.A	scale of the half-Cauchy prior for sqrtV_G_A (e.g. 5; used if V.G.A=NULL and T=1)
nu.G.A	degrees of freedom of the Wishart prior for V_G_A (used if V.G.A=NULL and T>1)
D	IxI matrix of dominant genetic relationships between genotypes (see estimGenRel with Vitezica's estimator); if D is singular (i.e. has a large condition number), the Choleski decomposition will be performed with pivoting
V.G.D	scalar (if T=1) or TxT matrix of dominant genetic variance-covariance between traits (e.g. 3 when T=1; used if D!=NULL)
scale.hC.G.D	scale of the half-Cauchy prior for sqrtV_G_D (e.g. 5; used if D!=NULL, V.G.D=NULL and T=1)
nu.G.D	degrees of freedom of the Wishart prior for V_G_D (used if D!=NULL, V.G.D=NULL and T>1)
V.E	scalar (if T=1) or TxT matrix of error variance-covariance between traits (used if T=1 and err.df=Inf)
scale.hC.E	scale of the half-Cauchy prior for sqrtV_E (e.g. 5; used if V.E=NULL and T=1 and err.df=Inf)
nu.E	degrees of freedom of the Wishart prior for V_E (used if V.E=NULL and T>1)
err.df	degrees of freedom of the Student's t-distribution of the errors (e.g. 3; Inf means Normal distribution; will be Inf if T>1)
perc.NA	percentage of missing phenotypes, at random
seed	seed for the pseudo-random number generator
nb.rows	number of rows, when assuming data come from a plant species with a regular spatial conformation
nb.cols	number of columns (nb.rows x nb.columns should be equal to the number of rows of A)
sigma2.ksi	variance of the spatially-correlated errors (only if T=1)
rho.rows	correlation between neighbor plants on the same row (only if T=1)
rho.cols	correlation between neighbor plants on the same column (only if T=1)

Value

list

Author(s)

Timothee Flutre

Examples

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

## estimate the additive genetic relationships
A <- estimGenRel(X, relationships="additive", method="vanraden1", verbose=0)

# 1: one trait with heritability h2=0.75, no covariate, Normal errors, no NA
model <- simulAnimalModel(T=1, Q=1, A=A, V.G.A=15, V.E=5, seed=1859)

# 2: one trait with heritability h2=0.75, three covariates, Normal errors, no NA
model <- simulAnimalModel(T=1, Q=3, A=A, V.G.A=15, V.E=5, seed=1859)

# 3: one trait with heritability h2=0.75, no covariate, Student errors, no NA
model <- simulAnimalModel(T=1, Q=1, A=A, V.G.A=15, err.df=3, seed=1859)

# 4: one trait with heritability drawn at random, no covariate, Normal errors, no NA
model <- simulAnimalModel(T=1, Q=1, A=A, scale.hC.G.A=10, V.E=5, seed=1859)

# 5: two traits with heritabilities drawn at random, no covariate, Normal errors, no NA
model <- simulAnimalModel(T=2, Q=1, A=A, nu.G.A=5, nu.E=3, seed=1859)

# 6: same as scenario 1 but with spatial coordinates
model <- simulAnimalModel(T=1, Q=1, A=A, V.G.A=15, V.E=5, seed=1859, nb.rows=20, nb.cols=10)
## if(require(sp))
##   plot(SpatialPoints(model$coords), axes=TRUE, las=1,
##        xlab="columns", ylab="rows", main="Map")

## End(Not run)

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