Author Response: Model-based Spatial-Temporal Mapping of Opisthorchiasis in Endemic Countries of Southeast Asia
Tingting Zhao,Yanping Feng,Pham Ngoc Doanh,Somphou Sayasone,Virak Khieu,Choosak Nithikathkul,Men‐Bao Qian,Yuantao Hao,Ying‐Si Lai
DOI: https://doi.org/10.7554/elife.59755.sa2
2021-01-01
Abstract:Article Figures and data Abstract Introduction Results Discussion Materials and methods Data availability References Decision letter Author response Article and author information Metrics Abstract Opisthorchiasis is an overlooked danger to Southeast Asia. High-resolution disease risk maps are critical but have not been available for Southeast Asia. Georeferenced disease data and potential influencing factor data were collected through a systematic review of literatures and open-access databases, respectively. Bayesian spatial-temporal joint models were developed to analyze both point- and area-level disease data, within a logit regression in combination of potential influencing factors and spatial-temporal random effects. The model-based risk mapping identified areas of low, moderate, and high prevalence across the study region. Even though the overall population-adjusted estimated prevalence presented a trend down, a total of 12.39 million (95% Bayesian credible intervals [BCI]: 10.10–15.06) people were estimated to be infected with O. viverrini in 2018 in four major endemic countries (i.e., Thailand, Laos, Cambodia, and Vietnam), highlighting the public health importance of the disease in the study region. The high-resolution risk maps provide valuable information for spatial targeting of opisthorchiasis control interventions. Introduction End of the epidemics of neglected tropical diseases (NTDs) by 2030 embodied in the international set of targets for the sustainable development goals (SDGs) endorsed by the United Nations empowers the efforts made by developing countries to combat the NTD epidemics (UN, 2015). To date, 20 diseases have been listed as NTDs, and opisthorchiasis is under the umbrella of food-borne trematodiasis (Ogorodova et al., 2015). Two species of opisthorchiasis are of public health significance, that is, Opisthorchis felineus (O. felineus), endemic in eastern Europe and Russia, and Opithorchis viverrini (O. viverrini), endemic in Southeast Asian countries (Petney et al., 2013). The later species is of our interest in the current article. According to WHO’s conservative estimation, an overall disease burden due to opisthorchiasis was 188,346 disability-adjusting life years (DALYs) in 2010 (Havelaar et al., 2015). Fürst and colleagues estimated that more than 99% of the burden worldwide attribute to O. viverrini infection in Southeast Asia (Fürst et al., 2012). Five countries in Southeast Asia, Cambodia, Lao PDR, Myanmar, Thailand, and Vietnam, are endemic for opisthorchiasis, with an estimated 67.3 million people at risk (Keiser and Utzinger, 2005). It is well documented that chronic and repeated infection with O. viverrini leads to the development of fatal bile duct cancer (cholangiocarcinoma) (International Agency for Research on Cancer, 1994). The life cycle of O. viverrini involves freshwater snails of the genus Bithynia as the first intermediate host, and freshwater cyprinoid fish as the second intermediate host. Humans and other carnivores (e.g., cats and dogs), the final hosts, become infected by consuming raw or insufficiently cooked infected fish (Andrews et al., 2008; Saijuntha et al., 2014). Behavioral, environmental, and socioeconomic factors affect the transmission of O. viverrini (Grundy-Warr et al., 2012, Phimpraphai et al., 2017, Phimpraphai et al., 2018, Prueksapanich et al., 2018). Raw or insufficiently cooked fish consumption is the cultural root in endemic countries, showing a strong relationship with the occurrence of the disease (Andrews et al., 2008; Grundy-Warr et al., 2012). Poorly hygienic conditions increase the risk of infection, especially in areas practicing raw-fish-eating habit (Grundy-Warr et al., 2012). In addition, environmental and climatic factors, such as temperature, precipitation, and landscape, affecting either snail/fish population or growth of the parasites inside the intermediate hosts, can potentially influence the risk of human infection (Forrer et al., 2012; Suwannatrai et al., 2017). Important control strategies of O. viverrini infection include preventive chemotherapy, health education, environmental modification, improving sanitation, as well as comprehensive approaches with combinations of the above (Saijuntha et al., 2014). For purposes of public health control, WHO recommends implementing preventive chemotherapy once a year or once every 2 years depending on the levels of prevalence in population, with complementary interventions such as health education and improvement of sanitation (WHO, 2009). Understanding the geographical distribution of O. viverrini infection risk at high spatial resolution is critical to prevent and control the disease cost-effectively in priority areas. Thailand conducted national surveys for O. viverrini prevalence in 1981, 1991, 2001, 2009, and 2014 (Echaubard et al., 2016; Suwannatrai et al., 2018), but the results of these surveys were presented at the province level, which is less informative for precisely targeting control interventions. Suwannatrai and colleagues, based on climatic and O. viverrini presence data, produced climatic suitability maps for O. viverrini in Thailand using the MaxEnt modeling approach (Suwannatrai et al., 2017). The maps brought insights for identifying areas with a high probability of O. viverrini occurrence; however, they did not provide direct information on prevalence of O. viverrini in population (Elith et al., 2011). A risk map of O. viverrini infection in Champasack province of Lao PDR was presented by Forrer and colleagues (Forrer et al., 2012). To our knowledge, high-resolution, model-based risk estimates of O. viverrini infection are unavailable in the whole endemic region of Southeast Asia. Bayesian geostatistical modeling is one of the most rigorous inferential approaches for high-resolution maps depicting the distribution of the disease risk (Karagiannis-Voules et al., 2015). Geostatistical modeling relates geo-referenced disease data with potential influencing factors (e.g., socioeconomic and environmental factors) and estimates the infection risk in areas without observed data (Gelfand and Banerjee, 2017). Common geostatistical models are usually based on point-referenced survey data (Banerjee et al., 2014). In practice, disease data collected from various sources often consists of point-referenced and area-aggregated data. Bayesian geostatistical joint modeling approaches provide a flexible framework for combining analysis of both kinds of data (Moraga et al., 2017; Smith et al., 2008). In this study, we aimed (1) to collect all available survey data on the prevalence of O. viverrini infection at point- or area-level in Southeast Asia through systematic review; and (2) to estimate the spatial-temporal distribution of the disease risk at a high spatial resolution, with the application of advanced Bayesian geostatistical joint modeling approach. Results A total of 2690 references were identified through systematically reviewing peer-review literatures, and 13 additional references were gathered from other sources. According to the inclusion and exclusion criteria, 168 records were included, resulted in a total of 580 ADM1-level surveys in 174 areas, 210 ADM2-level surveys in 142 areas, 53 ADM3-level surveys in 51 areas, and 251 point-level surveys at 207 locations in five endemic countries (i.e., Cambodia, Lao PDR, Myanmar, Thailand, and Vietnam) of Southeast Asia (Figure 1). Around 70% and 15% of surveys were conducted in Thailand and Lao PDR, respectively. Only two relevant records were obtained from Myanmar. To avoid large estimated errors, we did not include this data in the final geostatistical analysis. All surveys were conducted after 1970, with around 75% done after 1998. Most surveys (95%) are community based. Around 40% of surveys used the Kato–Katz technique for diagnosis, while another 42% did not specify diagnostic approaches. Mean prevalence calculated directly from survey data was 16.74% across the study region. A summary of survey data is listed in Table 1, and survey locations and observed prevalence in each period are shown in Figure 2. Area-level data cover all regions in Thailand and Lao PDR, and most regions in Cambodia and Vietnam, while point-referenced data are absent in most areas of Vietnam, the western part of Cambodia and southern part of Thailand. Around 70% of eligible literatures got a score equal or more than 7, indicating an overall good quality of eligible literatures in our study (Figure 2—figure supplement 1). Figure 1 with 1 supplement see all Download asset Open asset Data search and selection flow chart. Figure 2 with 1 supplement see all Download asset Open asset Survey locations and observed prevalence of O. viverrini infection in endemic countries of Southeast Asia. (A) 1978–1982, (B) 1983–1987, (C) 1988–1992, (D) 1993–1997, (E) 1998–2002, (F) 2003–2007, (G) 2008–2012, and (H) 2013–2018. Figure 2—source data 1 The original data of O. viverrini infection in endemic countries of Southeast Asia. https://cdn.elifesciences.org/articles/59755/elife-59755-fig2-data1-v2.xlsx Download elife-59755-fig2-data1-v2.xlsx Figure 2—source data 2 The results of the preferential sampling test. https://cdn.elifesciences.org/articles/59755/elife-59755-fig2-data2-v2.docx Download elife-59755-fig2-data2-v2.docx Table 1 Overview of opisthorchiasis survey data in Southeast Asia. CambodiaLao PDRMyanmarThailandVietnamTotalRelevant papers144329715168Total surveys/locations91/73156/996/6770/33571/641094/574Survey type (surveys/locations) School33/314/40/013/130/050/48 Community58/46152/946/6757/32571/641044/535Location type (surveys/locations) Point-level55/4363/513/3125/1055/5251/207 ADM3-level0/00/00/053/510/053/51 ADM2-level14/1135/270/0159/1022/2210/142 ADM1-level22/1958/183/3433/7764/57580/174Period1998–20161989–20162015–20161978–20181991–20151978–2018Year of survey (surveys/locations) 1978–19820/00/00/0123/1150/0123/115 1983–19870/00/00/07/60/07/6 1988–19920/02/20/097/891/1100/92 1993–19970/09/50/018/186/233/25 1998–200225/2228/220/0103/1032/2158/149 2003–20073/226/240/015/151/145/42 2008–201262/4875/540/0166/1539/8312/263 2013–20181/116/166/6241/20152/52316/276Diagnostic methods (surveys/locations) Kato–Katz86/70128/833/3212/1667/7436/329 FECT2/28/73/3109/990/0122/111 Stoll’s0/00/00/038/280/038/28 PCR0/05/40/01/10/06/5 Combined3/314/130/014/120/031/28 Others0/01/10/06/60/07/7 NS*0/05/50/0391/11164/57460/173Mean prevalence10.56%39.50%4.93%14.25%2.65%16.74% *NS: not stated or missing. Seven variables were selected for the final model through the Bayesian variable selection process (Table 2). The infection risk was 2.61 (95% BCI: 2.10–3.42) times in the community as much as that in school-aged children. Surveys using FECT (formalin-ethyl acetate concentration technique) as the diagnostic method showed a lower prevalence (OR 0.76, 95% BCI: 0.61–0.93) compared to that using Kato–Katz method, while no significant difference was found between Kato–Katz and the other diagnostic methods. Human influence index and elevation were negatively correlated with the infection risk. Each unit increase of the HII index was associated with 0.01 (95% BCI: 0.003–0.02) decrease in the logit of the prevalence. And increase in 1 m in elevation was associated with the 0.003 (95% BCI: 0.001–0.005) decrease in the logit of the prevalence. The spatial range was estimated as 83.55 km (95% BCI: 81.34–86.61), the spatial variance σϕ2 was 12.59 (95% BCI: 11.96–13.56), the variance of beta-likelihood σβ2 was 0.15 (95% BCI: 0.14–0.15), and the temporal correlation coefficient ρ was 0.66 (95% BCI: 0.65–0.67). Model validation showed that our model was able to correctly estimate 79.61% of locations within the 95% BCI, indicating the model had a reasonable capacity of prediction accuracy. The ME, MAE, and MSE were 0.24%, 9.06%, and 2.38%, respectively, in the final model, while they were −7.14%, 16.67%, and 5.09%, respectively, in the model only based on point-referenced data, suggesting that the performance of the final model was better than the model only based on point-referenced data. On the other hand, Monte Carlo test for preferential sampling suggested that preferential sampling may exist for survey locations in one third (6/18) of the survey years (Figure 2—source data 2). Table 2 Posterior summaries of model parameters. Estimated median (95% BCI)ORProb (%)*Intercept−4.51 (−5.08, –3.94)Survey type School-based surveyRefRef- Community-based survey0.96 (0.70, 1.23)2.61 (2.10, 3.42)>99.99Diagnostic methods Kato–KatzRefRef- FECT−0.28 (–0.49, –0.07)0.76 (0.61, 0.93)0.80 Other methods0.01 (–0.07, 0.10)1.01 (0.93, 1.12)64.20Land surface temperature (LST) in the daytime (°C) <30.65RefRef- 30.65–32.070.25 (–0.001, 0.50)1.28 (0.999, 1.65)97.40 >32.070.07 (–0.18, 0.33)1.07 (0.84, 1.39)73.40Human influence index−0.01 (–0.02, –0.003)0.99 (0.98, 1.00)0.80Distance to the nearest open water bodies (km)0.24 (–1.45, 1.94)1.27 (0.23, 6.96)60.20Elevation (m)−0.003 (–0.005,–0.001)0.997 (0.995, 0.999)<0.01Travel time to the nearest big city (min)0.0001 (–0.002, 0.002)1.00 (0.998, 1.002)56.60 *Posterior probability of OR > 1. The estimated risk maps of O. viverrini infection in different selected years (i.e., 1978, 1983, 1988, 1993, 1998, 2003, 2008, 2013, and 2018) are presented in Figure 3. In 2018, the high infection risk (with prevalence >25%) was mainly estimated in regions of the southern, the central, and the north-central parts of Lao PDR, some areas in the east-central parts of Cambodia, and some areas of the northeastern and the northern parts of Thailand. The southern part of Thailand, the northern part of Lao PDR, and the western part of Cambodia showed low risk estimates (with prevalence <5%) of O. viverrini infection. The central and several southern parts of Vietnam showed low to moderate risk of O. viverrini infection, while there was no evidence of O. viverrini in other parts of Vietnam. High estimation uncertainty was mainly present in the central part of Lao PDR, the northern and the eastern parts of Thailand, and the central part of Cambodia and Vietnam (Figure 4). Figure 3 with 1 supplement see all Download asset Open asset Model-based estimated risk maps of O. viverrini infection in endemic countries of Southeast Asia in different years. Estimated prevalence based on the median of the posterior estimated distribution of infection risk in (A) 1978, (B) 1983, (C) 1988, (D) 1993, (E) 1998, (F) 2003, (G) 2008, (H) 2013, and (I) 2018. Figure 3—source data 1 The sensitivity analysis results of model-based estimated risk maps in 2018. https://cdn.elifesciences.org/articles/59755/elife-59755-fig3-data1-v2.xlsx Download elife-59755-fig3-data1-v2.xlsx Figure 4 Download asset Open asset The estimation uncertainty in endemic countries of Southeast Asia in different years. (A) 1978, (B) 1983, (C) 1988, (D) 1993, (E) 1998, (F) 2003, (G) 2008, (H) 2013, and (I) 2018. Figure 4—source data 1 The results of the estimated uncertainty in endemic countries of Southeast Asia in different years. https://cdn.elifesciences.org/articles/59755/elife-59755-fig4-data1-v2.xlsx Download elife-59755-fig4-data1-v2.xlsx In addition, the infection risk varies over time across the study region (Figure 5). Areas of northern Thailand showed an increasing trend in periods 1978–1988 and 1993–2003, while most areas of the country presented a considerable decrease of infection risk after 2008. The infection risk first increased and then decreased in areas of the north, the central, and the southern parts of Lao PDR and the central parts of Vietnam. The east-central and western part of Cambodia showed an increasing trend in recent years. Figure 5 Download asset Open asset Changes of O. viverrini infection risk across time periods. Changes were calculated by the median of the posterior estimated distribution of infection risk for the latter time period minus that for the former time period divided by that for the former time period. The risk changes (A) between 1978 and 2018; (B) between 1978 and 1983; (C) between 1983 and 1988; (D) between 1988 and 1993; (E) between 1993 and 1998; (F) between 1998 and 2003; (G) between 2003 and 2008; (H) between 2008 and 2013; and (I) between 2013 and 2018 (source data: Figure 5—source data 1). Figure 5—source data 1 The results of the changes of O. viverrini infection risk across time periods. https://cdn.elifesciences.org/articles/59755/elife-59755-fig5-data1-v2.xlsx Download elife-59755-fig5-data1-v2.xlsx The population-adjusted estimated prevalence over the study region presents a trend down after 1995 (Figure 6 and Figure 6—figure supplements 1–9). At the country level, the estimated prevalence in Thailand showed a fast decline after 1995 and took on a gradually decreasing change in Cambodia. In Lao PDR, the overall prevalence maintained quite stable before 1990 and decreased slightly between 1990 and 1997, increased significantly after 1997, then decreased from 2006, and became stable after 2011. The prevalence is stable in Vietnam during the whole study period. We estimated that in 2018, the overall population-adjusted estimated prevalence of O. viverrini infection in the whole study region was 6.57% (95% BCI: 5.35–7.99%), corresponding to 12.39 million (95% BCI: 10.10–15.06) infected individuals (Table 3). Lao PDR showed the highest prevalence (35.21%, 95% BCI: 28.50–40.70%), followed by Thailand (9.71%, 95% BCI: 7.98–12.17%), Cambodia (6.15%, 95% BCI: 2.41–11.73%), and Vietnam (2.15%, 95% BCI: 0.73–4.40%). Thailand had the largest numbers of individuals estimated to be infected with O. viverrini (6.71 million, 95% BCI: 5.51–8.41), followed by Lao PDR (2.45 million, 95% BCI: 1.98–2.83), Vietnam (2.07 million, 95% BCI: 0.70–4.24), and Cambodia (1.00 million, 95% BCI: 0.39–1.90). Figure 6 with 9 supplements see all Download asset Open asset Trends in estimated prevalence of O. viverrini infection in Southeast Asia. Figure 6—source data 1 The results of the estimated prevalence of O. viverrini infection in Southeast Asia in different years. https://cdn.elifesciences.org/articles/59755/elife-59755-fig6-data1-v2.xlsx Download elife-59755-fig6-data1-v2.xlsx Table 3 Population-adjusted estimated prevalence and number of individuals infected with O. viverrini in endemic countries of Southeast Asia in 2018*. Population (×103)Prevalence (%)No. infected (×103)Cambodia16227.396.15 (2.41, 11.73)997.95 (390.46, 1903.46)Lao PDR6960.2835.21 (28.50, 40.70)2450.54 (1983.38, 2832.96)Thailand69112.649.71 (7.98, 12.17)6708.68 (5514.87, 8411.98)Vietnam96421.692.15 (0.73, 4.40)2073.72 (703.46, 4244.85)Total188722.016.57 (5.35, 7.98)12389.69 (10099.29, 15060.18) *Estimates were based on gridded population of 2018 and the median and 95% BCI of the posterior estimated distribution of the infection risk in 2018. Discussion In this study, we produced model-based, high-resolution risk estimates of opisthorchiasis across endemic countries of Southeast Asia. The disease is the most important foodborne trematodiasis in the study region (Sripa et al., 2010), taking into account most of the disease burden of opisthorchiasis in the world (Fürst et al., 2012). The estimates were obtained by systematically reviewing all possible geo-referenced survey data and applying a Bayesian geostatistical modeling approach that jointly analyzes point-referenced and area-aggregated disease data, as well as environmental and socioeconomic predictors. Our findings will be important for guiding control and intervention cost-effectively and serve as a baseline for future progress assessment. Our estimates suggested that there was an overall decrease of O. viverrini infection in Southeast Asia from 1995 onwards, which may be largely attributed to the decline of infection prevalence in Thailand. This decline was probably on account of the national opisthorchiasis control program launched by the Ministry of Public Health of Thailand from 1987 (Jongsuksuntigul and Imsomboon, 2003; Jongsuksuntigul et al., 2003). Our high-resolution risk estimates in Thailand in 2018 showed similar pattern as the climatic suitability map provided by Suwannatral and colleagues (Suwannatrai et al., 2017). In this case, we estimated the prevalence of the population instead of the occurrence probability of the parasite, which arms decision makers with more direct epidemiological information for guiding control and intervention. The national surveys in Thailand reported a prevalence of 8.7% and 5.2% in 2009 and 2014, respectively (MOPH, 2014; Wongsaroj et al., 2014). However, we estimated higher prevalence of 12.44% (95% BCI: 10.79–14.26%) and 9.34% (95% BCI: 7.88–11.02%) in 2009 and 2014, respectively. Even though the national surveys covered most provinces in Thailand, estimates were based on simply calculating the percentage of positive cases among all the participants (Wongsaroj et al., 2014), and the remote areas might not be included (Maipanich et al., 2004). Instead, our estimates were based on rigorous Bayesian geostatistical modeling of available survey data with environmental and socioeconomic predictors, accounting for heterogeneous distribution of infection risk and population density when aggregating country-level prevalence. Our findings suggested that the overall prevalence of O. viverrini remained high (>20%) in Lao PDR during the study periods, consistent with conclusions drawn by Suwannatrai et al., 2018. We estimated that a total number of 2.45 (95% BCI: 1.98–2.83) million people living in Lao PDR were infected with O. viverrini, equivalent to that estimated by WHO in 2004 (WHO, 2002). Besides, our risk mapping for Champasack province shares similarly risk map pattern produced by Forrer and colleagues (Forrer et al., 2012). A national-scaled survey in Cambodia during the period 2006–2011 reported infection rate of 5.7% (Yong et al., 2014), lower than our estimation of 8.34% (95% BCI: 5.25–14.95%) in 2011. The former may underestimate the prevalence because more than 77% of participants were schoolchildren (Yong et al., 2014). Another large survey in five provinces of Cambodia suggested a large intra-district variation, which makes the identification of endemic areas difficult (Miyamoto et al., 2014). Our high-resolution estimates for Cambodia help to differentiate the intra-district risk. However, the estimates should be taken cautious due to large district-wide variances and a relatively small number of surveys. Indeed, O. viverrini infection was underreported in Cambodia (Khieu et al., 2019), and further point-referenced survey data are recommended for more confirmative results. Although an overall low prevalence was estimated in Vietnam (2.15%, 95% BCI: 0.73–4.40%) in 2018, it corresponds to 2.07 million (95% BCI: 0.70–4.24 million) people infected, comparable to the number in Lao PDR, mainly due to a larger population in Vietnam. The risk mapping suggested moderate to high risk areas presented in central Vietnam, with a high risk in Phu Yen province for many years, particularly. This agreed with previous studies considering the province a ‘hotspot’ (Doanh and Nawa, 2016). Of note, even though there was no evidence of O. viverrini infection in the northern part of the country, Clonorchis sinensis, another important liver fluke species, is endemic in the region (Sithithaworn et al., 2012). We did not provide estimates for Myanmar in case of large estimated errors. Indeed, only two relevant papers were identified by our systematic review, where one shows low to moderate prevalence in three regions of Lower Myanmar (Aung et al., 2017), and the other found low endemic of O. viverrini infection in three districts of the capital city Yangon (Sohn et al., 2019). Nation-wide epidemiological studies are urgent for a more comprehensive understanding of the disease in Myanmar. We identified several important factors associated with O. viverrini infection in Southeast Asia, which may provide insights for the prevention and control of the disease. The infection risk was higher in the entire community than that in schoolchildren, consistent with multiple studies (Aung et al., 2017; Forrer et al., 2012; Miyamoto et al., 2014; Van De, 2004; Wongsaroj et al., 2014). A negative association was found between O. viverrini infection and elevation, suggesting the disease was more likely to occur in low altitude areas, which was consistent with a previous study (Wang et al., 2013). HII, a measure of human direct influence on ecosystems (Sanderson et al., 2002), showed a negative relationship with O. viverrini infection risk, indicating the disease was more likely to occur in areas with low levels of human activities, which were often remote and economically underdeveloped. The habit of eating raw or insufficiently cooked fish was more common in rural areas than that in economically developed ones, which could partially explain our findings (Grundy-Warr et al., 2012, Keiser, 2019). Indeed, this culturally rooted habit is one of the determinants for human opisthorchiasis (Kaewpitoon et al., 2008; Ziegler et al., 2011). However, the precise geographical distribution of such information is unavailable and thus we could not use it as a covariate in this study. Our estimate of the number of people infected with O. viverrini is higher than that of the previous study (12.39 million vs 8.6 million [Qian and Zhou, 2019]) emphasizing the public health importance of this neglected disease in Southeast Asia, and suggesting that more effective control interventions should be conducted, particularly in the high risk areas. The successful experience in the intervention of Thailand may be useful for reference by other endemic countries of the region. The national opisthorchiasis control program, supported by the government of Thailand, applied interrelated approaches, including stool examination and treatment of positive cases, health education aiming at the promotion of cooked fish consumption, and environmental sanitation to improve hygienic defecation (Jongsuksuntigul and Imsomboon, 2003). In addition, for areas with difficulties to reduce infection risk, a new strategy was developed by Sripa and colleagues, using the EcoHealth approach with anthelminthic treatment, novel intensive health education on both communities and schools, ecosystem monitoring, and active participation of the community (Sripa et al., 2015). This ‘Lawa model’ shows good effectiveness in Lawa Lake area, where the liver fluke was highly endemic (Sripa et al., 2015). Furthermore, common integrated control interventions (e.g., combination of preventive chemotherapy with praziquantel, improvement of sanitation and water sources, and health education) are applicable not only for opisthorchiasis but also for other NTDs, such as soil-transmitted helminth infection and schistosomiasis, which are also prevalent in the study region (Dunn et al., 2016; Gordon et al., 2019). Implementation of such interventions in co-endemic areas could be cost-effective (Linehan et al., 2011; WHO, 2012). Frankly, several limitations exist in our study. We collected data from different sources, locations of which might not be random and preferential sampling may exist. We performed a risk-preferential sampling test and the results showed that preferential sampling might exist for survey locations in one third (6/18) of the survey years (Figure 2—source data 2). The corresponding impacts might include improper variogram estimator, biased parameter estimation, and unreliable exposure surface estimates (Diggle et al., 2010; Pati et al., 2011; Gelfand et al., 2012). To avoid a more complex model, we did not take into account the preferential sampling issue for our final model, as the model validation showed a reasonable capacity of prediction accuracy. However, the disadvantage of this issue should be well aware. We set clear criteria for selection of all possible qualified surveys and did not exclude surveys that reported prevalence in intervals without exact observed values. Sensitivity analysis showed that the using the midpoint values of the intervals had little effects on the final results (Figure 3—figure supplement 1). For surveys across a large area, complex designs, such as randomly sampling from subgroups of the population under a well-designed scheme, are likely adopted, as it is impractical to draw simple random samples from the whole area. In such case, respondents may have unequal probabilities to be selected, thus weighting should be used to generalize results for the entire area. The observed disease data we collected were from surveys either at point-level (i.e., community or school) or aggregated over areas. For point-level data, as study areas were quite small, simple sampling design was mostly used in the corresponding surveys. And for areal-level data, particularly those aggregated across ADM1, complex designs were likely applied. However, most of the corresponding surveys were only reported raw prevalence or prevalence without clarifying whether weighting was applied. Thus, we did not have enough information to address the design effect for each single survey included. On the other hand, as population density across the study region was different, we calculated the estimated country- and provincial-level prevalence by averaging the estimated pixel-level prevalence weighted by population density. In this way, we took into account the diversity of population density across areas for regional summaries of the estimates. We assumed similar proportions of age and gender in different surveys, as most of which only reported prevalence aggreg