---
title: 'Edpsy 590BAY:  Gibbs Multiple Regression via jags'
author: "C. J. Anderson"
date: "September 19, 2021"
output: 
   html_document:
       toc: yes
---
This document illustrates multiple regression in jag using 1 school from the nels data.  This is an extension of linear regression for 1 predictor to multipl predictors.  


### Step up

The packages that will be used
```{r setup, message=FALSE, warning=FALSE, include=TRUE}
library(ellipse)
library(ggplot2)
library(coda)
library(runjags)
library(rjags)
library(jagsUI)
library(BayesFactor)

set.seed(19843)
```


Set my working directory and read in and check the data
```{r}
setwd("C:/Users/cja/Dropbox/edps 590BAY/Lectures/7 Multiple Linear Regression")
nels <- read.table("nels_school_62821.txt",header=TRUE)
head(nels)
```

# Exploratory and Descriptive Statistics

We should always do some exploratory data analysis (e.g., basic descriptive statistics).   I put them in a little table for future reference
```{r descriptive}
# Descriptive statistics
ybar<- mean(nels$math)
std <- sd(nels$math)
s2  <- var(nels$math)
n   <- length(nels$math)
descriptive <- c(n,ybar,std,s2)
names(descriptive) <- c("N","Ybar","std","s^2")
descriptive
summary(nels$math)
```

A histogram of math scores
```{r}
hist(nels$math, main="NELS Math Scores", xlab="Math Scores", breaks=15)
```

Since we're adding more predictor variables we should also compute correlations.  The following gets the correlation matrix 
```{r correlations}
# Correlations between potential predictors ---> see www.cookbook-r.com
xs <- data.frame(nels$math,nels$homework,nels$paredu,nels$ses)
ctab <- cor(xs)
round(ctab,2)     # round correlations to 2 decimials for easier viewing
```

Or if you prefer a figure of a correlation matrix
```{r corPlot}
# Figure of correlations --- grey scale
plotcorr(ctab,mar=c(0.1,0.1,0.1,0.1))
```

I like color, so
```{r corrPlotColor}
# Figures of correlations in color 
colorfun <- colorRamp(c("#CC0000","white","#3366CC"), space="Lab")
plotcorr(ctab, col=rgb(colorfun((ctab+1)/2), maxColorValue=255),
         mar = c(0.1, 0.1, 0.1, 0.1))
```

Another view of bi-variate relationships
```{r}
pairs(xs, cex=1.5)		 
```

The following a closer look at potential predictor variables, starting with how to treate "race".  
```{r}
# What about race?     
(r.tab <- table(nels$race))
tot <- sum(r.tab)
w <- r.tab[4]
```

Since `r r.tab[4]` out of `r tot` are whites, I decided to re-code this as a dichotomous variable.
```{r}
# Dichotomize race into white/not-white
nels$white <- ifelse(nels$race==4, 1, 0)
nels$white <- as.factor(nels$white)
```

We'll treat gender as a categorical/factor variable
```{r}
nels$gender <- as.factor(nels$sex)
```

# Ordinary Least Squares Regression

Just for comparison, I fit a couple of OLS regression models that turn out to be the more interesting ones. One with all predictors and one with a subject of them.

```{r}
# All predictors
nels$gender <- as.factor(nels$sex)
summary(lm(math~gender + ses + paredu + homework + white, data=nels))
```

Ones I think will be important
```{r}
summary(lm(math~ paredu + homework, data=nels))
```

# Bayesian Model 1:  all predictors

## The data
The data list with all the predictors 
```{r}
dataList <- list(
         y     =nels$math,
				 pared =nels$pared,
				 hmwk  =nels$homework,
         ses   =nels$ses,
				 gender=nels$gender,
				 white = nels$white,
                 N=length(nels$math),
                 sdY = sd(nels$math)
                 )
```

