---
title: 'Multilevel: Anorexia Data'
author: "Carolyn J. Anderson"
date: "10.14.2022"
output: 
  html_document:
    toc: yes
---

Below is code for a number of different models.  The results of some but not all models are run below.  You can do these on your own.

There is code for obtaining estimate of the posterior distribution of the icc (intra-class correlation), as well as the posterior probabilities for specific hypotheses. I also put in a better way to extract posterior values of parameter from samples (compared to what I did first time I did this).

# Setup
The packages  that we'll use are 
```{r setup, message=FALSE, warning=FALSE, include=TRUE}
library(lme4)
library(lmerTest)
library(rjags)
library(runjags)
library(coda)
library(texreg)
```
Set working directory to where the data live and set up data for initial analyses:
```{r}
setwd("D:/Dropbox/edps 590BAY/Lectures/8 Multilevel models")
ano.wide <- read.table("anorexia_wide.txt",header=TRUE)
ano.wide$change <- ano.wide$weight2 - ano.wide$weight1
(n <- nrow(ano.wide))
head(ano.wide)
```
# Examining the data
Let's do some graphics to look at the data.  First we'll just look at histograms of before treatment, after treatment, and change.  We'll add a 2nd graph of change but looking at is a bit differently.  For the 2nd graph of change we need to re-format the data, and create treatment dummy codes
```{r}
ano <- read.table("anorexia_long.txt",header=TRUE)
ano <- ano[order(ano$girl),]
ano$Rx1 <- ifelse(ano$Rx==1,1,0)
ano$Rx2 <- ifelse(ano$Rx==2,1,0)
ano$Rx3 <- ifelse(ano$Rx==3,1,0)
```
and now for the graphs
```{r}
par(mfrow=c(2,2))
hist(ano.wide$weight1, main='Weight Before Treatment',col="darkseagreen1")
hist(ano.wide$weight2, main='Weight After Treatment',col="darkseagreen1")
hist(ano.wide$change,  main='Change in Weight',col="darkseagreen1")
plot(ano$time,ano$weight,
     type="n", 
     ylim=c(65,104),
	 xlab="Before/After Treatment",
	 ylab="Weight in Pounds",
	 main="Change for Each Girl"
	 )
for (i in 1:n){
   g <- subset(ano,ano$girl==i)
   lines(g$time,g$weight,col=i,lwd=2)
}
```

This 4th plot show different intercepts and different slopes. To better see what's going on, we'll compute means for each treatment before and after treatment

```{r}
# means for each treatment 
byRx <- aggregate(ano$weight,list(ano$Rx),FUN="mean")

# mean for treatment x time
rx_by_time <- aggregate(ano$weight,list(ano$Rx,ano$time),FUN="mean")
names(rx_by_time) <-c("rx","time","weight")
rx1 <- subset(rx_by_time,rx==1)                 # Cognitive
rx2 <- subset(rx_by_time,rx==2)                 # Control 
rx3 <- subset(rx_by_time,rx==3)                 # Family
```
Now to look at these...what do you expect based on this?
```{r}
par(mfrow=c(1,1))
plot(rx2$time,rx2$weight,
     type='b',
     col="cyan",
     main="Weight Change by Treatment",
     ylim=c(80,95),
	 xlim=c(1,2),
	 xlab="When Weight Measured",
	 ylab="Weight",
	 xaxt="n"
	 	 )
axis(1, at = seq(1,2,by=1),labels=c("Before","After"))		 
lines(rx1$time,rx1$weight,type="b",col="blue")
lines(rx3$time,rx3$weight,type="b",col="red")
legend(1.02,95,title="Treatment",legend=c("Control", "Cognitive", "Family"), col=c("cyan","blue","red"),lty=c(1,1,1),cex=0.8)
```

# Maximum Likelihood Estimation

