##############################################################################
# #
# Supplementary R code accompanying Postma et al. Proc. Roy. Soc. B 2011 #
# #
# Disentangling the effect of genes, the environment and chance on sex ratio #
# variation in a wild bird population. #
# #
# This script (including any improvements, bug fixes etc.) can also be #
# downloaded from http://www.erikpostma.net/resources.html #
# #
##############################################################################
##############################################################################
# #
# 1. Fitting a generalised mixed model to sex ratio data using MCMCglmm #
# #
# NOTE: Analyses were run with a range of priors. The settings that we found #
# to perform best with our data are likely to be different from those that #
# work best for other data sets. #
# #
# In order to get reliable results with other data sets, the MCMCglmm #
# code needs to be adjusted accordingly. Therefore we strongly recommend #
# to read the excellent and in-depth documentation that accompanies #
# MCMCglmm, and which can be found at: #
# http://cran.r-project.org/web/packages/MCMCglmm/index.html #
# #
# Finally, there are several ways of running these analyses. This is just #
# one of those, and it is not necessarily the fastest or most elegant. #
# #
##############################################################################
# Load MCMCglmm library
library("MCMCglmm")
# Set working directory (i.e. the location of pedigree and data file)
setwd("/my/working/directory/")
# Preparing the data
# This script assumes the existence of two files.
# The data file contains a line for each nestling, with the sex of each
# nestling ('sex'), coded as 0/1, the identity of their mother, once to
# account for additive genetic effects ('animal') and once to account for
# permanent environment effects ('pe'), and finally a brood ID ('brood').
data <- read.delim("sexratio_data.txt", header=TRUE)
data$animal <- as.factor(data$mother)
data$pe <- as.factor(data$mother)
data$brood <- as.factor(data$brood)
# The pedigree file contains three columns, 'id', 'mother' and 'father'
pedigree <- read.delim("pedigree.txt", header=TRUE)
pedigree$id <- as.factor(pedigree$id)
pedigree$mother <- as.factor(pedigree$mother)
pedigree$father <- as.factor(pedigree$father)
# We can reduce the size of the pedigree by removing any individuals
# for which no sex ratio data are available, and that do not connect
# individuals for which sex ratio data is available.
# Create vector with all animals for which sex ratio data is available
keep <- unique(data$animal)
# Use this to create a pruned pedigree
pedigree.pruned <- prunePed(pedigree, keep)
# Specifying and running the model
# Specify prior
# For more information, see MCMCglmm reference manual
prior = list(R = list(V = 10, fix = 1),
G = list(G1 = list(V = 1, nu = 1, alpha.mu = 0, alpha.V = 1000),
G2 = list(V = 1, nu = 1, alpha.mu = 0, alpha.V = 1000),
G3 = list(V = 1, nu = 1, alpha.mu = 0, alpha.V = 1000)))
# Run the model with an intercept and three random effects
model <- MCMCglmm(molsex ~ 1, random = ~ brood + pe + animal, family = "categorical",
data = data, pedigree = pedigree.pruned, prior = prior,
verbose = TRUE, slice=TRUE,
burnin = 50000, nitt = 1050000, thin = 100)
# Process the output. Again, see the MCMCglmm documentation for more details
# Mean probability of producing a male (i.e. mean sex ratio)
effectiveSize(model$Sol)
autocorr(model$Sol)
plot(model$Sol)
c2 <- ((16*sqrt(3))/(15*pi))^2
intercept <- model$Sol/sqrt(1 + c2 * model$VCV[, "units"])
inv.logit(posterior.mode(intercept)) # Mean p
inv.logit(HPDinterval(intercept, 0.95)) # 95% credible interval
# Variance components
effectiveSize(model$VCV)
autocorr(model$VCV)
plot(model$VCV)
posterior.mode(model$VCV)
HPDinterval(model$VCV, 0.95)
# Heritability
brood <- model$VCV[,1]
pe <- model$VCV[,2]
animal <- model$VCV[,3]
units <- model$VCV[,"units"]
h2 <- animal/(brood + pe + animal + units + (pi^2)/3)
plot(h2)
posterior.mode(h2)
HPDinterval(h2, 0.95)