%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Purpose: Produce Reproducible output for Bayesian model fit in JAGS
%%%%%% Date: Sept 2018
%%%%%% Author: Ben Saville, Berry Consultants
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\documentclass{article}

\usepackage{setspace,booktabs,float}                % for latex(describe()) (Hmisc package) and latex.table
\usepackage{longtable}                                % for tables that break over multiple pages
\usepackage[pdftex]{lscape}
\usepackage[margin=0.5in,bottom=.5in,includehead,includefoot]{geometry}       % sets 1-inch margins
\setlength{\oddsidemargin}{.25in}
\setlength{\evensidemargin}{.25in}
\setlength{\textwidth}{6in}



\title{Bayesian Analysis via R/JAGS \\
       Reproducible Report (knitr)}
\author{Prepared by Ben Saville, Ph.D.}
\date{\today}



% 
 \usepackage{fancyhdr}
  \pagestyle{fancy}
 % \lhead{\begin{picture}(0,0) \put(-50,-25){\includegraphics[scale=0.4]{BerryPic2.png}} \end{picture}}
\fancyhead[R]{}

\renewcommand{\headrulewidth}{0pt}
\newcommand{\bi}{\begin{itemize}}
\newcommand{\ei}{\end{itemize}}
\newcommand{\be}{\begin{enumerate}}
\newcommand{\ee}{\end{enumerate}}
\setlength\parindent{0pt}


\begin{document}


\maketitle
\thispagestyle{fancy}


\tableofcontents

  <<opts, include=FALSE>>=
  library(knitr)
  opts_chunk$set(
    echo = FALSE,
    tidy = FALSE,
    dev="pdf",
    fig.path="graphics/plot",
    # out.width=".47\\textwidth",
    fig.keep="high",
    fig.show="hold",
    fig.align="center",
    # fig.lp="fig:",
    autodep = FALSE,
    comment=NA)
#knit_hooks$set(inline = function(x) {
#  if (is.numeric(x)) return(round(x, 3))})
par(las = 1)
options(width = 90, scipen = 6, digits = 3)
options(contrasts=c("contr.treatment","contr.treatment"))
@




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\section{Conjugate model}


<<conjugate, echo=FALSE,results='hide',message=FALSE, warning=FALSE>>=

## Simple linear regression model using JAGS
## Author: Ben Saville

library('rjags')


################################################################
# Data: allergy medication
# y contains the N hours of relief experienced
# x contains the N dosages.

    ## Input data as a list
df = list(
  y = c(9.1, 5.5, 12.3, 9.2, 14.2, 16.8, 22, 18.3, 
        24.5, 22.7), 
  x = c(3, 3, 4, 5, 6, 6, 7, 8, 8, 9), 
  N = 10
)
    
################## Classical Reg results ###################

freq.model <- lm(y~x,df)
summary(freq.model)
confint(freq.model)
  ## Variance estimate
(summary(freq.model)$sigma)^2
  
with(df,plot(x,y,xlab="dose",ylab="hours",main="Hours of relief by medication dose"))
abline(freq.model$coefficients)

################# Or input data from data set from a data frame (convert to list)
mydata =data.frame(cbind(df$y,df$x))
names(mydata) = c("y","x")    
    ## i.e., starting here
df2 = list(  'y'  = with(mydata, y),
            'x' = with(mydata, x),
            'N'  = nrow(mydata))

as.list(as.data.frame(mydata))

## You can fit model using either df or df2

############################# Fitting model with JAGS ############################################

################# Inital values
inits1 = list( beta0 = 0, beta1 = 0, tau = 1)                                                     # No seed
#inits1 = list( beta0 = 0, beta1 = 0, tau = 1, .RNG.seed=1, .RNG.name='base::Wichmann-Hill')      # Sets seed

################# Conjugate model ############################
    ## Create a txt file for use by JAGS
