# Structure induced by a multiple membership transformation on the Conditional
Autoregressive model

Marco Gramatica  School of Mathematical Sciences, Queen Mary University of
London, London, UK  Silvia Liverani  School of Mathematical Sciences, Queen
Mary University of London, London, UK  The Alan Turing Institute, The British
Library, London, UK  Peter Congdon  School of Geography, Queen Mary University
of London, London, UK

###### Abstract

The objective of disease mapping is to model data aggregated at the areal
level. In some contexts, however, (e.g. residential histories, general
practitioner catchment areas) when data is arising from a variety of sources,
not necessarily at the same spatial scale, it is possible to specify spatial
random effects, or covariate effects, at the areal level, by using a multiple
membership principle (MM) [42, 23]. A weighted average of conditional
autoregressive (CAR) spatial random effects embeds spatial information for a
spatially-misaligned outcome and estimate relative risk for both frameworks
(areas and memberships). In this paper we investigate the theoretical
underpinnings of these application of the multiple membership principle to the
CAR prior, in particular with regard to parameterisation, properness and
identifiability. We carry out simulations involving different numbers of
memberships as compared to number of areas and assess impact of this on
estimating parameters of interest. Both analytical and simulation study
results show under which conditions parameters of interest are identifiable,
so that we can offer actionable recommendations to researchers. Finally, we
present the results of an application of the multiple membership model to
diabetes prevalence data in South London, together with strategic implications
for public health considerations.

##  1 Introduction

The scope for geo-referenced data analysis is vast, though particular
methodological paradigms are more common. In the context of disease mapping,
the information available is typically aggregated at areal level, thus
constraining the range of viable modelling strategies. Bayesian hierarchical
models have long been a widely popular choice for areal disease mapping [32,
12, 5]. Our interest focuses on the conditional autoregressive prior (CAR) [2]
which has been employed extensively in the disease mapping literature to
account for spatial autocorrelation which is a common occurrence in areal data
[33].

The most common framework in this context arises from the observations of an
outcome of interest $\bm{y}=(y_{1},...,y_{n})$ and some explanatory
covariates, both observed on the same set of n areas that partition a certain
region of interest, modelled with a distribution that includes spatial random
effects such as the CAR prior. In the present work we will investigate a
different type of setup, where covariates are observed at the areal level but
outcomes refer to a different spatial framework.

An instance of this non-standard setup can be found in [23], where they
investigate general practitioner (GP) level data of diabetes prevalence
jointly modelled with areal level mortality data for the same condition. In
this instance covariates of interest are available by area (i.e. census
tracts) but not for practices. However, the patients of each practice do not
necessarily reside in a unique set of areas. This gives rise to two
misalignment problems: first, the population of each practice can be
attributed to several non-contiguous areas and second, residents of a certain
area can be registered at different practices. The former type of misalignment
is more related to the (non) spatial structure of the data, given that spatial
modelling strategies usually require continuity. On the other hand, the second
misalignment problem concerns the imperfect nesting of areas in GP practices.

Another example of the same ilk can be found in [42], where the impact of
residential history of people affected by a rare type of cancer caused by
exposure to asbestos (mesothelioma), is investigated. For each patient the
probability of contracting mesothelioma is modelled using age at death as a
covariate together with a spatial random effect. Clearly, given that each
individual might have resided in several different locations throughout the
course of their life, their residential history has to be taken into account
to correctly estimate the spatial risk. Hence, the random effect is specified
in terms of a weighted average of the random effects of the areas where the
individual has resided. The weights correspond to the share of time spent by
the individual living in a certain area, hence the multiple membership of each
area to different individuals’ residential history.

Both of the examples above can be framed as applying the multiple membership
principle (MM) [7] transformation to an areal framework. Hence, we will denote
this new prior as CAR-MM. Given that this model is a weighted average of areal
random effects, it can be seen as a linear transformation of a CAR prior. It
is well known that any covariance matrix of a normal distribution can be
represented as a CAR prior [50], as well as that the normal distribution is
closed under affine transformations [43]. Consequently, it is possible to use
tools from Gaussian distribution theory to investigate the properties and
assumptions underlying the CAR-MM prior, such as its identifiability.

Despite both [42] and [23] including simulation studies to verify the
performance of the model, they did not investigate the identifiability and the
structure of the prior. In this paper, we develop the underpinning theory for
this methodology, and provide a more general framework able to encompass a
wider range of applications. In the case of [23] as well as this paper, the
weights are obtained by computing the share of population of a certain GP that
resides in a specific area. This approach can take into account the fact that
a single area can contribute to more than one membership (where each GP
practice constitutes a membership) as well that each membership is usually
defined in terms of more than one area, thus addressing the misalignment
issues mentioned above.

Instead of a spatial misalignment problem, defined as “different spatial data
layers are collected at different scales” [1, p. 165], one could consider it
as a change of support (COSP) problem. For instance, [6] start from the
subregion $D\subset\mathbb{R}^{d}$, the domain of interest, and partition it
according to $L\geq 1$ potentially non-nesting areal frameworks called source
supports i.e. the spatial support for which data is available, as opposed to
target supports where inference should take place. Due to the uniqueness
requirement, this approach is not be feasible given the structure of the data
that we work with. In fact, it would require the region of interest to be
partitioned into smaller areas, each of which would have to be assigned
uniquely to a specific GP. If one were to follow this approach, they would
fail to take into account the actual practice population structure. Obviously,
with individual level data this problem could be easily overcome, however
these are rarely available in disease mapping applications.

The paper is organised as follows: in Section 2 we introduce the CAR model and
its variant, the ICAR. In Section 3 we review the CAR-MM model, its definition
and distributional form, we compare the two proposals by [42] and [23] and
introduce a more compact notation in matrix form. The main contribution of
this paper is in Section 4 which (i) introduces a more general framework for
the MM models and discusses results of both uniqueness of the covariance
matrix and CAR specification for the CAR-MM prior, (ii) presents two different
parameterisations for the CAR-MM model and (iii) uses these results to
investigate the issue of identifiability in a Bayesian framework. Section 5 is
mostly concerned with introducing the methodologies used for the simulation
study that will follow and the model fit and model comparison methods that
will be used in the data analysis. In Section 6 simulation results are
introduced in reference to the considerations of Section 4. Lastly, in Section
7 we consider an example on a real South London diabetes prevalence dataset
providing an application of the CAR-MM prior. Three approaches are compared:
one where no spatial prior accounting for spatial autocorrelation is added,
and two others employing the CAR-MM and ICAR-MM priors.

##  2 Generalised Linear Models (GLM) and CAR priors

When modelling area referenced data, the region of interest is generally
divided into n areas for which the outcome of interest is observed. For each
area $i=1,...,n$ a random variable $Y_{i}$ is used to model the observed data,
whose mean is generally specified as follows:

$\displaystyle
g(E(Y_{i}|\bm{\theta},\bm{X}))=\eta_{i}(\bm{\theta}_{i},\bm{X})\quad\quad\forall
i=1,...,n,$ (2.1)

where $\eta_{i}$ is a linear function of the parameters $\bm{\theta}_{i}$ and
the p covariates $\bm{X}$, and g is a continuous monotonic function called the
link-function [36]. The considerations in the remainder of this paper are
valid for all models such as in (2.1) however, given that for areal data of
prevalence or mortality counts Poisson distribution is the most commonly
assumed, for ease of presentation, we will work under the following framework:

$\displaystyle\begin{split}Y_{i}&\sim Poisson(E_{i}\rho_{i})\\\
\log(E_{i}\rho_{i})&=\log{E_{i}}+\gamma+\bm{x}_{i}^{T}\bm{\beta}+\phi_{i},\end{split}$
(2.2)

where the link function in (2.1) is the natural logarithm, $E_{i}$ is the
expected number of cases for area i and the relative risk $\rho_{i}$ is
modelled as a linear function of the parameters
$\bm{\theta}_{i}=(\gamma,\bm{\beta},\phi_{i})$. In turn, the random effects
$\phi_{i}$ are modelled using a CAR prior which is a random vector
$\bm{\phi}=\left(\phi_{1},...,\phi_{n}\right)^{T}\in\mathbb{R}^{n}$ defined on
an undirected graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$.

More specifically, using matrices $\mathbf{M}$ and $\mathbf{C}$ that are
specified below in terms of the $\mathcal{G}$, we define the prior as:

$\displaystyle\bm{\phi}\sim\mathcal{N}\left(\bm{0},\mathbf{\Sigma}=(\mathbf{I}-\mathbf{C})^{-1}\mathbf{M}\right).$
(2.3)

[50] state four conditions to ensure the positive definiteness of the
covariance matrix $\mathbf{\Sigma}$

  1. C1:

$(\mathbf{I}-\mathbf{C})$ has positive eigenvalues

  2. C2:

$\mathbf{M}$ is diagonal with positive diagonal elements

  3. C3:

$(\mathbf{C})_{ii}=0$ for all $i=1,...,n$

  4. C4:

$(\mathbf{C})_{ij}/(\mathbf{M})_{ii}=(\mathbf{C})_{ji}/(\mathbf{M})_{jj}$ for
all $i,j=1,...,n$

In practice, a very common specification that fulfils C1 \- C4, imposes
$\mathbf{C}=\alpha\mathbf{D}^{-1}\mathbf{W}$, with $\alpha\in\left[0,1\right)$
and $\mathbf{W}$ the adjacency matrix such that $(\mathbf{W})_{ij}=1$ iff i
and j are neighbours and 0 otherwise, and
$\mathbf{M}=\tau^{-1}\mathbf{D}^{-1}$ with $\tau\in\mathbb{R}^{+}$ and
$(\mathbf{D})_{ii}=\sum_{j=1}^{n}(\mathbf{W})_{ij}$ for all i thus indicating
the number of neighbours of area i. Finally, we can write the CAR prior as:

$\displaystyle\bm{\phi}\sim$ $\displaystyle
CAR(\alpha,\tau,\mathbf{W})\quad,\quad\mathbf{\Sigma}^{-1}=\mathbf{Q}=\tau(\mathbf{D}-\alpha\mathbf{W}),$
(2.4)

[50] note that the representation of the CAR covariance matrix in (2.3) is
rather general, and in fact they prove it can represent any positive-definite
covariance matrix:

######  Theorem 2.1.

Any positive-definite covariance matrix $\boldsymbol{\Sigma}$ can be expressed
as the covariance matrix of a CAR model
$(\mathbf{I}-\mathbf{C})^{-1}\mathbf{M}$, for a unique pair of matrices
$\mathbf{C}$ and $\mathbf{M}$.

For a proof of this result see [50]

A special case of the CAR prior in equation 2.4, arises by setting the
parameter $\alpha=1$. This prior is called intrinsic conditional
autoregressive (ICAR) [3] and is widely used in practical applications. We
will employ this prior for the data analysis in Section 7. Its precision
matrix $\mathbf{Q}$ is singular consequently leading to an improper density.
In our implementation of the ICAR prior we impose the constraint
$\sum_{i}\phi_{i}=0$ at every MCMC iteration, as is commonly done in the
literature [18] to improve convergence and avoid intercept identifiability
issues [22].

##  3 Multiple membership and the CAR-MM prior

Both [42] and [23] proposed a GLM model with a CAR spatial prior and MM to
deal with misaligned data. We will now lay down the model that they developed
independently.

Let us consider a region partitioned in n areas and define a CAR prior on it
as per equation (2.4) and let
$\tilde{\bm{y}}=(\tilde{y}_{1},...,\tilde{y}_{m})$ be an outcome observed on m
units called memberships. Typically no areal spatial framework would be
available for the memberships, as the main goal of this modelling approach is
to use the n areas as a latent spatial structure.

For each membership j we define a set of n weights $0\leq h_{ji}\leq 1$,
corresponding to each area, such that $\sum_{i=1}^{n}h_{ji}=1$ for all
$j=1,..,m$. For example, the $h_{ji}$ might be the proportions of the
population of GP practice j living in areas i. We can thus re-define the
relative risk term of the GLM in (2.2) for outcome $\tilde{\bm{Y}}$ so that it
incorporates the latent areal structure we defined earlier:

