Post date: Dec 05, 2015 3:2:4 AM
We decided to drop 2011 and combine 1996 and 2001 in a single analysis. We then fit clines based just on stripe vs. green or assuming intermediate individuals were equivalent to stripe or green individuals. We dropped 2011 because Patrik did this alone (the earlier data were collected with Cris Sandoval) and he and Cris scored intermediates differently. Clines with the combined data were steep. Here is the relevant part of the code from the clines.R script.
####### 1996 + 2001 ###########
# read data
hdat<-read.csv("refugio_by_year.csv")
## stripe allele frequency
p<-(hdat$striped_2001+hdat$striped_1996)/(hdat$green_2001+hdat$striped_2001+hdat$green_1996+hdat$striped_1996)
p<-sqrt(p)
## create data object
hobj<-hzar.doMolecularData1DPops(distance=hdat$distance_from_site1,pObs=p,nEff=2*(hdat$green_2001+hdat$striped_2001+hdat$green_1996+hdat$striped_1996))
## creat model
hmod<-hzar.makeCline1DFreq(data=hobj, scaling="none", tails="both",direction=NULL)
## generate and fit object
hfitC1<- hzar.first.fitRequest.old.ML(hmod, hobj, verbose = TRUE)
hfitC1$mcmcParam$chainLength<-2000000
hfitC1$mcmcParam$burnin<-1000000
houtC1<-hzar.doFit(hfitC1)
hfitC2<- hzar.first.fitRequest.old.ML(hmod, hobj, verbose = TRUE)
hfitC2$mcmcParam$chainLength<-2000000
hfitC2$mcmcParam$burnin<-1000000
houtC2<-hzar.doFit(hfitC2)
save(list=ls(),file="cline9601.Rdata")
## stripe (int.= stripe) allele frequency
p<-(hdat$striped_2001+hdat$intermediate_2001+hdat$striped_1996+hdat$intermediate_1996)/(hdat$green_2001+hdat$striped_2001+hdat$intermediate_2001+hdat$green_1996+hdat$striped_1996+hdat$intermediate_1996)
p<-sqrt(p)
## create data object
hobj<-hzar.doMolecularData1DPops(distance=hdat$distance_from_site1,pObs=p,nEff=2*(hdat$green_2001+hdat$striped_2001+hdat$intermediate_2001+hdat$green_1996+hdat$striped_1996+hdat$intermediate_1996))
## creat model
hmod<-hzar.makeCline1DFreq(data=hobj, scaling="none", tails="both",direction=NULL)
## generate and fit object
hfitC1<- hzar.first.fitRequest.old.ML(hmod, hobj, verbose = TRUE)
hfitC1$mcmcParam$chainLength<-2000000
hfitC1$mcmcParam$burnin<-1000000
houtC1<-hzar.doFit(hfitC1)
hfitC2<- hzar.first.fitRequest.old.ML(hmod, hobj, verbose = TRUE)
hfitC2$mcmcParam$chainLength<-2000000
hfitC2$mcmcParam$burnin<-1000000
houtC2<-hzar.doFit(hfitC2)
save(list=ls(),file="cline9601IntEqStr.Rdata")
## stripe (int.= green) allele frequency
p<-(hdat$striped_2001+hdat$striped_1996)/(hdat$green_2001+hdat$striped_2001+hdat$intermediate_2001+hdat$green_1996+hdat$striped_1996+hdat$intermediate_1996)
p<-sqrt(p)
## create data object
hobj<-hzar.doMolecularData1DPops(distance=hdat$distance_from_site1,pObs=p,nEff=2*(hdat$green_2001+hdat$striped_2001+hdat$intermediate_2001+hdat$green_1996+hdat$striped_1996+hdat$intermediate_1996))
## creat model
hmod<-hzar.makeCline1DFreq(data=hobj, scaling="none", tails="both",direction=NULL)
## generate and fit object
hfitC1<- hzar.first.fitRequest.old.ML(hmod, hobj, verbose = TRUE)
hfitC1$mcmcParam$chainLength<-2000000
hfitC1$mcmcParam$burnin<-1000000
houtC1<-hzar.doFit(hfitC1)
hfitC2<- hzar.first.fitRequest.old.ML(hmod, hobj, verbose = TRUE)
hfitC2$mcmcParam$chainLength<-2000000
hfitC2$mcmcParam$burnin<-1000000
houtC2<-hzar.doFit(hfitC2)
save(list=ls(),file="cline9601IntEqGrn.Rdata")
########################################
ho<-vector("list",3)
hp<-vector("list",3)
load("cline9601.Rdata")
ho[[1]]<-hzar.dataGroup.add(houtC1)
hp[[1]]<-p
load("cline9601IntEqStr.Rdata")
ho[[2]]<-hzar.dataGroup.add(houtC1)
hp[[2]]<-p
load("cline9601IntEqGrn.Rdata")
ho[[3]]<-hzar.dataGroup.add(houtC1)
hp[[3]]<-p
## 3-panel combined plot
hzar.plot.fzCline<-function (dataGroup, fzCline = hzar.getCredParamRed(dataGroup),
type = "p", pch = "+", col = "black", fzCol = "gray", ...) {
hzar.plot.obsData(dataGroup, col = "transparent", ...)
xSeries <- seq(from = par("usr")[1], to = par("usr")[2],
length.out = 109)
if (par("xaxs") == "r")
xSeries <- xSeries[2:108]
fzCoor <- fzCline$fzCline(xSeries)
polygon(x = c(fzCoor$x, rev(fzCoor$x)), y = c(fzCoor$yMin,
rev(fzCoor$yMax)), border = fzCol, col = fzCol)
lines(x = xSeries, y = dataGroup$ML.cline$clineFunc(xSeries),
col = "black")
#hzar.plot.obsData(dataGroup, col = col, type = type, pch = pch,
# add = TRUE)
}
xc<-rep(c(0,0,0),1)
tit<-c("(a) striped only","(b) int.=stripe","(c) int.=green")
pdf("transectPlots3p.pdf",width=4,height=12)
par(mfrow=c(3,1))
par(mar=c(4,5,2,1))
for(i in 1:3){
hzar.plot.fzCline(ho[[i]],fzCline=hzar.getCredParamRed(ho[[i]]),xlab="transect location (km)",ylab="stripe frequency",cex.lab=1.5,cex.axis=1.1,pch=20,col=c("blue","orange")[hdat$host])
points(hdat$distance_from_site1,hp[[i]],pch=20,col=c("blue","orange")[hdat$host])
center<-quantile(ho[[i]]$data.mcmc[,1],probs=c(0.5,0.025,0.975))
width<-quantile(ho[[i]]$data.mcmc[,2],probs=c(0.5,0.025,0.975))
text(xc[i],1,paste("center = ",round(center[1],2)," (",round(center[2],2),"-",round(center[3],2),")",sep=""),adj=0)
text(xc[i],.925,paste("width = ",round(width[1],2)," (",round(width[2],2),"-",round(width[3],2),")",sep=""),adj=0)
title(main=tit[i],adj=0)
}
dev.off()