fileForJAGS <- "allergy-bugConjugate.txt"
cat("
# This file is created dynamically by the main knitr file for this project.
# Make any changes there, not here!
model {
for(i in 1:N)  {
    y[i] ~ dnorm(mu[i], tau);
    mu[i] <- beta0 + beta1*x[i];
    }

  beta0 ~ dnorm(0, 1/tau * 0.001);
  beta1 ~ dnorm(0, 1/tau * 0.001);
  tau ~ dgamma(0.01, 0.01);
  sig <- 1/tau;
}
", file= fileForJAGS)

# n.chains = number of chains to monitor (with different starting values)
# n.adapt = initial sample phase to optimize the sampling mechanism (e.g. Metropolis Hastings step size)
# inits = initial values from above
# "allergy-bugConjugate.txt" is the filename which contains the model

jags.Conj <- jags.model("allergy-bugConjugate.txt",data=df,n.chains=1, n.adapt=1000, inits=inits1)
  ## burn-in
update( jags.Conj , n.iter=2000 )
  ## Store MCMC samples after burn-in
  ## Arugments are the parameters to be monitored and the number of samples to take
mcmc.samples.Conj <- coda.samples(jags.Conj,c('beta0','beta1','sig'),10000)
summary(mcmc.samples.Conj)
plot(mcmc.samples.Conj)

    ### Calculate probability that the slope is greater than 1
post.sam = as.matrix(mcmc.samples.Conj[1])
mean(post.sam[,2] > 0 )

@

\clearpage
\section{Non-conjugate model}

<<nonconjugate, echo=FALSE,results='hide'>>=

################### Non-conjugate model ###################################
    ## Create a txt file for use by JAGS
fileForJAGS <- "allergy-bug.txt"
cat("
# This file is created dynamically by the main knitr file for this project.
# Make any changes there, not here!
model {
for(i in 1:N)  {
    y[i] ~ dnorm(mu[i], tau);
    mu[i] <- beta0 + beta1*x[i];
    }

  beta0 ~ dnorm(0, 0.001);
  beta1 ~ dnorm(0, 0.001);
  tau ~ dgamma(0.01, 0.01);
  sig <- 1/tau;
}
", file= fileForJAGS)


jags <- jags.model("allergy-bug.txt",data=df,n.chains=1,n.adapt=1000, inits=inits1)
  ## burn-in
update(jags, n.iter=2000 )
mcmc.samples <- coda.samples(jags,c('beta0','beta1','sig'),10000)
summary(mcmc.samples)
#png('PostGraphs.png')
plot(mcmc.samples)
#dev.off()

@

\clearpage
\section{Three chains, non-conjugate model}

<<Threechains, echo=FALSE,results='hide'>>=


######################### 3 chains & convergence diagnostics, Non-conjugate model ############

inits2 = list( beta0 = 0.087, beta1 = -3.914, tau = 1.976 )
inits3 = list( beta0 = 0.568, beta1 = 2.440, tau = 0.841 )

jags3 <- jags.model("allergy-bug.txt",data=df,n.chains=3,n.adapt=1000, inits=list(inits1,inits2,inits3))
  ## burn-in
update(jags3, n.iter=2000 )
  ## Samples
mcmc.samples3 <- coda.samples(jags3,c('beta0','beta1','sig'),10000)
summary(mcmc.samples3)
plot(mcmc.samples3)
gelman.plot(mcmc.samples3)
gelman.diag(mcmc.samples3)
HPDinterval(mcmc.samples3)


@


\clearpage
\section{Two explanatory variables, Non-conjugate model}

<<nonconjugate2, echo=FALSE,results='hide'>>=

############ Two Parameters, Non-conjugate model ###################################
    ## Create a txt file for use by JAGS
fileForJAGS <- "allergy-bug2.txt"
cat("
# This file is created dynamically by the main knitr file for this project.
# Make any changes there, not here!
model {
for(i in 1:N)  {
    y[i] ~ dnorm(mu[i], tau);
    mu[i] <- beta0 + beta1*x[i] + beta2*tmt[i];
    }

  beta0 ~ dnorm(0, 0.001);
  beta1 ~ dnorm(0, 0.001);
  beta2 ~ dnorm(0, 0.001);
  tau ~ dgamma(0.01, 0.01);
  sig <- 1/tau;
}
", file= fileForJAGS)

    ## Additional variable: Treatment status (binary)
df$tmt <- rep(c(1,0),5)

inits2param = list( beta0 = 0, beta1 = 0, tau = 1)   
jags <- jags.model("allergy-bug2.txt",data=df,n.chains=1,n.adapt=1000, inits=inits2param)
  ## burn-in
update(jags, n.iter=2000 )
mcmc.samples <- coda.samples(jags,c('beta0','beta1','beta2','sig'),10000)
summary(mcmc.samples)
plot(mcmc.samples)

with(df,plot(x,y,xlab="dose",ylab="hours", ylim=c(0,25),
      main="Hours of relief by medication dose & tmt", col=c("red","blue")[df$tmt+1]))
legend("topleft",c("Tmt","PBO"), pch=c(1,1), col=c("blue","red"))


@



\end{document}