I started this example to help a student implement an epidemic model where we wanted to have the dynamics (Euler steps) occurring on a faster time scale/more frequently than the observations. I ended up doing this by coding missing data/intermediate time steps as negative values and telling Stan to skip the predicted-observed comparison for those cases. (I think this is actually fairly straightforward way to deal with missing data (there is various discussion of how to model missing data in Stan, I haven’t spent long enough to see if my approach is reasonable or not).

In the process I had fun writing some “tidying” methods (i.e., for the broom package) for rstan and rjags models.

library("R2jags")
library("rstan")
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
library("emdbook")  ## for perturb.params(): maybe unnecessary?
library("lattice")
library("plyr")
library("ggplot2"); theme_set(theme_bw())

## if necessary get broom package for tidying
## (BMB's development version has tidiers for 'rjags', 'rstan' objects)
if (FALSE) {
  library("devtools")
  install_github("bbolker/broom")
}
library("broom")

Simulation function:

simlogistdata <- function(seed=1001,
                          r=1,K=10,n0=1,t0=1,tot.t=10,dt=0.5,
                          sd.proc=1,sd.obs=1) {
    if (!is.null(seed)) set.seed(seed)
  tvec <- seq(1,tot.t,by=dt)
  n <- length(tvec)
  y <- numeric(n)
  ydet<-numeric(n)
  y[1] <- n0
  ydet[1] <- n0
  e.proc <- rnorm(n,mean=0,sd=sd.proc)
  e.obs <- rnorm(n,mean=0,sd=sd.obs)
  for (i in 2:n) {
      ydet[i] <- ydet[i-1]+r*ydet[i-1]*(1-ydet[i-1]/K)*dt
      y[i] <- y[i-1]+(r*y[i-1]*(1-y[i-1]/K))*dt
      y[i] <- y[i] + e.proc[i-1]
      ## note we have to be careful here if we are doing anything
      ## non-trivial with dt; nothing is done to ensure that
      ## variance of e.proc scales appropriately with dt ...
  }
    y.procobs <- y+e.obs
    return(cbind(tvec,y.procobs))
}
## now sim data ...
set.seed(102)
s1 <- simlogistdata(tot.t=10,dt=0.2)
plot(y.procobs~tvec,data=s1)

JAGS preliminaries.

model {
  ## initial values
  t[1] <- n0
  o[1] ~ dnorm(t[1],tau.obs)
  ## step through observations
  for (i in 2:N) {
     v[i] <- t[i-1]+r*t[i-1]*(1-t[i-1]/K)*dt
     t[i] ~ dnorm(v[i],tau.proc)
     o[i] ~ dnorm(t[i],tau.obs)
  }
  ## aux variables
  sd_proc <- 1/sqrt(tau.proc)
  sd_obs <- 1/sqrt(tau.obs)
  ## priors
  r ~ dunif(0.1,maxr)
  K ~ dgamma(0.005,0.005)
  tau.obs ~ dgamma(0.005,0.005)
  tau.proc ~ dgamma(0.005,0.005)
  n0 ~ dgamma(1,n0rate)
}
parameters <- c("r","K","sd_obs","sd_proc","n0")
inits <- perturb.params(
    list(n0=1,r=0.2,K=10,tau.obs=1,tau.proc=1),
    alt=list(r=c(0.1,0.4),tau.obs=3,tau.proc=3))

jagsdata <- list(o=s1[,2],dt=0.2,N=nrow(s1),n0rate=1,
                 maxr=5)

Run JAGS:

j1 <- jags(data=jagsdata,
           inits, param=parameters,
           model="logist.bug",
           n.chains=length(inits), n.iter=15000,
           progress.bar="none",  ## quiet ...
           n.burnin=5000,n.thin=37)
## module glm loaded

Diagnose (trace plots):

j1bugs <- as.mcmc.bugs(j1$BUGS)
xyplot(j1bugs)