# A Bayesian Bivariate Model for Spatially Correlated Binary Outcomes

Authors: Thierno Souleymane Barry; Oscar Ngesa; Nelson Owuor Onyango; Henry Mwambi
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Diseases have been studied separately, but two diseases have inherent dependencies on each other, modelling them separately negates practical reality. The authors’ modelling processes are based on univariate separate regressions, which connect each illness to covariates separately. Therefore, the focus of this article is to estimate the spatial correlation within geographic regions using latent variables. Individual and areal-level information, as well as spatially dependent random effects for each spatial unit, are incorporated into the models developed using a hierarchical structure. Simulation techniques provide to assess the models’ performance using Bayesian computing approaches (INLA and MCMC). The findings show a reasonable performance of the DIC and RMSE values of the proposed latent model. From that, the model can be considered as the best compared to the shared component model, multivariate conditional autoregressive model, and univariate models.


---

## Body

## 1. Introduction

The primary outcomes of many medical studies are binary. For example, the investigation of the prevalence of anemia in children under the age of five. Such an outcome can only take two values: 1 or 0, indicating the presence or absence of the disease, respectively. The logistic regression model is used in that framework to handle the effect of several explanatory variables [1, 2].The Multiple binary logistic regression is a more objective approach to study diseases effect. The goal is to determine whether or not there is a relationship between the observed phenomena. Even better, it is attempting to determine whether one phenomenon has a significant impact on the other. According to Davenport et al. [3], in the context of epidemiological studies, a disease complexity is often not adequately characterized by a single outcome, and several aspects of the patient’s response, either to a treatment or to a risk factor, must be considered.The mapping of disease prevalence and incidence is extremely important in the context of health policy. Its main goal is to smooth and predict some such health outcomes over a specific geographic domain. Areal data analysis can help in the investigation for scientific explanations. Furthermore, it aids in even general problem solving because observations in geographic areas are much more than often correlated [4]. Spatial modelling takes into account the spatial relationships and interactions between the diseases or diseases’ factors, and it aids in the identification and understanding of the mechanisms and processes underlying the phenomena. The CAR model must be used to account for spatial autocorrelation in classification and spatial regression problems. Since the surrounding areas have similar social, economic, and cultural characteristics, the effect may occur in social research [5]. Random effects are included in the nonspatial component to account for heterogeneity, overdispersion, and arbitrary choice of a spatial unit. These random effects would then correct and smooth the distribution [6].Due to the rapid development of computational techniques, the use of spatial models in a variety of fields has increased dramatically in recent years. Several authors have developed and applied spatial modelling in practical situations. To model the population and socioeconomic data, Zhang et al. [7] applied a Bayesian conditional autoregressive model to examine the association of meteorological and socioeconomic factors with Japanese encephalitis. Illian et al. [8] fitted complex spatial point process models using integrated nested Laplace approximation (INLA). In the framework of joint diseases modelling, Ngesa et al. [9] developed a model with HIV and HSV-2 applications to capture the covariance spatial random effects. A structured additive distributional biprobit model was fitted by Gayawan et al. [10]. The model developed was mainly focused on the spatial variability in the levels of, as well as the correlation between two outcomes in multiple West African countries. Ibanez et al. [11] used shared component models for joint spatial modelling in a Bayesian inference approach to identify shared patterns among chronic related potentially preventable hospitalizations. A multivariate Bayesian spatial modelling approach is used by Adeyemi et al. [12] to jointly model the counts of armed robbery and theft (stealing).This paper develops an appropriate multivariate approach, namely, a bivariate latent model, for spatial modelling while also accounting for the dependence of two outcomes via the possible linear correlation they may share. To compare the models, independent univariate models and Multivariate CAR model analysis (MCAR), as well as a shared component model, were developed. The development of a latent model is a critical component of our research.The article is organized as follows: Section2 provides a graphical representation of correlation coefficients on two-dimensional bivariate outcomes and the specification of the different models. Section 3 describes the simulation design. Section 4 covers the findings and discussions. Finally, Section 5 concludes the article by outlining future research directions.


## 2. Methods