$\displaystyle\begin{split}\tilde{Y}_{j}\sim
Poisson(&\tilde{E}_{j}\tilde{\rho}_{j}),\\\
\log(\tilde{\rho}_{j})=\sum_{i=1}^{n}h_{ji}\log(\rho_{i})=\sum_{i=1}^{n}&h_{ji}(\gamma+\bm{X}_{i}^{T}\bm{\beta}+\phi_{i}),\end{split}$
(3.1)

where $\tilde{E}$ indicate membership level offsets. This model defines each
membership level relative risk $\tilde{\rho}_{j}$ as a weighted average of the
areal log relative risks $\log(\bm{\rho})$. In this case each area can
contribute to multiple memberships with different weights.

We can more compactly represent the log relative risks using matrix notation,
for a $m\times n$ matrix $\mathbf{H}$ we define it element-wise using the
weights in (3.1) so that $h_{ji}=(\mathbf{H})_{ji}$. Additionally, we assume
$\mathbf{H}$ to be of full column rank and that for each area i there is at
least one membership j for which $h_{ji}>0$, i.e. all areas contribute to at
least one membership.

The second line in equation (3.1) can then be rewritten as:

$\displaystyle\log(\tilde{\bm{\rho}})=\mathbf{H}\log(\bm{\rho})=\gamma\mathbf{H}\mathbf{1}_{m}+\mathbf{H}\bm{X}\beta+\mathbf{H}\bm{\phi}=\gamma\mathbf{1}_{m}+\tilde{\bm{X}}\bm{\beta}+\tilde{\bm{\phi}},$
(3.2)

in line with [23] memberships corresponded to General Practitioner (GP)
practices, on which an outcome of interest was observed, and areas to English
census areal subdivisions (Middle Super-Output Area; MSOA) where the
covariates of interest were recorded. The MM matrix $\mathbf{H}$ was obtained
by computing the share of patients registered at a GP practice that resided in
a specific MSOA. [42] instead defined individual patients as memberships and
used membership level covariates, while the MM weights corresponded to the
share of time that patient j resided in a particular area i. Nevertheless,
both approaches can be represented with equation (3.2), by simply observing
that the intercept term $\gamma$ is unchanged by the linear transformation
$\mathbf{H}$ and so are coefficients $\bm{\beta}$.

##  4 Relationship between CAR and CAR-MM models

We will now focus on the CAR-MM prior and study its properties. Firstly, we
can observe that since the CAR-MM is a linear transformation of normally
distributed terms $\bm{\phi}$, $\bm{\tilde{\phi}}=\mathbf{H}\bm{\phi}$, it is
also normally distributed [43]. In fact:

$\displaystyle\bm{\tilde{\phi}}\sim\mathcal{N}_{m}(\bm{0},\tilde{\mathbf{\Sigma}}=\mathbf{H}\boldsymbol{\Sigma}\mathbf{H}^{T}),$
(4.1)

where $\mathbf{\Sigma}$ is as in equation (2.4). Additionally, in the case
where $m\leq n$, since we assumed full column rank for $\mathbf{H}$, the
resulting covariance matrix $\tilde{\mathbf{\Sigma}}$, with the corresponding
precision matrix $\tilde{\mathbf{Q}}$, will also be full rank (see Corollary
8.3.3 in [26]). Instead, when $m>n$, $\tilde{\mathbf{\Sigma}}$ will be of rank
$n$ and thus singular.

In the remainder of this section we explore the corresponding uniqueness
properties of the CAR-MM. As mentioned in Theorem 2.1, Section 2, [50] show
that any zero-mean Gaussian distribution on a finite set of points,
$\mathbf{Z}\sim N(0,\mathbf{\Sigma})$, can be written with a covariance matrix
parameterised as a CAR model.

Firstly, let us establish three conditions that will be common to the results
presented below:

  1. M1:

Let $\mathbf{\Sigma}$ be an $n\times n$ covariance matrix as defined in
Equation (2.4)

  2. M2:

Let $\mathbf{H}$ be a $m\times n$ matrix such that:

    1. (a)

for each row $j=1,...,m$ of $\mathbf{H}$: $\sum_{i=1}^{n}(H)_{ji}=1$,
$0\leq(H)_{ji}\leq 1$ for all i and j

    2. (b)

it is of full rank

    3. (c)

there is no column i such that $(\mathbf{H})_{\cdot i}=\bm{0}$

  3. M3:

Let $\tilde{\mathbf{\Sigma}}$ be the covariance matrix of the corresponding
CAR-MM model such that
$\tilde{\mathbf{\Sigma}}=\mathbf{H}\mathbf{\Sigma}\mathbf{H}^{T}$

where $(\mathbf{H})_{\cdot i}=\bm{0}$ indicates the i-th column of matrix
$\mathbf{H}$.

We remark that, although our approach is mainly focused on the CAR prior, the
following considerations directly apply to any Gaussian model. Therefore, the
spatial random effects in equation (3.1) could also be modelled through
Gaussian distributions with arbitrary covariance structure, for instance one
could use a Mátern covariance function using the Euclidean distances between
areas centroids [35, Example 6.1].

Above, following Theorem 2.1, we hinted at the fact that the CAR-MM can be
seen as a CAR prior. We now introduce a corollary (proved in the Supplemental
material) to that theorem, that draws a parallel between the two matrices
$\mathbf{M}$ and $\mathbf{C}$ that determine the CAR prior structure (C1-C4 in
Section 2), and the CAR-MM.

######  Corollary 1.

Let $\mathbf{\Sigma},\,\mathbf{H},\,\tilde{\mathbf{\Sigma}}$ be as defined in
M1-M3, then iff $m\leq n$, $\tilde{\mathbf{\Sigma}}$ can be expressed as the
covariance matrix of a CAR model
$(\mathbf{I}-\tilde{\mathbf{C}})^{-1}\tilde{\mathbf{M}}$ for a unique pair of
$m\times m$ matrices $\tilde{\mathbf{C}}$ and $\tilde{\mathbf{M}}$.

Since $\mathbf{C}$ contains the off-diagonal non-zero elements of the CAR
precision matrix, it determines the partial correlation structure of the
prior, i.e. its zero entries correspond to partially uncorrelated components
[45, Theorem 2.2]. Therefore, the same could be said for $\tilde{\mathbf{C}}$.
However, this matrix is not straightforward to obtain. In fact, the MM
transformation $\mathbf{H}$ is applied to the covariance matrix
$\mathbf{\Sigma}$ so that to obtain the precision matrix $\tilde{\mathbf{Q}}$
we could assume $m=n$, M1-M3 and, starting from (2.4) we could obtain:

$\displaystyle\tilde{\mathbf{Q}}=(\mathbf{H}(\mathbf{I}-\mathbf{C})^{-1}\mathbf{M}\mathbf{H}^{T})^{-1}=(\mathbf{H}^{T})^{-1}\mathbf{M}^{-1}(\mathbf{I}-(\mathbf{C}\mathbf{H}^{-1}-\mathbf{H}^{-1}+\mathbf{I})).$
(4.2)

Given that the resulting matrices $\tilde{\mathbf{M}}$ and
$\tilde{\mathbf{C}}$ would have to fulfil C1-C4 (Section 2) it is not
immediately clear how they would look in terms of $\mathbf{M}$, $\mathbf{C}$
and $\mathbf{H}$. In fact, if we took for example
$\tilde{\mathbf{M}}^{-1}=(\mathbf{H}^{T})^{-1}\mathbf{M}^{-1}$ and
$\tilde{\mathbf{C}}=\mathbf{C}\mathbf{H}^{-1}-\mathbf{H}^{-1}+\mathbf{I}$ we
would not have any guarantee that $\tilde{\mathbf{M}}^{-1}$ would be diagonal
as per C2, or $\tilde{\mathbf{C}}$ fulfils C3 (null diagonal elements) or C4
$(\tilde{\mathbf{C}})_{ij}/(\tilde{\mathbf{M}})_{ii}=(\tilde{\mathbf{C}})_{ji}/(\tilde{\mathbf{M}})_{jj}$
for all $i,j=1,...,n$. Therefore, differently from the spatial CAR model, we
have no way of obtaining analytically the partial correlation structure of the
CAR-MM prior.

Despite these difficulties regarding the precision matrix, we can consider the
covariance matrix, together with its transformation $\mathbf{H}$, and prove a
crucial result for model identifiability. The following lemma will help
evaluate under which circumstances the MM transformation of the CAR prior is
injective, and will also be used in the subsequent proofs.

######  Lemma 1.

Let $\mathbf{H}$ be a matrix as defined in M1-M3. $\mathbf{H}$ has a left
inverse $\mathbf{L}$ iff $m\geq n$, and this matrix is unique iff $m=n$.

It is worth mentioning that any left inverse $\mathbf{L}$ of $\mathbf{H}$ will
satisfy the same constraint, i.e the columns adding up to the $\mathbf{1}$
vector, in fact:
$\mathbf{H}\mathbf{1}_{n}=\mathbf{1}_{m}\,\Leftrightarrow\,\mathbf{L}\mathbf{H}\mathbf{1}_{n}=\mathbf{L}\mathbf{1}_{m}\,\Leftrightarrow\,\mathbf{1}_{n}=\mathbf{L}\mathbf{1}_{m}$
(proof in Section 1 of the Supplemental material.).

######  Theorem 4.1.

Let $\mathbf{\Sigma},\,\mathbf{H},\,\tilde{\mathbf{\Sigma}}$ be as defined in
M1-M3, then given $\tilde{\mathbf{\Sigma}}$ and $\mathbf{H}$,
$\mathbf{\Sigma}$ is uniquely defined iff $m\geq n$.

See Section 1 in the Supplemental material for a proof of both results.

###  4.1 Parameterisation

In light of the results presented above, we can devise two different
parameterisations for the random effects in equation (3.1): the post
parameterisation and the inverse parameterisation. In fact vectors $\bm{\phi}$
and $\tilde{\bm{\phi}}$ are linked deterministically by the linear
transformation $\mathbf{H}$ (see equation (3.2)), so the model has to specify
a prior only for one of the two.

The post-transform parameterisation consists of using the CAR prior directly
to sample random effects $\bm{\phi}^{*}$ from the spatial structure (equation
(2.4)), and then multiply them by $\mathbf{H}$ to obtain the MM equivalents
using $\mathbf{H}\bm{\phi}^{*}=\bm{\tilde{\phi}}^{*}$.

The inverse parameterisation is defined as follows. We specify a prior on
$\tilde{\bm{\phi}}^{\dagger}$ directly from the model in equation (4.1) and
obtain the corresponding spatial random effects using the inverse MM
transformation $\mathbf{H}^{-}\bm{\tilde{\phi}}^{\dagger}=\bm{\phi}^{\dagger}$

In terms of notation, all quantities associated with the post-transform
parameterisation will be denoted with an asterisk ($*$) superscript, while for
the inverse parameterisation a dagger ($\dagger$) will be used.

In order to be able to specify a proper CAR prior, we need the covariance
matrix to be positive definite. This means that due to Theorem 4.1, the
inverse parameterisation is available only in the case where $m\leq n$. Since
when m is strictly smaller than n by Lemma 1 there is no left inverse that
would allow us to obtain
$\mathbf{H}^{-}\bm{\tilde{\phi}}^{\dagger}=\bm{\phi}^{\dagger}$. Instead, the
post parameterisation can be implemented in all cases, given that it is only
necessary to specify a proper CAR prior on the areal structure as in Equation
(2.4).

It should be noted that whenever $m\neq n$ the specification of the CAR-MM
prior can be seen as a case of Low-rank representation of a spatial process
[19, Chapter 8], and specifically fixed rank kriging [14]. In fact, if we
consider the post parameterisation for instance, given that
$\mathbf{H}\bm{\phi}^{*}=\bm{\tilde{\phi}}^{*}$, for $m>n$ $\bm{\phi}^{*}$ is
a lower rank representation of $\tilde{\bm{\phi}}^{*}$, and vice versa for
$m<n$, as one random effect can be specified in terms of a deterministic
transformation of a lower rank distribution.

