Post date: Feb 17, 2016 6:14:38 PM
Originally I had compared change over the past 10 (or 11) generations to change over the whole time period (recent change in reversion lines). However, a reviewer pointed out that it would be preferable to compare the last 10 generations to the whole time excluding the last 10 generations (such that the two intervals would be independent). I did this, and while it reduced correlations some, the overall results were very similar. The new code is in rev2.R, but the relevant bit (that is, the bit the changed) is here:
pdf("reversion.pdf",width=4,height=12)
par(mfrow=c(3,1))
par(mar=c(4,5.5,2,0.5))
cl<-1.6
ca<-1.3
cm<-1.3
cx<-1.3
cs<-c("gray","#7570B3")
sym<-c(20,21)
## L1
## L1R-F35 - L1-F100
x<-l1[,4]-l1[,5]
## L1R-F46 - L1R-F35
y<-l1[,3]-l1[,4]
a<-which((x*y) > 0)
indx<-rep(1,length(x))
indx[a]<-2
indx2<-rep(1,length(x))
indx2[sigl1]<-2
plot(x,y,cex=cx,cex.lab=cl,cex.axis=ca,xlab=expression(paste(Delta[p]," L1R-F35 - L1-F100")),ylab=expression(paste(Delta[p]," L1R-F46 - L1R-F35")),pch=sym[indx2],col=cs[indx])
abline(h=0,lty=2)
abline(v=0,lty=2)
out<-cor.test(x,y)
r<-round(out$estimate,3)
title(main=paste("(a) L1R-F46 vs. L1R-F35, r = ",r," p = 0.013"),adj=0,cex.main=cm)
legend(-0.8,0.27,c("same, p < 0.05","same, p > 0.05","diff. p < 0.05","diff. p > 0.05"),cex=cx,pch=rep(rev(sym),2),col=rep(rev(cs),each=2))
## L2
## L2R-F35 - L2-F87
x<-l2[,4]-l2[,5]
## L2R-F45 - L2R-F35
y<-l2[,3]-l2[,4]
a<-which((x*y) > 0)
indx<-rep(1,length(x))
indx[a]<-2
indx2<-rep(1,length(x))
indx2[sigl2]<-2
plot(x,y,cex=cx,cex.lab=cl,cex.axis=ca,xlab=expression(paste(Delta[p]," L2R-F35 - L2-F87")),ylab=expression(paste(Delta[p]," L2R-F45 - L2R-F35")),pch=sym[indx2],col=cs[indx])
abline(h=0,lty=2)
abline(v=0,lty=2)
out<-cor.test(x,y)
r<-round(out$estimate,3)
title(main=paste("(b) L2R-F45 vs. L2R-F35, r = ",r," p < 0.001"),adj=0,cex.main=cm)
## L3
## L3R-F35 - L3-F85
x<-l3[,4]-l3[,5]
## L3R-F45 - L3R-F35
y<-l3[,3]-l3[,4]
a<-which((x*y) > 0)
indx<-rep(1,length(x))
indx[a]<-2
indx2<-rep(1,length(x))
indx2[sigl3]<-2
plot(x,y,cex=cx,cex.lab=cl,cex.axis=ca,xlab=expression(paste(Delta[p]," L3R-F35 - L3-F85")),ylab=expression(paste(Delta[p]," L3R-F45 - L3R-F35")),pch=sym[indx2],col=cs[indx])
abline(h=0,lty=2)
abline(v=0,lty=2)
out<-cor.test(x,y)
r<-round(out$estimate,3)
title(main=paste("(c) L3R-F45 vs. L3R-F35, r = ",r," p = 0.006"),adj=0,cex.main=cm)
dev.off()