#Run the RRAM first (model.LOCOM1) #Calculate age-specific variances (age 1-3 divided into 12 classes here) using pin_RR function varspec<-data.frame("variance"=NA,"se.variance"=NA,"lower.ci"=NA,"upper.ci"=NA,"age"=NA) for (i in 1:12){ VA.i<-paste("VA",as.character(i),sep=".") transform=as.formula(paste("V~",VA.i,sep="")) variance<-pin_RR(model.LOCOM1,transform,12)[1] se.variance<-pin_RR(model.LOCOM1,transform,12)[2] lower.ci<-variance-(1.96*se.variance) upper.ci<-variance+(1.96*se.variance) varspec[i,1:4]<-cbind(variance,se.variance,lower.ci,upper.ci) } varspec$age<-seq(0,3,length.out = 12) res.var=summary(model.LOCOM1)$varcomp[c(9:12),2] res.lower=res.var -1.96*summary(model.LOCOM1)$varcomp[c(9:12),3] res.upper=res.var +1.96*summary(model.LOCOM1)$varcomp[c(9:12),3] #Estimate age-specific heritability using pin_RR function h2<-c(as.numeric(pin_RR(model.LOCOM1,h2~VA.1/(VA.1+VR.1),4)[1]) ,as.numeric(pin_RR(model.LOCOM1,h2~VA.2/(VA.2+VR.2),4)[1]) ,as.numeric(pin_RR(model.LOCOM1,h2~VA.3/(VA.3+VR.3),4)[1]) ,as.numeric(pin_RR(model.LOCOM1,h2~VA.4/(VA.4+VR.4),4)[1])) h2_se<-c(as.numeric(pin_RR(model.LOCOM1,h2~VA.1/(VA.1+VR.1),4)[2]) ,as.numeric(pin_RR(model.LOCOM1,h2~VA.2/(VA.2+VR.2),4)[2]) ,as.numeric(pin_RR(model.LOCOM1,h2~VA.3/(VA.3+VR.3),4)[2]) ,as.numeric(pin_RR(model.LOCOM1,h2~VA.4/(VA.4+VR.4),4)[2])) #estimate correlations across ages using pin_RR function Correlations<-data.frame("Param"=NA,"Cor"=NA,"SE"=NA) for( i in 2:12){ C.1.i<-paste("C.1",as.character(i),sep=".") VA.i<-paste("VA",as.character(i),sep=".") transform=as.formula(paste("rA~",C.1.i,"/sqrt(VA.1*",VA.i,")",sep="")) Correlations[i,1]<-C.1.i Correlations[i,2]<-as.numeric(pin_RR(model.LOCOM1,transform, 12)[1]) Correlations[i,3]<-as.numeric(pin_RR(model.LOCOM1,transform, 12)[2]) } Correlations<-Correlations[-1,] Correlations$lower.ci<-Correlations[,2]-1.96*Correlations[,3] Correlations$upper.ci<-Correlations[,2]+1.96*Correlations[,3] Correlations$age<-seq(0,3,length.out = 12)[-1] #Plot the estimates par(mfrow=c(4,1),oma=c(1,1,1,1),mgp=c(4.5,1,0),mar=c(2,6,2,6)) layout(matrix(c(1,1,2,3), nrow = 4, ncol = 1, byrow = TRUE)) plot(varspec$age,varspec$variance,type="l",ylim=c(-33,87),ylab="Additive genetic variance", xlab="Age",xaxt="n",lwd=2,cex.axis=2,cex.lab=2,las=1) par(new=TRUE) plot(varspec$age,varspec$lower.ci,type="l",ylim=c(-33,87),ylab=" ", xlab=" ",xaxt="n", yaxt="n",lty=2,lwd=2,cex.axis=2,cex.lab=2) par(new=TRUE) plot(varspec$age,varspec$upper.ci,type="l",ylim=c(-33,87),ylab=" ", xlab=" ",xaxt="n", yaxt="n",lty=2,lwd=2,cex.axis=2,cex.lab=2) axis(side=1,at=c(0,1,2,3),labels=c("0","1","2","3+"),cex.axis=2,cex.lab=2) abline(0,0,lty=3) title("A", adj = 0.01, line = -1.5,cex.main=2) par(new=T) plot(c(0,1,2,3),res.var,ylim=c(-20,110),ylab=" ", xlab=" ",xaxt="n",pch=16,yaxt="n",cex=2) arrows(c(0,1,2,3),res.lower,c(0,1,2,3),res.upper,angle=90,length = 0.1,code=3,lwd=2) axis(side=4,at=c(0,20,40,60,80,100),cex.axis=2,las=1) axis(side=4,at=c(40),cex.lab=2,cex.axis=2,labels=c("Residual variance"),line=3.1,tick=FALSE) plot(c(0,1,2,3),h2,ylim=c(-0.45,1.2),ylab="Heritability", xlab="Age",xaxt="n",pch=16,cex=2,cex.axis=2,cex.lab=2,las=1) axis(side=1,at=c(0,1,2,3),labels=c("0-1","1-2","2-3","3+"),cex.axis=2,cex.lab=2) arrows(c(0,1,2,3),h2-1.96*h2_se,c(0,1,2,3),h2+1.96*h2_se,angle=90,length = 0.1,code=3,lwd=2) abline(0,0,lty=3) title("B", adj = 0.01, line = -1.5,cex.main=2) plot(Correlations$Cor~Correlations$age,type="l",ylim=c(-1.3,1.4),ylab="Correlation", xlab=" ",xaxt="n",lwd=2,cex.axis=2,cex.lab=2,las=1) par(new=TRUE) plot(Correlations$lower.ci~Correlations$age,type="l",ylim=c(-1.3,1.4),ylab=" ", xlab=" ",xaxt="n", yaxt="n",lty=2,lwd=2,cex.axis=2,cex.lab=2) par(new=TRUE) plot(Correlations$upper.ci~Correlations$age,type="l",ylim=c(-1.3,1.4),ylab=" ", xlab=" ",xaxt="n", yaxt="n",lty=2,lwd=2,cex.axis=2,cex.lab=2) abline(0,0,lty=3) axis(side=1,at=c(0,1,2,3),labels=c("0","0-1","0-2","0-3+"),cex.axis=2,cex.lab=2) title("C", adj = 0.01, line = -1.5,cex.main=2)