If we are taking a Bayesian approach, we would not look at MLE; however, for pedagogical purposes, we do it here.  Before using Bayesian estimation, we'll see what restricted maximum likelihood (REML) yields. I choose REML since these are unbiased; however, if we're only looking at lmer models from package lme4, I would normally choose MLE.   Using the lme4 package.
```{r  results=FALSE}
########################################################
# Simple Model using lme4 --- just a random intercept   #
########################################################
# Remember: treatment = 1 Cognitive
#                     = 2 Control
#                     = 3 Family
summary(model1reml <- lmer(weight ~ time + (1 | girl), data=ano, REML=TRUE))
```
Here are the fitted values from model 0
```{r}
girl <- data.frame(seq(1:n))
par(mfrow=c(1,1))

b1 <- data.frame(coef(model1reml)[[1]])
tmp1 <- cbind(b1,girl)
names(tmp1) <- c("b0","b1","girl")
ano <- merge(ano,tmp1,by="girl")
ano$y.1 <- ano$b0 + ano$b1*ano$time 
y <- subset(ano,ano$girl==1)
plot(y$time,y$y.1,
     type="n", 
     ylim=c(65,104),
	 xlab="Before/After Treatment",
	 ylab="Weight in Pounds",
	 main="weight ~ time + (1 | girl)",
	 xaxt="n"
	 	 )
axis(1, at = seq(1,2,by=1),labels=c("Before","After"))		 
for (i in 1:n){
   g <- subset(ano,ano$girl==i)
   lines(g$time, g$y.1, col=i)
}
```

Now make the model a bit more complex by adding in treatment

```{r  results=FALSE}

########################################################
#
# Model 2:  weight ~ time + treatment + (1 | girl)      #
########################################################

ano$treatment <- as.factor(ano$Rx)

model2reml <- lmer(weight ~ time + treatment + (1 | girl), data=ano, REML=TRUE)
summary(model2reml)
```
Now graphing conditional values
```{r}
b2 <- data.frame(coef(model2reml)[[1]])
tmp2 <- cbind(b2,girl)
names(tmp2) <- c("b0.2","b1.2","b2.2","b3.2","girl")
ano <- merge(ano,tmp2,by="girl")
ano$y.2 <- ano$b0.2 + ano$b1.2*ano$time + ano$b2.2*ano$Rx2 + ano$b3.2*ano$Rx3
y <- subset(ano,ano$girl==1)
plot(y$time,y$y.2,
     type="n", 
     ylim=c(65,104),
	 xlab="Before/After Treatment",
	 ylab="Weight in Pounds",
	 main="weight ~ time + treatment + (1 | girl)",
	 xaxt="n"
	 	 )
axis(1, at = seq(1,2,by=1),labels=c("Before","After"))		 
for (i in 1:n){
   g <- subset(ano,ano$girl==i)
   lines(g$time, g$y.2, col=i)
}
```

Doesn't look much different, so let's add an interaction between time and treatment.  We see this is in the data.

```{r  results=FALSE}
#######################################################
# Model 3: weight ~ time + treatment + time*treatment + ( 1 | girl) 
#
########################################################

model3reml <- lmer(weight ~ time + treatment + time*treatment + ( 1 | girl),data=ano,REML=TRUE)
summary(model3reml)
```
The fitted values are 
```{r}
b3 <- data.frame(coef(model3reml)[[1]])
tmp3 <- cbind(b3,girl)
names(tmp3) <- c("b0.3","b1.3","b2.3","b3.3","b4.3","b5.3","girl")
ano <- merge(ano,tmp3,by="girl")
ano$y.3 <- ano$b0.3 + ano$b1.3*ano$time + ano$b2.3*ano$Rx2 + ano$b3.3*ano$Rx3 + ano$b4.3*ano$time*ano$Rx2 + ano$b5.3*ano$time*ano$Rx3
y <- subset(ano,ano$girl==1)
plot(y$time,y$y.3,
     type="n", 
     ylim=c(65,104),
	 xlab="Before/After Treatment",
	 ylab="Weight in Pounds",
	 main="~time+treatment+time*treatment+( 1 | girl)",
	 xaxt="n"
	 	 )
axis(1, at = seq(1,2,by=1),labels=c("Before","After"))		 
for (i in 1:n){
   g <- subset(ano,ano$girl==i)
   lines(g$time, g$y.3, col=i)
}
```

Based on our initial look at the data, we might want a random slope but NO random intercept.

```{r  results=FALSE}
########################################################
# Model 4: weight ~ time + treatment + time*treatment + ( 0 +time|girl) 
########################################################

model4reml <- lmer(weight ~ time + treatment + time*treatment + ( 1 | girl) + ( 0 + time |girl),data=ano,REML=TRUE)
model4reml <- lmer(weight ~ time + treatment + time*treatment + ( 0 + time |girl),data=ano,REML=TRUE)
summary(model4reml)
```

What are fitted values now? These "look" like the data.