Low-rank representations can be particularly useful for geostatistical models
as they can significantly reduce the computational burden that is generally
required for spatial prediction. In this context matrix $\mathbf{H}$ takes the
name of basis and its selection, and in some cases estimation, is generally
arbitrary and depends on the problem at hand. This is different from the CAR-
MM case where the MM matrix $\mathbf{H}$ is entirely determined by the
population structure of the areas and memberships. Additionally, contrary to
most geostatistical applications in this context where the lower rank process
is simply used as a latent factor, for the CAR-MM both the membership and the
areal level, i.e. both the lower and higher rank process can be of inferential
interest.

###  4.2 Identifiability

We can now investigate the identifiability of the model by introducing two
definitions of non-identifiability in a Bayesian context and evaluate each one
in the case of CAR-MM random effects. We will not provide explicit
considerations for the hyperparameters $\alpha$ and $\tau$, as they are
generally not parameters of interest.

Firstly, according to [16]:

######  Definition 1 (Likelihood non-identifiability).

For a Bayesian model $M(\bm{\theta},\bm{y})$ with data $\bm{y}$, likelihood
$L$ and prior $f$, if the parameter vector $\bm{\theta}$ is partitioned into
$\bm{\theta}=\left(\bm{\theta}_{1},\bm{\theta}_{2}\right)$ then
$\bm{\theta}_{2}$ is non-identifiable if:

$\displaystyle
f(\bm{\theta}_{2}\mid\bm{\theta}_{1},\bm{y})=f(\bm{\theta}_{2}\mid\bm{\theta}_{1}),$
(4.3)

[20] remark that Definition 1 is equivalent to the classical notion of non-
identifiability of the likelihood (see Definition 1 Section 3.1.1 in [10])
since:

$\displaystyle f(\bm{\theta}_{2}\mid\bm{\theta}_{1},\bm{y})\propto
L(\bm{\theta}_{1},\bm{\theta}_{2}\mid\bm{y})f(\bm{\theta}_{2}\mid\bm{\theta}_{1})f(\bm{\theta}_{1}),$
(4.4)

is not identifiable according to Definition 1 iff
$L(\bm{\theta}_{1},\bm{\theta}_{2}\mid\bm{y})$ can be re-parameterised so that
it becomes free of $\bm{\theta}_{2}$.

In addition, it is worth recalling that [9] define parameter redundancy as

######  Definition 2 (Parameter redundancy).

A model $M(\bm{\theta},\bm{y})$ is parameter redundant if we can write
$\bm{\theta}$ as a function of just $\bm{\beta}$, where
$\bm{\beta}=g(\bm{\theta})$ and $dim(\bm{\beta})<dim(\bm{\theta})$, where dim
is length of the vector.

They also prove that parameter redundancy implies non-identifiability in the
likelihood which is then equivalent to Definition 1 (Theorem 4 in [9]). We
refer back to the model in equation (4.1) and conditions M1-M3 imposed on
$\mathbf{H}$ in Section 4 and evaluate the two parameterisations under the
definitions we just presented.

In the case of a post-parameterised model, it is evident that when $m>n$ the
MM random effects $\tilde{\bm{\phi}}^{*}$ are parameter redundant as they can
be written as a function of the sole areal random effects $\bm{\phi}^{*}$
which are n dimensional. More specifically, there is an injective function $g$
such that $g(\bm{\phi}^{*})=\mathbf{H}\bm{\phi}^{*}=\tilde{\bm{\phi}}^{*}$,
where $g:\mathbb{R}^{n}\to\mathbb{R}^{m}$.

This can also be seen more explicitly, by partitioning the model parameters in
equation (3.1) into $(\bm{\phi},\tilde{\bm{\phi}},\bm{\theta})$, where
$\bm{\theta}=(\gamma,\bm{\beta},\alpha,\tau)$, we can write the posterior
distribution, up to a constant, as

$\displaystyle\begin{split}L(\tilde{\bm{y}}|\tilde{\bm{\phi}}^{*},\bm{\phi},\bm{\theta},\mathbf{H},\bm{X})f(\tilde{\bm{\phi}}^{*},\bm{\phi}^{*},\bm{\theta})&=L(\tilde{\bm{y}}|g(\bm{\phi}^{*}),\bm{\phi},\bm{\theta},\mathbf{H},\bm{X})f(g(\bm{\phi}^{*})\bm{\phi}^{*},\bm{\theta})\\\
&=L(\tilde{\bm{y}}|g(\bm{\phi}^{*}),\bm{\theta},\mathbf{H},\bm{X})f(\bm{\phi}^{*},\bm{\theta})\end{split}$
(4.5)

where the second equality comes from the fact that the $\tilde{\bm{\phi}}^{*}$
are a deterministic function exclusively of the $\bm{\phi}^{*}$ and thus do
not depend on either $\bm{\theta}$ or $\\{\mathbf{H},\bm{X}\\}$. Consequently,
in (4.5) we can see that the MM random effects $\tilde{\bm{\phi}}$ are non-
identifiable in the likelihood (Definition 1). In the case of the post
parameterisation, this holds true for all m and n.

With regard to the inverse parameterisation, where a prior is specified on the
MM random effects $\tilde{\bm{\phi}}^{\dagger}$, in the case $m\geq n$ by
Lemma 1 there is a left inverse such that
$\bm{\phi}^{\dagger}=\mathbf{L}\tilde{\bm{\phi}}^{\dagger}$. Thus the same
considerations hold as for the post parameterisation made above. However in
this case the areal and MM random effects exchange statuses in terms of
identifiability and redundancy. To see this, suffices to swap $\bm{\phi}^{*}$
for $\tilde{\bm{\phi}}^{\dagger}$ and k for g in (4.5), where
$\mathbf{L}\tilde{\bm{\phi}}^{\dagger}=\bm{\phi}^{\dagger}=k(\tilde{\bm{\phi}}^{\dagger})$
and $k:\mathbb{R}^{m}\to\mathbb{R}^{n}$. In the $m<n$ case for the inverse
parameterisation, by Theorem 4.1 the areal random effects covariance matrix
$\mathbf{\Sigma}^{\dagger}$ is not uniquely defined and thus the
$\bm{\phi}^{\dagger}$ are not identifiable in the likelihood.

Despite these issues with the likelihood, we can adopt a different perspective
following a different definition, by [10] (Definition 12):

######  Definition 3 (Posterior non-identifiability).

A Bayesian model $M(\bm{\theta},\bm{y})$ described by the posterior
distribution, is globally identifiable if
$f(\bm{\theta}_{a}\mid\bm{y})=f(\bm{\theta}_{b}\mid\bm{y})$ implies
$\bm{\theta}_{a}=\bm{\theta}_{b}$. A model is locally identifiable if there
exists an open neighbourhood of any $\bm{\theta}$ such that this is true.
Otherwise the posterior is non-identifiable.

Under this perspective we can see that both membership and areal random
effects are identifiable in the posterior for the post parameterisation when
$m\geq n$. In fact, by setting
$\bm{\Delta}=\\{\tilde{\bm{y}},\mathbf{H},\bm{X}\\}$, we can investigate
identifiability by writing down the posterior distribution as:

$\displaystyle\begin{split}f(g(\bm{\phi}^{*}_{a})\mid\bm{\phi}^{*}_{a},\bm{\theta}_{a},\bm{\Delta})f(\bm{\phi}^{*}_{a},\bm{\theta}_{a}\mid\bm{\Delta})&=f(g(\bm{\phi}^{*}_{b})\mid\bm{\phi}^{*}_{b},\bm{\theta}_{b},\bm{\Delta})f(\bm{\phi}^{*}_{b},\bm{\theta}_{b}\mid\bm{\Delta}).\end{split}$
(4.6)

In the $m\geq n$ case, we can easily see that since the function $g$ is
injective, we have that if and only if $\bm{\phi^{*}}_{a}=\bm{\phi^{*}}_{b}$
then $g(\bm{\phi}^{*}_{a})=g(\bm{\phi}^{*}_{b})$. In turn, this entails that
for the purpose of verifying Definition 3, (4.6) reduces to:

$\displaystyle
f(\bm{\phi}^{*}_{a},\bm{\theta}_{a}\mid\bm{\Delta})=f(\bm{\phi}^{*}_{b},\bm{\theta}_{b}\mid\bm{\Delta})$
(4.7)

where the terms $\bm{\phi}^{*}$ are simple CAR random effects.

Therefore, for $m\geq n$ Equation 4.7 shows that a post parameterised CAR-MM
model as in (3.1), is identifiable if and only if the underlying CAR model
(2.2) is as well. Conversely, we remark that in the $m<n$ case, given that g
is not injective, (4.6) does not hold in general, thus making the random
effects non-identifiable in the posterior.

In the case of the inverse parameterisation for $m>n$, Equation 4.6 remains
the same with simply the $\bm{\phi}^{*}$ swapped for
$\tilde{\bm{\phi}}^{\dagger}$ and k for g. However, in this instance the
function k is not injective, thus leading to non-identifiability. Furthermore,
for $m<n$ there is no function k, given that $\mathbf{H}$ has no left inverse.
Only the $m=n$ case is identifiable under both parameterisations.

With regard to the covariates coefficients $\bm{\beta}$, it is apparent from
equation 3.2 that the MM transformation of the log relative risks
$\log(\bm{\rho})$ does not affect these parameters. Therefore, the question of
identifiability of both the areal and membership relative risks hinges on the
random effects $\bm{\phi}$ and $\tilde{\bm{\phi}}$ as discussed above.

To conclude, we recommend the use of the CAR-MM model only when a $m\geq n$
dataset is available. Only in this case, in fact, parameters of interest, such
as the relative risks as well as risk factor coefficients will be
identifiable.

##  5 Model validation and posterior checks

###  5.1 Simulation based Calibration (SBC)

In order to both check the software implementation of our posterior sampler
and verify that the model is able to correctly recover the true parameters, we
devised a simulation study following a simulation based calibration (SBC)
approach [46]. In an SBC study, instead of comparing posterior samples to some
fixed ground truth parameter, the procedure assesses the overall calibration
of the posterior samples with respect to the joint distribution
$\pi(\bm{y},\bm{\theta})$. In this way, it is possible to verify the
performance of the posterior sampler for multiple ground truths
$\tilde{\bm{\theta}}$ generated from the prior distribution.

In summary, the SBC procedure can be written as follows:

$\displaystyle\tilde{\bm{\theta}}\sim\pi(\bm{\theta}),\quad\tilde{\bm{y}}\sim
l(\bm{y}|\tilde{\bm{\theta}}),\quad\\{\bm{\theta}^{(1)},...,\bm{\theta}^{(B)}\\}\sim\pi(\bm{\theta}|\tilde{\bm{y}}).$
(5.1)

Since we have $\pi(\bm{y},\bm{\theta})=\pi(\bm{\theta})l(\bm{y}|\bm{\theta})$,
by sampling a parameter vector $\tilde{\bm{\theta}}$ from the prior and then
generate new data from the likelihood $\tilde{\bm{y}}\sim
l(\bm{y}|\tilde{\bm{\theta}})$, this would be equivalent to drawing a sample
from the joint distribution
$(\tilde{\bm{y}},\tilde{\bm{\theta}})\sim\pi(\bm{y},\bm{\theta})$. Therefore,
whenever we implement an algorithm that correctly samples from a posterior
distribution, the $\\{\bm{\theta}^{(1)},...,\bm{\theta}^{(B)}\\}$ are going to
have the same marginal distribution as $\tilde{\bm{\theta}}$ sampled from the
prior [13].

On the basis of this, we can straightforwardly check that the implemented
algorithm performs inference correctly by verifying that prior and posterior
samples are distributed according to the same distribution [46]: for any one-
dimensional random variable $f:\bm{\Theta}\to\mathbb{R}$, the rank statistic
$r(\cdot)$ of the prior sample relative to the posterior samples:

$\displaystyle
r(\left\\{f(\bm{\theta}^{(1)}),...,f(\bm{\theta}^{(B)})\right\\},f(\tilde{\bm{\theta}}))=\sum_{b=1}^{B}\mathbb{I}\left[f(\bm{\theta}^{(b)})<f(\tilde{\bm{\theta}})\right]\in\left[0,B\right],$
(5.2)

will be uniformly distributed across the integers $\left[0,B\right]$.
Therefore, by sampling from the Bayesian joint distribution N times, and for
each sample drawing B posterior samples, computing the rank statistic in
equation (5.2) and finally checking it is actually uniformly distributed, we
can build a validation procedure for any Bayesian model. We mention that in
practice, given that we have both N and B very large, we divide the sampled
ranks by $B+1$ and verify calibration with the Uniform distribution on $[0,1]$

###  5.2 Posterior checks and model comparison measures

Before illustrating the simulation results we introduce the methodologies used
to assess both model fits and comparisons. Posterior predictive p-values (PPP)
are a widely used tool in Bayesian analysis to compare the estimated posterior
distribution with the observed data [21]. In the case of marginal predictive
checks, model calibration can be verified in two steps: first generating new
observations $\bm{y}^{rep}$ using posterior samples and secondly by computing
the marginal PPP’s, for each observation
$\bm{y}=\left(y_{1},...,y_{m}\right)^{T}$:

$\displaystyle
p_{j}=P\left(T(y_{j}^{rep})<T(y_{j})\mid\bm{y}\right)+\frac{1}{2}P\left(T(y_{j}^{rep})=T(y_{j})\mid\bm{y}\right),$
(5.3)

if $\bm{y}$ is a discrete outcome, where $T(\cdot)$ is a generic test
quantity; in our case we will take the identity function $T(y_{j})=y_{j}$.

A computationally inexpensive alternative is the one from [34] where they
propose a method to approximate this “leave-one-out” cross-validatory model
fit measure: $p_{j|-\bm{j}}=P\left(y_{j}^{rep}\leq
y_{j}\mid\bm{y}_{-j}\right)$, which for continuous data will have a uniform
distribution if the model is well calibrated. Referring to Equation 2.2, we
can simulate a new random effect for area i using $\phi_{i}\sim
p(\phi^{rep}_{i}\mid\alpha,\tau,\bm{\phi}_{-i})$ using posterior samples of
$\alpha$ and $\tau$ and thus approximate the $p_{j|-\bm{j}}$’s above. This
approach can be directly extended to the CAR-MM model, by simply post-
multiplying the quantities in the relative risk by the MM matrix $\mathbf{H}$
as per equation (3.1).

With regard to model comparison, in the data analysis Section we compare
different models fit using three criteria. The first one is Pareto smoothed
importance-sampling leave-one-out cross-validation (PSIS-LOO) [47]. Similarly
to the approach described above for PPP, this method make use of a conditional
independence assumption of data given model parameters $\bm{\theta}$, which in
our case is implicit in Equations (2.2) and (3.1). The PSIS-LOO approach
allows to compute the expected log pointwise predictive density (elpd) to
estimate out-of-sample predictive fit: $elpd_{loo}=\sum_{j=1}^{m}\log
p(y_{j}\mid\bm{y}_{-j})$. Without the need to fit the model m times, estimates
of $elpd_{loo}$ can be obtained using

$\displaystyle\widehat{elpd}_{psis-
loo}=\sum_{j=1}^{m}\log\left(\frac{\sum_{s=1}^{S}w_{j}^{s}p(y_{j}\mid\bm{\theta}^{s})}{\sum_{s=1}^{S}w_{j}^{s}}\right),$
(5.4)

where $s=1,...,S$ indexes the posterior sample iterations, and $w_{j}^{s}$ the
Pareto smoothed importance weights. Model comparison can be carried out by
computing the difference in $\widehat{elpd}_{loo}$ between the two models, so
for models A and B:
$\text{diff}(\widehat{elpd}_{loo}^{A,B})=\widehat{elpd}_{loo}^{A}-\widehat{elpd}_{loo}^{B}$,
where a larger $\widehat{elpd}_{loo}$ indicates better estimated predictive
accuracy of the model. The standard error of the difference is thus:

$\displaystyle se(\text{diff}(\widehat{elpd}_{loo-
psis}^{A,B}))=\sqrt{mV_{j=1}^{m}\left(\text{diff}(\widehat{elpd}_{loo-
psis}^{A,B})\right)}.$ (5.5)

where $V_{j=1}^{m}(a_{j})$ is the sample variance of m values with mean
$\bar{a}$:
$V_{s=1}^{S}(a_{s})=\frac{1}{m-1}\sum_{j=1}^{m}(a_{j}-\bar{a})^{2}$. The
paired estimate of equation (5.5) relies on the two models being fitted on the
same set of m data points.

The $\widehat{elpd}_{loo}$ (we drop “psis” in the subscript for conciseness)
is easily computable using the R package `loo` [48] in conjunction with the
RStan package [8] for sampling from the posterior distribution. The authors
suggest checking the estimated shape parameters $\hat{k}$ of the Pareto
weights for all observation to assess reliability of estimates: when
$\hat{k}>0.7$ there might be stability issues.

Additionally, we also considered the mean Ranked Probability Score for count
data ($\overline{RPS}$) [15]

$\displaystyle\overline{RPS}=\frac{1}{m}\sum_{j=1}^{m}\left[\frac{1}{B}\sum_{i=1}^{B}\absolutevalue{y^{rep(i)}_{j}-y_{j}}-\frac{1}{B}\sum_{i=1}^{B/2}\absolutevalue{y^{rep(i)}_{j}-y^{rep(i+B/2)}_{j}}\right],$
(5.6)

where $y^{rep(i)}_{j}$ is the i-th replicate of observation $y_{j}$. Finally,
the Dawid-Sebastiani Score ($\overline{DSS}$) [17] was also included in our
comparisons

$\displaystyle\overline{DSS}=\frac{1}{m}\sum_{j=1}^{m}\left[\left(\frac{y_{j}-\overline{y^{rep}_{j}}}{\sigma_{y_{j}^{rep}}}\right)^{2}+2\log\sigma_{y_{j}^{rep}}\right],$
(5.7)

where $\overline{y^{rep}_{j}}$ and $\sigma_{y_{j}^{rep}}$ are respectively the
mean and standard deviations of the replicates for the j-th observation. For
both $\overline{RPS}$ and $\overline{DSS}$ a lower value indicates better fit
for the model.

##  6 Simulation study

We now turn to a simulation study designed to assess the considerations on
identifiability of Section 4. We designed the simulations to evaluate the two
parameterisations (post and inverse) and the number of memberships relative to
areas, i.e. $m<n$, $m=n$ and $m>n$. We will first describe the implementation
of the posterior sampler for each model, then the simulation setup in detail
and finally the results.

###  6.1 Model Implementation

All models are implemented in RStan with significant differences between the
two parameterisations. For the post parameterised covariance matrices, we
followed [30]’s sparse representation combined with the approach suggested by
[29] for the determinant. The former relies on a sparse representation of the
matrix $\mathbf{W}$ whereby the terms corresponding to non-neighbouring areas
are from the computation. With regard to the determinant, following [29] we
can avoid inverting the entire precision matrix $\mathbf{Q}$ at each
iteration. Taking $\tau=1$, without loss of generality, it can be shown that:
$\mathbf{Q}=\mathbf{D}-\alpha\mathbf{W}=\mathbf{D}^{1/2}(\mathbf{I}-\alpha\mathbf{D}^{-1/2}\mathbf{W}\mathbf{D}^{-1/2})\mathbf{D}^{1/2}$
then $\det(\mathbf{Q})=\det(\mathbf{D})\prod_{i=1}^{n}(1-\alpha\lambda_{i})$,
where the $\lambda_{i}$’s are the eigenvalues of
$\mathbf{D}^{-1/2}\mathbf{W}\mathbf{D}^{-1/2}$. This representation
significantly speeds up the computation of $det(\mathbf{Q})$ since the
$\lambda_{i}$’s do not depend on either $\alpha$ and $\tau$ and thus have to
be calculated only once.

Conversely, for the inverse parameterisation we specified the entire $m\times
m$ precision matrix $\tilde{\mathbf{Q}}$. We remark that in this instance, the
sampler can only be applied to the $m\leq n$ case, as it is the only case
where the covariance matrix is proper. Under this parameterisation, in order
to obtain the areal random effects $\bm{\phi}^{\dagger}$ we used the Moore-
Penrose inverse of the $\mathbf{H}$ matrix.

The code for both the simulation study and the RStan models can be found at
https://github.com/markgrama/CAR-MM.

With regard to the model, we used the following priors (see equation (3.1)): a
uniform on $(0,1)$ for $\alpha$, a $Gamma(2,0.2)$, parameterised in terms of
shape and rate, for the marginal precision $\tau$, while for $\gamma$,
$\beta_{1}$ and $\beta_{2}$ three independent $\mathcal{N}(0,0.7^{2})$
distributions were used.

In terms of computational time, with the same number of iterations (see
Section 6.2 for details), the post parameterisation performed significantly
better than the inverse, with an average runtime of 2 minutes across
membership sizes, against the 10 minutes for the $m=70$ case and 29 minutes
for $m=100$ (Table 1 in the Supplemental Material). Our hypothesis is that the
post parameterisation runs more quickly due to the fact that the sampling of
the CAR random effect is done through the sparse representation mentioned
above, and the linear transformation $\mathbf{H}$ is only evaluated in the
likelihood with RStan optimised matrix operations.

Conversely, the inverse parameterisation requires the evaluation of the entire
precision matrix, since no adjacency matrix is directly available. This is due
to the problem in identifying the partial correlation structure of the CAR-MM
mentioned in Section 4. Additionally, even if we were able to determine
partial independence relationships between components, there is no guarantee
that the resulting precision matrix will be sparse.

###  6.2 Setup

In this simulation study we use the post and inverse parameterisations to both
generate the data and sample from the posterior, thus obtaining 4 different
scenarios. In order to refer to a particular scenario we might write post-
inverse, where the first parameterisation refers to the data generation
(post), while the second indicates the MCMC sampler parameterisation
(inverse). For each scenario we consider three membership sizes ($70,100,130$)
so that we can correspondingly evaluate the three cases $m<n$, $m=n$, and
$m>n$. However, as we remark in Section 4.1, the inverse parameterisation is
available only for $m\leq n$, therefore we simulate the $m=130$ case only for
the post parameterisation.

The spatial structure used in the study is a $10\times 10$ grid, hence
$n=100$. Expected rates for each memberships ($\tilde{E}_{j}$ in equation
(3.1)) were simulated from $m$ independent $Poisson(20)$ distributions. For
both the MM matrix and the $\tilde{\bm{E}}$ vector, we simulated the case
$m=130$ and removed the extra rows or elements to obtain the other membership
sizes.

With regard to the MM matrix $\mathbf{H}$, this was simulated starting from
each area’s neighbourhood structure, encoded in the adjacency matrix
$\mathbf{W}$, up to second order neighbours. For each area a new membership
was created by first generating random weights from an uniform distribution to
all first and second order neighbours. Then those quantities were normalised
so that the area itself would amount to 50% of the total weights, the sum of
all first order neighbours would be attributable 35% of the total weight, and
the remaining 15% would be taken up by the second order neighbours. More
formally, the simulated $m\times n$ matrix $\mathbf{H}$ is obtained by first
randomly drawing samples:

$\displaystyle\begin{split}(G)_{ji}\overset{iid}{\sim}\mbox{Uniform}(0,1)\quad&\mbox{for}\quad
i=1,...,n\text{ and }j=1,...,m\\\
\mbox{s.t.}\quad\\{w_{ji}=1\\}&\vee\\{w_{il}=1\wedge w_{lj}=1\\},\end{split}$
(6.1)

for $l=1,...,n$. Therefore the elements of the matrix $\mathbf{H}$ are defined
as follows