# The Model
The model with all predictors is basically the same as simple linear regression, but there are just more variables.  Rather than just use b1, b2, etc, I added more descriptive names to the coefficients.
```{r}
model1 = "model { for (i in 1:N){
          y[i] ~ dnorm(mu[i] , precision)
          mu[i] <- b0  + bpared*pared[i] + bhmwk*hmwk[i] + bses*ses[i]+ bsex*gender[i]
                   + bwhite*white[i]        
          }
          b0 ~ dnorm(0 , 1/(100*sdY^2) )   
          bpared ~ dnorm(0 , 1/(100*sdY^2) )   
          bhmwk ~ dnorm(0 , 1/(100*sdY^2) )   
          bses ~ dnorm(0 , 1/(100*sdY^2) )   
          bsex ~ dnorm(0 , 1/(100*sdY^2) )   
          bwhite ~ dnorm(0 , 1/(100*sdY^2) )   
          sigma ~ dunif( 1E-3, 1E+30 )
          precision <- 1/sigma^2           
          } "

writeLines(model1, con="model1.txt")
```

## Initial Values (Starting values)
Need some starting values --- on set per chain that I plan to run.  I used sample mean and standard deviation for one, but used random draws from normal or uniform for rest of them.
```{r}
initsList = list(
list("b0"=mean(nels$math), "bpared"=0, "bhmwk"=0, "bses"=0, "bsex"=0, "bwhite"=0, "sigma"=sd(nels$math)),
list("b0"=rnorm(1,50,5),"bpared"=rnorm(1,-2,1),"bhmwk"=rnorm(1,2,1),"bses"=rnorm(1,0,1), "bsex"=rnorm(1,1,1),  "bwhite"=rnorm(1,0,1),"sigma"=sd(nels$math)), 
list("b0"=rnorm(1,50,5), "bpared"=rnorm(1,1,1),"bhmwk"=rnorm(1,1,1), "bses"=rnorm(1,1,1), 	"bsex"=rnorm(1,1,1), "bwhite"=rnorm(1,1,1), "sigma"=runif(1,0,10)), 
list("b0"=rnorm(1,50,6), "bpared"=rnorm(1,0,1),"bhmwk"=rnorm(1,0,1), "bses"=rnorm(1,0,1), "bsex"=rnorm(1,0,1),"bwhite"=rnorm(1,0,1),"sigma"=runif(1,0,10)) 
 )
```

## Get some samples
Run simulations of drawing samples from posterior using runjags
```{r}
mod1.runjags <- run.jags(model=model1, 
    monitor=c("b0","bpared","bhmwk","bses","bsex","bwhite","sigma","dic"),
    data=dataList, 
		n.chains=4,
		inits=initsList)
```

No error messages but should check convergence.  

## Convergence checks 
```{r fig.show="hide", results="hide", echo=TRUE}
plot(mod1.runjags)
```

Without doing anything else there are problems:  

1.     Chains for all parameters except sigma are not mixing well
2.     High auto-correlations for b0, bpared, bses, bwhite

Before trying to get better convergence, we will look at the results, which has further diagnostic information.
```{r}
print(mod1.runjags)
```

## Adding Iterations 

We can draw more samples; however, rather than re-run what we already have, you can use the following.  Note that we don't need to burnin because we alreay done this.  Since I bumped up the number of iterations to 200,000, it will take some time (it takes 30 seconds on my desktop).  
```{r eval=FALSE}
mod1.extend <- extend.jags(mod1.runjags,  
                           burnin=0, 
						               sample=200000)
plot(mod1.extend)
print(mod1.extend)      
```

Mixing is better, but still have large auto-correlations.

## Thinning

When auto-correlations are large, thinning will often help. Only every 10th sample will be kept.  This will increase computation time.
```{r eval=FALSE}
mod1.thin10 <- run.jags(model=model1, 
            monitor=c("b0","bpared","bhmwk","bses","bsex","bwhite","sigma","dic"),
            data=dataList, 
						n.chains=4, 
						inits=initsList,
						sample=100000,
						thin=10)

plot(mod1.thin10)

gelman.plot(mod1.thin10)
```


