Genetic basis of aggression-size relationship in voles

We now know vole size is heritable as well as being under positive selection and so evolution of larger size is predicted under simple models (univariate breeder’s equation). However, selection on body size may have evolutionary consequences for other aspects of phenotype, depending on whether or not they are genetically correlated. Here we will test whether there is a genetic basis to (hypothesised) covariance of body size and aggressiveness. A priori we might expect these traits to positively covary at the phenotypic level if (i) aggressive individuals are able to win more food resource in competition (and so grow faster), and/or (ii) expression of aggression towards conspecifics is mediated by relative size (ie bigger individuals can afford to act aggressively towards smaller rivals as they have lower risk of injury if a conflict escalates).

We’ll return to the data set we’ve already seen, extending our previous animal model analysis to the bivariate case. Since the free REML programs we’ve used already don’t do multivariate models we will switch to an MCMC approach using the package mcmcglmm. Note that we will be estimating a G matrix here, but for two traits only. Extension to any number of traits is possible in principle. In practice models with more traits do get harder to fit and take longer to run.

Getting ready!

Before we do anything else, set your working directory and reload the data files if needed.

#1. read in the main data file
voles<-read.table("voles2.txt", header=T)
voles$ID<-as.factor(voles$ID)
voles$sex<-as.factor(voles$sex)
voles$SIRE<-as.factor(voles$SIRE)
voles$DAM<-as.factor(voles$DAM)

#read in the pedigree file
volesPED<-read.table("volesPED.txt",header=T)

Then load up some packages we are going to use

#2. Automatically installs packages in this list that user does not already have

list_pkgs <- c("lme4","lmerTest","MCMCglmm","pedantics","pedigreemm")
new_pkgs <- list_pkgs[!(list_pkgs %in% installed.packages()[,"Package"])]
if(length(new_pkgs) > 0){ install.packages(new_pkgs) }

library(pedantics)
library(MCMCglmm)
library(lme4)
library(lmerTest)
library(pedigreemm)

An animal model with mcmcglmm

MCMCglmm is a Bayesian method - we need to specify priors or let it use its own defaults. In a way we might want to avoid using priors here, because we would like all of the information in our analysis to come from our data. So we are not really being Bayesian in the philosophical sense, just taking advantage of MCMC as a way to get our parameter estimates. Prior specification is a complex issue and we won’t get into it here. Jarrod Hadfield’s MCMCglmm course notes give a lot more information and the definitive explanation (a MUST read if you want to work with this package)

We are going to run a bivariate model, but just to see how mcmcglmm works let’s run a univariate model of aggression first to check it is heritable (if it isn’t then there’s little sense in estimating a genetic correlation with size).

#4. Univariate animal models in mcmcglmm

#first we'll set a prior for the variance components
prior_uni <- list(G = list(G1 = list(V = 1, nu = 0.002)), 
            R = list(V = 1, nu = 0.002))

#I think the package likes to reserve ID for an identity matrix so we’ll duplicate 
#our identity column and call it “animal” in the data file 
voles$animal<-as.factor(voles$ID)

#then run the model

model_uni <- MCMCglmm(aggression ~ 1+sex+forage, random = ~animal, pedigree = volesPED,
            data = voles, prior = prior_uni)
## 
##                       MCMC iteration = 0
## 
##                       MCMC iteration = 1000
## 
##                       MCMC iteration = 2000
## 
##                       MCMC iteration = 3000
## 
##                       MCMC iteration = 4000
## 
##                       MCMC iteration = 5000
## 
##                       MCMC iteration = 6000
## 
##                       MCMC iteration = 7000
## 
##                       MCMC iteration = 8000
## 
##                       MCMC iteration = 9000
## 
##                       MCMC iteration = 10000
## 
##                       MCMC iteration = 11000
## 
##                       MCMC iteration = 12000
## 
##                       MCMC iteration = 13000

We can use the summary function just as before to see the parameter estimates:

summary(model_uni)
## 
##  Iterations = 3001:12991
##  Thinning interval  = 10
##  Sample size  = 1000 
## 
##  DIC: 1996.729 
## 
##  G-structure:  ~animal
## 
##        post.mean l-95% CI u-95% CI eff.samp
## animal     1.589   0.9682    2.297    251.6
## 
##  R-structure:  ~units
## 
##       post.mean l-95% CI u-95% CI eff.samp
## units     2.354     1.79    2.856    451.3
## 
##  Location effects: aggression ~ 1 + sex + forage 
## 
##             post.mean  l-95% CI  u-95% CI eff.samp  pMCMC    
## (Intercept)  6.512771  6.181638  6.812078     1000 <0.001 ***
## sex2         0.310464  0.009252  0.662258     1000  0.056 .  
## forage      -0.184156 -0.349216 -0.029575     1000  0.022 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We now a set of fixed and random effect estimates, along with credible intervals based on the posteriors. If 95% CI for the fixed effects don’t overlap zero so we can conclude they are “significant” (but of course true Bayesians don’t care about significance).

