\documentclass{article}
\usepackage{graphics}
\evensidemargin 0cm
\oddsidemargin 0cm
\topmargin -2cm
\textwidth 15cm
\textheight 24cm
\columnsep 2cm
\newcommand{\E}{{\rm E}}
\newcommand{\var}{{\rm var}}
\newcommand{\covar}{{\rm covar}}
\newcommand{\bias}{{\rm b}}
\newcommand{\MSE}{{\rm MSE}}
\newcommand{\muu}{\underline{\mu}}
\newcommand{\Mu}{\underline{M}}
\newcommand{\eps}{\varepsilon}
\newcommand{\thetu}{\underline{\theta}}
\newcommand{\betu}{\underline{\beta}}
\newcommand{\epsu}{\underline{\eps}}
\newcommand{\phiu}{\underline{\phi}}
\newcommand{\Xu}{\underline{X}}
\newcommand{\xu}{\underline{x}}
\newcommand{\Yu}{\underline{Y}}
\newcommand{\yu}{\underline{y}}
\newcommand{\piu}{\underline{\pi}}
\newcommand{\tauu}{\underline{\tau}}
\newcommand{\bu}{\underline{b}}
\newcommand{\nnu}{\underline{n}}
\newcommand{\wu}{\underline{w}}
\newcommand{\binco}[2]{\left( \begin{array}{c} #1 \\ #2 \end{array} \right) }

\begin{document}
\title{Durham/Newcastle Statistics Graduate Course \\
Introduction to OpenBUGS {\it via} BRugs}
\author{M.\ Farrow \\
School of Mathematics and Statistics \\
Newcastle University}
\date{January 2011}
\maketitle
\section{Introduction}
\subsection{Introduction}
The software BUGS ({\bf B}ayesian Inference {\bf U}sing {\bf G}ibbs {\bf S}ampling) was developed to allow users to specify models and priors, connect these with data and compute samples of unknowns from the posterior distribution using a Gibbs sampler (Spiegelhalter {\it et al}, 1995). Later a menu-driven version to run under MS Windows, called WinBUGS (Lunn {\em et al}, 2000) was developed. This eventually incorporated new features not found in the original, or `Classic'', BUGS. There is now also OpenBUGS, developed at the University of Helsinki, and various other implementations of the basic ``BUGS'' idea. In particular we will be using BRugs which is a R package which implements OpenBUGS within R. All of these use (apart from a few small differences) the same {\em Model Specification Language} and, in this part of this week's course, it is this language which is of particular interest.

Various resources are available at the following Web page.

\begin{verbatim}
                  http://www.mas.ncl.ac.uk/~nmf16/teaching/mas8303/
\end{verbatim}
This is the Web page for a module taught at Newcastle in which similar material is covered. References to the ``Module Web Page'' in what follows refer to this page.

The WinBUGS manual is available from the Module Web Page. The details of how you tell BRugs to do things are different from WinBUGS but the model specification language and many other features are the same.

A good overview is given by Lunn {\it et al} (2009).

\subsubsection{Potted history}
\begin{description}
\item[1989] BUGS Project starts at MRC Biostatistics Unit, Cambridge.
\item[1993] Unix version available. At first it used only Gibbs sampling (with adaptive rejection sampling where necessary) -- hence the name. Metropolis-Hastings was added later.
\item[1996] Project moved to Imperial College.
\item[1997] Windows version: WinBUGS.
\item[2004] Work began at Helsinki University on the open-source version, OpenBUGS. This led to BRugs -- interface which allows BUGS to run from R. (There are also interfaces to other statistical software now).
\item[2007] WinBUGS development work stopped. The future lies with OpenBUGS.
\item[2008] Competition: JAGS (``{\bf J}ust {\bf A}nother {\bf G}ibbs {\bf S}ampler'').
\end{description}


\subsection{Graphical models and the BUGS model specification language}
We will be using BRugs to do practical work. This is a R package which implements OpenBUGS. Models and priors are specified using the BUGS language. We will need to make a {\em model specification} file using the BUGS language. The BUGS language is closely linked to the idea of {\em directed acyclic graphs} (DAG). At each node, we only need to specify the conditional distribution of the variable given the variables at the parent nodes (direct predecessors).

\begin{description}
\item[Example 1] (Adapted from Gelfand {\it et al}, 1990).

Suppose we have $ k $ samples of observations and observation $ j $ in sample $ 
i $ is $ Y_{ij} $ where, given the values of $ \theta_{i} $ and $ \tau_{\eps}, $

\[ Y_{ij} \sim N(\theta_{i},\ \tau^{-1}_{\eps}) \]

for $ j = 1, \ldots, J, $  and, given the values of $ \mu $ and $ \tau^{-1}_{\theta}), $ 
\[ \theta_{i} \sim N(\mu,\ \tau^{-1}_{\theta}) \]
 (all independent).

(This could be a repeated measure model, for example).

Suppose we have independent priors for the three parameters:

\begin{eqnarray*}
\mu & \sim & N(\mu_{0}, \tau^{-1}_{0}) \\
\tau_{\theta} & \sim & {\rm Ga}(a_{1}, b_{1}) \\
\tau_{\eps} & \sim & {\rm Ga}(a_{2}, b_{2})
\end{eqnarray*}

where Ga$(a,b)$ stands for  a gamma 
distribution.

Figure \ref{fig:repmeas} shows the model with $ k=2 $ and $ J=3. $

\begin{figure}

\setlength{\unitlength}{0.1cm}

\begin{picture} (150,80)

\put (80,70) {\circle{20}}
\put (80,70) {\makebox(0,0){$ \tau_{\eps} $}} 

\put (30,50) {\circle{20}}
\put (30,50) {\makebox(0,0){$ Y_{11} $}} 
\put (50,50) {\circle{20}}
\put (50,50) {\makebox(0,0){$ Y_{12} $}} 
\put (70,50) {\circle{20}}
\put (70,50) {\makebox(0,0){$ Y_{13} $}} 
\put (90,50) {\circle{20}}
\put (90,50) {\makebox(0,0){$ Y_{21} $}} 
\put (110,50) {\circle{20}}
\put (110,50) {\makebox(0,0){$ Y_{22} $}} 
\put (130,50) {\circle{20}}
\put (130,50) {\makebox(0,0){$ Y_{23} $}}

\put (50,30) {\circle{20}}
\put (50,30) {\makebox(0,0){$ \theta_{1} $}} 
\put (110,30) {\circle{20}}
\put (110,30) {\makebox(0,0){$ \theta_{2} $}} 

\put (70,10) {\circle{20}}
\put (70,10) {\makebox(0,0){$ \mu $}}  
\put (90,10) {\circle{20}}
\put (90,10) {\makebox(0,0){$ \tau_{\theta} $}}

\put (65,15) {\vector(-1,1){9}}
\put (85,15) {\vector(-2,1){28}}
\put (75,15) {\vector(2,1){28}}
\put (95,15) {\vector(1,1){9}}

\put (45,35) {\vector(-1,1){9}}
\put (50,38) {\vector(0,1){4}}
\put (55,35) {\vector(1,1){9}}
\put (105,35) {\vector(-1,1){9}}
\put (110,38) {\vector(0,1){4}}
\put (115,35) {\vector(1,1){9}}

\put (77,63) {\vector(-1,-2){3}}
\put (75,65) {\vector(-3,-2){18}}
\put (73,68) {\vector(-3,-1){39}}

\put (83,63) {\vector(1,-2){3}}
\put (85,65) {\vector(3,-2){18}}
\put (87,68) {\vector(3,-1){39}}

\end{picture}

\caption{DAG: Repeated measures example.} 
\label{fig:repmeas}
\end{figure}    

The BUGS model specification is shown in Figure \ref{fig:repmeasBUG}.
 



\begin{figure}
\begin{verbatim}
model Gelfand

{for (i in 1:k)
    {for (j in 1:J)
        { y[i,j]~dnorm(theta[i],p.eps) 
         }
     theta[i]~dnorm(mu,p.theta)
     }

 mu~dnorm(mu0,p0)
 p.theta~dgamma(a1,b1)
 p.eps~dgamma(a2,b2)

 sigmasq.theta<-1/p.theta
 sigmasq.eps<-1/p.eps

 }
\end{verbatim}

\caption{BUGS model specification, repeated measures example.}
\label{fig:repmeasBUG}
\end{figure}


\item[Example 2] Patients in four groups are observed for various lengths of time. During this 
time tumours may develop. The dependent variable is the number of tumours 
observed for each patient. The mean number of tumours for a patient in group $ g 
$ is $ \lambda_{g} t = \exp(\beta_{g} + \ln t) $ where $ t $ is the time 
observed in weeks. Thus the parameters are $ \beta_{1}, \ldots , \beta_{4} $ 
where $ \beta_{g}= \ln(\lambda_{g}). $ There is no intercept here and the coefficient of
 $ \ln t $ is known to be 1.

 Suppose that we observe $ n $ patients (written as \verb@ N @ in the BUGS code). For each patient we have a group number $ g $ (\verb@ group@), the time $ t $ for which the patient was observed (\verb@ t @) and the number $ y $ of tumours observed (\verb@ y @).

Figure \ref{fig:tumours} shows some suitable BUGS code. 
\end{description}

Note that the code in Figures \ref{fig:repmeasBUG} and \ref{fig:tumours} are {\em not}  programs with commands to be executed. They are model specifications. We are defining the joint distribution of the unknowns and the data, mostly by specifying conditional distributions. For example
\begin{verbatim}
                 y[i]~dpois(mean[i])
\end{verbatim}
might be written in standard mathematical notation as
\[ Y_{i} \mid m_{i} \sim {\rm Po}(m_{i}). \]
The code \verb@ dpois @ represents the Poisson distribution and the symbol \verb@~@ has its usual meaning of ``has the following distribution.'' Similarly \verb@ dnorm @ stands for a normal distribution. Note however that the parameters are mean and {\em precision}, not mean and variance. So, in the example, we are saying that

\[ \beta_{g} \mid \mu \sim N(\mu,\ 0.1). \]


\begin{figure}
\begin{verbatim}

model tumours

{
 for (i in i:N)
     {y[i]~dpois(mean[i])
      mean[i]<-lambda[group[i]]*t[i]
      }

 for (g in 1:4)
      {lambda[g]<-exp(beta[g])
       beta[g]~dnorm(mu,10)
       }
 mu~dnorm(0,5)
 }
\end{verbatim}

\caption{BUGS code for the tumours example}
\label{fig:tumours}
\end{figure}

 

\newpage

\section{BUGS functions}

\begin{tabular}{lll}
Function & Usage & Definition \\ \hline
Complementary         & \verb@ cloglog(p)<-a+b*x @ & $ \log[-\log(1-p)]=a+bx $ \\
 \ \ \ \ \ log log    & \verb@ y<-cloglog(p) @     & $ y = \log[-\log(1-p)] $ \\
Logical equals        & \verb@ y<-equals(x,z) @    & $ y=1 $ if $ x=z $ \\
                      &                            & $ y=0 $ if $ x \neq z $ \\
Exponential           & \verb@ y<-exp(x) @         & $ y = e^{x} $ \\
Inner product         & \verb@ y<-inprod(a[],b[]) @ & $ y=\sum_{i} a_{i}b_{i} $ \\
Matrix inverse        & \verb@ y[,]<-inverse(x[,]) @ & $ y=x^{-1} $ \\
                      &                             & $ y,\ x $ both $ n \times n $ matrices \\
Natural logarithm     & \verb@ log(lambda)<-a+b*x @ & $ \log(\lambda)=a+bx $ \\
                      & \verb@ y<-log(x) @         & $ y=\log x $ \\
Log determinant       & \verb@ y<-logdet(x[,]) @   & $ y=\log|x| $ \\
                      &                            & $ x $ is a $ n \times n $ matrix \\
Log factorial         & \verb@ y<-logfact(x) @     & $ y=\log(x!) $ \\
Log(gamma function)   & \verb@ y<-loggam(x) @      & $ y=\log[ \Gamma(x)] $ \\
Logit                 & \verb@ y<-logit(p) @       & $ y=\log[p/(1-p)] $ \\
                      & \verb@ logit(p)<-a+b*x @   & $ \log[p/(1-p)]=a+bx $ \\
Maximum               & \verb@ c<-max(a,b) @       & $ c=\max(a,b) $ \\
Mean                  & \verb@ x.bar<-mean(x[]) @  & $ \bar{x}=\sum_{i} x_{i}/n $ \\
Minimum               & \verb@ c<-min(a,b) @       & $ c=\min(a,b) $ \\
Standard Gaussian     & \verb@ p<-phi(x) @         & $ p=\int_{-\infty}^{x} (2 \pi)^{-1/2} 
e^{ -t^{2}/2 }\ dt $ \\
\ distribution function &                           & i.e. $ p = \Phi(x) $  \\
Power                 & \verb@ z<-pow(x,y) @       & $ z=x^{y} $ \\
Probit                & \verb@ y<-probit(p) @      & $ y=\Phi^{-1}(p) $ \\
                      & \verb@ probit(p)<-a+b*x @  & $ \Phi^{-1}(p)=a+bx $ \\
Standard deviation    & \verb@ s<-sd(x[]) @        & $ s= \sqrt{ \sum_{i}(x_{i}-\bar{x})^{2}/n} $ \\
Square root           & \verb@ sigma<-sqrt(tau) @  & $ \sigma=\sqrt{\tau} $ \\
Unit step             & \verb@ y<-step(x) @        & $ y=0 $ if $ x<0 $ \\
                      &                            & $ y=1 $ if $ x \geq 0 $ \\
Sum                   & \verb@ x.sum<-sum(x[]) @   & $ x_{\rm sum}=\sum_{i} x_{i} $
\end{tabular}

\newpage

\section{BUGS distributions}

\begin{tabular}{lll}
Distribution & Usage & Definition\\ \hline
Bernoulli    & \verb@r~dbern(p)@ & $ f(r \mid p) = p^{r}(1-p)^{1-r}; $ \\
             &                     & $ r=0,1 $ \\
beta         & \verb@p~dbeta(a,b)@ & $ f(p \mid a,b) = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} p^{a-1}(1-p)^{b-1}; $ \\
             &                       & $ 0<p<1 $ \\
binomial     & \verb@r~dbin(p,n)@  & $ f(r \mid p,n) = \frac{n!}{r!(n-r)!} p^{r}(1-p)^{n-r}; $ \\
             &                       & $ r=0, \ldots ,n $ \\
categorical  & \verb@r~dcat(p[])@  & $ f(r \mid p_{1}, \ldots , p_{R})=p_{r}; $ \\
             &                       & $ r=1,2, \ldots , R $ where $ R={\rm dim}({\bf p})$ \\
chi-squared  & \verb@x~dchisq(k)@  & $ f(x \mid k)=2^{-k/2} x^{k/2-1} e^{-x/2}/\Gamma(\frac{k}{2}); $ \\
             &                       & $ x>0 $ \\
double       & \verb@x~ddexp(mu,tau)@ & $ f(x \mid \mu, \tau) = \frac{\tau}{2}e^{-\tau|x-\mu|}; $ \\
\ \ exponential &                       & $ -\infty<x<\infty $ \\
Dirichlet    & \verb@p[]~ddirch(alpha[])@ & $ f({\bf p} \mid \alpha) = \frac{ \Gamma(\sum_{i} \alpha_{i})}{\prod_{i} \Gamma(\alpha_{i})} \prod_{i} p_{i}^{\alpha_{i}-1}; $ \\
             &                          & $ 0<p_{i}<1,\ \ \sum_{i}p_{i}=1 $ \\
exponential  & \verb@x~dexp(lambda)@ & $ f(x \mid \lambda) = \lambda e^{-\lambda x}; $ \\
             &                         & $ x>0 $ \\
gamma        & \verb@x~dgamma(r,mu)@ & $ f(x \mid r,\mu) = \mu^{r}x^{r-1}e^{-\mu x}/\Gamma(r); $ \\
             &                         & $ x>0 $ \\
lognormal    & \verb@x~dlnorm(mu,tau)@ &  \\
             & \multicolumn{2}{c}{$ f(x \mid \mu,\tau)= \sqrt{\frac{\tau}{2 \pi}} x^{-1} \exp[-\frac{\tau}{2}(\log x - \mu)^{2}]; $ \ \  $ x>0 $ } \\
logistic     & \verb@x~dlogis(mu,tau)@ & \\
 & \multicolumn{2}{c}{ $ f(x \mid \mu, \tau)=\tau e^{\tau(x-\mu)}/(1+e^{\tau(x-\mu)})^{2}; $ \ \ $ -\infty<x<\infty $} \\
multivariate & \verb@x[]~dmnorm(mu[],T[,])@ & \\
\ \ normal   & \multicolumn{2}{c}{$ f( {\bf x} \mid \mu, {\bf T}) = (2 \pi)^{-N/2} | {\bf T} |^{1/2} \exp[-\frac{1}{2}({\bf x}- \mu)'{\bf T}({\bf x}-\mu)]; $} \\  
             &                           & $ -\infty< x_{i} < \infty $      \\
multinomial  & \verb@x[]~dmulti(p[],N)@ & $ f({\bf x} \mid {\bf p},N)= \frac{ ( \sum_{i} x_{i})!}{\prod_{i} x_{i} !} \prod_{i} p_{i}^{x_{i}} ; $ \\
   &  $ 0 < p_{i}<1, \ \sum_{i} p_{i}=1 $ & $ \sum_{i} x_{i}=N $  \\
negative     & \verb@x~dnegbin(p,r)@  & $ f(x \mid p,r)=\frac{(x+r-1)!}{x!(r-1)!}p^{r}(1-p)^{x}; $ \\
\ \ binomial &                        & $ x=0,1,2, \ldots $ \\
normal       & \verb@x~dnorm(mu,tau)@ & \\
 & \multicolumn{2}{c} {$ f(x \mid \mu,\tau)=\sqrt{ \frac{\tau}{2 \pi}} \exp[-\frac{\tau}{2}(x-\mu)^{2}]; $ \ \ $ -\infty<x<\infty $} \\
Pareto       & \verb@x~dpar(alpha,c)@  & $ f(x \mid \alpha,c)=\alpha c^{\alpha} x^{-(\alpha+1)}; $ \\
             &                         & $ x>c $ \\
Poisson      & \verb@r~dpois(lambda)@  & $ f(r \mid \lambda)= e^{-\lambda} \frac{\lambda^{r}}{r!}; $ \\
             &                          & $ r=0,1, \ldots $ \\
Student's $ t $ & \verb@x~dt(mu,tau,k)@  & \\
   & \multicolumn{2}{c}{ $ f(x \mid \mu,\tau,k)= \frac{ \Gamma([k+1]/2)}{\Gamma(k/2)} \sqrt{ \frac{\tau}{k \pi}}[1+ \frac{\tau}{k}(x-\mu)^{2}]^{-(k+1)/2} ; $ } \\
   & & $ -\infty<x<\infty $  \\
uniform         & \verb@x~dunif(a,b)@  & $ f(x \mid a,b)=1/(b-a); $ \\
                &                      & $ a<x<b $ \\
Weibull         & \verb@x~dweib(v,lambda)@ & \\
 & \multicolumn{2}{c}{$ f(x \mid v, \lambda)=v \lambda x^{v-1} \exp(-\lambda x^{v}); $}  \\
                &                       & $ x>0 $ \\
Wishart         & \verb@x[,]~dwish(R[,],k)@ & \\
 & \multicolumn{2}{c}{$ f({\bf x} \mid {\bf R},k) \propto | {\bf R} |^{k/2} | {\bf x} |^{(k-p-1)/2} \exp(-\frac{1}{2}{\rm tr}[ {\bf Rx}]); $} \\
         &  \multicolumn{2}{c}{ $ {\bf x} $ symmetric and positive definite}
\end{tabular}

\newpage

\section{Some BRugs (and other) Commands}

\begin{verbatim}
library(BRugs)


modelCheck("modelfile.txt")

modelData("datafile.txt")

modelCompile()
modelCompile(2)

modelInits("initsfile.txt")

modelGenInits()

modelUpdate(5000)

samplesSet("theta")
samplesSet(c("alpha","beta","gamma"))

samplesHistory("theta")

samplesDensity("theta")

samplesStats("theta")
samplesStats(c("alpha","beta","gamma"))

thetavalues<-samplesSample("theta")

iter<-1:length(thetavalues)
plot(iter,thetavalues,type="l",xlab="Iteration",ylab=expression(theta))

thetadens<-density(theta,adjust=1.5)
plot(thetadens$x,thetadens$y,type="l",xlab=expression(theta),ylab="Density")

\end{verbatim}



\section*{References}
\begin{list}{ } {\setlength{\itemindent}{-1cm}}
\item Box, G.E.P.\ and Jenkins, G.M., 1976. {\em Time Series Analysis: Forecasting and Control}, Revised edition. San Francisco: Holden-Day.

\item Cooper, C., Shah, S., Hand, D.J., Adams, J., Compston, J., Davie, M.\ and Woolf,A., 1991. Screening for vertebral osteoporosis using individual risk factors. {\em Osteoporosis International}, {\bf 2}, 48-53.

\item Cowburn, G.J., 2003. 
{\em Bayesian mixture modelling with application to road traffic flow.} 
Ph.D. thesis, University of Sunderland.

\item Cowburn, G.J.\ and Farrow, M., 2007.
Mixtures and diagnostic plots in modelling road traffic vehicle headways.
Statistical Modelling, \underline{7}, 73-89.

\item Davies, O.L.\ and Goldsmith, P.L.\ (eds), 1972. {\em Statistical Methods in Research and Production}, 4th edition. Edinburgh: Oliver and Boyd.

\item Freeman, D.H., 1987. {\em Applied Categorical Data Analysis}. New York: Marcel Dekker.

\item Gelfand, A.E., Hills, S.E., Racine-Poon, A.\ and Smith, A.F.M., 1990. Illustration of Bayesian inference in normal data models using Gibbs sampling. {\em Journal of the American Statistical Association}, {\bf 85}, 972-85.

\item Krzanowski, W.J., 1988. {\em Principles of Multivariate Analysis}. Oxford: Oxford University Press.

\item Lunn, D.J., Thomas, A., Best, N.G., and Spiegelhalter, D.,J., 2000. WinBUGS -- a Bayesian modelling framework: concepts, structure, and extensibility. {\em Statistics and Computing}, {\bf 10}, 325-337.

\item Lunn, D.J., Spiegelhalter, D.J., Thomas, A.\ and Best, N.G., 2009.
The BUGS project: Evolution, critique and future directions.
{\em Statistics in Medicine}, {\bf 28}, 3049-3067.

\item Phillips, D.P., 1978. Airplane accident fatalities increase just after newspaper stories about murder and suicide. {\em Science}, {\bf 201}, 748-750.

\item Spiegelhalter, D.J., Thomas, A., Best, N.G.\ and Gilks, W.R., 1995. {\em BUGS: Bayesian inference Using Gibbs Sampling, Version 0.50}. MRC Biostatistics Unit, Cambridge.

\end{list}


%\item Azzalini, A.\ and Bowman, A.W., 1990. A look at some data on the Old Faithful geyser. {\it Journal of the Royal Statistical Society, Series C}, {\bf 39}, 357-366.

%\item Breslow, N.E.\ and Clayton, D.G., 1993. Approximate inference in generalised linear mixed models. {\em Journal of the American Statistical Association}, {\bf 88}, 9-25.

%\item Hand, D.J., Daly, F., Lunn, A.D., McConway, K.J.\ and Ostrowski, E., 1994. {\it A Handbook of Small Data Sets}. Chapman and Hall, London.

%\item Open University, 1983. {\em MDST242 Statistics in Society, Unit C3: Is my child normal?}. Milton Keynes: The Open University.

%\item Mosteller, F.\ and Tukey, J.W., 1977. {\em Data analysis and regression}. Reading, Massachusetts: Addison-Wesley.

%\item Pearl,J., Geiger,D.\ and Verma,T., 1990, The logic of influence diagrams, In {\em Influence Diagrams, Belief Nets and Decision Analysis} (R.M.Oliver, J.Q.Smith eds.) New York: Wiley.


%\item Snedecor, G.W.\ and Cochran, G.C., 1967. {\em Statistical Methods}, 6th edition. Iowa State University Press, Ames, Iowa.

%\item Sokal, R.R.\ and Rohlf, F.J., 1981. {\em Biometry}, 2nd edition. San Francisco: W.H.\ Freeman.


%\item Thall, P.F.\ and Vail, S.C., 1990. Some covariance models for longitudinal count data with overdispersion. {Biometrics}, {\bf 46}, 657-71.

%\item Till, R., 1974. {\em Statistical Methods for the Earth Scientist}. London: MacMillan.





\end{document}
