Alahmadi, H et al

*Hanan Alahmadi.

# Joint Quantile Disease Mapping with Application to Malaria and G6PD
Deficiency

Hanan Alahmadi  Håvard Rue  Janet van Niekerk  CEMSE Division, King Abdullah
University of Science and Technology, Kingdom of Saudi Arabia  King Saud
University, Kingdom of Saudi Arabia<EMAIL_ADDRESS>
(<day> <Month>, <year>; <day> <Month>, <year>; <day> <Month>, <year>)

###### Abstract

[Abstract]Statistical analysis based on quantile regression methods is more
comprehensive, flexible, and less sensitive to outliers when compared to mean
regression methods. When the link between different diseases are of interest,
joint disease mapping is useful for measuring directional correlation between
them. Most studies study this link through multiple correlated mean
regressions. In this paper we propose a joint quantile regression framework
for multiple diseases where different quantile levels can be considered. We
are motivated by the theorized link between the presence of Malaria and the
gene deficiency G6PD, where medical scientist have anecdotally discovered a
possible link between high levels of G6PD and lower than expected levels of
Malaria initially pointing towards the occurrence of G6PD inhibiting the
occurrence of Malaria. This link cannot be investigated with mean regressions
and thus the need for flexible joint quantile regression in a disease mapping
framework. Our joint quantile disease mapping model can be used for linear and
non-linear effects of covariates by stochastic splines, since we define it as
a latent Gaussian model. We perform Bayesian inference of this model using the
INLA framework embedded in the R software package INLA. Finally, we illustrate
the applicability of model by analyzing the malaria and G6PD deficiency
incidences in $21$ African countries using linked quantiles of different
levels.

###### keywords:

Bayesian analysis; Disease mapping; INLA; Joint quantile regression

††articletype: Article Type

##  1 Introduction

Malaria is considered a leading cause of mortality worldwide, and the disease
is most prominent in Africa. It has been estimated that malaria in 2017
affected about 219 million people and causing around 435,000 deaths 1. The
Malaria Atlas Project 2 provides a global database on malaria risk in order to
solve critical questions. This project disseminates free, accurate, and up-to-
date geographical on malaria and associated topics. One of their research
outputs points out a relationship between malaria and Glucose 6 phosphate
dehydrogenase (G6PD) deficiency, a genetic disorder that affects red blood
cells. The G6PD is a gene that provides instructions for making the
glucose-6-phosphate dehydrogenase enzyme. The research by the Malaria Atlas
Project found that G6PD deficiency is common in populations that have a high
level of malaria infection 3. Studies dating back to the early 1960s, 4, 5,
postulated that G6PD deficiency inhibits the occurrence of malaria. The
reasoning was that G6PD deficiency leads to the accumulation of oxygen
radicals inside the red blood cells ($H_{2}O_{2}$). This accumulation offers
resistance against malaria infection because the Plasmodium falciparum
parasite (the parasite that causes malaria) does not have any antioxidant
mechanism, which makes them more vulnerable to oxygen radicals 6, 7. The
hypothesis that G6PD deficiency provides some protection against Plasmodium
falciparum malaria was further supported by a review by Greene 8, published in
1993, based on experimental and population studies. At the same time, it was
acknowledged that there is not enough data in population studies, due to
limited sample sizes, to produce concluding evidence 9, 8. However, there are
opposing arguments, also based on limited population studies, stating that
G6PD deficiency by itself is unlikely to produce a significant protection
against malaria, see 10. In 1995, Ruwende et al. 11 suggest, from two case-
control studies of more than 2,000 African children, that G6PD deficiency
reduced the risk of severe malaria by around $50\%$. In 2017, a systematic
review by Mbanefo et al. 12 based on a selection of 28 various studies arrived
at that G6PD deficiency could potentially offer some protection against
uncomplicated malaria, but less likely so for severe malaria.

Following the results of 5, 9, it is of interest to perform a statistical
inference of such a relationship between diseases and quantify the
uncertainties involved. For this case, we propose using a quantile-based joint
model, instead of the standard joint that models correlation of the means of
the two diseases, since the G6PD deficiency may act as a resistance factor
against malaria, while not the other way around. Thus, to identify possible
directional correlation, this study looks at the joint quantiles between the
two diseases by modeling the high quantile of G6PD deficiency and the low
quantile of malaria. The joint quantile model can be applied to other disease
mapping problems. Quantile regression was introduced by Koenker and Basset
art0. After that the quantile regression has been widely used, in particular
for Bayesian spatial analysis 13. The R package bayesQR proposed by 14 can be
used to estimate the parameters in quantile regression using a Bayesian
approach with the asymmetric Laplace distribution. This package supports both
continuous-dependent and binary-dependent variables. In 15, the authors
proposed using a negative-binomial regression $\alpha$-quantiles approach with
an ecological regression model with application to disease mapping of lip
cancer.

The main difference to our work is that we consider joint quantile regression
with two diseases, instead of a single one. In joint quantile regression, one
can model spatial dependence through a Gaussian or t-copula process of the
quantile levels 16, which could provide certain benefits for cases with heavy-
tailed spatial data. One of the approaches to spatial quantile regression is
to use the Asymmetric Laplace Process (ALP) for modeling the data 17. However,
this assumes the data is coming from the ALP, regardless of the actual
generating distribution of the data. A quantile regression-based Bayesian
joint modeling analysis of longitudinal-survival data has been proposed in 18
and it extends the use of the asymmetric Laplace process as in 17 to joint
quantile regression. Markov chain Monte Carlo (MCMC) methods have been used
for parameter estimation in Bayesian quantile regression models, for example
in 19 for multivariate quantile regression. However, we advocate the use of
INLA over MCMC for practical disease mapping due to its computational
advantages. Spatial quantile regression is widely used with applications
ranging from modeling of wildfire risk 16 to studying healthy life years
expectancy 20 to economics 17. In 21, a Bayesian multiple quantile regression
method is proposed for linear models, and they used the working likelihood
instead of the likelihood of the generating distribution. In contrast, the
quantile regression in 22 was developed such that the likelihood of the
generating distribution is respected. For a comprehensive introduction to
quantile regression for spatial data, see 23, and for multivariate disease
mapping modeling we refer to 24, which includes many practical exercises and
examples, often provided with R-code implementation. As far as the authors are
aware, there is no available literature on joint quantile disease mapping is
available, which we aim to contribute in this study.

##  2 Disease Mapping

Disease mapping, also known as spatial epidemiology, analyzes the incidence of
disease using geographical information. In other words, Disease mapping
describe the spatial variation of disease. The two characteristics of disease
mapping are the location of the events, which is called spatial or
geographical distribution, and the disease.25. The Poisson distribution is
well representing the disease count for the data that have low disease count
for a relatively large population 25. For the region that consists of $n$ non-
overlapping areas 26, let $y_{i}$ denote the number of cases in regions $i$.
Often $y_{i}$ is assumed to be distributed as :

$\displaystyle y_{i}\sim\text{Poisson}(\mu_{i})$ (1)

