Simple univariate animal model

Fitting a simple univariate model in R.

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 ' ' 1

Using 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

Feedback

Was this page helpful?

Glad to hear it! Please tell us how we can improve.

Sorry to hear that. Please tell us how we can improve.

Last modified February 13, 2023: Merge pull request #43 from wamwiki/simple (f33c080)