$\displaystyle(H)_{ji}=\begin{cases}0.5&\mbox{if }i=j,\\\
\frac{(G)_{ji}}{\sum_{j:w_{ji}=1}(G)_{ji}}\cdot 0.35&\mbox{if }w_{ji}=1,\\\
\frac{(G)_{ji}}{\sum_{j:\\{w_{il}=1\wedge w_{lj}=1\\}}(G)_{ji}}\cdot
0.15&\mbox{if }w_{il}=1\mbox{ and }w_{lj}=1.\end{cases}$ (6.2)

In order to obtain a 130 membership matrix, the extra rows were created by
simply repeating the same process and generating different random weights.
Before deleting the bottom rows to obtain the varying number of memberships,
in order to ensure that all areas were represented as uniformly as possible in
the final matrix its rows were randomly shuffled. Finally, the obtained 70
membership matrix was checked to ensure that each area was represented with at
least one weight. In our case the area with the least weights had 2, so it
contributed to 2 memberships.

A total of $10^{4}$ datasets were simulated for the memberships structure
using the post parameterisation, each membership size was obtained by simply
truncating the simulated outcome vector accordingly, thus increasing the total
of datasets used to 9000. Covariates $\bm{X}$ were kept fixed across all
datasets and they were generated using two independent standard normal samples
and min-max normalised. Both the MM matrix and offsets were also the same in
all simulated datasets.

Following the approach described in Section 5.1, hyperparameters for the CAR
prior where sampled from the same priors described in Section 6.1 as well as
Hamiltonian Monte Carlo (HMC) [4] sampling which was run using RStan [8] on
Queen Mary’s Apocrita High Performance Computing facility, supported by QMUL
Research-IT [31], using five $10^{4}$ iterations chains and 50% warm-up. To
reduce autocorrelation, as suggested in [46], samples were thinned by 5
leading to $L=5000+1$ posterior samples (see Section 5.1).

Additionally, in order to ensure convergence of the posterior sampler, we only
keep those simulations for which no parameter report a $\hat{R}$ statistic
over $1.01$, following [49]. In Table 2 in the Supplemental material we report
the number of simulations that were excluded. Convergence was satisfactory
overall, since we always retained at least 98% of the simulations.

Finally we evaluate bias and root mean square error (RMSE) across simulations
by computing the following quantities. For a generic parameter $\theta$ with
true value $\hat{\theta}_{k}$ for simulated dataset $k$, $k=1,...,10^{4}$ we
compute

$\displaystyle\text{Bias}_{k}=\frac{1}{B}\sum_{b=1}^{B}\theta^{(b)}_{k}-\hat{\theta}_{k}\,\,;\,\,\text{Abs
bias}_{k}=\frac{1}{B}\sum_{b=1}^{B}\absolutevalue{\theta^{(b)}_{k}-\hat{\theta}_{k}}\,\,;\,\,\text{RMSE}_{k}=\sqrt{\frac{1}{B}\sum_{b=1}^{B}(\theta^{(b)}_{k}-\hat{\theta}_{k})^{2}},$
(6.3)

where $\theta^{(b)}_{k}$ is the $b^{th}$ HMC sample for the $k^{th}$ dataset
and $B=2.5\cdot 10^{4}$.

###  6.3 Results

| Data generation: Post | Data generation: Inverse
---|---|---
| MCMC: Post | MCMC: Inverse | MCMC: Post | MCMC: Inverse
| $70$ | $100$ | $130$ | $70$ | $100$ | $70$ | $100$ | $70$ | $100$
$\alpha$ | 87.35 | 66.62 | 80.61 | 37.75 | 73.45 | 43.16 | 76.44 | 95.80 | 71.72
$\tau$ | 68.52 | 93.80 | 94.43 | 8.92 | 93.36 | 30.82 | 90.85 | 91.47 | 91.54
$\gamma$ | 98.38 | 99.42 | 98.60 | 99.70 | 99.23 | 99.96 | 99.52 | 99.95 | 99.57
$\beta_{1}$ | 96.15 | 96.14 | 83.34 | 99.79 | 97.76 | 99.99 | 96.41 | 99.96 | 97.25
$\beta_{2}$ | 99.97 | 99.96 | 99.75 | 99.96 | 99.77 | 99.95 | 99.52 | 99.81 | 99.46
$\bm{\phi}$ | 94.63 | 94.33 | 94.22 | 50.72 | 94.53 | 57.72 | 94.53 | 93.83 | 94.54
$\bm{\rho}$ | 94.45 | 94.37 | 93.51 | 26.35 | 94.68 | 59.30 | 94.46 | 27.95 | 94.17
$\tilde{\bm{\rho}}$ | 95.37 | 95.78 | 94.72 | 94.32 | 95.78 | 94.38 | 95.72 | 95.16 | 95.88
Table 1: Coverage of the $95\%$ expected variation intervals of the SBC
obtained posterior samples of the CAR-MM parameters for each scenario and
membership size

In Table 1 we report the coverage of the $95\%$ intervals of the expected
variation of the SBC rank statistics against the uniform distribution 111The
plots for all parameters can be found in Section 2.2 in the Supplemental
material. For the $m=100,130$ cases, we observe high coverage for all the
parameters in all scenarios, with the exception of the $\alpha$ parameter,
which is generally not of inferential interest. We believe that the relatively
lower coverage of $\beta_{1}$ in the $m=130$ scenario is not a cause of
concern for two reasons: first, as will be seen below, posterior means are
unbiased across the simulation study, and second, this fact does not seem to
hinder the good calibration of the relative risks $\bm{\rho}$ and
$\tilde{\bm{\rho}}$.

Overall, in the $m=70$ case, we have satisfactory coverage for the linear term
parameters $\gamma$, $\beta_{1}$ and $\beta_{2}$ in all scenarios. Despite
this, we observe poor calibration for the spatial random effects $\bm{\phi}$
and the areal relative risks $\bm{\rho}$, the only exception being the post-
post case. These results are consistent with the considerations on
identifiability that we discussed in Section 4 where we remarked the
identifiability only of the $m=n$ case under both parameterisation and the
$m>n$ for the post one. Although theoretically the $m=70$ should not be
identifiable, the calibration in the post-post case could be due to the
coincidence between data generation and the MCMC parameterisation. Conversely,
in the inverse-inverse case, despite the two parameterisations being the same,
we remark that the generalised inverse used to obtain areal relative risks
($\bm{\rho}$) from the sampled membership risks ($\tilde{\bm{\rho}}$) is not a
right inverse (Lemma 1). This leads to the poor calibration of the $\bm{\rho}$

All in all, these simulations suggest that in the post parameterisation case
the posterior samples are well calibrated and support the considerations on
identifiability of the CAR-MM model discussed above. Overall these results
lead to us recommend using the post parameterisations in all applications, not
only due to it computational efficiency as mentioned in Section 6.1, but also
due to its overall better calibration.

However, despite the positive results for the $m<n$ case, we remind the reader
that in Section 4.2 we showed its non identifiability in the posterior, under
both parameterisations. The good calibration results in the post-post (both
data and MCMC sampler are post parameterised) case could simply an artefact of
the simulation setup. In fact, in that scenario the function that maps the
areal relative risks to the membership risks, i.e. the MM matrix $\mathbf{H}$
is the same both when simulating the data and sampling from the posterior
distribution. Conversely, whenever an inverse parameterisation is involved,
either in the data simulation stage or during MCMC sampling, given that a left
inverse does not exist (Lemma 1), the data does not contain enough information
to correctly recover the spatial random effects at the areal level. Finally,
high calibration for all parameters across different parameterisations for the
two $m\geq n$ cases is evidence of the reliability of the CAR-MM model.

##  7 Data analysis

###  7.1 Data structure and modelling

We present an application of this methodology to primary care catchment areas.
As argued above, the CAR-MM approach is of utility in situations where data
are observed for non-spatial units (memberships) but where an underlying
spatial framework can be conceived, and where an area to membership population
cross-referencing is available [28]. In particular, UK primary care data on
disease prevalence (which is recorded at GP practice level) offers a suitable
application. Specifically, we consider diabetes prevalence for a set of GP
practices in South East London (SEL), and to the MSOAs (small areas) defining
the catchment areas of these GP practices.

The data pertain to 153 GP practices which together constitute what are known
as Clinical Commissioning Groups or CCGs: the 153 practices are those for
Bexley, Lewisham, Greenwich and Bromley CCGs. Based on the choice of GPs, a
total of 152 MSOAs were included so that at least 90% coverage of each GPs
population was reached.

Prevalence data for GPs refers to October 2015, as does the cross reference
information providing the number of each MSOA residents registered at each GP
practice [27]. Age structure information was also available at GP level, so
that expected rates could be easily computed using England and Wales national
rates of diabetes prevalence.

For MSOA level covariates, two were included: from the 2011 Census the share
of South Asian residents $(\bm{X})_{\cdot 1}$ i.e. Indian, Pakistani and
Bangladeshi, see Table KS201EW in the NOMIS database [41], and the 2015 Index
of Multiple deprivation (IMD) $(\bm{X})_{\cdot 2}$. The two variables were
included as a control of socio-economic status and due to the higher
prevalence of diabetes among the South Asian population [25]. Both variables
were min-max normalised.

###  7.2 Results

Initial runs using the Poisson likelihood model (Eq. 3.1) returned concerning
PSIS-LOO weights $\hat{k}$ which might be a symptom of misspecification of the
model. In fact, [47] note that high $\hat{k}$ mean that the full posterior
distribution is a poor approximation of the leave-one-out posterior.

Since the possible misspecification might be attributable to overdispersion,
i.e. when the variance of the outcome is larger than its average [11] we
therefore turned to a Negative Binomial likelihood for the model in equation
3.1, that allows for direct modelling of overdispersion.

The density of the parameterisation we used is the following:

$\displaystyle f_{Y}(y_{j};\mu_{j},\psi)=\begin{pmatrix}y_{j}+\psi-1\\\
y_{j}\end{pmatrix}\left(\frac{\mu_{j}}{\mu_{j}+\psi}\right)^{y_{j}}\left(\frac{\psi}{\mu_{j}+\psi}\right)^{\psi},$
(7.1)

where $y_{j}\in\mathbb{R}$ for all $j=1,...,m$;
$E(Y_{j})=\mu_{j}\in\mathbb{R}^{+}$, $\psi\in\mathbb{R}^{+}$ and
$Var(Y_{j})=\mu_{j}+\mu_{j}^{2}/\psi$. Therefore the mean of each $Y_{j}$ can
be modelled directly through the $\mu_{j}$’s. Equation 3.1 then becomes:

$\displaystyle\tilde{Y}_{j}\sim
NegBin(\tilde{E}_{j}\tilde{\rho}_{j},\psi)\quad,\quad\log(\tilde{\bm{\rho}})=\log(\bm{\mathbf{H}\rho})=\mathbf{H}\left(\gamma\mathbf{1}_{n}+\bm{X}\bm{\beta}+\bm{\phi}\right),$
(7.2)

so the overdispersion parameter $\psi$ is assumed common across all GP
practices (i.e. memberships). The priors used were the same as for the
simulation study in Section 6.

In the analysis we compared three models: the CAR-MM as described above, the
ICAR-MM and a GLM-MM that simply consists in removing the spatial random
effects $\bm{\phi}$ from equation 7.2. We include the ICAR since it is one of
the most widely used prior in the disease mapping literature [35]. All models
were coded in RStan (see Section 6.1 and [38] for the ICAR prior) and run for
five $10^{5}$ iterations chains with 50% warm-up. With regard to the prior for
the added overdispersion parameter $\psi$ for the Negative Binomial
likelihood, $Gamma(2,0.2)$, parameterised in terms of shape and rate, was
used.

Table 2: Posterior results for SEL data analysis comparing Negative Binomial
GLM-MM models (white columns), CAR-MM (light grey columns) and ICAR-MM (dark
grey columns). $\beta_{1}$ is the coefficient, associated to the share of
South Asian population, and $\beta_{2}$ for the IMD (both min-max normalised)

