Bayesian implementation of Aster models

[crush...@gmail.com]

I'  trying to implement an Aster-type model in JAGS and am struggling to understand how to estimate the unconditional effects of each predictor variable on the terminal nodes. For example, considering a simple graphical model with three nodes (z1, z2, and y), where z1 is Bernoulli, z2 is also Bernoulli (conditional on z1 = 1) and y is zero-truncated Poison (conditional on z2 = 1), I have implemented the following model in JAGS:

model {

 C <- 10000

 ## PRIORS
 psi0 ~ dnorm(0, 0.0001)
 psi ~ dnorm(0, 0.0001)    

 gamma0 ~ dnorm(0, 0.0001)
 gamma ~ dnorm(0, 0.0001)  

 beta0 ~ dnorm(0, 0.0001)
 beta ~ dnorm(0, 0.0001)          


# Likelihoods
for(i in 1:N){

  ## Hurdle 1
  logit(phi1[i]) <- psi0 + psi * X[i]

  ## Hurdle 2
  logit(phi2[i]) <- gamma0 + gamma * X[i]

  ## Count model
  count[i] <- (pow(mu[i], y[i]) * exp(-mu[i])) / (y_fact[i] * (1 -  exp(-mu[i])))
  log(mu[i]) <- beta0 + beta * X[i]


  # Ones trick
  p[i] <- Lik[i]/C 
  ones[i] ~ dbern(p[i])

  # Full likelihood
  Lik[i] <- (1 - phi1[i]) * (1 - z1[i]) +
            (1 - phi2[i]) * phi1[i] * (1 - z2[i]) * z1[i] + 
            count[i] * phi2[i] * phi1[i] * z1[i] * z2[i]
}
} # Model
 

Using simulated data, this model does estimate psi, gamma, and beta (i.e. the conditional effect of X on z1, z2, and y). For example, when X is simulated to have a negative infuence on z1 (phi1) but no effect on z2  (phi2, conditional on z1 = 1) or y (f, conditional on z2 = 1), the model correctly estimates gamma and beta as 0.

However, even in this case, individuals with large values of X will have lower counts because they are more likely to have z1 = 0. Can anyone explain how to estimate the overall influence (i.e. unconditional effect) of X on the count process and how to implement those estimates in JAGS?

Data and code to run model:

### Function to generate counts from zero-truncated Poisson
rpoisson0 <- function(n, lambda){
  U <- runif(n)   # the uniform sample
  t <- -log(1 - U*(1 - exp(-lambda))) # the "first" event-times
  T1 <- (lambda - t)   # the set of (T-t)
  X <- rpois(n, T1) + 1 # the final truncated Poisson sample
  return(X)
}

n <- 200
X <- rnorm(150)

psi0 <- -0.5
psi <- -2
exp.z1<- exp(psi0 + psi*X)/(1 + exp(psi0 + psi*X))

z1 <- rbinom(n = n, size = 1, prob = exp.z1)

gamma0 <- 1.5
gamma <- 0
exp.z2 <- exp(gamma0 + gamma*X)/(1 + exp(gamma0 + gamma*X))

z2 <- rbinom(n = n, size = 1, prob = exp.z2) * z1

beta0 <- 1.1
beta <- 0
exp.y <- exp(beta0 + beta * X) 

y <- rpoisson0(n = n, lambda = exp.y) * z2

### Prepare data for JAGS
jags_data <- list(N = length(y), z1 = z1, z2 = z2, y = y, y_fact = factorial(y), ones = rep(1, length(y)), X = X)

### Parameters to monitor
jags_param <- c("psi", "gamma", "beta")

[geyer]

I don't have any experience with JAGS or BUGS.  Most of the models I want to do, they don't handle.  What we are now calling the aster transform, equation (5) in Geyer, Wagenius, and Shaw (2007), which is described more fully in the course notes for my special topics course on aster models slide deck 2, slides 28 through 37, is really important.  If you don't know how to implement it in JAGS, then I don't either.  To me this just shows how limited the JAGS/BUGS notion of "model" is.

If I were going to do a Bayesian analysis of an aster model, I would use the mlogl function in the aster package to calculate the log likelihood and the metrop function in the mcmc package to do the sampling.  See attached example.

Admittedly this is slow because it doesn't do everything in C/C++.  Still it manages to do 10^5 iterations in 10 minutes on a pretty slow laptop.

I used a flat prior for my example.  Your use of normal priors is debatable.  A normal prior has lighter tails than any conjugate prior.  It says you have prior "knowledge" that is more than could be obtained by observing any arbitrary amount of data (no matter how much).  It would not be easy to use conjugate priors in my analysis because the mlogl function insists on checking that the data have integer values for integer variables (like Bernoulli and zero-truncated Poisson).  So I would have to remove these checks to allow for conjugate priors, which give unnormalized posteriors that are just like log likelihoods except for having some "made-up" data added to the real data with the "made-up" data allowed to have any possible mean value and any possible strictly positive real number sample size.  If I were going to do a Bayes function for the aster package, that is what I would do.

But no one before you has asked me about Bayesian analysis of aster models.  If there really is a call for it, I will have to modify the mlogl function to allow non-integer data.