```{r}
b4 <- data.frame(coef(model4reml)[[1]])
tmp4 <- cbind(b4,girl)
names(tmp4) <- c("b0.4","b1.4","b2.4","b3.4","b4.4","b5.4","girl")
ano <- merge(ano,tmp4,by="girl")
ano$y.4 <- ano$b0.4 + ano$b1.4*ano$time + ano$b2.4*ano$Rx2 + ano$b3.4*ano$Rx3 + ano$b4.4*ano$time*ano$Rx2 + ano$b5.4*ano$time*ano$Rx3
plot(ano$time,ano$y.4,
     type="n", 
     ylim=c(65,104),
	 xlab="Before/After Treatment",
	 ylab="Weight in Pounds",
	 main="~ ... +time*treatment+(0+time|girl)",
	 xaxt="n"
	 	 )
axis(1, at = seq(1,2,by=1),labels=c("Before","After"))		 
for (i in 1:n){
   g <- subset(ano,ano$girl==i)
   lines(g$time, g$y.4, col=i)
}
```
## Summary of model fit via MLE

```{r}
screenreg(c(model1reml, model2reml,  model3reml, model4reml))
```

# Bayesian estimation

To illustrate how to fit models, we start simple and work to more complex.

## Model 0:  The random intercept model

This model is also called the Null HLM and is the same as random effects ANOVA.

```{r}
########################################################

# Bayesian Model 0     : random intercept no predictors (null hlm)     
########################################################

dataList <- list(
      y = ano$weight,          # sd of weight (response)
		 sdY = sd(ano$weight),     # sd of response
		 n = length(ano$weight),   # number obs
		 ng= length(ano$weight)/2, # number of girls
		 girl= ano$girl
		)
```
The model part
```{r}
model0 <- "model { 
  for (i in 1:n) {                   # likelihood
	     y[i] ~ dnorm(mu[i],precision)
       mu[i] <-  U0j[girl[i]]        # note index 
  }		
  for (j in 1:ng) {                  # random intercepts
	       U0j[j] ~ dnorm(g0,ptau)
	}
	g0 ~ dnorm(0,1/(100*sdY^2))    # other priors, mean of random
  tau ~ dunif(0.0001,200)        # sd of random
  ptau <- 1/tau^2                # precision random
  sigma ~ dunif(0.0001,2000)     # sd of response (weight)
  precision <- 1/sigma^2         # precision of weight/response
total.var <- sigma^2 + tau^2
icc <- tau^2/total.var         
}"

writeLines(model0, con="model0.txt")
```
I used different random number generators to show you what is available.  This shouldn't make a difference in results.
```{r}
start1 = list("g0"=mean(ano$weight), "sigma"=sd(ano$weight), "tau"=.5, .RNG.name="base::Wichmann-Hill", .RNG.seed=523) 
start2 = list("g0"=rnorm(1,0,3),     "sigma"=5,              "tau"=1, .RNG.name="base::Marsaglia-Multicarry", .RNG.seed=57)
start3 = list("g0"=rnorm(1,3,4),     "sigma"=10,             "tau"=5, .RNG.name="base::Super-Duper", .RNG.seed=24)
start4 = list("g0"=rnorm(1,-3,10),   "sigma"=50,             "tau"=20, .RNG.name="base::Mersenne-Twister", .RNG.seed=72100)
                   
start <- list(start1,start2,start3,start4)
```
And run jags
```{r}
model0.runjags <- run.jags(model=model0,
          method="parallel",  
				  monitor=c("g0", "sigma", "tau","total.var", "icc"),
          data=dataList, 
				  n.chains=4, 
				  sample=50000, 
				  burnin=5000, 
				  inits=start, 
				  thin=15)          
plot(model0.runjags)
print(model0.runjags)
summary(model1reml)
```

### Compute ICC (intra-class correlation) and credible intervals

There are multiple ways that we can do this: (1) within jags (see above), (2) monte carlo after obtaining posteriors, (3) multiple posteriors which is illustrated below.


```{r}
samps <-combine.mcmc(mcmc.objects = model0.runjags$mcmc)
# I need to know which columns I need
colnames(samps)
sigma2 <- samps[,2]**2
tau2   <- samps[,3]**2
icc <- tau2/(sigma2+tau2)
summary(icc)
```

## Model 1:  random intercept 
Now for our next model where we add in a predictor, "time":


