vignettes/chytrid_dede_example.R

## ----opts, echo=FALSE, results='hide'------------------------------------
knitr::opts_chunk$set(dev="png", dpi=150)

## ----install1, eval=FALSE------------------------------------------------
#  install.packages("deBInfer")

## ----install2, eval=FALSE------------------------------------------------
#  if (require("devtools")){
#    #install deBInfer from github
#    devtools::install_github("pboesu/debinfer")
#  }

## ----loadlib, message=FALSE----------------------------------------------
library(deBInfer)

## ----dde-def-------------------------------------------------------------
#dede version
CSZ.dede<-function(t,y,p){

  sr    <- p["sr"]
  fs    <- p["fs"]
  ds    <- p["ds"]
  eta   <- p["eta"]
  Tmin  <- p["Tmin"]
  muz <- p["muz"]

  Rs <- Ms <- 0
  lag1 <- lag2 <- 0

  if (t>Tmin) {
    lag1 <- lagvalue(t - Tmin)
    Rs <- sr * fs * lag1[1]
    }

  phiZ <- eta * y[2]
  dy1 <- -(muz + sr) * y[1]
  dy2 <- Rs - Ms - ds * y[2]
  dy3 <- phiZ - (muz + sr) * y[3]

  if(y[1]<0) dy1<-0
    if (y[2] < 0) {
    dy2 <- Rs - Ms
    dy3 <- -(muz + sr) * y[3]
    }
    if (y[3] < 0) {
    dy3 <- dy3 + (muz + sr) * y[3]
    }
    
    list(c(dy1,dy2,dy3))
    }

## ----data----------------------------------------------------------------
#load chytrid data
data(chytrid)
#have a look at the variables
head(chytrid)
#plot the data
plot(chytrid, xlab = "Time (days)", ylab = "Zoospores x 10e4", xlim = c(0,10))

## ----obs-model-----------------------------------------------------------
# observation model
chytrid_obs_model <- function(data, sim.data, samp) {
  
  ec <- 0.01
  llik.Z <- 0
  for(i in unique(data$time)){
    try(llik.Z <- llik.Z + sum(dpois(data$count[data$time == i], 
                                lambda = (sim.data[,"Z"][sim.data[,"time"] == i] + ec),
                                log = TRUE)))
  }
  llik <- llik.Z
  return(llik)
}

## ----vars----------------------------------------------------------------
sr <- debinfer_par(name = "sr", var.type = "de", fixed = FALSE,
                   value = 2, prior = "gamma", hypers = list(shape = 5, rate = 1),
                   prop.var = c(3,4), samp.type = "rw-unif")

fs <- debinfer_par(name = "fs", var.type = "de", fixed = FALSE,
                   value = 0.5, prior = "beta", hypers = list(shape1 = 1, shape2 = 1),
                   prop.var = 0.01, samp.type = "ind")

ds <- debinfer_par(name = "ds", var.type = "de", fixed = FALSE,
                   value = 2, prior = "gamma", hypers = list(shape = 1, rate = 1),
                   prop.var = 0.1, samp.type = "rw")

muz <- debinfer_par(name = "muz", var.type = "de", fixed = FALSE,
                    value = 1, prior = "gamma", hypers = list(shape = 5, rate = 1),
                    prop.var = c(4,5), samp.type = "rw-unif")

eta <- debinfer_par(name = "eta", var.type = "de", fixed = FALSE,
                    value = 10, prior = "gamma", hypers = list(shape = 1, rate = 0.25),
                    prop.var = 5, samp.type = "rw")

Tmin <- debinfer_par(name = "Tmin", var.type = "de", fixed = FALSE,
                     value = 3, prior = "unif", hypers = list(min = 2, max = 6),
                     prop.var = 0.05, samp.type = "rw")


# ----inits---------------------------------------------------------------
C <- debinfer_par(name = "C", var.type = "init", fixed = TRUE, value = 120)
S <- debinfer_par(name = "S", var.type = "init", fixed = TRUE, value = 0)
Z <- debinfer_par(name = "Z", var.type = "init", fixed = TRUE, value = 0)

## ----mcmc-setup----------------------------------------------------------
# ----setup---------------------------------------------------------------
mcmc.pars <- setup_debinfer(sr, fs, ds, muz, eta, Tmin, C, S, Z)
pboesu/debinfer documentation built on April 20, 2018, 4:34 a.m.