# DOWNLOAD http://dr-k-lo.blogspot.com/2013/03/problems-loading-jags-into-r.html, RUN INSTALL SCRIPT
# MAKE SURE R IS VERSION 3.1.0
# LOAD LIBRARIES
lapply(c("rjags","arm","coda","superdiag","R2WinBUGS","R2jags","lme4"),library, character.only=TRUE)

source("/Users/jgill/Article.JPART/asap.rjags.dat")  # READS IN THE LIST VERSION OF THE DATA FOR asap.jags.list
names(asap.jags.list) <- c( "state.id", "contracting", "gov.influence", "leg.influence",
    "elect.board", "years.tenure", "education", "party.ID", "category2", "category3",
    "category4", "category5", "category6", "category7", "category8", "category9",
    "category10", "category11", "category12", "med.time", "medt.contr", "grp.influence",
    "gov.ideology", "lobbyists", "nonprofits", "STATES", "SUBJECTS")

# DEFINE THE MODEL
asap.model2.rjags  <- function()  {
for (i in 1:SUBJECTS) {
    mu[i] <- alpha[state.id[i]]
             + beta[1]*contracting[i] + beta[2]*gov.influence[i] + beta[3]*leg.influence[i] 
             + beta[4]*elect.board[i] + beta[5]*years.tenure[i]  + beta[6]*education[i]
             + beta[7]*party.ID[i]    + beta[8]*category2[i]     + beta[9]*category3[i] 
             + beta[10]*category4[i]  + beta[11]*category5[i]    + beta[12]*category6[i] 
             + beta[13]*category7[i]  + beta[14]*category8[i]    + beta[15]*category9[i] 
             + beta[16]*category10[i] + beta[17]*category11[i]   + beta[18]*category12[i] 
             + beta[19]*med.time[i]   + beta[20]*medt.contr[i]
      grp.influence[i] ~ dnorm(mu[i],tau.y)
    }
    for (j in 1:STATES) {
      eta[j]  <- gamma[1]*gov.ideology[j] + gamma[2]*lobbyists[j] + gamma[3]*nonprofits[j]
      alpha[j] ~ dnorm(eta[j],tau.alpha)
    }
    beta[1]   ~ dnorm(0.070,1);     # PRIOR MEANS FROM KELLEHER AND YACKEE 2009, MODEL 3
    beta[2]   ~ dnorm(-0.054,1);
    beta[3]   ~ dnorm(0.139,1);
    beta[4]   ~ dnorm(0.051,1);      
    beta[5]   ~ dnorm(0.017,1);
    beta[6]   ~ dnorm(0.056,1);
    beta[7]   ~ dnorm(0.039,1);
    beta[8]   ~ dnorm(0.0,1);       # DIFFUSE PRIORS 
    beta[9]   ~ dnorm(0.0,1);      
    beta[10]  ~ dnorm(0.0,1);      
    beta[11]  ~ dnorm(0.0,1);      
    beta[12]  ~ dnorm(0.0,1);      
    beta[13]  ~ dnorm(0.0,1);      
    beta[14]  ~ dnorm(0.0,1);      
    beta[15]  ~ dnorm(0.0,1);       
    beta[16]  ~ dnorm(0.0,1);        
    beta[17]  ~ dnorm(0.0,1);      
    beta[18]  ~ dnorm(0.0,1);    
    beta[19]  ~ dnorm(0.184,1);     # PRIOR MEANS FROM KELLEHER AND YACKEE 2009, MODEL 3
    beta[20]  ~ dnorm(0.156,1);
    gamma[1]  ~ dnorm(0.0,1);       # DIFFUSE PRIORS 
    gamma[2]  ~ dnorm(0.0,1);
    gamma[3]  ~ dnorm(0.0,1);
    tau.y     ~ dgamma(1.0,1);
    tau.alpha ~ dgamma(1.0,1);
}

# SAVE MODEL TO A FILE
write.model(asap.model2.rjags, "Article.JPART/asap.model2.rjags")  

# SETUP INITIAL VALUES AND PARAMETER NAMES
asap.inits <- function(){ list("tau.y" = 10, "tau.alpha" = 10,
    "beta" = rep(1,20), "gamma" = c(1,1,1)) }
asap.params <- c("beta","gamma","tau.y","tau.alpha")

# RUN THE SAMPLER AND COLLECT coda SAMPLES
asap.out <- jags(data=asap.jags.list, inits=asap.inits, asap.params, n.iter=5000, 
    model="Article.JPART/asap.model2.rjags", DIC=TRUE)

# USE THE autojags() FUNCTION TO RUN UNTIL G&R < 1.1
#asap.out2 <- autojags(asap.out) 

# OR update() TO RUN LONGER
asap.out3 <- update(asap.out, n.iter=200000)
print(asap.out3)
asap.mcmc3 <- as.mcmc(asap.out3)

# CHECK CONVERGENCE
superdiag(as.mcmc.list(asap.mcmc3), burnin=0)