inst/doc/R2admb.R

In R2admb: 'ADMB' to R Interface Functions

## ----opts,echo=FALSE,message=FALSE--------------------------------------------
library("knitr")
opts_chunk$set(tidy=FALSE)

## ----do_admb_fake,eval=FALSE--------------------------------------------------
#  fitted.model <- do_admb(fn,data,params,
#                          run.opts=run.control(checkparam="write",
#                          checkdata="write"))

## ----libs,message=FALSE, results = "hide"-------------------------------------
library("R2admb")
library("ggplot2") ## for pictures
theme_set(theme_bw())  ## cosmetic
zmargin <- theme(panel.spacing=grid::unit(0,"lines"))
library("bbmle")

## ----dat1---------------------------------------------------------------------
ReedfrogSizepred <- 
  data.frame(TBL = rep(c(9,12,21,25,37),each=3),
             Kill = c(0,2,1,3,4,5,0,0,0,0,1,0,0,0,0L))

## ----mlefit-------------------------------------------------------------------
m0 <- mle2(Kill~dbinom(c*((TBL/d)*exp(1-TBL/d))^g,size=10),
           start=list(c=0.45,d=13,g=1),data=ReedfrogSizepred,
           method="L-BFGS-B",
           lower=c(c=0.003,d=10,g=0),
           upper=c(c=0.8,d=20,g=60),
           control=list(parscale=c(c=0.5,d=10,g=1)))

## ----predvals-----------------------------------------------------------------
TBLvec = seq(9.5,36,length=100)
predfr <- 
  data.frame(TBL=TBLvec,
             Kill=predict(m0,newdata=data.frame(TBL=TBLvec)))

## ----fig1,echo=FALSE----------------------------------------------------------
g1  <- ggplot(ReedfrogSizepred,
              aes(x=TBL,y=Kill/10))+
    geom_point()+stat_sum(aes(size=after_stat(n)))+
    geom_smooth(method="loess",formula=y~x)+
    labs(size="n",x="Size (total body length",
         y="Proportion killed")+
    coord_cartesian(ylim=c(-0.05,0.55))
startest <- stat_function(fun = function(x) { 0.45*((x/13)*exp(1-x/13)) },
                          lty=2,colour="red")
print(g1+startest+
      geom_line(data=predfr,colour="purple",lty=2))

## ----admbfit_getruns,echo=FALSE-----------------------------------------------
load_doc <- function(x) {
    ##    Define a bit of magic to load results from previous runs: 
    ## would like parent.frame(2) above instead of
    ## hacking around with global environment
    ## ... but doesn't work as I thought?
    load(system.file("doc",x,package="R2admb"),envir=.GlobalEnv) ##
}
zz <- load_doc("Reedfrog_runs.RData")

## ----setup_admb,eval=FALSE----------------------------------------------------
#  setup_admb()

## ----rfs_setup----------------------------------------------------------------
rfs_params <- list(c=0.45,d=13,g=1) ## starting parameters
rfs_bounds <- list(c=c(0,1),d=c(0,50),g=c(-1,25)) ## bounds
rfs_dat <- c(list(nobs=nrow(ReedfrogSizepred),
                  nexposed=rep(10,nrow(ReedfrogSizepred))),
             ReedfrogSizepred)

## ----admbfit_fake,eval=FALSE--------------------------------------------------
#  m1 <- do_admb("ReedfrogSizepred0",
#                data=rfs_dat,
#                params=rfs_params,
#                bounds=rfs_bounds,
#                run.opts=run.control(checkparam="write",
#                checkdata="write",clean=FALSE))
#  unlink(c("reedfrogsizepred0.tpl",
#           "reedfrogsizepred0_gen.tpl",
#           "reedfrogsizepred0")) ## clean up leftovers

## ----rffit2,eval=FALSE--------------------------------------------------------
#  rfs_params$g <- 2
#  m2 <- do_admb("ReedfrogSizepred1",
#                data=rfs_dat,
#                params=rfs_params,
#                bounds=rfs_bounds,
#                run.opts=run.control(checkparam="write",
#                checkdata="write"))

## ----basic--------------------------------------------------------------------
m1

## ----coef---------------------------------------------------------------------
coef(m1)

## ----summary------------------------------------------------------------------
summary(m1)

## ----vcov---------------------------------------------------------------------
vcov(m1)

## ----others-------------------------------------------------------------------
c(logLik(m1),deviance(m1),AIC(m1))

## ----profrun,eval=FALSE,tidy=FALSE--------------------------------------------
#  m1P <- do_admb("ReedfrogSizepred0",
#                 data=c(list(nobs=nrow(ReedfrogSizepred),
#                 nexposed=rep(10,nrow(ReedfrogSizepred))),
#                 ReedfrogSizepred),
#                 params=rfs_params,
#                 bounds=rfs_bounds,
#                 run.opts=run.control(checkparam="write",
#                 checkdata="write"),
#                 profile=TRUE,
#                 workdir=".",
#                 profile.opts=list(pars=c("c","d","g")))

## ----mleprof,cache=TRUE-------------------------------------------------------
m0prof <- profile(m0)

