Modelling collinear and spatially correlated data
Silvia Liverani
Department of Mathematics
Brunel University London
UK
Medical Research Centre Biostatistics Unit
Cambridge
UK
Department of Epidemiology and Biostatistics
Imperial College London
UK
Aurore Lavigne
University of Lille 3
France
Marta Blangiardo
MRC-PHE Centre for Environment and Health
Department of Epidemiology and Biostatistics
Imperial College London
UK
In this work we present a statistical approach to distinguish and interpret the complex relationship between several predictors and a response variable at the small area level, in the presence of i) high correlation between the predictors and ii) spatial correlation for the response. Covariates which are highly correlated create collinearity problems when used in a standard multiple regression model. Many methods have been proposed in the literature to address this issue. A very common approach is to create an index which aggregates all the highly correlated variables of interest. For example, it is well known that there is a relationship between social deprivation measured through the Multiple Deprivation Index (IMD) and air pollution; this index is then used as a confounder in assessing the effect of air pollution on health outcomes (e.g. respiratory hospital admissions or mortality). However it would be more informative to look specifically at each domain of the IMD and at its relationship with air pollution to better understand its role as a confounder in the epidemiological analyses. In this paper we illustrate how the complex relationships between the domains of IMD and air pollution can be deconstructed and analysed using profile regression, a Bayesian non-parametric model for clustering responses and covariates simultaneously. Moreover, we include an intrinsic spatial conditional autoregressive (ICAR) term to account for the spatial correlation of the response variable. \vspace{0.3cm {Keywords: profile regression, Bayesian clustering, spatial modelling, collinearity, index of multiple deprivation, pollution.
SSTEarxiv_v1.tex
\documentclass[11pt]{article}
\usepackage{fullpage}

\RequirePackage{amsthm,amsmath}
\RequirePackage{natbib}
\RequirePackage[colorlinks,citecolor=blue,urlcolor=blue]{hyperref}
\usepackage{subfigure}
% settings
%\pubyear{2005}
%\volume{0}
%\issue{0}
%\firstpage{1}
%\lastpage{8}
%\arxiv{arXiv:0000.0000}

%\startlocaldefs
%\numberwithin{equation}{section}
%\theoremstyle{plain}
%\newtheorem{thm}{Theorem}[section]
%\endlocaldefs

\usepackage{fullpage}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm,graphicx,bm}

\usepackage{natbib}

\newcommand{\todo}[1]{\textbf{\hl{To Do}}\hl{ - #1}}
\newcommand{\question}[1]{\textbf{\hl{Qn}}\hl{ - #1}}

\newcommand{\bTheta}{\boldsymbol{\Theta}}
\newcommand{\bV}{\boldsymbol{V}}
\newcommand{\bU}{\boldsymbol{U}}
\newcommand{\bD}{\boldsymbol{D}}
\newcommand{\bX}{\boldsymbol{X}}
\newcommand{\bY}{\boldsymbol{Y}}
\newcommand{\bW}{\boldsymbol{W}}
\newcommand{\bZ}{\boldsymbol{Z}}
\newcommand{\tbTheta}{\tilde{\boldsymbol{\Theta}}}
\newcommand{\tbV}{\tilde{\boldsymbol{V}}}
\newcommand{\tbZ}{\tilde{\boldsymbol{Z}}}

\newtheorem{proposition}{Proposition}

\begin{document}
\title{Modelling collinear and spatially correlated data}
%\author{}
\author{Silvia Liverani\\ Department of Mathematics, Brunel University London, UK\\ Medical Research Centre Biostatistics Unit, Cambridge, UK\\ Department of Epidemiology and Biostatistics, Imperial College London, UK  \and Aurore Lavigne\\ University of Lille 3, France \and Marta Blangiardo \\MRC-PHE Centre for Environment and Health, \\Department of Epidemiology and Biostatistics, Imperial College London, UK}
\date{}
\maketitle

\begin{abstract}
In this work we present a statistical approach to distinguish and interpret the complex relationship between several predictors and a response variable at the small area level, in the presence of i) high correlation between the predictors and ii) spatial correlation for the response. 

Covariates which are highly correlated create collinearity problems when used in a standard multiple regression model. Many methods have been proposed in the literature to address this issue. A very common approach is to create an index which aggregates all the highly correlated variables of interest. For example, it is well known that there is a relationship between social deprivation measured through the Multiple Deprivation Index (IMD) and air pollution; this index is then used as a confounder in assessing the effect of air pollution on health outcomes (e.g. respiratory hospital admissions or mortality). However it would be more informative to look specifically at each domain of the IMD and at its relationship with air pollution to better understand its role as a confounder in the epidemiological analyses.

In this paper we illustrate how the complex relationships between the domains of IMD and air pollution can be deconstructed and analysed using profile regression, a Bayesian non-parametric model for clustering responses and covariates simultaneously. Moreover, we include an intrinsic spatial conditional autoregressive (ICAR) term to account for the spatial correlation of the response variable. 


\vspace{0.3cm}

{\bf Keywords}: profile regression, Bayesian clustering, spatial modelling, collinearity, index of multiple deprivation, pollution.
\end{abstract}

\section{Introduction}
In many statistical applications a common challenge arises when trying to assess meaningful relationships between explanatory variables and outcomes through regression models, due to the potential collinearity of the explanatory variables. This issue is well known in epidemiological or social studies, for instance where questionnaires or surveys collect information on a large number of potential risk factors for particular end points; in this context a simplistic approach consists in examining  each variable in turn to avoid the instability in the estimates due to the collinearity, making it impossible to judge the more realistic complex relationship involving several risk factors at the same time. A different approach combines all the relevant variables into summary scores or indexes and assesses the relationship of these with the outcome of interest, which is free from the collinearity issue, but loses information on the single variables included in the summary. 

Recently, Dirichlet process mixture models have been used as an alternative to regression models \citep{Dunson08,Bigelow09}. In this paper we focus on the model known as profile regression and proposed by \cite{Molitor10}. Profile regression is a Bayesian non--parametric method which assesses the link between potentially collinear variables and a response through cluster membership. This allows to formally take into account the correlation between the variables without the need to create a summary score, giving more flexibility to the inferential process. Profile regression has been used on several applications in environmental and social epidemiology and the R package PReMiuM \citep{Liverani15} makes it readily available to any applied researcher. For instance \cite{Molitor10} considered the National Survey of Children's Health and in particular investigated a large number of health and social related variables on mental health of children age 6--17, while \cite{Papathomas11} focussed on profiles of exposure to environmental carcinogens and lung cancer in the EPIC European cohort.  Profile regression has also been used in environmental epidemiology \citep{Pirani15}, for studying risk functions associated with multi-dimensional exposure profiles \citep{Hastie13,Molitor14} as well as for looking for gene-gene interactions \citep{Papathomas12}.%Finally, \cite{Molitor11} used again the National Survey of Children's Health  to study the link between air pollution and poverty in Los Angeles and how this affects mental health of children living in that area.

In its present formulation, profile regression has only been used for studies based on cohorts or surveys where information on the predictors/outcomes is available on each individual; in this paper we extend the method to fit small area studies,  commonly used in epidemiological surveillance (see for instance \citealt{Elliott04}) or in studies where the interest lies on the spatial variability of an outcome \citep{Barcelo09} or on cluster detection \citep{Abellan08,Li12}. In this types of studies information is available at the area level rather than at the individual level and space is used as a proxy for any unmeasured variable; the common assumption is that areas which are close to each other are more similar than those further apart, suggesting that an additional source of correlation, namely \textit{spatial correlation} needs to be accommodated in the models. We incorporate it in the model  through a conditional autoregressive structure \citep{Besag91} based on a neighborhood definition, thus assuming that conditional on the neighborhood structure, two areas are independent from each other if they do not share boundaries. We apply the \textit{spatial profile regression} to the problem of environmental and social inequalities in London, jointly modelling social deprivation and air pollution to highlight the presence of environmental justice.

The paper is structured as follows. In Section \ref{sect:casestudy} we present the motivating example for our methodological development of the spatial profile regression, introducing the context of social and environmental inequalities and how they are related; we also describe the available data. In Section \ref{sect:methods} we provide a brief summary of the profile regression and present how to extend it to include spatial correlation.  In Section \ref{sect:results} we illustrate how the model works on evaluating the relationship between social deprivation and air pollution. Section \ref{sect:conclusion} presents some discussion points and ideas for future work. 
 
\section{Example: social deprivation and air pollution in London}\label{sect:casestudy}

The scientific literature reports mixed evidence on the link between socio-economic status and air pollution. Recent studies indicated that air pollution tends to influence most deprived groups, suggesting that people with lower socio-economic status are more likely to live in a more hazardous and polluted living environment, accidentally or deliberately \citep{Brown95,Oneill03,Blowers94,Morello02}. In particular, ecological studies using small areas such as neighbourhoods, census tracts and post codes, report this association, while studies carried out at a lower spatial resolution (e.g. region, country), thus characterised by more aggregate measurements of socio-economic characteristics, showed either non-existent or negative associations \citep{Davidson00,Laurent07}, presumably due to the large within-area variability not taken into account, or even an inverse association, with higher exposures in less deprived groups \citep{Perlin95}. In the UK several studies reported positive or non-linear correlation between environmental pollution and the deprivation index at both small area level and country level. However the results varied depending on the selection of environmental hazards and scale of analysis (Briggs et al., 2008), calling for some more research on the topic.

Understanding environmental and social inequalities is a key issue as growing health disparities appear between people with socially disadvantaged and privileged social classes, which can translate into increased mortality or morbidity for the low socio-economic groups across a wide range of diseases \citep{Brulle06,Benach01}, including lung cancer \citep{Pope11}, cardiovascular events \citep{Tonne07,Peters04}, and childhood respiratory diseases \citep{Morgenstern07}. 

In addition, the exposure of air pollution can lead to negative health outcomes acutely or chronically \citep{Chen08}. Previous studies reported possible mechanisms to explain how environmental exposures result in greater health impact among socially disadvantaged groups, who may have increased susceptibility to the effect of these exposures because of limited access to health care and psychosocial stress; underlying health conditions such as cardiovascular diseases and respiratory diseases that increase susceptibility to the effect of these exposure may also vary between deprived and privileged populations \citep{Oneill03,Morello06}. These environmental exposure inequalities are increasingly considered as a potential determinant of health disparities \citep{Morello06}. In addition, it has been suggested that the disparities grow in more deprived areas as health improves faster in high socio-economic groups \citep{Leyland07,Higgs98}. 

Although individual determinants (such as smoking) or individual risk responses (such as closing windows to avoid exposure)  may frequently contribute to these health inequities, only a fraction of the overall disparities are attributed to individual factors \citep{Lantz01}. In fact, human health is not only influenced by individual health behaviours but also by contextual and ecological factors \citep{Marmot07}. Furthermore, socio-economic status plays a potential role of confounding or effect modification in epidemiological studies investigating the relationship between environmental variables and health outcomes, especially at  aggregated level \citep{Blakely00,Blakely04}. The further effect of confounding and effect modification will potentially lead to bias of the results, whose level depends on the relationship between environmental pollutions and socio-economic status. Hence it is extremely important to study this association, which, at the moment still remains uncertain and subjected to the fundamental methodological issue of correlation between variables. 

To study this relationship in the present work we consider the following data:
\begin{itemize}
\item nitrogen oxides (NO$_x$), which is generated mainly through combustion, thus is a good proxy for traffic related air pollution. The data were obtained from the environmental research group at Kings College as annual mean for the period 2003-2010 at the Lower Super Output Area geographical level in Greater London (LSOA, 4,767 in Greater London) as part of the TRAFFIC project (\url{http://www.kcl.ac.uk/lsm/research/divisions/aes/research/ERG/research-projects/traffic/index.aspx}).

\item Index of Multiple Deprivation (IMD), publicly available from the Department for Communities and Local Government (\texttt{data.gov.uk}). It is commonly used at the small area level to synthesize multiple aspects of deprivation. It is originally built at LSOA level and is formed by 38 indicators collapsed into seven domains: Income, Employment, Health, Education, Crime, Access to Services (Housing) and Living Environment. As we want to evaluate the relationship between the domains of the IMD and air pollution (NOx) we have not considered the living environment domain, which includes air quality. IMD is available for 2004, 2007 and 2010 and we have considered the most recent one in this work (correlation between Index at different years ranges from 0.94 and 0.97).
\end{itemize}

\begin{figure}[ht]
\centering
\subfigure[NO$_x$ concentration for 2003-2010]{\includegraphics[scale=0.4]{NOX.png}}
\qquad \subfigure[IMD score, 2010]{\includegraphics[scale=0.4]{IMD2004.png}}
\caption{Quintilesof the NO$_x$ concentration (average 2003-2010) and of IMD score (2010) at LSOA level in Greater London.\label{Fig:NOX-IMD}}
\end{figure}


Figure \ref{Fig:NOX-IMD} shows the map for NO$_x$ (left) and IMD score (right) and a clear spatial pattern is visible in both: air pollution concentration increases steadily going from outer to inner London, while IMD shows the highest deprived areas in the northeastern part of London and most of the central southern part. However looking at the maps of each of the six domains highlights a different picture (Figure \ref{Fig:IMD}): Crime shows the absence of a clear pattern, with scattered areas of high crime (dark grey) next to areas of low crime (light grey); on the other hand income, employment and health/disability are in agreement with the total IMD score, while barriers to housing and services are more pronounced in central London and education shows more deprivation in East London. This suggests how simplistic is the approach that considers the total IMD score and highlights the importance of including all the domains to disclose the relationship between social characteristics and environmental pollution at a small area level. 

If we want to investigate the relationship between each domain and air pollution we cannot include all the domains in a regression model due to their  collinearity issues: the pairwise Pearson correlation between domains (Table \ref{tab:cor}) shows high values for income and employment (0.91), income and health (0.77), income and education (0.68), employment and health (0.81) and employment and education (0.64).

\begin{table}[ht]
\caption{Correlation between IMD domains. In bold correlation higher than 0.6. \label{tab:cor}}
\centering
\begin{tabular}{l|llllll}
\hline
	&      Income & Employ. &   Health & Educ. &  Hous.  &   Crime\\
\hline
Income  &   1.0 & \textbf{0.91} & \textbf{0.77} & \textbf{0.68} & 0.48 & 0.52\\
Employ. & 0.91 & 1.0 & \textbf{0.81} & \textbf{0.64} & 0.42 & 0.53\\
Health    & 0.77 & 0.81 & 1.0 & 0.55 & 0.41 & 0.59\\
Educ. & 0.68 & 0.64 & 0.55 & 1.0 & 0.16 & 0.36\\
Hous.   & 0.48 & 0.42 & 0.41 & 0.16 & 1.0 & 0.29\\ 
Crime    &  0.52 & 0.53 & 0.59 & 0.37 & 0.29 & 1.0\\
\hline
\end{tabular}
\end{table}


\begin{figure}[h]\centering
\subfigure[Income]{\includegraphics[scale=0.4]{Income.png}}\qquad
\subfigure[Employment]{\includegraphics[scale=0.4]{Employment.png}}

\subfigure[Crime]{\includegraphics[scale=0.4]{Crime.png}}\qquad
\subfigure[Barriers to Housing and Services]{\includegraphics[scale=0.4]{Housing.png}}

\subfigure[Education]{\includegraphics[scale=0.4]{Educ.png}}\qquad
\subfigure[Health and Disability]{\includegraphics[scale=0.4]{Health.png}}
\caption{Maps of the six IMD domains considered: quintiles of the scores (note that higher positive values means higher deprivation.\label{Fig:IMD}}
\end{figure}

\clearpage
\section{Modelling highly correlated covariates with profile regression}\label{sect:methods}
To include all the domains into the same statistical model we use the profile regression, a model that non-parametrically links a response vector $\bY$ to covariate data $\bX$ through cluster membership. It was proposed by \cite{Molitor10} and it has been implemented in the R package PReMiuM \citep{Liverani15}. 

Profile regression implements a Bayesian clustering model through a Dirichlet process mixture model. The data $\bD=(\bY,\bX,\bW)$ contain the response $\bY$, covariate $\bX$ and fixed effects $\bW$ if they are available. The fixed effects are potentially confounding variables. In our running example the response is the nitrogen oxides, the covariates are the selected six domains of IMD and we do not include fixed effects. 

For each individual $i$, for $i=1,\ldots,n$, the response is given by $y_i$, the covariate vector $\mathbf{x}_i$ and the fixed effect vector $\mathbf{w}_i$. The data are then jointly modelled as the product of a response model and a covariate model, leading to the following likelihood:
\begin{equation*}
f(\mathbf{x}_i,y_i | \bTheta_{Z_i},\boldsymbol\Lambda,W_i,\psi) = \sum_c \psi_c f(\mathbf{x}_i | z_i = c,\phi_c)f(y_i|z_i=c,\theta_c,\boldsymbol\Lambda,\mathbf{w}_i) 
\end{equation*}
where $z_i = c$, the allocation variable, indicates that individual $i$ belongs to cluster $c$. The parameters $\boldsymbol\Lambda=(\boldsymbol{\theta},\boldsymbol\phi)$ are cluster specific and represent the contribution of the response and the covariates to the mixture model. There is also the possibility to include additional fixed effects $\mathbf{w}_i$ for each individual, which are constrained to only have a global (i.e., non-cluster specific) effect on the response $y_i$. The parameters $\psi$ are the mixture weights. 

Multicollinearity arises when regression models of the response with respect to highly correlated covariates are implemented, due to identifiability issues. However, as here the response is conditionally independent from the covariates, we do not encounter such issues, but we can explore in depth the potentially complex relationship between response and covariates. 

The prior model for the mixture weights is given by the stick-breaking priors (constructive definition of the Dirichlet process), that is,
\begin{eqnarray*}
\psi_c &=& V_c \prod_{l<c} (1-V_l) \quad \mbox{for all $c$},\\
\psi_1 &=& V_1,\\
 V_c &\sim& \mbox{Beta} (1,\alpha) \quad \mbox{i.i.d.}
\end{eqnarray*}
The parameter $\alpha$ can be fixed or can have a Gamma$(s_\alpha,r_\alpha)$ distribution with $s_\alpha$ and $r_\alpha$ as shape and rate parameters respectively. Other prior models for the mixture weights are possible, and, for example, the Pitman-Yor construction is also available in the R package PReMiuM. 

The covariate model $f(\mathbf{x}_i | z_i = c,\phi_c)$ can be defined as continuous or discrete. In the continuous case $\bX$ assumes a mixture of Gaussian distributions. In the discrete case for each individual $i$, $\mathbf{x}_i$ is a vector of $J$ locally independent discrete categorical random variables, where the number of categories for covariate $j=1,2,\ldots,J$ is $K_j$. Then we can write $\Phi_c=(\Phi_{c,1},\Phi_{c,2}\ldots,\Phi_{c,J})$ with $\Phi_{c,j} = (\phi_{c,j,1},\phi_{c,j,2},\ldots,\phi_{c,j,K_j})$ and
\begin{equation}
   f(\mathbf{x}_i | z_i = c,\phi_c)=\prod_{j=1}^J\phi_{{Z_i},j,X_{i,j}}. 
\end{equation}
We let $a = (a_1 , a_2 ,\ldots , a_J )$, where for $j = 1, 2,\ldots , J$, $a_j = (a_{j,1} , a_{ j,2} ,\ldots , a_{ j,K_{ j}} )$ and $\Phi_{ c,j} \sim \mbox{Dirichlet}(a_j )$. The covariate model can also be defined as a mixture of continuous and discrete covariates.  

The response model $f(y_i|z_i=c,\theta_c,\boldsymbol\Lambda,\mathbf{w}_i) $ can be defined as binary, categorical, count (modelled as Binomial or Poisson) or Gaussian. For example, for Gaussian response the mixture model is extended to contain $\theta_c$ for each $c$ and the global parameters $\boldsymbol\Lambda=(\boldsymbol\beta, \sigma^2_Y)$. These parameters allow us to write the response model as:
\[
f(y_i|z_i=c,\theta_c,\boldsymbol\Lambda,W_i)=f(y_i|z_i=c,\theta_c,\boldsymbol\beta,\sigma^2_Y,W_i)=\frac{1}{\sqrt{2\pi\sigma^2_Y}}\exp\left\{-\frac{1}{2\sigma^2_Y}(Y_i-\lambda_i)^2\right\},
\]
where $\lambda_i=\theta_{Z_i}+\boldsymbol\beta^\top W_i$ and $\boldsymbol\beta$ represent the effect of the counfounding variables, the fixed effects, on the response. For each cluster $c$, we adopt a $t$ location-scale distribution on $\theta_c$, with hyperparameters $\mu_\theta$ and $\sigma_\theta$ with 7 degrees of freedom. Similarly, we adopt the same prior for the fixed effect $\beta_k$, but with hyperparameters $\mu_\beta$ and $\sigma_\beta$. We set $\tau_Y = 1/\sigma_{Y}^2$ to Gamma$(s_{\tau_Y},r_{\tau_Y})$, where $s_{\tau_Y}$ and $r_{\tau_Y}$ are the shape and rate hyperparameters. 

More details on the Markov chain Monte Carlo (MCMC) algorithm for this model are provided in \cite{Liverani15}. When the signal in the data is strong the MCMC results for different runs, with different initial values and chain lengths, give stable results. However, this is not the case when the signal is not strong. \cite{Hastie15} discuss strategies to identify convergence issues. They recommend starting the MCMC with a large number of clusters as the algorithm can struggle to explore the partition space when starting with a small number of clusters if the signal is not strong. When many clusters are identified it is more challenging to characterise each cluster meaningfully and interpretation of the results can be facilitated by the posterior predictive distrbution. Moreover, they suggest the use of the posterior distribution of predictive profiles for the assessment of convergence instead of the posterior distribution of the parameters. This is because the posterior distribution of the parameters can appear to have converged when the model as a whole has not (often the case for $\boldsymbol\beta$), or cannot be used for this scope because they are cluster-specific and the number of clusters changes at each iteration (such as $\theta_c$).
 
It is often useful to characterise the partition which is most supported by the data. However, as at each iteration of the sampler individual profiles are assigned to clusters, the MCMC output is very rich. \cite{Molitor10} developed methods to process this output to make useful and interpretable inference.  Several methods for this are available in the R package PReMiuM but we find the most robust method is to process the similarity matrix using partitioning around medoids (PAM), which is available in the R package cluster. First of all, a score matrix is constructed, where each element of the matrix is set equal to 1 if individuals $i$ and $j$ belong to the same cluster and 0 otherwise. Then a similarity matrix $S$ is computed by dividing each element of the score matrix by the number of iterations, so that $S_{ij}$  denotes the probability that individuals i and j are assigned to the same cluster. PAM then assigns individuals to clusters in a way consistent with matrix $S$. 

%The R package PReMiuM also includes two methods for variable selection \citep{Papathomas12} and postprocessing functions to identify the optimal partition in the rich MCMC output.  








\subsection{The Spatial Conditional Autoregressive Model}

When clustering data from small area studies, we need to modify the model to account for spatial correlation. In this paper we propose to extend the response model described above to include an intrinsic spatial conditional autoregressive (ICAR) term \citep{Besag91} as follows. The likelihood component for the Gaussian response becomes 
\[
f(y_i|z_i=c,\theta_c,\boldsymbol\Lambda,W_i)=f(y_i|z_i=c,\theta_c,\boldsymbol\beta,\sigma^2_Y,u_i,W_i)=\frac{1}{\sqrt{2\pi\sigma_Y^2}}\exp(-\frac{1}{2\sigma_Y^2}(Y_i-\lambda_i)^2)
\]
where $\lambda_i= \theta_{Z_i}+  W_i\boldsymbol\beta + u_i $ and $u=(u_1, \cdots, u_n) \sim N(0, \mathbf{\tau P})$ with $\mathbf{P}=\{P_{ij}\}$ a precision matrix such that 
\[
P_{ij} = \begin{cases} 
n_i & \mbox{if} \;  i = j\\
-I\{i\sim j\}  &\mbox{if} \; i \ne j \\
\end{cases} 
\]
where $n_i$ is the number of neighbours of subject $i$, $I$ is the indicator function and $i \sim j$ indicates that regions $i$ and $j$ are neighbours. The prior of $\tau$ is given by 
\[
\tau \sim \mbox{Gamma}(a_{\tau}, b_{\tau})
\]
such that 
\[
E(\tau)=\frac{a_{\tau}}{b_{\tau}} \quad and \quad \mbox{Var}(\tau)=\frac{a_{\tau}}{b_{\tau}^2}.
\] 
Details of the sampling strategy for the ICAR parameters are given in Appendix \ref{app1}. We have implemented this model in the R package PReMiuM for Gaussian and Poisson responses. 

\section{Results}\label{sect:results}

We have fit profile regression to the data using the R package PReMiuM \citep{Liverani15}. The pollution data are modelled with a Gaussian distribution including a spatial ICAR term. The covariate profiles, given by the selected IMD domains, are modelled with a discrete distribution, as we have transformed each IMD domain into quintiles. We have not included any additional fixed effects. The results were very robust on several MCMC runs using a range of initial values and different chain lengths. We present here the results obtained with 5,000 iterations after a burn in of 5,000 with the following hyperparameter settings. 
\begin{eqnarray*}
&&s_\alpha = 2, \quad \quad r_\alpha=1,\\
&&a_1=\ldots=a_6=1,\\
&&\mu_\theta=0, \quad \sigma_\theta=2.5,\\
&&\mu_\beta=0, \quad \sigma_\beta=2.5,\\
&&s_{\tau-Y}=2.5, \quad r_{\tau-Y}=2.5.
\end{eqnarray*}

\begin{figure}[h]
\centering
\includegraphics[height=7cm]{clustering.png}
\caption{Geographical representation of the eleven clusters of the areas in Greater London identified by profile regression. The colours reflect the mean of the observed pollution levels, with dark grey identifying the most polluted clusters and light grey the least polluted clusters. \label{fig:clusters}}
\end{figure}

In this section we aim to illustrate how profile regression can help shed light on complex patterns between highly covariates covariates and response. We propose several ways to explore the results.  

The main output of profile regression are the clusters, given by regions with similar covariate and response characteristics. We have used different types of plots, which highlight specific features, at different level of details. In Figure \ref{fig:clusters} we represent the clusters geographically. The eleven clusters identified are plotted using colours that reflect their observed pollution levels. As expected, the most central areas have higher pollution levels. In particular we see that areas that belong to clusters with higher observed pollution levels are mostly located in North-East London, where areas with higher levels of deprivation are found. In contrast, the less deprived areas in South-West London are clustered together and have lower mean for the observed pollution levels. 


%\begin{table}[!htbp]
%\caption{Cluster sizes, mean observed pollution per cluster ($\bar y_c$) and average IMD per cluster.\label{tab:profiles}}
%\centering
%\begin{tabular}{lllllllllllll}
%  Cluster &  & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 \\ 
%   \hline
%Size &  & 153 & 948 & 649 & 294 & 614 & 199 & 178 & 604 & 722 & 274 & 132 \\ 
%   \hline
%$\bar y_c $ &  & 61.9 & 64.7 & 66.0 & 66.2 & 69.1 & 70.8 & 77.3 & 82.6 & 82.8 & 82.9 & 84.6 \\ 
%   \hline
%
%IMD &  & 3.9 & 1.14 & 1.86 & 3.34 & 2.71 & 4.48 & 2.64 & 4.23 & 4.95 & 3.55 & 2.32 \\ 
%   \hline
%\end{tabular}
%\end{table}
%
%\iffalse
% & 1 & 0.01 & 0.9 & 0.22 & 0 & 0.02 & 0 & 0.06 & 0 & 0 & 0.01 & 0.65 \\ 
%   & 2 & 0.01 & 0.09 & 0.64 & 0.02 & 0.24 & 0.01 & 0.61 & 0 & 0 & 0.08 & 0.3 \\ 
%  Income & 3 & 0.07 & 0.01 & 0.13 & 0.41 & 0.67 & 0.02 & 0.31 & 0.08 & 0 & 0.77 & 0.04 \\ 
%   & 4 & 0.8 & 0 & 0 & 0.53 & 0.07 & 0.19 & 0.02 & 0.82 & 0.07 & 0.13 & 0.01 \\ 
%   & 5 & 0.12 & 0 & 0 & 0.04 & 0 & 0.78 & 0 & 0.1 & 0.92 & 0.01 & 0.01 \\ 
%   \hline
% & 1 & 0.01 & 0.9 & 0.21 & 0 & 0.02 & 0 & 0.09 & 0 & 0 & 0.02 & 0.79 \\ 
%   & 2 & 0.01 & 0.09 & 0.74 & 0.07 & 0.51 & 0.01 & 0.67 & 0.01 & 0 & 0.19 & 0.19 \\ 
%  Emloyp & 3 & 0.07 & 0 & 0.05 & 0.38 & 0.37 & 0.03 & 0.2 & 0.09 & 0 & 0.37 & 0.01 \\ 
%   & 4 & 0.79 & 0 & 0 & 0.52 & 0.09 & 0.55 & 0.03 & 0.73 & 0.16 & 0.41 & 0.01 \\ 
%   & 5 & 0.12 & 0 & 0 & 0.03 & 0 & 0.42 & 0.01 & 0.16 & 0.84 & 0.02 & 0.01 \\ 
%   \hline
% & 1 & 0.01 & 0.8 & 0.22 & 0 & 0.02 & 0 & 0.05 & 0 & 0 & 0.01 & 0.29 \\ 
%   & 2 & 0.02 & 0.17 & 0.6 & 0.09 & 0.33 & 0.01 & 0.31 & 0.01 & 0 & 0.05 & 0.53 \\ 
%  Health & 3 & 0.2 & 0.02 & 0.17 & 0.53 & 0.55 & 0.14 & 0.54 & 0.11 & 0.01 & 0.33 & 0.16 \\ 
%   & 4 & 0.54 & 0 & 0.01 & 0.32 & 0.1 & 0.64 & 0.09 & 0.53 & 0.12 & 0.55 & 0.02 \\ 
%   & 5 & 0.24 & 0 & 0 & 0.06 & 0.01 & 0.2 & 0.01 & 0.35 & 0.86 & 0.06 & 0.01 \\ 
%   \hline
% & 1 & 0.02 & 0.6 & 0.31 & 0.03 & 0.1 & 0.01 & 0.21 & 0.01 & 0 & 0.07 & 0.37 \\ 
%   & 2 & 0.04 & 0.23 & 0.32 & 0.17 & 0.31 & 0.02 & 0.39 & 0.1 & 0.02 & 0.35 & 0.35 \\ 
%  Educ & 3 & 0.11 & 0.11 & 0.22 & 0.21 & 0.29 & 0.1 & 0.24 & 0.33 & 0.09 & 0.33 & 0.19 \\ 
%   & 4 & 0.14 & 0.06 & 0.12 & 0.26 & 0.18 & 0.3 & 0.11 & 0.37 & 0.35 & 0.17 & 0.07 \\ 
%   & 5 & 0.69 & 0 & 0.03 & 0.33 & 0.12 & 0.58 & 0.06 & 0.18 & 0.54 & 0.08 & 0.02 \\ 
%   \hline
% & 1 & 0.17 & 0.48 & 0.35 & 0.19 & 0.17 & 0.05 & 0.14 & 0.03 & 0.01 & 0.04 & 0.2 \\ 
%   & 2 & 0.33 & 0.28 & 0.29 & 0.29 & 0.27 & 0.15 & 0.25 & 0.08 & 0.03 & 0.11 & 0.18 \\ 
%  Housing & 3 & 0.21 & 0.12 & 0.21 & 0.28 & 0.27 & 0.23 & 0.28 & 0.22 & 0.1 & 0.37 & 0.18 \\ 
%   & 4 & 0.2 & 0.06 & 0.09 & 0.18 & 0.19 & 0.28 & 0.21 & 0.39 & 0.3 & 0.23 & 0.17 \\ 
%   & 5 & 0.1 & 0.06 & 0.06 & 0.06 & 0.09 & 0.3 & 0.13 & 0.28 & 0.57 & 0.25 & 0.27 \\ 
%   \hline
% & 1 & 0.1 & 0.51 & 0.28 & 0.13 & 0.17 & 0.1 & 0.15 & 0.04 & 0.03 & 0.1 & 0.22 \\ 
%   & 2 & 0.32 & 0.24 & 0.29 & 0.22 & 0.25 & 0.12 & 0.23 & 0.13 & 0.08 & 0.14 & 0.19 \\ 
%  Crime & 3 & 0.3 & 0.13 & 0.22 & 0.3 & 0.26 & 0.29 & 0.24 & 0.17 & 0.16 & 0.22 & 0.18 \\ 
%   & 4 & 0.16 & 0.1 & 0.12 & 0.19 & 0.22 & 0.26 & 0.21 & 0.3 & 0.25 & 0.28 & 0.23 \\ 
%   & 5 & 0.13 & 0.03 & 0.08 & 0.17 & 0.1 & 0.23 & 0.17 & 0.36 & 0.49 & 0.26 & 0.19 \\ 
%   \hline
%\fi

\begin{figure}[!htbp]
\centering
\includegraphics[angle=90,height=18cm]{plotsProfRegrEmpiricalRemoved.png}
\caption{Summary plot of the posterior distribution of the parameters $\Phi_c$, for $c=1,\ldots11$. Each column $j$ in the figure represents $\Phi_{c,j,k}$, for $c=1,\ldots,11$ and $k=0,1,\ldots,4$. Within a column $j$, each row $k$ is a visualisation of the boxplots for $\phi_{c,j,k}$ for each cluster $c$.\label{fig:summary}}
\end{figure}

The cluster data are high-dimensional and complex. Figure \ref{fig:summary} provides a visual representation of the parameters $\Phi_c$. Through the boxplots of the MCMC samples of the parameters $\Phi_c$, this figure provides also a representation of the uncertainty around these parameters. Each column $j$ in the figure represents $\Phi_{c,j,k}$, for $c=1,\ldots,11$ and $k=0,1,\ldots,4$ (the five quintiles of each covariate). Within a column $j$, each row $k$ is a visualisation of the boxplots for $\phi_{c,j,k}$ for each cluster $c$. This visual representation provides a further insight into the eleven clusters. We can identify patterns in the relationships between pollution and IMD domains at a glance. For example, Housing and Crime appear to generally increase as pollution levels increase, while the other domains show less linear patterns. We can also see here the details of the distributions of the levels of covariates that define the different clusters. Cluster 2 has the lowest levels of deprivation, although not the lowest levels of mean pollution. In contrast, cluster 9 has the highest levels of deprivation, and high levels of pollution, although not the highest among all clusters. The mean IMD per cluster highlights the complex relationship between IMD, pollution and deprivation. For example, for cluster 6, which has a high mean IMD, there is strong deprivation for the first four domains. However, the domains of Crime and Housing are rather evenly spread among all levels of the covariates, suggesting that they do not contribute to the deprivation that characterises these areas. This is an example of a complex pattern that cannot be identified when the domains are simply collapsed into the IMD. 

\begin{figure}[!htbp]
\centering
\includegraphics[height=12cm]{heatmap.png}
\caption{Summary table of the clusters. The quintiles for pollution and each domain of the IMD are shown for each cluster. \label{fig:heatmap}}
\end{figure}

In Figure \ref{fig:heatmap} we provide a summary of the posterior means for each cluster. Each row represents a cluster. The columns represents, respectively, the mean observed pollution and each domain of IMD. The colour of each column of the matrix corresponds to a quintile of the distribution of that variable. As before, the clusters are ordered by their observed pollution level. Note that the colours in the matrix do not become darker (or lighter) in a smooth manner. Together, Figure \ref{fig:heatmap} and Figure \ref{fig:clusters} suggest that the areas of low pollution and low deprivation are in outer London. As we get closer to the centre, pollution increases and many of the deprivation variables increase levels. However, there are many notable exceptions to this. For example, cluster 11 has the highest levels of pollution, but among the lowest levels of deprivation. On the contrary, cluster 1 has the lowest level of pollution, but rather high levels of deprivation on all domains except Housing. 

\begin{figure}[!htbp]
\centering
\includegraphics[height=6cm]{u.png}
\caption{Posterior mean of the spatial conditional autoregressive term.\label{fig:u}}
\end{figure}

Figure \ref{fig:u} shows the posterior mean of the spatial term $\exp(u_i)$ for each area, which accounts for the residual spatial variation in NOx after having adjusted for the cluster assignment. The map presents a clear pattern going from central London (darker) to outer London (lighter) with values ranging from -30 to 48 $\mu g/m^3$, thus suggesting that the model picks up the spatial dependence in air pollution concentration which is not explained by deprivation.


\begin{figure}[!htbp]
\centering
\includegraphics[height=8cm]{preds12.pdf}
\caption{Beanplots of the posterior predictive distributions for these four pseudo profiles: $(0,0,0,0,0,0), (0,\text{NA},\text{NA},\text{NA},\text{NA},\text{NA}), (4,4,4,4,4,4), (4,4,4,4,4,0)$ where the elements of each vector represent respectively the IMD domains (Income, Employment, Health, Education, Housing, Crime). \label{fig:preds1}}
\end{figure}

We can explore the relationships between covariates and response further by looking at the posterior predictive distributions. The profile regression model allows us to predict the pollution level for specific combinations of the IMD domains. If we wish to understand the role of a particular covariate or group of covariates, we can specify a number of predictive scenarios (pseudo-profiles), that capture the range of possibilities for the covariates that we are interested in \citep{Hastie13}. For each of these pseudo-profiles we can see how these would have been allocated in our mixture model to understand the level of pollution associated with them once we have accounted for the spatial residuals. 

In Figure \ref{fig:preds1} we show beanplots of four pseudo-profiles: $(0,0,0,0,0,0)$, $(0,\text{NA},\text{NA},\text{NA},\text{NA},\text{NA})$, $(4,4,4,4,4,4)$, $(4,4,4,4,4,0)$. The elements of each vector represent the IMD domains in the following order: Income, Employment, Health, Education, Housing, Crime. Each of these beanplots shows the pseudo-profile corresponding to particular values of the IMD domain of interest for an average area (ie. including the mean spatial residual, which is 0). When `NA' is set, this allows that domain to vary, de facto marginalising for it, i.e. capturing all possible values it can take. For example, the beanplot for the first pseudo-profile on the left presents the posterior predictive distribution of NOx for areas with the lowest levels of deprivation and shows values between 68 and 75$\mu$g/m$^3$. Using this as benchmark we can make some comparisons: areas characterised by low Income and marginalising for the remaining domains (second beanplot from the left) have a much wider posterior predictive distribution, with values going from 67 to 78$\mu$g/m$^3$ while areas with in the highest quintile of deprivation for all the domains (third beanplot from the left) present consistently highest level of pollution (ranging between 71 to 77$\mu$g/m$^3$). The last beanplot shows the posterior predictive distribution of NOx in an area where crime has decreased to the first quintile (for instance through the implementation of a policy) while the other domains remains in the last quintile. Comparing it with the previous one it an be seen a similar distribution, but with a lower tail which could be a consequence of the policy implementation. 

All predictive profiles can be computed and we provide here only these examples to show how these can be interpreted, if there was an interest in the posterior predictive distribution of specific combinations of deprivation levels, these could be explored in depth. Moreover, if there was interest in a specific area, the pseudo-profiles could be adjusted by adding the spatial residual. 

\section{Discussion}\label{sect:conclusion}

In this paper we have considered a spatially-correlated response variable and a set of highly correlated covariates. We have extended the profile regression model, a Bayesian clustering method used to deal with collinearity in the predictors, to account for spatial correlation adding a spatial conditionally autoregressive term. 

We have applied our method to explain the relationship between air pollution and social deprivation in Greater London. The Index of Multiple Deprivation is commonly used as a proxy for deprivation, as its domains usually cannot be analysed individually due to the high correlation between them. We have illustrated how profile regression can produce meaningful and useful results which shed light on the complex non linear relationship between pollution and the different domains.  

We want to stress that we are not framed in a standard regression approach, where the interest is to estimate the effect of each predictor on the outcome, as we do not attempt to explain the level of air pollution through the IMD domains. On the other hand through cluster assignment the profile regression  is able to disentangle the complex relationship between IMD and air pollution; this method has the added benefit of providing readily available prediction estimates which can be used to evaluate how the response could change for specific combinations of the predictors, and which could be used to evaluate the effect of policies.

A limitation of our model is that in its present formulation the spatial structure is not included on the cluster allocation, thus it accounts for local spatial dependency in the response, but not in the covariates, which is an extension we are going to work on in the future. 

%The code and the data to replicate our figures and results are available on the journal's webpage. 


\section*{Acknowledgements}
Marta Blangiardo acknowledges support from the UK NERC-MRC funded project “Traffic pollution and health in London” (NE/I00789X/1). Silvia Liverani acknowledges support from the Leverhulme Trust (ECF-2011-576). 

\appendix
\section{Sampling for the spatial ICAR parameters} \label{app1}
We include details of the sampling algorithm for the spatial ICAR parameters for Gaussian and Poisson distributed response. We have implemented both in the R package PReMiuM. 



\subsection{Gaussian response}

The conditional distribution for $u_i$ is given by 
\begin{eqnarray*}
\log(p(u_i|u_{-i},\boldsymbol\beta,\boldsymbol\theta,\sigma_Y^2,Z_i,\tau,Y_i,W_i,T_i)) & \propto  & \log p(Y_i|u_i,\boldsymbol\beta,\theta,Z_i,\sigma^2_Y,W_i,T_i) + \log p(u_i|u_{-i},\tau)\\
& \propto & -\frac{1}{2\sigma_Y^2}(Y_i-(\theta_{Z_i}+  W_i\boldsymbol\beta + u_i))^2 -\frac{1}{2}\tau n_i (u_i-\bar{u}_i)^2\\
& \propto & -\frac{1}{2\sigma^2_i}(u_i - m_i)^2\\
\end{eqnarray*}
with
\[
\left\{\begin{array}{l}
m_i=\frac{\frac{1}{\sigma_Y^2}(Y_i-\theta_{Z_i}-W_i\boldsymbol\beta)-\tau n_i \bar{u}_i}{\frac{1}{\sigma_Y^2}+\tau n_i}\\
\sigma^2_i=\frac{1}{\frac{1}{\sigma_Y^2}+\tau n_i}
\end{array}\right.
\]
with , $\bar{u}_i=\frac{1}{n_i}\sum_{j \in \rho_i}u_j$ $\rho_i$ is the set of neighbours of $i$. Thus, for Normal response the prior is conjugated, the conjugated complete conditional distribution is Normal with mean $m_i$ and variance $\sigma_i^2$. The conditional distribution for $\tau$ is given by 
\[
\log(p(\tau|u))= (a_{\tau}+\frac{n-1}{2}-1)\log(\tau)-\tau(b_{\tau}+\frac{1}{2}u^T\mathbf{P}u)
\]
Thus $\tau \sim \mbox{Gamma}(a_{\tau}+\frac{n-1}{2},b_{\tau}+\frac{1}{2}u^T\mathbf{P}u)$. 

\subsection{Poisson response} For Poisson response, suitable for count data, the likelihood is given by 
\[
f_Y(y_i|z_i=c,\theta_c,\boldsymbol\Lambda,W_i)=p(Y_i|z_i=c,\theta_c,\boldsymbol\beta,u_i,W_i)=\frac{\mu_i^{Y_i}}{Y_i!}\exp\{-\mu_i\},
\]
where each individual $i$ is associated with an expected offset $E_i$, 
\[ 
\mu_i=E_i \exp\{\lambda_i\},\;\;\;\textrm{for } \lambda_i=\theta_{Z_i} + \boldsymbol\beta^\top W_i.
\]
As for the Gaussian response, the parameters $u=(u_1, \cdots, u_n) \sim N(0, \mathbf{\tau P})$ with $\mathbf{P}=\{P_{ij}\}$ a precision matrix such that 
\[
P_{ij} = \begin{cases} 
n_i & \mbox{if} \;  i = j\\
-I\{i\sim j\}  &\mbox{if} \; i \ne j \\
\end{cases} 
\]
where $n_i$ is the number of neighbours of subject $i$, $I$ is the indicator function and $i \sim j$ indicates that regions $i$ and $j$ are neighbours. The prior of $\tau$ is given by 
\[
\tau \sim \mbox{Gamma}(a_{\tau}, b_{\tau})
\]
such that 
\[
E(\tau)=\frac{a_{\tau}}{b_{\tau}} \quad and \quad \mbox{Var}(\tau)=\frac{a_{\tau}}{b_{\tau}^2}.
\] 

The conditional distribution for $u_i$ is given by 
\[
\log(p(u_i|u_{-i},\boldsymbol\beta,\theta,Z,\tau,Y))=Y_iu_i-E_i\exp(X_i\boldsymbol\beta+\theta_{Z_i}+u_i)-\frac{1}{2}\tau n_i (u_i-\bar{u}_i)^2
\]
with $\bar{u}_i=\frac{1}{n_i}\sum_{j \in \rho_i}u_j$ and $\rho_i$ is the set of neighbours of $i$. We implemented an adaptive rejection sampler for $u_i$. The conditional distribution for $\tau$ is given by 
\[
\log(p(\tau|u))= (a_{\tau}+\frac{n-1}{2}-1)\log(\tau)-\tau(b_{\tau}+\frac{1}{2}u^T\mathbf{P}u).
\]
Thus $\tau \sim \mbox{Gamma}(a_{\tau}+\frac{n-1}{2},b_{\tau}+\frac{1}{2}u^T\mathbf{P}u)$. 


\bibliographystyle{chicago}
\bibliography{biblio}

\end{document}
SSTEarxiv_v1.bbl
\begin{thebibliography}{}

\bibitem[\protect\citeauthoryear{Abellan, Richardson, and Best}{Abellan
  et~al.}{2008}]{Abellan08}
Abellan, J., S.~Richardson, and N.~Best (2008).
\newblock Use of space-time models to investigate the stability of patterns of
  disease.
\newblock {\em Environmental Health Perspectives\/}~{\em 116}, 1111--1119.

\bibitem[\protect\citeauthoryear{Barcelo, Saez, and C.}{Barcelo
  et~al.}{2009}]{Barcelo09}
Barcelo, M., M.~Saez, and S.~C. (2009).
\newblock Spatial variability in mortality inequalities, socioeconomic
  deprivation, and air pollution in small areas of the barcelona metropolitan
  region, spain.
\newblock {\em Science of the Total Environment\/}~{\em 407\/}(21), 5501--5523.

\bibitem[\protect\citeauthoryear{Benach, Yasui, Borrell, Saez, and
  Pasarin}{Benach et~al.}{2001}]{Benach01}
Benach, J., Y.~Yasui, C.~Borrell, M.~Saez, and M.~I. Pasarin (2001).
\newblock Material deprivation and leading causes of death by gender: evidence
  from a nationwide small area study.
\newblock {\em Journal of Epidemiology and Community Health\/}~{\em 55},
  239--245.

\bibitem[\protect\citeauthoryear{Besag, York, and Molli{\'e}}{Besag
  et~al.}{1991}]{Besag91}
Besag, J., J.~York, and A.~Molli{\'e} (1991).
\newblock Bayesian image restoration, with two applications in spatial
  statistics.
\newblock {\em Annals of the institute of statistical mathematics\/}~{\em
  43\/}(1), 1--20.

\bibitem[\protect\citeauthoryear{Bigelow and Dunson}{Bigelow and
  Dunson}{2009}]{Bigelow09}
Bigelow, J.~L. and D.~B. Dunson (2009, March).
\newblock {Bayesian Semiparametric Joint Models for Functional Predictors}.
\newblock {\em Journal of the American Statistical Association\/}~{\em
  104\/}(485), 26--36.

\bibitem[\protect\citeauthoryear{Blakely, Kawachi, Atkinson, and
  Fawcett}{Blakely et~al.}{2004}]{Blakely04}
Blakely, T., I.~Kawachi, J.~Atkinson, and J.~Fawcett (2004).
\newblock Income and mortality: the shape of the association and confounding
  new zealand census-mortality study.
\newblock {\em international Journal of Epidemiology\/}~{\em 33}, 874--883.

\bibitem[\protect\citeauthoryear{Blakely and Woodward}{Blakely and
  Woodward}{2000}]{Blakely00}
Blakely, T. and A.~Woodward (2000).
\newblock Ecological effects in multi-level studies.
\newblock {\em Journal of Epidemiology Community and Health\/}~{\em 54},
  367--374.

\bibitem[\protect\citeauthoryear{Blowers and Leroy}{Blowers and
  Leroy}{1994}]{Blowers94}
Blowers, A. and P.~Leroy (1994).
\newblock Power, politics and environmental inequality: A theoretical and
  empirical analysis of the process of peripheralisation.
\newblock {\em Environmental Politics\/}~{\em 3}, 197--228.

\bibitem[\protect\citeauthoryear{Brown}{Brown}{1995}]{Brown95}
Brown, P. (1995).
\newblock Race, class, and environmental health: a review and systematization
  of the literature.
\newblock {\em Environmental Research\/}~{\em 69}, 15--30.

\bibitem[\protect\citeauthoryear{Brulle and Pellow}{Brulle and
  Pellow}{2006}]{Brulle06}
Brulle, R. and D.~Pellow (2006).
\newblock Environmental justice: human health and environmental inequalities.
\newblock {\em Annual Review of Public Health\/}~{\em 27}, 103--124.

\bibitem[\protect\citeauthoryear{Chen, Goldberg, and Villeneuve}{Chen
  et~al.}{2008}]{Chen08}
Chen, H., M.~Goldberg, and P.~Villeneuve (2008).
\newblock A systematic review of the relation between long-term exposure to
  ambient air pollution and chronic diseases.
\newblock {\em Reviews on Environmental Health\/}~{\em 23}, 243--297.

\bibitem[\protect\citeauthoryear{Davidson and Anderton}{Davidson and
  Anderton}{2000}]{Davidson00}
Davidson, P. and D.~Anderton (2000).
\newblock Demographics of dumping ii: a national environmental equity survey
  and the distribution of hazardous materials handlers.
\newblock {\em Demography\/}~{\em 37}, 461--466.

\bibitem[\protect\citeauthoryear{Dunson, Herring, and Siega-Riz}{Dunson
  et~al.}{2008}]{Dunson08}
Dunson, D.~B., A.~B. Herring, and A.~M. Siega-Riz (2008, December).
\newblock {Bayesian Inference on Changes in Response Densities Over Predictor
  Clusters}.
\newblock {\em Journal of the American Statistical Association\/}~{\em
  103\/}(484), 1508--1517.

\bibitem[\protect\citeauthoryear{Elliott and Wartenberg}{Elliott and
  Wartenberg}{2004}]{Elliott04}
Elliott, P. and D.~Wartenberg (2004).
\newblock Spatial epidemiology: current approaches and future challenges.
\newblock {\em Environmental Health Perspectives\/}~{\em 112\/}(9), 998--1006.

\bibitem[\protect\citeauthoryear{Hastie, Liverani, Azizi, Richardson, and
  St{\"u}cker}{Hastie et~al.}{2013}]{Hastie13}
Hastie, D.~I., S.~Liverani, L.~Azizi, S.~Richardson, and I.~St{\"u}cker (2013).
\newblock {A semi-parametric approach to estimate risk functions associated
  with multi-dimensional exposure profiles: application to smoking and lung
  cancer}.
\newblock {\em BMC Medical Research Methodology\/}~{\em 13\/}(1), 129.

\bibitem[\protect\citeauthoryear{Hastie, Liverani, and Richardson}{Hastie
  et~al.}{2015}]{Hastie15}
Hastie, D.~I., S.~Liverani, and S.~Richardson (2015).
\newblock {Sampling from Dirichlet process mixture models with unknown
  concentration parameter: mixing issues in large data implementations}.
\newblock {\em Statistics and Computing\/}~{\em 25\/}(5), 1023--1037.

\bibitem[\protect\citeauthoryear{Higgs, Senior, and Williams}{Higgs
  et~al.}{1998}]{Higgs98}
Higgs, G., M.~Senior, and H.~Williams (1998).
\newblock . spatial and temporal variation of mortality and deprivation. 1:
  widening health inequalities.
\newblock {\em Environment and Planning A\/}~{\em 30\/}(1), 661--682.

\bibitem[\protect\citeauthoryear{Lantz, Lynch, House, Lepkowski, Mero, Musick,
  and Williams}{Lantz et~al.}{2001}]{Lantz01}
Lantz, P., J.~Lynch, J.~House, J.~Lepkowski, R.~Mero, M.~Musick, and
  D.~Williams (2001).
\newblock . socioeconomic disparities in health change in a longitudinal study
  of us adults: the role of health-risk behaviors.
\newblock {\em Social Science and Medicine\/}~{\em 53}, 29--40.

\bibitem[\protect\citeauthoryear{Laurent, Bard, Filleul, and Segala}{Laurent
  et~al.}{2007}]{Laurent07}
Laurent, O., D.~Bard, L.~Filleul, and C.~Segala (2007).
\newblock Effect of socioeconomic status on the relationship between
  atmospheric pollution and mortality.
\newblock {\em Journal of Epidemiology and Community Health\/}~{\em 61},
  665--675.

\bibitem[\protect\citeauthoryear{Leyland, Dundas, Mcloone, and Boddy}{Leyland
  et~al.}{2007}]{Leyland07}
Leyland, A., R.~Dundas, P.~Mcloone, and F.~Boddy (2007).
\newblock Cause-specific inequalities in mortality in scotland: two decades of
  change. a population-based study.
\newblock {\em BMC Public Health\/}~{\em 7}, 1--12.

\bibitem[\protect\citeauthoryear{Li, Best, Hansell, Ahmed, and Richardson}{Li
  et~al.}{2012}]{Li12}
Li, G., N.~Best, A.~Hansell, I.~Ahmed, and S.~Richardson (2012).
\newblock Baystdetect: detecting unusual temporal patterns in small area data
  via bayesian model choice.
\newblock {\em Biostatistics\/}~{\em 13\/}(4), 695--710.

\bibitem[\protect\citeauthoryear{Liverani, Hastie, Azizi, Papathomas, and
  Richardson}{Liverani et~al.}{2015}]{Liverani15}
Liverani, S., D.~I. Hastie, L.~Azizi, M.~Papathomas, and S.~Richardson (2015).
\newblock {PReMiuM: An R package for Profile Regression Mixture Models using
  Dirichlet Processes}.
\newblock {\em Journal of Statistical Software\/}~{\em 64\/}(7), 1--30.

\bibitem[\protect\citeauthoryear{Marmot}{Marmot}{2007}]{Marmot07}
Marmot, M. (2007).
\newblock Achieving health equity: from root causes to fair outcomes.
\newblock {\em Lancet\/}~{\em 370}, 1153--1163.

\bibitem[\protect\citeauthoryear{Molitor, Brown, Chan, Papathomas, Liverani,
  Molitor, Richardson, Van~Horn, Daviglus, Dyer, Stamler, Elliott, and
  Group}{Molitor et~al.}{2014}]{Molitor14}
Molitor, J., I.~J. Brown, Q.~Chan, M.~Papathomas, S.~Liverani, N.~Molitor,
  S.~Richardson, L.~Van~Horn, M.~L. Daviglus, A.~Dyer, J.~Stamler, P.~Elliott,
  and I.~R. Group (2014).
\newblock {Blood Pressure Differences Associated With Optimal Macronutrient
  Intake Trial for Heart Health (OMNIHEART)--Like Diet Compared With a Typical
  American Diet}.
\newblock {\em Hypertension\/}~{\em 64\/}(6), 1198--1204.

\bibitem[\protect\citeauthoryear{Molitor, Papathomas, Jerrett, and
  Richardson}{Molitor et~al.}{2010}]{Molitor10}
Molitor, J., M.~Papathomas, M.~Jerrett, and S.~Richardson (2010, July).
\newblock {Bayesian profile regression with an application to the National
  Survey of Children's Health.}
\newblock {\em Biostatistics\/}~{\em 11\/}(3), 484--498.

\bibitem[\protect\citeauthoryear{Morello-Frosch and Jesdale}{Morello-Frosch and
  Jesdale}{2006}]{Morello06}
Morello-Frosch, R. and B.~Jesdale (2006).
\newblock Separate and unequal: residential segregation and estimated cancer
  risks associated with ambient air toxics in u.s. metropolitan areas.
\newblock {\em Environmental Health Perspectives\/}~{\em 114}, 386--393.

\bibitem[\protect\citeauthoryear{Morello-Frosch, Pastor, Porras, and
  Sadd}{Morello-Frosch et~al.}{2002}]{Morello02}
Morello-Frosch, R., M.~Pastor, C.~Porras, and J.~Sadd (2002).
\newblock . separate and unequal: residential segregation and estimated cancer
  risks associated with ambient air toxics in u.s. metropolitan areas.
\newblock {\em Environmental Health Perspectives\/}~{\em 110\/}(2), 149--154.

\bibitem[\protect\citeauthoryear{Morgenstern, Zutavern, Cyrys, Brockow,
  Gehring, Koletzko, Bauer, Reinhardt, Wichmann, and Heinrich}{Morgenstern
  et~al.}{2007}]{Morgenstern07}
Morgenstern, V., A.~Zutavern, J.~Cyrys, I.~Brockow, U.~Gehring, S.~Koletzko,
  C.~Bauer, D.~Reinhardt, H.~Wichmann, and J.~Heinrich (2007).
\newblock A case-control analysis of exposure to traffic and acute myocardial
  infarction.
\newblock {\em Occupational Environmental Medicine\/}~{\em 64}, 8--16.

\bibitem[\protect\citeauthoryear{O'Neill, Jerrett, Kawachi, Levy, Cohen,
  Gouveia, Wilkinson, Fletcher, Cifuentes, and Schwartz}{O'Neill
  et~al.}{2003}]{Oneill03}
O'Neill, M., M.~Jerrett, I.~Kawachi, J.~Levy, A.~Cohen, N.~Gouveia,
  P.~Wilkinson, T.~Fletcher, L.~Cifuentes, and J.~Schwartz (2003).
\newblock Health, wealth, and air pollution: advancing theory and methods.
\newblock {\em Environmental Health Perspectives\/}~{\em 111}, 1861--1870.

\bibitem[\protect\citeauthoryear{Papathomas, Molitor, Hoggart, Hastie, and
  Richardson}{Papathomas et~al.}{2012}]{Papathomas12}
Papathomas, M., J.~Molitor, C.~Hoggart, D.~Hastie, and S.~Richardson (2012).
\newblock Exploring data from genetic association studies using bayesian
  variable selection and the dirichlet process: application to searching for
  gene$\times$ gene patterns.
\newblock {\em Genetic epidemiology\/}~{\em 36\/}(6), 663--674.

\bibitem[\protect\citeauthoryear{Papathomas, Molitor, Richardson, Riboli, and
  Vineis}{Papathomas et~al.}{2011}]{Papathomas11}
Papathomas, M., J.~Molitor, S.~Richardson, E.~Riboli, and P.~Vineis (2011).
\newblock Examining the joint effect of multiple risk factors using exposure
  risk profiles: lung cancer in nonsmokers.
\newblock {\em Environmental health perspectives\/}~{\em 119\/}(1), 84.

\bibitem[\protect\citeauthoryear{Perlin, Setzer, Creason, and Sexton}{Perlin
  et~al.}{1995}]{Perlin95}
Perlin, S., R.~Setzer, J.~Creason, and K.~Sexton (1995).
\newblock Distribution of industrial air emissions by income and race in the
  united states: an approach using the toxic release inventory.
\newblock {\em Environmental Science and Technology\/}~{\em 29}, 69--80.

\bibitem[\protect\citeauthoryear{Peters, Von~Klot, Heier, Trentinaglia,
  Hormann, Wichmann, and Lowel}{Peters et~al.}{2004}]{Peters04}
Peters, A., S.~Von~Klot, M.~Heier, I.~Trentinaglia, A.~Hormann, H.~Wichmann,
  and H.~Lowel (2004).
\newblock A case-control analysis of exposure to traffic and acute myocardial
  infarction.
\newblock {\em New England Journal of Medicine\/}~{\em 351}, 1721--1730.

\bibitem[\protect\citeauthoryear{Pirani, Best, Blangiardo, Liverani, Atkinson,
  and Fuller}{Pirani et~al.}{2015}]{Pirani15}
Pirani, M., N.~Best, M.~Blangiardo, S.~Liverani, R.~W. Atkinson, and G.~W.
  Fuller (2015).
\newblock {Analysing the health effects of simultaneous exposure to physical
  and chemical properties of airborne particles}.
\newblock {\em Environment International\/}~{\em 79}, 56--64.

\bibitem[\protect\citeauthoryear{Pope, Burnett, Turner, Cohen, Krewski,
  Jerrett, Gapstur, and Thun}{Pope et~al.}{2011}]{Pope11}
Pope, C., R.~Burnett, M.~Turner, A.~Cohen, D.~Krewski, M.~Jerrett, S.~Gapstur,
  and M.~Thun (2011).
\newblock Lung cancer and cardiovascular disease mortality associated with
  ambient air pollution and cigarette smoke: shape of the exposure-response
  relationships.
\newblock {\em Environmental Health Perspectives\/}~{\em 119}, 1616--1621.

\bibitem[\protect\citeauthoryear{Tonne, Melly, Mittleman, Coull, Goldberg, and
  Schwartz}{Tonne et~al.}{2007}]{Tonne07}
Tonne, C., S.~Melly, M.~Mittleman, B.~Coull, R.~Goldberg, and J.~Schwartz
  (2007).
\newblock A case-control analysis of exposure to traffic and acute myocardial
  infarction.
\newblock {\em Environmental Health Perspectives\/}~{\em 115}, 53--57.

\end{thebibliography}
