---
title: 'Gibbs Sampling using jags: Mean and variance of Anorexia data'
author: "C.J. Andersion"
date: "9/25/2022"
output: 
   html_document:
       toc: yes
---

In this session, we'll use rjags, runjags and jagsUI to estimate the mean and variance of a weight change of girls being treated for Anorexia. Before you can run this you need to 

1.     Install the program jags
2.     Install the r packages:  coda, rjags, runjags, and jagsUI.

There is at least one other r package that is a wrapper for jags, but these are more than enough (and the ones illustrated below)


## Set-up

Below are the four packages that we'll be using.  The coda has fucntions that compute/create various diagnostics. The multiple jags packages have slightly different input and output, but all are wrappers for the program jags.
```{r}
library(coda)
library(rjags)
library(runjags)
library(jagsUI)
```

As always, our first step is to set your working directory, read in the data, create a change in weight score, and take a look at the data:

```{r}
setwd("C:/Users/cja/Dropbox/edps 590BAY/Lectures/5 Gibbs Sampling")
ano <- read.table("anorexia_data.txt",header=T)
ano$change <- ano$weight2 - ano$weight1
head(ano)
```

I would like to get the same results each time I run this code, so I'll set a seed for the random number generaters.

```{r}
set.seed(234590)
```

## Create data list

The data list will be the same for rjags, runjags and jagsUI.  We list and define variables that will be used to create our data model (i.e., the likelihood). For our problem this is 

```{r}
dataList <- list(y=ano$change, 
                 Ntotal=length(ano$change),
                 meanY = mean(ano$change),
                 sdY = sd(ano$change)
                 )
```

The dataLlist object now contains the following
```{r}
dataList
```

## Create the data model 

The model will define the likelihood and also give the priors for the parameters being estimated.  Note that jags "dnorm" is different from base "dnorm".  For jags, the parameters for jags "dnorm" are the mean and preicsion, which is 1/variance of the data.

```{r}
Model1 = "model {                                    
    for (i in 1:Ntotal){
       y[i] ~ dnorm( mu, precision )                
        }
    mu ~ dnorm( meanY , 1/(100*sdY^2) )  # mean (y) and precision N(2.763889,1.568927e-06)
    precision <- pow(sigma, -2)         # tau is presicion = 1/var(y)
    sigma ~ dunif( sdY/1000, sdY*1000 )  # variance (y) and uniform(.00798,7983.598)        
   }
   "
```

Note the quotes around this code.  We will save the model to a text file that will be called by jags.  The following command writes the model to a file that I'm calling "Model1.txt".  You can find this on your hard drive in the current working directory that you set.

```{r}
writeLines(Model1, con="Model1.txt")
```
Go to your working memory and look for "Model1.txt" and take a look at it.



## Set Initial values (rjags)

Regardless of which implemenations you use for jags, they all require that intital values are a list object.   In this case, I simpley used the same mean and variance.

```{r}
thetaInit = mean(ano$change)
sigmaInit = sd(ano$change)
initsList = list(mu=thetaInit,sigma=sigmaInit)
```


## Run rjags

```{r}
jagsModel1 <- jags.model(file="Model1.txt",     # compiles and intializes model
                        data=dataList,
                        inits=initsList,
                        n.chains=4,
                        n.adapt=500)            # adaptive phase for 500 iterations
```

Model model has successfully complieds; that is, it has be initialized and we can do some sampling.

```{r}
update (jagsModel1, n.iter=500)                 # gets samples from posterior with a 
                                                #  burn in of 500 iterations
``` 
And now for getting the samples for many iterations.

```{r}
# contains samples from all chains with 500 iterations
Samples <- coda.samples(jagsModel1, variable.names=c("mu"," precision","sigma"), n.iter=4000)
```

## Did our model converge?

There are many diagnostics that we can look at including densities, trace plots, geweke statistics, Rhat/psrf/Gelman statistics, auto-correlations, effective sample sizes and more.

We will start with trace plots and densities for each of our parameters that we've estimated.  Note that these look better (nicer formatting) when I run them in R.

```{r}
plot(Samples)
```

Although we can look at the trace plots for stationary and nice mixing of the chains, to back what we think we see we can get the Gelman (proportional scale reduction factors or Rhats).  These should be less than .01.

```{r}
gelman.diag(Samples)
```
Can't get much better than these numbers.  We can also get a plot of them
```{r}
gelman.plot(Samples)
```