print(mod1.thin10)
```

Much better.

# Auto-run
When you don't want to think about what you're doing and want to let the software make more decisions, you can use autorun
```{r} 
mod1b.runjags <- autorun.jags(model=model1, 
                 monitor=c("b0","bpared","bhmwk","bses","bsex","bwhite","sigma","dic"),
                 data=dataList, 
							   n.chains=4, 
							   inits=initsList,
							   thin=10)

plot(mod1b.runjags)

print(mod1b.runjags)
```

# Model 2:  Refinement 

I will drop non-significant predictors, which means change data list, model, starting values, and call to run jags.
```{r}
dataList <- list(
         y     =nels$math,
				 pared =nels$paredu,
				 hmwk  =nels$homework,
         N=length(nels$math),
         sdY = sd(nels$math)
         )
model2 = "model {  for (i in 1:N){     
                       y[i] ~ dnorm(mu[i] , precision)
                       mu[i] <- b0  + bpared*pared[i] + bhmwk*hmwk[i]     
 					          }
                     b0 ~ dnorm(0 , 1/(100*sdY^2) )   
                     bpared ~ dnorm(0 , 1/(100*sdY^2) )   
                     bhmwk ~ dnorm(0 , 1/(100*sdY^2) )    
                     sigma ~ dunif( 1E-3, 1E+30 )
                     precision <- 1/sigma^2           
              } "
writeLines(model2, con="model2.txt")

initsList = list(list("b0"=mean(nels$math), "bpared"=0, "bhmwk"=0,            "sigma"=sd(nels$math)), 
                 list("b0"=rnorm(1,50,5),   "bpared"=rnorm(1,-2,1), "bhmwk"=rnorm(1,2,1), "sigma"=sd(nels$math)), 
                 list("b0"=rnorm(1,50,5),   "bpared"=rnorm(1,1,1), "bhmwk"=rnorm(1,1,1),  "sigma"=runif(1,0,10)), 
                 list("b0"=rnorm(1,50,6),   "bpared"=rnorm(1,0,1), "bhmwk"=rnorm(1,0,1),  "sigma"=runif(1,0,10)) 
                 )

mod2.runjags <- run.jags(model=model2, 
                monitor=c("b0","bpared","bhmwk","sigma","dic"),
                data=dataList, 
						    n.chains=4, 
						    inits=initsList)
plot(mod2.runjags)

print(mod2.runjags)
```

This model look pretty good. Also according to DIC statistics (smaller is better), the simpler model is preferred.

But could use some thinning\ldots
```{r}
mod2.runjags.thin5 <- run.jags(model=model2, 
                monitor=c("b0","bpared","bhmwk","sigma","dic"),
                data=dataList, 
						    n.chains=4, 
						    inits=initsList,
						    thin=5)
						    
plot(mod2.runjags.thin5)

print(mod2.runjags.thin5)
```

# Model 3:  Interactions

Since adding interactions, often leads to multicolinearity, I centered the predictors between adding the interaction.
```{r}
nels$chomework <- nels$homework -mean(nels$homework)
nels$cparedu  <-  nels$paredu  -mean(nels$paredu )
```

Below is the rest of the code for this model.
```{r}
dataList <- list(
          y     =nels$math,
				 cpared =nels$cparedu,
				 chmwk  =nels$chomework,
         N      =length(nels$math),
         sdY    = sd(nels$math)
                 )
```

```{r}
model3 = "model {  for (i in 1:N){
                       y[i] ~ dnorm(mu[i] , precision)
                       mu[i] <- b0  + bpared*cpared[i] + bhmwk*chmwk[i]  
                                + bintact*cpared[i]*chmwk[i]    
                     }
                     b0 ~ dnorm(0 , 1/(100*sdY^2) )   
                     bpared ~ dnorm(0 , 1/(100*sdY^2) )   
                     bhmwk ~ dnorm(0 , 1/(100*sdY^2) )    
                     bintact ~ dnorm(0 , 1/(100*sdY^2) )    
                     sigma ~ dunif( 1E-3, 1E+30 )
                     precision <- 1/sigma^2           
              } "
