#################################################################### ## a linear regression example ## section 6.3, pp256ff library(pscl) data(absentee) attach(absentee) ## create variables for regression analysis y <- (absdem - absrep)/(absdem + absrep)*100 x <- (machdem - machrep)/(machdem + machrep)*100 ## exploratory graph, regression plot(y~x, las=1, xlab="Democratic Margin, Machine Ballots", ylab="Democratic Margin, Absentee Ballots") points(x[22],y[22], col="red", cex=3,pch=16) abline(v=0,lty=2) # ols linear regression, subset to 1st 21 observations ols <- lm(y~x, subset=c(rep(TRUE,21),FALSE)) abline(ols,col="blue",lwd=3) ## overlay regression fit ## environment for passing to JAGS forJags <- list(y=y[1:21], x=x[1:21], n=21, xstar=x[22]) ## initial values, list containing one list, one chain inits <- list(list(beta=c(0,0), sigma=5, ystar=0, .RNG.seed=1234, ## seed the RNG explicitly .RNG.name="base::Mersenne-Twister")) ## p259 book library(rjags) ## compilaton, no Gibbs sampler yet foo <- jags.model(file="regression.bug", n.chains=1, data=forJags, inits=inits) update(foo,5E3) # now we get Gibbs sampling out <- coda.samples(model=foo, variable.names=c("beta","sigma","ystar"), n.iter=50E3) summary(out) ## posterior predictive distribution for ## disputed election ystar <- out[[1]][,"ystar"] hist(ystar,51, probability=TRUE, main="Posterior Predictive Density,\nDisputed Absentee Ballot Result") abline(v=y[22],col="red",lwd=3) ## Bayesian p-value table(ystar>y[22]) mean(ystar>y[22])