```{r}
####################################################################
# Bayesian Model 1     : random intercept and fixed time           #
####################################################################
dataList <- list(
         y = ano$weight,
				 time = ano$time,        
				 sdY = sd(ano$weight),
				 n = length(ano$weight),
				 ng= length(ano$weight)/2,
				 girl= ano$girl
				 )
				 
model1 <- "model { 
       for (i in 1:n) {
	     y[i] ~ dnorm(mu[i],precision)
		   mu[i] <-  b0j[girl[i]] + g1*time[i] 
		}
		
	   for (j in 1:ng) {
	     b0j[j] ~ dnorm(g0,ptau)
	    }
		
	    g0 ~ dnorm(0,1/(100*sdY^2))
	  	g1 ~ dnorm(0,1/(100*sdY^2))

        tau ~ dunif(0.0001,200)
        ptau <- 1/tau^2
		sigma ~ dunif(0.0001,2000)
        precision <- 1/sigma^2	
		}"

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

start1 = list("g0"=mean(ano$weight),"g1"=rnorm(1,0,3),      "sigma"=sd(ano$weight), "tau"=.5, 
             .RNG.name="base::Wichmann-Hill", .RNG.seed=523) 
start2 = list("g0"=rnorm(1,0,3),    "g1"=rnorm(1,0,3),      "sigma"=5,              "tau"=1, 
             .RNG.name="base::Marsaglia-Multicarry", .RNG.seed=57)
start3 = list("g0"=rnorm(1,3,4),    "g1"=rnorm(1,3,4),      "sigma"=50,             "tau"=3, 
             .RNG.name="base::Super-Duper", .RNG.seed=24)
start4 = list("g0"=rnorm(1,-3,10),  "g1"=rnorm(1,-3,10),    "sigma"=10,             "tau"=10, 
             .RNG.name="base::Mersenne-Twister", .RNG.seed=72100)

start <- list(start1,start2,start3,start4)

start <- list(start1,start2,start3,start4)

model1.runjags <- run.jags(model=model1,      
                  method="parallel",  
                  monitor=c("g0", "g1", "sigma", "tau"),
                  data=dataList,
				          sample=20000,				  
				          n.chains=4,
                  thin=10,				  
				         inits=start)

# And checking a bit more whether model1.runjags converged
plot(model1.runjags)

print(model1.runjags)
# For comparison
summary(model1reml)	
```

After fitting model, I decided I want an ICC, so...
```{r}
# Compute residual ICC (intra-class correlation) and credible interval
samps <-combine.mcmc(mcmc.objects = model1.runjags$mcmc)
colnames(samps)
sigma2 <- samps[,3]**2
tau2   <- samps[,4]**2
icc1 <- tau2/(sigma2+tau2)
summary(icc1)
```

This model is too simple for the data, so we'll go on to our next model, model2 which adds treatment as a fixed effect.

### Model 2: Random intercept + fixed time and Rx
```{r}
##################################################################
# Bayesian Model 2 : random intercept + fixed time   + Rx        #
##################################################################
dataList <- list(
         y = ano$weight,
				 time = ano$time,
				 rx1 = ano$Rx1,
				 rx3 = ano$Rx3,
				 sdY = sd(ano$weight),
				 n = length(ano$weight),
				 ng= length(ano$weight)/2,
				 girl= ano$girl
				 )
				 
model2 <- "model { 
       for (i in 1:n) {
	     y[i] ~ dnorm(mu[i],precision)
		   mu[i] <- b0j[girl[i]] + g1*time[i] + g2*rx1[i] + g3*rx3[i]
		}
		
	   for (j in 1:ng) {
	     b0j[j] ~ dnorm(g0,ptau)
	    }
		
	    g0 ~ dnorm(0,1/(100*sdY^2))
		  g1 ~ dnorm(0,1/(100*sdY^2))
	    g2 ~ dnorm(0,1/(100*sdY^2))
		  g3 ~ dnorm(0,1/(100*sdY^2))

      tau ~ dunif(0.0001,200)
      ptau <- 1/tau^2
		  sigma ~ dunif(0.0001,2000)
      precision <- 1/sigma^2	
		}"

writeLines(model2, con="model2.txt")

start1 = list("g0"=mean(ano$weight),"g1"=rnorm(1,0,3),  "g2"=rnorm(1,0,3),"g3"=rnorm(1,0,3),     
              "sigma"=sd(ano$weight), "tau"=.5,   .RNG.name="base::Wichmann-Hill", .RNG.seed=523) 
			  
start2 = list("g0"=rnorm(1,0,3),    "g1"=rnorm(1,0,3), "g2"=rnorm(1,-1,3), "g3"=rnorm(1,0,3),        
              "sigma"=5,            "tau"=1,      .RNG.name="base::Marsaglia-Multicarry", .RNG.seed=57)
			  
start3 = list("g0"=rnorm(1,3,4),    "g1"=rnorm(1,3,4),   "g2"=rnorm(1,3,4), "g3"=rnorm(1,0,3),   
              "sigma"=50,           "tau"=3,      .RNG.name="base::Super-Duper", .RNG.seed=24)
			  
start4 = list("g0"=rnorm(1,-3,10),  "g1"=rnorm(1,-3,10),   "g2"=rnorm(1,10,5), "g3"=rnorm(1,0,3),
               "sigma"=10,          "tau"=10,     .RNG.name="base::Mersenne-Twister", .RNG.seed=72100)

start <- list(start1,start2,start3,start4)

model2.runjags <- run.jags(model=model2,      
                  method="parallel",  
                  monitor=c("g0", "g1", "g2","g3","sigma", "tau"),
                  data=dataList,
				          sample=20000,				  
				          n.chains=4,
				          thin=10,
				          inits=start)

# And we should do some more model checking...
plot(model2.runjags)


print(model2.runjags)
# For comparison
summary(model2reml)
```
### Compute residual ICC (intra-class correlation) and credible interval
```{r}
samps <-combine.mcmc(mcmc.objects = model2.runjags$mcmc)
colnames(samps)        # to see which columns I need
sigma2 <- samps[,5]**2
tau2   <- samps[,6]**2
icc2 <- tau2/(sigma2+tau2)
summary(icc2)
```


