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()