inst/doc/mixstock.R

In mixstock: functions for mixed stock analysis

## ----Rver,echo=FALSE-----------------------------------------------------
Rver <- paste(R.version$major,R.version$minor,sep=".")

## ----knitopts,echo=FALSE,message=FALSE-----------------------------------
library("knitr")
opts_chunk$set(fig.width=6,fig.height=4)
knit_hooks$set(base.figure=function() { par(las=1,bty="l") })

## ----instpkgs,eval=FALSE-------------------------------------------------
#  install.packages("mixstock")

## ----loadpkg-------------------------------------------------------------
library(mixstock)

## ----readdata1,eval=FALSE------------------------------------------------
#  mydata = read.table("mydata.dat",header=TRUE)

## ----readdata2,eval=FALSE------------------------------------------------
#  mydata = read.csv("mydata.csv")

## ----lahanas98raw--------------------------------------------------------
mydata = lahanas98raw

## ----printdata,eval=FALSE------------------------------------------------
#  mydata

## ----head----------------------------------------------------------------
head(mydata)

## ----convert1------------------------------------------------------------
mydata = as.mixstock.data(mydata)

## ----data1,echo=FALSE,fig.scap="Basic plot of turtle mtDNA haplotype data",fig.cap="Basic plot of turtle mtDNA haplotype data, using \\code{plot(mydata,mix.off=2)} (\\code{mix.off=2} leaves a slightly larger space between the rookery and mixed stock data)",base.figure=TRUE----
plot(mydata,mix.off=2)

## ----echo=FALSE,eval=FALSE-----------------------------------------------
#  require(ggplot2)
#  require(reshape2)
#  tmpf <- function(m0) {
#    mm1 <- as.matrix(m0)
#    mm2 <- sweep(mm1,2,colSums(mm1),"/")
#    transform(melt(mm2),
#              Var1=factor(Var1,levels=rownames(m0)))
#  }
#  m2 <- matrix(mydata$mixsamp,ncol=1,dimnames=list(names(mydata$mixsamp),"Mixed"))
#  g1 <- ggplot(tmpf(mydata$sourcesamp),aes(x=Var2,y=value,fill=Var1))+geom_bar()+coord_flip()+
#    scale_fill_brewer(palette="Set3")+theme_bw()+labs(x="")
#  g2 <- ggplot(tmpf(m2),aes(x=Var2,y=value,fill=Var1))+geom_bar()+coord_flip()+
#    scale_fill_brewer(palette="Set3")+theme_bw()+labs(x="")+opts(legend.position="none")
#  ## library(ggExtra)
#  ## align.plots(g1,g2,heights=unit(c(0.1,0.9)))
#  print()

## ----condense------------------------------------------------------------
mydata = markfreq.condense(mydata)

## ----condensed,echo=FALSE,fig.scap="Condensed haplotype data from Lahanas 1998",fig.cap="Condensed haplotype data from Lahanas 1998 (\\code{plot(lahanas98, mix.off=2, leg.space=0.4)}; \\code{leg.space=0.4} leaves more room for the legend)."----
plot(lahanas98,mix.off=2,leg.space=0.4)

## ----cml0----------------------------------------------------------------
mydata.cml = cml(mydata)
mydata.cml

## ----cmlplot1,eval=FALSE-------------------------------------------------
#  plot(mydata.cml)

## ----cml1,echo=FALSE,fig.scap="CML estimates for Lahanas 1998 data",fig.cap="CML estimates for Lahanas 1998 data; \\code{plot(mydata.cml)}",base.figure=TRUE,fig.width=8----
plot(mydata.cml)

## ----uml0----------------------------------------------------------------
mydata.uml = uml(mydata)

## ----eval=FALSE----------------------------------------------------------
#  par(ask=TRUE)
#  plot(mydata.uml,plot.freqs=TRUE)
#  par(ask=FALSE)

## ----eval=FALSE----------------------------------------------------------
#  mydata.umlboot = genboot(mydata,"uml")

## ----echo=FALSE----------------------------------------------------------
if (!file.exists("umlboot.RData")) {
  t1 <- system.time(mydata.umlboot <- genboot(mydata,"uml"))
  save("t1","mydata.umlboot",file="umlboot.RData",compress=TRUE)
} else load("umlboot.RData")

## ----confint1------------------------------------------------------------
confint(mydata.umlboot)

## ----umlboot,echo=FALSE,fig.scap="UML with bootstrap CI",fig.cap="UML estimates with bootstrap confidence limits for Lahanas 1998 data: \\code{plot(mydata.umlboot)}",base.figure=TRUE,fig.width=8----
plot(mydata.umlboot,ylim=c(0,1))

## ----runmcmc,eval=FALSE--------------------------------------------------
#  mydata.mcmc = tmcmc(mydata)

## ----getmcmcdata,echo=FALSE----------------------------------------------
if (!file.exists("mcmc.RData")) {
  mydata.mcmc <- tmcmc(mydata)
  save("mydata.mcmc",file="mcmc.RData",compress=TRUE)
} else load("mcmc.RData")

