Bivariate spatial clustering in differential time trends of related tropical diseases: Application to diarrhea and intestinal parasite infections

There has been a rapid development of space–time multivariate disease mapping methods that focus on space–time variation of risks. Examining the posterior estimates of the space–time random effects can provide compelling epidemiological information that is necessary for public health monitoring. In this study, we propose and evaluate the posterior estimates of the random effects to examine spatial and temporal trends of two tropical diseases. Our model is a multivariate Bayesian space– time model with common spatial and temporal trends, and a space–time interaction term that allows different time trends for different areas. When applied to diarrhea and intestinal parasites data from Ghana, the model that consolidates all random effects as a multivariate conditional autoregressive prior was the best fit. The implementation is based on the notion that diarrheal and intestinal parasite infections share common risk factors Our novel contribution concerns the posterior joint evaluations of the spatial and temporal random effects into a 5 × 5 cross-tabulation table. The method presented is useful for developing and implementing joint epidemiological control strategies, especially in countries where resources are scarce. © 2023TheAuthors.PublishedbyElsevierB.V.Thisisanopen accessarticleundertheCCBYlicense


Introduction
Disease mapping plays a significant role in the geographical surveillance of diseases.There has been significant development and benefit of spatial statistical methods in disease mapping.Current literature on disease mapping has mostly focused on small-area processes.In western countries where effective health surveillance and information systems exist, the protection of patient privacy leads to disease data mostly being available in aggregated formats.In developing countries, data are mostly aggregated because (1) surveillance and health information systems are weak to capture individual-level data, and (2) ethical and legal restrictions limit the acquisition of individual-level data.In both ways, there is the advantage of obtaining estimates of standard incidence ratios (SIR) for small areas.Small areas are often the implementation units for targeting interventions and policy developments.Area-specific estimates of the SIR can therefore be resourceful for allocating public health resources and may further provide clues to the etiology, or generate hypotheses on the etiology (Wakefield, 2007).
The aims of disease mapping are twofold.First, spatial random effects are incorporated to globally or locally smooth the SIRs which are further mapped.Second, the exploration of associations between disease incidences and environmental or demographic exposures can be achieved.The recent focus on dynamic disease mapping has necessitated the development of space-time disease mapping methods that account for spatial, temporal, and spatio-temporal interactions.Spatiotemporal disease mapping by incorporating the temporal dynamics of each observation location can extend the objectives of disease mapping to evaluate the evolution of the disease clusters.The application of Bayesian modeling for disease mapping provides the opportunities to use exceedance posterior probabilities to detect and map higher than expected (hot-spots) and lower than expected (cold-spots) space-time clusters of the SIR (Richardson et al., 2004).
Bernardinelli and colleagues (Bernardinelli et al., 1995) introduced separable parametric spacetime modeling to accommodate the linear temporal changes in disease risk.Using Poisson sampling, they decomposed the log of the unknown space-time risk into a spatial trend common to all periods, a linear temporal trend common to all areas, and a space-time interaction term that allows different time trends for the different spatial units.Knorr-Held (2000) extended the separable space-time models to include nonparametric space-time interaction effects that may include unstructured temporal and unstructured spatial effects (Type I), structured temporal and unstructured spatial effects (Type II), unstructured temporal and structured spatial effects (Type III), and structured temporal and structured spatial effects (Type IV).So far, the focus has been mainly on providing space-time variation estimates of single diseases.The detection of small areas whose time trends depart from the common trend would have significant epidemiological relevance.For instance, epidemiologists may be interested in monitoring small areas with unexpected changes in the disease trend.Building on the work of Abellan et al. (2008) and Li et al. (2012) developed a novel detection method for short-time series of small area data using Bayesian model choice between two competing space-time models.Despite their significant contribution to the spatial epidemiology literature, their method considers unusual time trends in all areas.For resource-scarce countries, the detection of unusual time trends within hot-spots or higher than expected risk areas can reduce the number of areas that need urgent intervention.In our previous study, a adopted the Bayesian space-time model (Bernardinelli et al., 1995) to analyze the spatial clustering of the differential time trends of the incidences.In this study, we aim to develop a methodology to study the temporal dynamics of the joint clustering of related tropical diseases.There are two important questions whose answers would be of strategic importance to epidemiologists and public health managers.First, how do the identified clusters, i.e. hot-spots, cold-spots, or neutral-spots, evolve?Second, do areas that are not yet hot-spots, but e.g.neutral-spots, show the tendency to become hot-spots?These questions are adapted from Li et al. (2014) posing similar questions for univariate space-time variability of burglary risk.The univariate classification proceeds as follows.First, areas are classified as hot-spots, cold-spots, or neutral-spots.Next, the temporal behaviors of these areas are classified as showing a decreasing trend, an increasing trend, or neither.The simultaneous evaluation of multiple diseases that share common environmental factors to answer the questions raised above can serve as a valuable tool for disease surveillance.In this study, we will in addition answer the following three questions in respect of two diseases: -Are the area-specific risks jointly diverging from the same level or converging to the same level as they evolve?-How do their joint hot-spots and cold-spots evolve?-Do the neutral-spots show the tendency to become hot-spots?
To do so, we extend the Bayesian space-time model (Bernardinelli et al., 1995) into a joint model for two diseases.Similarly, we extend the univariate two-stage classification described above into a bivariate classification.We apply our models to study the space-time variability of diarrhea and intestinal parasite infections in Ghana.These two tropical diseases remain a major public health burden in developing countries, especially in sub-Saharan Africa.The methodological contribution is the extension of the univariate two-stage classification into a multivariate two-stage classification.The application of the model and posterior classifications will also have practical and epidemiological value to public health managers.

