A Systematic Review of Joint Spatial and Spatiotemporal Models in Health Research

With the advancement of spatial analysis approaches, methodological research addressing the technical and statistical issues related to joint spatial and spatiotemporal models has increased. Despite the benefits of spatial modelling of several interrelated outcomes simultaneously, there has been no published systematic review on this topic, specifically when such models would be useful. This systematic review therefore aimed at reviewing health research published using joint spatial and spatiotemporal models. A systematic search of published studies that applied joint spatial and spatiotemporal models was performed using six electronic databases without geographic restriction. A search with the developed search terms yielded 4077 studies, from which 43 studies were included for the systematic review, including 15 studies focused on infectious diseases and 11 on cancer. Most of the studies (81.40%) were performed based on the Bayesian framework. Different joint spatial and spatiotemporal models were applied based on the nature of the data, population size, the incidence of outcomes, and assumptions. This review found that when the outcome is rare or the population is small, joint spatial and spatiotemporal models provide better performance by borrowing strength from related health outcomes which have a higher prevalence. A framework for the design, analysis, and reporting of such studies is also needed.


Introduction
Recent advances in geographic information systems in medicine have led to the development of advanced spatial analysis of geocoded data in health research [1]. Disease mapping, also defined as the spatial analysis of disease risk, is an important area of public health research [2]. Commonly georeferenced data used in epidemiological investigations have information about space and perhaps time [3]. In spatial epidemiology, Tobler's first law of geography is considered as the foundation of spatial statistics, which asserts that everything is related to everything else, although proximate things are more closely related than distant things [4]. Nearby areas are more likely to share similar geographic characteristics linked to the disease and, likewise, the temporal dependence is greater for succeeding years than for years apart from one another [5,6].
Spatial models have been in use in the field of public health research for decades [7], and important progress over decades has enabled the development of complex models to examine a potential correlation between disease patterns and covariates that are geographically and/or temporally varied [8]. In disease risk mapping, the Standardised Incidence Rate (SIR) and Standardised Mortality Ratio (SMR) are frequently used to measure spatial Therefore, we conducted a systematic review of joint spatial and spatiotemporal models applied to health outcomes to provide insightful recommendations for future researchers as to when and how to fit a joint spatial and spatiotemporal model. This systematic review helps to improve the decision-making process through the joint spatial and spatiotemporal modelling of two or more health outcomes. In addition, the results could inform researchers in terms of providing insights about advanced joint spatial and spatiotemporal statistical methods and related issues.

Data Source and Search Strategy
We performed a systematic review of peer-reviewed published health research that employed joint spatial and spatiotemporal methods. For the formulation of the systematic review methodology, we used the Preferred Reporting Items for Systematic Reviews and Meta-analysis (PRISMA) checklist [66]. We registered the systematic review on the PROSPERO international prospective register of systematic reviews (registration number: CRD42022365445). A comprehensive search strategy was carried out for joint spatial and spatiotemporal models (joint spatial autocorrelation, joint spatiotemporal autocorrelation, joint spatial model, and joint spatiotemporal models) applied to any health or health-related outcomes with no geographic limits.
In our review, a spatial model incorporates a geo-spatial index, a temporal model includes a time index, a spatiotemporal model involves both a geospatial and time index, and joint spatial and spatiotemporal models accommodate a geospatial and/or time index of two or more health outcomes [67,68]. The search was conducted on 19 September 2022. Databases such as PubMed, Medline, Scopus, PsycINFO, Emcare, and Embase were searched. The reference lists of retrieved studies were further searched on Google Scholar and advanced Google to identify more papers. Search terms for the joint spatial and spatiotemporal studies are detailed in Table S1. Studies published between January 2011 and October 2022 without geographic restrictions were considered in our review.
Retrieved articles from each database were exported to Endnote version 20 reference citation software and stored as a single file name and then exported to Covidence for further processing (Covidence systematic review software (Veritas Health Innovation, Melbourne, Australia. Available at www.covidence.org, accessed on 20 August 2022)). Duplicates were deleted manually during the title and abstract and full-text screening and automatically in Endnote and Covidence software. Searches were conducted using the following terms: "multivariate spatiotemporal" OR "bivariate spatiotemporal" OR "multivariate spatiotemporal" OR "bivariate spatio-temporal" OR "joint shared spatial model*" OR "joint spacetime model" OR "multivariate space-time model*" OR "bivariate space-time model" OR "small area analys*" AND "shared component model" OR "disease mapping" AND "shared component model" OR "space-time mixture model" OR "shared component model" OR "spatial analys*" AND "joint model*" OR "joint spatial model*" OR "joint spatial analys*" OR "shared latent component model" OR "joint model*" AND "spatial model*" OR "spatial factor analys*" OR "risk map*" AND "shared component model*" OR "shared spatial model*" OR "multivariate spatial analys*" OR "bivariate spatial analys*" OR "bivariate conditional autoregressive model" OR "multivariate conditional autoregressive model" OR "joint conditional autoregressive model" OR "joint spatial autocorrelation" OR "bivariate spatial autocorrelation" OR "multivariate spatial autocorrelation" OR "spatial co-cluster*" OR "spatio-temporal co-cluster*" (Table S1).

