%\documentclass{article}
\documentclass[article,shortnames,nojss]{jss}
\usepackage{amssymb}
\usepackage{natbib}
\usepackage{graphicx}
\usepackage{color}
\usepackage{verbatim}
\usepackage{hyperref}
\usepackage{url}
%\usepackage{ae}



\setlength{\unitlength}{1cm}
\newcommand{\ex}[1]{\ensuremath{\mathbb{E}[#1]}}
\newcommand{\var}[1]{\ensuremath{\mathrm{Var}[#1]}}
\newcommand{\cov}[1]{\ensuremath{\mathrm{Cov}[#1]}}
\newcommand{\corr}[1]{\ensuremath{\mathrm{Corr}[#1]}}
\newcommand{\bd}[1]{\ensuremath{\mbox{\boldmath $#1$}}}

%\VignetteIndexEntry{Vignette for CARBayes  package.}


%% almost as usual
\author{Duncan Lee\\University of Glasgow }
\title{\pkg{CARBayes} version 4.6: An \proglang{R} Package for Spatial Areal Unit Modelling  with Conditional Autoregressive Priors}

%% for pretty printing and a nice hypersummary also set:
\Plainauthor{Duncan Lee} %% comma-separated
\Plaintitle{CARBayes version 4.6: An R package for spatial areal unit modelling  with conditional autoregressive priors} %% without formatting
\Shorttitle{\pkg{CARBayes}: Bayesian Conditional Autoregressive modelling} %% a short title (if necessary)

%% an abstract and keywords
\Abstract{
This is a vignette for the \proglang{R} package \pkg{CARBayes} version 4.6, and is an updated version of a paper in the Journal of Statistical Software in 2013 Volume 55 Issue 13 by the same author. This vignette describes the class of models that can be implemented in CARBayes, and gives 3 worked examples of spatial data analysis using the package. Version 4.6 has one minor change from version 4.5, namely the removal of errors in the R code in the vignette. However, version 4.5 has more major changes from version 4.4, namely the introduction of an MCAR model for multivariate spatial data, for details see Section 2 of this vignette. Additionally, when fitting a model with no spatial autocorrelation (using the \code{S.CARleroux()} or \code{MVS.CARleroux()} functions and setting the arguments, \code{fix.rho=TRUE, rho=0}) the neighbourhood matrix $\mathbf{W}$ is not part of the model. However, for compatability with the case that $\rho$ is not fixed at zero, a fake $\mathbf{W}$ matrix, that is a $K\times K$ matrix of zeros, should be entered, where $K$ is the number of spatial units.
}
\Keywords{Bayesian inference, conditional autoregressive priors, \proglang{R} package \pkg{CARBayes}}
\Plainkeywords{Bayesian inference, conditional autoregressive priors, R package CARBayes}

\Address{
  Duncan Lee\\
  School of Mathematics and Statistics\\ 
  15 University Gardens\\  
  University of Glasgow\\
  Glasgow\\ 
  G12 8QQ, Scotland\\
  E-mail: \email{Duncan.Lee@glasgow.ac.uk}\\
  URL: \url{http://www.gla.ac.uk/schools/mathematicsstatistics/staff/duncanlee/}
}
%% The address of (at least) one author should be given
%% in the following format:

\begin{document}
\SweaveOpts{concordance=TRUE}

%%%%%%%%%%%%%%
%%%% Section 1
%%%%%%%%%%%%%%

<<echo=FALSE>>=
options(prompt = "R> ")
@



\section{Introduction}
Data relating to a set of non-overlapping spatial areal units are prevalent in many fields, including agriculture (\cite{besag1999}), ecology (\cite{brewer2007}), education (\cite{wall2004}), epidemiology (\cite{lee2011}) and image analysis (\cite{gavin1997}).  There are numerous motivations for modelling such data, including ecological regression  (see \cite{wakefield2007} and \cite{lee2009}), disease mapping (see \cite{green2002} and \cite{lee2011}) and Wombling (see \cite{lu2007}, \cite{ma2007}). The set of areal units on which data are recorded can form a regular lattice or differ largely in both shape and size, with examples of the latter including the set of electoral wards or census tracts corresponding to a city or country. In either case such data typically exhibit spatial autocorrelation, with observations from areal units close together tending to have similar values. A proportion of this spatial autocorrelation may be modelled by known covariate risk factors in a regression model, but it is common for  spatial structure to remain in the residuals after accounting for these covariate effects. This residual spatial autocorrelation can be induced by a number of factors, and violates the assumption of independence that is common in many regression models. One possible cause is unmeasured confounding, which occurs when an important spatially autocorrelated covariate is either unmeasured or unknown. The spatial structure in this covariate induces spatial autocorrelation into the response, which hence cannot be accounted for in a regression model. Other possible causes of residual spatial autocorrelation are neighbourhood effects, where subjects behaviour is influenced by that of neighbouring subjects, and grouping effects, where subjects choose to be close to similar subjects.\\

The most common remedy for this residual autocorrelation is to augment the linear predictor with a set of spatially autocorrelated random effects, as part of a Bayesian hierarchical model. These random effects are typically represented with a conditional autoregressive (CAR, \cite{besag1991}) prior,  which induces spatial autocorrelation through the adjacency structure of the areal units. A number of CAR priors have been proposed in the literature, including the intrinsic and Besag-York-Molli\'{e} (BYM) models (both \cite{besag1991}), as well as alternatives  developed by \cite{leroux1999} and \cite{stern1999}. However, the CAR priors listed above force the random effects to exhibit a single global level of spatial autocorrelation, ranging from independence through to strong spatial smoothness. Such a uniform level of spatial autocorrelation for the entire region maybe unrealistic for real data, which instead may exhibit sub-regions of spatial autocorrelation separated by discontinuities. Such localised spatial autocorrelation may occur where rich and poor communities live side-by-side, and in this context the response variable is likely to evolve smoothly within each community with a sudden change in its value at the border where the two communities meet. However, covariate data quantifying this localised strucure may not be available, meaning that it has to be modelled by the random effects. A number of approaches have been proposed for extending the class of CAR priors to deal with localised spatial smoothing amongst the random effects, including papers by \cite{lawson2002} (combining the intrinsic model with a `jump' component for discontinuities), \cite{brewer2007} (variable smoothing via a spatially varying variance),  \cite{lu2007} (modelling the adjacency structure of the areal units using logistic regression),  \cite{reich2008} (varible smoothing via a spatially varying variance in a spatio-temporal setting), \cite{lee2012} (modelling the partial correlation between random effects in adjacent areal units as a function of their dissimilarity), and \cite{lee2014} (localised smoothing in an ecological regression context).\\


The models described above are typically implemented in a Bayesian setting, where inference is based on Markov chain Monte Carlo (McMC) simulation. The most commonly used software to implement this class of models is the \proglang{BUGS} project (\cite{lunn2009}, \proglang{WinBUGS} and \proglang{OpenBUGS}), which has in-built functions \code{car.normal} and \code{car.proper} to implement the intrinsic, BYM and \cite{stern1999} models, as well as allowing users to write code to implement their own spatial random effects models. The intrinsic and BYM models can also be implemented in \proglang{BayesX} (\cite{belitz2009}), while the \proglang{R} software (\citealp{R}) has packages \pkg{hSDM} (\citealp{hSDM}), \pkg{spatcounts} (\citealp{schabenberger2009}) and \pkg{spdep} (\citealp{bivand2013}) can also implement a restricted set of CAR models. CAR models can also be implemented in \proglang{R} using Integrated Nested Laplace Approximations (INLA, \emph{http://www.r-inla.org/}), using the package \pkg{INLA} (\citealp{rue2009}).\\


However, each of these software packages either can only fit a limited set of CAR models or require a degree of programming to implement them, which is the motivation for creating the \proglang{R} package \pkg{CARBayes} (original reference \citealp{lee2013}). The main advantage of this package is its ease of use in fitting CAR models, because: (1) the spatial adjacency information is easy to specify as a binary neighbourhood matrix; and (2) given the neighbourhood matrix the models can be implemented by a single function call in \proglang{R}. In addition, \pkg{CARBayes} can implement a much wider class of spatial areal unit models than is possible using the \proglang{R} packages listed above, as the response data can follow binomial, Gaussian or Poisson distributions, while a range of CAR priors have been implemented for the random effects.\\

The aim of this vignette is to present the software \pkg{CARBayes}, by outlining the class of models that it can implement and illustrating its use by means of 3 worked examples. The remainder of this paper is organised as follows. Section two outlines the general Bayesian hierarchical model that can be implemented in the \pkg{CARBayes} package, while Section three gives details about the software.  Sections four to six give three worked examples of using the software, including how to create the neighbourhood matrix and produce spatial maps of the results. Finally, Section 7 contains a concluding discussion, and outlines areas for future development.


%%%%%%%%%%%%%%
%%%% Section 2
%%%%%%%%%%%%%%
\section{Spatial models for areal unit data}
This section outlines the class of spatial generalised linear mixed models for areal unit data that can be implemented in the \pkg{CARBayes} software. Inference for all models is set in a Bayesian framework, and is based on Markov chain Monte Carlo (McMC) simulation. The majority of the models in \pkg{CARBayes} relate to univariate spatial data and are described in Section 2.1, but the function \code{MVS.CARleroux()} fits a model to multivariate spatial data and is described in Section 2.2.


\subsection{Univariate spatial data models}
The study region $\mathcal{S}$ is partitioned into $K$ non-overlapping areal units $\mathcal{S}=\{\mathcal{S}_{1},\ldots,\mathcal{S}_{K}\}$, which are linked to a corresponding set of responses $\mathbf{Y}=(Y_{1},\ldots,Y_{K})$, and a vector of known offsets $\mathbf{O}=(O_{1},\ldots,O_{K})$. Missing, \code{NA}, values are allowed in the repsonse $\mathbf{Y}$ except for the \code{S.CARlocalised()} function which does not allow them due to model complexity and corresponding poor predictive performance. The spatial pattern in the response is modelled by a matrix of covariates $\mathbf{X}=(\mathbf{x}_{1},\ldots,\mathbf{x}_{K})$ and  a spatial structure component $\bd{\psi}=(\psi_{1},\ldots,\psi_{K})$, the latter of which is included to model any spatial autocorrelation that remains in the data after the covariate effects have been accounted for. The vector of covariates for areal unit $\mathcal{S}_{k}$ are denoted by $\mathbf{x}_{k}=(1, x_{k1},\ldots,x_{kp})$, the first of which corresponds to an intercept term. The general model class that \pkg{CARBayes} can implement is a generalised linear mixed model for spatial areal unit data which is given by


\begin{eqnarray}
Y_{k}|\mu_{k}&\sim&f(y_{k}|\mu_{k},\nu^{2})~~~~\mbox{for }k=1,\ldots,K,\label{equation likelihood}\\
g(\mu_{k})&=& \mathbf{x}_{k}^{\top}\bd{\beta} + O_k + \psi_{k} ,\nonumber\\
\bd{\beta}&\sim&\mbox{N}(\bd{\mu}_{\beta}, \bd{\Sigma}_{\beta}).\nonumber
\end{eqnarray}


The responses $Y_{k}$ come from an exponential family of distributions $f(y_{k}|\mu_{k}, \nu^{2})$, and in  \pkg{CARBayes} these can be the binomial, Gaussian or Poisson families. The expected value of $Y_{k}$ is denoted by $\E(Y_{k})=\mu_{k}$, while $\nu^{2}$ is an additional scale parameter that is required if the Gaussian family is used. The expected values of the responses are related to the linear predictor via an invertible link function $g(.)$, which in this software is one of the logit (binomial family), identity (Gaussian family) and natural log (Poisson family) functions. Thus \pkg{CARBayes} can fit the following likelihood models:

\begin{itemize}
\item \textbf{Binomial - } $Y_{k}~\sim~\mbox{Binomial}(n_{k}, \theta_{k})$ and $\ln(\theta_{k}/(1-\theta_{k}))~=~\mathbf{x}_{k}^{\top}\bd{\beta} + O_{k} + \psi_{k}$. 

\item \textbf{Gaussian - } $Y_{k}~\sim~\mbox{N}(\mu_{k}, \nu^{2})$ and $\mu_{k}~=~\mathbf{x}_{k}^{\top}\bd{\beta} + O_{k} + \psi_{k}$. 

\item \textbf{Poisson - } $Y_{k}~\sim~\mbox{Poisson}(\mu_{k})$ and $\ln(\mu_{k})~=~\mathbf{x}_{k}^{\top}\bd{\beta} + O_{k} + \psi_{k}$. 
\end{itemize}

In the binomial model above $n_k$ is the number of trials in the $k$th area, while $\theta_k$ is the probability of success in a single trial. The vector of regression parameters are denoted by $\bd{\beta}=(\beta_{1},\ldots,\beta_{p})$, and non-linear covariate effects can be incorporated into the above model by including natural cubic spline or polynomial basis functions of the covariates in $\mathbf{X}$. A multivariate Gaussian prior is assumed for $\bd{\beta}$, and the mean $\bd{\mu}_{\beta}$ and diagonal variance matrix $\bd{\Sigma}_{\beta}$ can be chosen by the user. Default values specified by the software are a constant zero-mean vector and diagonal elements of  $\bd{\Sigma}_{\beta}$ equal to 1000. The scale parameter $\nu^2$ for the Gaussian likelihood is assigned a conjugate inverse-gamma prior distribution,  where the default specification is $\nu^{2}\sim\mbox{Inverse-Gamma}(1, 0.01)$. \pkg{CARBayes} can implement a number of different spatial random effects models for $\bd{\psi}$, and they are summarised below.


\begin{itemize}
\item \code{S.CARbym()} - fits the convolution or Besag-York-Mollie (BYM) CAR model outlined in \cite{besag1991}.

\item \code{S.CARleroux()} - fits the CAR model proposed by \cite{leroux1999}. This model can also fit the intrinsic CAR model proposed by \cite{besag1991}, as well as a model with independent random effects. 

\item \code{S.CARdissimilarity()} - fits the localised spatial autocorrelation model proposed by \cite{lee2012}.

\item \code{S.CARlocalised()} - fits the localised spatial autocorrelation model proposed by \cite{lee2015}.
\end{itemize}


The spatial structure component $\psi$ in each of the above models includes a set of random effects $\bd{\phi}=(\phi_{1},\ldots,\phi_{K})$, which come from a conditional autoregressive model. These models are a special case of a Gaussian Markov Random Field (GMRF), and can be written in the general form $\bd{\phi}\sim\mbox{N}(\mathbf{0},\tau^{2}\mathbf{Q}^{-1})$, where $\mathbf{Q}$ is a precision matrix that may be singular (intrinsic model). This matrix controls the spatial autocorrelation structure of the random effects, and is based on a non-negative symmetric $K\times K$ neighbourhood or weight matrix $\mathbf{W}$. A binary specification based on geographical contiguity is most commonly used, where $w_{kj}=1$ if areal units $(\mathcal{S}_k, \mathcal{S}_j)$ share a common border (denoted $k\sim j$), and is zero otherwise. This specification forces $(\phi_k, \phi_j)$ relating to geographically adjacent areas (that is where $w_{kj}=1$) to be autocorrelated, whereas random effects relating to non-contiguous areal units are conditionally independent given the values of the remaining random effects. A binary specification is not necessary in \pkg{CARBayes}, as the only requirement is that $\mathbf{W}$ is non-negative and symmetric. However, each area must have at least one positive element $w_{kj}$, so that the row sums of the $\mathbf{W}$ matrix are not allowed to contain zeros. CAR priors are commonly specified as a set of $K$ univariate full conditional distributions $f(\phi_{k}|\bd{\phi}_{-k})$ for $k=1,\ldots,K$ (where $\bd{\phi}_{-k}=(\phi_{1},\ldots,\phi_{k-1},\phi_{k+1},\ldots,\phi_{K})$), rather than via the multivariate specification described above. We now outline the four models that \pkg{CARBayes} can fit below.\\


\subsubsection{Globally smooth CAR models}
\code{S.CARbym()}\\
The  convolution or Besag-York-Mollie (BYM) CAR model outlined in \cite{besag1991} was the first to be proposed, and contans two sets of random effects, spatially autocorrelated and independent. It has the following form. 


\begin{eqnarray}
\psi_k&=&\phi_{k} + \theta_k,\label{equation intrinsic}\\
\phi_{k}| \bd{\phi}_{-k},\mathbf{W},\tau^2&\sim&\mbox{N}\left(\frac{\sum_{i=1}^{K}w_{ki}\phi_{i}}{\sum_{i=1}^{K}w_{ki}},~
\frac{\tau^{2}}{\sum_{i=1}^{K}w_{ki}}\right),\nonumber\\
\theta_k&\sim&\mbox{N}(0, \sigma^2),\nonumber\\
\tau^{2}, \sigma^2&\sim&\mbox{Inverse-Gamma}(a, b).\nonumber
\end{eqnarray}

The $\bd{\theta}=(\theta_{1},\ldots,\theta_{K})$ random effects are independent with zero mean and a constant variance, while the spatial autocorrelation is modelled via $\bd{\phi}$. For this latter set the conditional expectation is the average of the random effects in neighbouring areas, while the conditional variance is inversely proportional to the number of neighbours. The latter is appropriate because if the random effects are strongly spatially autocorrelated, then the more neighbours an area has the more information there is from its neighbours about the value of its random effect. In common with the other variance parameters the default prior specification for ($\tau^{2},\sigma^2$) has $(a=1, b=0.01)$. This model is known as the BYM or convolution model, and is the most commonly used conditional autoregressive model in practice. However, it requires two random effects to be estimated for each data point, whereas only their sum is identifiable from the data. Therefore in \pkg{CARBayes} only McMC samples for the sum $\psi_k=\phi_k+\theta_k$ are returned to the user. \\

\code{S.CARleroux()}\\
As a result of this limitation \cite{leroux1999} proposed an alternative CAR prior for modelling varying strengths of spatial autocorrelation, using only a single set of random effects. Their model  is given by

\begin{eqnarray}
\psi_k&=&\phi_k\label{equation leroux}\\
\phi_{k}| \bd{\phi}_{-k},\mathbf{W},\tau^2, \rho&\sim&\mbox{N}\left(\frac{\rho\sum_{i=1}^{K}w_{ki}\phi_{i}}{\rho\sum_{i=1}^{K}w_{ki} + 1-\rho},~
\frac{\tau^{2}}{\rho\sum_{i=1}^{K}w_{ki} + 1-\rho}\right),\nonumber\\
\tau^{2}&\sim&\mbox{Inverse-Gamma}(a, b),\nonumber\\
\rho&\sim&\mbox{Uniform}(0,1),\nonumber
\end{eqnarray}

where $\rho$ is a spatial autocorrelation parameter.  A uniform prior on the unit interval is specified for $\rho$, while the usual inverse gamma prior is adopted for $\tau^{2}$. The value of $\rho$ can be fixed in the interval $[0,1]$ when fitting the model rather than being estimated from the data, which in particular allows 2 special cases to be fitted. Fixing $\rho=1$ results in the intrinsic CAR model (corresponding to the model for $\bd{\phi}$ in (\ref{equation intrinsic})), which was available via a separate function (\code{S.CARiar()}) in older versions of \pkg{CARBayes}. In contrast, fixing $\rho=0$ in the above corresponds to independence (corresponding to the model for $\bd{\theta}$ in (\ref{equation intrinsic})), which was available via a separate function (\code{S.independent()}) in older versions of \pkg{CARBayes}. In this case the neighbourhood matrix $\mathbf{W}$ is not part of the model, so enter a fake $\mathbf{W}$ matrix that is a $K\times K$ matrix of zeros, where $K$ is the number of spatial areal units. These global CAR models were compared in a review by \cite{lee2011}, who concluded  that the model proposed by \cite{leroux1999} was the most appealing from both theoretical and practical  standpoints.\\


\subsubsection{Locally smooth CAR models}
The CAR priors described above enforce a single global level of spatial smoothing for the set of random effects, which for model (\ref{equation leroux}) is controlled by $\rho$. This is illustrated by the partial autocorrelation structure implied by that model, which for $(\phi_k, \phi_j)$ is given by

\begin{equation}
\mbox{COR}(\phi_{k},\phi_{j}|\bd{\phi}_{-kj}, \mathbf{W}, \rho)~=~\frac{\rho w_{kj}}{\sqrt{(\rho\sum_{i=1}^{K}w_{ki} + 1-\rho)(\rho\sum_{i=1}^{K}w_{ji} + 1-\rho)}}\label{equation partialcorrelation}.
\end{equation}


For non-neighbouring areal units (where $w_{kj}=0)$ the random effects are conditionally independent, while for neighbouring areal units (where $w_{kj}=1$) their partial autocorrelation is controlled by $\rho$. However, this representation of spatial smoothness is likely to be overly simplistic in practice, as the random effects surface is likely to include sub-regions of smooth evolution as well as boundaries where abrupt step changes occur.  A number of approaches have been proposed for modelling this localised spatial autocorrelation, and \pkg{CARBayes} can implement the proposals of \cite{lee2012} and \cite{lee2015} which are both summarised below.\\


\code{S.CARdissimilarity()}\\ 
The paper by \cite{lee2012} proposes a method for capturing such localised spatial autocorrelation, including the identification of boundaries in the random effects surface. The underlying idea is to model the elements of $\mathbf{W}$ corresponding to geographically adjacent areal units as binary random quantities, rather than assuming they are fixed at one. Conversely, if areal units $(\mathcal{S}_{k},\mathcal{S}_{j})$ do not share a common border then $w_{kj}$ is fixed at zero. From (\ref{equation partialcorrelation}), it is straightforward to see that if $w_{kj}$ is estimated as one then $(\phi_k, \phi_j)$ are spatially autocorrelated, and are smoothed over in the modelling process. In contrast, if $w_{kj}$ is estimated as zero then no smoothing is imparted between $(\phi_k, \phi_j)$, as they are modelled as conditionally independent. In this case a boundary is said to exist in the random effects surface between areal units $(\mathcal{S}_{k},\mathcal{S}_{j})$. We note that if covariates are excluded from (\ref{equation likelihood}) then any boundaries identified also relate to the mean surface $\bd{\mu}=(\mu_1,\ldots,\mu_K)$ in the absence of an offset term, because it has the same spatial structure as the random effects as $g(\mu_k)=\beta_0 + \phi_k$.\\

The model proposed by \cite{lee2012}  is based on the Poisson log-linear specification of (\ref{equation likelihood}) and the CAR prior (\ref{equation leroux}), with the restriction that $\rho$ is fixed at 0.99 (\pkg{CARBayes} can also fit binomial and Gaussian specifications). This restriction on $\rho$ was made by \cite{lee2012} to ensure that the random effects exhibit strong spatial smoothing globally, which can be altered locally by estimating $\{w_{kj}|k\sim j\}$. They model each $w_{kj}$ as a function of the dissimilarity between areal units $(\mathcal{S}_{k},\mathcal{S}_{j})$, because large differences in the response are likely to occur where neighbouring populations are very different. This dissimilarity is captured by $q$ non-negative dissimilarity metrics $\mathbf{z}_{kj}=(z_{kj1},\ldots,z_{kjq})$, which could include social or physical factors, such as the absolute difference in smoking rates, or the proportion of the shared border that is blocked by a physical barrier (such as a river or railway line) and cannot be crossed. Using these measures of dissimilarity, $\{w_{kj}|k\sim j\}$ are collectively modelled as

\begin{eqnarray}
w_{kj}(\bd{\alpha})&=&\left\{ \begin{array}{ll}
1&\mbox{if }\exp(-\sum_{i=1}^{q}z_{kji}\alpha_{i})\geq0.5 \mbox{ and } k\sim j\\
0&\mbox{otherwise}\\\end{array}\right.,\label{equation binary neighbourhood}\\
\alpha_{i}&\sim&\mbox{Uniform}(0, M_{i})\hspace{1cm}\mbox{for $i=1,\ldots,q$.}\nonumber
\end{eqnarray}

The $q$ regression parameters $\bd{\alpha}=(\alpha_{1},\ldots,\alpha_{q})$ determine the effects of the dissimilarity metrics on $\{w_{kj}|k\sim j\}$, and if $\alpha_i<-\ln(0.5)/\max\{z_{kji}\}$, then the $i$th dissimilarity metric has not solely identified any  boundaries because $\exp(-\alpha_i z_{kji})>0.5$ for all $k\sim j$. The upper limits $M_{i}$ for the priors for $\alpha_i$ depend on the distribution of $z_{kji}$, and are chosen to be weakly informative and fixed in the software. The aim of \cite{lee2012} was to identify the locations of any boundaries (abrupt step changes) in disease risk surfaces, so the available covariates were used to construct dissimilarity metrics rather than being incorporated into the linear predictor. In contrast, if the aim of the analysis was to explain the spatial pattern in the response, then covariates would be included in (\ref{equation likelihood}), and only metrics directly describing the dissimilarity between two areas, such as the existence of a physical boundary or the distance between the areas centroids, would be included in (\ref{equation binary neighbourhood}).\\

\code{S.CARlocalised()}\\
The above approach allows for localised spatial autocorrelation by modelling the elements in $\mathbf{W}$, which can reduce the partial autocorrelations between certain pairs of adjacent random effects. An alternative is to augment the set of spatially smooth random effects with a piecewise constant intercept or cluster model, thus allowing large jumps in the mean surface between adjacent areal units in different clusters. One such model was proposed by \cite{lee2015}, and partitions the $K$ areal units into a maximum of $G$ clusters each with their own intercept term $(\lambda_{1},\ldots,\lambda_{G})$. The  model is given by


\begin{eqnarray}
\psi_k & = & \phi_{k} + \lambda_{Z_{k}},\label{equation cluster}\\
\phi_{k}| \bd{\phi}_{-k}, \mathbf{W}, \tau^2&\sim&\mbox{N}\left(\frac{\sum_{i=1}^{K}w_{ki}\phi_{i}}{\sum_{i=1}^{K}w_{ki}},\frac{\tau^{2}}{\sum_{i=1}^{K}w_{ki}}\right),\nonumber\\
\tau^{2}&\sim&\mbox{Inverse-Gamma}(a, b),\nonumber\\
\lambda_i & \sim & \mbox{Uniform}(\lambda_{i-1}, \lambda_{i+1})\hspace{1cm}\mbox{for }i=1,\ldots,G,\nonumber\\
f(Z_{k})&=&\frac{\exp(-\delta(Z_{k}-G^{*})^{2})}{\sum_{r=1}^{G}\exp(-\delta(r-G^{*})^{2})},\nonumber\\
\delta&\sim&\mbox{Uniform}(1, M).\nonumber
\end{eqnarray}

The cluster means $(\lambda_{1},\ldots,\lambda_{G})$ are ordered so that $\lambda_{1}<\lambda_{2}<\ldots<\lambda_{G}$, which prevents the label switching problem common in mixture models. Thus in the above $\lambda_{0}=-\infty$ and $\lambda_{G+1}=\infty$.  Area $k$ is assigned to one of the $G$ intercepts by $Z_{k}\in\{1,\ldots,G\}$, and $G$ is the maximum number of different intercept terms. Here we penalise $Z_{k}$ towards the middle intercept value, so that the extreme intercept classes (e.g. $1$ or $G$) may be empty. This is achieved by the penalty term $\delta(Z_{k}-G^{*})^{2}$ in the prior for $Z_{k}$, where $G^{*}=(G+1)/2$ if $G$ is odd and $G^{*}=G/2$ if $G$ is even is the middle of the intercept terms. A weakly informative uniform prior is specified for the penalty parameter $\delta\sim\mbox{Uniform}(1,M)$ (by default $M=10$), so that the data play the dominant role in estimating its value. Again $(a=1, b=0.01)$ are the default values in the Inverse-Gamma prior for $\tau^2$. We note that we do not put any spatial smoothing constraints on the indicator vector $\mathbf{Z}$, because similar random effect values can occur at opposite ends of the study region. The clustering is thus inherently non-spatial, and the spatial autocorrelation is accounted for by the CAR prior for $\bd{\phi}$.  The \pkg{CARBayes} package can fit binomial and Poisson variants of this model. Note, a Gaussian likelihood is not allowed because of a lack of identifiability among the parameters.


\subsection{Multivariate spatial data models}
The study region $\mathcal{S}$ is again partitioned into $K$ non-overlapping areal units $\mathcal{S}=\{\mathcal{S}_{1},\ldots,\mathcal{S}_{K}\}$, which are linked to a corresponding set of responses $\mathbf{Y}=(\mathbf{Y}_{1},\ldots,\mathbf{Y}_{K})$, 
where $\mathbf{Y}_{k}=(Y_{k1},\ldots,Y_{kJ})$. Thus for each areal unit there are $J$ different response variables being modelled, leading to both spatial autocorrelation and between variable correlation. The general likelihood model allowed for these data is

\begin{eqnarray}
Y_{kj}|\mu_{kj}&\sim&f(y_{kj}|\mu_{kj})~~~~\mbox{for }k=1,\ldots,K,~~~~ j=1,\ldots,J,\label{equation MV likelihood}\\
g(\mu_{kj})&=& \mathbf{x}_{kj}^{\top}\bd{\beta} + O_{kj} + \phi_{kj} ,\nonumber\\
\bd{\beta}&\sim&\mbox{N}(\bd{\mu}_{\beta}, \bd{\Sigma}_{\beta}).\nonumber
\end{eqnarray}

In common with the univariate models $(\mathbf{x}_{kj}^{\top}, O_{kj}, \phi_{kj})$ are respectively a vector of $p$ covariates, an offset, and a random effect for area $k$ and variable $j$. Gaussian priors are again assumed for the regression parameters 
$\bd{\beta}$ in common with the univariate models described above. At present only binomial and Poisson data likelihood models are allowed, that is:

\begin{itemize}
\item \textbf{Binomial - } $Y_{kj}~\sim~\mbox{Binomial}(n_{kj}, \theta_{kj})$ and $\ln(\theta_{kj}/(1-\theta_{kj}))~=~\mathbf{x}_{kj}^{\top}\bd{\beta} + O_{kj} + \phi_{kj}$. 

\item \textbf{Poisson - } $Y_{kj}~\sim~\mbox{Poisson}(\mu_{kj})$ and $\ln(\mu_{kj})~=~\mathbf{x}_{kj}^{\top}\bd{\beta} + O_{kj} + \phi_{kj}$. 
\end{itemize}

The set of random effects are denoed by $\bd{\phi}=(\bd{\phi}_1,\ldots,\bd{\phi}_{K})$, where $\bd{\phi}_{k}=(\phi_{k1},\ldots,\phi_{kJ})$. Thus, when fitting this model all data vectors (e.g. response, offset, covariates) should be ordered so that all $J$ variables for areal unit 1 come first and are followed by all $J$ variables for areal unit 2 and so on.  The random effects need to model both spatial autocorrelation and between variable correlation, and this is achieved using a multivariate conditional autorgressive (MCAR) model, for details see \cite{gelfand2003}. At present \pkg{CARBayes} can only fit the following MCAR model.\\

\code{MVS.CARleroux()}\\
The random effects $\bd{\phi}$ are modelled by:

\begin{equation}
\bd{\phi}~\sim~\mbox{N}\left(\mathbf{0},\left[\mathbf{Q}(\mathbf{W},\rho) \otimes \bd{\Sigma}^{-1}\right]^{-1}\right)
\label{MCAR}
\end{equation}


Here $\mathbf{Q}(\mathbf{W},\rho)=\rho[\mbox{diag}(\mathbf{W}\mathbf{1}) - \mathbf{W}] + (1-\rho)\mathbf{I}$ is the precision matrix for the joint distribution corresponding to the CAR prior proposed by \cite{leroux1999} and described above in the univariate models section, while $\bd{\Sigma}$ is a cross variable covariance matrix. In common with the univariate models, the correlation structure imposed by (\ref{MCAR}) is more easily seen by its full conditional form, that is:


\begin{eqnarray}
\bd{\phi}_{k}| \bd{\phi}_{-k},\mathbf{W},\bd{\Sigma}, \rho&\sim&\mbox{N}\left(\frac{\rho\sum_{i=1}^{K}w_{ki}\bd{\phi}_{i}}{\rho\sum_{i=1}^{K}w_{ki} + 1-\rho},~
\frac{\bd{\Sigma}}{\rho\sum_{i=1}^{K}w_{ki} + 1-\rho}\right),\nonumber\\
\bd{\Sigma}&\sim&\mbox{Inverse-Wishart}(df, \bd{\Omega}),\nonumber\\
\rho&\sim&\mbox{Uniform}(0,1),\nonumber
\end{eqnarray}

where $\bd{\phi}_{-k}$ denotes the vector of random effects except those relating to the $k$th areal unit. Here $df$ is the degrees of freedom for the Inverse-Wishart prior for $\bd{\Sigma}$ and a default value of $df=J+1$ is selected by the software. Similarly,  $\bd{\Omega}$ is the $J\times J$ scale matrix, with the default value being the identity matrix. In common with the univariate model \code{S.CARleroux} the spatial autocorrelation parameter $\rho$ can be fixed to any value in the unit interval using the arguments (\code{fix.rho, rho}).

\subsection{Inference}
All models in this package are fitted in a Bayesian setting using McMC simulation. A combination of Gibbs sampling (when the appropriate full conditional distributions are proportional to standard statistical distributions) and Metropolis / Metropolis-Hastings steps are used, and the majority of the latter use simple random walk proposals. The overall functions that implement the McMC algorithms are written in \proglang{R}, while the computationally intensive updating steps are written as computationally efficient  \proglang{C++} routines using the \proglang{R} package \pkg{Rcpp} (\citealp{eddelbuettel2011}). Additionally, the sparsity of the neighbourhood matrix $\mathbf{W}$ is utilised via its triplet form when updating the random effects  within the algorithms, which increases the computational efficiency of the software.



%%%%%%%%%%%%%%
%%%% Section 3
%%%%%%%%%%%%%%
\section{Loading and using the software}

\subsection{Loading the software}
\pkg{CARBayes} is an add-on package to the statistical software \proglang{R}, and is freely available to download from the Comprehensive R Archive Network (CRAN, http://cran.r-project.org/) for Windows, Linux and Apple platforms. The package requires \proglang{R} ($\geq$ 3.0.0) and depends on packages \pkg{MASS} (\citealp{mass2002}), and \pkg{Rcpp} ($\geq$ 0.11.5). Additionally, it imports functionality from the following other packages: \pkg{CARBayesdata} (\citealp{carbayesdata}), \pkg{coda} (\citealp{coda2006}), \pkg{matrixcalc} (\citealp{matrixcalc}), \pkg{MCMCpack} (\citealp{mcmcpack}), \pkg{spam} (\citealp{spam2010}), \pkg{sp} (\citealp{bivand2013b}), \pkg{spdep}, \pkg{stats}, \pkg{truncdist} (\citealp{truncdist2012}) and \pkg{utils}. Once installed it can be loaded using the command \code{library(CARBayes)}.\\

Note, that certain functionality from the packages listed in the previous paragraph are automatically loaded upon loading \pkg{CARBayes}, but only for use within the package. However, a complete spatial analysis will typically also include the creation of the neighbourhood matrix $\mathbf{W}$ from a shapefile, the production of spatial maps of the fitted values and residuals, and tests for the presence of spatial autocorrelation. To achieve these tasks the following packages should be loaded separately into \proglang{R}: \pkg{maptools}, \pkg{shapefiles}, \pkg{sp} and \pkg{spdep}.



\subsection{Using the software}
The software can fit  five main models: \code{S.CARbym()}, \code{S.CARleroux()}, \code{S.CARdissimilarity()} and \code{S.CARlocalised()} for univariate spatial data and  \code{MVS.CARleroux()} for multivariate data, and full details of the arguments required for each function are given in the helpfiles. However, the main arguments common to all the functions  that are required for a baseline analysis (for example using default priors) are as follows.

\begin{itemize}
\item \code{formula} - A formula for the covariate part of the model using the same syntax used in the \code{lm()} function. Offsets can be included here using the \code{offset()} function.

\item \code{family} - The likelihood model which must  be one of \code{`binomial'}, \code{`Gaussian'} or \code{`Poisson'}.

\item \code{trials} - This is only needed if  \code{family=`binomial'}, and is a vector the same length as the response containing the total number of trials  for each area. 

\item \code{W} - A $K \times K$ symmetric and non-negative neighbourhood matrix whose row sums must all be positive. Typically a binary specification is used, where the $kj$th element $w_{kj}$ equals one if areas $(\mathcal{S}_j, \mathcal{S}_k)$ are spatially close (e.g. share a common border) and is zero otherwise. This matrix can be created from a shapefile and data frame using the \code{combine.data.shapefile()} function in the \pkg{CARBayes} package, which is illustrated in the next section.

\item \code{burnin} - The number of McMC samples to discard as the burn-in period.
    
\item \code{n.sample} - The number of McMC samples to generate.
\end{itemize}

When a model has run \pkg{CARBayes} has its own \code{print()} function that prints a summary of the fitted model to the screen, including both parameter summaries and convergence diagnostics for the McMC run. Additionally, it has functions \code{summarise.samples()} and \code{summarise.lincomb()} to summarise the results, and both functions are illustrated in the examples that follow. The software updates the user on its progress to the \proglang{R} console, which allows the user to monitor the function's progress. However, using the \code{verbose=FALSE} option will disable this feature. Once run, each model in this package returns a list object with the following components. 


\begin{itemize}
\item \code{summary.results} - A summary table of selected parameters that is presented when using the \code{print()} function. The table includes the posterior median (\code{Median}) and 95$\%$ credible interval (\code{2.5\%, 97.5\%}), the number of samples generated (\code{n.sample}), the acceptance rate for the Markov chain (\code{\% accept}), the effective number of independent samples using the  \code{effectiveSize()} function from the \pkg{coda} package (\code{n.effective}), and the convergence  diagnostic proposed by \cite{geweke1992} and implemented in the \pkg{coda} package (\code{Geweke.diag}). This diagnostic takes the form of a Z-score, so that convergence is suggested by the statistic being within the range (-1.96, 1.96).

\item \code{samples} - A list containing the McMC samples from the model, where each element in the list is a matrix. The names of these matrix objects correspond to the parameters defined in Section 2 of this paper, and each column contains the set of samples for a single parameter. This list includes samples from the posterior distribution of the fitted values for each data point (\code{fitted}).

\item \code{fitted.values} - A vector of fitted values based on posterior medians for each area.

\item \code{residuals} - A vector of raw residuals based on posterior medians for each area.

\item \code{modelfit} - Model fit criteria including the Deviance Information Criterion (DIC, \cite{spiegelhalter2002}) and its corresponding estimated effective number of parameters (p.d), the Watanabe-Akaike Information Criterion (WAIC, \cite{watanabe2010}) and its corresponding estimated number of effective parameters (p.w), and the Log Marginal Predictive Likelihood (LMPL, \cite{congdon2005}). The best fitting model is one that minimises the DIC and WAIC but maximises the LMPL.

\item \code{accept }{The acceptance probabilities for the parameters.}

\item \code{localised.structure} - This element is \code{NULL} except for the models \code{S.CARdissimilarity()} and \code{S.CARlocalised()}. For \code{S.CARdissimilarity} it is a list containing two matrices. \code{W.posterior} contains posterior medians for each element $w_{kj}$ of the $K \times K$ neighbourhood matrix $\mathbf{W}$, while \code{W.border.prob} contains posterior probabilities that each $w_{kj}$ element equals zero, which corresponds to the posterior probability of a boundary in the random effects surface. In both cases elements which correspond to two non-neighbouring areas have NA values. For \code{S.CARlocalised()} this element is a vector of length $K$, and gives the posterior median class ($Z_{k}$ value) that each data point is assigned to.
        
\item \code{formula} - The formula (as a text string) for the covariate and offset part of the model.

\item \code{model}- A text string describing the model that has been fitted.

\item \code{X} - The design matrix of covariates inherited from the \code{formula} argument.
\end{itemize}



The remainder of this vignette illustrates the \pkg{CARBayes} software via 3 worked examples.




\section{Example 1 - Scottish lip cancer data}
The first example is the famous Scottish lip cancer data set, which is included purely to illustrate how to combine a data frame and shapefile together into a \code{SpatialPolygonsDataFrame} object. The creation of this object allows spatial maps to be produced of variables of interest, as well as allowing the neighbourhood matrix $\mathbf{W}$ to be created for use in the models implemented in \pkg{CARBayes}. The Scottish lip cancer data are contained in the \pkg{CARBayesdata} package, which can be loaded using the following command:


<<>>=
library(CARBayesdata)
library(shapefiles)
library(sp)
data(lipdata)
data(lipdbf)
data(lipshp)
@

To create a \code{SpatialPolygonsDataFrame} object you essentially need two types of data. The first is a \code{data.frame}, containing the data you wish to model. If this is a  comma separated variable (csv) file then it can be read into \proglang{R} using the command \code{read.csv()}. The second data type is a shapefile, which comprises two components: \code{shapefile.shp} containing the polygons, and \code{shapefile.dbf} containing a unique identifier linking each row in the \code{data.frame} to a polygon in the \code{shapefile.shp} file. The shapefiles can be read in to \proglang{R} using the \code{read.shp()} and \code{read.dbf()} functions. In the above example the \code{data(lipdata)} loads the \code{data.frame} object, \code{data(lipdbf)} loads the \emph{.dbf} component of the shapefile while \code{data(lipshp)} loads the \emph{.shp} component of the shapefile. These three data sets can be combined together to create a \code{SpatialPolygonsDataFrame} object using the \code{combine.data.shapefile()} function from \pkg{CARBayes}.

<<>>=
library(CARBayes)
lipdbf$dbf <- lipdbf$dbf[ ,c(2,1)]
data.combined <- combine.data.shapefile(data=lipdata, shp=lipshp, dbf=lipdbf)
@

For this function to work the rownames of the dataframe (\code{lipdata}) must be contained in the first column of the \emph{.dbf} (\code{lipdbf\$dbf}) object, which is the reason for re-ordering the columns in the second line of the above code. The \code{data.combined} object is a \code{SpatialPolygonsDataFrame} object, and the data it contains can be mapped using the \code{spplot()} function. A binary neighbourhood matrix $\mathbf{W}$ based on sharing a common border can be constructed from this object using the code:


<<>>=
library(spdep)
W.nb <- poly2nb(data.combined, row.names = rownames(lipdata))
W.mat <- nb2mat(W.nb, style="B")
@

This example was included purely to illustrate this process, and in the next two examples full spatial analyses using the \pkg{CARBayes} software will be illustrated.





\section{Example 2 - property prices in Greater Glasgow}
The utility of the \pkg{CARBayes} software is illustrated by modelling the spatial pattern in average property prices across Greater Glasgow, Scotland, in 2008. This is an ecological regression analysis, whose aim is to  identify the factors that affect property prices and quantify their effects.


\subsection{Data and exploratory analysis}
The data come from the Scottish Statistics database (\emph{http://statistics.gov.scot}), but are also included in the \pkg{CARBayesdata} \proglang{R} package. The study region is the Greater Glasgow and Clyde health board (GGHB), which is split into 271 intermediate geographies (IG). These IGs are small areas that have a median area of 124 hectares and a median population of 4,239. These data can be loaded into \proglang{R} using the code shown below:


<<>>=
library(CARBayesdata)
library(sp)
data(GGHB.IG)
data(pricedata)
@

The \code{GGHB.IG} object is a \code{SpatialPolygonsDataFrame} object containing the spatial information for the GGHB, which is used to plot the data, construct the neighbourhood matrix $\mathbf{W}$, and conduct a test for spatial autocorrelation. The \code{pricedata} object is a  \code{data.frame} containing the property price data. However, the latter only has data for 270 of the 271 IGs in GGHB, as one area had outlying values and was hence removed. Thus \code{GGHB.IG} needs to have the appropriate row removed and then be combined with the \code{pricedata} \code{data.frame}. This is achieved using the code below:

<<>>=
missing.IG <- setdiff(rownames(GGHB.IG@data), pricedata$IG)
missing.IG.row <- which(missing.IG==rownames(GGHB.IG@data))
propertydata.spatial <- GGHB.IG[-missing.IG.row, ]
propertydata.spatial@data <- data.frame(propertydata.spatial@data, pricedata)
@

The first two lines identify the missing row and its row number, while the third line subsets the \code{GGHB.IG} object by removing that row to create the \code{propertydata.spatial} object. Finally, the last line adds the data contained in \code{pricedata} to the \code{propertydata.spatial} object.\\

The variables in \code{pricedata} are summarised in Table \ref{table_scotland data}, which displays the percentiles of their distribution (with the exception of the categorical variable \code{type}). The response variable in this study is the median price (in thousands, \code{price}) of all properties sold in 2008 in each IG. The table shows large variation in this variable, with average prices ranging between $\pounds 50,000$ and $\pounds 372,800$ across the study region. The first covariate in this study is the crime rate (\code{crime}) in each IG, because areas with higher crime rates are likely to be less desirable to live in. Crime rate is measured as the total number of recorded crimes in each IG per 10,000 people that live there, and the values range between 85 and 1994. Other covariates included in this study are the median number of rooms in a property (\code{rooms}), the number of properties that sold in a year (\code{sales}), and the average time taken to drive to the nearest shopping centre (\code{driveshop}). The latter is a proxy measure of access to services which may affect property prices. Finally, a categorical variable measuring the  most prevalent property type in each area is available (\code{type}), with levels; `flat' (68$\%$ of areas), `terraced' (7$\%$), `semi-detached' (13$\%$) and `detached' (12$\%$).\\ 


\begin{table}
\centering\begin{tabular}{lrrrrr}
  \hline
\textbf{Variable}&\multicolumn{5}{c}{\textbf{Percentiles}}\\
&0$\%$&25$\%$&50$\%$&75$\%$&100$\%$\\\hline
Property price (in thousands)&50.0& 95.0 &121.8 &159.2 &372.8\\
Crime rate (per 10,000)&85.0&  303.2&  517.0&  728.0& 1994.0\\
Number of rooms (median)&3 &3 &4 &4 &6 \\
Property sales &4 & 46 & 58 & 85 & 266\\
Drive time to a shop (minutes)&0.3 & 0.8 & 1.2 & 1.9 & 8.5\\\hline
\end{tabular}
\caption{\label{table_scotland data} Summary of the distribution of the data.}
\end{table}



A spatial map of the \code{price} response variable can be plotted using the functionality of the \pkg{sp} package, using the following \proglang{R} code.



<<>>=
northarrow <- list("SpatialPolygonsRescale", layout.north.arrow(), 
    offset = c(220000,647000), scale = 4000)
scalebar <- list("SpatialPolygonsRescale", layout.scale.bar(), 
    offset = c(225000,647000), scale = 10000, fill=c("transparent","black"))
text1 <- list("sp.text", c(225000,649000), "0")
text2 <- list("sp.text", c(230000,649000), "5000 m")
breakpoints <- seq(min(pricedata$price)-1, max(pricedata$price)+1, length.out=8)
spplot(propertydata.spatial, c("price"), sp.layout=list(northarrow, scalebar, 
    text1, text2), scales=list(draw = TRUE), at=breakpoints,
    col.regions=terrain.colors(n=length(breakpoints)-1), col="transparent")
@


\begin{figure}
\centering 
\scalebox{1}{\includegraphics{pricemap.png}}
\caption{Map showing the average property price in each IG (in thousands).\label{pricemap}}
\end{figure} 


The map is shown in Figure \ref{pricemap} and is achieved by the \code{spplot()} function, with the preceding lines adding a north arrow, a scale bar and accompanying text. The figure suggests that Glasgow has a number of property sub-markets, whose prices are not related to those in neighbouring areas. An example of this is the two groups of higher priced regions north of the river Clyde (the thin white line running south east), which are the highly sought after Westerton / Bearsden (northerly cluster) and Dowanhill / Hyndland (central cluster) districts.\\ 


\subsection{Non-spatial modelling}
The natural log of the median property price variable \code{price} is treated as the response and assumed to be Gaussian, and an initial covariate only model is built in a frequentist framework using linear models. Initial plots of the data using the \code{pairs()} command suggest that the natural log of  drive time to a shopping centre (\code{driveshop}) is linearly related to the response, and that crime rate (\code{crime}) has a non-linear relationship to the response. The natural log transformation of \code{price} and \code{driveshop} are created using the following commands.

<<>>=
propertydata.spatial@data$logprice <- log(propertydata.spatial@data$price)
propertydata.spatial@data$logdriveshop <- log(propertydata.spatial@data$driveshop)
@
    
A model with all the covariates is fitted to the data, where the crime rate variable is modelled as non-linear using a natural cubic spline with 3 degrees of freedom. This is achieved using the following \proglang{R} code:
    
<<>>=
library(splines)
form <- logprice~ns(crime,3)+rooms+sales+factor(type) + logdriveshop
model <- lm(formula=form, data=propertydata.spatial@data)
@

From fitting this model all of the numeric covariates are significantly related to the response at the 5$\%$ level, suggesting they all play an important role in explaining the spatial pattern in median property price. The predominant property type variable also appears to be important, with areas where the level is `detached' (the baseline level) having significantly higher property prices than the other three levels.\\ 

To quantify the presence of spatial autocorrelation in the residuals from this model we can compute Moran's I statistic (\citealp{moran1950}) and conduct a permutation test. The permutation test has the null hypothesis of no spatial autocorrelation and an alternative hypothesis of spatial autocorrelation (either positive or negative), and is conducted using the \code{moran.mc()} function from the \pkg{spdep} package. The test can be implemented using the code below. Lines 2 and 3 turn the \code{SpatialPolygonsDataFrame} object \code{propertydata.spatial} into an \code{nb} and then a \code{listw} object, which is required by the \code{moran.mc()} function.

<<>>=
library(spdep)
W.nb <- poly2nb(propertydata.spatial, row.names = rownames(propertydata.spatial@data))
W.list <- nb2listw(W.nb, style="B")
resid.model <- residuals(model)
moran.mc(x=resid.model, listw=W.list, nsim=1000)
@
    
    
The Moran's I statistic equals 0.2733 with a corresponding p-value of much less than 0.05, which suggests that the residuals contain substantial positive spatial autocorrelation.

\subsection[Spatial modelling with CARBayes]{Spatial modelling with \pkg{CARBayes}}
The residual spatial autocorrelation can be accounted for by adding a set of random effects to the model, using the functions outlined in the previous section. We illustrate this by applying model (\ref{equation likelihood}) and  (\ref{equation leroux}) to the data. The code to implement this model in \pkg{CARBayes} is shown below, where the first line creates the binary neighbourhood matrix \code{W} from the \code{W.nb} object.


\begin{CodeInput}
R> library(CARBayes)
R> W <- nb2mat(W.nb, style="B")
R> model.spatial <- S.CARleroux(formula=form, data=propertydata.spatial@data, 
+   family="gaussian", W=W, burnin=20000, n.sample=120000, thin=10)
\end{CodeInput}

In the above code the covariate and offset component defined by \code{formula} is the same as for the simple linear model fitted earlier, and the neighbourhood matrix \code{W} is binary and defined by whether or not two areas share a common border. Inference for this model is based on 10,000 McMC samples, which were obtained following a burn-in period of 20,000 and thinning the remaining 100,000 samples by 10 to reduced their autocorrelation. A summary of the fitted model can be obtained using the following code.

\begin{CodeInput}
R> print(model.spatial)
\end{CodeInput}

\begin{CodeOutput}
#################
#### Model fitted
#################
Likelihood model - Gaussian (identity link function) 
Random effects model - Leroux CAR
Regression equation - logprice ~ ns(crime, 3) + rooms + sales + factor(type) + 
logdriveshop
Number of missing observations - 

############
#### Results
############
Posterior quantities and DIC

                     Median    2.5%   97.5% n.sample % accept n.effective Geweke.diag
(Intercept)          4.2419  3.9630  4.5179    10000    100.0     10000.0         0.8
ns(crime, 3)1       -0.2459 -0.3967 -0.0968    10000    100.0      7406.5         0.7
ns(crime, 3)2       -0.4010 -0.7049 -0.1046    10000    100.0      8004.3         0.1
ns(crime, 3)3       -0.2007 -0.4073  0.0102    10000    100.0     10000.0         1.2
rooms                0.2198  0.1693  0.2711    10000    100.0     10000.0        -0.5
sales                0.0022  0.0016  0.0029    10000    100.0     10000.0        -1.3
factor(type)flat    -0.2488 -0.3680 -0.1299    10000    100.0     10000.0        -0.2
factor(type)semi    -0.1622 -0.2632 -0.0611    10000    100.0      9601.3         0.9
factor(type)terrace -0.2943 -0.4222 -0.1684    10000    100.0      9635.0        -2.0
logdriveshop        -0.0050 -0.0611  0.0500    10000    100.0      7013.5         1.0
nu2                  0.0245  0.0143  0.0340    10000    100.0      2244.5         0.7
tau2                 0.0433  0.0188  0.0826    10000    100.0      1893.7        -1.0
rho                  0.9325  0.7402  0.9912    10000     45.4      3662.0         0.9

DIC =  -134.9729       p.d =  98.37034       LMPL =  151.1882 
\end{CodeOutput}


The Summary is presented in two parts, the first of which describes the model that has been fit. The second summarises the results, and includes the posterior median (\code{Median}) and 95\% credible intervals (\code{2.5\%, 97.5\%}) for selected parameters (not the random effects), the convergence diagnostic proposed by \cite{geweke1992} (\code{Geweke.diag}) as a Z-score and the effective number of independent samples (\code{n.effective}). Also displayed are the actual number of samples kept from the McMC run  (\code{n.sample}) as well as the acceptance rate for each parameter (\code{\% accept}). Note, parameters that have an acceptance rate of 100\% have been Gibbs sampled due to their full conditional distributions being a standard distribution. Finally, the DIC and LMPL overall model fit criteria are displayed.\\

In addition to producing the summary above, fitting the model returns a list object with the following components:

\begin{CodeInput}
R> summary(model.spatial)
\end{CodeInput}

\begin{CodeOutput}
                    Length Class   Mode     
summary.results       91   -none-  numeric  
samples                7   -none-  list     
fitted.values        270   -none-  numeric  
residuals            270   -none-  numeric  
modelfit               5   -none-  numeric  
accept                 5   -none-  numeric  
localised.structure    0   -none-  NULL     
formula                3   formula call     
model                  2   -none-  character
X                   2700   -none-  numeric  
\end{CodeOutput}



The first element of this list is a summary table of results, which is what is printed by the \code{print()} function. The next element is a list containing matrices of the thinned and post burn-in McMC samples for each set of parameters.  For example, \code{model.spatial$samples$beta} is a matrix containing the McMC samples for all the regression parameters. The next two elements in the list \code{fitted.values}  and \code{residuals} are vectors of fitted values and residuals from the model, while \code{modelfit} gives a selection of model fit criteria. These criteria include the Deviance Information Criterion (DIC), the log Marginal Predictive Likelihood (LMPL) and the Watanabe-Akaike Information Criterion (WAIC). For further details about Bayesian modelling and model fit criteria see \cite{gelman2003}. The item \code{accept} contains the acceptance rates for the model, while \code{localised.structure} is \code{NULL} for this model and is used for compatabiltiy with the other functions in the package. Finally, the \code{formula} and \code{model} are text strings describing the formula used and the model fit, while \code{X} gives the design matrix corresponding to the \code{formula} object.



\subsection{Inference}
The summary table above gives posterior medians and 95$\%$ credible intervals for a selection of model parameters, but these can be recreated (or similar summarise created for other parameters) using the function

\begin{CodeInput}
R> summarise.samples(model.spatial$samples$beta, quantiles=c(0.5, 0.025, 0.975))
\end{CodeInput}

\begin{CodeOutput}
$quantiles
               0.5        0.025       0.975
 [1,]  4.241896323  3.962957346  4.51787943
 [2,] -0.245877968 -0.396698998 -0.09679782
 [3,] -0.401030811 -0.704894560 -0.10461861
 [4,] -0.200704600 -0.407327946  0.01018534
 [5,]  0.219841199  0.169342794  0.27108763
 [6,]  0.002237906  0.001611683  0.00287336
 [7,] -0.248767056 -0.368024984 -0.12990882
 [8,] -0.162194165 -0.263154252 -0.06109645
 [9,] -0.294254065 -0.422182354 -0.16841391
[10,] -0.005002887 -0.061142928  0.04995092

$exceedences
NULL
\end{CodeOutput}


which here has summarised posterior medians and 95$\%$ credible intervals for the covariate effects. However, for the crime variable its relationship is non-linear and summarised by the results for all 3 basis functions \code{ns(crime, 3)1, ns(crime, 3)2, ns(crime, 3)3}. Therefore to summarise the entire linear relationship we can use the \code{summarise.lincomb()} function, which allows us to compute the posterior distribution and quantiles of a linear combination of the covariates. This can be achieved and then plotted using the code:


\begin{CodeInput}
R> crime.effect <- summarise.lincomb(model=model.spatial, columns=c(2,3,4), 
+    quantiles=c(0.5, 0.025, 0.975), distribution=FALSE)
R> plot(propertydata.spatial@data$crime, crime.effect$quantiles[ ,1], pch=19, 
+    ylim=c(-0.55,0.05), xlab="Number of crimes", ylab="Effect of crime")
R> points(propertydata.spatial@data$crime, crime.effect$quantiles[ ,2], pch=19, col="red")
R> points(propertydata.spatial@data$crime, crime.effect$quantiles[ ,3], pch=19, col="red")
\end{CodeInput}


\begin{figure}
\centering 
\scalebox{1}{\includegraphics{crimeplot.png}}
\caption{Plot showing the estimated non-linear relationship between crime rate and log-price.\label{crimeplot}}
\end{figure} 

The first line creates the crime effect, while the remaining lines plot the graph which is shown in Figure \ref{crimeplot} and shows a decreasing relationship as crime rate increases as expected.



\section{Example 3 - identifying high-risk disease clusters}
The third example illustrates the utility of the localised spatial autocorrelation model proposed by \cite{lee2012}, which can identify boundaries that represent step changes in the (random effects) response surface between geographically adjacent areal units. The aim of this analysis is to identify boundaries in the risk surface of respiratory disease in Greater Glasgow, Scotland in 2010, so that the spatial extent of high-risk clusters can be identified. The identification of boundaries in spatial data is affectionately known as \emph{Wombling}, after the seminal paper by \cite{womble1951}. 


\subsection{Data and exploratory analysis}
The data again relate to the Greater Glasgow and Clyde health board, and are also freely available to download from  Scottish Statistics (\emph{http://statistics.gov.scot}). However, the river Clyde partitions the study region into a northern and a southern sub-region, and no areal units on opposite banks of the river border each other. This means that boundaries could not be identified across the river, and therefore here we only consider those areal units that are on the northern side of the study region. This leaves 134 areal units in the new smaller study region, and the data on respiratory disease risk are included with the \pkg{CARBayesdata} package and can be loaded with the command:


<<>>=
library(CARBayesdata)
library(sp)
data(GGHB.IG)
data(respiratorydata)
@

The \code{GGHB.IG} \code{spatialPolygonsDataFrame} object can then be subset to just include IGs in the \code{respiratorydata} data.frame using the following code.

<<>>=
missing.IG <- setdiff(rownames(GGHB.IG@data), respiratorydata$IG)
missing.IG.row <- rep(NA, length(missing.IG))
    for(i in 1:length(missing.IG))
    {
    missing.IG.row[i] <- which(missing.IG[i] == rownames(GGHB.IG@data))     
    }
respiratorydata.spatial <- GGHB.IG[-missing.IG.row, ]
respiratorydata.spatial@data <- data.frame(respiratorydata.spatial@data, 
    respiratorydata)
@

This code is similar to that in the previous example, except this time more than one area needs to be removed, hence the \code{for} loop. The the first 6 rows of the data can be viewed using the code below.

<<>>=
head(respiratorydata.spatial@data)
@

In addition to the unique identifier IG codes (\code{IG}) and the geographical coordinates of each area's centroids (\code{easting, northing}), the data contain 4 variables.  \code{observed} is the  numbers of hospital admissions in 2010 in each IG due to respiratory disease (International Classification of Disease tenth revision codes J00-J99). However, these observed numbers will depend on the size and demographic structure of the populations living in each IG, and these factors need to be adjusted for before estimating disease risk. This is typically achieved  by computing the expected numbers of hospital admissions in each IG based on this demographic information, using indirect standardisation. These expected numbers are stored in the \code{expected} column, and the simplest measure of disease risk is the Standardised Morbidity Ratio (SMR), which is the ratio of the observed to the expected numbers of hospital admissions and stored in the \code{SMR} column. Finally, the percentage of people defined to be income deprived (in receipt of means tested benefits) is stored in the \code{incomedep} column. A map showing the SMR is displayed in the Figure \ref{smrmap}, and is created using the code below that has a similar structure to that used in the previous analysis



<<>>=
northarrow <- list("SpatialPolygonsRescale", layout.north.arrow(), offset = 
    c(220000,647000), scale = 4000)
scalebar <- list("SpatialPolygonsRescale", layout.scale.bar(), offset = 
    c(225000,647000), scale = 10000, fill=c("transparent","black"))
text1 <- list("sp.text", c(225000,649000), "0")
text2 <- list("sp.text", c(230000,649000), "5000 m")
breakpoints <- seq(min(respiratorydata.spatial@data$SMR)-0.05, 
    max(respiratorydata.spatial@data$SMR)+0.05, length.out=8)
spplot(respiratorydata.spatial, c("SMR"), sp.layout=list(northarrow, scalebar, 
    text1, text2), scales=list(draw = TRUE), at=breakpoints, 
    col.regions=terrain.colors(n=length(breakpoints)-1), col="transparent")
@


\begin{figure}
\centering 
\scalebox{1}{\includegraphics{smrmap.png}}
\caption{Map displaying the SMR for each area.\label{smrmap}}
\end{figure} 


Values of the SMR in the figure above one relate to areas exhibiting above average risks, while values below one correspond to below average risks. For example, an SMR of 1.2 corresponds to a 20\% increased risk relative to the expected numbers of respiratory disease cases. The figure shows evidence of localised spatial structure in these disease data, with numerous different locations where high and low risk areas border each other. This in turn suggests that boundaries are likely to be present in these data, and their identification is the goal of this analysis. The method proposed by \cite{lee2012} identifies these boundaries using dissimilarity metrics, which are non-negative measures of the dissimilarity between all pairs of adjacent areas. In this example we use the absolute difference in the percentage of people in each IG who are defined to be income deprived (\code{incomedep}), because it is well known that socio-economic deprivation plays a large role in determining people's health. However, before fitting the model the spatial neighbourhood matrix $\mathbf{W}$ based on sharing a common border is computed using the following code.



<<>>=
W.nb <- poly2nb(respiratorydata.spatial, row.names = 
    rownames(respiratorydata.spatial@data))
W <- nb2mat(W.nb, style="B")
@



\subsection[Spatial modelling with CARBayes]{Spatial modelling with \pkg{CARBayes}}
Let the  observed and expected numbers of hospital admissions be denoted by $\mathbf{Y}=(Y_{1},\ldots,Y_{K})$ and $\mathbf{E}=(E_{1},\ldots,E_{K})$ respectively. Then as the observed numbers of hospital admissions are counts, a Poisson likelihood model given by $Y_{k}\sim\mbox{Poisson}(E_{{k}}R_{k})$ is appropriate, where $R_{k}$ represents disease risk in areal unit $\mathcal{S}_{k}$. A log-linear model is specified for $R_{k}$, that is, $\ln(R_{k})=\beta_{0}+\phi_{k}$, and for a general review of disease mapping see \cite{wakefield2007}. We note that in fitting this model in \pkg{CARBayes}, the offset is specified on the linear predictor scale rather than the expected value scale, so in this analysis the offset is $\log(\mathbf{E})$ rather than $\mathbf{E}$. The dissimilarity metric used here is the absolute difference in the level of income deprivation between each pair of areal units, which can be created from the  vector of area level income deprivation scores using the following code.

<<>>=
income <- respiratorydata.spatial@data$incomedep
Z.incomedep <- as.matrix(dist(cbind(income, income), method="maximum", diag=TRUE, 
    upper=TRUE))
@
    
The first line pulls out the income deprivation covariate while the second line computes the $K\times K$ matrix of absolute differences in income deprivation levels between each pair of areal units. The function to implement the localised CAR model is called \code{S.CARdissimilarity()}, and it takes the same arguments as the global CAR models except that it additionally requires the dissimilarity metrics. These are required in the form of a list of $K\times K$ matrices, and for this example we only have a single dissimilarity metric. The model is run using the following code.




\begin{CodeInput}
R> formula <- observed ~ offset(log(expected))
R> model.dissimilarity <- S.CARdissimilarity(formula=formula, 
+    data=respiratorydata.spatial@data, family="poisson", W=W, 
+    Z=list(Z.incomedep=Z.incomedep), burnin=20000, n.sample=60000, thin=10)
R> print(model.dissimilarity)
\end{CodeInput}


\begin{CodeOutput}
#################
#### Model fitted
#################
Likelihood model - Poisson (log link function) 
Random effects model - Localised CAR 
Dissimilarity metrics -  Z.incomedep 
Regression equation - observed ~ offset(log(expected))
Number of missing observations - 

############
#### Results
############
Posterior quantities and DIC

             Median    2.5%   97.5% n.sample % accept n.effective Geweke.diag alpha.min
(Intercept) -0.2196 -0.2421 -0.1985     4000     31.9      3652.8        -0.7        NA
tau2         0.1207  0.0745  0.2380     4000    100.0       259.6        -2.1        NA
Z.incomedep  0.0535  0.0378  0.0577     4000     45.9       197.8         2.6    0.0159

DIC =  1072.698       p.d =  101.9844       LMPL =  -451.082 

The number of stepchanges identified in the random effect surface
     no stepchange stepchange
[1,]           215        145
\end{CodeOutput}

The first line of the above code specifies the formula with an offset (the natural log of the expected numbers of cases) but no covariates, the latter being required so that boundaries identified in the random effects surface can also be interpreted as boundaries in the risk surface (that is $\mathbf{R}=(R_{1},\ldots,R_{n})$). The above model bases inference on 4,000 post burn-in and thinned McMC samples. When the model has been fit the \code{print()} function produces the summary output above which is similar to that produced for the property price data in the previous example. The main difference between this and the corresponding output from the property price analysis is the addition of a column in the parameter summary table headed \code{alpha.min}. This column only applies to the dissimilarity metrics, which is why it is \code{NA} for the remaining parameters. The value of \code{alpha.min} is the threshold value for the regression parameter $\alpha$, below which the dissimilarity metric has had no effect in identifying boundaries in the response (random effects) surface. A brief description is given in Section 2.1, while full details are given in \cite{lee2012}. For these data the posterior median and 95$\%$ credible interval lie completely above this threshold, suggesting that the income deprivation dissimilarity metric has identified a number of boundaries.\\

The number and locations of these boundaries are summarised in the element of the output list called \code{model.dissimilarity$localised.structure$W.posterior}, which is a $K\times K$ symmetric matrix containing the posterior median for the set $\{w_{kj}|k\sim j\}$. Values equal to zero represent a boundary, values equal to one correspond to no boundary, while \code{NA} values correspond to non-adjacent areas. The locations of these boundaries can be overlaid on a map of the estimated disease risk (that is the posterior median of $\mathbf{R}$) using the following code. 

\begin{CodeInput}
R> border.locations <- model.dissimilarity$localised.structure$W.posterior
R> respiratorydata.spatial@data$risk <- model.dissimilarity$fitted.values / 
+   respiratorydata.spatial@data$expected
R> boundary.final <- highlight.borders(border.locations=border.locations, 
+   spdata=respiratorydata.spatial)
R> spplot(respiratorydata.spatial, c("risk"), sp.layout=list(northarrow, scalebar, 
+   text1, text2, boundary.final), scales=list(draw = TRUE), at=breakpoints, 
+   col.regions=terrain.colors(n=length(breakpoints)-1), col="transparent")
\end{CodeInput}

\begin{figure}
\begin{center}
\scalebox{1.2}{\includegraphics{boundarymap.png}}
\end{center}
\vspace{-1cm}\caption{Map displaying estimated risk and locations of the boundaries for the northern part of Greater Glasgow.}\label{figure_boundary}
\end{figure}


The first line saves the matrix of border locations, while the second adds the estimated risk values to the \code{respiratorydata.spatial} object. The next line identifies the boundary points (using the \pkg{CARBayes} function \code{highlight.borders()}) and formats them to enable plotting. The remaining line plots the map and is similar to that used to produce the earlier spatial maps.  
The result of these commands are displayed in Figure \ref{figure_boundary}, which shows the fitted risk surface and the locations of the boundaries (denoted by blue crosses).  The model has identified 145 boundaries in the risk surface. The majority of these visually seem to correspond to sizeable changes in the risk surface, suggesting that the model has the power to distinguish between boundaries and non-boundaries. The notable boundaries are the demarkation between the low risk (shaded green) city centre / west end of Glasgow in the middle of the region and the deprived neighbouring areas on both sides, which include Easterhouse / Parkhead in the east and Knightswood / Drumchapel in the west. The other interesting feature of this map is that the boundaries are not closed, suggesting that the spatial pattern in risk is more complex than being partitioned into groups of non-overlapping areas of similar risk.


\section{Discussion}
This vignette has illustrated the \proglang{R} package \pkg{CARBayes}, which can fit a number of commonly used univariate and multivariate  conditional autoregressive models to spatial areal unit data, as well as the localised spatial smoothing models proposed by \cite{lee2012} and \cite{lee2015}. The response data can be binomial, Gaussian or Poisson, with the canonical link functions logit, identity and natural log respectively. The availability of areal unit data has grown dramatically in recent times, due to the launch of freely available online databases such as in the UK (see \emph{http://www.neighbourhood.statistics.gov.uk} and \emph{http://statistics.gov.scot}), and Surveillance Epidemiology and End Results (SEER, \emph{http://seer.cancer.gov/}) in the USA, the latter being a cancer registry. This increased availability of spatial data has fuelled a growth of modelling in this area, leading to the need for user friendly software such as \pkg{CARBayes} for use by both statisticians and non-statisticians alike.\\

A number of other software packages can also  fit conditional autoregressive models to spatial data, including \proglang{BUGS}, \proglang{BayesX} and \proglang{R} packages \pkg{hSDM}, \pkg{INLA}, \pkg{spatcounts} and \pkg{spdep}. However, these software packages either can only fit a limited selection of CAR models, or require a degree of programming which may be beyond some users of spatial data. Thus a gap in the market exists for user friendly software that can fit a wide class of CAR models, which was the motivation behind the \pkg{CARBayes} software. The user friendly features of \pkg{CARBayes} have been illustrated by the three worked examples presented in Sections 4 to 6, which include: (i) models can be implemented using a single function call; (ii) the spatial information required by the models is straightforward to create from a shapefile; (iii) only a small number of arguments are required to run a \emph{default} analysis; and (iv) the software reports on the progress of model fitting, and produces a summary table of the results when it has finished. Finally, this software now has a sister spatio-temporal modelling package called \pkg{CARBayesST}, which can fit a range of spatio-temporal areal unit models based on CAR priors. These models include similar models to those proposed by \cite{bernardinelli1995} and \cite{knorrheld2000b}, as well as other alternatives.

\section*{Acknowledgements}
The data and shapefiles used in sections 5 and 6 of this vignette were provided by the Scottish Government.


\bibliography{jss996}
\end{document}