### Is there a difference between Rx1 and Rx3 (cognitive and family)?
```{r}
g2 <- samps[,3]
g3 <- samps[,4]
diff <- g3-g2
summary(diff)
hist(diff, main="Posterior of Differences between  family & cognitive", col="cyan")

### What is the probability that family is better than cognitive?
prob <- mean(diff>0)
```
These results should not be treated as final, because we have more models to fit.  However, this does illustrate how to be get postiors samples, compute variaous posteriors, get probabilities for hypothses, etc.


## Model 3: Random intercept + fixed time, Rx & time*Rx

Next model (this will run fine but it is not run here)

```{r eval=FALSE}
##################################################################################
# Bayesian Model 3     : random intercept + fixed time   + Rx    + time*Rx       #
##################################################################################
dataList <- list(
       y = ano$weight,
			time = ano$time,
			rx1 = ano$Rx1,
			rx3 = ano$Rx3,
			sdY = sd(ano$weight),
			n = length(ano$weight),
			ng= length(ano$weight)/2,
			girl= ano$girl
			)
				 
model3 <- "model { 
       for (i in 1:n) {
	     y[i] ~ dnorm(mu[i],precision)
		   mu[i] <- b0j[girl[i]]  + g1*time[i] + g2*rx1[i] + g3*rx3[i] + g4*time[i]*rx1[i] +  g5*time[i]*rx3[i]
		}
		
	   for (j in 1:ng) {
	     b0j[j] ~ dnorm(g0,ptau)
	    }
		
	    g0 ~ dnorm(0,1/(100*sdY^2))
		  g1 ~ dnorm(0,1/(100*sdY^2))
	    g2 ~ dnorm(0,1/(100*sdY^2))
		  g3 ~ dnorm(0,1/(100*sdY^2))
	    g4 ~ dnorm(0,1/(100*sdY^2))
		  g5 ~ dnorm(0,1/(100*sdY^2))

      tau ~ dunif(0.0001,200)
      ptau <- 1/tau^2
	  	sigma ~ dunif(0.0001,2000)
      precision <- 1/sigma^2	
		}"

writeLines(model3, con="model3.txt")

start1 = list("g0"=mean(ano$weight),"g1"=rnorm(1,0,3),  "g2"=rnorm(1,0,3),"g3"=rnorm(1,0,3),  
              "g4"=rnorm(1,0,3),"g5"=rnorm(1,0,3),      "sigma"=sd(ano$weight), "tau"=.5,   
			  .RNG.name="base::Wichmann-Hill", .RNG.seed=523) 
			  
start2 = list("g0"=rnorm(1,0,3),  "g1"=rnorm(1,0,3), "g2"=rnorm(1,-1,3), "g3"=rnorm(1,0,3),
              "g4"=rnorm(1,-1,3), "g5"=rnorm(1,0,3),   "sigma"=5,        "tau"=1,     
			  .RNG.name="base::Marsaglia-Multicarry", .RNG.seed=57)
			  
start3 = list("g0"=rnorm(1,3,4),  "g1"=rnorm(1,3,4), "g2"=rnorm(1,3,4), "g3"=rnorm(1,0,3),   
              "g4"=rnorm(1,3,4), "g5"=rnorm(1,0,3),  "sigma"=50,        "tau"=5,     
			  .RNG.name="base::Super-Duper", .RNG.seed=24)
			  
start4 = list("g0"=rnorm(1,-3,10),"g1"=rnorm(1,-3,10), "g2"=rnorm(1,10,5), "g3"=rnorm(1,0,3),
              "g4"=rnorm(1,2,5), 
			  "g5"=rnorm(1,3,3),  "sigma"=10,          "tau"=10,     
			  .RNG.name="base::Mersenne-Twister", .RNG.seed=72100)

start <- list(start1,start2,start3,start4)

model3.runjags <- run.jags(model=model3,      
          method="parallel",  
          monitor=c("g0","g1","g2","g3","g4","g5","sigma","tau"),
          data=dataList,
				  sample=20000,				  
				  n.chains=4,
				  thin=20,
				  inits=start)

print(model3.runjags)
summary(model3.reml)		
plot(model3.runjags)
```

