Post date: Dec 01, 2015 11:14:53 PM
Patrik has data on morph frequencies across a transect where host plants transition from C to A rather abruptly. I used hzar (a R package) to fit 6-paramter geogrphic clines for these data. For now I have only considered green and striped individuals and converted these to stripe allele frequencies based on the assumption that stripe is recessive and populations are in HWE (the latter is a bit dubious of course). He has data from three years: 1996, 2001, and 2011. I fit each year separately and all three in a single analysis. All analyses suggest a steep cline, though the width varies a bit (steeper for 1996 and 2011 on their own than all three combined; 1996 is missing a key data point so is not very relevant on its own). The R code, results, and infiles are all in gompert:/home/zgompert/Local/timemaCline/. And the R script (clines.R) is included below.
# read data
hdat<-read.csv("refugio_by_year.csv")
####### 2011 ############
## stripe allele frequency
p<-hdat$striped_2011/(hdat$green_2011+hdat$striped_2011)
p<-sqrt(p)
## create data object
hobj<-hzar.doMolecularData1DPops(distance=hdat$distance_from_site1,pObs=p,nEff=(2*hdat$green_2011+hdat$striped_2011))
## 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)
## extract parameters and make plots
pdf("transectPlot2011.pdf",width=6,height=6)
par(mar=c(5,5,0.5,0.5))
plot(hdat$distance_from_site1,p,pch=20,col=c("blue","orange")[hdat$host],xlab="transect location (km)",ylab="stripe frequency",cex.lab=1.4,cex.axis=1.1)
hcline<-hzar.get.ML.cline(fitRequest=houtC1)
hzar.plot.cline(cline=hcline,add=TRUE)
dev.off()
pdf("transectPlot2011CI.pdf",width=6,height=6)
par(mar=c(5,5,0.5,0.5))
hzar.plot.fzCline(houtC1,fzCline=hzar.getCredParamRed(houtC1),xlab="transect location (km)",ylab="stripe frequency",cex.lab=1.4,cex.axis=1.1)
points(hdat$distance_from_site1,p,pch=20,col=c("blue","orange")[hdat$host])
dev.off()
save(list=ls(),file="cline2011.Rdata")
####### 1996 ############
## stripe allele frequency
p<-hdat$striped_1996/(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_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)
## extract parameters and make plots
pdf("transectPlot1996.pdf",width=6,height=6)
par(mar=c(5,5,0.5,0.5))
plot(hdat$distance_from_site1,p,col=c("orange","blue")[hdat$host],xlab="transect location",ylab="stripe frequency",cex.lab=1.4,cex.axis=1.1)
hcline<-hzar.get.ML.cline(fitRequest=houtC1)
hzar.plot.cline(cline=hcline,add=TRUE)
dev.off()
save(list=ls(),file="cline1996.Rdata")
hzar.plot.fzCline(houtC1,fzCline=hzar.getCredParamRed(houtC1),add=TRUE)
####### 2001 ############
## stripe allele frequency
p<-hdat$striped_2001/(hdat$green_2001+hdat$striped_2001)
p<-sqrt(p)
## create data object
hobj<-hzar.doMolecularData1DPops(distance=hdat$distance_from_site1,pObs=p,nEff=(2*hdat$green_2001+hdat$striped_2001))
## 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)
## extract parameters and make plots
pdf("transectPlot2001.pdf",width=6,height=6)
par(mar=c(5,5,0.5,0.5))
plot(hdat$distance_from_site1,p,col=c("orange","blue")[hdat$host],xlab="transect location",ylab="stripe frequency",cex.lab=1.4,cex.axis=1.1)
hcline<-hzar.get.ML.cline(fitRequest=houtC1)
hzar.plot.cline(cline=hcline,add=TRUE)
dev.off()
save(list=ls(),file="cline2001.Rdata")
hzar.plot.fzCline(houtC1,fzCline=hzar.getCredParamRed(houtC1),add=TRUE)
####### Combined years ############
## stripe allele frequency
p<-(hdat$striped_2011+hdat$striped_2001+hdat$striped_1996)/(hdat$green_2011+hdat$striped_2011+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_2011+hdat$striped_2011+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="clineAll.Rdata")
#### MAKE FINAL PLOTS #######
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)
}
load("clineAll.Rdata")
pdf("transectPlotCombinedCI2.pdf",width=6,height=6)
par(mar=c(5,5,0.5,0.5))
hzar.plot.fzCline(hout,fzCline=hzar.getCredParamRed(hout),xlab="transect location (km)",ylab="stripe frequency",cex.lab=1.4,cex.axis=1.1)
center<-quantile(houtC2$mcmcRaw[,1],probs=c(0.5,0.025,0.975))
width<-quantile(houtC2$mcmcRaw[,2],probs=c(0.5,0.025,0.975))
text(0,1,paste("center = ",round(center[1],2)," (",round(center[2],2),"-",round(center[3],2),")",sep=""),adj=0)
text(0,.925,paste("width = ",round(width[1],2)," (",round(width[2],2),"-",round(width[3],2),")",sep=""),adj=0)
dev.off()
pdf("transectPlotCombinedCI.pdf",width=6,height=6)
par(mar=c(5,5,0.5,0.5))
hzar.plot.fzCline(hout,fzCline=hzar.getCredParamRed(hout),xlab="transect location (km)",ylab="stripe frequency",cex.lab=1.4,cex.axis=1.1,pch=20,col=c("blue","orange")[hdat$host])
center<-quantile(hout$data.mcmc[,1],probs=c(0.5,0.025,0.975))
width<-quantile(hout$data.mcmc[,2],probs=c(0.5,0.025,0.975))
text(0,1,paste("center = ",round(center[1],2)," (",round(center[2],2),"-",round(center[3],2),")",sep=""),adj=0)
text(0,.925,paste("width = ",round(width[1],2)," (",round(width[2],2),"-",round(width[3],2),")",sep=""),adj=0)
dev.off()
load("cline2011.Rdata")
pdf("transectPlot2011CI.pdf",width=6,height=6)
par(mar=c(5,5,0.5,0.5))
hzar.plot.fzCline(houtC1,fzCline=hzar.getCredParamRed(houtC1),xlab="transect location (km)",ylab="stripe frequency",cex.lab=1.4,cex.axis=1.1)
points(hdat$distance_from_site1,p,pch=20,col=c("blue","orange")[hdat$host])
center<-quantile(houtC1$mcmcRaw[,1],probs=c(0.5,0.025,0.975))
width<-quantile(houtC1$mcmcRaw[,2],probs=c(0.5,0.025,0.975))
text(0,1,paste("center = ",round(center[1],2)," (",round(center[2],2),"-",round(center[3],2),")",sep=""),adj=0)
text(0,.925,paste("width = ",round(width[1],2)," (",round(width[2],2),"-",round(width[3],2),")",sep=""),adj=0)
dev.off()
load("cline1996.Rdata")
pdf("transectPlot1996CI.pdf",width=6,height=6)
par(mar=c(5,5,0.5,0.5))
hzar.plot.fzCline(houtC1,fzCline=hzar.getCredParamRed(houtC1),xlab="transect location (km)",ylab="stripe frequency",cex.lab=1.4,cex.axis=1.1)
points(hdat$distance_from_site1,p,pch=20,col=c("blue","orange")[hdat$host])
center<-quantile(houtC1$mcmcRaw[,1],probs=c(0.5,0.025,0.975))
width<-quantile(houtC1$mcmcRaw[,2],probs=c(0.5,0.025,0.975))
text(0,1,paste("center = ",round(center[1],2)," (",round(center[2],2),"-",round(center[3],2),")",sep=""),adj=0)
text(0,.925,paste("width = ",round(width[1],2)," (",round(width[2],2),"-",round(width[3],2),")",sep=""),adj=0)
dev.off()
load("cline2001.Rdata")
pdf("transectPlot2001CI.pdf",width=6,height=6)
par(mar=c(5,5,0.5,0.5))
hzar.plot.fzCline(houtC1,fzCline=hzar.getCredParamRed(houtC1),xlab="transect location (km)",ylab="stripe frequency",cex.lab=1.4,cex.axis=1.1)
points(hdat$distance_from_site1,p,pch=20,col=c("blue","orange")[hdat$host])
center<-quantile(houtC1$mcmcRaw[,1],probs=c(0.5,0.025,0.975))
width<-quantile(houtC1$mcmcRaw[,2],probs=c(0.5,0.025,0.975))
text(2,1,paste("center = ",round(center[1],2)," (",round(center[2],2),"-",round(center[3],2),")",sep=""),adj=0)
text(2,.925,paste("width = ",round(width[1],2)," (",round(width[2],2),"-",round(width[3],2),")",sep=""),adj=0)
dev.off()