Post date: Aug 30, 2013 7:59:28 PM
I think bgsr is working properly based on my analysis of simulated data. Most importantly, true and estimated values of allele frequency change were highly correlated (~0.9) regardless of the selection model. Ne is higher than I might expect (~90, but my simulations didn't really include a Ne). In some ways Ne will only be meaningful when lower than N in the next generation. There is no evidence of selection (95% CI's for s overlap zero) with assumed sd on betas of 0.01 or 0.05 (the former yields very low values of s), but point estimates of s are correlated with allele frequency change (r ~0.49). But estimates of s are also higher for loci with intermediate allele frequencies and this is particularly true of loci with highest s (i.e., the highest s are for loci with intermediate allele frequencies, not those that were actually selected). This is important, suggests that only very strong selection can be reliably detected (note, selection on rare variants is of course weaker at the population-level, so this is kind of ok), and that caution and restraint will be needed when interpreting the analysis of the real data... in other words, we need to pay attention to the credible intervals.
Here is the main R code I used to summarize the output:
load("~/Documents/programs/bgsr/sims/sim.rwsp")
## surviving allele freq
psurv<-matrix(NA,nrow=nloci,ncol=npops)
for(i in 1:nloci){
psurv[i,]<-tapply(INDEX=popsurv,X=gsurv[i,],mean)/2
}
dpfull05<-read.table("dpfull05",header=FALSE,sep=",")
dpnull01<-read.table("dpnull01",header=FALSE,sep=",")
dpcons05<-read.table("dpcons05",header=FALSE,sep=",")
dpfull01<-read.table("dpfull01",header=FALSE,sep=",")
dpcons01<-read.table("dpcons01",header=FALSE,sep=",")
dpobs<-matrix(NA,nrow=nloci,ncol=npops)
for(i in 1:nloci){
dpobs[i,]<-psurv[i,]-pi[i]
}
dpestFull05<-matrix(dpfull05[,2],nrow=nloci,ncol=npops,byrow=TRUE)
dpestNull01<-matrix(dpnull01[,2],nrow=nloci,ncol=npops,byrow=TRUE)
p1<-read.table("p1null",header=FALSE,sep=",")
p0<-read.table("p0null",header=FALSE,sep=",")
sfull05<-read.table("sfull05",header=FALSE,sep=",")
scons05<-read.table("scons05",header=FALSE,sep=",")
sfull01<-read.table("sfull01",header=FALSE,sep=",")
scons01<-read.table("scons01",header=FALSE,sep=",")