Inclusion and Exclusion Criteria
This systematic review covered peer-reviewed studies published in the English language between 2011 and 2022 modelled using joint spatial and spatiotemporal models. The year 2011 was chosen as the starting point because the joint spatial and spatiotemporal analysis of two or more health outcomes was widely implemented as a new area of study over the past decades and because of the need for most recent evidence due to the rapid changes in analytical techniques due to the ongoing advancement in science and technology. Studies before that date were either outdated or superseded by newer methods included in our review. Every article retrieved from the databases (Medline, PubMed, PsycINFO, Emcare, Scopus, and Embase) was exported to Endnote. After excluding duplicates from Endnote, we transferred to Covidence for additional article screening and extraction.
Title and abstracts and full-texts were screened by two authors (GAT and ZTT) independently to identify eligible studies based on the inclusion and exclusion criteria. When conflicts emerged over the inclusion or exclusion of studies, a consensus was reached through discussion, and if the conflicts were not resolved a third reviewer (AE) was consulted. Research articles performed their analyses using the joint spatial approaches; joint spatial and spatiotemporal autocorrelations, joint spatial models, and joint spatiotemporal models, and the outcomes analysed (could be on health outcomes among humans (not animal study)) were eligible for this review. No exclusion was made based on geography and types of health outcomes studied. Conference abstracts, reviews, texts published in a language other than English, non-human studies, and those not considering joint spatial and spatiotemporal models were considered as exclusion criteria.
Methods for spatially and temporally modelling two or more health outcomes, or the same health outcome in two or more subsets of the population at risk, are referred to as joint spatial and spatiotemporal methods [69].

Data Extraction
A data extraction template was developed in Microsoft Excel. The tool was developed considering the review question. Two authors (GAT and ZTT) extracted the data independently. When a disagreement appeared, it was resolved by a third author (AE). The data abstraction tools contained key information such as bibliographic information, research study objectives, the nature and type of data, covariate type, and data analysis methods, i.e., modelling approaches, key findings, and methodological gaps.
For each study, data such as the last name of the author, article title, name of the journal, year of publication, country, data source, spatial data type, outcomes of interest, number of outcomes, incidence/prevalence of outcomes, inference approaches, estimation techniques, study design, spatial unit, number of the spatial unit, temporal unit, number of temporal units, objectives, sample size, spatial model type, spatial structure, temporal structure, space-time interaction term, assigned priors, the reason for using joint spatial modelling, covariates used in the model, variable selection approach, number of covariate considered in the model, standardisation, spatial neighbourhood structure, temporal adjacency, software used, analysis method, model validation, model comparison measures, effect measure reported, key findings, map reported, script provided, and methodological gaps were extracted.