## ----confint2------------------------------------------------------------
mydata.mcmc
confint(mydata.mcmc)

## ----plot0,eval=FALSE----------------------------------------------------
#  plot(mydata.mcmc)

## ----plot_mcmc,echo=FALSE,base.figure=TRUE,message=FALSE,fig.width=8,base.figure=TRUE----
plot(mydata.mcmc,ylim=c(0,1),axes=FALSE)
axis(side=1,at=1:mydata.mcmc$R,labels=gsub("^contrib\\.","",
                               names(coef(mydata.mcmc)$input.freq)))
axis(side=2)

## ----eval=FALSE----------------------------------------------------------
#  diag1=calc.RL.0(mydata)

## ----echo=FALSE----------------------------------------------------------
if (!file.exists("RL.RData")) {
  diag1=calc.RL.0(mydata)
  save("diag1",file="RL.RData")
} else load("RL.RData")

## ------------------------------------------------------------------------
head(diag1$current)

## ------------------------------------------------------------------------
diag1$history

## ----eval=FALSE----------------------------------------------------------
#  diag2=calc.GR(mydata)

## ----echo=FALSE----------------------------------------------------------
if (!file.exists("GR.RData")) {
  diag2=calc.GR(mydata)
  save("diag2",file="GR.RData")
} else load("GR.RData")

## ----plot_sim1,fig.keep="none"-------------------------------------------
Z = simmixstock2(nsource=4,nmark=5,nmix=3,
                 sourcesize=c(4,2,1,1),
                 sourcesampsize=rep(25,4),
                 mixsampsize=rep(30,3),rseed=1001)               
Z
plot(Z)

## ----mmplot0,fig.cap="Simulated multiple-mixed-stock data",base.figure=TRUE,echo=FALSE,fig.scap="multiple-mixed-stock simulation"----
plot(Z)

## ----echo=FALSE,cache=TRUE-----------------------------------------------
if (!file.exists("wbugs_cache.RData")) {
  ## FIXME: would be more efficient to run Zfit0 and then convert
  ##  the results rather than fitting the whole model twice!
  st0 <- system.time(Zfit0 <- mm.wbugs(Z,sourcesize=c(4,2,1,1),
                                       returntype="bugs",pkg="JAGS"))
  st1 <- system.time(Zfit <- mm.wbugs(Z,sourcesize=c(4,2,1,1),
                                      pkg="JAGS"))
  st2 <- system.time(Zfit2 <- mm.wbugs(Z,sourcesize=c(4,2,1,1),
                                       bugs.code="BB",pkg="JAGS"))
  save("st0","st1","st2","Zfit","Zfit2","Zfit0",file="wbugs_cache.RData")
} else load("wbugs_cache.RData")

## ----eval=FALSE,echo=FALSE-----------------------------------------------
#  ## tests
#  st3 <- system.time(Zfit0 <- mm.wbugs(Z,sourcesize=c(4,2,1,1),
#                                       bugs.code="BB",pkg="JAGS"))
#  Zfit2 <- mm.wbugs(Z,sourcesize=c(4,2,1,1),bugs.code="BB",n.iter=200)
#  Zfit3 <- mm.wbugs(Z,sourcesize=c(4,2,1,1),bugs.code="BB",n.iter=1000,pkg="JAGS")

## ----bugsfit2,eval=FALSE-------------------------------------------------
#  Zfit0 = mm.wbugs(Z,sourcesize=c(4,2,1,1),returntype="bugs")

## ----eval=FALSE----------------------------------------------------------
#  plot(Zfit0)

## ----plot_coda,fig.keep="none",eval=FALSE--------------------------------
#  plot(as.mcmc.bugs(Zfit0))

## ----plot_res1-----------------------------------------------------------
plot(Zfit)

## ----plot_res2-----------------------------------------------------------
plot(Zfit,sourcectr=TRUE)

## ----sumZfit-------------------------------------------------------------
summary(Zfit)

## ----plot_g2,message=FALSE,tidy=FALSE,fig.height=3-----------------------
Zfit.d <- as.data.frame(Zfit)
library(ggplot2)
library(grid)
theme_set(theme_bw())
zmargin <- theme(panel.margin=unit(0,"lines"))
ggplot(subset(Zfit.d,type=="input.freq"),
       aes(x=from,y=est,ymin=lwr,ymax=upr))+
    geom_pointrange()+facet_wrap(~to)+zmargin

## ----quick_install,eval=FALSE--------------------------------------------
#  install.packages("mixstock")

## ----quick_condense,eval=FALSE-------------------------------------------
#  mydata = hapfreq.condense(as.mixstock.data(read.csv("myfile.dat")))

## ----quick_analyze,eval=FALSE--------------------------------------------
#  mydata.mcmc = tmcmc(mydata)
#  mydata.mcmc
#  intervals(mydata.mcmc)
#  plot(mydata.mcmc)