Here we demonstrate how to fit the simplest univariate animal model to calculate the additive genetic variance component of a phenotypic variable, for example body size. In later sections we show how to build more complex models with fixed or random effects and how to compute heritability.
is the simplest of all animal models and assumes no repeated measures across years or individuals.
Example data
For this tutorial we will use the simulated gryphon dataset (download zip file).
phenotypicdata <- read.csv("data/gryphon.csv")
pedigreedata <- read.csv("data/gryphonped.csv")The univariate animal model
Definitions
Sample code in the following examples includes the following variables:
- Response variable: Size
- Fixed effects: Intercept (\(\mu\))
- Random effects: Additive genetic variance
- Data containing phenotypic information: phenotypicdata
- Data containing pedigree data:pedigreedata
We first fit the simplest possible animal model: no fixed effects apart from the interecept, a single random effect (the breeding values, associated with the additive genetic variance), and Gaussian redisuals.
Here it is not too difficult to describe the model in plain English, but for more complex models a mathematical equation may prove easier, much shorter and less ambiguous. Here we could write the model as (y_i = \mu + a_i + \epsilon_i), where $y_i $ is the response for individual (i), with the residuals assumed to follow a normal distribution of mean zero and variance (\sigma_R^2), which we can write as (\mathbf{\epsilon} \sim N(0,\mathbf{I}\sigma_R^2)), where (\mathbf{I}) is an identity matrix; and with the breeding values (\mathbf{a} \sim N(0,\mathbf{A} \sigma_A^2), where (\mathbf{A}) is the pairwise additive genetic relatedness matrix.
Using MCMCglmm
Here is the simplest implementation of an animal model in MCMCglmm:
library(MCMCglmm) #load the package
inverseAmatrix <- inverseA(pedigree = pedigreedata)$Ainv # compute the inverse relatedness matrix
model1.1<-MCMCglmm(birth_weight~1, #Response and Fixed effect formula
random=~id, # Random effect formula
ginverse = list(id=inverseAmatrix), # correlations among random effect levels (here breeding values)
data=phenotypicdata)# data set
Note the use of the argument ginverse
summary(model1.1)##
## Iterations = 3001:12991
## Thinning interval = 10
## Sample size = 1000
##
## DIC: 3908.779
##
## G-structure: ~id
##
## post.mean l-95% CI u-95% CI eff.samp
## id 3.454 2.209 4.72 170
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 3.822 2.755 4.787 209.8
##
## Location effects: birth_weight ~ 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 7.594 7.300 7.863 1000 <0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Using ASReml
From R:
model1<-asreml(fixed=SIZE~ 1 #1
, random= ~ped(ANIMAL,var=T,init=1) #2
, data=phenotypicdata #3
,ginverse=list(ANIMAL=pedigreedata_inverse) #4
, na.method.X="omit", na.method.Y="omit") #5