writeLines(model3, con="model3.txt")
```

```{r}
initsList = list(list("b0"=mean(nels$math), "bpared"=0, "bhmwk"=0, "bintact"=0, "sigma"=sd(nels$math)), 
                 list("b0"=rnorm(1,50,5),   "bpared"=rnorm(1,-2,1), "bhmwk"=rnorm(1,2,1), "bintact"=.3, "sigma"=sd(nels$math)), 
                 list("b0"=rnorm(1,50,5),   "bpared"=rnorm(1,1,1), "bhmwk"=rnorm(1,1,1),  "bintact"=.5,  "sigma"=runif(1,0,10)), 
                 list("b0"=rnorm(1,50,6),   "bpared"=rnorm(1,0,1), "bhmwk"=rnorm(1,0,1),  "bintact"=-.3,  "sigma"=runif(1,0,10)) 
                 )
```

```{r}
mod3.runjags <- run.jags(model=model3,
                monitor=c("b0","bpared","bhmwk","bintact","sigma","dic"),
                data=dataList, 
						    n.chains=4, 
						    thin=5, 
						    inits=initsList)
```

A quick convergence check looks good:
```{r eval=FALSE}
plot(mod3.runjags)
print(mod3.runjags)
```

Since the 95% interval for the interaction include 0, so we will delete it.  This takes us back to model 2; however, in what follows I will go back to un-centered variables.  

# Model 4:  Final model with various derived statistics

I will add some values that I want computed and some derived variables.  For some of the analyses done here, I find it more convenient to use rjags than runjags, so I make the switch here

Since I'm switching back to UNCENTERED I 
```{r}
dataList <- list(
         y     =nels$math,
				 pared =nels$paredu,
				 hmwk  =nels$homework,
         N     =length(nels$math),
         sdY   = sd(nels$math)
                 )

nels_final = "model {  for (i in 1:N){
            y[i] ~ dnorm(mu[i] , precision)
            mu[i] <- b0  + bpared*pared[i] + bhmwk*hmwk[i]    

            res[i] <- y[i] - mu[i]               # Data residuals
            emp.new[i] ~ dnorm(mu[i],precision)  # New Predicted y 
            res.new[i] <- emp.new[i] - mu[i]     # New Predicted residuals
            }
            b0 ~ dnorm(0 , 1/(100*sdY^2) )   
            bpared ~ dnorm(0 , 1/(100*sdY^2) )   
            bhmwk ~ dnorm(0 , 1/(100*sdY^2) )    
            sigma ~ dunif( 1E-3, 1E+30 )
            precision <- 1/sigma^2 
# Derived quantities			
            fit <- sum(res[])                    # sum data residuals this iterations
            fit.new <- sum(res.new[])            # sum predicted residuals 
            fit.mean <- mean(emp.new[])          # mean of posterior predictions
              } "
writeLines(nels_final, con="nels_final.txt")
```
I added res[ ] as raw residuals, emp.new[ ] are predicted y, and res.new[ ] are the residuals of these new predictions.  The predicted y will have the same values of cpared[ ] and chmwk[ ].  

I used the same initial or starting values as above, so we don't need to re-run this here.  Which brings us to the simulation itself.
```{r}
b0Init = mean(nels$math)
bparedInit = 0
sigmaInit = sd(nels$math)
initsList = list(b0=b0Init, bpared=bparedInit, sigma=sigmaInit )

nels_final <- jags.model(file="nels_final.txt", 
                        data=dataList,
                        inits=initsList,
                        n.chains=4,
                        n.adapt=1000)           

update (nels_final, n.iter=1000)             

# gets samples from all chains for 2000 iterations
Samples <- coda.samples(nels_final, 
             variable.names=c("b0","bpared","bhmwk", "sigma", "fit", "fit.new",
                          "fit.mean", "emp.new"), 
             thin=5,
			       n.iter=10000)
