unexpected invalid model spec when hoping to use Polya-gamma

[paciorek]

I thought my error trapping for model structure would allow this but it doesn't seem to. There are easy workarounds, but would be nice if this worked.

set.seed(1)
n_sites <- 100
n_visits <- 5
rho <- 0.7
x1 <- rnorm(n_sites)
x2 <- rnorm(n_sites) + rho*x1
x3 <- rnorm(n_sites)

beta <- c(1, 0.5, 0, -0.5)
alpha <- c(-.5, 0.25, 1, 0)

prob_present <- expit(beta[1] + beta[2]*x1 + beta[3]*x2 + beta[4]*x3)
prob_detected <- expit(alpha[1] + alpha[2]*x1 + alpha[3]*x2 + alpha[4]*x3)

z <- rbinom(n_sites, 1, prob_present)
y <- matrix(rbinom(n_sites*n_visits, 1, prob_detected*z), n_sites, n_visits)
visited <- apply(y, 1, max)
z_data <- rep(NA, n)
z_data[visited == 1] <- 1

z_init <- z_data
z_init[[is.na](http://is.na/)(z_init)] <- 0

code <- nimbleCode({
   for(j in 1:4)
      beta[j] ~ dnorm(0, sd=10)  # use all three covariates
   for(j in 1:3)
      alpha[j] ~ dnorm(0, sd=10)  # only use first two covariates
   for(i in 1:n_sites) {
      logit(p[i]) <- beta[1] + beta[2]*x1[i]+beta[3]*x2[i] + beta[4]*x3[i]
      z[i] ~ dbern(p[i])
      p_eff[i] <- z[i]*expit(alpha[1]+alpha[2]*x1[i] + alpha[3]*x2[i])
      for(j in 1:n_visits) {
         y[i,j] ~ dbern(p_eff[i])
      }
    }
 })      

model <- nimbleModel(code, data = list(y=y, z=z_data),
   constants = list(x1=x1, x2=x2, x3=x3, n_sites = n_sites, n_visits=n_visits),
   inits = list(beta = rep(0,4), alpha = rep(0,3), z=z_init))

conf <- configureMCMC(model)
conf$removeSamplers(c('alpha','beta'))
conf # binary samplers assigned to z[i]
conf$addSampler('beta','polyagamma')
conf$addSampler('alpha','polyagamma')
conf
mcmc <- buildMCMC(conf)
cmodel <- compileNimble(model)
cmcmc <- compileNimble(mcmc, project = model)
samples <- runMCMC(cmcmc, niter = 5000)

[paul-vdb]

@paciorek I'm also surprised this doesn't work. I thought we checked to make sure p_eff[i] <- z[i]*expit(alpha[1]+alpha[2]*x1[i] + alpha[3]*x2[i]) is okay if z[i] is a bernoulli. I suppose we checked a different case of that.

This minor change version is allowed:

code <- nimbleCode({
   for(j in 1:4)
      beta[j] ~ dnorm(0, sd=10)  # use all three covariates
   for(j in 1:3)
      alpha[j] ~ dnorm(0, sd=10)  # only use first two covariates
   for(i in 1:n_sites) {
      logit(p[i]) <- beta[1] + beta[2]*x1[i]+beta[3]*x2[i] + beta[4]*x3[i]
      z[i] ~ dbern(p[i])
      p_eff[i] <- expit(alpha[1]+alpha[2]*x1[i] + alpha[3]*x2[i])
      for(j in 1:n_visits) {
         y[i,j] ~ dbern(z[i]*p_eff[i])
      }
    }
 })