Barbara Class, Jon. E. Brommer, “Senescence of personality in a wild bird”
In Behavioral Ecology and Sociobiology
Corresponding author:
Text S3: R script for the simulation comparing model output when using within-individual vs. grand-mean centering of data prior to random-regression analysis.
#Create the data frame which will store the estimates of the K matrix
setwd("~/GXA in BT")
K.estimates<-data.frame("V1"=NA, "V2"=NA, "V3"=NA, "Conv"=NA, "Chisq"=NA,
"V1b"=NA, "V2b"=NA, "V3b"=NA, "Conv.b"=NA, "Chisq.b"=NA)
library(asreml)
ped<-read.table("Pedigree03_14.txt",header=T, na.strings="0");
ainv<-asreml.Ainverse(ped)$ginv
######Start of the simulation
for (z in 1:1000){
#Use individuals that were measured more than once
setwd("~/GXA in BT")
adults.full<- read.table("GAI_HA_full_old.txt", header=T, na.strings=c("-999"))
adults.full<-adults.full[is.na(adults.full$HA)==FALSE,]
table1<-aggregate(adults.full$HA, list(adults.full$Ring),length)
colnames(table1)<-c("Ring","n")
multiple<-table1[table1$n>=2,]
ring.multi<-multiple$Ring #406 individuals
adults.full_multi<-adults.full[adults.full$Ring %in% ring.multi,c(1,2)]
#Simulate 406 intercepts and slopes with their variance and covariance derived from the #RRAM analysis
library(mvtnorm)
K <- matrix(c(0.5,-0.3 ,-0.3,0.4), ncol=2)
int.slope<- rmvnorm(n=406, mean=c(3,-0.1), sigma=K)
colMeans(int.slope)
var(int.slope)
colnames(int.slope)<-c("Intercept", "Slope")
#Calculate min and max age of measurement for each individual
table2<-aggregate(adults.full_multi$Age, list(adults.full_multi$Ring), max)
table2$mean.age<-(table2[,2]+1)/2
colnames(table2)<-c("Ring","max.age","mean.age")
table3<-aggregate(adults.full_multi$Age, list(adults.full_multi$Ring), min)
colnames(table3)<-c("Ring","min.age")
table4<-merge(table2, table3)
#Associate each individual with an intercept and a slope
Indiv.values<-cbind(table4,int.slope)
adults.full_multi<-merge(adults.full_multi,Indiv.values,all=T)
#Add a column with individual centered age and grand mean centered age
adults.full_multi$Age_F<-as.factor(adults.full_multi$Age)
adults.full_multi$Icent<-adults.full_multi[,2]-adults.full_multi[,4]
adults.full_multi$Gcent<-adults.full_multi[,2]-4
#Calculate individual centered age scaled by asreml
scale<-data.frame("Icent"=seq(from=-3, to=+3, by=0.5), "asreml.scale"=c(-1, -0.83466, -0.66931, -0.50397, -0.33862, -0.17328, -0.00794, 0.15741, 0.32275, 0.4881, 0.65344, 0.81878 ,0.98413))
adults.full_multi<-merge(adults.full_multi,scale,all=T)
adults.full_multi<-adults.full_multi[,c(2,3,4,5,6,7,8,9,1,10,11)]
#Fill the HA column for each individual for each age at which it was measured
adults.full_multi$HA<-NA
adults.full_multi$res<-rnorm(1042,0,0.6)
#HA is calculated for each age (intercept is HA for mean age (Icent=0)based on asreml scaled age)
for(i in 1:406){
a<-adults.full_multi[adults.full_multi$Ring==as.character(ring.multi[i]),]
n<-a$Age
for (y in 1:length(n)){
adults.full_multi[adults.full_multi$Ring==as.character(ring.multi[i]) & adults.full_multi$Age== n[y], 12]<- a[a$Age==n[y],6]+(a[a$Age==n[y],11] *a[a$Age==n[y],7])+a[a$Age==n[y],13]
}
}
####Analyze it using RRAM model
model.1<-asreml(fixed=HA~ 1+ Gcent#Grand-mean centered age, Model I
,random=~ide(Ring)
,data=adults.full_multi
,ginverse=list(Ring=ainv)
, na.method.X="include", na.method.Y ="include"
,maxiter=1000)
model.2<-asreml(fixed= HA~ 1+ Gcent#Grand-mean centered age,Model IXA
,random=~ us(pol(Gcent)):ide(Ring)
,data=adults.full_multi
,ginverse=list(Ring=ainv)
, na.method.X="include", na.method.Y ="include"
,maxiter=1000)
model.1b<-asreml(fixed=HA~ 1+ Icent #Individual centered age, Model I
,random=~ide(Ring)
,data=adults.full_multi
,ginverse=list(Ring=ainv)
, na.method.X="include", na.method.Y ="include"
,maxiter=1000)
model.2b<-asreml(fixed= HA~ 1+ Icent #Individual centered age, Model IXA
,random=~ us(pol(Icent)):ide(Ring)
,data=adults.full_multi
,ginverse=list(Ring=ainv)
, na.method.X="include", na.method.Y ="include"
,maxiter=1000)
K.estimates[z,1]<-summary(model.2)$varcomp[1,2]
K.estimates[z,2]<-summary(model.2)$varcomp[2,2]
K.estimates[z,3]<-summary(model.2)$varcomp[3,2]
K.estimates[z,4]<-model.2$conv
K.estimates[z,5]<-2*(model.2$loglik-model.1$loglik)
K.estimates[z,6]<-summary(model.2b)$varcomp[1,2]
K.estimates[z,7]<-summary(model.2b)$varcomp[2,2]
K.estimates[z,8]<-summary(model.2b)$varcomp[3,2]
K.estimates[z,9]<-model.2b$conv
K.estimates[z,10]<-2*(model.2b$loglik-model.1b$loglik)
}
write.table(K.estimates,"K_estimates_sim1000.txt")
################End of the simulation
#Calculate bias
K.estimates[,11]<-K.estimates[,1]-0.5
K.estimates[,12]<-K.estimates[,2]+0.3
K.estimates[,13]<-K.estimates[,3]-0.4
B1.1<-posterior.mode(as.mcmc((K.estimates[,11])))
B1.1.CI<-HPDinterval(as.mcmc((K.estimates[,11])),0.95)
B2.1<-posterior.mode(as.mcmc((K.estimates[,12])))
B2.1.CI<-HPDinterval(as.mcmc((K.estimates[,12])),0.95)
B3.1<-posterior.mode(as.mcmc((K.estimates[,13])))
B3.1.CI<-HPDinterval(as.mcmc((K.estimates[,13])),0.95)
K.estimates[,14]<-K.estimates[,6]-0.5
K.estimates[,15]<-K.estimates[,7]+0.3
K.estimates[,16]<-K.estimates[,8]-0.4
B1.2<-posterior.mode(as.mcmc((K.estimates[,14])))
B1.2.CI<-HPDinterval(as.mcmc((K.estimates[,14])),0.95)
B2.2<-posterior.mode(as.mcmc((K.estimates[,15])))
B2.2.CI<-HPDinterval(as.mcmc((K.estimates[,15])),0.95)
B3.2<-posterior.mode(as.mcmc((K.estimates[,16])))
B3.2.CI<-HPDinterval(as.mcmc((K.estimates[,16])),0.95)
#Plot the bias
point.estimates<-data.frame("points"=c(as.numeric(B1.1), as.numeric(B1.2),
as.numeric(B2.1), as.numeric(B2.2),
as.numeric(B3.1), as.numeric(B3.2)))
point.estimates$cor<-c("V1","V1","V2", "V2","V3", "V3")
CI.inf<-c(B1.1.CI[1],B1.2.CI[1],B2.1.CI[1],B2.2.CI[1],B3.1.CI[1],B3.2.CI[1] )
CI.sup<-c(B1.1.CI[2],B1.2.CI[2],B2.1.CI[2],B2.2.CI[2],B3.1.CI[2],B3.2.CI[2] )
plot(point.estimates$points,ylim=c(-0.5, 0.6),cex=1.5,cex.axis=1,ylab="Variance",
xlab=" ",
cex.lab=1.15,pch=c(1,2,1,2,1,2),
xaxt= "n")
axis(1, at =c(1.5,3.5,5.5),labels=c('V1', "V2", "V3"))
x<-c(1,2,3,4,5,6)
arrows(x, CI.inf, x, CI.sup, angle=90, code=3, length=0.1)
abline(0,0)
breaks<-c(2.5,4.5)
arrows(breaks, -0.8, breaks, 0.8, length=0,lty=2)