RJAGS - How to pass more complex functions in BUGS file

126 Views Asked by At

My goal is to basically migrate this code to R. All the preprocessing wrt datasets has been already done, now however I am stuck in writing the "model" file. As a first attempt, and for the sake of clarity, I wrote the code which is shown below in R language.

What I want to do is to run an MCMC to have an estimate of the parameter R_t, given the daily reported data for Italian Country. The main steps that have been pursued are:

  1. Sample an array parameter, namely the log(R_t), from a Gaussian RW distribution
Gauss_RandomWalk <- function(N, x0, mu, variance) {   
z <- cumsum(rnorm(n=N, mean=mu, sd=sqrt(variance)))   
t <- 1:N   
x <- (x0 + t*mu + z)   
return(x) 
}

log_R_t <- Gauss_RandomWalk(tot_dates, 0., 0., 0.035**2)
R_t_candidate <- exp(log_R_t)
  1. Compute some quantities, that are function of this sampled parameters, namely the number of infections. This dependence is quite simple, since it is linear algebra:
infections <- rep(0. , tot_dates)
infections[1] <- exp(seed)
for (t in 2:tot_dates){
  infections[t] <- sum(R_t_candidate * infections * gt_to_convolution[t-1,])
}

  1. Convolve the array I have just computed with a delay distribution (onset+reporting delay), finally rescaling it by the exposure variable:
test_adjusted_positive <- convolve(infections, delay_distribution_df$density, type = "open")
test_adjusted_positive <- test_adjusted_positive[1:tot_dates]
positive <- round(test_adjusted_positive*exposure) 
  1. Compute the Likelihood, which is proportional to the probability that a certain set of data was observed (i.e. daily confirmed cases), by sampling the aforementioned log(R_t) parameter from which the variable positive is computed.
likelihood <- dnbinom(round(Italian_data$daily_confirmed), mu = positive, size = 1/6) 

Finally, here we come to my BUGS model file:

model {

  #priors as a Gaussian RW
  log_rt[1] ~ dnorm(0, 0.035)
  log_rt[2] ~ dnorm(0, 0.035)
  for (t in 3:tot_dates) {
    log_rt[t] ~ dnorm(log_rt[t-1] + log_rt[t-2], 0.035)
    R_t_candidate[t] <- exp(log_rt[t])
  }
  
    # data likelihood
    for (t in 2:tot_dates) {
        infections[t] <- sum(R_t_candidate * infections * gt_to_convolution[t-1,])
    }
    
    test_adjusted_positive <- convolve(infections, delay_distribution)
  test_adjusted_positive <- test_adjusted_positive[1:tot_dates]
  positive <- test_adjusted_positive*exposure

  for (t in 2:tot_dates) {
        confirmed[t] ~ dnbinom( obs[t], positive[t], 1/6) 
    }
}

where gt_to_convolution is a constant matrix, tot_dates is a constant value and exposure is a constant array.

When trying to compile it through:

data <- NULL
data$obs <- round(Italian_data$daily_confirmed)
data$tot_dates <- n_days
data$delay_distribution <- delay_distribution_df$density
data$exposure <- exposure
data$gt_to_convolution <- gt_to_convolution

inits <- NULL
inits$log_rt  <- rep(0, tot_dates)

library (rjags)
library (coda)

set.seed(1995)

model <- "MyModel.bug"
jm <- jags.model(model , data, inits)

It raises the following raising error:

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Deleting model

Error in jags.model(model, data, inits) : RUNTIME ERROR:
Compilation error on line 19.
Possible directed cycle involving test_adjusted_positive

Hence I am not even able to debug it a little, even though I'm pretty sure there is something wrong more in general but I cannot figure out what and why.

At this point, I think the best choice would be to implement a Metropolis Algorithm myself according to the likelihood above, but obviously, I would way much more prefer to use an already tested framework that is BUGS/JAGS, this is the reason why I am asking for help.

0

There are 0 best solutions below