Estimating h2 for aggression

A nice bonus of MCMC is that we can used the posteriors of estimated parameters to get a valid posterior for derived parameters. So here we might want to estimate h2 as VA / VP, and get 95% CI’s for that too.

#5.Get a posterior for h2
posterior.heritability<-model_uni$VCV[,1]/(model_uni$VCV[,1]+ model_uni$VCV[,2])

You can look at this posterior distribution in a histogram, pull out the mean or median as our point estimate, and get 95% CI. Be careful though, MCMCglmm constrains variance components to be positive, so the 95% CI’s for VA and h2 will never span zero (i.e. we cannot “test” these parameters in the frequentist sense).

mean(posterior.heritability)   
## [1] 0.4003053
HPDinterval(posterior.heritability)
##          lower     upper
## var1 0.2586752 0.5425361
## attr(,"Probability")
## [1] 0.95

Is it really that simple?

Sort of. However, using MCMC we really need to do some post-model fitting work to see whether we should believe the model output. For our posterior to be good we want to make sure the MCMC chain has converged and the posteriors of our estimates are stationary. If not we probably need to run the chain longer than the default, and/or change other options controlling the model fitting process. We would probably run teh model severl times to ensure our answers remained consistent (note they will not be identical as they would under REML as there is MCMC sampling variation here, not a single maximum likelihood solution) Being good Bayesians we should probably also try different prior specifications just so we can be sure our prior is not overly informing our posterior (i.e. what we get out is really from the data, not our prior belief of what we expect to get out).

Noting that the above points are very important for real analyses you want to publish, here we’ll just plot the posteriors to make a quick visual check on the MCMC run!

#6. Some diagnostic model plots
plot(model_uni$Sol)   #plots traces and posterior distributions of fixed effects

In the trace plots we are looking for consistent variation around a stable mean. Ideally the density plots show a clean unimodal (and roughly Gaussian) distribution. Take a look at the plots for the variance component estimates too.

plot(model_uni$VCV)   #plots traces and posterior distributions of random effect (co)variances