```

A long summary table...
```{r}
summary(Samples)
```

When you get plots and print Samples, you will get something for every value in the variable names list. Note that this include `r n` values of fit.  This is one reason that I don't include all the these "extras" when working on model refinement.


To be able to use the samples, I need to do some data manipulation. In "Samples"", the data are in separate chains.  The object "Samples"" is an mcmc object which is difficult for me to manipulate, so I changed the object type.
```{r}
ffit <- as.array(Samples)                          # change from mcmc to array object                                 
b <- rbind(ffit[,,1],ffit[,,2],ffit[, , 3], ffit[, , 4])  # combine the chains
```
The following takes Samples and breaks out specific estimates into different objects.  To know what the "numbers" for the indices, I had to use what I know about the sample size of nels data and the number of other parameters (looking at the samples helps).
```{r}
b0 <- b[,1]
bpared <- b[,2]
bhmwk <- b[,3]
pred <- b[,4:70]     
fit <- b[,71]
fit.mean <- b[,72]
fit.new <- b[,73]
sigma <- b[,74]
```

A little checking to see if things make sense.  I expect the following two values to be very similar, because one is the mean of the data and the other is the mean of the postierior distribution of the fitted values.
```{r}
mean(nels$math)
mean(pred)
```

I wanted to look at the data and posterior distribution in different ways.  In the first case, I collapsed over students so that I have one fitted value for each replications.  This is like finding means for pred and nels$math, but doesn't collapse as much.
```{r}
fitted <- apply(pred,2,mean)
```
And some graphics
```{r}
par(mfrow=c(2,2))
hist(nels$math,breaks=10,main="Distribution of Data")
hist(fitted,breaks=10,main="Posterior Predictions (8000 iterations)")
R2 <- round(cor(nels$math, fitted)**2, digit=4)
plot(nels$math, fitted, 
     main=paste("Fitted by Observed Math  R2=",R2),
     xlab="Observed Math Scores",
     ylab="Mean over Samples of Fitted Math")
```

```{r}
plot(nels$homework, nels$math,main="Data Points")
plot(nels$homework, fitted, ylim=c(58,70))
# Get the means to add to plot
schMathHmwk <- aggregate(nels$math, list(nels$homework), mean)
lines(schMathHmwk$Group.1,schMathHmwk$x,type='p',pch=19,col="blue")
```

I know collapse over replications (which makes more sense), and get a single mean value for each student.
```{r}
# Compute statistics over rows (students)
ymeans <- apply(pred,1,"mean")
ymin  <- apply(pred,1,"min")
ymax  <- apply(pred,1,"max")
ysd   <- apply(pred,1,"sd")
```
And the Bayesian p-values which gives another way to look at how well the model is doing for different aspects of the data.
```{r}
pmean <- mean(mean(nels$math) < ymeans)
pmin  <- mean(min(nels$math) < ymin)
pmax  <- mean(max(nels$math) < ymax)
psd   <- mean(sd(nels$math) < ysd)
```
Understanding what these are is a bit easier by looking at them relative to data.
```{r}
par(mfrow=c(2,2))
 
hist(ymin,breaks=15,main=paste("Simulated within Jags: Minimums", "\nBayesian P-value =", round(pmin, 2)))
lines(c(min(nels$math),min(nels$math)),c(0,2000),col="blue",lwd=2)
2
hist(ymax,breaks=15,main=paste("Simulated within Jags: Maximums", "\nBayesian P-value =", round(pmax, 2)))
lines(c(max(nels$math),max(nels$math)),c(0,2000),col="blue",lwd=2)

hist(ymeans,breaks=15,main=paste("Simulated within Jags:  Means", "\nBayesian P-value =", round(pmean, 2)))
lines(c(mean(nels$math),mean(nels$math)),c(0,2000),col="blue",lwd=2)

hist(ysd,breaks=15,main=paste("Simulated within Jags:  SDs", "\nBayesian P-value =", round(psd, 2)))
lines(c(sd(nels$math),sd(nels$math)),c(0,2000),col="blue",lwd=2)
```

# Robust Regression:  Student's t-distribution

We swap out "dnorm" by "dt", add priors for degrees of freedom (nu), and add nu to list of parameters to monitor.  Only the parts of the code that change are given below.  The model,

```{r}
dataList <- list(
         y     =nels$math,
				 pared =nels$paredu,
				 hmwk  =nels$homework,
         N     =length(nels$math),
         sdY   =sd(nels$math)
                 )