$\alpha$ $\tau$ $\gamma$ $\beta_{1}$ $\beta_{2}$ $\psi$ mean - 0.52 - - 18.42
26.02 -0.75 -0.74 -0.69 0.75 0.69 0.54 1.02 1.02 0.99 31.61 33.85 34.71 sd -
0.28 - - 7.89 9.01 0.05 0.06 0.07 0.14 0.17 0.20 0.07 0.09 0.11 4.05 4.59 4.63
$2.5\%$ - 0.03 - - 7.13 12.26 -0.84 -0.86 -0.83 0.49 0.34 0.13 0.88 0.85 0.78
24.35 25.73 26.35 $50\%$ - 0.54 - - 17.04 24.71 -0.75 -0.74 -0.69 0.75 0.69
0.54 1.02 1.02 0.99 31.37 33.55 34.46 $97.5\%$ - 0.96 - - 37.64 47.06 -0.66
-0.62 -0.56 1.02 1.02 0.94 1.16 1.19 1.19 40.20 43.68 44.58 width CI - 0.93 -
- 30.51 34.80 0.18 0.24 0.28 0.53 0.68 0.80 0.28 0.34 0.41 15.85 17.94 18.23
$\hat{R}$ - 1.00 - - 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
1.00 1.00 1.00 ESS - 3414 - - 14048 13034 10821 7586 5325 12741 9850 6496
12570 14058 5937 14984 27963 28552

Table 2 shows posterior estimates for the three models estimated on the SEL
dataset. In general estimates are very similar across the three models, with
the exceptions of marginal CAR precision $\tau$ and $\beta_{1}$ coefficient,
the former larger and the latter smaller in the ICAR-MM. As expected, given
the fact that no term is directly accounting for residual spatial
autocorrelation, the GLM-MM model estimates slightly smaller credible
intervals for $\beta_{1}$ and $\beta_{2}$ than the two CAR based models.

Model comparisons in Table 3 do not seem to show substantial differences
between model fits: the ICAR shows the highest $\widehat{elpd}_{loo}$ although
with no significant difference from the other models. Despite the GLM-MM
having the lowest $\overline{RPS}$ and $\overline{DSS}$, the small differences
and the result from other model comparison criteria, there is not sufficient
evidence to claim it provides a better fit. In fact, in Figure 10
(Supplemental material), we can see no discernible difference between
posterior calibration of marginal PPP’s. Therefore, we focus on the CAR-MM
model as it is the one on which we focussed our attention in this paper.

Table 3: Comparison of model check measures. ICAR-MM model is used as reference (model A) when computing differences $\text{diff}(\widehat{elpd}_{loo}^{A,B})$ | $\widehat{elpd}$ | $se(\widehat{elpd}_{loo})$ | $\text{diff}(\widehat{elpd}_{loo}^{A,B})$ | $se\left[\text{diff}(\widehat{elpd}_{loo}^{A,B})\right]$ | $\overline{RPS}$ | $\overline{DSS}$
---|---|---|---|---|---|---
ICAR-MM | -847.8 | 14.2 | 0 | 0 | 35.93 | 9.39
CAR-MM | -849.0 | 13.8 | -1.2 | 1.1 | 35.11 | 9.32
GLM-MM | -849.1 | 13.0 | -1.3 | 2.5 | 34.86 | 9.28
Figure 1: (a) Posterior means for each areal relative risk ($\bm{\rho}$ in
equation 7.2) for the CAR-MM model (b) Probabilities $P(\rho_{i}>1|\bm{y})$
that the relative risk in each exceeds 1 among posterior samples under the
CAR-MM prior

(a)

(b)

###  7.3 Strategic implications

Comparison of posterior relative risk estimates between the 152 MSOAs in the
catchment area shows a wide contrast in estimated prevalence and hence in the
need for diabetes care. Such contrasts show substantial ecological variation –
that is, variation due to area context. Such variations support locality based
or place-based interventions and allocation of resources for tackling diabetes
[39].

Considering the CAR-MM model, the ratio of the $90^{th}$ percentile relative
risk to the $10^{th}$ percentile relative risk has posterior mean ($95\%$
Credible Interval) of $2.31$ $(2.06,2.55)$, an over 2-fold variation in risk.

Figure 1b maps out the probabilities $P(\rho_{i}>1|\tilde{\bm{y}})$ that the
relative risk exceeds 1 under the CAR-MM in the posterior samples [44]. There
are $24$ areas with over $95\%$ probability of elevated risk, indicating
highly elevated diabetes prevalence, and $14$ areas with this probability
exceeding $99\%$. By contrast, a considerably larger minority of $61$ areas
have under $5\%$ probability of elevated risk, with $49$ areas having
exceedance probabilities below $0.01$, indicating significantly low
prevalence.

Figures 1a and 1b also show the spatial features of elevated diabetes levels
in a minority of areas, particularly a spatial concentration in the north of
the region, reinforcing findings such as [24]. The posterior mean ($95\%$ CI)
of the Moran’s I for the relative risks is $0.44$ $(0.35,0.49)$. The
clustering of high risk reflects residential concentration in observed area
risk factors, namely deprived communities and in the South Asian community, as
well as spatial clustering in unobserved influences on diabetes risk. Hence
both area risk factors should figure in any area rankings for locality based
intervention.

Table 4 disaggregates the 152 MSOAs into quintiles according to their
deprivation scores and their percentages South Asian. If we consider the mean
of posterior relative risks in the 30 MSOAs which are most deprived we see an
average ($95\%$ CI) of $1.32$ $(1.15,1.55)$, as compared to a much lower
relative risk in the least deprived MSOAs. The posterior relative risk in the
30 MSOAs with the relatively largest Asian communities has mean ($95\%$ CI)
$1.07$ $(0.64,1.68)$ as compared to significantly lower relative risks in
MSOAs with relatively small Asian communities. Such findings add to existing
UK-based findings that area deprivation and ethnic mix significantly affects
area diabetes prevalence and incidence variations [37, 40].

| IMD | South Asian population
---|---|---
| mean | $2.5\%$ | $97.5\%$ | mean | $2.5\%$ | $97.5\%$
$1^{st}$ | 0.62 | 0.54 | 0.79 | 0.92 | 0.54 | 1.39
$2^{nd}$ | 0.74 | 0.62 | 0.89 | 0.91 | 0.56 | 1.38
$3^{rd}$ | 0.90 | 0.75 | 1.10 | 0.87 | 0.55 | 1.45
$4^{th}$ | 1.13 | 0.96 | 1.64 | 0.93 | 0.58 | 1.45
$5^{th}$ | 1.32 | 1.15 | 1.55 | 1.07 | 0.64 | 1.68
Table 4: Relative risk profiles of posterior means by covariate quintiles

Finally, so far we analysed the areal relative risks $\bm{\rho}$ however, as
mentioned in Section 4, the membership relative risks $\tilde{\bm{\rho}}$ are
also identifiable in the posterior (since $m>n$ in our case) and can be of
inferential interest in many applications. In fact, in our case each
$\tilde{\rho_{j}}$ corresponds to the relative risk of a single GP practice
and can be of significant epidemiological and public health interest.

Therefore, we can replicate the analysis of the probabilities of a relative
risk exceeding 1 as we did above. In the case of GP practices, we can see
that: $41$ practices have over $95\%$ probability, $29$ with over $99\%$, and
$68$ are under $5\%$, with $62$ under $1\%$. In Figure 1b we also show the
locations of those practices on the map.

As expected, in the SEL region, these results are similar to the ones observed
for areal relative risks, given that we have a similar number of practices and
areas. However, these membership level relative risk can give insight in which
GP’s could be facing significant public health challenges. In fact, although
most of the practices with either high or low probability are in accordance
located in areas with similar risk profiles, there are some GP’s for which the
opposite is true, i.e. high or low risk practices in respectively low or high
risk areas. A purely areal analysis of prevalence would have not been able to
identify these practices.

##  8 Discussion

In this work we have studied and evaluated the theoretical underpinnings of
the CAR-MM model devising a general framework to evaluate its structure and
implications to spatial and relative risk modelling. We compared our approach
to that of [42], provided a more general framework able to encompass a
potentially wide range of applications, and showed how this could apply to the
parameterisation and evaluation of a CAR-MM model. Building on that, we
considered identifiability in both the likelihood and the posterior for
different membership sizes compared to areas and both parameterisations.

The simulation study highlights how this approach leads to well calibrated
posterior samples and small bias for estimated parameters. We showed how the
two parameterisations, although theoretically identical in the $m=n$ case, can
lead to significantly different results in terms of both computational time
and identifiability when $m\neq n$. More specifically, we recommend the use of
the post parameterisation in all cases although we warn of the possible
identifiability problems in the $m<n$ case. In this regard, the present paper
supports the previous analyses and findings of both [23] and [42], where $m>n$
datasets were considered, and where the models were post parameterised.

Finally, the data analysis allowed us to show how this approach can be used in
a real world application where the MM principle and CAR frameworks merge
naturally, thus allowing for a direct implementation of the CAR-MM prior.
Datasets from different institutions, in this case the National Health Service
(NHS) and the Census, can be directly integrated leading to a modelling
strategy that not only is able to pinpoint areas with increased relative
risks, but also specific practices, by making use of the membership level
modelling. These evaluations can be of interest for both epidemiological
research and public health considerations.

