Hierarchical Poisson Regression

This code fits a hierarchical Poisson regression to the women's
cosponsorship pooled data.   The first-level covariates in the model are
the member's D-NOMINATE first dimension score (a measure of
conservatism), the party of the member (coded one for Democrats and zero
for Republicans), and the member's gender (coded one for women and zero
for men).   The data are discussed by Wolbrecht (2000). This version
includes a set of  posterior predictive distributions to test the
plausibility of the model.

Each first level parameter is independently drawn from a common
distribution where the mean and variance are estimated.  Here the
constant is modeled as a function of three eras -- pre-ERA passage,
post-ERA passage, and post-ERA passage with unlimited cosponsorship. 
The preference measure is modeled as a function of measures of
heterogeneity -- standard deviation of preferences in the House, and a
measure of party distance.

There are four datasets loaded by the program.  The first contains the
level one data as a N by J matrix, and has Y in the first column,
followed by the X matrix (with the constant in the file).  The second
dataset contains a K by 2 matrix that contains the upper and lower
limits for the data matrix.  This matrix defines the clustering scheme. 
The third matrix contains the dataset of covariates for beta[1] in the
second level of the hierarchy, and the fourth contains covariates for
beta[2] in the second level of the hierarchy.

model hier;
{
# Model Level One
for (k in 1:K)
   {
   for (i in lower[k]:upper[k])
      {
      log(mu[i]) <-
         (beta[1,k] * cons[i] +
          beta[2,k] * dnom[i] + 
          beta[3,k] * party[i] + 
          beta[4,k] * gender[i]);
      cospon[i] ~ dpois(mu[i]);
      }
   }

# Model Level Two
for (k in 1:K)
   {
   nu[1,k] <- alpha1[1] + alpha1[2] * era1[k] + alpha1[3] * era2[k];
   nu[2,k] <- alpha2[1] + alpha2[2] * hsdnom[k] + alpha2[3] * ptydist[k];
   nu[3,k] <- alpha3[1];
   nu[4,k] <- alpha4[1];
   for (j in 1:J)
      {
      beta[j,k] ~ dnorm(nu[j,k],gamma[j]);
      }
   }

# Priors
for(m1 in 1: M1)
   {
   alpha1[m1] ~ dnorm(0.0, 0.0001);
   }
for(m2 in 1: M2)
   {
   alpha2[m2] ~ dnorm(0.0, 0.0001);
   }
for(m3 in 1: M3)
   {
   alpha3[m3] ~ dnorm(0.0, 0.0001);
   }
for(m4 in 1: M4)
   {
   alpha4[m4] ~ dnorm(0.0, 0.0001);
   }

for (j in 1:J)
   {
   gamma[j] ~ dgamma(0.001,0.001);
   igamma[j] <- 1/gamma[j];
   }

# Posterior Predictive Distributions
for (k in 1:K)
   {
   pred[1,k] <- exp(beta[1,k] + -0.627 * beta[2,k] + 0 * beta[3,k] + 1 * beta[4,k]); # LibRW
   pred[2,k] <- exp(beta[1,k] + -0.627 * beta[2,k] + 1 * beta[3,k] + 1 * beta[4,k]); # LibDW
   pred[3,k] <- exp(beta[1,k] + -0.627 * beta[2,k] + 0 * beta[3,k] + 0 * beta[4,k]); # LibRM
   pred[4,k] <- exp(beta[1,k] + -0.627 * beta[2,k] + 1 * beta[3,k] + 0 * beta[4,k]); # LibDM
   pred[5,k] <- exp(beta[1,k] + -0.027 * beta[2,k] + 0 * beta[3,k] + 1 * beta[4,k]); # ModRW
   pred[6,k] <- exp(beta[1,k] + -0.027 * beta[2,k] + 1 * beta[3,k] + 1 * beta[4,k]); # ModDW
   pred[7,k] <- exp(beta[1,k] + -0.027 * beta[2,k] + 0 * beta[3,k] + 0 * beta[4,k]); # ModRM
   pred[8,k] <- exp(beta[1,k] + -0.027 * beta[2,k] + 1 * beta[3,k] + 0 * beta[4,k]); # ModDM
   }
}