In this study, we developed a new model for multivariate spatial data recorded on a lattice. The proposed latent model estimates the spatial correlation within geographical regions. To illustrate this, we considered a lattice with only three regions and two binary outcomes (Figure1). The regions are represented by the circles, and for each region, two binary outcomes are recorded Y1 and Y2. Three different correlation parameters are indicated in the figure. Spatial autocorrelation parameters α1 and α2 are related to the observations of variables Y1 and Y2 across the regions on the lattice. There are also α3 and α4 coefficients, which link variable Y1 in region i with variable Y2 in region j and variable Y1 in region j with variable Y2 in region i, respectively, (spatial cross-correlations). There is a possibility that the two variables Y1 and Y2, which are recorded in each region, are correlated, and the correlation coefficients are ρ1, ρ2, and ρ3. The latter correlations are referred to as within correlation coefficients in this study.Figure 1
A representation of correlation coefficients on two-dimensional bivariate binary outcomes.


### 2.1. Specification of the Models

LetYijk be the disease status of individual i for disease k in region j. Especially, k=1 for disease 1 and k=2 for disease 2, j=1,2,…,p and i=1,2,…,nj, where nj is the number of individuals in region j. The response variables are binary and we have the following equation:(1)Yijk=1,if individualiinregionjis positive for diseasek,0,otherwise.Yijk is assumed to be Bernoulli distributed, i.e., Yijk|pijk∼ Bernoulli pijk with an unknown mean EYijk=pijk, being related to the predictor variables as follows:(2)fEYijk=logPYijk=1|x,ukj,vkj1−PYijk=1|x,ukj,vkj,=XTβk+ukj+vkj,i=1,…,nj,j=1,…,p.


#### 2.1.1. Univariate Models

Univariate disease mapping approaches based on Poisson or Binomial distributions have traditionally been used in standard analyses. For instance, Lawson [13] developed a hierarchical model in spatial epidemiology using Poisson distribution to estimate the set of disease risks. Ayalew et al. [14] used a Binomial distribution to compare three spatial smoothing models. In the two studies, spatial and nonspatial random effects are introduced into the models. Spatial random effects represent the spatial autocorrelation, which occurs when adjacent regions are more related to each other than more distant regions. The nonspatial component accounts for the heterogeneity, over-dispersion, and arbitrary choice of a spatial unit.For the separate univariate models, the independent variables and random effects are introduced as follows equations (3)–(5).Model M1: Only the non-spatial random effectsvj are considered in the model,(3)fEYij=XTβ+vj.Model M2: Only the spatial random effectsuj are included in the model,(4)fEYij=XTβ+uj.Model M3: Both types of random effectsvj and uj are present (convolution model),(5)fEYij=XTβ+uj+vj,XT in equations (3)–(5) represents a d− dimensional vector of covariates, which is the same for the three separate models, with β the corresponding vector of regression model coefficients.


#### 2.1.2. Bivariate Models