```

```{r}
ymean <- mean(nels$math) 
sdy <- sd(nels$math)

initsList = list(list("b0"=ymean, "bpared"=0, "bhmwk"=0, "sigma"=sdy), 
list("b0"=rnorm(1,50,5), "bpared"=rnorm(1,-2,1), "bhmwk"=rnorm(1,2,1), "sigma"=sdy), 
 list("b0"=rnorm(1,50,5),   "bpared"=rnorm(1,1,1), "bhmwk"=rnorm(1,1,1), "sigma"=runif(1,0,10)),
list("b0"=rnorm(1,50,6),   "bpared"=rnorm(1,0,1), "b"=rnorm(1,0,1),  "sigma"=runif(1,0,10)) )
```


```{r}
tmr = "model { for (i in 1:N){
# Likelihood          
              y[i] ~ dt(mu[i] , precision, nu)
               mu[i] <- b0  + bpared*pared[i] + bhmwk*hmwk[i]   
             }
            b0 ~ dnorm(0 , 1/(100*sdY^2) )   
            bpared ~ dnorm(0 , 1/(100*sdY^2) )   
            bhmwk ~ dnorm(0 , 1/(100*sdY^2) )    
            sigma ~ dunif( 1E-3, 1E+30 )
            precision <- 1/sigma^2  
            nuMinusOne  ~ dexp(1/29)  
            nu <- nuMinusOne+1  					 
           }"

writeLines(tmr, con="tmr.txt")
```

and the call for samples,
```{r}
tmr.runjags <- run.jags(model=tmr, 
               monitor=c("b0","bpared","bhmwk","sigma","nu","dic"),
               data=dataList, 
               n.chains=4, 
               thin=5, 
               inits=initsList)
```
```{r}
plot(tmr.runjags)

gelman.diag(tmr.runjags)

gelman.plot(tmr.runjags)

cumuplot(tmr.runjags,probs=c(.25,.50,.75),lwd=c(1,2),lty=c(2,1),col=c("blue","red"))

print(tmr.runjags)


```

# Bayes factor

If you feel compelled to compute these, it's pretty if you use the BayesFactor package. This package does more then just linear regressions.   See http://bayesfactorpcl.r-forge.r-project.org/ for more information about the package.

You basically, ask for a regression using the regressionBF command with (in this case) all our possible variables.  The following code gives Bayes factor for all possible models where the comparison model (default) is model with only an intercept.
```{r}
library(BayesFactor)
nels$xwhite <- as.numeric(nels$white)
(bf <- regressionBF(math ~ paredu + homework + ses + xwhite + sex, data=nels))
```

If you have lots and lots of possible models, you may want just the top ones (again compared against model with only an intercept), the add the line
```{r}
#The best 5
head(bf,n=5)
```

or compare the best with the top competitors and put the best model in the denominator.
```{r}
# Compare the best with top 5
(top_models <- head(bf)/max(bf) )
```


# Speeding up convergence

As models become more complex and with large data sets, the time it take for an algorithm to converge can be an issue.  Possible step to recude computational time incude

*     Turn of real time progress indicator off.  
*     Use parallel procressing which can be done in jagsUI using the parallel=TRUE option.  This is illustrated below.  It does take a little time to set up parallel but for large problems, the cpu savings is likely to be worth it.
*    In terms of wrappers (serial processing), from me mini study the fastest is jagsUI and rjags and the slowest was runjags.  There is a "parallel.jags" comand in runJags, but it is buggy and I couldn't get it to run (trouble with how many intital values are needed).
*     BLAS matters (Basic Linear ALgebra Subprograms).  For PC, the easiest thing to do is install Microsoft Open R which uses the Intel MKL libraries for linear algebra.  For Apple, you might be able to install for example OpenBLAS or ATLAS.  For more one this see https://csantill.github.io/RPerformanceWBLAS.
*     User faster computer or cluster computing.
*     Write you underlying C, but this is more than we'll get into.


## Parallel Processing with GIBSS


Times vary from run to run.  Expect to see bigger differences 

*     Expect to see bigger differences with more iterations
*     Expect to see bigger differences more parameters (complex models)
*     More chains (should be <= number of cores (I didn't check but                           think this should be so).
*     Larger data sets


## Set up 

```{r}
dataList <- list(y=nels$math,
                 x=nels$homework,
                 x2=nels$pared,
                 N=length(nels$math),
                 sdY = sd(nels$math)
                 )