Statistical modeling
We consider the bivariate spatio-temporal counts y ijk for i = 1, . . ., N s districts, k = 1, 2 diseases, and j = 1, . . ., N T time steps.An extension to k > 2 will be discussed in the application section.In the Bayesian hierarchical modeling framework, we incorporate spatial processes as an intermediate model (process model) between the prior distributions (prior model) and the data likelihood (data model).For the data model, we consider the bivariate spatio-temporal counts y ijk as conditionally independent samples from the Poisson distribution, y ijk |r ijk ∼ P(E ijk r ijk ), where r ijk is the relative risk and E ijk is the expected number of cases.The expected number of cases E ijk = n ij rjk after accounting for district-level population differentials n ij , where rjk = ∑ i y ijk / ∑ i n ij is the overall risk within the study area for time step j and disease k.
We propose a variant of Bernardinelli's (Bernardinelli et al., 1995) method for the process model.Random effect terms are used to uncover any joint spatial, temporal, and spatio-temporal trends.First, we use a monotonically and differentiable log-link function to match the relative risk r ijk and the systematic component, log r ijk ∼ N(η ijk , σ 2 jk ), where the variance parameter σ 2 jk accounts for time and disease-specific over-dispersion in the counts.The variance is time and disease-specific to ensure that log r ijk is exchangeable only within the spatial domain.Next, the linear parameter is decomposed into disease-specific spatial patterns common to all periods, disease-specific linear temporal trends common to all areas, and disease-specific spatio-temporal interaction terms that allow different time trends for the different spatial units x iq γ q + u 0ik .
We refer to this modeling structure as a log-Gaussian parameterization and has minimal Markov chain Monte Carlo (MCMC) convergence issues compared with its log-linear counterpart (Osei et al., 2022).Here, the parameter β 0k denotes the disease-specific overall risk on the log scale, and β1k is a fixed effect parameter for the overall rate of growth for disease k.The time steps are centered at the mid-observation period, t The parameter u 1ik accounts for district-specific spatially structured differential time trends for disease k.Inferentially, β 1ik = β1k + u 1ik specifies the district and disease-specific temporal trends.The coefficients γ q , q = 1, . . ., Q are the fixed effects of Q vulnerability or exposure variables x iq .The disease-specific spatially structured intercepts u 0ik are residual spatial variations that account for unobserved ecological factors which may give rise to spatial clustering.From our parameterization, the unstructured space-time variation is deduced as v ijk = log r ijk − η ijk .Therefore, the disease-specific temporal relative risk is rjk = exp ) .To accommodate for the association between the two diseases via within-area correlations, the assumption of a bivariate conditional autoregressive (BCAR) prior for the structured intercepts u 0i1 and u 0i2 is expedient.The univariate conditional autoregressive model (Besag et al., 1991) extends naturally to the BCAR specification by replacing the univariate normal conditional distribution with a bivariate normal (BN) conditional distribution.Let u 0i = (u 0i1 , u 0i2 ) T be a vector of the random intercepts for k = 1, 2 and u 0,−i = ( u 0,−i1 , u 0,−i,2 ) be the vector of all random intercepts except The weights w il are fixed constants that measure the adjacency of districts i and l.Let the set of boundary points on a district say i be denoted by (i).We defined w il as a binary connectivity weight matrix such that w il = 1 if (i) ∩ (l) ̸ = ∅, and w il = 0 otherwise.The covariance matrix Γ 0 is of dimension 2 × 2. The diagonal elements of Γ 0 are equal to the conditional variances of the random effects.The off-diagonal elements Γ 0 [1, 2] = Γ 0 [2, 1] capture the within-area correlation between them.We add the constraints ∑ i u 0i1 = ∑ i u 0i2 = 0 to ensure identifiability of the mean and proper posteriors.Likewise, for the structured random slopes u 1i1 and u 1i2 , the conditional density ) . Here too, u 1i = (u 1i1 , u 1i2 ) T is the vector of random slopes, ) T is the mean vector of the random slopes, Γ 1 is a 2 × 2 matrix covariance matrix, and is the vector of all random slopes except u 1i .The off-diagonal elements of the covariance matrix Γ 1 [1, 2] = Γ 1 [2, 1] capture the within-area correlations between the random slopes in this case, whereas the diagonal elements capture their conditional variances.
The incorporation of correlation between the random slopes and intercepts could answer critical etiological questions such as how the population has reacted towards environmental factors introduced at the reference time of the study.When the possible correlation between the random intercepts and slopes of the time trends is avoided, the area-specific time trends tend to be pulled towards the overall mean trend (Bernardinelli et al., 1995).In a previous study (Osei and Stein, 2019) we explored different mechanisms of inducing correlations and found the BCAR model to be superior to the other models for a single disease.Considering two diseases in this study, we have four spatial random effects, u = (u 0i , u 1i ) T .The BCAR model naturally extends to the MCAR model as et al., 2005).Now, Γ is a 4 × 4 covariance matrix with diagonal elements equal to the conditional variances of the random effects.The off-diagonal elements capture the within-area correlation between the 4 random effects.

