Text S4. R Code for Data Example Using the Between-Individual Covariance Matrix ID As

Text S4. R code for data example using the between-individual covariance matrix ID as a proxy for G

The code below reads in the data file and runs a Bayesian multivariate mixed model on one morphometric trait (tarsus length tarsus_avg) and three behavioral traits (activity sqrt_act, aggression aggres and breath rate br.rate). The design was that a large number of adults were repeatedly measured (in different seasons and different years) for these traits (see Kluen et al. 2014). The between-individual covariance matrix (ID matrix) can then be separated from the residual covariance matrix. The second part sums up some summary 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 ID matrix. That is, matrix ID is used as a proxy for G in the evaluation of average autonomy.

#get the data

dat<-read.table("ASremlIG_file.txt", sep="\t", header=TRUE,na.strings="-999")

require(MCMCglmm)

pheno.var<-apply(cbind(dat$tarsus_avg, dat$sqrt_act, dat$aggres, dat$br_rate),2,function(m) var(m,na.rm=T))

prior4T<-list(R=list(V=diag(c(pheno.var))*0.7, nu=3), G=list(G1=list(V=diag(c(pheno.var))*0.3, nu=3) ))

mod.ID <- MCMCglmm(cbind(tarsus_avg, sqrt_act, aggres, br_rate) ~ (trait - 1)+ trait*time + trait*as.factor(sex) + trait* as.factor(age) + trait*as.factor(year_int)+trait*as.factor(season)+trait*who, random = ~us(trait):ring ,rcov = ~us(trait):units, family = c("gaussian", "gaussian", "gaussian", "gaussian"), data=dat, prior = prior4T, verbose = FALSE,burnin=1000,nitt=51000,thin=50)

#some model investigations

plot(mod.ID$VCV)

autocorr(mod.ID$VCV[,1:16])

#between-individual correlations

ID.corr.mat<-matrix(posterior.mode(posterior.cor(mod.ID$VCV[,1:16])),4,4,dimnames=list(c("Tars","Act","Agg","BR"),c("Tars","Act","Agg","BR")))

lower.ID.corr.mat<-matrix(HPDinterval(posterior.cor(mod.ID$VCV[,1:16]))[,1],4,4,dimnames=list(c("Tars","Act","Agg","BR"),c("Tars","Act","Agg","BR")))

upper.ID.corr.mat<-matrix(HPDinterval(posterior.cor(mod.ID$VCV[,1:16]))[,2],4,4,dimnames=list(c("Tars","Act","Agg","BR"),c("Tars","Act","Agg","BR")))

#residual correlations

R.corr.mat<-matrix(posterior.mode(posterior.cor(mod.ID$VCV[,17:(2*16)])),4,4,dimnames=list(c("Tars","Act","Agg","BR"),c("Tars","Act","Agg","BR")))

lower.R.corr.mat<-matrix(HPDinterval(posterior.cor(mod.ID$VCV[,17:(2*16)]))[,1],4,4,dimnames=list(c("Tars","Act","Agg","BR"),c("Tars","Act","Agg","BR")))

upper.R.corr.mat<-matrix(HPDinterval(posterior.cor(mod.ID$VCV[,17:(2*16)]))[,2],4,4,dimnames=list(c("Tars","Act","Agg","BR"),c("Tars","Act","Agg","BR")))

#################################################

# observed and expected average autonomy

#################################################

#initialise

setwd("C:/Users/joegbr/Documents/data/Tutkimus/projects/MethodsPapers/autonomy")

gen.a <- dget("gena.R")

autonomy.ID<-vector("numeric",dim(mod.ID$VCV)[1])

Null.a.ID<-vector("numeric",dim(mod.ID$VCV)[1])

ID.stand<-matrix(NA,dim(mod.ID$VCV)[1],36)

#calculate

for (i in 1:dim(mod.ID$VCV)[1]){

## STANDARDISATION by trait.specific SD

id <- matrix(mod.ID$VCV[i,1:16], 4, 4)

res <- matrix(mod.ID$VCV[i,17:(2*16)], 4, 4)

#vector with phenotypic standard deviation for the i'th posterior

#calculated as the sum of all variance components in the model

Psd=sqrt(diag(id)+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 ID matrix of the i'th posterior

ID.stand[i,1:16]<-as.vector(id/SquaredPMatrix)

#autonomy for the standardized ID matrix of the i'th posterior

autonomy.ID[i]<-gen.a(id/SquaredPMatrix)

Null.a.ID[i]<-gen.a((id/SquaredPMatrix)*diag(4))

}

#point estimates for the variance standardised ID matrix with its 95% lower and upper credible interval

ID.obs.stand<-matrix(posterior.mode(as.mcmc(ID.stand[,1:16])),4,4,dimnames=list(c("Tars","Act","Agg","BR"),c("Tars","Act","Agg","BR")))

lo.ID.obs.stand<-matrix(HPDinterval(as.mcmc(ID.stand[,1:16]))[,1],4,4,dimnames=list(c("Tars","Act","Agg","BR"),c("Tars","Act","Agg","BR")))

hi.ID.obs.stand<-matrix(HPDinterval(as.mcmc(ID.stand[,1:16]))[,2],4,4,dimnames=list(c("Tars","Act","Agg","BR"),c("Tars","Act","Agg","BR")))

#stats regarding the average autonomy

posterior.mode(as.mcmc(autonomy.ID)) #point estimate observed autonomy

HPDinterval(as.mcmc(autonomy.ID),prob=0.9) #and its 90% CRI

posterior.mode(as.mcmc(Null.a.ID)) #expected autonomy in case cov's are 0

HPDinterval(as.mcmc(Null.a.ID),prob=0.9) #and its 90% CRI

# ¤¤ plotting

tiff(file="FigS2.tif", width = 8, height = 5, units = 'in', res = 300,compression = "lzw" )

hist(autonomy.ID,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,200),freq=TRUE)

lines(c(HPDinterval(as.mcmc(Null.a.ID),prob=0.9)[1],HPDinterval(as.mcmc(Null.a.ID),prob=0.9)[1]),c(0,200),lwd=2,col="red",lty="dashed")

lines(c(HPDinterval(as.mcmc(autonomy.ID),prob=0.9)[2],HPDinterval(as.mcmc(autonomy.ID),prob=0.9)[2]),c(0,200),lwd=2)

dev.off()