## References

  * [1] Sudipto Banerjee, Bradley P. Carlin and Alan E. Gelfand  “Hierarchical modeling and analysis for spatial data”  ChapmanHall/CRC, 2014, pp. 562  URL: https://www.crcpress.com/Hierarchical-Modeling-and-Analysis-for-Spatial-Data/Banerjee-Carlin-Gelfand/p/book/9781439819173
  * [2] Julian Besag  “Spatial Interaction and the Statistical Analysis of Lattice Systems”  In _Journal of the Royal Statistical Society. Series B (Methodological)_ 36 WileyRoyal Statistical Society, 1974, pp. 192–236  URL: https://www.jstor.org/stable/2984812
  * [3] Julian Besag, Jeremy York and Annie Mollié  “Bayesian image restoration, with two applications in spatial statistics”  In _Annals of the Institute of Statistical Mathematics_ 43.1 Kluwer Academic Publishers, 1991, pp. 1–20
  * [4] Michael Betancourt  “A Conceptual Introduction to Hamiltonian Monte Carlo”, 2018  arXiv:1701.02434v2
  * [5] Marta Blangiardo and Michela Cameletti  “Spatial and Spatio-temporal Bayesian Models with R - INLA”  In _Spatial and Spatio-temporal Bayesian Models with R - INLA_ John Wiley & Sons, 2015, pp. 1–308
  * [6] Jonathan R. Bradley, Christopher K. Wikle and Scott H. Holan  “Bayesian Spatial Change of Support for Count-Valued Survey Data With Application to the American Community Survey”  In _Journal of the American Statistical Association_ 111.514 American Statistical Association, 2016, pp. 472–487  URL: https://www.tandfonline.com/doi/abs/10.1080/01621459.2015.1117471
  * [7] William J Browne, Harvey Goldstein and Jon Rasbash  “Multiple membership multiple classification (MMMC) models”  In _Statistical Modelling: An International Journal_ 1.2 Sage PublicationsSage CA: Thousand Oaks, CA, 2001, pp. 103–124  URL: http://journals.sagepub.com/doi/10.1177/1471082X0100100202
  * [8] Bob Carpenter et al.  “Stan : A Probabilistic Programming Language”  In _Journal of Statistical Software_ 76.1, 2017, pp. 1–32  DOI: 10.18637/jss.v076.i01
  * [9] E.. Catchpole and B..T. Morgan  “Detecting parameter redundancy”  In _Biometrika_ 84.1 Oxford University Press, 1997, pp. 187–196  DOI: 10.1093/biomet/84.1.187
  * [10] D Cole  “Parameter Redundancy and Identifiability”, Chapman $\backslash$& Hall/CRC Interdisciplinary Statistics  CRC Press, 2020  URL: https://books.google.co.uk/books?id=JrnjDwAAQBAJ
  * [11] Sylvain Coly, Myriam Garrido, David Abrial and Anne-Françoise Yao  “Bayesian hierarchical models for disease mapping applied to contagious pathologies”  In _bioRxiv_ Cold Spring Harbor Laboratory, 2019, pp. 766071  DOI: 10.1101/766071
  * [12] Peter D. Congdon  “Bayesian Hierarchical Models”  In _Bayesian Hierarchical Models_ ChapmanHall/CRC, 2019  DOI: 10.1201/9780429113352
  * [13] Samantha R Cook, Andrew Gelman and Donald B Rubin  “Validation of Software for Bayesian Models Using Posterior Quantiles”  In _Journal of Computational and Graphical Statistics_ 15.3 Taylor & Francis, 2006, pp. 675–692  DOI: 10.1198/106186006X136976
  * [14] Noel Cressie and Gardar Johannesson  “Fixed rank kriging for very large spatial data sets”  In _J. R. Statist. Soc. B_ 70, 2008, pp. 209–226
  * [15] Claudia Czado, Tilmann Gneiting and Leonhard Held  “Predictive model assessment for count data”  In _Biometrics_ 65.4 Blackwell Publishing Inc., 2009, pp. 1254–1261  DOI: 10.1111/j.1541-0420.2009.01191.x
  * [16] A P Dawid  “Conditional Independence in Statistical Theory”  In _J. R. Statist. Soc. B_ 41.1, 1979, pp. 1–31
  * [17] A. Dawid and Paola Sebastiani  “Coherent dispersion criteria for optimal experimental design”  In _The Annals of Statistics_ 27.1 Institute of Mathematical Statistics, 1999, pp. 65–81  DOI: 10.1214/AOS/1018031101
  * [18] Lynn E. Eberly and Bradley P. Carlin  “Identifiability and convergence issues for Markov chain Monte Carlo fitting of spatial models”  In _Statistics in Medicine_ 19.17-18, 2000, pp. 2279–2294  DOI: 10.1002/1097-0258(20000915/30)19:17/18¡2279::AID-SIM569¿3.0.CO;2-R
  * [19] A E Gelfand, P Diggle, P Guttorp and M Fuentes  “Handbook of Spatial Statistics”, Chapman & Hall/CRC Handbooks of Modern Statistical Methods  CRC Press, 2010  URL: https://books.google.co.uk/books?id=Xf4leslPDzsC
  * [20] Alan E. Gelfand and Sujit K. Sahu  “Identifiability, Improper Priors, and Gibbs Sampling for Generalized Linear Models”  In _Journal of the American Statistical Association_ 94.445, 1999, pp. 247–253  DOI: 10.1080/01621459.1999.10473840
  * [21] Andrew Gelman et al.  “Bayesian data analysis”  CRC Press, 2013, pp. 667  URL: https://books.google.co.uk/books/about/Bayesian%7B%5C_%7DData%7B%5C_%7DAnalysis%7B%5C_%7DThird%7B%5C_%7DEdition.html?id=ZXL6AQAAQBAJ%7B%5C&%7Dredir%7B%5C_%7Desc=y
  * [22] T. Goicoa, A. Adin, M.. Ugarte and J.. Hodges  “In spatio-temporal disease mapping models, identifiability constraints affect PQL and INLA results”  In _Stochastic Environmental Research and Risk Assessment_ 32.3 Springer New York LLC, 2018, pp. 749–770  DOI: 10.1007/S00477-017-1405-0
  * [23] Marco Gramatica, Peter Congdon and Silvia Liverani  “Bayesian modelling for spatially misaligned health areal data: A multiple membership approach”  In _Journal of the Royal Statistical Society. Series C: Applied Statistics_ 70.3 Blackwell Publishing Ltd, 2021, pp. 645–666  DOI: 10.1111/rssc.12480
  * [24] Chris Green, Robert D. Hoppa, T. Young and J.. Blanchard  “Geographic analysis of diabetes prevalence in an urban area”  In _Social Science & Medicine_ 57.3 Pergamon, 2003, pp. 551–560  DOI: 10.1016/S0277-9536(02)00380-5
  * [25] Lesley M.. Hall et al.  “Fat Oxidation, Fitness and Skeletal Muscle Expression of Oxidative/Lipid Metabolism Genes in South Asians: Implications for Insulin Resistance?”  In _PLoS ONE_ 5.12 Public Library of Science, 2010, pp. e14197  DOI: 10.1371/journal.pone.0014197
  * [26] D A Harville  “Matrix Algebra From a Statistician’s Perspective”  Springer, 2008  URL: https://books.google.it/books?id=kZGBQijgGV8C
  * [27] Health and Social Care Information Centre  “Quality and Outcomes Framework (QOF) - 2015-16”, 2016  URL: https://digital.nhs.uk/data-and-information/publications/statistical/quality-and-outcomes-framework-achievement-prevalence-and-exceptions-data/quality-and-outcomes-framework-qof-2015-16
  * [28] Health and Social Care Information Centre  “Patients Registered at a GP Practice - January 2019”, 2019  URL: https://digital.nhs.uk/data-and-information/publications/statistical/patients-registered-at-a-gp-practice/january-2019
  * [29] Xiaoping Jin, Bradley P. Carlin and Sudipto Banerjee  “Generalized hierarchical multivariate CAR models for areal data”  In _Biometrics_ 61.4, 2005, pp. 950–961  DOI: 10.1111/j.1541-0420.2005.00359.x
  * [30] Max Joseph  “Exact sparse CAR models in Stan”  In _Stan Case studies_ , 2016  URL: https://mc-stan.org/users/documentation/case-studies/mbjoseph-CARStan.html
  * [31] Thomas King, Simon Butcher and Lukasz Zalewski  “Apocrita - High Performance Computing Cluster for Queen Mary University of London”, 2017  DOI: 10.5281/zenodo.438045
  * [32] A.B. Lawson  “Bayesian Disease Mapping: Hierarchical Modeling in Spatial Epidemiology, Third Edition”  CRC Press, 2018, pp. 652
  * [33] Ying C MacNab  “On Gaussian Markov random fields and Bayesian disease mapping:”  In _https://doi.org/10.1177/0962280210371561_ 20.1 SAGE PublicationsSage UK: London, England, 2010, pp. 49–68  DOI: 10.1177/0962280210371561
  * [34] E.. Marshall and D.. Spiegelhalter  “Approximate cross-validatory predictive checks in disease mapping models”  In _Statistics in Medicine_ 22.10 Stat Med, 2003, pp. 1649–1660  DOI: 10.1002/sim.1403
  * [35] Miguel A. Martinez-Beneito and Paloma Botella-Rocamora  “Disease Mapping : From Foundations to Multidimensional Modeling”  In _Disease Mapping_ ChapmanHall/CRC, 2019  DOI: 10.1201/9781315118741
  * [36] P McCullagh and J.A. Nedler  “Generalized Linear Models”  CRC Press, 1983  URL: https://books.google.co.uk/books?id=SPSDDwAAQBAJ
  * [37] Alison Moody, Giovanna Cowley, Linda Ng Fat and Jennifer S. Mindell  “Social inequalities in prevalence of diagnosed and undiagnosed diabetes and impaired glucose regulation in participants in the Health Surveys for England series”  In _BMJ Open_ 6.2 British Medical Journal Publishing Group, 2016, pp. e010155  DOI: 10.1136/BMJOPEN-2015-010155
  * [38] Mitzi Morris et al.  “Bayesian hierarchical spatial models: Implementing the Besag York Mollié model in stan”  In _Spatial and Spatio-temporal Epidemiology_ 31 Elsevier Ltd, 2019, pp. 100301  DOI: 10.1016/j.sste.2019.100301
  * [39] National Institute for Health and Care Excellence  “Preventing Type 2 Diabetes: Population and Community-level Interventions”, 2014
  * [40] Yoshitaka Nishino, Stuart Gilmour and Kenji Shibuya  “Inequality in Diabetes-Related Hospital Admissions in England by Socioeconomic Deprivation and Ethnicity: Facility-Based Cross-Sectional Analysis”  In _PLOS ONE_ 10.2 Public Library of Science, 2015, pp. e0116689  DOI: 10.1371/JOURNAL.PONE.0116689
  * [41] Northern Ireland Statistics and Research Agency Office for National Statistics  “2011 Census aggregate data”  UK Data Service, 2016  DOI: http://dx.doi.org/10.5257/census/aggregate-2011-1
  * [42] Oana Petrof et al.  “On the impact of residential history in the spatial analysis of diseases with a long latency period: A study of mesothelioma in Belgium”  In _Statistics in Medicine_ 39.26 John WileySons Ltd, 2020, pp. 3840–3866  DOI: 10.1002/sim.8697
  * [43] “Linear Statistical Inference and its Applications”, Wiley Series in Probability and Statistics  Hoboken, NJ, USA: John Wiley & Sons, Inc., 1973  DOI: 10.1002/9780470316436
  * [44] Sylvia Richardson, Andrew Thomson, Nicky Best and Paul Elliott  “Interpreting Posterior Relative Risk Estimates in Disease-Mapping Studies”  In _Environmental Health Perspectives_ 112.9 National Institute of Environmental Health Sciences, 2004, pp. 1016  DOI: 10.1289/EHP.6740
  * [45] Håvard Rue and Leonhard Held  “Gaussian Markov random fields : theory and applications”  Chapman & Hall/CRC, 2005, pp. 263  URL: https://www.crcpress.com/Gaussian-Markov-Random-Fields-Theory-and-Applications/Rue-Held/p/book/9781584884323
  * [46] Sean Talts et al.  “Validating Bayesian Inference Algorithms with Simulation-Based Calibration”, 2018  arXiv: https://arxiv.org/abs/1804.06788v2
  * [47] Aki Vehtari, Andrew Gelman and Jonah Gabry  “Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC”  In _Statistics and Computing_ 27.5 Springer US, 2017, pp. 1413–1432  DOI: 10.1007/s11222-016-9696-4
  * [48] Aki Vehtari et al.  “loo: Efficient leave-one-out cross-validation and WAIC for Bayesian models” R package version 2.4.1, 2020  URL: https://mc-stan.org/loo/
  * [49] Aki Vehtari et al.  “Rank-Normalization, Folding, and Localization: An Improved Rˆ for Assessing Convergence of MCMC”  In _Bayesian Analysis_ 16.2 Institute of Mathematical Statistics, 2021, pp. 667–719  DOI: 10.1214/20-ba1221
  * [50] Jay M. Ver Hoef, Ephraim M. Hanks and Mevin B. Hooten  “On the relationship between conditional (CAR) and simultaneous (SAR) autoregressive models”  In _Spatial Statistics_ 25 Elsevier B.V., 2018, pp. 68–85  DOI: 10.1016/j.spasta.2018.04.006


Charmonia photo-production in ultra-peripheral and peripheral PbPb collisions
with LHCb

Weisong Duan on behalf of the LHCb collaboration

E-Mail<EMAIL_ADDRESS>

South China Normal University, Guangzhou, China

Presented at the Low-$x$ Workshop, Elba Island, Italy, September 27–October 1
2021

###### Abstract

The LHCb recorded ∼ 210 $\mu b^{-1}$ integrated luminosity of PbPb collisions
at $\sqrt{s_{NN}}$ = 5.02 TeV in 2018. With an increase of the luminosity by a
factor of 20 compared to the previous 2015 PbPb dataset, precise measurements
on photo-produced charmonia in ultra-peripheral collisions are now possible.
Moreover, the great momentum resolution of the detector allows photo-produced
$J/\psi$ in collisions with a nuclear overlap to be studied. This new type of
probe is sensitive to the geometry of the collisions but also to the
electromagnetic field of the Pb nuclei. In this contribution, we present the
latest results on $J/\psi$ photo-production measured by LHCb in peripheral and
ultra-peripheral PbPb collisions.

##  1 Introduction

The LHCb detector is a single-arm forward spectrometer fully instrumented in
the pseudorapidity range 2 < $\eta$ < 5 Collaboration_2008 . It has a high
precision tracking system, which provides excellent vertex and momentum
resolution, and full particle identification. Compared to the ALICE, CMS and
ATLAS, LHCb covers the forward rapidity region, providing better access to the
gluon distribution at small $x$.

