Quick reference code
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
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
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)
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
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)
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
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.
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
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.
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)
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!