Risk of Bias Assessment
Two authors carried out a thorough assessment of the included studies' methodological quality (GAT and ZTT). All included studies' risk of bias was evaluated using a quality assessment tool that has an 8-point scoring system updated and modified to evaluate each study's quality based on its aims and objectives, model validity, overall results, and study conclusion [70,71]. A standardised item list was used to grade the quality and risk of bias of included studies (Table S2). The checklist consists of 8 questions with possible answers ranging from 0 to 2, with a maximum overall score of 16. Low-, medium-, high-, and very high-quality levels were used to classify the overall score (low = score < 8, medium = score 8-10, high = score 11-13, and very high >13). Two authors (GAT and ZTT) independently assessed each study to determine its score and to determine the overall quality of the included studies. Discussion between the two authors attempted to settle any differences, and for those that could not be settled, a third reviewer (AE) was engaged.

Data Synthesis and Analysis
Microsoft Excel and STATA version 17 software were used for data entry and analysis, respectively. The review findings were summarised into texts, tables, and figures. The descriptive analysis was presented using the proportion, mean, medians, and ranges.

Search Results and Characteristics of Included Studies
A comprehensive search of international peer-reviewed journals yielded a total of 4077 published articles. Of these, 4071 articles (PubMed: 2813, Embase: 103, PsycINFO: 24, Emcare: 45, Medline: 82, and Scopus: 1004) were obtained from 6 databases, and 6 more studies were found by manual searches in Google and Google Scholar. All these articles were published from 2011 onwards. During title and abstract screening, about 168 duplicates, 1278 non-relevant studies, and 2544 non-joint spatial studies were discarded, and only 87 studies were left for the full-text screening. A total of 44 studies were excluded from the 87 articles that met the inclusion and exclusion criteria for full-text screening; of these, 32 were not joint spatial models, 4 were methodological reviews, 3 were multivariable analyses, 2 were animal studies, 1 article was a duplication, 1 was a genetic study, and 1 was not a spatial study. The systematic review included 43 research publications that met the eligibility criteria (Table S3 and Figure 1).   Since 2011, the number of publications has fluctuated; it increased in 2016 and reached a peak in 2019 and 2020, then sharply dropped in 2021 before rising again in 2022. More ongoing research in this area is also anticipated with an increase in the number of clinical registries being set up. The majority of studies (n = 30, 69.78%) have been published after 2016 ( Figure 2). Nearly one-fifth (18.6%) and 9.3% of the studies were conducted in Iran and the United States of America (USA), respectively ( Figure 3). Of the 43 total studies, 15 (34.88%) were applied to infectious diseases, such as HIV/AIDS, herpes simplex virus-2, malaria, Zika, Leishmania, and hookworm, and 11 (25.58%) were applied to cancer, respectively. The International Journal of Environmental Research and Public Health (16.67%) and Spatial and Spatio-temporal Epidemiology (11.63%) were the most common journals where the articles were published (Table 1).

Data Source, Study Design, and Unit of Analysis
Seven studies (26.6%) used data from a national health survey or the Demographic and Health Survey (DHS) [43,46,51,60,64,72,75], while seven other studies (26.6%) used data from cancer registries. An estimated 11.63% of the studies used data from multiple surveys. More than one-third (41.86%) and 14 articles (32.56%) were ecological and crosssectional studies, respectively. The majority (55.81%) of the studies analysed two outcomes simultaneously in the model. Twenty-one studies reported the prevalence/incidence of outcomes, and nearly half of them had a prevalence of less than 10%.

Data Source, Study Design, and Unit of Analysis
Seven studies (26.6%) used data from a national health survey or the Demographic and Health Survey (DHS) [43,46,51,60,64,72,75], while seven other studies (26.6%) used data from cancer registries. An estimated 11.63% of the studies used data from multiple surveys. More than one-third (41.86%) and 14 articles (32.56%) were ecological and crosssectional studies, respectively. The majority (55.81%) of the studies analysed two outcomes simultaneously in the model. Twenty-one studies reported the prevalence/incidence of outcomes, and nearly half of them had a prevalence of less than 10%.