If you get unstable posterior means and/or strange looking distribution, the normal first step would be to change the mcmcglmm run options from the default values (total run length of 13000 iterations, the first 3000 of which are discarded as a `burn-in’, the remaining 10000 are sampled every 10 iterations (the thinning interval), to give a posterior size of 1000). For instance you could use a longer burn-in period and more iterations, and sample the non-discarded part of the chain less frequently which to reduce autocorrelation. Of course this will take longer to run.

#7. A longer run. Will you get the same answers?
model_uni2 <- MCMCglmm(aggression ~ 1+sex+forage, random = ~animal, pedigree = volesPED,
            data = voles, nitt = 65000, thin = 50, burnin = 15000, prior = prior_uni)

Fitting the bivariate model animal model

Now you’ve seen mcmcglmm let’s extend to the multivariate case to estimate the genetic variance-covariance matrix G for size and aggression. We’ll keep the extended length of the MCMC chain from above (indeed this may not be enough for a good solution here in which case you’ll all get quite different answers!).

#8. Fitting a bivariate model, first specify a prior

prior_bi <- list(G = list(G1 = list(V = diag(2), n = 1.002)),
                          R = list(V = diag(2), n = 1.002))


model_bi <- MCMCglmm(cbind(size, aggression) ~ trait + trait:sex +trait:forage- 1, 
                    random = ~us(trait):animal,
                    rcov = ~us(trait):units, family = c("gaussian", "gaussian"),
                    pedigree = volesPED, data = voles, prior = prior_bi, 
                    nitt = 65000, thin = 50, burnin = 15000, verbose = TRUE)
            #

A SLIGHT DIGRESSION - This probably took a minute (or several?) to run. A fast REML program would have found the maximum likelihood solution for the model in seconds. This illustrates that computing time is a relative disadvantage of MCMC. For small data sets and simple models like this it’s rarely an issue, for large data sets and complex models (e.g. lots af traits and or random effects) MCMC analyses can take days or weeks. However, you end up with posteriors that are really useful for carrying uncertainty forward into subsequent analyses. mcmcglmm also has much more advanced (and robust) capabilities for fitting generalised models than likelihood based methods. So there are advantages and disadvantages to consider!

Let’s look at the plots of the variance component esimates to see if we think the run will have generated reliable results. To reiterate this is not sufficient model checking for a formal analysis to generate publishable results Read Hadfield’s mcmcglmm vignette and course notes for what should properly be done. We are taking major short cuts here to save time.

#9. Plot variance covariance estimates for diagnostics 
plot(model_bi$VCV)   #plots traces and posterior distributions of random effect (co)variances

If you are happy with the plots then summarise the model to get parameter estimates.

#10. Summarise the bivariate model
summary(model_bi)
## 
##  Iterations = 15001:64951
##  Thinning interval  = 50
##  Sample size  = 1000 
## 
##  DIC: 3997.535 
## 
##  G-structure:  ~us(trait):animal
## 
##                              post.mean l-95% CI u-95% CI eff.samp
## size:size.animal                 2.342    1.596    3.229    849.1
## aggression:size.animal           1.668    1.058    2.246    908.0
## size:aggression.animal           1.668    1.058    2.246    908.0
## aggression:aggression.animal     1.619    1.017    2.332    795.3
## 
##  R-structure:  ~us(trait):units
## 
##                             post.mean l-95% CI u-95% CI eff.samp
## size:size.units                 2.566   2.0556   3.1471    844.2
## aggression:size.units          -0.261  -0.6723   0.1109   1000.0
## size:aggression.units          -0.261  -0.6723   0.1109   1000.0
## aggression:aggression.units     2.328   1.7923   2.7990    906.1
## 
##  Location effects: cbind(size, aggression) ~ trait + trait:sex + trait:forage - 1 
## 
##                        post.mean l-95% CI u-95% CI eff.samp  pMCMC    
## traitsize               20.88482 20.51353 21.21895   1000.0 <0.001 ***
## traitaggression          6.53254  6.20822  6.83346   1000.0 <0.001 ***
## traitsize:sex2           0.77219  0.45172  1.14197   1000.0 <0.001 ***
## traitaggression:sex2     0.32811  0.01642  0.65004    934.5  0.042 *  
## traitsize:forage         0.24815  0.07808  0.41458   1000.0  0.002 ** 
## traitaggression:forage  -0.17189 -0.33516 -0.01625   1000.0  0.038 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

This output contains estimates of the elements in G, the corresponding residual (environment) variance-covariance matrix (which we could call R), and the fixed effects. In all cases the posterior mean is provided as the point estimate.

Making a G matrix and calculating the genetic correlation (rG)

In the G-structure table section of the model summary, the elements are written in the order VA(size), COVA(size.agg), COVA(size.agg), VA(agg). Writing the covariance term twice seems a bit unecessary but actually makes it easier to read values into a matrix form, which we may want to do if using G in subsequent analyses.

#11. Make a G matrix and calcular rG

G_elements<-c(mean(model_bi$VCV[,1]), mean(model_bi$VCV[,2]), 
              mean(model_bi$VCV[,3]), mean(model_bi$VCV[,4]))
                        
              #model_bi$VCV contains the posteriors (each 1000 values) of the 
              #eight var-covar parameters (VA1,COVA12, COVA12, VA2, VR1, COVR12, COVR12, VR2).
              #Here we get the means of columns 1:4 to get point estimates of genetic 
              #parameters then arrange them in a matrix

G<-matrix(G_elements, nrow=2)

#Turn G into the correspoding correlation matrix
Gcor<-cov2cor(G)

Of course noting that rG =COVA/(VA1.VA2)0.5 we can also apply this calculation to the full posteriors of the genetic parameters, then use the posterior mean and 95 percentiles of the resulting posterior for rG. Since a correlation can be either positive or negative, whether the 95% CI span zero can be used to infer significance (or lack thereof) at the standard 0.05 level.

#11. Make a G matrix and calculate rG

rG_posterior<-model_bi$VCV[,2]/(model_bi$VCV[,1]*model_bi$VCV[,4])^0.5
rG_posterior.mean<-mean(rG_posterior)

rG_posterior.mean
## [1] 0.8606355
HPDinterval(rG_posterior)
##          lower     upper
## var1 0.7355148 0.9542086
## attr(,"Probability")
## [1] 0.95

Note that the rG posterior mean may not be exactly the same as the point estimate of the correlation as calculated using the posterior means of VA1, VA2 and COVA. However, it should be very close.