# This template assumes that you have installed JAGS # and the R packages car, rjags and pscl. # MCMC settings n.adapt <- 100 n.burn <- 100 n.sample <- 500 # load rjags library(rjags) library(car) library(pscl) # convert the data format of the provided 109th # US Senate roll call object to a matrix of # 1s (yea), 0s (nay), and NA (abstention, absence, etc). Y <- recode(s109$votes,"1:3=1;4:6=0;7:9=NA;0=NA") N <- nrow(Y) # extract number of legislators M <- ncol(Y) # extract number of votes SenatorParty <- s109$legis.data$party SenatorName <- rownames(s109$legis.data) jagsModelSpecification <- ' model{ # specify likelihood of vote conditional # on vote parameters, ideal points (theta) and legislator specific variances (sigma^2) for (i in 1:N){ for(j in 1:M){ Y[i,j] ~ dbern(pi[i,j]) pi[i,j] <- phi((alpha[j] + beta[j]*theta[i])/sigma[i]) } } # specify priors for legislator ideal points and variances for (i in 1:N){ theta[i] ~ dnorm(0,1) sigma[i] ~ dgamma(nu,nu) } # specify priors for vote parameters for (j in 1:M){ alpha[j] ~ dnorm(0,0.25) beta[j] ~ dnorm(0,0.25) } nu ~ dunif(0,1000) } ' # specify data to be passed to JAGS jagsData <- list( Y = Y, N = N, M = M ) # specify initial values of ideal points with Republicans at 1, Dem/Ind at -1 jagsInits <- list( theta = 2*(SenatorParty == "R")-1 ) # create the JAGS model object jagsModel <- jags.model( file = textConnection(jagsModelSpecification), # specification of the model data= jagsData, # data for the model inits= jagsInits, # initial values for ideal points n.chains=1, # number of MCMC chains n.adapt=n.adapt # number of iterations of MCMC adaptation, not saved ) # run the MCMC for a while without saving the samples update(jagsModel,n.iter=n.burn) # number of iterations of MCMC burn in, not saved # save sample from the posterior distribution jagsSample <- jags.samples( jagsModel, # model to sample from variable.names=c("theta","sigma","nu"), # names of all variables to be saved n.iter = n.sample # number of iterations to sample ) # extract the samples for theta and sigma (the final index # corresponds to multiple MCMC chain, we just # want the first/only one here.) theta.mcmc <- jagsSample$theta[,,1] sigma.mcmc <- jagsSample$sigma[,,1] # calculate the posterior mean theta.est <- apply(theta.mcmc,1,mean) sigma.est <- apply(sigma.mcmc,1,mean) # calculate the posterior standard deviation (i.e. standard error) theta.se <- apply(theta.mcmc,1,sd) sigma.se <- apply(sigma.mcmc,1,sd) # calculate the 95% posterior interval # (we should really get a larger sample in # order to get stable estimates for this) theta.int <- apply(theta.mcmc,1,quantile,c(0.025,0.975)) sigma.int <- apply(sigma.mcmc,1,quantile,c(0.025,0.975)) nu.est <- mean(jagsSample$nu) nu.se <- sd(jagsSample$nu) nu.int <- quantile(jagsSample$nu,c(0.025,0.975)) # create an data frame with stats for each legislator legislator.summary <- data.frame( name=SenatorName, party=SenatorParty, colour=as.character(c("blue","green","red")[as.numeric(SenatorParty)]), theta.est=theta.est, theta.se=theta.se, theta.lo=theta.int[1,], theta.hi=theta.int[2,], sigma.est=sigma.est, sigma.se=sigma.se, sigma.lo=sigma.int[1,], sigma.hi=sigma.int[2,], stringsAsFactors=FALSE) par(mfrow=c(1,2)) # side-by-side caterpiller plots for ideal points and variances # sort data frame on legislator ideal point estimate legislator.summary <- legislator.summary[order(legislator.summary$theta.est),] # create a caterpiller plot for legislator ideal points plot(0,0,type="n",xlim=c(-2.5,2.5),ylim=c(0,N+1),xlab="Ideal Point",ylab="",axes=FALSE) axis(1) for (i in 1:N){ lines(c(legislator.summary$theta.lo[i], legislator.summary$theta.hi[i]),rep(i,2),lwd=1,col= legislator.summary$colour[i]) points(legislator.summary$theta.est[i],i,pch=16,col= legislator.summary$colour[i],cex=0.5) if (legislator.summary$theta.est[i] > 0){ text(legislator.summary$theta.lo[i],i, legislator.summary$name[i],pos=2,cex=0.25) } else { text(legislator.summary$theta.hi[i],i, legislator.summary$name[i],pos=4,cex=0.25) } } # sort data frame on legislator variance estimate legislator.summary <- legislator.summary[order(legislator.summary$sigma.est),] # create a caterpiller plot for legislator variances plot(0,0,type="n",xlim=c(0,3),ylim=c(0,N+1),xlab="Legislator sqrt(Variances)",ylab="",axes=FALSE) axis(1) for (i in 1:N){ lines(c(legislator.summary$sigma.lo[i], legislator.summary$sigma.hi[i]),rep(i,2),lwd=1,col= legislator.summary$colour[i]) points(legislator.summary$sigma.est[i],i,pch=16,col= legislator.summary$colour[i],cex=0.5) if (legislator.summary$sigma.est[i] > 1){ text(legislator.summary$sigma.lo[i],i, legislator.summary$name[i],pos=2,cex=0.25) } else { text(legislator.summary$sigma.hi[i],i, legislator.summary$name[i],pos=4,cex=0.25) } }