Data Source, Study Design, and Unit of Analysis
Seven studies (26.6%) used data from a national health survey or the Demographic and Health Survey (DHS) [43,46,51,60,64,72,75], while seven other studies (26.6%) used data from cancer registries. An estimated 11.63% of the studies used data from multiple surveys. More than one-third (41.86%) and 14 articles (32.56%) were ecological and crosssectional studies, respectively. The majority (55.81%) of the studies analysed two outcomes simultaneously in the model. Twenty-one studies reported the prevalence/incidence of outcomes, and nearly half of them had a prevalence of less than 10%.
Given that spatial analysis can be conducted at different spatial scales, about 11 studies performed analyses at the provincial level, and seven studies performed analyses at the county level. The mean number of spatial units was 428, ranging from 11 to 3577 spatial units. Among 43 joint spatial studies considered for the systematic review, 17 employed a joint spatiotemporal model. Of them, the vast majority (n = 15, 88.23%) used year as the unit of analysis. The mean temporal period was 21.71, ranging from 5 to 260 (Table 2).

Items Number Percentage (%) References
Types of spatial data

Key Implications of Applying Joint Spatial Modelling, Findings, and Methodological Gaps
The justifications provided in the included studies for fitting the joint spatial and spatiotemporal model, shared component spatial and spatiotemporal model, or multivariate spatial and spatiotemporal model varied. Fourteen (38.89%) studies applied the shared component spatial and spatiotemporal model to consider the spatial dependence of interrelated outcome variables and to better explore their overlapping epidemiology [43,44,47,52,57,58,[60][61][62]65,72,73,75,78]. Of the 36 joint spatial and spatiotemporal studies, 12 (33.33%) studies used the joint model for the ease of interpretation and to improve the precision of estimation [38,39,43,44,49,[53][54][55]57,60,65,76]. Different reasons were provided for using the joint spatial and spatiotemporal analysis. Nine studies (25%) applied the joint spatial and spatiotemporal model to borrow strength between diseases and to incorporate data from a more common and related disease when interest was in a relatively rare disease, thereby strengthening the relevant results of the rare disease [36,37,42,44,49,53,55,60,61].
The studies included in this systematic review have self-reported methodological gaps. Seven studies acknowledged that in aggregated data, ecological fallacies are introduced, and some relevant information may be concealed by using large geographical units of study [40,[47][48][49]53,65,80]. Therefore, using smaller units of analysis as a methodological gap may be a preferred approach. Four of the studies revealed that a meaningful number of temporal units is required to efficiently detect the temporal effect [38,41,44,80], and assuming the shared and specific components as independent ignores the possibility of interactions between the true covariates [38,44,62,73]. Three of the studies reported that MCMC has computational problems, model fitting, and convergence issues [42,43,56] ( Table 4). Table 4. Summary of the purpose of fitting joint spatial model, key findings, and reported methodological gaps of selected studies.