### Compute Residuals ICC (intra-class correlation) when using runjags
```{r eval=FALSE}
samps <-combine.mcmc(mcmc.objects = model3.runjags$mcmc)
colnames(samps)        # to see which columns I need
sigma2 <- samps[,7]**2
tau2   <- samps[,8]**2
icc3 <- tau2/(sigma2+tau2)
summary(icc3)
```


## Model 4:  Random slope, fixed time, Rx & time*Rx

Now we'll switch to random slope and drop the random intercept.  This might be OK.  This code will run but not run here.
```{r eval=FALSE}
##################################################################################
# Bayesian Model  4   : random slope + fixed time   + Rx    + time*Rx            #
##################################################################################
dataList <- list(
      y = ano$weight,
			time = ano$time,
			rx1 = ano$Rx1,
			rx3 = ano$Rx3,
			sdY = sd(ano$weight),
			n = length(ano$weight),
			ng= length(ano$weight)/2,
			girl= ano$girl
			)
				 
model4 <- "model { 
       for (i in 1:n) {
	      y[i] ~ dnorm(mu[i],precision)
		    mu[i] <- g0 + g1*time[i] + b1j[girl[i]]*time[i] + g2*rx1[i] + g3*rx3[i] 
		              + g4*time[i]*rx1[i] +  g5*time[i]*rx3[i]
		}
		
	   for (j in 1:ng) {
	     b1j[j] ~ dnorm(0,ptau)
	    }
		
	    g0 ~ dnorm(0,1/(100*sdY^2))
	    g1 ~ dnorm(0,1/(100*sdY^2))
	    g2 ~ dnorm(0,1/(100*sdY^2))
		  g3 ~ dnorm(0,1/(100*sdY^2))
	    g4 ~ dnorm(0,1/(100*sdY^2))
		  g5 ~ dnorm(0,1/(100*sdY^2))

      tau ~ dunif(0.0001,200)
      ptau <- 1/tau^2
		  sigma ~ dunif(0.0001,2000)
      precision <- 1/sigma^2	
		}"

writeLines(model4, con="model4.txt")

start1 = list("g0"=mean(ano$weight), "g1"=rnorm(1,1,3),  "g2"=rnorm(1,0,3),"g3"=rnorm(1,0,3),  
              "g4"=rnorm(1,0,3),     "g5"=rnorm(1,0,3),      "sigma"=sd(ano$weight), "tau"=.5,   
			  .RNG.name="base::Wichmann-Hill", .RNG.seed=523) 
			  
start2 = list("g0"=rnorm(1,0,3), "g1"=rnorm(1,-2,3), "g2"=rnorm(1,-1,3), "g3"=rnorm(1,0,3),
              "g4"=rnorm(1,-1,3),"g5"=rnorm(1,0,3),   "sigma"=5,        "tau"=1,     
			  .RNG.name="base::Marsaglia-Multicarry", .RNG.seed=57)
			  
start3 = list("g0"=rnorm(1,3,4), "g1"=rnorm(1,0,3), "g2"=rnorm(1,3,4), "g3"=rnorm(1,0,3),   
              "g4"=rnorm(1,3,4), "g5"=rnorm(1,0,3),  "sigma"=50,        "tau"=5,     
			  .RNG.name="base::Super-Duper", .RNG.seed=24)
			  
start4 = list("g0"=rnorm(1,-3,10), "g1"=rnorm(1,5,3),"g2"=rnorm(1,10,5), "g3"=rnorm(1,0,3),
              "g4"=rnorm(1,2,5),   "g5"=rnorm(1,3,3),  "sigma"=10,          "tau"=10,     
			  .RNG.name="base::Mersenne-Twister", .RNG.seed=72100)

start <- list(start1,start2,start3,start4)

model4.runjags <- run.jags(model=model4,      
          method="parallel",  
          monitor=c("g0","g1","g2","g3","g4","g5", "sigma", "tau"),
          data=dataList,
				  sample=20000,				  
				  n.chains=4,
				  thin=5,
				  inits=start)

plot(model4.runjags)
print(model4.runjags)
summary(model4reml)		

```