Application to diarrhea and intestinal parasites infections in Ghana
More than 1.7 billion episodes of diarrhea are recorded globally every year with the majority of these occurring in low and middle-income countries (Black et al., 2010(Black et al., , 2003;;Boschi-Pinto, 2008;Fischer Walker et al., 2012).Intestinal parasites infection on the other hand is estimated to infect more than 1 billion people, although there are concerns about the actual number of infections (Bethony et al., 2006;de Silva et al., 2003).Available literature indicates that diarrhea is an important symptom of intestinal parasites infection.Hence, symptomatic diagnostics which is the common practice of diagnosis in developing countries could attribute some proportion of intestinal parasites cases to diarrhea.Understanding their joint burden through joint modeling plays a fundamental public health role to achieve the sustainable development goals (SDGs).
To demonstrate the methods described above, we used yearly diarrhea and intestinal parasites morbidities from 2010 to 2014 from outpatient records.They were obtained from the District Health Information Management System (DHIMS) managed by the Centre for Health Information and Management (CHIM) of the Ghana Health Services (GHS).The data exist in aggregated format per administrative district.The geographical scale of analysis is restricted to the 170 administrative districts of which data had been recorded for the period 2010 to 2014.The population data for the 170 districts were obtained from the Ghana Statistical Service (GSS).The summary statistics and the spatial distribution of the relative risks are shown in Table 1 and Fig. 1, respectively.Within the study period, the relative risk for diarrhea ranged from 0 in 2010 to 12.50 in 2012.That of intestinal parasites infection also ranged from 0 in 2010 to 19.97 in 2012.The spatial distribution between the two diseases appears similar within the southern belt, but diverse within the northern belt (Fig. 1).
We fitted two models.Model 1 imposes a correlation between the random intercepts and random slopes separately for each of the two diseases via the BCAR formulation.Model 2 is an extension of model 1 where all the random effects are simulated from the MCAR density.Both models do not include covariates.Model 2, therefore, has six correlation coefficients that  can be estimated between the random effects, i.e. within-area correlation between the random intercepts of intestinal parasites and diarrhea; within-area correlation between the random slopes of intestinal parasites and diarrhea; within-area correlation between the random intercepts and slopes of diarrhea; within-area correlation between the random intercepts and slopes of intestinal parasites; within-area correlation between the random intercepts of diarrhea and random slopes of intestinal parasites; and within-area correlation between the random slopes of diarrhea and random intercepts of intestinal parasites.The last four of these correlations express the correlations between random intercepts and slopes.The structure of the models is then as follows: Model 1: Model 2: The models can be extended to more than two diseases.For k = 1, ..K diseases, the random components of Model 1 become {u 0i1 , u 0i2 , . . ., u 0iK } ∼ MCAR (Γ 0 ) and {u 1i1 , u 1i2 , . . ., u 1iK } ∼ MCAR (Γ 1 ) where Γ 0 and Γ 1 are of dimension K ×K .That of Model 2 becomes {u 0i1 , u 0i2 , . . ., u 0iK , u 1i1 , u 1i2 , . . .,

Model estimation
In our study, we estimated the model parameters by generating samples from the posterior density via Gibbs sampling within a Bayesian framework.Let ψ 1 be the full Gaussian latent (unobservable) field and let ψ 2 be a vector of hyper-parameters.Then

}
. Our model is then summarized within a three-stage hierarchical framework: a data model, a process model, and a parameter model.Stage 1 is equivalent to the data model y|ψ 1 , ψ 2 ∼ p (y|ψ 1 , ψ 2 ); Stage 2 is the process model ψ 1 |ψ 2 ∼ p (ψ 1 |ψ 2 ), while Stage 3 is equivalent to the parameter model ψ 2 ∼ p (ψ 2 ).The joint posterior by the chain rule equals p (ψ 1 , ψ 2 |y) ∝ p (y|ψ 1 , ψ 2 ) × p (ψ 1 , |ψ 2 ) × p (ψ 2 ).Before estimation, we assign prior distributions to the variance parameters σ 2 jk , the disease-specific fixed effects β 0k , β1k and the covariance, for instance, Γ as in Model 2. In our study, we assigned a non-informative flat prior distribution for the intercept, p(β 0k ) ∝ 1, to ensure that the data drive the inference.Non-informative priors result in posterior inference similar to maximum likelihood inference (Waller and Gotway, 2004).For the disease-specific slopes β1k , we assigned vague normal priors β1k ∼ N(0, 10 5 ), being the equivalent of a non-informative Gaussian prior.We assigned a Wishart prior to the inverse of the covariance matrix, thus Γ −1 ∼ W(Ω, df ).This prior is a conjugate of the inverse of the covariance parameters in a multivariate normal distribution (Gelman et al., 2013;Press, 2005).The scale matrix is Ω and the number of degrees of freedom df equals the number of random effects for a weakly informative distribution.We set Ω as a scaled identity matrix of dimension equal to the number of random effects, a specification (Moraga and Lawson, 2012) utilized in simulation studies of MCAR modeling.We fitted the models using the R2WinBUGS package (Spiegelhalter et al., 2008b) together with the R software (R Core Team, 2016).

Model evaluation and comparison
We evaluated the adequacy of the models using cross-validating predictive checks.We used the Chi-square goodness-of-fit statistic based on the discrepancy function, and Spiegelhalter, 2003).Similarly, we obtained χ 2 pred for the predicted datasets y pred ijk .The two quantities were compared using Bayesian p-values, as the probability Pr . Bayesian p-values close to 0.5 commonly suggest that the generated data are compatible with the model.Since their distribution is not symmetrical, we considered values <0.1 and >0.9 as extreme values and, hence, as an indication of poor fit (Tzala and Best, 2008).For model comparison, we  (Spiegelhalter et al., 2002).Negative twice the log-likelihood of the deviance informs the model fit, while the effective number of parameters informs the model complexity.Like the AIC, a low DIC value correspond with a good predictive performance of the model.