## ----profcalcs2,echo=FALSE,warning=FALSE--------------------------------------
tmpf <- function(p,w="prof") {
  pp <- log(m1P$prof[[p]][[w]][,2])
  pp <- max(pp)-pp
  data.frame(param=p,z=sqrt(2*pp),
             par.vals.c=NA,par.vals.d=NA,par.vals.g=NA,
             focal=m1P$prof[[p]][[w]][,1])
}
quadf <- function(p) {
    m <- coef(m0)[p]
    se <- stdEr(m0)[p]
    pvec <- seq(m-3*se,m+3*se,length=31)
    data.frame(param=p,z=abs(pvec-m)/se,
                focal=pvec,
                par.vals.c=NA,par.vals.d=NA,par.vals.g=NA)
}
proflist <- do.call(rbind,lapply(list("c","d","g"),tmpf))
profnlist <- do.call(rbind,lapply(list("c","d","g"),tmpf,w="prof_norm"))
quadlist <- do.call(rbind,lapply(list("c","d","g"),quadf))
pdat <- rbind(cbind(as.data.frame(m0prof),method="mle2"),
              cbind(proflist,method="ADMB"),
              cbind(profnlist,method="ADMB_norm"),
              cbind(quadlist,method="Wald"))

## ----profpic,echo=FALSE,fig.width=8,fig.height=3.2,warning=FALSE--------------
ggplot(pdat,aes(x=focal,y=abs(z),group=method,colour=method))+geom_line()+
      geom_point(alpha=0.5)+
  facet_grid(.~param,scale="free_x")+ylim(0,3)+xlab("")+
      ylab(expression(Delta(sqrt(-2*L))))+
      geom_hline(yintercept=1.96,lty=2)+zmargin

## ----admbfakemc,eval=FALSE,tidy=FALSE-----------------------------------------
#  m1MC <- do_admb("ReedfrogSizepred0",
#                  data=rfs_dat,
#                  params=rfs_params,
#                  bounds=rfs_bounds,
#                  run.opts=run.control(checkparam="write",
#                    checkdata="write"),
#                  mcmc=TRUE,
#                  mcmc.opts=mcmc.control(mcmcpars=c("c","d","g")))
#  ## clean up leftovers:
#  unlink(c("reedfrogsizepred0.tpl",
#           "reedfrogsizepred0_gen.tpl",
#           "reedfrogsizepred0"))

## ----mchistplot,message=FALSE-------------------------------------------------
plot(m1MC$hist)

## ----coda---------------------------------------------------------------------
library("coda")
mmc <- as.mcmc(m1MC$mcmc)

## ----mctraceplot--------------------------------------------------------------
library("lattice")
xyplot(mmc)

## ----rgdiags------------------------------------------------------------------
raftery.diag(mmc)
geweke.diag(mmc)

## ----gdiag,echo=FALSE---------------------------------------------------------
gd <- geweke.diag(mmc)
gd1 <- gd[["z"]][1]

## ----effsize------------------------------------------------------------------
effectiveSize(mmc)

## ----hpd----------------------------------------------------------------------
HPDinterval(mmc)

## ----mcdensplot---------------------------------------------------------------
densityplot(mmc)

## ----lme4,message=FALSE-------------------------------------------------------
library(lme4)
if (as.numeric(R.version$major)<3) {
## FIXME
    gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
                 family = binomial, data = cbpp)
}

## ----toymats------------------------------------------------------------------
X <- model.matrix(~period,data=cbpp)
Zherd <- model.matrix(~herd-1,data=cbpp)

## ----toydat-------------------------------------------------------------------
tmpdat <- list(X=X,Zherd=Zherd,
                 incidence=cbpp$incidence,size=cbpp$size,
                 nobs=nrow(cbpp))

## ----loadtoy1,echo=FALSE------------------------------------------------------
zz2 <- load_doc("toy1_runs.RData")

## ----fakerun2,eval=FALSE------------------------------------------------------
#  d1 <- do_admb("toy1",
#                data=tmpdat,
#                params=list(beta=rep(0,ncol(X)),sigma_herd=0.1),
#                bounds=list(sigma_herd=c(0.0001,20)),
#                re=list(u_herd=ncol(Zherd)),
#                run.opts=run.control(checkdata="write",checkparam="write"),
#                mcmc=TRUE,
#                mcmc.opts=mcmc.control(mcmc=20,mcmcpars=c("beta","sigma_herd")))

## ----testprofinput,echo=FALSE,eval=FALSE--------------------------------------
#  do_admb("toy1", data=tmpdat,
#          params=list(beta=rep(0,ncol(X)),sigma_herd=0.1),
#          bounds=list(sigma_herd=c(0.0001,20)),
#          re=list(u_herd=ncol(Zherd)),
#          run.opts=run.control(checkdata="write",checkparam="write",
#            clean=FALSE))
#  run_admb("toy1_gen",profile=TRUE)
#  read_admb("toy1_gen",profile=TRUE)

## ----coefsumlmer--------------------------------------------------------------
## FIXME
## coef(summary(gm1))

## ----coefsumadmb--------------------------------------------------------------
coef(summary(d1))[1:5,]

## ----ranefs-------------------------------------------------------------------
## FIXME
## plot(ranef(gm1)$herd[,1],coef(d1)[6:20]*coef(d1)["sigma_herd"],
##     xlab="glmer estimate",ylab="ADMB estimate")
## abline(a=0,b=1)

## ----mcmcconf-----------------------------------------------------------------
detach("package:lme4") ## HPDinterval definition gets in the way
HPDinterval(as.mcmc(d1$mcmc[,6:20]))