Below are dome other diagnostics as can look at and these done for each chain for each parameter.
```{r}
autocorr.plot(Samples,auto.layout=TRUE)
geweke.diag(Samples,frac1=0.1,frac2=0.5)
cumuplot(Samples,probs=c(.25,.50,.75),lwd=c(1,2),lty=c(2,1),col=c("blue","red"))
```


Lastly, we might want to know what the effective sample size is for each of our parameters in the posterior distribution.
```{r}
effectiveSize(Samples)
```
Recall that the effective sample size is the number of independent Monte Carlo samples needed to give the same precision as the MCMC samples. 

## Tables of Summary Statistics from our Posterior

Dropping the burn-ins or warm-ups, we can compute various sample statistics that describe our posterior distributions.  Note that these are computed over all chains.

```{r}
summary(Samples)
```

# Working With Samples: 

If we put our posterior sample into a more useable format, we can compute any statistic that we want to need.   This includes what is in the summary table as well as anything else that you might want to know.

The main reason for doing this is to show you how the three measures of variability are computed.  For these I will give formulas but also show what they are with "Samples"; that is, how they are obtained.

## Converting "Samples" into Data Frame

Let's first looks are the object Samples
```{r}
class(Samples)

#-- What's in First chain
dim(Samples[[1]])
class(Samples[[1]])
head(Samples[[1]])
```

We need to vertically concatenate the samples from each chain into a single matrix or array.   The following will do this
```{r concatenatw}
post_samples <- as.data.frame(rbind(Samples[[1]], Samples[[2]],
                      Samples[[3]], Samples[[4]]))

class(post_samples)
dim(post_samples)
head(post_samples)
```

## Sample statistics of posterior

### Mean
Compare the following with what is in summary table

```{r xbar}
(Ntotal <- dim(post_samples)[1])   #--- number iterations X number chains
(Nchains <- 4)
(Niter.per.chain <- 4000)          

#--- These are means are same as those above
(mu.mean <- colMeans(post_samples)) 
```


## Variability Measures

### SD
The column label "SD" in the summary tables for "Samples" is usualy way we would compute the standard deviation of a sample of data:
$$ SD = \sum_{i=1}^I (\hat{\theta}_i - \bar{\hat{\theta}})^2 / N, $$
where the summation is over all iterations and chains.

The R to get SDs,
```{r}
(sd.samples <- sapply(post_samples, sd))
```

We can also compute the variance, which we'll be using in the next measures
```{r}
# --- for below 
(var.samples <- sapply(post_samples, var))
```

### SE of MCMC Estimates (i.e. MCerr)

This is given in the run.jags output. We might want a measure of the precision (or variance) of the MCMC samples.  This is the SE of MCMC estimates but uses a different sample size; namely,
\[ \mbox{var}_{mcmc(\theta)} = \frac{\mbox{var}(\theta)}{\mbox{effective sample size}} \]

and in R
```{r}
(eff.size <- effectiveSize(Samples))
(var.mcmc <- var.samples/eff.size)
(MCerr <- sqrt(var.mcmc))
```


### Naive SE
The next two sample statistics is based on information in https://stats.stackexchange.com/questions/74450/naive-se-vs-time-series-se-which-statistics-should-i-report-after-bayesian-esti

The MCerr (above), the naive se, and the times series are measures of possible computational errors of the MCMC's posterior **expected values**.  

The naive SE are standard errors of expected values similar to the standard deviation of (hypothetical) sampling distributions in frequentest frame work (but ours are based on data).  

The naive standard error (se.naive) ignores that possible auto-correlation (dependency) of sample estimates; therefore, it is not a good estimator of the of variability of the expected values (i.e., the "Means"). The naive SE is even worse than the MCerr.  The MCerr takes into account the effective sample size, which in turn is related to auto-correlation.  If the auto-correlation is large, then effective sample size is small and naieve SE is too small.

The naive values are computed as
$$ \mbox{se.naive } = \sqrt{\frac{\mbox{var}(\theta)}{N}},$$
were $N$ equals the number of chains times iterations.

```{r naive}
(se.naive <- sqrt(var.samples/Ntotal))
```

Note that these are smaller for precisio and sigma, but are the same for mu (the effective sample size for mu is 16,000 --- the sames as $N$).

### Time-series SE

The time-series estimate corrects for possible auto-correlation. An auto-correlation is for specific "lags is the correlation between successive values.

An auto-correlation is the correlation between a set of numbers (ordered over time or space) and same set but shifted.  For example, the correlation between the first two columns below is a lag 1 correlation, the correlation between the first and third column is a lag 2 correlation, and so on.

