lmerAM: Animal model

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

Animal model

Description

Given I genotypes, Q covariates and N=I*Q phenotypes for the trait, fit an "animal model" with the lme4 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

lmerAM(
  formula,
  data,
  relmat,
  REML = TRUE,
  na.action = stats::na.exclude,
  ci.meth = NULL,
  ci.lev = 0.95,
  nb.boots = 10^3,
  parallel = "no",
  ncpus = 1,
  cl = NULL,
  verbose = 1
)

Arguments

formula	formula (see lmer)
data	data.frame containing the data corresponding to formula and relmat (see lmer); the additive genotypic effect should be named "geno.add" and the dominance genotypic effect, if any, should be named "geno.dom"
relmat	list containing the matrices of genetic relationships (A is compulsory but D is optional); the list should use the same names as the colnames in data (i.e., "geno.add" and "geno.dom") to compute heritability properly; the matrices can be in the "matrix" class (base) or the "dsCMatrix" class (Matrix package); see estimGenRel
REML	default is TRUE (use FALSE to compare models with different fixed effects)
na.action	a function that indicates what should happen when the data contain NAs (see lmer)
ci.meth	method to compute confidence intervals (profile/boot)
ci.lev	level to compute confidence intervals
nb.boots	number of bootstrap replicates; used only if ci.meth="boot"
parallel	the type of parallel operation to be used, if any (no/multicore/snow)
ncpus	number of processes to be used in parallel operation
cl	an optional object of class "cluster" returned by makeCluster; if not supplied, a cluster on the local machine is created for the duration of the call
verbose	verbosity level (0/1)

Value

list with a merMod object, a vector of point estimates of variance components, a point estimate of h2, a thpr object if ci.meth="profile", a boot object if ci.meth="profile", a data frame with confidence intervals (if ci.meth is not NULL)

Note

If A is not positive definite, an error will be raised (via chol); in such cases, using the nearPD function from the Matrix package can be useful.

Author(s)

Timothee Flutre (inspired by Ben Bolker at http://stackoverflow.com/q/19327088/597069)

See Also

inlaAM, jagsAM, 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)
kappa(A)
A <- as.matrix(nearPD(A)$mat)
kappa(A)
modelA <- simulAnimalModel(T=1, Q=3, A=A, V.G.A=15, V.E=5, seed=1859)

## infer with lme4
modelA$data$geno.add <- modelA$data$geno
fitA <- lmerAM(formula=response1 ~ year + (1|geno.add),
               data=modelA$data,
               relmat=list(geno.add=A), verbose=0)
summary(fitA$merMod)
REMLcrit(fitA$merMod)
extractAIC(fitA$merMod)
summary(residuals(fitA$merMod)) # "deviance residuals"
summary(residuals(fitA$merMod) / sigma(fitA$merMod)) # "scaled/Pearson residuals"
c(modelA$C); modelA$V.G.A; modelA$V.E
fixef(fitA$merMod)
coefficients(summary(fitA$merMod))[, "Std. Error"]
vc <- as.data.frame(VarCorr(fitA$merMod))
c(vc[vc$grp == "geno.add", "vcov"], vc[vc$grp == "Residual", "vcov"])
blups.geno <- ranef(fitA$merMod, condVar=TRUE, drop=TRUE)$geno.add
var.blups.geno <- setNames(attr(blups.geno, "postVar"), names(blups.geno))

## compute confidence intervals in parallel
library(parallel)
(nb.cores <- max(1, detectCores() - 1))
cl <- makeCluster(spec=nb.cores, type="PSOCK")
fitA <- lmerAM(formula=response1 ~ year + (1|geno.add),
               data=modelA$data,
               relmat=list(geno.add=A), verbose=1,
               ci.meth="profile", parallel="snow", ncpus=nb.cores, cl=cl)
stopCluster(cl)

## simulate phenotypes with additive and dominant parts of genotypic values
D <- estimGenRel(X, relationships="dominant", method="vitezica", verbose=0)
kappa(D)
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 lme4
modelAD$data$geno.add <- modelAD$data$geno
modelAD$data$geno.dom <- modelAD$data$geno; modelAD$data$geno <- NULL
fitAD <- lmerAM(formula=response1 ~ year + (1|geno.add) + (1|geno.dom),
                data=modelAD$data, relmat=list(geno.add=A, geno.dom=D),
                verbose=0)
summary(fitAD$merMod)
c(modelAD$C); modelAD$V.G.A; modelAD$V.E; modelAD$V.G.D
fixef(fitAD$merMod)
vc <- as.data.frame(VarCorr(fitAD$merMod))
c(vc[vc$grp == "geno.add", "vcov"], vc[vc$grp == "Residual", "vcov"],
  vc[vc$grp == "geno.dom", "vcov"])

## End(Not run)

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