Spatial clustering of time trends
Our prime interest concerns the spatial clustering of the residual relative risk and the differential time trends.To evaluate this, we adapt the univariate spatial clustering classification rule by Richardson et al. (2004) and Li et al. (2014) and extend it towards a bivariate situation as a twostep procedure.First, we classified districts as either hot-spot, cold-spot, or neutral-spot based upon the posterior probabilities of the residual risks.Following Richardson et al. (2004), we define a district as a univariate hot-spot if for each disease k = 1, 2 the posterior probability that p (exp (u 0ik ) > 1|data) > 0.8.A district is defined as a cold-spot if the posterior probability that p (exp (u 0ik ) > 1|data) < 0.2.Neutral-spots are those districts with posterior probabilities 0.2 ≤ p (exp (u 0ik ) > 1|data) ≤ 0.8.Second, for the classification of time trends into decreasing, stable and decreasing trends, we used the same probability cutoffs applied on the local slopes, β 1ik .Thus, districts with increasing time trends are defined as those trends with p (β 1ik > 0|data) > 0.8, decreasing time trends as those with p (β 1ik > 0|data) < 0.2, and stable time trends as those with 0.2 ≤ p (β 1ik > 0|data) ≤ 0.8.
For the joint or bivariate clustering of the two diseases, we defined five categories of spatial clusters and time trends.These definitions are similar to the univariate case, except that they are based upon two posterior probability estimates.Table 2 summarizes the classifications of districts into 5 × 5 clusters and time trends.For instance, we defined a district as a bivariate hot-spot (Hothot-spot) if the posterior probabilities p (exp (u 0i1 ) > 1|data) > 0.8 and p (exp (u 0i2 ) > 1) > 0.8.If p (exp (u 0i1 ) > 1|data) > 0.8 and p (exp (u 0i2 ) > 1|data) < 0.2 then the district is defined as a hot-spot for disease 1 and a cold-spot for disease 2 (Hot-cold-spot).
We report our results based upon the posterior samples from two MCMC chains, each with 20 000 iterations and a burn-in of 10 000.The convergence checks were by visual examination of the autocorrelation trace plots.We also formally assessed convergence using the Brooks and Gelman plots and shrink factor (Brooks and Gelman, 1998).Here, convergence is said to be attained when the shrink factor approaches 1.