$$ \begin{array}{rrr}
\theta_1 &  -- & --       \\
\theta_2 & \theta_1 & --   \\
\theta_3 & \theta_2 & \theta_1 \\
\theta_4 & \theta_3 & \theta_2 \\
\vdots   &      &         \\
\theta_s & \theta_{2-1} & \theta_{s-2} \\
\end{array} $$

The time-series SE involves fitting auto-correlation models.  The definition of time-series SE is 
\[  \mbox{SE}_{ts} = \sqrt{\frac{\mbox{var}_{ts}(\theta)}{\mbox{Ntotal}}}, \]
where $\mbox{var}_{ts}(\theta)$ is the mean of $\mbox{var}_{ts}(\theta^(c))^(c))$ which are obtained by fitting an auto-regressive models to each chain c.  The auto-regressive model that is fit is 
\[ \theta_t^(c) - \mu^(c) = \sum_{k=1}^K \rho_k(\theta_k^(c) - \mu^(c)) +\epsilon,\]
where $\epsilon \sim (N(0,\sigma^2_{\epsilon})$, K is the order or largest lag needed, $\theta_t^{(c)}$ is the sampled value from posterior for chain c.  Using this we can compute
\[ \mbox{var}_{ts}(\theta^{(c)})^{(c)}= \frac{\sigma^2_{\epsilon}}{(1- \sum_{k=1}^{K} \rho_k)^2 }\]

Since this has to be done repeatably (i.e., for each chain and each parameter), I wrote function.  This is was my way to make sure that this gives the same results as in the summary table. 



```{r}
#########################################
# function that computes time-series SE #
#########################################
# parameter = 1 first parameter
#             2 for 2nd parameter
#			        ...
#			        P for Pth parameter	
#
# Samples is an mcmc.list 
#

se.ts <- function(Samples, parameter) {
	
	Nchains <- length(Samples)
 	Ntotal <- nrow(Samples[[1]])*Nchains
	
  var.tsx <- matrix(NA, nrow=Nchains, ncol=1)

  for (c in 1:Nchains) {
	  chain.c <- Samples[[c]]
      fit.ar <- ar(chain.c[ ,parameter], order.max=10, aic=TRUE)
	  var.tsx[c,1] <- fit.ar$var.pred/((1-sum(fit.ar$ar))**2)
  }
	var.tsc<- colMeans(var.tsx)
	var.ts <- var.tsc/Ntotal
	sd.ts <- sqrt(var.ts)
	return(sd.ts)
}
```
For mu, precision and sigma

```{r}
se.ts(Samples, 1)
se.ts(Samples, 2)
se.ts(Samples, 3)
```


# runjags run!


Different front ends onto jags are a bit different in terms of default output and input.  All these use the same datalist and model but differ in terms of initial or starting values and getting the samples. First we'll look at runjag and then jagsUI.  The latter is overly verbose.

So to we set up out starting values.  I'll use the sample mean and standard deviation of change values for one of the chains and for the others use different means and variances. For the other 3 chains, I sampled from normal distributions with different means and standard deviations.

```{r}
(meanY = mean(ano$change))
(sdY = sd(ano$change))
initsList = list(list("mu"=meanY, "sigma"=sdY),
                 list("mu"=rnorm(1,2,4), "sigma"=1 ),
                 list("mu"=rnorm(1,4,1), "sigma"=8 ),
                 list("mu"=rnorm(1,-4,2),"sigma"=.5 )
                 )
```
Go ahead and take a look at this.

The next comand gets the samples (i.e., runs jags)
```{r}
out.runjags <- run.jags(model=Model1, 
                        monitor=c("mu","sigma","precision","dic"),
                        data=dataList, 
                        n.chains=4, 
                        inits=initsList)
```
runjag will easly product diagnostic plots, for example
```{r}
plot(out.runjags)
```

You can also use coda to produce the same graphics as I did with rjags.

The table of summary statistics from runjag includes more than rjags gives.
```{r}
print(out.runjags)
```


# jagsUI

If you only use jags once in a while and don't remember what all this various statistics are, the jagsUI is for you.  We use the same data list, model, and intial values list as we did for runjags, but a little different command to get the samples.

```{r}
out.jagsUI <- jags(model.file="Model1.txt",
              data=dataList,inits=initsList, 
              parameters.to.save=c("mu","sigma","precision"), 
              n.iter=2000,
              n.burnin=500,
              n.chains=4)
```
and the verbose table of summary statistics
```{r}
print(out.jagsUI)
```