# Quick reference code

### Overview

This univariate model is the simplest of all animal models and assumes no repeated measures across years or individuals. Here we provide the code to calculate the additive genetic variance component of a phenotypic variable, for example body size.

If you havenâ€™t used these programs before, this information will not make a huge amount of sense. We strongly recommend that you read the additional background information and step-by-step instructions which can be found in the freely available ecologists guide to the animal model, published in the Journal of Animal Ecology and available here. Have a look at our quick presentation of the software we use.

# 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

## Using MCMCglmm

Here is the simplest implementation of an animal model:

```
#set the prior
prior1.1<-list(G=list( G1=list(V=1,n=0.002)), R=list(V=1,n=0.002))
#model specification
model1.1<-MCMCglmm(SIZE~1,random=~ANIMAL,
pedigree=pedigreedata,data=phenotypicdata,prior=prior1.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
#1: SIZE is the response variable and the only fixed effect is the mean(denoted as1)
#2: fit random effect of ANIMAL Va with an arbitrary starting value of 1
#3: use data file phenotypic data
#4: connect the individual in the data file to the pedigree
#5: omit any rows where the response or predictor variables are missing
#to see the estimates of the fixed effects
summary(model)$coef.fixed
#and the estimates of the random effects
summary(model)$varcomp
```

Or with the standaloneversion:

```
# provide a title - failure to do this will prevent the program from running properly
ASReml analysis of size
# next section describes contents of data file
# note that column headings are indented (use at least one space)
# !P associates a term with the pedigree file
ANIMAL !P
SIZE
# then specify the pedigree and data files
# for the data file we tell ASReml to skip the first line since these are headers not data
# Note that the pedigree and data file names below should not be indented
pedigreedata.ped !skip 1
phenotypicdata.dat !skip 1
# Then we specify the model
# mu is the mean & any term after !r is treated as random
SIZE ~ mu !r ANIMAL
```

# Adding fixed effects to a model

## Using MCMCglmm

```
model1.2<-MCMCglmm(BWT~SEX,random=~animal,
pedigree=Ped,data=Data,prior=prior1.1,
nitt=65000,thin=50,burnin=15000,verbose=FALSE)
posterior.mode(model1.2$Sol[,"SEX2"])
HPDinterval(model1.2$Sol[,"SEX2"],0.95)
posterior.mode(model1.2$VCV)
posterior.heritability1.2<-model1.2$VCV[,"animal"]/
(model1.2$VCV[,"animal"]+model1.2$VCV[,"units"])
posterior.mode(posterior.heritability1.2)
HPDinterval(posterior.heritability1.2,0.95)
```

## Using ASReml

```
ASReml analysis of size
ANIMAL !P
SIZE
SEX !A #denotes a factor
AGE #here treated as a linear effect
pedigreedata.ped !skip 1
phenotypicdata.dat !skip 1 !DDF 1 !FCON #specifies method of sig testing
SIZE ~ mu AGE SEX AGE.SEX !r ANIMAL
```

# Testing significance of random effects

## Using MCMCglmm

MCMCglmm, and more in general Bayesian methods, do not provide a simple, consensual way to test the statistical significance of a variance parameters. Indeed, variances parameters are constrained to be positive, and their credible intervals (e.g., as returned by HPDinterval()) cannot include exactly zero (although it may look like it due to rounding. Covariance and correlation parameters do not have this issue because they are not constrained to be positive and their credible interval can be used to estimate the probability that they are positive or negative. The old WAMBAM website recommended to compare DIC (Deviance Information Criterion, analog to AIC) across models with and without a random effect. However, DIC may be focused at different levels of a mixed model, and is calculated for the lowest level of the hierarchy in MCMCglmm, which may not be appropriate for comparing different random effect structures.

## Using ASReml

In ASReml statistical the significance of a variance parameter can be tested using a Likelihood Ratio Test. Fit a model with and without a particular random effect. Then use log likelihoods reported in the primary results file to perform a ratio test.

In ASReml standalone:

```
ASReml analysis of size
ANIMAL !P
SIZE
SEX !A #sex as a factor
pedigreedata.ped !skip 1
phenotypicdata.dat !skip 1 !dopart 1 #change part to be run required
!part 1
SIZE ~ mu SEX !r ANIMAL RANDOMEFFECT
!part 2
SIZE ~ mu SEX !r ANIMAL
```

From R:

```
model1<-asreml(fixed=SIZE~1+SEX
,random=~ped(ANIMAL,var=T,init=1)+RANDOMEFFECT
,data=phenotypicdata
,ginverse=list(ANIMAL=ainv), na.method.X="omit', na.method.Y="omit')
model2<-asreml(fixed=SIZE~1+SEX
,random=~ped(ANIMAL,var=T,init=1)
,data=phenotypicdata
,ginverse=list(ANIMAL=ainv), na.method.X="omit', na.method.Y="omit')
#calculate the chi-squared stat for the log-likelihood ratio test
2*(model1$loglik-model2$loglik)
#calculate the associated significance
1-pchisq(2*(model1$loglik-model2$loglik),df=1)
```

However, this test is conservative with 1 degree of freedom. Using df=0.5 gives a better (but still a bit conservative) test.

# Calculating heritability

## Using MCMCglmm

```
posterior.heritability1.1<-model1.1$VCV[,"animal"]/
(model1.1$VCV[,"animal"]+model1.1$VCV[,"units"])
HPDinterval(posterior.heritability1.1,0.95)
posterior.mode(posterior.heritability1.1)
plot(posterior.heritability1.1)
```

## Using ASReml

In ASReml standalone: In ASReml a second command file (with extension .pin) is used to caculate functions of estimated variance components ad their associated standard errors. So for a model in the .as file such as

```
SIZE ~ mu ! ANIMAL
```

the primary output file (.asr) will contain two variance components. The first will be the ANIMAL (i.e. additive genetic component), the second will be the residual variance. A .pin file to calculate heritability from these components migt be

```
F VP 1+2 #adds components 1 and 2 to make a 3rd variance denoted VP
H h2 1 3 #divides 1 (VA) by 3 (VP) to calculate h2
```

NOTE - if you change the random effects stucture of your model in .as you need to modify the .pin file accordingly or you will get the wrong answer!

From R:

```
summary(model)$varcomp[1,3]/sum(summary(model)$varcomp[,3])
```