SSVS and Update Coefficients en bloc

[John]

Hey all,

I'm trying to do variable selection within a zero-truncated poisson glm. Some of the variables are higher order terms and I want to enforce marginality, such that a higher order term can't enter into the model unless its lower order analogs are also included. Obviously these variables will be correlated with one another, so I've scaled them and want to update them en bloc within JAGS. I've done this using the following code (I've included year as a random effect and only showed the code with one co-variate, although I have many):

data{
    C <- 10000

    for(i in 1:n){
        zeros[i] <- 0
    }

    for(i in 1:3){
        b[i] <- 0
    }
}
model{

    for(i in 1:n){

        log(mu[i]) <- log(offset[i]) + alpha[year[i]] + inprod(beta[year[i], ],x[i,])

        #Zero-truncated Poisson
        ll[i] <-  y[i]*log(mu[i])-mu[i]-loggam(y[i]+1)-log(1-exp(-mu[i]))

        phi[i] <- -ll[i] + C
        zeros[i] ~ dpois(phi[i])

    }

    for(i in 1:3){       
        w_pre[i] ~ dbern(0.5)
    }

    #Co-variate 1
    w[1] <- w_pre[1]
    #Co-variate 1^2
    w[2] <- w_pre[1]*w_pre[2]
    #Co-variate 1^3
    w[3] <- w_pre[1]*w_pre[2]*pre[3]

    #Co-variate 1
    B1[1,1] <- tau[w[1]+1]
    B1[1,2] <- 0
    B1[1,3] <- 0
    B1[2,1] <- 0
    B1[2,2] <- tau[w[2]+1]
    B1[2,3] <- 0
    B1[3,1] <- 0
    B1[3,2] <- 0
    B1[3,3] <- tau[w[3]+1]

    for(i in 1:years){ 
        alpha[i] ~ dnorm(0,tau[2])
        beta[i,1:3] ~ dmnorm(b[1:3],B1[1:3,1:3])
    }

    tau[1] <- 10000
    tau[2] <- tauBeta
    tauBeta <- pow(sigmaBeta,-2)
    sigmaBeta ~ dunif(0.0001,10)

}

By updating the coefficients en bloc, it appears that all terms are included equally (for certain variables) such that w[1] = w[2] = w[3] = 1 and the posterior marginal probability is always the same for w[1], w[2], and w[3] (this seems odd). If I don't update the coefficients en bloc this doesn't happen, but mixing is poor. Am I coding something incorrectly, I can't understand why updating coefficients en bloc forces all terms within a block to be included at any iteration. Obviously, this logic doesn't work for exclusion because if a low-order term doesn't make it into the model no higher-order term will either. I also realize the I'm penalizing models with higher order terms from the get go.

Thanks!

[Martyn Plummer]

Hi John

You can fit the zero-truncated Poisson model in JAGS without resorting to the zeros trick

y[i] ~ dpois(mu[i]) T(1, )

Your plan is to force block updating of the beta parameters by defining them as a multivariate normal, even though they are a priori independent. You should check what sampling method is being used (using list.samplers in rjags, or "samplers to <filename>" in the command line interface). I suspect it is a random-walk Metropolis album, which is a lousy method and may not be any better than Gibbs sampling.

I think your problem is that you still have poor mixing, but now the poor mixing is manifesting itself in the indicator variables w[1], w[2], and w[3]. To solve this, you should adopt the approach of Carlin and Chib (Bayesian Model Choice via Markov Chain Monte Carlo Methods, JRSS B, 1995). Instead of having a prior for beta that is concentrated around zero when the corresponding term is out of the model, you should create a pseudo-prior that resembles the posterior distribution when the term is in the model. This allows the Markov chain to jump back and forth more easily between the two states. Calibrating the pseudo-priors requires some test runs but should pay off in terms of mixing.

[John]

Hey Dr. Plummer,

Thanks for the advice, I very much appreciate it. Although, I hesitate to use Carlin and Chib's method as they require you to take a two-stage approach. First fitting all the models separately, then using the posteriors to approximate the pseudo-priors for other models, so as to adjust the prior probability on each model when fitting all models together to ensure they're visited frequently. I have a lot of covariates, and this method is impractical. Do you have any other suggestions?

Thanks for you time.