# Edps 590BAY
# Spring 2019
# C.J.Anderson
#
# HLM of anorexia data
#
# Graphics of data 
#
# Multilevel models with lmer & graphing
#
# Bayes models:
#    Model 0:  Random intercept and no predictors
#    Model 1:  Random intercept:  time
#    Model 2:  Random intercept:  time + treatment
#    Model 3:  Random intercept:  time + treatment + time*treatment 
#    Model 4:  Random  slope
#
library(lme4)
library(lmerTest)
library(rjags)
library(runjags)
library(coda)

#setwd("C:\\Users\\cja\\Documents\\Dropbox\\edps 590BAY\\Lectures\\10 Multilevel models")

#setwd("C:\\Users\\cja\\Dropbox\\edps 590BAY\\Lectures\\10 Multilevel models")

# How we looked at the data ordinal and now

ano.wide <- read.table("anorexia_wide.txt",header=TRUE)
ano.wide$change <- ano.wide$weight2 - ano.wide$weight1

n <- length(ano.wide$change)

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")

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)

par(mfrow=c(1,1))
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)
}
   
# 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

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)

############################################################
# Simple Model using lme4                                  #
############################################################

# Remember: treatment = 1 Cognitive
#                     = 2 Control
#                     = 3 Family

model1.mle <- lmer(weight ~ time + (1 | girl),data=ano,REML=FALSE)
summary(model1.mle)

girl <- data.frame(seq(1:n))
par(mfrow=c(2,2))

b1 <- data.frame(coef(model1.mle)[[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)
}

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

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

model2.mle <- lmer(weight ~ time + treatment + (1 | girl),data=ano,REML=FALSE)
summary(model2.mle)

b2 <- data.frame(coef(model2.mle)[[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)
}

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

model3.mle <- lmer(weight ~ time + treatment + time*treatment + ( 1 | girl),data=ano,REML=FALSE)
summary(model3.mle)

b3 <- data.frame(coef(model3.mle)[[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)
}


#####################################################################
# Model 4: weight ~ time + treatment + time*treatment + ( 0 +time|girl) #
#####################################################################
#model4.mle <- lmer(weight ~ time + treatment + time*treatment + ( 1 | girl) + ( 0 + time |girl),data=ano,REML=FALSE)

model4.mle <- lmer(weight ~ time + treatment + time*treatment + ( 0 + time |girl),data=ano,REML=FALSE)
summary(model4.mle)

b4 <- data.frame(coef(model4.mle)[[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)
}


##################################################################################
# Bayesian Model 0     : random intercept no predictors                          #
##################################################################################

model0.reml <- lmer(weight ~ (1 | girl),data=ano,REML=TRUE)
summary(model0.reml)

dataList <- list(
                 y = ano$weight,
				 sdY = sd(ano$weight),
				 n = length(ano$weight),
				 ng= length(ano$weight)/2,
				 girl= ano$girl
				 )
				 
model0 <- "model { 
       for (i in 1:n) {
	     y[i] ~ dnorm(mu[i],precision)
		 mu[i] <-  U0j[girl[i]]  
		}
		
	   for (j in 1:ng) {
	    U0j[j] ~ dnorm(g0,ptau)
	    }
			 
	    g0 ~ 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(model0, con="model0.txt")

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)

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

print(model0.runjags)
summary(model0.reml)		
plot(model0.runjags)

##################################################################################
# Bayesian Model 1     : random intercept and fixed time                         #
##################################################################################

model1.reml<- lmer(weight ~ time + (1 | girl),data=ano,REML=TRUE)
summary(model1.reml)

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=19,				  
				  inits=start)

print(model1.runjags)
summary(model1.reml)		
plot(model1.runjags)


##################################################################################
# Bayesian Model 2     : random intercept + fixed time   + Rx                    #
##################################################################################

model2.reml<- lmer(weight ~ time + Rx1 + Rx3 +(1 | girl),data=ano,REML=TRUE)
summary(model2.reml)

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)

print(model2.runjags)
summary(model2.reml)		
plot(model2.runjags)




##################################################################################
# Bayesian Model 3     : random intercept + fixed time   + Rx    + time*Rx       #
##################################################################################

model3.reml<- lmer(weight ~ time + Rx1 + Rx3 + time*Rx1 + time*Rx3 + (1 | girl),data=ano,REML=TRUE)
summary(model3.reml)

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)



##################################################################################
# Bayesian Model  4   : random slope + fixed time   + Rx    + time*Rx            #
##################################################################################

model4.reml<- lmer(weight ~ time + Rx1 + Rx3 + time*Rx1 + time*Rx3 + (0 + time | girl),data=ano,REML=TRUE)
summary(model4.reml)

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)

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



##################################################################################
# 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",  "tau1"),
                  data=dataList,
				  sample=20000,				  
				  n.chains=4,
				  thin=20,
				  inits=start)

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

#############  Probably don't need both random slope and intercept #######################
#JAGS model summary statistics from 80000 samples (thin = 20; chains = 4; adapt+burnin = 5000):
#                                                                                           
#       Lower95 Median Upper95   Mean      SD Mode     MCerr MC%ofSD SSeff     AC.400   psrf
#g0      74.503 77.678  80.821 77.675  1.6103   -- 0.0083055     0.5 37589 -0.0047055 1.0002
#g1     0.98671 2.7887  4.5941 2.7846 0.92155   -- 0.0046962     0.5 38508  0.0012152 1.0001
#g2    -0.44021 2.3134  5.1091    2.3  1.4138   -- 0.0070581     0.5 40126 -0.0026701 1.0001
#g3     0.91706 4.2921  7.6316 4.2688  1.7173   -- 0.0087671     0.5 38369  0.0033132 1.0001
#sigma   4.2966 5.1457  6.1231 5.1802 0.47197   -- 0.0025477     0.5 34320  0.0087721 1.0001
#ptau0 0.032431 11.653  659.89 126.71  330.29   --    3.7299     1.1  7841    0.18711 1.0018
          ^^      ^^      ^^                            
#tau1    1.4083 2.5115  3.5236 2.4871 0.53629   -- 0.0036089     0.7 22082   0.016333 1.0009
#
#Total time taken: 49.7 seconds
#
# tau0 = 1/sqrt(ptau0) = 1/sqrt(126.71) = 0.0888





##################################################################################
# Bayesian Model  4 refinement : random slope + fixed time   + rx3+ time*Rx3     #
# 1st: drop time*Rx1
# 2nd: drop Rx1
##################################################################################

model4.reml<- lmer(weight ~ time  + Rx3  + time*Rx3 + (0 + time | girl),data=ano,REML=TRUE)
summary(model4.reml)

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(model4.reml)		
plot(model4.runjags)


# Takes output and put in format so can use coda plots and diagnostics on it.
samp4 <- as.mcmc.list(model4.runjags)                               # change to mcmc object so that can use coda.
samples <- as.array(samp4)                                           # change from mcmc to array object, one from each chain.
sampled <- rbind(samples[,,1],samples[,,2],samples[, , 3], samples[, , 4]) 
write(sampled,"sampled_model4.txt")

g0 <- sampled[,1]
g1 <- sampled[,2]
g3 <- sampled[,3]
g5 <- sampled[,4]
sigma <- sampled[,5]
tau <- sampled[,6]
fitted <- sampled[,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)
}