Post date: Jun 25, 2015 5:21:32 PM
I ran 500,000 simulations and then estimated the parameters with the local linear regression correction with a tolerance of 0.002. I also ran this with a tolerance of 0.001 which gave similar but stronger results. The take home is that the migration rate estimates are consistent with reduced gene flow in parapatry at selection candidate and stripe regions, but not in allopatry, but the credible intervals are huge (so the results must be interpreted cautiously).
The results are in /labs/evolution/projects/timema_speciation/demog/sims/NmEstComb*csv and TEstComb*csv. The version with tolerance of 0.001 has 'Tol001' in the name.
## R code to summarize abc results, in /labs/evolution/projects/timema_speciation/demog/sims/summarizePost.R
library(abc)
## pops by meadian, 95CI and quartiles
pops<-c("HVAxHVC","MR1AxMR1C","R12AxR12C","LAxPRC")
sums<-c("ran_median","ran_se","stripe_median","stripe_se","pr(str<ran)","scan_median","scan_se","pr(scan<ran)")
npar<-10
nsim<-500000
lgb<-matrix(c(rep(0,7),rep(0.2,7)),nrow=7,ncol=2,byrow=FALSE)
nms<-c("N0","N1","T","M","Navg","Nm","T/N","fst","priv0","priv1")
NmEstComb<-matrix(NA,nrow=4,ncol=8)
TEstComb<-matrix(NA,nrow=4,ncol=8)
rownames(NmEstComb)<-pops
colnames(NmEstComb)<-sums
rownames(TEstComb)<-pops
colnames(TEstComb)<-sums
## random
d<-matrix(scan("abcSimCat_pRanHVAxHVC.txt",n=npar*nsim,sep=" "),nrow=nsim,ncol=npar,byrow=TRUE)
colnames(d)<-nms
obs<-read.table("../obsSS_pRanHVAxHVC.txt",header=FALSE)
RanHvaHvc<-abc(targe=obs,param=d[,1:7],sumstat=d[,8:10],method="loclinear",transf=c("log","log","log","logit","log","log","log"),logit.bounds=lgb,tol=0.002)
NmEstComb[1,1]<-median(RanHvaHvc$adj.values[,6])
NmEstComb[1,2]<-sd(RanHvaHvc$adj.values[,6])
TEstComb[1,1]<-median(RanHvaHvc$adj.values[,3])
TEstComb[1,2]<-sd(RanHvaHvc$adj.values[,3])
d<-matrix(scan("abcSimCat_pRanMR1AxMR1C.txt",n=npar*nsim,sep=" "),nrow=nsim,ncol=npar,byrow=TRUE)
colnames(d)<-nms
obs<-read.table("../obsSS_pRanMR1AxMR1C.txt",header=FALSE)
RanMr1aMr1c<-abc(targe=obs,param=d[,1:7],sumstat=d[,8:10],method="loclinear",transf=c("log","log","log","logit","log","log","log"),logit.bounds=lgb,tol=0.002)
NmEstComb[2,1]<-median(RanMr1aMr1c$adj.values[,6])
NmEstComb[2,2]<-sd(RanMr1aMr1c$adj.values[,6])
TEstComb[2,1]<-median(RanMr1aMr1c$adj.values[,3])
TEstComb[2,2]<-sd(RanMr1aMr1c$adj.values[,3])
d<-matrix(scan("abcSimCat_pRanR12AxR12C.txt",n=npar*nsim,sep=" "),nrow=nsim,ncol=npar,byrow=TRUE)
colnames(d)<-nms
obs<-read.table("../obsSS_pRanR12AxR12C.txt",header=FALSE)
RanR12aR12c<-abc(targe=obs,param=d[,1:7],sumstat=d[,8:10],method="loclinear",transf=c("log","log","log","logit","log","log","log"),logit.bounds=lgb,tol=0.002)
NmEstComb[3,1]<-median(RanR12aR12c$adj.values[,6])
NmEstComb[3,2]<-sd(RanR12aR12c$adj.values[,6])
TEstComb[3,1]<-median(RanR12aR12c$adj.values[,3])
TEstComb[3,2]<-sd(RanR12aR12c$adj.values[,3])
d<-matrix(scan("abcSimCat_pRanLAxPRC.txt",n=npar*nsim,sep=" "),nrow=nsim,ncol=npar,byrow=TRUE)
colnames(d)<-nms
obs<-read.table("../obsSS_pRanLAxPRC.txt",header=FALSE)
RanLaPrc<-abc(targe=obs,param=d[,1:7],sumstat=d[,8:10],method="loclinear",transf=c("log","log","log","logit","log","log","log"),logit.bounds=lgb,tol=0.002)
NmEstComb[4,1]<-median(RanLaPrc$adj.values[,6])
NmEstComb[4,2]<-sd(RanLaPrc$adj.values[,6])
TEstComb[4,1]<-median(RanLaPrc$adj.values[,3])
TEstComb[4,2]<-sd(RanLaPrc$adj.values[,3])
## stripe
d<-matrix(scan("abcSimCat_pStrHVAxHVC.txt",n=npar*nsim,sep=" "),nrow=nsim,ncol=npar,byrow=TRUE)
colnames(d)<-nms
obs<-read.table("../obsSS_pStrHVAxHVC.txt",header=FALSE)
StrHvaHvc<-abc(targe=obs,param=d[,1:7],sumstat=d[,8:10],method="loclinear",transf=c("log","log","log","logit","log","log","log"),logit.bounds=lgb,tol=0.002)
NmEstComb[1,3]<-median(StrHvaHvc$adj.values[,6])
NmEstComb[1,4]<-sd(StrHvaHvc$adj.values[,6])
NmEstComb[1,5]<-mean(StrHvaHvc$adj.values[,6] < RanHvaHvc$adj.values[,6])
TEstComb[1,3]<-median(StrHvaHvc$adj.values[,3])
TEstComb[1,4]<-sd(StrHvaHvc$adj.values[,3])
d<-matrix(scan("abcSimCat_pStrMR1AxMR1C.txt",n=npar*nsim,sep=" "),nrow=nsim,ncol=npar,byrow=TRUE)
colnames(d)<-nms
obs<-read.table("../obsSS_pStrMR1AxMR1C.txt",header=FALSE)
StrMr1aMr1c<-abc(targe=obs,param=d[,1:7],sumstat=d[,8:10],method="loclinear",transf=c("log","log","log","logit","log","log","log"),logit.bounds=lgb,tol=0.002)
NmEstComb[2,3]<-median(StrMr1aMr1c$adj.values[,6])
NmEstComb[2,4]<-sd(StrMr1aMr1c$adj.values[,6])
NmEstComb[2,5]<-mean(StrMr1aMr1c$adj.values[,6] < RanMr1aMr1c$adj.values[,6])
TEstComb[2,3]<-median(StrMr1aMr1c$adj.values[,3])
TEstComb[2,4]<-sd(StrMr1aMr1c$adj.values[,3])
d<-matrix(scan("abcSimCat_pStrR12AxR12C.txt",n=npar*nsim,sep=" "),nrow=nsim,ncol=npar,byrow=TRUE)
colnames(d)<-nms
obs<-read.table("../obsSS_pStrR12AxR12C.txt",header=FALSE)
StrR12aR12c<-abc(targe=obs,param=d[,1:7],sumstat=d[,8:10],method="loclinear",transf=c("log","log","log","logit","log","log","log"),logit.bounds=lgb,tol=0.002)
NmEstComb[3,3]<-median(StrR12aR12c$adj.values[,6])
NmEstComb[3,4]<-sd(StrR12aR12c$adj.values[,6])
NmEstComb[3,5]<-mean(StrR12aR12c$adj.values[,6] < RanR12aR12c$adj.values[,6])
TEstComb[3,3]<-median(StrR12aR12c$adj.values[,3])
TEstComb[3,4]<-sd(StrR12aR12c$adj.values[,3])
d<-matrix(scan("abcSimCat_pStrLAxPRC.txt",n=npar*nsim,sep=" "),nrow=nsim,ncol=npar,byrow=TRUE)
colnames(d)<-nms
obs<-read.table("../obsSS_pStrLAxPRC.txt",header=FALSE)
StrLaPrc<-abc(targe=obs,param=d[,1:7],sumstat=d[,8:10],method="loclinear",transf=c("log","log","log","logit","log","log","log"),logit.bounds=lgb,tol=0.002)
NmEstComb[4,3]<-median(StrLaPrc$adj.values[,6])
NmEstComb[4,4]<-sd(StrLaPrc$adj.values[,6])
NmEstComb[4,5]<-mean(StrLaPrc$adj.values[,6] < RanLaPrc$adj.values[,6])
TEstComb[4,3]<-median(StrLaPrc$adj.values[,3])
TEstComb[4,4]<-sd(StrLaPrc$adj.values[,3])
## selection candidates
d<-matrix(scan("abcSimCat_pScanHVAxHVC.txt",n=npar*nsim,sep=" "),nrow=nsim,ncol=npar,byrow=TRUE)
colnames(d)<-nms
obs<-read.table("../obsSS_pScanHVAxHVC.txt",header=FALSE)
ScanHvaHvc<-abc(targe=obs,param=d[,1:7],sumstat=d[,8:10],method="loclinear",transf=c("log","log","log","logit","log","log","log"),logit.bounds=lgb,tol=0.002)
NmEstComb[1,6]<-median(ScanHvaHvc$adj.values[,6])
NmEstComb[1,7]<-sd(ScanHvaHvc$adj.values[,6])
NmEstComb[1,8]<-mean(ScanHvaHvc$adj.values[,6] < RanHvaHvc$adj.values[,6])
TEstComb[1,6]<-median(ScanHvaHvc$adj.values[,3])
TEstComb[1,7]<-sd(ScanHvaHvc$adj.values[,3])
d<-matrix(scan("abcSimCat_pScanMR1AxMR1C.txt",n=npar*nsim,sep=" "),nrow=nsim,ncol=npar,byrow=TRUE)
colnames(d)<-nms
obs<-read.table("../obsSS_pScanMR1AxMR1C.txt",header=FALSE)
ScanMr1aMr1c<-abc(targe=obs,param=d[,1:7],sumstat=d[,8:10],method="loclinear",transf=c("log","log","log","logit","log","log","log"),logit.bounds=lgb,tol=0.002)
NmEstComb[2,6]<-median(ScanMr1aMr1c$adj.values[,6])
NmEstComb[2,7]<-sd(ScanMr1aMr1c$adj.values[,6])
NmEstComb[2,8]<-mean(ScanMr1aMr1c$adj.values[,6] < RanMr1aMr1c$adj.values[,6])
TEstComb[2,6]<-median(ScanMr1aMr1c$adj.values[,3])
TEstComb[2,7]<-sd(ScanMr1aMr1c$adj.values[,3])
d<-matrix(scan("abcSimCat_pScanR12AxR12C.txt",n=npar*nsim,sep=" "),nrow=nsim,ncol=npar,byrow=TRUE)
colnames(d)<-nms
obs<-read.table("../obsSS_pScanR12AxR12C.txt",header=FALSE)
ScanR12aR12c<-abc(targe=obs,param=d[,1:7],sumstat=d[,8:10],method="loclinear",transf=c("log","log","log","logit","log","log","log"),logit.bounds=lgb,tol=0.002)
NmEstComb[3,6]<-median(ScanR12aR12c$adj.values[,6])
NmEstComb[3,7]<-sd(ScanR12aR12c$adj.values[,6])
NmEstComb[3,8]<-mean(ScanR12aR12c$adj.values[,6] < RanR12aR12c$adj.values[,6])
TEstComb[3,6]<-median(ScanR12aR12c$adj.values[,3])
TEstComb[3,7]<-sd(ScanR12aR12c$adj.values[,3])
d<-matrix(scan("abcSimCat_pScanLAxPRC.txt",n=npar*nsim,sep=" "),nrow=nsim,ncol=npar,byrow=TRUE)
colnames(d)<-nms
obs<-read.table("../obsSS_pScanLAxPRC.txt",header=FALSE)
ScanLaPrc<-abc(targe=obs,param=d[,1:7],sumstat=d[,8:10],method="loclinear",transf=c("log","log","log","logit","log","log","log"),logit.bounds=lgb,tol=0.002)
NmEstComb[4,6]<-median(ScanLaPrc$adj.values[,6])
NmEstComb[4,7]<-sd(ScanLaPrc$adj.values[,6])
NmEstComb[4,8]<-mean(ScanLaPrc$adj.values[,6] < RanLaPrc$adj.values[,6])
TEstComb[4,6]<-median(ScanLaPrc$adj.values[,3])
TEstComb[4,7]<-sd(ScanLaPrc$adj.values[,3])
write.csv(NmEstComb,"NmEstComb.csv")
write.csv(TEstComb,"TEstComb.csv")
## plots of posterior densities
myc<-c("red","forestgreen","darkgray")
pdf("NmHVAxHVC.pdf",width=6,height=6)
par(mar=c(5,5,0.5,0.5))
plot(density(ScanHvaHvc$adj.values[,6],from=0,to=500),xlab="number of migrants",ylab="posterior density",cex.lab=1.4,col=myc[1],main="",xlim=c(0,500),lwd=1.5)
lines(density(strHvaHvc$adj.values[,6],from=0,to=500),col=myc[2],lwd=1.5)
lines(density(RanHvaHvc$adj.values[,6],from=0,to=500),col=myc[3],lwd=1.5)
legend(275,0.05,c("neutral","stripe","selection candidate"),lty=1,lwd=1.5,col=myc[c(3,2,1)],bty='n')
dev.off()
pdf("NmMR1AxMR1C.pdf",width=6,height=6)
par(mar=c(5,5,0.5,0.5))
plot(density(ScanMr1aMr1c$adj.values[,6],from=0,to=500),xlab="number of migrants",ylab="posterior density",cex.lab=1.4,col=myc[1],main="",ylim=c(0,0.05),lwd=1.5)
lines(density(strMr1aMr1c$adj.values[,6],from=0,to=500),col=myc[2],lwd=1.5)
lines(density(RanMr1aMr1c$adj.values[,6],from=0,to=500),col=myc[3],lwd=1.5)
legend(275,0.05,c("neutral","stripe","selection candidate"),lty=1,lwd=1.5,col=myc[c(3,2,1)],bty='n')
dev.off()
pdf("NmR12AxR12C.pdf",width=6,height=6)
par(mar=c(5,5,0.5,0.5))
plot(density(ScanR12aR12c$adj.values[,6],from=0,to=500),xlab="number of migrants",ylab="posterior density",cex.lab=1.4,col=myc[1],main="",xlim=c(0,500),lwd=1.5)
lines(density(strR12aR12c$adj.values[,6],from=0,to=500),col=myc[2],lwd=1.5)
lines(density(RanR12aR12c$adj.values[,6],from=0,to=500),col=myc[3],lwd=1.5)
legend(275,0.027,c("neutral","stripe","selection candidate"),lty=1,lwd=1.5,col=myc[c(3,2,1)],bty='n')
dev.off()
pdf("NmLAxPRC.pdf",width=6,height=6)
par(mar=c(5,5,0.5,0.5))
plot(density(ScanLaPrc$adj.values[,6],from=0,to=500),xlab="number of migrants",ylab="posterior density",cex.lab=1.4,col=myc[1],main="",xlim=c(0,500),lwd=1.5,ylim=c(0,0.08))
lines(density(strLaPrc$adj.values[,6],from=0,to=500),col=myc[2],lwd=1.5)
lines(density(RanLaPrc$adj.values[,6],from=0,to=500),col=myc[3],lwd=1.5)
legend(275,0.075,c("neutral","stripe","selection candidate"),lty=1,lwd=1.5,col=myc[c(3,2,1)],bty='n')
dev.off()