where $\mu_{i}$ is the mean of $y_{i}$. The mean function often consists of
two components. The first component is usually called the relative risk, which
represents the risk within a region, it is unknown and the purpose of this
work to estimate these values. The second component is usually called
standardization, which represents the expected local count. The expected local
count is the value that represents our expectation if the population locally
behaved the way the standard population behaves. The expectation of the cases
in region $i$ can be written as follows:

$\displaystyle E(y_{i})=\mu_{i}=E_{i}\lambda_{i}$ (2)

where $E_{i}$ is the expected number for the $i$th area, which is usually
assumed to be a fixed quantity 25. The expected number can be obtained by
using indirect standardization as follows:

$\displaystyle E_{i}=\sum_{j=1}^{m}r_{j}^{(s)}n_{j}^{(i)}$ (3)

here, $r_{j}^{(s)}$ denotes the disease rate of the standard population in
stratum j, the rate is the number of cases divided by the population,
$n_{j}^{(i)}$ is the size in stratum j of area $i$ 27, and $\lambda_{i}$ is
the relative risk for $i$th area. Here $\lambda_{i}=1$ means there is no
augmented risk in comparison with the whole study area; $\lambda_{i}>1$ ,
$\lambda_{i}<1$ indicates higher risk and lower risk than the average
respectively 28. The maximum likelihood estimator of $\lambda_{i}$ is
$\hat{\lambda_{i}}=y_{i}/E_{i}$ which is correspond to the standardised
mortality ratio (SMR). However, mapping SMRs directly are misleading and
insufficient for counties with small populations. Therefore the covariates
need to be incorporated in order to smooth extreme values because of the small
sample sizes by borrowing information from neighboring counties. The model
considered in this work for disease mapping is formulated as follow:

$\displaystyle y_{i}\sim\text{Poisson}(E_{i}\lambda_{i}),i=1,...,n$ (4)
$\displaystyle\log(\lambda_{i})=\log(\eta_{i})=m_{0}+\sum_{f=1}^{F}\beta_{f}X_{if}+\sum_{r=1}^{R}\rho^{r}(u_{ir})+b_{i}$
(5)

where $\lambda_{i}$ is the mean of unit $i$, $m_{0}$ is the intercept that
follow a weakly informative Gaussian prior with mean zero and large variance,
$\sum_{f=1}^{F}\beta_{f}$ is the fixed effect of the covariates $X_{if}$.
Random effects such us splines for non linear effect of covariates
$\boldsymbol{u_{i}}$ is included through the functions
$\\{\rho^{r}\\}_{r=1}^{R}$, $\boldsymbol{b}$ is the spatial effects.
For the spatial effects, $\boldsymbol{b}$, different spatial models for areal
data can be assumed such as the Besag model 29, or the extended Besag-York-
Mollie model 29, the Leroux model 30, or the Dean’s model 31.

$\boldsymbol{x}=\\{\eta_{1},...,\eta_{n},m_{0},\beta_{1},...,\beta_{F},\boldsymbol{\rho},\boldsymbol{b}\\}$
is called a latent field. with hyperparameters
$\boldsymbol{\theta}=\\{\boldsymbol{\theta}_{\boldsymbol{\rho}},\boldsymbol{\theta}_{\boldsymbol{b}}\\}$,
then the data $\boldsymbol{y}$ is conditionally independent given the latent
field and the hyperparameters such that the likelihood function is

$\pi(\boldsymbol{x},\boldsymbol{\theta}|\boldsymbol{y})=\prod_{i=1}^{n}f(y_{i}|x_{i},\boldsymbol{\theta})$
(6)

###  2.1 Prior specification and posterior propriety

We assume prior independence amongst the parameters and as such we assign
Gaussian priors to the latent field elements and various other prior to the
hyperparameters as set out next.
For the latent field elements assume the following:

$\displaystyle m_{0}$ $\displaystyle\sim
N(0,\tau_{m}^{-1}),\quad\boldsymbol{\beta}|\tau_{\beta}\sim
N(\boldsymbol{0},\tau_{\beta}^{-1}\boldsymbol{I}),\quad$ (7)
$\displaystyle\boldsymbol{\rho}|\boldsymbol{\theta}_{\rho}$ $\displaystyle\sim
N(\boldsymbol{0},\boldsymbol{Q}^{-1}_{\rho}),\quad\boldsymbol{b}|\boldsymbol{\theta}_{b}\sim
N(\boldsymbol{0},\boldsymbol{Q}^{-1}_{b})\quad$

so that the joint prior for this part of the latent field is

$\boldsymbol{x}\sim N(\boldsymbol{0},\boldsymbol{Q}^{-1})$

where $\boldsymbol{Q}^{-1}$ has a block diagonal structure as formed from (7).
The vector of hyperparameters, $\boldsymbol{\theta}$ is assigned a joint prior
$\pi(\boldsymbol{\theta})$ which is composed of independent marginal proper
priors of any shape (not necessarily Gaussian).

The joint posterior of the unknown parameters, $\boldsymbol{x}$ and
$\boldsymbol{\theta}$ from (6) and (7) is

$\pi(\boldsymbol{x},\boldsymbol{\theta}|\boldsymbol{y})\propto\pi(\boldsymbol{y}|\boldsymbol{x},\boldsymbol{\theta})\pi(\boldsymbol{x}|\boldsymbol{\theta})\pi(\boldsymbol{\theta}),$

and based on the prior structures the posterior propriety holds.

##  3 Quantile Regression

Quantile regression describes the conditional quantile of the response
variable given the explanatory variables, instead of the conditional mean. Let
$Y$ be a real valued random variable. The $\alpha^{\text{th}}$ quantile of $Y$
is given by

$Q(\alpha)=F^{-1}(\alpha)=\inf\\{y:F(y)\geq\alpha\\}\quad\text{ for }\quad
0\leq\alpha\leq 1,$

where $F(y)=P(Y\leq y)$ is the Cumulative Distribution Function (CDF) of the
random variable $Y$. Like the mean regression, a loss function is used in
order to infer the parameters. The loss function of the quantile regression is
the check loss function. Given that $0\leq\alpha\leq 1$, $\forall
x\in\mathbb{R}$ the quantile loss function is defined as

