Pure birth process inference using JAGS

Question

I am trying to use JAGS to infer the birth rate in a (Stochastic) pure birth process.

In the language of chemistry, this model is equivalent to the reaction: X->2X with rate alpha*X (also can be seen as a model of a chain reaction)

This is the R code I'm using to generate the process (at fixed times) and the jags code to make inference of the parameter alpha.

library(rjags)

y <- 1; # Starting number of "individuals"
N <- 25 # number of time samplings
alpha <- 0.2 # per-capita birth rate
# Generate the time series
for(i in 2:N) {
  y<- c(y,y[i-1]+rpois(1,alpha*y[i-1]))
};

# The jags code
model_string <- "model{
  for(i in 2:N) {
    New[i] ~ dpois(alpha*y[i-1])
    y[i] <- y[i-1] + New[i]
  }
  alpha ~ dunif(0, 2)
}"

# Create and run the jags model
model <- jags.model(textConnection(model_string), data = list(y = y,N = N), n.chains = 3, n.adapt= 10000)
update(model, 5000); # Burnin for 10000 samples
mcmc_samples <- coda.samples(model, variable.names=c("alpha"), n.iter=5000)

When I run the code, I get the following error:

Error in jags.model(textConnection(model_string), data = list(y = y, N = N),  : 
  RUNTIME ERROR:
Compilation error on line 4.
y[2] is a logical node and cannot be observed

I have tried different things like putting alpha*y[i-1] in a new variable (say, lambda[i]) or changing the calls to New[i] by New[i-1] but nothing worked. Any idea why is this failing? Another smarter way to do this?

Thank you in advance.

I have found the answer somewhere else. You can read it here: https://sourceforge.net/p/mcmc-jags/discussion/610036/thread/6b159634/ — Donald Lockwood, Sep 29 '16 at 13:42

score 1 · Answer 1 · answered Sep 29 '16 at 14:01

Another solution would be to change how you simulate the data and use a link function with the model.

N <- 25 # number of time samplings
alpha <- 0.2 # log per-capita birth rate

# Generate the time series
steps <- 1:N # simulating 25 steps
log.y<- alpha*steps # the log-scale expected count
expected.y <- exp(log.y) # back to the real scale
y <- rpois(N, expected.y) # add Poisson noise to your expected.y

# The jags code
model_string <- "model{
  for(i in 1:N) {
    y[i] ~ dpois(lambda[i])
    log(lambda[i]) <- log.lambda[i]
    log.lambda[i] <- alpha * i
  }
  alpha ~ dunif(-10, 10)
}"

# Create and run the jags model
model <- jags.model(textConnection(model_string),inits = list(alpha = 1), data = list(y = y,N = N), n.chains = 1, n.adapt= 10000)
update(model, 5000); # Burnin for 10000 samples
mcmc_samples <- coda.samples(model, variable.names=c("alpha"), n.iter=5000)

You can see below that this model also retrieves the parameter correctly (alpha = 0.2).

Taking the exponential of that would give you the birth rate (i.e. exp(0.2) = 1.22), or you could do it within the model and track a derived parameter which is just the exponential of alpha. The model would then be:

model_string <- "model{
  for(i in 1:N) {
    y[i] ~ dpois(lambda[i])
  log(lambda[i]) <- log.lambda[i]
  log.lambda[i] <- alpha * i
  }
  alpha ~ dunif(-10, 10)
  birth.rate <- exp(alpha)
}"

And you would just track birth.rate in the variable.names argument.

This is definitely, a better solution than what I expected. But I don't fully understand why the rate lambda is the exponential of alpha*i, instead of simply alpha*i. — Donald Lockwood, Sep 30 '16 at 10:25
Because this model uses the log-link. The inverse of the log link is the exponential, which takes it back to a scale that is easier to interpret (a rate). The Wikipedia page on [Poisson regression](http://en.wikipedia.org/wiki/Poisson_regression) has sufficient information on this. So , for example, if you wanted to simulate a rate of increase of 1.2 (20% increase from one time step to the next), you would take its log `alpha <- log(1.2)` and use that in the simulation. This model also lets you include covariates as well, while the other solution would not. — mfidino, Sep 30 '16 at 13:50

Pure birth process inference using JAGS

1 Answers1