Results and analysis
Fig. 2 shows the Brooks and Gelman (1998) plots for the fixed and variance parameters (β 0k , β1k , σ 2 tk ).For each of these parameters, the shrink factor approached 1 as the number of iterations increased, suggesting convergence.The estimated Bayesian p-values of the Chi-square goodness of fit discrepancy statistics also suggest that the models adequately fit the data.
Model comparison was based on the DIC values.DIC differences of less 1-2 are treated as negligible, 3-6 as moderate, and 7 as substantial (Spiegelhalter et al., 2002).Table 3 shows that Model 1 has the highest DIC value as compared to Model 2. The difference in DIC equals 19.4, and it suggests the support for Model 2 over Model 1.The indication is that there are significant spacetime correlated residuals within and between the diseases that need to be considered.Exponents of the intercepts represent the overall relative risks after accounting for fixed and random effects.
Thus, the relative risks of diarrhea and intestinal parasites are exp (β 01 ) ≈ 1 and exp (β 02 ) ≈ 0.77 after controlling for the spatial trend common to all periods, a linear temporal trend common to all areas, and a space-time interaction term that allows different time trends for the different areas.The time trends common to all districts show positive slopes for both diseases.But intestinal parasites infection has a slightly higher slope (exp Table 4 shows the covariances and correlations between the random effects.For both diseases, the variances of the random intercepts exceed that of the random slopes.The variance of the random intercepts of intestinal parasites is also nearly twice that of diarrhea, i.e. σ 2 u 02 = 3.085 for intestinal parasites against σ 2 u 01 = 1.469 for diarrhea.Similarly, for the random slopes, the variance of intestinal parasites equals σ 2 u 12 = 0.254, and that of diarrhea equals σ 2 u 11 = 0.149. Fig. 3 (left panel) shows the spatial variabilities of the exponentiated residual relative risks (i.e. the random intercepts) and the residual time trends (i.e. the random slopes).These can be interpreted as multiplicative effects.There are obvious spatial similarities between the residual relative risks of both diseases.We observe spatial similarities also for their residual time trends.We can affirm this from Table 4; we notice that the two diseases are positively correlated in terms of their residual risk (R 31 = 0.66) and time trends (R 42 = 0.64).Thus, districts of similar risks and time trends of both diseases agglomerate, indicating that the two diseases share similar environmental or demographic determinants.We observe a similar within-area negative correlation between the residual risks and the time trends for both diarrhea R 21 = −0.48,and intestinal parasites R 43 = −0.47.This suggests that the district-specific disease rates are jointly converging to the same level.Fig. 3 (right panel) also shows the estimated common temporal relative risks and their 95% Bayesian credible intervals (shaded regions).The common temporal relative risks show similar general increasing trends despite the occurrence of a large cluster with decreasing time trends within the south-central areas.
Table 5 shows the results of the two-stage univariate classification of clusters into increasing, stable, and decreasing time trends.For only diarrhea, 46.5% of districts were classified as hot-spots, 39.4% as cold-spots, and 14.1% as neutral-spots.We observed a comparable pattern for intestinal parasites.Figs. 4,5,and 6 show the spatial distribution of the time trends within the hot-spots, cold-spots, and neutral-spots, respectively.Associated with each spatial cluster are the plots of district-specific relative risks for those with increasing time trends (for only two districts each) together with the common temporal relative risks.Within each clustering category, the stable time trends predominate except for the cold-spots.For diarrhea, 11 out of the 79 hot-spots are detected as increasing time trends; for intestinal parasites, 6 out of the 84 hot-spots showed increasing time trends.
The bivariate or joint cross-classification of clusters and their time trends resulted in 25 categories, consisting of the 5 cluster categories and the 5 time trend categories (Table 6).The   and 8 show the spatial distribution of the time trends for each of the 5 cluster categories.Amongst the 54 districts which are hot-hot-spots, i.e. hot-spots for both diseases, only one district shows a bivariate increasing time trend, while 34 districts are stable.Of the 34 districts identified as cold-cold-spots, i.e. cold-spots for both diseases, 12 show bivariate increasing time  trends.Amongst the 48 neutral districts, the majority (31 districts) are stable, while six districts have increasing time trends.Since districts with increasing time trends are critical, in Fig. 7 (right panel), we also show the district-specific temporal relative risks for the hot-hot-spot and two cold-coldspots with bivariate increasing time trends.Their deviations from the common temporal relative risks are widely noticeable.In Fig. 8, we also observe that the bivariate neutral-spots have temporal relative risks, as shown for two of the districts, that deviate widely from the common temporal relative risks.

Discussion
In this paper, we considered a multivariate space-time statistical modeling framework and developed a bivariate posterior analysis of time trends within spatial clusters.We have adapted and extended the development proposed by Li et al. (2014) towards the joint counts on two diseases.Our study is specifically applied to analyze diarrhea and intestinal parasites morbidities in Ghana.Space-time parametrizations of the random effects have been utilized to capture variances across space, time, and space-time.The fitted models, especially Model 2, serve as a multivariate version of the general space-time framework proposed by Knorr-Held (2000) for space-time disease mapping.There are two primary points of discussion regarding the methods.The first is the modeling framework.Knorr-Held (2000) proposed 4 types (type I, II, III, and IV) of inseparable spacetime interaction mechanisms for the space-time modeling framework.This modeling framework has gained attention in disease mapping literature.Our modeling framework can be seen as a multivariate version of the type IV interaction, where a linear trend with a local slope is specified for each spatial unit and the local slopes are smoothed over space.Our model further smooths the time trends and the residual risks across the two diseases.This allowed us to quantify the correlation between the two diseases in terms of their residual risks and time trends.Consolidation of all random effects via MCAR modeling led to a substantial reduction in the DIC (Spiegelhalter et al., 2002), implying that the inclusion of within-area correlations amongst random effects enjoys a better fit than when excluded.There are alternatives to the space-time structure used in this study, in particular when one is interested in making spatial and temporal inferences at a resolution higher than the level at which the data have been observed.The space and time domains can be viewed as either discrete or continuous.There is rich literature on the integration of Gaussian random fields (Gaussian processes) for either the spatial domain or temporal domain or both domains, e.g.Gneiting (2002) and Stein (2005).In the future, we intend to extend the modeling structure to be based upon dynamic MCAR models that evolve continuously over time so that investigators can have the opportunity to make local inferences of temporal effects at a resolution (temporal) finer than that at which the data are observed.This is achievable if the temporal process is assumed to occur within a continuous domain and modeled via a Gaussian random field process.
The second primary point of discussion concerns the strength of the approach to detect disease clusters and how it is set apart from current methods.Often, there is a focus on either the spatial, temporal, or space-time variabilities of the disease data.In some occasions, the exceedance posterior probabilities of the residual risks establish the hot-spots, cold-spots, or the neutralspots.The exceedance probability cutoffs dwell on the seminal work of Richardson and colleagues (Richardson et al., 2004) and so is our methodology.Our study extends their approach by exploring the tendency of the different types of spots to evolve within different time trends.We do not question the sensitivity of the probability boundaries for these classifications since Richardson et al. (2004) already concluded on their robustness.As suggested by Li et al. (2014), the performance of this classification criterion and the sensitivity of the probability cutoffs for the hot-spots and cold-spots need to be assessed through simulations.Since using a stricter cutoff (i.e.0.1-0.9)did not show different results, we conclude that our results are reliable.The Bayesian model selection method, BaySTDetect (Li et al., 2012), is an alternative to detect areas of unusual departures from the common time trend.This method may be applicable if extended toward a multivariate framework.Yet, the added value of our method is that we can detect the areas of unusual times trends within the different types of spots.Our method efficiently scales down the number of implementation units that may require urgent public health attention in the presence of scarce resources.More flexible time trends like random walk models or splines may be fitted to the common and area-specific time trends if there is data with many time steps.Data used in this study only included time steps which are arguably not adequate for nonlinear trend modeling.Also, for data with more time steps, our model can be modified for space-time-series modeling of multiple diseases.
We now turn to the epidemiological relevance of our results in terms of public health monitoring and evaluation.We found no existing study with which to compare our results.To the best of our knowledge, our study is the first on the small-area spatial similarity between diarrhea and intestinal parasites infections.The epidemiological significance of the correlations cannot be overlooked.Here, we observed a negative correlation between the random intercepts and slopes, suggesting that the district-specific disease rates, at different levels in the past, are converging to the same level (Bernardinelli et al., 1995).The possible explanation is that demographic or environmental determinants which were highly variable in the past are becoming homogeneous with time.Considering this observation together with the yearly increasing trends for both diseases, we argue that the trends of such possible determinants are deteriorating.
The increasing common time trends for both diseases imply that interventions so far, if any, have not been fully effective.This is worrying, especially when increasing trends are predominant even within cold-spots.If we juxtapose this finding with the district-specific disease rates that are jointly converging to the same level, it becomes highly urgent to call for a review of the existing health intervention strategies such as health promotion.Based on our hypothesis that increasing time trends within cold-spots and neutral-spots have a higher tendency to become hot-spots, we make the following remarks.There are 12 cold-cold-spots that have a higher tendency to migrate into hot-hot-spots.There are 6 bivariate neutral-spots that have the tendency to migrate into hothot-spots.A large proportion of the univariate and bivariate hot-spots show the tendency to migrate into neutral or cold-spots since they have decreasing time trends.While urgent attention is required at the hot-spots with increasing time trends, continuous monitoring should be undertaken at the cold-spots and neutral-spots with increasing time trends.
The study also unveils significant spatial within-area correlation for residual intercepts and residual time trends.This indicates spatial and temporal similarities between diarrhea and intestinal parastates infections, suggesting that the two diseases have similar transmission pathways.Diarrhea and intestinal parasites are influenced by similar environmental and demographic determinants.Since these determinants are spatially continuous, joint clustering of their relative risks and time trends is expected.Diarrhea infection is mainly through water and food contaminated with feces, due to poor hygiene (Vesikari and Torun, 1994).Its persistence has been attributed to socioeconomic inequalities such as low-income levels, illiteracy, and inadequate safe water and sanitation (Dasgupta, 2008;Mekasha and Tesfahun, 2003;Woldemicael, 2001).Similarly, infection with intestinal parasites occurs mainly by contact with infected environments, hand-to-hand contact, or contaminated food or water (fecal-oral) (Alum et al., 2010) which are also indicators of socioeconomic inequalities.The common intestinal parasites thrive under climatic and environmental conditions such as high temperatures, severe precipitation, and adequate soil moisture (Brooker et al., 2006;Brooker and Michael, 2000;Chammartin et al., 2014).In our previous studies, the risk factors observed for diarrhea were similarly observed for intestinal parasites infections (Osei et al., 2019;Osei and Stein, 2017).

Conclusions
We have presented an extension of a univariate cross-classification of clusters and time trends into a multivariate setting of multiple diseases.It expands the current state of univariate modeling of differential time trends into bivariate modeling.We focused on the joint spatial clustering of time trends rather than on the estimation of only space-time variances.The bivariate model applied to diarrhea and intestinal parasite infections in Ghana classified time trends within the different spatial clusters.Our findings generally imply that our approach can reduce the number of spatial units that may require location-specific interventions in the presence of scarce resources.The method was able to evaluate whether the two diseases are converging or diverging together.In our application, we have observed that diarrhea and intestinal parasites infections are converging to the same level.
Epidemiologically, the study has demonstrated the applicability of our bivariate model to evaluate the spatial clustering of tropical diseases.The observed correlation between diarrhea and intestinal parasites across various years suggests that they share similar predisposing determinants.Future work will consider elucidating separate correlations for the different time steps using nonparametric random walk priors.The greater impact of this study concerns both the findings related to the specific disease data as well as the subsequent application of the methods by public health professionals.There may, however, be a capacity gap in terms of the technical staff at public health services to reproduce these results and the application to new data.Our next step is to facilitate the development of the capacity of local health officials, especially those in sub-Saharan Africa, with simplified workflow applications.

Fig. 2 .
Fig. 2. Brooks and Gelman plots for the fixed and variance parameters.The shrink factor for all the parameters approaches 1 as the number of iterations increases.

Fig. 3 .
Fig. 3. Spatial distribution of the posterior means of the common residual relative risks and district-specific residual time trends (Left), and the posterior means of the common temporal relative risks (Right): A: spatially correlated common residual relative risks for diarrhea.B: spatially correlated district-specific residual time trends for diarrhea.C: spatially correlated common residual relative risks for intestinal parasites.D: Spatially correlated district-specific residual time trends for intestinal parasites.

Fig. 4 .
Fig. 4. Spatial distribution of the time trends within the univariate hot-spots (Left).Increasing district-specific time trends (for two districts whose IDs are also shown on the maps) are also plotted together with the common time trends (Right).

Fig. 5 .
Fig. 5. Spatial distribution of the time trends within the univariate cold-spots (Left).Increasing district-specific time trends (for two districts whose IDs are also shown on the maps) are also plotted together with the common time trends (Right).

Fig. 6 .
Fig. 6.Spatial distribution of the time trends within the univariate neutral-spots (Left).Increasing district-specific time trends (for two districts whose IDs are also shown on the maps) are also plotted together with the common time trends (Right).

Fig. 7 .
Fig. 7. Spatial distribution of time trends within bivariate hot-hot-spots and cold-cold-spots.Increasing district-specific time trends (for districts whose ID are also shown on the maps) are also plotted together with the common time trends (Right).

Fig. 8 .
Fig.8.Spatial distribution of time trends within bivariate neutral-spots, hot-cold-spots, and cold-hot-spots.Increasing district-specific time trends (for two districts whose ID are also shown on the maps) for the neutral-spots are also plotted together with the common time trends.

0i |u 0,−i ∼ BN
The conditional density is the bivariate normal density u

Table 1
Summary statistics of the yearly relative risk of diarrhea and intestinal parasites from 2010 to 2014.

Table 2
Criteria for the two-stage bivariate classification of clusters and time trends based on the posterior probabilities of the residual risk p (exp (u 0ki ) > 1) and time trends p (β 1ik > 0).
used the deviance information criterion (DIC), a generalization of Akaike's information criterion (AIC).The DIC = D + p D is the sum of the model fit D and model complexity p D

Table 3
The posterior means of fixed effects and model fit parameters of Model 1 and Model 2. The 95% posterior credible intervals (CI) are shown in brackets.

Table 4
The covariance and correlation matrix for the random effects in Model 1 and Model 2.
a Not applicable.

Table 5
Results from the univariate cross-classification (Clustering and time trends) of the 170 districts.
in the category 'neutral-spots', either or both of the two diseases have neutral-spots.In addition toTable 6, Figs.7

Table 6
Results from the bi-variate cross-classification (Clustering and time trends) of the 170 districts.