## Model 5:  Random intecept & Slope with 2 time points!

We'll try a random intercept and slope.  With only 2 time points, you can't fit both random intercept and slope as well as correlation between the random intercept and slope. 
In this example, $$\tau_01=0$$. This seems to be just enough to get the model to fit.  Let's see what happens when you try to get estimate both. 
```{r eval=FALSE}
##################################################################################
# Bayesian Model  5: random intercept & slope + fixed time  + Rx1 + Rx3          #
#   density of tau_0 doesn't look good                                           #
##################################################################################
dataList <- list(
      y = ano$weight,
			time = ano$time,
			rx1 = ano$Rx1,
			rx3 = ano$Rx3,
			sdY = sd(ano$weight),
			n = length(ano$weight),
			ng= length(ano$weight)/2,
			girl= ano$girl
				 )
				 
model5 <- "model { 
       for (i in 1:n) {
	     y[i] ~ dnorm(mu[i],precision)
		   mu[i] <- g0 + g1*time[i] + g2*rx1[i] + g3*rx3[i] + U0j[girl[i]] + U1j[girl[i]]*time[i]
		}
		
	   for (j in 1:ng) {
	     U0j[j] ~ dnorm(0,ptau0)
	     U1j[j] ~ dnorm(0,ptau1)
	    }
		
	    g0 ~ dnorm(0,1/(100*sdY^2))
	    g1 ~ dnorm(0,1/(100*sdY^2))
	    g2 ~ dnorm(0,1/(100*sdY^2))
	    g3 ~ dnorm(0,1/(100*sdY^2))
		
      ptau0 ~ dgamma(0.001,0.001)
      tau0 <- 1/sqrt(ptau0)
		
      tau1 ~ dunif(0.0001,200)
      ptau1 <- 1/tau1^2
		
		sigma ~ dunif(0.0001,2000)
        precision <- 1/sigma^2	
		}"

writeLines(model5, con="model5.txt")

start1 = list("g0"=mean(ano$weight), "g1"=rnorm(1,1,3), "g2"=rnorm(1,1,3), "g3"=rnorm(1,4,3), "sigma"=sd(ano$weight), "ptau0"=.005,    "tau1"=2, 
			  .RNG.name="base::Wichmann-Hill", .RNG.seed=523) 
			  
start2 = list("g0"=dnorm(1,0,3), "g1"=rnorm(1,-2,3), "g2"=rnorm(1,0,3), "g3"=rnorm(1,-1,3), "sigma"=5,   "ptau0"=.1, "tau1"=3,    
			  .RNG.name="base::Marsaglia-Multicarry", .RNG.seed=57)
			  
start3 = list("g0"=dnorm(1,3,4), "g1"=rnorm(1,0,3), "g2"=rnorm(1,-1,3), "g3"=rnorm(1,0,4),  "sigma"=50,  "ptau0"=.045,  "tau1"=10,     
			  .RNG.name="base::Super-Duper", .RNG.seed=24)
			  
start4 = list("g0"=dnorm(1,-3,10), "g1"=rnorm(1,5,3), "g2"=rnorm(1,5,3), "g3"=rnorm(1,3,4), "sigma"=10, "ptau0"=.10, "tau1"=1,      
			  .RNG.name="base::Mersenne-Twister", .RNG.seed=72100)

start <- list(start1,start2,start3,start4)

model5.runjags <- run.jags(model=model5,      
                  method="parallel",  
                  monitor=c("g0","g1","g2","g3","sigma","ptau0","ptau1","tau0",  "tau1"),
                  data=dataList,
				  sample=20000,				  
				  n.chains=4,
				  thin=20,
				  inits=start)

print(model5.runjags)		
plot(model5.runjags)
```