Items Number Percentage (%) References
Reasons for using joint modelling (n = 36) To borrow strength between diseases and to incorporate data from a more common and related disease when interest is in a relatively rare disease strengthens the relevant results of the rare disease
Joint spatial and spatiotemporal modelling of two or more cancers has become increasingly frequent over the past few decades [24,29,84] to examine the shared and differing trends of cancers regarding geographic patterns and shared risk factors. Contrary to univariate analysis, joint spatial models include shared components as various groups of cancers share common risk factors [85,86]. The majority of studies in spatial studies were based on single health outcomes, even though diseases such as cancer have common risk factors. Joint spatial and spatiotemporal modelling has recently become popular. It can model rare and common cancers to improve the estimates by borrowing strengths. It is feasible to modify behaviour common to cancers, and this has huge potential for preventing cancer. Major cancer risk factors can be altered by applying behavioural strategies including ceasing smoking, getting more exercise, managing weight, improving diet, limiting alcohol consumption, getting regular cancer screening tests, and limiting sun exposure [87]. Given that many of these cancer-preventive techniques lower the risk of many cancers, these might be supported by generating evidence through shared spatial and spatiotemporal studies.
The number of publications on joint spatial studies decreased dramatically during the coronavirus disease (COVID-19) period. It could be explained by the fact that since 2019, journals have prioritised COVID-19 research since little was known about the disease [88]. Evidence on the route of transmission, typical clinical features, underlying risk factors, and pathogenesis was limited, which is why the spatial analysis of COVID-19 with other infectious diseases was infrequent. In addition, advanced statistical approaches such as joint spatial and spatiotemporal models were not commonly applied as little was known about COVID-19's clinical manifestation, route of transmission, etiology, and treatment.
More than one-third of the studies used ecological data for analysis. This demonstrates how disease mapping frequently used aggregated data at a specific geographic level to generate area-level estimates to guide healthcare decisions and effective resource allocation [89]. Compared with other studies, ecological data can be accessed or retrieved from reports quite easily. The main data sources were DHS or national health survey data [43,46,51,60,64,72,75] and cancer registry data [35][36][37][38][39]41,42]. This could be because of DHS and other national health surveys having location data (GPS) and geolocated covariates including environmental, pollution, demographic, and socio-economic covariates available in these surveys [90].
Some of the studies were exploratory spatial analyses, including joint spatial autocorrelation techniques such as Moran Index (MI) statistics, Local Indicator Spatial Analysis (LISA), Getis Ord Gi * statistics, or Kulldroff spatial or spatiotemporal scan statistical tests [34,72,74,[76][77][78][79]. Unlike univariate spatial autocorrelations, multivariate spatial autocorrelations can determine how interrelated health outcomes such as TB and HIV, HSV-2 and HIV, and cancers influence each other spatially or spatiotemporally. They can explore the overlapping spatial and/or temporal distribution of two or more interrelated diseases. To detect the local clusters (hotspot and cold spot areas), spatial and space-time scanning statistical analysis and Getis-Ord Gi statistic are commonly used for cluster detections [91]. The hotspot and coldspot cluster detection is sensitive to the change in the size of the spatial and temporal units of analysis in which the data are aggregated; thus, analysing at a small spatial scale is preferable to identify hotspot areas efficiently. Multivariate spatial autocorrelation methods can investigate the spatial dependence of two or more interrelated health outcomes. In multivariate spatial autocorrelation analysis, the presence of a disease in a particular area is not only influenced by the prevalence of diseases in neighbouring areas, but also influenced by the presence of diseases that are related to one another in the area. However, they are unable to look at how the existence of one health outcome in one area may affect the spread of someone else in nearby areas, how one health outcome may affect the spread of another in adjacent areas, or, further, how it is affected by the spatial risk factors [92].
When variables were accounted for, the majority of the joint spatial and spatiotemporal models showed model improvement [32,34,37,38,43,[45][46][47][48][51][52][53][54][55]57,[59][60][61]64,73,75,79,80]. If the covariates were considered, joint spatial and spatiotemporal models performed much better than others and offered more insights than univariate models [93]. When the outcome is rare, the joint spatial models could improve the model performance by borrowing strength from interrelated diseases, neighbourhood areas, and/or time [10]. Besides, they can capture the spatial and spatiotemporal effects unexplained by the observed covariates by introducing random effects. This is well-suited to the data that have limited predictors and models that capture few covariates.
Only five studies underwent standardisation for common demographic variables [38,40,48,55,61]. Demographic covariates such as age, sex, race, etc. are the most obvious risk factors for almost all health and health-related conditions. In the univariate spatial analysis, the difference in incidence or mortality of given diseases because of age, sex, or other variables can be addressed by standardisation, which ignores the effects of these variables in the analysis. Age, sex, and other demographic variables are determinants for a multitude of health-related and other outcomes. The impacts of age, sex, and other common risk factors should be estimated using these variables in the model, as well as their interactions with other factors, such as spatial and temporal effects.
The vast majority of the studies were conducted based on a Bayesian modelling framework [32,33,73,79]. This was in line with the advancement of Bayesian statistics in the field of disease mapping. Bayesian spatial models are now commonly applied because of advances in knowledge of advanced statistics, programming, coding skills, and improved computer power resources that overcome computational problems. This methodology can provide more reliable area-level estimates specifically when the cases are rare or the population is small [90]. It can smoothen the observed extreme estimates towards the global mean or the neighbourhood values and provide more robust estimates, specifically when studied areas have sparse populations, by assigning prior distributions to define spatial structure to ensure the closer areas have more contribution than distant areas [22,94].
The majority of the Bayesian spatial and spatiotemporal studies undertook model fitting and inference based on MCMC [32,33,[35][36][37][38][39][40][41][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60] followed by INLA [42,43,[61][62][63][64][65]. The most popular models for spatial studies were Bayesian hierarchical models with structured and unstructured random effects. The MCMC approach is computationally very demanding and has convergence issues [95]. This is especially true for hierarchical models, which by their very nature make MCMC convergence unpredictable and slow. No matter the model, it is necessary to assess the convergence of posterior samples because there is no guarantee that such models can be easily fitted, in which additional simulations and model simplification will be necessary. In contrast, INLA has recently emerged as a reliable alternative method for fitting Bayesian spatial models that overcomes MCMC's drawbacks [96,97]. The package uses INLA to estimate Bayesian models without the need for posterior sampling techniques. In practical terms, numerical integration is performed via this approximation; therefore, it does not require a lot of iterative processing [98]. Bayesian estimation utilising the INLA methodology typically takes substantially less time than MCMC, which is the reason why this package was developed for spatial statistics in the first place.
The majority of studies applied joint spatial models [36,[38][39][40]42,43,[45][46][47][48][49][50][51][52]54,56,57,[59][60][61][62]64,73,75] followed by joint spatiotemporal models [32,33,35,37,41,44,55,58,63,65,80]. Unlike spatial autocorrelation studies, the spatial and spatiotemporal models can investigate the covariates influencing the distribution of the outcome across space and/or time [33,37,42,43,53,57,60,65,75,80]. The joint spatial and spatiotemporal models are applied for spatial analysis of two or more health outcomes or one outcome across different population groups. They typically rely on Generalised Linear Mixed Models (GLMMs) that consider shared and specific spatial, temporal, and spatiotemporal random effects. The data type, distribution, nature, and incidence of the outcomes determine the type of joint spatial and spatiotemporal models, e.g., poisson or negative binomial GLMM for count data and logistic regression for categorical outcomes. The dependence between diseases with similar spatial or temporal patterns is captured by prior distributions [13]. The multivariate CAR models smoothen noisy estimates and leverage information from nearby areas and interrelated diseases to predict spatially autocorrelated area-level disease risks [99]. However, multivariate CAR models are unable to show how the correlation of outcomes varies across space. The Copula geoadditive model overcomes this limitation and can demonstrate the change in the association between outcomes across geographic locations [100]. However, it is unable to detect the geographic areas contributing to higher or lower risk of simultaneous occurrence of multiple outcomes. Recently, the majority of studies applied the shared spatial and spatiotemporal models that decompose the spatial effect into shared and disease-specific spatial effects [24].
This review pointed out several recommendations for the development of improved joint spatial and spatiotemporal models. Some studies acknowledged that when data are aggregated, ecological fallacies are introduced, and some relevant information may be concealed by using large geographical units of study [40,[47][48][49]53,65,80]. Thus, using smaller units of analysis increases the precision of the estimates. Assuming the shared and specific components as independent denies the possibility of interactions between the true covariates [38,41,44,80], and a relevant number of temporal units is necessary to efficiently identify the temporal effect [38,44,62,73]. Moreover, some of the studies showed that MCMC has computational problems, model fitting, and convergence issues [42,43,56]. INLA is better for developing statistical models to obtain efficient risk estimations and direct the efficient distribution of medical interventions.
DIC was the most commonly used model comparison criteria to measure and compare the model goodness of fit and model complexity [32,33,[37][38][39]41,43,44,[46][47][48][49][52][53][54][60][61][62][63][64]73,78], followed by RMSPE [35,37,44,58,59,78]. Most of the studies used a combination of goodnessof-fit measures for model assessment. DIC and WAIC are model performance measures that are calculated by combining the model likelihood function (deviance (-2LLR)) and a model complexity term (number of effective parameters) [102,103]. In addition to model performance assessment, model accuracy assessment such as mean absolute prediction error and mean square prediction error are considered for model comparison [104]. Moreover, local measures of fit such as the conditional predictive ordinate are also used for making a model comparison. Apart from three studies [34,44,77], most studies reported maps for the visualisation of risk estimates. Maps offer epidemiologists enough evidence to display spatial risk and/or risk factors across time and/or space. It can give decision-makers motivation, insight, and assist potential health interventions in high-risk areas.
The scientific community may benefit from the epidemiological and statistical insights this systematic study provides in terms of joint spatial and spatiotemporal model applications in health research. First, the utility of joint spatial and spatiotemporal models is more pronounced in large data registries and when multiple interrelated diseases are fitted simultaneously. It is therefore crucial to estimate the smoothed relative risk of lower-prevalence cases through the borrowing of strength from the related cases and neighbourhood areas. Although there were joint spatial and spatiotemporal studies, the systematic review found heterogeneity in methods of estimation technique, statistical models, prior selections, defining adjacencies, and model complexities. This showed that a consistent framework for undertaking joint spatial and spatiotemporal models is needed. This framework is currently a focus of our research program. This systematic review provides insight suggesting that jointly modelling two or more cases that have shared characteristics is better to detect clusters of cases specifically when the number of cases is rare, such as in rare cancers and orphan diseases, or when the population is small. Another important finding was that the most complex models (joint spatial and spatiotemporal models incorporating covariates and interaction) performed very well. Overall, the systematic review identified several areas of improvement in joint spatial studies such as providing data, maps, scripts, and methodological gaps.
This review has some strengths and limitations, including an extensive search of six electronic databases to retrieve studies in an area without a previous systematic review. Careful title/abstract and full-text screening was carried out with predefined inclusion and exclusion criteria. One of the limitations of this review was that the majority of the studies were from a few countries, which might be because of spatial data being limited (GPS data and software packages for joint spatial and spatiotemporal model) or insufficient funding or statistical skills, indicating the possibility of publication bias or a focus of research effort on countries included in research publications. Another limitation was that only articles published in English were considered, so we may have excluded valuable contributions.

Conclusions
Multivariate disease mapping is crucial for understanding the burden of interrelated health outcomes over space and/or time. Numerous joint spatial and spatiotemporal methodologies aiming to explore the spatial risk of two or more health outcomes simultaneously were reviewed. The majority of studies used Bayesian methods, which handled a wider range of variance components at different levels in the model and could consider model uncertainties to provide reliable estimates. The most often utilised covariates in joint spatial and spatiotemporal models were socio-economic and demographic. Most of the reviewed studies used shared component spatial and spatiotemporal models with a Poisson-based and negative binomial modelling approach. Relatively few studies have been published on the applications of joint spatial and spatiotemporal models since the COVID-19 pandemic. Reviewed studies have acknowledged that aggregated data are liable to ecological fallacies and some relevant information may be concealed by using large geographical units of study. Therefore, this systematic review highlighted the need for future joint spatial and spatiotemporal models to analyse correlated health outcomes to guide decision-making for effective prevention and control strategies.
Supplementary Materials: The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/ijerph20075295/s1, Table S1: Databases and searching terms; Table S2: Risk of bias tool for assessment; Table S3: List of publications included for systematic review. Funding: This systematic review did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: This study is a systematic review of already published articles and the extracted data are provided as a supplementary file.