$\rho_{\alpha}(x)=\left\\{\begin{array}[]{ll}x\alpha&x\geq 0\\\
x(\alpha-1)&x<0\end{array}\right.$

An estimate of the $\alpha^{\text{th}}$ quantile of the random variable $Y$
can be obtained by minimizing the following risk function:

$\underset{q_{\alpha}\in\mathbb{R}}{\operatorname{minimize}}\,\mathbb{E}\left[\rho_{\alpha}\left(Y-q_{\alpha}\right)\right]$
(8)

When $q_{\alpha}$ depends on the explanatory variables X, then it is called a
conditional quantile. The estimate of the conditional quantile is called
quantile regression. The quantile regression summarizes the relationship
between X and the quantile of Y. The estimate of the quantile regression can
be written as

$\hat{q}_{\alpha}=\underset{q_{\alpha}\in\mathbb{R}}{\operatorname{argmin}}\,\mathbb{E}\left[\rho_{\alpha}\left(Y-q_{\alpha}\right)\right]$
(9)

Then $F(\hat{q}_{\alpha})=\alpha$, where $F(y)$ is the CDF of the random
variable $Y$.

###  3.1 Model-based Quantile Regression

The goal of the statistical analysis based on the Bayesian methods is to make
inference from the posterior distribution for unknown parameters. Model-based
quantile regression is an approach for quantile regression that considers the
quantiles of the generating distribution proposed by 22. This approach extends
the Generalized Linear Mixed Model (GLMM) framework from modelling means to
modelling quantiles. Two steps can do this extension. The first step is
modeling the quantile; in this step, the quantile of the distribution is
linked to the linear predictor through an invertible function $g$. The second
step is mapping the quantile; in this step, the quantile is mapped to the
parameter of the distribution through a map function $h$. This approach can be
applied to both frequentist and Bayesian frameworks. The resulting parameters
of the Bayesian framework are all identifiable, making model-based quantile
regression appealing in the Bayesian inference. To see these steps, let
$F(y_{i};\lambda_{i})$ be the distribution of $Y_{i}|X_{i}$, where
$\lambda_{i}$ is the parameter of the distribution. Given $0\leq\alpha\leq 1$,
the $\alpha^{th}$ quantile of $Y_{i}|X_{i}$ is
$q_{i,\alpha}={Q}_{\alpha}(Y_{i}|X_{i})$. The two steps can be written as
follows:

$\textbf{Step 1 - Modelling}.$
The quantile $q_{i,\alpha}$ of the distribution $F(y_{i},\lambda_{i})$ is
modeled as follows:

$q_{i,\alpha}=g\left(\eta_{i,\alpha}\right),$

where $g$ is an invertible function, and $\eta_{i}^{\alpha}$ is the linear
predictor for the level $\alpha$ quantile for $i=1,\ldots,n$. The linear
predictor can include fixed effects, random effects, or both. Moreover,
parametric or semi-parametric models can be included in this approach in order
to study the impact of the covariates at different levels of the distribution
and non-parametric models can be used for prediction.

Step 2 - Mapping.
The quantile $q_{i,\alpha}$ is mapped to the parameter $\lambda_{i}$ of the
distribution $F(y_{i};\lambda_{i})$ as

$\lambda_{i}=h\left(q_{i,\alpha}\right),$ (10)

where $h$ is an invertible map function. The map $h$ can be obtained by two
steps. First, taking the inverse of the CDF ($F(y_{i};\lambda_{i})$) which
give you the quantile function $Q(\alpha,\lambda_{i})$. Then, we write the
parameter $\lambda_{i}$ as a function of the quantile, and that function is
the map $h$. In this approach, the parameter $\lambda_{i}$ is modeled
indirectly by the link between the quantiles of the generating distribution
and $\lambda_{i}$. Unlike mean regression, when the parameter of the
generating distribution links to the linear predictor through a function
$\lambda_{i}=g(\eta_{i})$, in model-based quantile regression the parameter of
the generating distribution is linked to the linear predictor through a
composition function $\lambda_{i}=h(g(\eta_{i}))$. In other words, in the mean
regression, the (GLMM) have a link function $\lambda_{i}=g(\eta_{i})$ to link
the parameter of the generating distribution to the linear predictor.

###  3.2 Model-based quantile regression for count data

The extension of model-based quantile regression for discrete random variables
is not straight-forward since the objective function in (9) is non-
differentiable for discrete random variables. The positive mass of the points
for the discrete variable prevent the sample quantile from having an
asymptotic distribution. Additionally, it is not easy to apply the modeling
and mapping steps of model-based quantile regression to discrete data. First,
in the modeling step, the common models for $g$ are the log for count data and
the logit for binary data, and they are continuous functions. Therefore, the
model $q_{i,\alpha}=g\left(\eta_{i,\alpha}\right)$ is not appropriate, since
the quantile which is on the left hand side is discrete whereas the function
$g$ is continuous. The second reason, in the mapping step, it is hard to get
the map $h$ because the CDF of the discrete is non-invertible, which implies
that there is no unique $\lambda_{i}$ to generate each quantile, as one can be
seen in Figure 1.

To address these issues, 22 approximated discrete distributions by continuous
counterparts, and then model the quantile for the continuous version instead
of the discrete. The continuous counterpart is obtained by interpolating the
cumulative distribution function (CDF) of the discrete random variable. The
model-based quantile method can be applied to discrete variables if their CDF
can be expressed as

$F_{X}(x;\lambda)=\mathbb{P}(X\leq x)=k(\lfloor x\rfloor,\lambda)$

where $k$ is a continuous function, and $X$ is a discrete random variable. The
interpolation can be obtained by removing the floor operator, so that
$k(x,\lambda)$ is the CDF of the continuous version of $X$, assigned
$X^{\prime}$. By definition of the floor, for all integers $x$

$F_{X}(x)=k(\lfloor x\rfloor,\lambda)=k(x,\lambda)=F_{X^{\prime}}(x).$

The continuous distribution of $X^{\prime}$ is considered as a continuous
generalization of the original variable because the two CDFs are equal for all
integer values $x$. The advantage of working with the continuous version of a
discrete distribution in the Bayesian framework is that a likelihood function
can be obtained by using the model-based quantile method, since the sample
quantiles for a discrete random variable are generally not asymptotically
normal 22.

Figure 1: (Top) The CDF of the discrete Poisson (dashed line), and the CDF of
the continuous Poisson (solid line). (Bottom) Quantile function of the
discrete (dashed line) and continuous (continuous line) Poisson distributions

22

####  3.2.1 Continuous Poisson

Here we present the details on the approximation of the discrete Poisson
distribution with a continuous Poisson counterpart.

The CDF of a Poisson distribution can be expressed as the ratio of an
incomplete and regular Gamma function as follows:

$X\sim\operatorname{Poisson}(\lambda)\quad F_{X}(x)=\mathbb{P}(X\leq
x)=\frac{\Gamma(\lfloor x\rfloor+1,\lambda)}{\Gamma(\lfloor x\rfloor+1)}\quad
x\geq 0$ (11)

where $\Gamma(x,\lambda)=\int_{\lambda}^{\infty}e^{-s}s^{x-1}\mathrm{~{}d}s$
is the upper incomplete Gamma function. Following Section 3.3, the Continuous
Poisson is then defined from (11) as

$X^{\prime}\sim\text{ Continuous Poisson }(\lambda)\quad
F_{X^{\prime}}(x)=\mathbb{P}\left(X^{\prime}\leq
x\right)=\frac{\Gamma(x+1,\lambda)}{\Gamma(x+1)}\quad x>-1$

The reason for changing the support from $x\geq 0$ to $x>-1$ is to avoid mass
at $0$, so there will be no jump on the CDF of the Continuous Poisson (CP)as
illustrated in Figure 1. If the support remains the same, then the value of
the CP will be $0$ if $x<0$ and about $0.4$ at $x=0$, which introduces a jump
at zero. However, if the support is $x>-1$, then there will be no jump at zero
because the CDF of the CP will be $0$ at $x=-1$, then an interpolation will be
applied from $x=-1$ to $x=0$. The Continuous and discrete Poisson random
variables can be related as $X=\left\lceil X^{\prime}\right\rceil$.

The model-based quantile regression model for Poisson data is then defined for
$Y_{i}\mid\eta_{i}$ a continuous Poisson random variable with parameter
$\lambda_{i}$ as

$\displaystyle q_{i,\alpha}$ $\displaystyle=$ $\displaystyle
g\left(\eta_{i,\alpha}\right)=\exp\left\\{\eta_{i,\alpha}\right\\}$
$\displaystyle\lambda_{i}$ $\displaystyle=$ $\displaystyle
h\left(q_{i,\alpha}\right)=\frac{\Gamma^{-1}\left(q_{i,\alpha}+1,1-\alpha\right)}{\Gamma\left(q_{i,\alpha}+1\right)}.$
(12)

###  3.3 Model-based quantile regression for disease mapping

From Sections 2 and 3.3 we can define a model-based quantile regression model
for disease mapping. One issue that remains is how to decompose the expected
number of cases into the local expectation, $E_{i}$ and the relative risk
$\lambda_{i}$. In the case of modeling the quantile instead of the mean there
are two options as discussed by 22:

  * •

Include $E_{i}$ in the linear model as an offset

$\displaystyle q_{i,\alpha}$ $\displaystyle=$
$\displaystyle\exp\left\\{\eta_{i,\alpha}+\log\left(E_{i}\right)\right\\}=E_{i}\exp\left\\{\eta_{i,\alpha}\right\\}$
$\displaystyle\lambda_{i}$ $\displaystyle=$
$\displaystyle\frac{\Gamma^{-1}\left(q_{i,\alpha}+1,1-\alpha\right)}{\Gamma\left(q_{i,\alpha}+1\right)}$
(13)

  * •

Consider it as a scaling of the parameter of the distribution

$\displaystyle q_{i,\alpha}$ $\displaystyle=$
$\displaystyle\exp\left\\{\eta_{i,\alpha}\right\\}$ $\displaystyle\lambda_{i}$
$\displaystyle=$ $\displaystyle
E_{i}\frac{\Gamma^{-1}\left(q_{i,\alpha}+1,1-\alpha\right)}{\Gamma\left(q_{i,\alpha}+1\right)}$
(14)

These two approaches are equivalent in the Poisson mean regression, but not
equal in the Poisson quantile regression and the choice of approach depends on
the purpose of the analysis. If the focus of the study is to infer a quantile-
specific model then (13) is more appropriate whereas (14) can be considered as
a model for the parameter $\lambda_{i}$.

##  4 Bayesian Joint Quantile Disease Mapping

The main goal of disease mapping is to estimate the relative risk of diseases
across regions. Sometimes specific diseases have similar spatial patterns due
to sharing the same etiologies. In this case, these diseases have some
dependence, and it would be more appropriate to model them jointly rather than
separately. Moreover, sometimes the dependence might be in different quantiles
between the diseases or some diseases could inhibit the occurence of another
disease. The proposed joint quantile disease mapping model links different
quantiles of multiple diseases by a more general framework by considering
dependence not in the mean, but in the quantiles.

###  4.1 Model specification

The joint quantile model for two diseases can be formulated as:

$\displaystyle y_{i1}$ $\displaystyle\sim$ $\displaystyle\text{ Poisson
}\left(\lambda_{i1}\right)$ $\displaystyle y_{i2}$ $\displaystyle\sim$
$\displaystyle\text{ Poisson }\left(\lambda_{i2}\right)$
$\displaystyle\log(q_{i1,\alpha_{1}})$ $\displaystyle=$
$\displaystyle\log(\eta_{i1,\alpha_{1}})=m_{1}+\sum_{f=1}^{F_{1}}\beta_{f}X_{if}+\sum_{r=1}^{R_{1}}\rho^{r}(u_{ir})+b_{i1}+S_{i}$
(15) $\displaystyle\log(q_{i2,\alpha_{2}})$ $\displaystyle=$
$\displaystyle\log(\eta_{i2,\alpha_{2}})=m_{2}+\sum_{f=1}^{F_{2}}\gamma_{f}Z_{if}+\sum_{r=1}^{R_{2}}\xi^{r}(v_{ir})+b_{i2}+c\,S_{i}$
(16)

where $\lambda_{ik}$ is the mean of unit $i$ for disease $k$ and it is mapped
to the $\alpha_{k}$ level quantile $q_{ik,\alpha_{k}}$ as in (3.2.1). In the
modeling part, $m_{k}$ is a disease-specific intercept that follows a weakly
informative Gaussian prior with zero mean and large variance, $b_{ik}$ is a
spatial random effect and $S_{i}$ is the shared spatial component. The model
also incorporates fixed effects of covariates $\boldsymbol{X}_{i}$ and
$\boldsymbol{Z}_{i}$ by $\sum_{f=1}^{F_{1}}\beta_{f}X_{if}$ in (15) and
$\sum_{f=1}^{F_{2}}\gamma_{f}Z_{if}$ in (16), respectively. Various random
effects such as splines for non-linear effects of covariates
$\boldsymbol{u}_{i}$ and $\boldsymbol{v}_{i}$ is included through functions
$\\{\rho^{r}\\}_{r=1}^{R_{1}}$ in (15) and $\\{\xi^{r}\\}_{r=1}^{R_{2}}$ in
(16), respectively.

For the disease-specific spatial effects, $\boldsymbol{b}_{k}$, various
spatial models for areal data can be assumed such as the Besag model 29, or
the extended Besag-York-Mollie model 29, the Leroux model 30, or the Dean’s
model 31.
The shared spatial component, $\boldsymbol{S}$, that links the two diseases
through their quantiles is assumed to be a Besag effects with precision matrix
$Q=\left(Q_{ij}\right)$, where for $j\neq i$,
$Q_{ii}=\tau\left(n_{i}+d\right)$, $Q_{ij}=-\tau$, and $n_{i}$ is the number
of neighbours of node $i$. The parameter $c\in\Re$ is used to scale the shared
component and correlate the two diseases in space.

Now we can collect
$m_{1},m_{2},\boldsymbol{\beta},\boldsymbol{\rho},\boldsymbol{b}_{1},\boldsymbol{S},\boldsymbol{\gamma},\boldsymbol{\xi},\boldsymbol{b}_{2}$
together with the linear predictors
$\eta_{11,\alpha_{1}},...,\eta_{n1,\alpha_{1}},\eta_{12,\alpha_{2}},...,\eta_{n2,\alpha_{2}}$
and form the latent field

$\boldsymbol{x}=\\{\eta_{11,\alpha_{1}},...,\eta_{n1,\alpha_{1}},\eta_{12,\alpha_{2}},...,\eta_{n2,\alpha_{2}},m_{1},m_{2},\beta_{1},...,\beta_{F_{1}},\boldsymbol{\rho},\boldsymbol{b}_{1},\boldsymbol{S},\boldsymbol{\gamma},\boldsymbol{\xi},\boldsymbol{b}_{2}\\}$

with hyperparameters
$\boldsymbol{\theta}=\\{c,\tau,d,\boldsymbol{\theta}_{\boldsymbol{\rho}},\boldsymbol{\theta}_{\boldsymbol{\xi}},\boldsymbol{\theta}_{\boldsymbol{b}_{1}},\boldsymbol{\theta}_{\boldsymbol{b}_{2}}\\}$,
then we have that the data
$\boldsymbol{y}=\\{\boldsymbol{y}_{1},\boldsymbol{y}_{2},\\}$ is conditionally
independent given the latent field and the hyperparameters such that the
likelihood function is

$\pi(\boldsymbol{x},\boldsymbol{\theta}|\boldsymbol{y})=\prod_{i=1}^{n}f(y_{i}|x_{i},\boldsymbol{\theta})$
(17)

###  4.2 Prior specification and posterior propriety

We assume prior independence amongst the parameters and as such we assign
Gaussian priors to the latent field elements and various other prior to the
hyperparameters as set out next.
For the latent field elements assume the following:

$\displaystyle m_{k}$ $\displaystyle\sim$ $\displaystyle
N(0,\tau_{m}^{-1}),\quad\boldsymbol{\beta}|\tau_{\beta}\sim
N(\boldsymbol{0},\tau_{\beta}^{-1}\boldsymbol{I}),\quad\boldsymbol{\gamma}|\tau_{\gamma}\sim
N(\boldsymbol{0},\tau_{\gamma}^{-1}\boldsymbol{I})$
$\displaystyle\boldsymbol{\rho}|\boldsymbol{\theta}_{\rho}$
$\displaystyle\sim$ $\displaystyle
N(\boldsymbol{0},\boldsymbol{Q}^{-1}_{\rho}),\quad\boldsymbol{\xi}|\boldsymbol{\theta}_{\xi}\sim
N(\boldsymbol{0},\boldsymbol{Q}^{-1}_{\xi}),\quad\boldsymbol{S}|\tau,d\sim
N(\boldsymbol{0},\boldsymbol{Q}^{-1}_{S})$
$\displaystyle\boldsymbol{b}_{1}|\boldsymbol{\theta}_{b_{1}}$
$\displaystyle\sim$ $\displaystyle
N(\boldsymbol{0},\boldsymbol{Q}^{-1}_{b_{1}}),\quad\boldsymbol{b}_{2}|\boldsymbol{\theta}_{b_{2}}\sim
N(\boldsymbol{0},\boldsymbol{Q}^{-1}_{b_{2}})$ (18)

so that the joint prior for these elements of the latent field is

$\boldsymbol{x}\sim N(\boldsymbol{0},\boldsymbol{Q}^{-1})$

where $\boldsymbol{Q}^{-1}$ has a block diagonal structure as formed from
(18).
The vector of hyperparameters, $\boldsymbol{\theta}$ is assigned a joint prior
$\pi(\boldsymbol{\theta})$ which is composed of independent marginal proper
priors of any shape (not necessarily Gaussian).

The shared spatial field is assumed to follow Besag model but with an
additional parameter $d$, to ensure a proper prior of $\boldsymbol{S}$ as
follows:

$\displaystyle\boldsymbol{S}|\tau,d$ $\displaystyle\sim$ $\displaystyle
N(\boldsymbol{0},\boldsymbol{Q}^{-1}_{S})$

with the entries of $\boldsymbol{Q}_{S}$ as follows:

$Q_{S,ii}=\tau(n_{i}+d)\quad\text{and}\quad Q_{S,ij}=-\tau,$ (19)

for $i\neq j$, and $j$ is in the neighbourhood of $i$.

The disease-specific spatial fields $\boldsymbol{b}_{1}$ and
$\boldsymbol{b}_{2}$ are assumed to follow uncorrelated BYM/CAR models where
we reparameterize the precision matrix similar to 32 to have more orthogonal
parameters resulting in useful practical interpretation of the weight
parameter, $\phi$. One issue with the proper CAR parameterization proposed by
32 is that the weight parameter is still not practically a weight since the
unstructered effect and the Besag field might have different generalized
variances. To alleviate this issue, we scale both the unstructured and Besag
components to have the same geometric mean and define the proper scaled BYM
field as

$\boldsymbol{b}_{k}=\frac{1}{\sqrt{\tau_{b_{k}}}}\left(\sqrt{1-\phi_{b_{k}}}\boldsymbol{b}_{k}^{1}+\sqrt{\phi_{b_{k}}}\boldsymbol{b}_{k}^{2}\right)$
(20)

with $\boldsymbol{b}_{k}^{1}$ a scaled IID effect and $\boldsymbol{b}_{k}^{2}$
a scaled Besag effect as in (19). With this formulation, $\phi_{b_{k}}$ can be
interpreted as the proportion of the marginal variance explained by the
spatial effect, and $1-\phi_{b_{k}}$ is the proportion of the marginal
variance explained by the unstructured effect.

The joint posterior of the unknown parameters, $\boldsymbol{x}$ and
$\boldsymbol{\theta}$ from (17) and (18) is

$\pi(\boldsymbol{x},\boldsymbol{\theta}|\boldsymbol{y})\propto\pi(\boldsymbol{y}|\boldsymbol{x},\boldsymbol{\theta})\pi(\boldsymbol{x}|\boldsymbol{\theta})\pi(\boldsymbol{\theta}),$

and based on the prior structures the posterior propriety holds.

###  4.3 Approximate inference using INLA

Computational Bayesian inference can be achieved largely in one of two ways,
either through sampling-based methods like Markov Chain Monte Carlo (MCMC) and
deviants or approximately using approximate methods like Variational methods
or Laplace approximations like Integrated Nested Laplace Approximation (INLA).
INLA, as introduced by 33, has been shown to be widely applicable to various
statistical models; in particular, to the latent Gaussian models class of
which disease mapping models are included 34, 35, 36, 37

INLA employs a series of Laplace approximations and numerical integration to
perform approximate Bayesian inference through numerically approximating the
posterior densities of the latent field and hyperparameters. For data
$\boldsymbol{y}$, latent field $\boldsymbol{x}$ and hyperparameters
$\boldsymbol{\theta}$, INLA can be summarized as follows:

  1. 1.

Find the $m$-variate Gaussian approximation of
$\pi(\boldsymbol{x}|\boldsymbol{\theta},\boldsymbol{y})$ at the mode
$\boldsymbol{\mu}(\boldsymbol{\theta})$ with matching curvature using the
Hessian of $\pi(\boldsymbol{x}|\boldsymbol{\theta},\boldsymbol{y})$ at the
mode $\boldsymbol{\mu}(\boldsymbol{\theta})$.

  2. 2.

Let

$\tilde{\pi}(\boldsymbol{\theta}|\boldsymbol{y})\propto\frac{\pi(\boldsymbol{x}^{*},\boldsymbol{\theta}|\boldsymbol{y})}{\pi_{G}(\boldsymbol{x}^{*}|\boldsymbol{\theta},\boldsymbol{y})}|_{\boldsymbol{x}^{*}=\boldsymbol{\mu}(\boldsymbol{\theta})}$
(21)

and locate the mode of $\tilde{\pi}(\boldsymbol{\theta}|\boldsymbol{y})$ and
find a set of integration points $\boldsymbol{\theta}_{k},k=1,2,...,T$ in the
area of the highest probability mass.

  3. 3.

Calculate

$\tilde{\pi}(\theta_{j}|\boldsymbol{y})=\int_{\boldsymbol{\theta}_{-j}}\tilde{\pi}(\boldsymbol{\theta}|\boldsymbol{y})d\boldsymbol{\theta}_{-j}$
(22)

where we note that this is a low-dimensional integral since $p$ is generally
small.

  4. 4.

Now define

$\displaystyle\tilde{\pi}(x_{i}|\boldsymbol{\theta}_{k},\boldsymbol{y})\approx\frac{\pi(\boldsymbol{x}^{*},\boldsymbol{\theta}_{k}|\boldsymbol{y})}{\pi_{G}(\boldsymbol{x}^{*}_{-i}|x_{i},\boldsymbol{\theta}_{k},\boldsymbol{y})}|_{\boldsymbol{x}^{*}_{-i}=\boldsymbol{\mu}_{-i}(\boldsymbol{\theta}_{k})}$
(23)

with
$\pi_{G}(\boldsymbol{x}^{*}_{-i}|x_{i},\boldsymbol{\theta},\boldsymbol{y})$
the $(m-1)$-variate Gaussian approximation at the mode
$\boldsymbol{\mu}_{-i}(\boldsymbol{\theta})$ for the $T$ configuration points
$\boldsymbol{\theta}_{k},k=1,2,...,T$, and calculate

$\tilde{\pi}(x_{i}|\boldsymbol{y})\approx\sum_{k=1}^{T}\tilde{\pi}(x_{i}|\boldsymbol{\theta}_{k},\boldsymbol{y})\tilde{\pi}(\theta_{k}|\boldsymbol{y})\Delta_{k}$
(24)

where $\tilde{\pi}(\theta_{k}|\boldsymbol{y})$ is from step 3, with
$\Delta_{k}$ the step size.

Various simplifications to the approximations have been proposed as well in
order to achieve increased computational efficiency such as an empirical Bayes
approach where the integration points $\boldsymbol{\theta}_{k}$ are all set to
the mode of $\tilde{\pi}(\boldsymbol{\theta}|\boldsymbol{y})$, which is named
a Simplified Laplace approximation strategy.

###  4.4 Simulation Study

The code for this simulation study is available at
https://github.com/JanetVN1201/Code_for_papers/tree/main/Joint%20quantile%20disease%20mapping%20.

In this part, simulated independent and correlated data was added to the
Pennsylvania map, which is considered as a connected graph of size $67$.
Figure 2 shows one realization of the correlated data that was added to the
Pennsylvania map.

Figure 2:  One realization of cases of disease 1 (left) and disease 2 (right).

The correlated data was generated as follows.

$\displaystyle y_{i1}$ $\displaystyle\sim$ $\displaystyle\text{ Poisson
}\left(\lambda_{i1}\right)\quad\text{with}\quad\log(q_{i1,0.2})=1+S_{i}$
$\displaystyle y_{i2}$ $\displaystyle\sim$ $\displaystyle\text{ Poisson
}\left(\lambda_{i2}\right)\quad\text{with}\quad\log(q_{i2,0.8})=1+0.7\,S_{i}$

where $S_{i}$ is the shared component that follow a Besag proper model (19)
with precision matrix $Q$, where $Q_{ii}=\left(n_{i}+1\right)$ and $Q_{ij}=-1$
for $j\neq i$. In other words, the extra term added on the diagonal is $d=1$,
and the precision parameter $\tau=1$.
There are two hyperparameters for this model: a parameter $d$ for controlling
the properness, and a parameter $\tau>0$ that is a scaling parameter. The
hyperpriors defined on the $\log$ scale of these hyperparameters are as
follows.

$\displaystyle\log(d)$ $\displaystyle\sim$ $\displaystyle\text{loggamma}(1,1)$
$\displaystyle\log(\tau)$ $\displaystyle\sim$
$\displaystyle\text{loggamma}(1,5e-04)$

The estimated values of the parameters obtained by R-INLA are similar to the
true values as seen in Tables 1 and 2, and Figures 3 and 4.

| mean | sd | 0.025quant | 0.975quant
---|---|---|---|---
$m_{1}$ | 1.137 | 0.042 | 1.053 | 1.218
$m_{2}$ | 1.003 | 0.039 | 0.923 | 1.079
Table 1: The estimated intercepts and the 95% credible intervals. Figure 3: The posterior distribution for the intercepts | mean | 0.025quant | 0.975quant | mode
---|---|---|---|---
$\tau$ | 1.337 | 0.77 | 2.074 | 1.257
$d$ | 1.409 | 0.472 | 3.612 | 0.887
$c$ | 0.838 | 0.636 | 1.034 | 0.845
Table 2: The estimated hyperparameters and the 95% credible intervals.

Figure 4: The posterior distribution for the hyperparameters

Model selection criteria are presented in Tables 3 and 4 and it shows a
preference for the joint quantile model when the data are correlated and a
preference for the separate models, which are (15) and (16) without the shared
components, when the data are independent. This indicates stable estimation
and the model’s ability to distinctly estimate an associated joint model if
needed.

| DIC | WAIC
---|---|---
Separate 1 | 3311 | 3350
Separate 2 | 2439 | 2452
Sum of Separates | 5750 | 5802
Joint quantile | 5642 | 5660
Table 3: The DIC and WAIC for correlated data. | DIC | WAIC
---|---|---
Separate 1 | 2951 | 2967
Separate 2 | 2925 | 2953
Sum of Separates | 5876 | 5919
Joint quantile | 6002 | 6100
Table 4: The DIC and WAIC for independent data.

##  5 Joint quantile disease mapping model for Malaria and G6PD

In this section, we fit the Bayesian model by using R-INLA to estimate the
risks of malaria and G6PD deficiency in some African countries by using
separate and joint quantile models. The code for this analysis is available at
https://github.com/JanetVN1201/Code_for_papers/tree/main/Joint%20quantile%20disease%20mapping%20.

###  5.1 Exploratory data analysis

The malaria cases and G6PD cases per region is obtained from
https://malariaatlas.org/. Various country-level covariates can be used in our
model but for the motivating example the emphasis is placed on the joint
component, even though for a thorough analysis of the data itself, various
fixed and random effects might be considered.

We selected only the countries for which information for both Malaria and G6PD
is available, as indicated in Figure 5.

Figure 5:  Countries where both G6PD deficiency and malaria cases exist (left)
on the African continent (right)

According to Figure 5, the countries are distributed around the world. Since
we want to investigate the spatial correlation we consider the African
continent so that most countries have some neighbours as in Figure 5. In
Figure 6, 6(b) and 6(e) the SMRs for malaria and G6PD deficiency. It can be
seen that, in general, the risk of G6PD deficiency is higher than the risk of
malaria because G6PD deficiency has a higher SMR. Some countries like Abidjan
and Madagascar that considered to have the highest risk of G6PD deficiency,
they have the lowest risk of malaria according to the SMR values, which could
indicate a prohibitive relationship between these two diseases. The observed
cases for malaria and G6PD can be seen in Figures 6(a) and 6(d) respectively.
Kenya has the highest number of malaria cases, while Nigeria has the highest
number of G6PD deficiency cases.

((a)) Observed cases

((b)) SMR

((c)) Expected cases

((d)) Observed cases

((e)) SMR

((f)) Expected cases

Figure 6: Summaries from exploratory data analysis for Malaria (top row) and
G6PD (bottom row)

###  5.2 Results

It is believed that G6PD deficiency limits the occurrence of malaria 9, 5. We
expect a correlation between a high quantile of G6PD and a low quantile of
malaria. Therefore we applied the joint quantile model that was discussed in
section 4 where $y_{i1}$ and $y_{i2}$ represent the cases of malaria and G6PD
deficiency, respectively. The quantile levels are $\alpha_{1}=0.2$ and
$\alpha_{2}=0.8$.

Table 5 shows the overall means for the malaria and the G6PD deficiency. Table
6 shows the estimation of the hyperparameters. The precision of the random
effects indicate that most of the spatial variability comes from the shared
component. There is a significant correlation between a high quantile of G6PD
deficiency and a low quantile of malaria as observed from the point estimate
and the credible interval of $c$, which is the coefficient of the shared
component. This finding is consistent with the studies 9, 5.

In contrast, Table 7 shows the estimation of the hyperparameters for a high
quantile of malaria with a low quantile of G6PD deficiency. As can be seen
based on the credible interval of $c$, there is no significant correlation
between these two quantiles. This is expected because having G6PD deficiency
protects you from having malaria, but malaria does not influence G6PD
deficiency.

Figure 10 presents maps of the spatial effects. Figures 10(a) and 10(d) show
the shared spatial effect for malaria and G6PD deficiency, respectively. The
structure of the effects appear similar. However, the shared spatial effect
for the G6PD deficiency is much lower. This is expected because the value of
$c$ is smaller than 1. The disease-specific spatial effects in Figures 10(b)
and 10(e) are very low compared to the shared spatial effect. This can also be
seen from the posterior precision estimates for $b_{1}$ and $b_{2}$ compared
to that of the shared component in Table 6. Because the disease-specific
spatial effects are smaller than the shared spatial effect, the total spatial
effect (that is, the sum of the shared and the specific-spatial effects) is
very similar to those for the shared effects, see Figures 10(c) and 10(f). The
iid random effect for G6PD deficiency, see Figures 11(b) and 11, is higher
than the iid for malaria. This is understandable because the value of Phi for
malaria is bigger than the one for G6PD deficiency, which means G6PD
deficiency accounts for more iid effect.

The relative risks in Figures 9(a) and 9(b) show similar relative risks as
obtained by the separate models. However, observe that the relative risk for
G6PD deficiency from the joint model is higher than the one from the separate.
This difference between the relative risks is due to borrowing strength from
the spatial pattern of malaria through the shared component. The joint
quantile model predicts the cases well for both diseases, as can be seen in
Figures 11(c) and 11(d). The model comparison shows a preference for the joint
quantile model over the separate models because the values of DIC and WAIC for
the joint quantile are less than the sum of the tests for the separate models.

| mean | sd | 0.025quant | 0.975quant
---|---|---|---|---
$m_{1}$ | 7.852 | 0.634 | 6.58 | 9.12
$m_{2}$ | 4.245 | 0.285 | 3.669 | 4.81
Table 5: The estimated intercepts and the 95 % credible intervals. Figure 7:
The posterior distributions for the intercepts

The estimation results of the hyperparameters are given in Tables 6 and 7.

| mean | 0.025quant | 0.975quant | mode
---|---|---|---|---
$\tau$ | 0.107 | 0.044 | 0.224 | 0.082
$d$ | 1.665 | 0.379 | 4.434 | 0.953
$\tau_{b_{1}}$ | 39.1 | 0.831 | 244.1 | 1.862
$\phi_{b_{1}}$ | 0.321 | 0.021 | 0.827 | 0.056
$\tau_{b_{2}}$ | 1.226 | 0.625 | 2.091 | 1.096
$\phi_{b_{2}}$ | 0.186 | 0.008 | 0.658 | 0.019
$c$ | 0.291 | 0.066 | 0.521 | 0.285
Table 6: The estimation of the hyperparameters for low quantile of malaria and high quantile of G6PD deficiency. | mean | 0.025quant | 0.975quant | mode
---|---|---|---|---
$\tau$ | 3.648 | 1.437 | 7.976 | 2.688
$d$ | 1.819 | 0.423 | 5.258 | 0.977
$\tau_{b_{1}}$ | 11.31 | 5.485 | 24.16 | 8.088
$\phi_{b_{1}}$ | 0.056 | 0 | 0.311 | 0
$\tau_{b_{2}}$ | 131.5 | 55.62 | 326.2 | 82.55
$\phi_{b_{2}}$ | 0.14 | 0 | 0.853 | 0
$c$ | 0.002 | -0.008 | 0.015 | -0.002
Table 7: The estimation of the hyperparameters for high quantile of malaria
and low quantile of G6PD deficiency.

The values for the model choice criteria, DIC and WAIC, are given in Table 8.

| DIC | WAIC
---|---|---
G6PD Deficiency | 168 | 164.8
Malaria | 246.8 | 241.4
Sum | 414.8 | 406.2
Joint quantile | 413.6 | 402.2
Table 8: Model choice criteria.

Figure 8: The posterior distributions for the hyperparameters

((a)) The relative risk estimate for malaria

((b)) The relative risk estimate for G6PD

((c)) The relative risk estimate for malaria

((d)) The relative risk estimate for G6PD

Figure 9: The estimates from the joint model (top), and the separate models
(bottom)

((a)) Shared effect

((b)) specific effect

((c)) The total effect

((d)) Shared effect

((e)) specific effect

((f)) The total effect

Figure 10: Posterior estimates of the spatial effects from the joint model for
malaria (top) and for G6PD deficiency (bottom)

((a)) The iid effect for Malaria

((b)) The iid effect for G6PD

((c)) The predicted cases for Malaria

((d)) The predicted cases for G6PD

Figure 11: The iid effect and predicted cases

##  6 Concluding Remarks

The motivation stemmed from estimating the relative risk of Malaria and G6PD
deficiency, jointly, on the African continent. The G6PD deficiency is
considered as a resistance against malaria based on anecdotal medical studies
(see 9 and 5). In this case, joint mean disease mapping will not provide the
information needed to investigate these initial findings. Therefore, we
considered a joint quantile disease mapping of different quantiles for the
diseases. The approach is successful since considering the joint quantile
model allows a possible investigation of the correlation between any level of
the conditional distributions of the random variables that represents the
number of cases, not only the correlation between their means. An advantage of
the proposed approach is that the computationally efficient INLA method is
used for statistical inference, such as estimating the relative risk.

Our main contribution is two-fold. Firstly, we propose a very general joint
quantile disease mapping model where the correlation between different
quantiles can be inferred and multiple diseases can be considered, together
with an efficient computational framework for the inference thereof. Secondly,
the significant correlation between a high quantile of G6PD cases and a low
quantile of Malaria cases encourages further investigation based on expanded
data collection efforts as already underway at the Malaria Atlas Project. This
analysis provides a solid statistical framework to the anecdotal findings as
remarked by medical professionals, and could underpin future studies in this
direction.

## References

  * 1 Talapko J, Škrlec I, Alebić T, Jukić M, Včev A. Malaria: the past and the present. Microorganisms 2019; 7(6): 179.
  * 2 Hay SI, Snow RW. The Malaria Atlas Project: developing global maps of malaria risk. PLoS medicine 2006; 3(12): e473–e473.  doi: 10.1371/journal.pmed.0030473
  * 3 Howes RE, Piel FB, Patil AP, et al. G6PD deficiency prevalence and estimates of affected populations in malaria endemic countries: a geostatistical model-based map. PLoS Medicine 2012; 9(11): e1001339.
  * 4 Allison A. Glucose-6-Phosphate Dehydrogenase Deficiency in Red Blood Cells of East Africans. Nature 1960; 186(4724): 531–532.  doi: 10.1038/186531a0
  * 5 Allison A, Clyde D. Malaria in African children with deficient erythrocyte glucose-6-phosphate dehydrogenase. British Medical Journal 1961; 1(5236): 1346.
  * 6 Aziz MSS. Biochemistry, Hexose Monophosphate Pathway. StatPearls Publishing 2021\.
  * 7 Laslett N, Hibbs J, Hallett M, Ghaneie A, Zemba-Palko V. Glucose-6-Phosphate Dehydrogenase Deficiency-Associated Hemolytic Anemia and Methemoglobinemia in a Patient Treated With Hydroxychloroquine in the Era of COVID-19. Cureus 2021; 13(5).
  * 8 Greene LS. G6PD deficiency as protection against falciparum malaria: An epidemiologic critique of population and experimental studies. American Journal of Physical Anthropology 1993; 36(S17): 153-178.
  * 9 Beutler E. G6PD deficiency. Blood 1994; 84(11): 3613–3636.
  * 10 Martin S, Miller L, Alling D, et al. Severe malaria and glucose-6-phosphate-dehydrogenase deficiency: a reappraisal of the malaria/G-6-PD hypothesis. The Lancet 1979; 313(8115): 524–526.
  * 11 Ruwende C, Khoo SC, Snow RW, et al. Natural selection of hemi- and heterozygotes for G6PD deficiency in Africa by resistance to severe malaria. Nature 1995; 376(6537): 246–249.  doi: 10.1038/376246a0
  * 12 Mbanefo EC, Ahmed AM, Titouna A, et al. Association of glucose-6-phosphate dehydrogenase deficiency and malaria: a systematic review and meta-analysis. Scientific reports 2017; 7: 45963–45963.  doi: 10.1038/srep45963
  * 13 Reich BJ, Fuentes M, Dunson DB. Bayesian Spatial Quantile Regression. Journal of the American Statistical Association 2011; 106(493): 6-20.
  * 14 Benoit DF, Poel V. dD. bayesQR: A Bayesian Approach to Quantile Regression. Journal of Statistical Software 2017; 76(7): 1–32.
  * 15 Chambers R, Dreassi E, Salvati N. Disease mapping via negative binomial regression M-quantiles. Statistics in Medicine 2014; 33(27): 4805-4824.
  * 16 Chen X, Tokdar ST. Joint quantile regression for spatial data. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 2021; 83(4): 826-852.
  * 17 Lum K, Gelfand AE. Spatial Quantile Multiple Regression Using the Asymmetric Laplace Process. Bayesian Analysis 2012; 7(2): 235 – 258\.
  * 18 Zhang H, Huang Y. Quantile regression-based Bayesian joint modeling analysis of longitudinal–survival data, with application to an AIDS cohort study. Lifetime Data Analysis 2020; 26(2): 339–368.
  * 19 Santos B, Kneib T. Noncrossing structured additive multiple-output Bayesian quantile regression models. Statistics and Computing 2020; 30(4): 855–869.
  * 20 Trzpiot G, Orwat-Acedańska A. Spatial quantile regression in analysis of healthy life years in the European Union countries. Comparative Economic Research 2016; 19(5): 179-199.
  * 21 Wu T, Narisetty NN. Bayesian Multiple Quantile Regression for Linear Models Using a Score Likelihood. Bayesian Analysis 2021; 1(1): 1–29.
  * 22 Padellini T, Rue H. Model-based Quantile Regression for Discrete Data. 2018.  arXiv no: 1804.03714.
  * 23 McMillen DP. Quantile Regression for Spatial Data.  Springer Science & Business Media .  2013\.
  * 24 Botella-Rocamora P. Disease Mapping: From Foundations to Multidimensional Modeling.  Chapman and Hall/CRC .  2019\.
  * 25 Lawson AB. Bayesian Disease Mapping: Hierarchical Modeling in Spatial Epidemiology.  CRC press .  2018\.
  * 26 Riebler A, Sørbye SH, Simpson D, Rue H. An intuitive Bayesian spatial model for disease mapping that accounts for scaling. Statistical Methods in Medical Research 2016; 25(4): 1145–1165.
  * 27 Moraga P. Small Area Disease Risk Estimation and Visualization Using R. R Journal 2018; 10(1): 495.
  * 28 Blangiardo M, Cameletti M, Baio G, Rue H. Spatial and spatio-temporal models with R-INLA. Spatial and Spatio-Temporal Epidemiology 2013; 4: 33–49.
  * 29 Besag J, York J, Mollié A. Bayesian image restoration, with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics 1991; 43(1): 1–20.
  * 30 Leroux BG, Lei X, Breslow N. Estimation of disease rates in small areas: a new mixed model for spatial dependence. In: Springer. 2000 (pp. 179–191).
  * 31 Dean C, Ugarte M, Militino A. Detecting interaction between random region and fixed age effects in disease mapping. Biometrics 2001; 57(1): 197–202.
  * 32 Banerjee S, Carlin BP, Gelfand AE. Hierarchical modeling and analysis for spatial data.  CRC press .  2014\.
  * 33 Rue H, Martino S, Chopin N. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 2009; 71(2): 319–392.
  * 34 Martínez-Minaya J, Conesa D, López-Quílez A, Vicent A. Climatic distribution of citrus black spot caused by Phyllosticta citricarpa. A historical analysis of disease spread in South Africa. European journal of plant pathology 2015; 143(1): 69–83.
  * 35 Moraga P. Geospatial health data: Modeling and visualization with R-INLA and shiny.  Chapman and Hall/CRC .  2019\.
  * 36 Moraga P, Dean C, Inoue J, Morawiecki P, Noureen SR, Wang F. Bayesian spatial modelling of geostatistical data using INLA and SPDE methods: A case study predicting malaria risk in Mozambique. Spatial and Spatio-temporal Epidemiology 2021; 39: 100440.
  * 37 Ugarte MD, Adin A, Goicoa T, Militino AF. On fitting spatio-temporal disease mapping models using approximate Bayesian inference. Statistical methods in medical research 2014; 23(6): 507–530.