The photonuclear production of vector mesons such as $J/\psi$ is sensitive to
the gluon parton distribution function in the nucleus at small Bjorken-$x$,
which is estimated by $x\approx(m_{J/\psi}\cdot e^{-y})/\sqrt{s_{NN}}$, where
$m_{J/\psi}$ is mass of $J/\psi$ and $y$ is its rapidity. Coherent
photoproduction of $J/\psi$ meson provides a way to study the nuclear
shadowing effects at small Bjorken-$x$ ranging from $10^{-5}$ to $10^{-2}$ at
LHC energies.

##  2 Study of coherent $J/\psi$ production in ultra-peripheral lead-lead
collisions at $\sqrt{s_{NN}}$ = 5 TeV

The ultra-peripheral collisions, UPCs, are $\rm{PbPb\to Pb+Pb+X}$ in which two
ions interact via their cloud of virtual photons. If the photon couples
coherently to the nucleus as a whole, it is called coherent production. If the
photon couples with one nucleon leading to the breakup of the target nucleus,
it is called incoherent production.

In UPCs, coherent $J/\psi$ meson production can be described by the
interaction between photons and gluons, according to the Regge theory
PhysRevC.86.014905 ; PhysRevC.93.055206 , gluons are considered as a single
object with vacuum quantum numbers, which is called pomeron ($\enspace\rm
P\kern-10.00002ptI$ ). The cross-section for photoproduction gives constraints
on the gluon parton distribution functions. This process has low
multiplicities and very low transverse momentum $p_{T}$.

The $J/\psi$ mesons are reconstructed through the $J/\psi\to\mu^{+}\mu^{-}$
decay channel, using 2015 Pb-Pb data samples, corresponding to an integrated
luminosity of 10 $\mu b^{-1}$. An ultra-peripheral electromagnetic interaction
could occur simultaneously with the hadron collision. The HeRSCheL detector
Akiba_2018  is used to reject backgrounds from hadronic interactions.

The number of candidates are obtained by fitting the di-muon spectrum as shown
in Fig. 1. The $J/\psi$ and $\psi(2S)$ mass spectrum are modeled by a double-
sided Crystal ball function, and the non-resonant background is modeled by
exponential function multiplied by an first-order polynomial function. Thus we
determine the number of $J/\psi$ candidates within the $J/\psi$ mass window
3040 MeV/$c^{2}$ $\sim$ 3165 MeV/$c^{2}$ and the number of $\psi(2S)$
candidates within the $\psi(2S)$ mass window 3608 MeV/$c^{2}$ $\sim$ 3763
MeV/$c^{2}$.

Figure 1: The di-muon invariant mass distribution in the rapidity range
between 2.0 and 4.5. The solid blue line corresponds to the $J/\psi$ meson.
The solid green line represents the $\psi(2S)$ meson. The dashed black line
represents non-resonant background.

To determine the coherent $J/\psi$ production, a fit to the ${\rm
log}(p^{2}_{T})$ is performed to extract the coherent $J/\psi$ mesons within
the $J/\psi$ mass window. The ${\rm log}(p^{2}_{T})$ distribution of the
$J/\psi$ mesons is shown in Fig. 2. In short, the number of inclusive $J/\psi$
mesons is obtained by the invariant mass fit, and the number of coherent
$J/\psi$ mesons is obtained by the ${\rm log}(p^{2}_{T})$ fit.

Figure 2: The ${\rm log}(p^{2}_{T})$ distribution of $J/\psi$ candidates in
the interval 2.5 < $y$ < 3.0, where the unit of $p_{T}$ is GeV/c. The solid
blue line represents coherent $J/\psi$ distribution. The dashed red line
represents the incoherent distribution. The dashed green line represents feed
down from $\psi(2S)$. The dashed black line represents non-resonant
background.

The results of the differential cross section are calculated in five rapidity
bins as shown in Fig. 3. We compared the results between experimental results
and theoretical predictions PhysRevC.93.055206 ; PhysRevC.84.011902 ;
PhysRevC.97.024901 ; PhysRevD.96.094027 ; MANTYSAARI2017832 . The coherent
$J/\psi$ cross-section production is given by:

$\frac{\mathrm{d}\sigma_{coh,J/\psi}}{\mathrm{d}y}=\frac{N_{coh,J/\psi}}{\varepsilon_{t}\cdot\mathcal{L}\cdot\Delta
y\cdot\mathcal{B}(J/\psi\to\mu^{+}\mu^{-})}$ (1)

where $\varepsilon_{t}$ is the total efficiency, $\mathcal{L}$ is an
integrated luminosity of the Pb-Pb data sample, and
$\mathcal{B}=(5.961±0.033)\%$ is the $J/\psi\to\mu^{+}\mu^{-}$ branching
ratio.

Figure 3: Differential cross-section of coherent $J/\psi$ production as a
function of rapidity, with comparisons to phenomenological models.

Similarly, ALICE measured the coherent $J/\psi$ production cross-section
2019134926 , so we also compared the coherent $J/\psi$ production cross-
section between ALICE and LHCb, as shown in Fig. 4. The LHCb result is
slightly lower than the ALICE measurement by around 1.3 $\sigma$. Measurements
of the coherent $J/\psi$ and $\psi(2S)$ are currently underway using 2018 Pb-
Pb data, corresponding to an integrated luminosity of 210 $\mu b^{-1}$. Fig. 5
shows the di-muon invariant mass distribution in the range between 2.7 and 4.0
GeV. The final results are expected in the near future.

Figure 4: Differential cross-section of coherent $J/\psi$ production as a
function of rapidity, with comparison to ALICE measurements. Figure 5: The
invariant mass distribution of $J/\psi$ and $\psi(2S)$ candidates. The solid
green line represents $J/\psi$ signal, the solid red line represents
$\psi(2S)$ signal, and the yellow regions represents non-resonant background.

##  3 Study of $J/\psi$ photo-production in lead-lead peripheral collisions at
$\sqrt{s_{NN}}$ = 5 TeV

The second result in this contribution is photo-production of $J/\psi$ at low
$p_{T}$, studied in peripheral Pb-Pb collisions at $\sqrt{s_{NN}}$ = 5 TeV,
using the data sample collected by LHCb in 2018, with an integrated luminosity
of about 210 $\mu b^{-1}$ lhcbcollaboration2021study .

The $J/\psi$ candidates are selected through the $J/\psi\to\mu^{+}\mu^{-}$
decay channel. The di-muon invariant mass spectrum of the selected candidates
in the range between 3.0 and 3.2 GeV is shown in the left of Fig. 6, for the
$J/\psi$ mesons with $\rm{p_{T}}$ <15.0 GeV/c and the number of participants
$\left<N_{part}\right>$ = 10.6 $\pm$ 2.9 in full LHCb rapidity coverage 2 <
$y$ <4.5.

The inclusive $J/\psi$ candidates consists of photo-produced and hadronically
produced $J/\psi$ mesons, which are separated by an unbinned maximum
likelihood fit to the ${\rm log}(p^{2}_{T})$ distribution, as shown in the
right of Fig. 6. In this figure, the transverse momentum of photo-produced
$J/\psi$ yields (red dotted line) are visible in the range between 0 and 250
MeV/c.

Fig. 7 shows the photo-produced $J/\psi$ meson yields as a function of
$p_{T}$(right), and $\left<N_{part}\right>$ (left). The mean $p_{T}$ of the
coherent $J/\psi$ is estimated to be $\left<p_{T}\right>$ = 64.9 $\pm$ 2.4
MeV/c. Theoretical predictions PhysRevC.97.044910 ; PhysRevC.99.061901  are
drawn in open circles, and are qualitatively in agreement with the
experimental results in shape.

Figure 6: The left plot is the invariant mass distribution of $J/\psi$
candidates in $p_{T}$ < 15.0 GeV/c and 2 < $y$ <4.5, with
$\left<N_{part}\right>$ = 10.6 $\pm$ 2.9. The right plot is the ${\rm
ln}(p^{2}_{T})$ distribution of $J/\psi$ candidates after background
subtraction for the same kinematic interval.

Figure 7: Illustration of invariant yields of photo-produced $J/\psi$ mesons
as a function of $N_{part}$ and $p_{T}$. Note that the blue error bars
represent the statistical uncertainty and the red error bars the total
uncertainty. Theoretical model predictions PhysRevC.97.044910 ;
PhysRevC.99.061901  are shown in open circles in black and green.

##  4 Conclusions

LHCb detector’s unique geometry acceptance allows us to study the nuclear
shadowing in small $x$ region through the UPC collisions. The $J/\psi$
photoproduction in UPC using PbPb collisions at 5TeV is measured by LHCb, and
is compared to results by ALICE. The photoproduced $J/\psi$ production in
peripheral PbPb collisions is also measured with high precision at very low
$p_{T}$, and is compared to theoretical calculations. These results
demonstrate the capabilities of the LHCb detector in studying nuclear effects.
More results from the large 2018 PbPb dataset are expected in the future.

## References

  * (1) LHCb Collaboration, A. A. Alves Jr. et al.  The LHCb detector at the LHC.  JINST 3(08):S08005–S08005, aug 2008.
  * (2) S. P. Jones et al.  Exclusive $J/\psi$ and $\Upsilon$ photoproduction and the low $x$ gluon.  Journal of Physics G: Nuclear and Particle Physics, 43(3):035002, Jan 2016.
  * (3) A. Cisek, W. Schäfer, and A. Szczurek.  Exclusive coherent production of heavy vector mesons in nucleus-nucleus collisions at energies available at the cern large hadron collider.  Phys. Rev. C, 86:014905, Jul 2012.
  * (4) V. Guzey, E. Kryshen, and M. Zhalov.  Coherent photoproduction of vector mesons in ultraperipheral heavy ion collisions: Update for run 2 at the cern large hadron collider.  Phys. Rev. C, 93:055206, May 2016.
  * (5) K. Carvalho Akiba et.al.  The HeRSCheL detector: high-rapidity shower counters for LHCb.  JINST, 13(04):P04017–P04017, apr 2018.
  * (6) V. P. Gonçalves and M. V. T. Machado.  Vector meson production in coherent hadronic interactions: Update on predictions for energies available at the bnl relativistic heavy ion collider and the cern large hadron collider.  Phys. Rev. C, 84:011902, Jul 2011.
  * (7) J. Cepila, J. G. Contreras, and M. Krelina.  Coherent and incoherent $J/\psi$ photonuclear production in an energy-dependent hot-spot model.  Phys. Rev. C, 97:024901, Feb 2018.
  * (8) V. P. Gonçalves, M. V. T. Machado, B. D. Moreira, F. S. Navarra, and G. Sampaio dos Santos.  Color dipole predictions for the exclusive vector meson photoproduction in $pp$, $p\mathrm{Pb}$, and pbpb collisions at run 2 lhc energies.  Phys. Rev. D, 96:094027, Nov 2017.
  * (9) Heikki Mäntysaari and Björn Schenke.  Probing subnucleon scale fluctuations in ultraperipheral heavy ion collisions.  Physics Letters B, 772:832–838, 2017.
  * (10) S. Acharya et al. ALICE collaboration.  Coherent $J/\psi$ photoproduction at forward rapidity in ultra-peripheral pb–pb collisions at $\sqrt{s_{NN}}$=5.02 tev.  Physics Letters B, 798:134926, 2019.
  * (11) LHCb collaboration and R. Aaij et.al.  Study of $J/\psi$ photo-production in lead-lead peripheral collisions at $\sqrt{s_{\rm NN}}$ = 5 TeV, arXiv: 2108.02681, 2021.
  * (12) W. Zha, S. R. Klein, R. Ma, L. Ruan, T. Todoroki, Z. Tang, Z. Xu, C. Yang, Q. Yang, and S. Yang.  Coherent $J/\psi$ photoproduction in hadronic heavy-ion collisions.  Phys. Rev. C, 97:044910, Apr 2018.
  * (13) W. Zha, L. Ruan, Z. Tang, Z. Xu, and S. Yang.  Double-slit experiment at fermi scale: Coherent photoproduction in heavy-ion collisions.  Phys. Rev. C, 99:061901, Jun 2019.