```
```{r}
nelsLR1 = "model {                                    
                 for (i in 1:N){
                  y[i] ~ dnorm(mu[i] , precision)
                  mu[i] <- b0 + b1*x[i] +b2*x2[i]        
                 }
                 b0 ~ dnorm(0 , 1/(100*sdY^2) )   
                 b1 ~ dnorm(0 , 1/(100*sdY^2) )   
                 b2 ~ dnorm(0 , 1/(100*sdY^2) )   
                 precision ~ dgamma(.01, .01)
               } 
"
writeLines(nelsLR1, con="nelsLR1.txt")
```

## rjags
```{r}
initsList <- list("b0"=mean(nels$math), "b1"=sd(nels$math), "b2"=sd(nels$pared), "precision"=1/sd(nels$math)**2) 

start_time <- Sys.time()
jagsNelsLR1 <- jags.model(file="nelsLR1.txt", 
                        data=dataList,
                        inits=initsList,
                        n.chains=4,
                        n.adapt=4000)   
						
Samples1 <- coda.samples(jagsNelsLR1, variable.names=c("b0", "b1", "b2", "precision"), n.iter=10000)
end_time <- Sys.time()
end_time - start_time
```

## runjags

```{r}
initsList = list(list("b0"=mean(nels$math), "b1"=sd(nels$math), "b2"=sd(nels$pared), "precision"=1/sd(nels$math)**2), 
                 list("b0"=rnorm(1,-2,4), "b1"=rnorm(1,0,4), "b2"=rnorm(1,0,4), "precision"=1/2),
                 list("b0"=rnorm(1,0,1),  "b1"=rnorm(1,6,4),"b2"=rnorm(1,0,4), "precision"=1/5),
                 list("b0"=rnorm(1,2,2),  "b1"=rnorm(1,-2,4), "b2"=rnorm(1,2,4) ,"precision"=1/1515 )   
                 )
```
```{r}
start_time <- Sys.time()
out.runjags <- run.jags(model=nelsLR1, 
              monitor=c("b0","b1","precision","dic"),
              data=dataList, 
					  	n.chains=4, 
						  inits=initsList)
end_time <- Sys.time()
end_time - start_time
```

## jagsUI -- serial 

```{r}
start_time <- Sys.time()				 
out.jagsUI <- jags(model.file="nelsLR1.txt",
              data=dataList,
			        inits=initsList, 
              parameters.to.save=c("b0", "b1", "precision"), 
              n.iter=10000,
              n.burnin=4000,
              n.chains=4,
              verbose=FALSE)
end_time <- Sys.time()
end_time - start_time
```

## parallel via jagsUI

```{r}
start_time <- Sys.time()				 
parallel.jagsUI <- jags(model.file="nelsLR1.txt",
              data=dataList,
			        inits=initsList, 
              parameters.to.save=c("b0", "b1", "precision"), 
              n.iter=10000,
              n.burnin=4000,
              n.chains=4,
	            parallel=TRUE,
	     )
end_time <- Sys.time()
end_time - start_time
```

## Output from jagsUI

```{r}
whiskerplot(parallel.jagsUI, parameters=c("b0", "b1", "precision"),
            quantiles=c(0.025,0.975), zeroline=TRUE)
```

To use coda plotting functions input the samples
For example

```{r}
gelman.diag(parallel.jagsUI$samples)
gelman.plot(parallel.jagsUI$samples)

geweke.diag(parallel.jagsUI$samples, frac1=.10, frac2=.50)

HPDinterval(parallel.jagsUI$samples)
				 
autocorr.plot(parallel.jagsUI$samples)
```