Text S3. R code for data example using the animal model
The code below reads in the data file and the pedigree file and runs a Bayesian multivariate mixed model on three behavioral traits. The design of this experiment was such that nestlings were reciprocally cross-fostered between pairs of nests (see Brommer and Kluen 2012). Hence, there are a number of variance components estimable which are included to get an as precise as possible estimate of the additive genetic variance. The second part sums up some the aspects of the different (co)variance matrices. The third part applies the above outlined approach for calculating observed and expected average autonomy on the basis of the posterior matrices, while carrying out variance standardization of the G matrix.
############################################################
## DATA EXAMPLE
############################################################
gen.a <- dget("gena.R") #reads autonomy function, see supplement
require(MCMCglmm)
#read data
data2=read.table("young_CF_07_09.txt",header=T,na.strings = "-999.00")
#read pedigree
pedigree.data=read.table("prunedPed.txt",na.strings = "0")
names(pedigree.data)=c("animal","sire","dam")
pedigree.data$animal<-as.factor(pedigree.data$animal)
pedigree.data$sire<-as.factor(pedigree.data$sire)
pedigree.data$dam<-as.factor(pedigree.data$dam)
##
pheno.var<-apply(cbind(data2$Agg_d16, data2$DOC, data2$BR),2,function(m) var(m,na.rm=T))
#define priors
prior3T<-list(R=list(V=diag(c(pheno.var))*0.4, nu=2), G=list(G1=list(V=diag(c(pheno.var))*0.2, nu=2),G2=list(V=diag(c(pheno.var))*0.2, nu=2),G3=list(V=diag(c(pheno.var))*0.2, nu=2) ) )
#run model
mod <- MCMCglmm(cbind(Agg_d16, DOC, BR) ~ (trait - 1) + trait*as.factor(Year) + trait*Sex + trait*obs_d16, random = ~us(trait):RearID + us(trait):GeneticID + us(trait):animal, rcov = ~us(trait):units, family = c("gaussian", "gaussian", "gaussian"), pedigree=pedigree.data,data=data2, prior = prior3T, verbose = FALSE,burnin=10000,nitt=510000,thin=500)
#below can be used to explore stats of the MCMC model see also the MCMCglmm manuals
plot(mod$VCV)
autocorr(mod$VCV[,1:9])
autocorr(mod$VCV[,10:18])
autocorr(mod$VCV[,19:27])
autocorr(mod$VCV[,28:36])
###################################################################################
###### Part 2: construct matrices ###################################################################################
#the point estimates for the correlation matrices with their lower and upper 95%CI
cor.matrix.rear<-matrix(posterior.mode(posterior.cor(mod$VCV[,1:9])),3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
lo.cor.matrix.rear<-matrix(HPDinterval(posterior.cor(mod$VCV[,1:9]))[,1],3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
hi.cor.matrix.rear<-matrix(HPDinterval(posterior.cor(mod$VCV[,1:9]))[,2],3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
cor.matrix.origin<-matrix(posterior.mode(posterior.cor(mod$VCV[,10:18])),3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
lo.cor.matrix.origin<-matrix(HPDinterval(posterior.cor(mod$VCV[,10:18]))[,1],3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
hi.cor.matrix.origin<-matrix(HPDinterval(posterior.cor(mod$VCV[,10:18]))[,2],3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
cor.matrix.gen<-matrix(posterior.mode(posterior.cor(mod$VCV[,19:27])),3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
lo.cor.matrix.gen<-matrix(HPDinterval(posterior.cor(mod$VCV[,19:27]))[,1],3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
hi.cor.matrix.gen<-matrix(HPDinterval(posterior.cor(mod$VCV[,19:27]))[,2],3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
cor.matrix.res<-matrix(posterior.mode(posterior.cor(mod$VCV[,28:36])),3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
lo.cor.matrix.res<-matrix(HPDinterval(posterior.cor(mod$VCV[,28:36]))[,1],3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
hi.cor.matrix.res<-matrix(HPDinterval(posterior.cor(mod$VCV[,28:36]))[,2],3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
#the point estimates for the covariance matrices with their lower and upper 95%CI
cov.matrix.rear<-matrix((posterior.mode(mod$VCV[,1:9])),3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
lo.cov.matrix.rear<-matrix(HPDinterval((mod$VCV[,1:9]))[,1],3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
hi.cov.matrix.rear<-matrix(HPDinterval((mod$VCV[,1:9]))[,2],3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
cov.matrix.origin<-matrix(posterior.mode((mod$VCV[,10:18])),3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
lo.cov.matrix.origin<-matrix(HPDinterval((mod$VCV[,10:18]))[,1],3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
hi.cov.matrix.origin<-matrix(HPDinterval((mod$VCV[,10:18]))[,2],3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
cov.matrix.gen<-matrix(posterior.mode((mod$VCV[,19:27])),3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
lo.cov.matrix.gen<-matrix(HPDinterval((mod$VCV[,19:27]))[,1],3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
hi.cov.matrix.gen<-matrix(HPDinterval((mod$VCV[,19:27]))[,2],3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
cov.matrix.res<-matrix(posterior.mode((mod$VCV[,28:36])),3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
lo.cov.matrix.res<-matrix(HPDinterval((mod$VCV[,28:36]))[,1],3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
hi.cov.matrix.res<-matrix(HPDinterval((mod$VCV[,28:36]))[,2],3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
#################################################
# Part 3: observed and expected average autonomy
#################################################
#initialise
autonomy<-vector("numeric",dim(mod$VCV)[1])
Null.a<-vector("numeric",dim(mod$VCV)[1])
G.stand<-matrix(NA,dim(mod$VCV)[1],9)
#calculate
for (i in 1:dim(mod$VCV)[1]){
## STANDARDISATION by trait.specific SD
rear <-matrix(mod$VCV[i,1:9],3,3)
origin <- matrix(mod$VCV[i,10:18], 3, 3)
gen <- matrix(mod$VCV[i,19:27], 3, 3)
res <- matrix(mod$VCV[i,28:36], 3, 3)
#vector with phenotypic standard deviation for the i'th posterior
#calculated as the sum of all variance components in the model
Psd=sqrt(diag(rear)+diag(origin)+diag(gen)+diag(res))
#turn the phenotypic SD into cross-product matrix form
SquaredPMatrix<-matrix(rep(Psd,length(Psd)),length(Psd),length(Psd),byrow=TRUE)*matrix(rep(Psd,length(Psd)),length(Psd),length(Psd),byrow=FALSE)
#save the standardised G matrix of the i'th posterior
G.stand[i,1:9]<-as.vector(gen/SquaredPMatrix)
#autonomy for the standardized G matrix of the i'th posterior
autonomy[i]<-gen.a(gen/SquaredPMatrix)
Null.a[i]<-gen.a((gen/SquaredPMatrix)*diag(3))
}
#point estimates for the variance standardised G matrix with its 95% lower and upper credible interval
G.obs.stand<-matrix(posterior.mode(as.mcmc(G.stand[,1:9])),3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
lo.G.obs.stand<-matrix(HPDinterval(as.mcmc(G.stand[,1:9]))[,1],3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
hi.G.obs.stand<-matrix(HPDinterval(as.mcmc(G.stand[,1:9]))[,2],3,3,dimnames=list(c( "AGG", "DOC", "BR"),c( "AGG", "DOC", "BR")))
#stats regarding the average autonomy
posterior.mode(as.mcmc(autonomy)) #point estimate observed autonomy
HPDinterval(as.mcmc(autonomy),prob=0.9) #and its 90% CRI
mean(Null.a) #expected autonomy in case cov's are 0
HPDinterval(as.mcmc(Null.a),prob=0.9) # and its 90%CI
#¤¤¤¤¤¤
#¤¤¤¤¤¤ plotting the supplementary figure for illustration
tiff(file="FigS1.tif", , width = 8, height = 5, units = 'in', res = 300,compression = "lzw" )
hist(autonomy,main=NA,xlab=list("Posterior obs. average autonomy", cex=1.5),ylab=list("Count",cex=1.5),breaks=seq(0,1,0.02),col='grey',xlim=c(0,1),ylim=c(0,80),freq=TRUE)
lines(c(HPDinterval(as.mcmc(Null.a),prob=0.9)[1],HPDinterval(as.mcmc(Null.a),prob=0.9)[1]),c(0,80),lwd=2,col="red",lty="dashed")
lines(c(HPDinterval(as.mcmc(autonomy),prob=0.9)[2],HPDinterval(as.mcmc(autonomy),prob=0.9)[2]),c(0,80),lwd=2)
dev.off()