## Refinement of Model 4
So, we'll take Model 4 as our final model and get some additional statistics from the code (i.e., draws from posterior distribution).
```{r}
###############################################################################
# Bayesian Model  4 refinement : random slope + fixed time   + rx3+ time*Rx3  #
# 1st: drop time*Rx1                                                          #
# 2nd: drop Rx1                                                               #
###############################################################################
dataList <- list(
         y = ano$weight,
				 time = ano$time,
				 rx3 = ano$Rx3,
				 sdY = sd(ano$weight),
				 n = length(ano$weight),
				 ng= length(ano$weight)/2,
				 girl= ano$girl
				 )
				 
model4 <- "model { 
       for (i in 1:n) {
	     y[i] ~ dnorm(mu[i],precision)
		   mu[i] <- g0 + g1*time[i] + b1j[girl[i]]*time[i]  + g3*rx3[i]  +  g5*time[i]*rx3[i]
		   yhat[i] <- g0 + g1*time[i] + b1j[girl[i]]*time[i]  + g3*rx3[i]  +  g5*time[i]*rx3[i]
		}
		
	   for (j in 1:ng) {
	     b1j[j] ~ dnorm(0,ptau)
	    }
		
	   g0 ~ dnorm(0,1/(100*sdY^2))
	   g1 ~ dnorm(0,1/(100*sdY^2))
		 g3 ~ dnorm(0,1/(100*sdY^2))
		 g5 ~ dnorm(0,1/(100*sdY^2))

     tau ~ dunif(0.0001,200)
     ptau <- 1/tau^2
		 sigma ~ dunif(0.0001,2000)
     precision <- 1/sigma^2	
		}"

writeLines(model4, con="model4.txt")

start1 = list("g0"=mean(ano$weight), "g1"=dnorm(1,1,3), "g3"=dnorm(1,0,3),  
                "g5"=dnorm(1,0,3),      "sigma"=sd(ano$weight), "tau"=.5,   
			  .RNG.name="base::Wichmann-Hill", .RNG.seed=523) 
			  
start2 = list("g0"=dnorm(1,0,3), "g1"=dnorm(1,-2,3), "g3"=dnorm(1,0,3),
              "g5"=dnorm(1,0,3),   "sigma"=5,        "tau"=1,     
			  .RNG.name="base::Marsaglia-Multicarry", .RNG.seed=57)
			  
start3 = list("g0"=dnorm(1,3,4), "g1"=dnorm(1,0,3), "g3"=dnorm(1,0,3),   
               "g5"=dnorm(1,0,3),  "sigma"=50,        "tau"=5,     
			  .RNG.name="base::Super-Duper", .RNG.seed=24)
			  
start4 = list("g0"=dnorm(1,-3,10), "g1"=dnorm(1,5,3),"g3"=dnorm(1,0,3),
              "g5"=dnorm(1,3,3),  "sigma"=10,          "tau"=10,     
			  .RNG.name="base::Mersenne-Twister", .RNG.seed=72100)

start <- list(start1,start2,start3,start4)

model4.runjags <- run.jags(model=model4,      
          method="parallel",  
          monitor=c("g0","g1","g3","g5", "sigma", "tau","yhat"),
          data=dataList,
				  sample=20000,				  
				  n.chains=4,
				  thin=5,
				  inits=start)

print(model4.runjags)
summary(model4reml)	
```

Why didn't print(model4.runjags) give you the table of statistics?   There are many of them; in particular, there are an extra 72 X 2 = 143 yhat values, 2 values per girl (weight before and weight after). 

```{r eval=FALSE}
plot(model4.runjags)
```

Let's do some other things with posterior samples.  This is a better (more straightforward) than previous version that I had here.
```{r}
# Takes output and put in format so can use coda plots and diagnostics on it.
samps <- combine.mcmc(model4.runjags)
colnames(samps)

g0 <- samps[,1]
g1 <- samps[,2]
g3 <- samps[,3]
g5 <- samps[,4]
sigma <- samps[,5]
tau <- samps[,6]
fitted <- samps[,7:150]

ano$yhat <- apply(fitted,2,"mean")
plot(ano$time,ano$yhat,
     type="n", 
     ylim=c(65,104),
	 xlab="Before/After Treatment",
	 ylab="Weight in Pounds",
	 main="Bayesian refined Model 4",
	 xaxt="n"
	 	 )
axis(1, at = seq(1,2,by=1),labels=c("Before","After"))		 
for (i in 1:n){
   g <- subset(ano,ano$girl==i)
   lines(g$time, g$yhat, col=i)
}
```

This graph of fitted values resembles pattern in the data.

More methods for checked convergence should be done.