We now extend the notation to two diseases, assuming thatpij1 and pij2 represent the probabilities from the two diseases in region j: j=1,2,…,p. Both of the outcomes follow a Bernouilli distribution.Traditionally, the shared component model and multivariate conditional autoregressive model (MCAR) are two different approaches of implementing a bivariate model to measure the risks of two diseases. A bivariate latent model can also be used to estimate spatial correlation.(1) Shared Component Model. Knorr–Held and Best [15] introduced a new joint modelling paradigm known as the shared component model. For the case of two diseases, the relative risk of each disease depends on an unobserved spatial component shared by both diseases and a disease-specific latent component. The relative risk for each disease is modelled as follows:(6)ModelM4logitpij1=α1+XTβ1+u1jδ+u2j,logitpij2=α2+XTβ2+u1jδ,where u1j is the shared component, common to both diseases while u2j is the component specific to the first disease only. These two components of random effects are modelled using conditional autoregressive priors. The unknown parameter δ>0 is introduced to allow for a differential gradient of the shared component for the two diseases. The ratio of the scaling parameters δ to 1/δ compares the weight of disease 1 relative to disease 2 associated with the shared component.(2) MCAR Model. The multivariate CAR model is a spatial model often used to smooth across space while accounting for associations between the diseases. The fundamental theory for multivariate Gaussian Markov random field (GMRF) has developed by Mardia [16]. Based on his results, Carlin et al. [17] extended multivariate conditional autoregressive (MCAR) models. In the MCAR model, spatial random effects, uj=u1j,u2jT are assigned a multivariate condition autoregressive prior,(7)uj∼MCAR1,Σu,where Σu is the covariance matrix inducing correlation. The resultant MCAR model is expressed by the following equation:(8)ModelM5logitpij1=α1+XTβ1+u1j,logitpij2=α2+XTβ2+u2j.(3) Latent Model. For the full explicit model, the correlation between regions and within a region, we used data augmentation (latent variables) such that,(9)Yijk=0,Wijk≤0,1,Wijk>0.An individual is identified with disease 1 or disease 2 if and only if the underlying latent variableWijk has a positive value. Wijk is a continuous variable that is considered to be normally distributed. The approach adopted offers the opportunity of unravelling the linear relationship between two diseases.The bivariate regression binary augmented model is specified as follows:(10)ModelM6⇒Wij1Wij2∼Nθij1θij2,Σ,Σ=σ112ρwσ1σ2ρwσ1σ2σ222,with(11)θij1θij2∼Nγj1γj2,Λ,Λ=λ112ρbλ1λ2ρbλ1λ2λ222,where Σ and Λ are the within-region and between-region covariance matrices, respectively. The inclusion of within-region correlations (i.e., the ρw and between-region correlations (i.e., the ρb distinguishes bivariate random effects from two independent univariate random effects. We noted that for two independent univariate random effects, ρw=ρb=0 i.e., there is zero correlation.According to Knorr–Held and Best [15], many diseases share common risk factors, which can provide more compelling evidence of true clustering in the underlying risk surface if similar patterns of geographical variation of related diseases are eventually realised. We considered zj as unobserved covariate common to both diseases. The log relative risk γj1,γj2 for disease 1 and 2 are given by the following equation:(12)γj1=α1+ζ1zj,(13)γj2=α2+ζ2zj.αb,b=1,2 in equations (12) and (13) represents individual specific disease intercept and ζ1 and ζ2 are the different risk gradients associated with the covariate for the two diseases.Now, suppose we specify the following model for the log relative risk,(14)γj1=logηjδ,γj2=logηjδ.Here, ηj is the share component.In a general case, the models are expressed as follows:(15)γj1=α1+ζ1zj+φ1qj,(16)γj2=α2+ζ2zj+φ2sj.The contribution of the shared component to overall relative risk is scaled by theδ scaling parameter. qj and sj are diseases’ specific risk factors that are relevant to one or other of diseases only. The cluster model for shared component ηj will capture the spatial distribution of the zj. Whereas the cluster model for specific components ω1 and ω2 will account for the underlying distributions of the qj and sj, respectively. φ1 and φ2 are the different risk gradients associated with the specific risk factors for the two diseases.


### 2.2. Hierarchical Structure

The posterior distributions of the parameters of interest are obtained by updating the priors with the observed data. The first stage is the process model that relates the observed data to the explanatory variables and the spatial component. The second stage is the spatial process model for the spatial random effectsuj and nonspatial random effects v. The parameters β, φ, α and ζ are then assigned prior distributions. β, φ, α, and ζ represent the regression coefficients, different risk gradients associated with the specific risk factors for the two diseases, specific disease intercept, and distinct risk gradients related with the covariate for both the two diseases, respectively. The hyperprior distributions are denoted by the inverse variances of the nonspatial and spatial random effects, 1/σv2 and 1/τu2, respectively. The entire model, as well as its hierarchical structure, is illustrated as follows:(17)level1Yijk|pijk∼BernoullipijkLogitpijk=XTβk+ukj+vkj,i=1,…,nj,j=1,…,p,level2uj∣u−j∼N1nj∑r∼jur,τu2nj,v∼N0,σv2,Priorsβ∼Nμ,Σ,φ∼Nμφ,σφ2,α∼Nμα,σα2,ζ∼Nμζ,σζ2,hyperpriors1τu2∼Gατ,βτ,1σv2∼Gασ,βσ.The Bayesian MCMC simulation involves estimating the posterior distribution of all parameters by integrating prior information with the likelihood for the respective model and sampling each parameter sequentially from its conditional distribution. The three componentszj,qj, and sj are assumed to be independent. The v capture region-wide heterogeneity via an exchangeable normal prior. Then, the uj are the parameters that make this a truly spatial model by capturing regional clustering. They assumed that the random and fixed effects and hyperpriors were mutually independent.In the case of data augmentation for a binary response, the probability is expressed as follows:(18)pi=PWi>0.We can now regard the problem as a data augmentation problem for whichWi’s are not observed but only the indicators are, Yi:=1Wi>0.The likelihood expression is now given by the following equation:(19)∏j=1j∏i=1njpiyij1−pi1−yij=∏j=1J∏i=1nj1Wij>01yij=1+1Wij≤01yij=0ϕWij;μ,σ2.


## 2.1. Specification of the Models

LetYijk be the disease status of individual i for disease k in region j. Especially, k=1 for disease 1 and k=2 for disease 2, j=1,2,…,p and i=1,2,…,nj, where nj is the number of individuals in region j. The response variables are binary and we have the following equation:(1)Yijk=1,if individualiinregionjis positive for diseasek,0,otherwise.Yijk is assumed to be Bernoulli distributed, i.e., Yijk|pijk∼ Bernoulli pijk with an unknown mean EYijk=pijk, being related to the predictor variables as follows:(2)fEYijk=logPYijk=1|x,ukj,vkj1−PYijk=1|x,ukj,vkj,=XTβk+ukj+vkj,i=1,…,nj,j=1,…,p.


### 2.1.1. Univariate Models

Univariate disease mapping approaches based on Poisson or Binomial distributions have traditionally been used in standard analyses. For instance, Lawson [13] developed a hierarchical model in spatial epidemiology using Poisson distribution to estimate the set of disease risks. Ayalew et al. [14] used a Binomial distribution to compare three spatial smoothing models. In the two studies, spatial and nonspatial random effects are introduced into the models. Spatial random effects represent the spatial autocorrelation, which occurs when adjacent regions are more related to each other than more distant regions. The nonspatial component accounts for the heterogeneity, over-dispersion, and arbitrary choice of a spatial unit.For the separate univariate models, the independent variables and random effects are introduced as follows equations (3)–(5).Model M1: Only the non-spatial random effectsvj are considered in the model,(3)fEYij=XTβ+vj.Model M2: Only the spatial random effectsuj are included in the model,(4)fEYij=XTβ+uj.Model M3: Both types of random effectsvj and uj are present (convolution model),(5)fEYij=XTβ+uj+vj,XT in equations (3)–(5) represents a d− dimensional vector of covariates, which is the same for the three separate models, with β the corresponding vector of regression model coefficients.


### 2.1.2. Bivariate Models

We now extend the notation to two diseases, assuming thatpij1 and pij2 represent the probabilities from the two diseases in region j: j=1,2,…,p. Both of the outcomes follow a Bernouilli distribution.Traditionally, the shared component model and multivariate conditional autoregressive model (MCAR) are two different approaches of implementing a bivariate model to measure the risks of two diseases. A bivariate latent model can also be used to estimate spatial correlation.(1) Shared Component Model. Knorr–Held and Best [15] introduced a new joint modelling paradigm known as the shared component model. For the case of two diseases, the relative risk of each disease depends on an unobserved spatial component shared by both diseases and a disease-specific latent component. The relative risk for each disease is modelled as follows:(6)ModelM4logitpij1=α1+XTβ1+u1jδ+u2j,logitpij2=α2+XTβ2+u1jδ,where u1j is the shared component, common to both diseases while u2j is the component specific to the first disease only. These two components of random effects are modelled using conditional autoregressive priors. The unknown parameter δ>0 is introduced to allow for a differential gradient of the shared component for the two diseases. The ratio of the scaling parameters δ to 1/δ compares the weight of disease 1 relative to disease 2 associated with the shared component.(2) MCAR Model. The multivariate CAR model is a spatial model often used to smooth across space while accounting for associations between the diseases. The fundamental theory for multivariate Gaussian Markov random field (GMRF) has developed by Mardia [16]. Based on his results, Carlin et al. [17] extended multivariate conditional autoregressive (MCAR) models. In the MCAR model, spatial random effects, uj=u1j,u2jT are assigned a multivariate condition autoregressive prior,(7)uj∼MCAR1,Σu,where Σu is the covariance matrix inducing correlation. The resultant MCAR model is expressed by the following equation:(8)ModelM5logitpij1=α1+XTβ1+u1j,logitpij2=α2+XTβ2+u2j.(3) Latent Model. For the full explicit model, the correlation between regions and within a region, we used data augmentation (latent variables) such that,(9)Yijk=0,Wijk≤0,1,Wijk>0.An individual is identified with disease 1 or disease 2 if and only if the underlying latent variableWijk has a positive value. Wijk is a continuous variable that is considered to be normally distributed. The approach adopted offers the opportunity of unravelling the linear relationship between two diseases.The bivariate regression binary augmented model is specified as follows:(10)ModelM6⇒Wij1Wij2∼Nθij1θij2,Σ,Σ=σ112ρwσ1σ2ρwσ1σ2σ222,with(11)θij1θij2∼Nγj1γj2,Λ,Λ=λ112ρbλ1λ2ρbλ1λ2λ222,where Σ and Λ are the within-region and between-region covariance matrices, respectively. The inclusion of within-region correlations (i.e., the ρw and between-region correlations (i.e., the ρb distinguishes bivariate random effects from two independent univariate random effects. We noted that for two independent univariate random effects, ρw=ρb=0 i.e., there is zero correlation.According to Knorr–Held and Best [15], many diseases share common risk factors, which can provide more compelling evidence of true clustering in the underlying risk surface if similar patterns of geographical variation of related diseases are eventually realised. We considered zj as unobserved covariate common to both diseases. The log relative risk γj1,γj2 for disease 1 and 2 are given by the following equation:(12)γj1=α1+ζ1zj,(13)γj2=α2+ζ2zj.αb,b=1,2 in equations (12) and (13) represents individual specific disease intercept and ζ1 and ζ2 are the different risk gradients associated with the covariate for the two diseases.Now, suppose we specify the following model for the log relative risk,(14)γj1=logηjδ,γj2=logηjδ.Here, ηj is the share component.In a general case, the models are expressed as follows:(15)γj1=α1+ζ1zj+φ1qj,(16)γj2=α2+ζ2zj+φ2sj.The contribution of the shared component to overall relative risk is scaled by theδ scaling parameter. qj and sj are diseases’ specific risk factors that are relevant to one or other of diseases only. The cluster model for shared component ηj will capture the spatial distribution of the zj. Whereas the cluster model for specific components ω1 and ω2 will account for the underlying distributions of the qj and sj, respectively. φ1 and φ2 are the different risk gradients associated with the specific risk factors for the two diseases.


## 2.1.1. Univariate Models

Univariate disease mapping approaches based on Poisson or Binomial distributions have traditionally been used in standard analyses. For instance, Lawson [13] developed a hierarchical model in spatial epidemiology using Poisson distribution to estimate the set of disease risks. Ayalew et al. [14] used a Binomial distribution to compare three spatial smoothing models. In the two studies, spatial and nonspatial random effects are introduced into the models. Spatial random effects represent the spatial autocorrelation, which occurs when adjacent regions are more related to each other than more distant regions. The nonspatial component accounts for the heterogeneity, over-dispersion, and arbitrary choice of a spatial unit.For the separate univariate models, the independent variables and random effects are introduced as follows equations (3)–(5).Model M1: Only the non-spatial random effectsvj are considered in the model,(3)fEYij=XTβ+vj.Model M2: Only the spatial random effectsuj are included in the model,(4)fEYij=XTβ+uj.Model M3: Both types of random effectsvj and uj are present (convolution model),(5)fEYij=XTβ+uj+vj,XT in equations (3)–(5) represents a d− dimensional vector of covariates, which is the same for the three separate models, with β the corresponding vector of regression model coefficients.


## 2.1.2. Bivariate Models

We now extend the notation to two diseases, assuming thatpij1 and pij2 represent the probabilities from the two diseases in region j: j=1,2,…,p. Both of the outcomes follow a Bernouilli distribution.Traditionally, the shared component model and multivariate conditional autoregressive model (MCAR) are two different approaches of implementing a bivariate model to measure the risks of two diseases. A bivariate latent model can also be used to estimate spatial correlation.(1) Shared Component Model. Knorr–Held and Best [15] introduced a new joint modelling paradigm known as the shared component model. For the case of two diseases, the relative risk of each disease depends on an unobserved spatial component shared by both diseases and a disease-specific latent component. The relative risk for each disease is modelled as follows:(6)ModelM4logitpij1=α1+XTβ1+u1jδ+u2j,logitpij2=α2+XTβ2+u1jδ,where u1j is the shared component, common to both diseases while u2j is the component specific to the first disease only. These two components of random effects are modelled using conditional autoregressive priors. The unknown parameter δ>0 is introduced to allow for a differential gradient of the shared component for the two diseases. The ratio of the scaling parameters δ to 1/δ compares the weight of disease 1 relative to disease 2 associated with the shared component.(2) MCAR Model. The multivariate CAR model is a spatial model often used to smooth across space while accounting for associations between the diseases. The fundamental theory for multivariate Gaussian Markov random field (GMRF) has developed by Mardia [16]. Based on his results, Carlin et al. [17] extended multivariate conditional autoregressive (MCAR) models. In the MCAR model, spatial random effects, uj=u1j,u2jT are assigned a multivariate condition autoregressive prior,(7)uj∼MCAR1,Σu,where Σu is the covariance matrix inducing correlation. The resultant MCAR model is expressed by the following equation:(8)ModelM5logitpij1=α1+XTβ1+u1j,logitpij2=α2+XTβ2+u2j.(3) Latent Model. For the full explicit model, the correlation between regions and within a region, we used data augmentation (latent variables) such that,(9)Yijk=0,Wijk≤0,1,Wijk>0.An individual is identified with disease 1 or disease 2 if and only if the underlying latent variableWijk has a positive value. Wijk is a continuous variable that is considered to be normally distributed. The approach adopted offers the opportunity of unravelling the linear relationship between two diseases.The bivariate regression binary augmented model is specified as follows:(10)ModelM6⇒Wij1Wij2∼Nθij1θij2,Σ,Σ=σ112ρwσ1σ2ρwσ1σ2σ222,with(11)θij1θij2∼Nγj1γj2,Λ,Λ=λ112ρbλ1λ2ρbλ1λ2λ222,where Σ and Λ are the within-region and between-region covariance matrices, respectively. The inclusion of within-region correlations (i.e., the ρw and between-region correlations (i.e., the ρb distinguishes bivariate random effects from two independent univariate random effects. We noted that for two independent univariate random effects, ρw=ρb=0 i.e., there is zero correlation.According to Knorr–Held and Best [15], many diseases share common risk factors, which can provide more compelling evidence of true clustering in the underlying risk surface if similar patterns of geographical variation of related diseases are eventually realised. We considered zj as unobserved covariate common to both diseases. The log relative risk γj1,γj2 for disease 1 and 2 are given by the following equation:(12)γj1=α1+ζ1zj,(13)γj2=α2+ζ2zj.αb,b=1,2 in equations (12) and (13) represents individual specific disease intercept and ζ1 and ζ2 are the different risk gradients associated with the covariate for the two diseases.Now, suppose we specify the following model for the log relative risk,(14)γj1=logηjδ,γj2=logηjδ.Here, ηj is the share component.In a general case, the models are expressed as follows:(15)γj1=α1+ζ1zj+φ1qj,(16)γj2=α2+ζ2zj+φ2sj.The contribution of the shared component to overall relative risk is scaled by theδ scaling parameter. qj and sj are diseases’ specific risk factors that are relevant to one or other of diseases only. The cluster model for shared component ηj will capture the spatial distribution of the zj. Whereas the cluster model for specific components ω1 and ω2 will account for the underlying distributions of the qj and sj, respectively. φ1 and φ2 are the different risk gradients associated with the specific risk factors for the two diseases.


## 2.2. Hierarchical Structure

The posterior distributions of the parameters of interest are obtained by updating the priors with the observed data. The first stage is the process model that relates the observed data to the explanatory variables and the spatial component. The second stage is the spatial process model for the spatial random effectsuj and nonspatial random effects v. The parameters β, φ, α and ζ are then assigned prior distributions. β, φ, α, and ζ represent the regression coefficients, different risk gradients associated with the specific risk factors for the two diseases, specific disease intercept, and distinct risk gradients related with the covariate for both the two diseases, respectively. The hyperprior distributions are denoted by the inverse variances of the nonspatial and spatial random effects, 1/σv2 and 1/τu2, respectively. The entire model, as well as its hierarchical structure, is illustrated as follows:(17)level1Yijk|pijk∼BernoullipijkLogitpijk=XTβk+ukj+vkj,i=1,…,nj,j=1,…,p,level2uj∣u−j∼N1nj∑r∼jur,τu2nj,v∼N0,σv2,Priorsβ∼Nμ,Σ,φ∼Nμφ,σφ2,α∼Nμα,σα2,ζ∼Nμζ,σζ2,hyperpriors1τu2∼Gατ,βτ,1σv2∼Gασ,βσ.The Bayesian MCMC simulation involves estimating the posterior distribution of all parameters by integrating prior information with the likelihood for the respective model and sampling each parameter sequentially from its conditional distribution. The three componentszj,qj, and sj are assumed to be independent. The v capture region-wide heterogeneity via an exchangeable normal prior. Then, the uj are the parameters that make this a truly spatial model by capturing regional clustering. They assumed that the random and fixed effects and hyperpriors were mutually independent.In the case of data augmentation for a binary response, the probability is expressed as follows:(18)pi=PWi>0.We can now regard the problem as a data augmentation problem for whichWi’s are not observed but only the indicators are, Yi:=1Wi>0.The likelihood expression is now given by the following equation:(19)∏j=1j∏i=1njpiyij1−pi1−yij=∏j=1J∏i=1nj1Wij>01yij=1+1Wij≤01yij=0ϕWij;μ,σ2.


## 3. Simulation Design

We generated 1,600 observations (in 8 areas) from a Bernoulli distribution using 2018 Guinea DHS results to approximate the probability of having anemia and malaria [18]. The six models M1,M2,M3,M4,M5, and M6 described in Section 2.1 were simulated to provide a fair comparison between them using R-INLA (Integrated Nested Laplace Approximation) developed by Rue et al. [19] and R2WinBUGS. Two different scenarios were adopted to assess the models M1,M2,M3, and M4. We used INLA with the default options prior in the first scenario, and for the second scenario, we set prior distributions for the parameters (to analyze the sensitivity). The precisions’ default prior is 0.00005. In terms of sitting, according to Neyens et al. [20], precisions of 0.5 and 5 were used in the models. Instead of specifying a prior distribution with INLA for each hyperparameter, a joint prior distribution is considered for the precision matrix in the MCAR model M5. Thus, the precision matrix follows a Wishart distribution. We assigned a uniform prior to the hyperparameter α, which defines the spatial autocorrelation parameter, α∼Un0,1. However, the latent model M6 was assessed using R2WinBUGS, with 10,000 iterations, a burn-in period of 5,000 and 4 chains.To compare the models,M1, M2, M3, M4, and M5 were assessed using the marginal likelihood (MLIK), conditional predictive ordinate (CPO), Deviance information criterion (DIC), and Root Mean Squared Error (RMSE). These criteria can be computed with INLA by setting the appropriate control option when calling inla() [21]. As formal tests for the latent model M6, the DIC, RMSE, and Gelman Robin statistics R^ were used [22]. Following that, the trace plot and density plot are used for visual inspection.


## 4. Results and Discussions

The disease mapping models are simulated to 1,600 observations inn = 8 areas. To assess the performance of the models, various criteria (DIC, CPO, MLIK, and RMSE) were explored. The traditional models are compared to the proposed latent model. We used R-INLA and R2WinBUGS for the simulation. For the latent model, 10,000 iterations with a burn-in period of 5,000 were run, and the MCMC convergence was monitored using the trace plot.According to the results generated with default priors (Table1), the DIC and CPO of the model with nonspatial random effects M1 are smaller than the models M2, M3, and M4, 1610.20 and 1610.28, respectively. These two criteria show that the model fits better with the simulated data. When both random effects are present (convolution model), the model generates the highest MLIK value (−812.95), indicating that the model performs reasonably well. Whereas if we only consider the spatial random effects in the model (M2), it generates a small RMSE (0.2010), which is a good indicator of the precision.Table 1
Summary of values of DIC, CPO, MLIK and RMSE for different models fit (default priors).
ModelsDICCPOMLIKRMSEModel with non-spatial random effectsM11610.201610.28−816.410.2011Model with spatial random effectsM21610.631610.70−819.730.2010Convolution modelM31610.781610.87−812.950.2012Shared component modelM42217.662217.76−1109.500.5000In contrast, the model with spatial random effectsM2 generates small values of DIC and CPO for the precision values of 0.5 and 5. The results in Table 2 show that the DIC and CPO values obtained are 1613.40 and 1613.47, respectively. Table 3 shows 1615.19 for the DIC and 1615.32 for the CPO. In this case, M2 is considered as the best model. As observed in Table 1, when both random effects are present in the model, the highest CPO values generated in Tables 2 and 3 are −819.72 and −827.09, respectively. Similarly, the smallest RMSE values observed in Tables 2 and 3 are 0.2022 and 0.2034, respectively.Table 2
Summary of values of DIC, CPO, MLIK and RMSE for different models fit (precision value of 0.5.).
ModelsDICCPOMLIKRMSEModel with non-spatial random effectsM11614.011614.12−821.130.2029Model with spatial random effectsM21613.401613.47−823.050.2022Convolution modelM31614.951615.07−819.720.2033Shared component modelM42216.392216.41−1111.890.5013Table 3
Summary of values of DIC, CPO, MLIK and RMSE for the different models fit (precision value of 5.).
ModelsDICCPOMLIKRMSEModel with non-spatial random effectsM11615.691615.84−827.890.2037Model with spatial random effectsM21615.191615.32−828.380.2034Convolution modelM31615.901616.05−827.090.2038Shared component modelM42217.442217.49−1116.890.5017As can be seen in Table4, the DIC (−12056.50) and RMSE (0.0002) generated from the latent model M6 are smaller than the ones generated by all other models. Furthermore, the trace plot of inference and convergence show that the R^ value is close to 1 for the correlation parameter, which indicates a good mixing of the four chains and no convergence issues are detected (Figures 2(a), and 3). Also, the density plot shows the normality (Figure 2(b)).Table 4
Summary of values of DIC, CPO, MLIK and RMSE for the MCAR model, and DIC, RMSE for latent model.
ModelsDICCPOMLIKRMSEMultivariate conditional autoregressive (MCAR) modelM53771.803771.80−1731.740.4418Latent modelM6−12056.50——0.0002Figure 2
(a) Trace plot of overall correlation between two outcomes and (b) Density plot of overall correlation between two outcomes from the latent modelM6.
(a)(b)Figure 3
(a) 80% interval for each chain and (b) Gelman Rubin statisticsR^ (overall correlation between two outcomes) from the latent model M6.
(a)(b)Although the latent modelM6 does not generate CPO and MLIK values, the large difference values of DIC and RMSE seen in the table show that this model performs fairly well compared to the models M1, M2, M3, M4, and M5.This study is not exhaustive as it does not cover all aspects of spatial disease modelling. Further research could be conducted to incorporate time effects into the models. Spatial-temporal models, when applied to real-world data, can aid in the discovery of new predictors or explanatory variables. Cressie et al. [23] proposed temporal extensions to hierarchical spatial models based on a parametric description of trends over time, independent risk estimates for each time period, or defining the joint covariance matrix including all periods as a Kronecker product of matrices. Martinez et al. [24] also developed an autoregressive approach to spatio-temporal disease mapping by fusing ideas from autoregressive time series and spatial modelling to link information in time and space. The spatial-temporal models that were developed offered information with a high epidemiological value. Another potential limitation of the proposed latent model is the correlation between variable Y1 at site i and variable Y2 at site j. This was not considered and could be interpreted as spatial cross-correlation of latent outcomes.


## 5. Conclusion

In this study, a bivariate latent model was developed to investigate the relationship between two outcomes within a geographic region. The performance of the proposed latent model was compared to existing models (univariate models, MCAR model, and Shared component model) using the DIC, CPO, MLIK, and RMSE. According to the results, the DIC and RMSE of the latent model developed were smaller than the modelsM1, M2, M3, M4, and M5. Furthermore, the inference and convergence trace plot shows that the R^ value is close to 1 for the correlation parameter, indicating that no convergence issues were detected. From this, it can be concluded that the latent model is considered as the better fit model. In the future, the correlation between variable Y1 at site i and variable Y2 at site j will be investigated.



---

*Source: 7852042-2022-08-11.xml*