Spatiotemporal Trends and Distributions of Malaria Incidence in the Northwest Ethiopia

Understanding and extracting noticeable patterns of malaria surveillance data at the district level are crucial for malaria prevention, control, and elimination progress. This study aimed to analyze spatiotemporal trends and nonparametric dynamics of malaria incidences in northwest Ethiopia, considering spatial and temporal correlations. The data were analyzed using count regression spatiotemporal models under the Bayesian setups, and parameters were estimated using integrated nested Laplace approximations (INLA). The region had a declining linear trend, and the average annual malaria incidence rate was 24.8 per 1,000 persons between 2012 and 2020. The malaria incidence rate was decreased by 0.984 (95% CI: 0.983, 0.986) per unit increase in months between July 2012 and June 2020. Districts found in the western and northwestern parts of the region had a steeper trend, while districts in the eastern and southern parts had a less steep trend than the average trend of the region. Compared to the regional level trend, the decreasing rate of malaria incidence trends was lower in most town administrations. The nonparametric dynamics showed that the monthly malaria incidence had a sinusoidal wave shape that varied throughout study periods. Malaria incidence had a decreasing linear trend changed across districts of the study region, and the steepness of trends of districts might not depend on incidences. Thus, an intervention and controlling mechanism that considers malaria incidences and district-specific differential trends would be indispensable to mitigate malaria transmission in the region.


Introduction
Malaria surveillance is a continuous, systematic collection, analysis, and interpretation of malaria data and is used as input for planning, implementation, and evaluation of public health practice [1,2]. Understanding the trend of malaria surveillance data and extracting its noticeable patterns of prevalence across districts are crucial to inform concerned organizations about the state of malaria control and elimination [3].
Even though malaria morbidity and mortality have significantly declined in Ethiopia since 2001, 9% and 8% of the global Plasmodium vivax cases occurred in 2017 and 2018, respectively [4]. According to the world malaria report, the estimated number of malaria cases was 4,231,328 in 2020, significantly declining compared to 2010 [5]. Malaria is the major public health challenge in the Amhara region of Ethiopia [6]. e average malaria incidence rate was 23.51 per 1000 persons in the region between 2004 and 2014, and it was the fifth-largest incidence in the country [7]. e region has a decreased malaria burden, and the estimated prevalence was 4.6% in 2006, 0.6% in 2007, and 0.8% in 2011 [8]. More than 1.1 million cases occurred in 2012, with an annual incidence rate of 60 per 1000 persons, accounting for 19% of the national malaria cases. A cross-sectional survey in the 19 districts of more vulnerable revealed that malaria prevalence was 1.9% in 2013, and 40% of confirmed cases were asymptomatic [9].
Malaria transmission exhibits spatial and temporal variability due to diverse ecology, climate, altitude, topography, and human settlement patterns [10]. Even under lower transmission settings, districts have substantial malaria incidence variations in the region. Studies revealed that malaria risks have spatial, temporal, and spatiotemporal heterogeneity in the south and north Gondar zones [3,6,9]. Malaria incidence has spatial variations across districts of the region, influenced by trends or the number of malaria cases in their neighborhood districts [11]. Exploring and estimating malaria trends at the district level in considering the space-time interaction effect would reflect the entire picture of malaria transmission that would be crucial for planning, intervention, elimination, and control. Previous studies did not cover districts in the region, both with lower and higher malaria transmission scenarios. e regional level trend, differential trends of districts, and nonparametric dynamics of malaria incidence were not studied with the account of space-time interaction effects through spatial and temporal autocorrelations [6,9]. In this study, malaria surveillance data from July 2012 and June 2020 were utilized to investigate and estimate global linear effects of time in a month, differential trends of districts, and nonparametric dynamics of malaria incidence.

Study Setting.
is study was conducted in the Amhara National Regional State (ANRS), Ethiopia ( Figure 1). e region is located in northwestern Ethiopia between 9°20′ and 14°20′ north latitude and 36°20′ and 40°20′ east longitude. e region has a monsoon climate, and its elevation ranges from 506 meters at the bottom of the Blue Nile gorge to 4,533 meters at Ras-Dajen mount [12]. e study area encompasses 152 rural and town districts that have had weekly malaria surveillance reports to the Amhara Public health institute (APHI) since 2012.

Study Design and Period.
Using the weekly reported malaria surveillance data of Amhara Public Health Institute (APHI) from July 2012 to June 2020, a retrospective study design was used to investigate monthly spatiotemporal trends and space-time dynamics of malaria incidence.

Data.
e weekly malaria surveillance data of rural districts, town administrations, general and specialized hospitals of the study region were obtained from Amhara public health institute (APHI), within the Public Health Emergency Management (PHEM) directorate. e weekly malaria data were reported using WHO EPI week that starts on Sunday and ends on Saturday.
e weekly malaria surveillance data encompass various information about districts, EPI-week, budget year, total malaria cases (clinical and confirmed), number of cases in different age groups, and types of Plasmodium species. is study used monthly malaria cases between July 2012 and June 2020, obtained by aggregating weekly surveillance data based on a month that had four or more days from a given WHO EPI-week. e estimated monthly population size of districts between 2012 and 2020, obtained from the central statistical agency of Ethiopia and Amhara National Regional State Planning Commission, is used as an offset in the model.

Inclusion and Exclusion Criteria.
e districts encompass various healthcare institutions such as health posts, health stations, and different levels of hospitals. Especially town administrations have general, referral, or specialized hospitals where patients get treatment with or without referral letters from the rural districts. Hence, malaria cases from the referral and specialized hospitals were not included in the analysis to preclude the overestimation of malaria incidence trends of town administrations.

Statistical
Analysis. Spatiotemporal disease mapping models are widely used in disease surveillance studies [13] when the interest is identifying the underlying spatial and temporal patterns of diseases. Count spatiotemporal models were used to investigate and estimate malaria incidence trends at the district level in considering space-time interaction effects. Let Y it be the number of malaria cases in the i th district at time t and n it be the size of the corresponding population at risk, for i � 1, 2, . . . , I � 152 and t � 1, 2, . . . , T � 96. Suppose μ c it be the conditional expectations of the Y it given random effects, then the disease mapping model is as follows: where μ is the mean log-number of cases in overall areas, u i and v i are structured and unstructured spatial effects, respectively, and T t represents temporal effects and can be specified in a parametric or nonparametric structure [14].

Parametric
Trend. e parametric time trend on the temporal components of the spatiotemporal model was introduced by Bernardinelli et al. [15], and the linear predictor given in equation (1) is written as follows: here, the equation includes the structured and unstructured spatial effects that represent between area-specific lognumber of cases and overall mean rate, β is mean linear time trend, which represents global time effect, and δ i is a differential trend (DT) that represents the difference between district-specific trends and regional level trends, which measures time and space interaction effects. For identifiability purposes, a sum-to-zero constraint is imposed on δ(δ 1 , δ 2 , . . . , δ n ) and u(u 1 , u 2 , . . . , u n ) and a value of δ i < 0 implies that an area-specific trend is less steep than a mean trend, whilst a value of δ i > 0 implies that the area-specific trend is steeper than the mean trend. e presence of a statistically significant differential trend in a given area is tested using posterior probability (PP), the Bayesian equivalent of the p-value [16]. A value of posterior probability, which is greater than 0.90, indicates some evidence of differential trend is greater than zero (greater than one in the exponential scales of parameters); that is, the area or districtspecific trend is as steep as or steeper than the mean trend [15]. In disease mapping, spatial effects were considered random [15], and a Bayesian approach estimates parameters through assignments of prior distributions to the fixed parameters and spatial effect u i . e prior distribution for the intercept is assumed to be zero-mean Gaussian with variance σ 2 ρ . e conditional autoregressive model proposed by Besag et al. [17] is used as a prior distribution for the spatially structured effects, that is, where μ i � (1/ j≠i w ij ) j≠i w ij u j , σ 2 i � σ 2 / j≠i w ij , and w ij � 1, if i and j are geographically adjacent districts and w ij � 0 otherwise [18]. e uncorrelated spatial effect prior is defined as v i ∼ N(0, σ 2 v ) and is used to allow for uncorrelated extra variation. e differential trend, which accounts spacetime interaction, is assumed to follow normal distributions with mean zero and variance σ 2 δ , that is, δ it ∼ N(0, σ 2 δ ). In addition, noninformative prior distributions are considered for hyperparameters, so that the inference is based on the assumed model and observed data [19].

Nonparametric Dynamic Trend.
Bernardinelli et al. [15] evaluated spatiotemporal interactions on the disease risk by imposing a linearity constraint on differential temporal trends. e temporal trends in disease risk may be different for different spatial locations with a noticeable spatial variation. However, disease risk might not have linear temporal trends across areas. Knorr-Held [20] extended the linear time trend using a dynamic nonparametric formulation for the linear predictors underlying inseparable spacetime effects.
e log-link of the mean number of cases where μ is the overall malaria risk level, u i and v i represent unspecified features of district i that, respectively, do and do not display spatial structure. Similarly, α t and c t represent the unspecified features of month t that, respectively, do and do not display temporal structure. e interaction between space and time effect is incorporated using δ it (δ 11 , . . . , δ nT ), which would explain differences in the time trend of malaria cases for different areas or districts. e parameter δ it is assumed to be Gaussian with precision matrix λ δ K δ , where λ δ is unknown scalar and K δ � I v ⊗ I c � I is a prespecified structure matrix that is Kronecker's product of the structured matrices of the main effects, which are assumed to interact [21]. Here, the unobserved covariate effect at each pixel or area (i, t) assumed does not have any structure in space-time interactions, and all interaction parameters are prior independent [14,20]. e spatial and interaction effects are assigned conditional autoregressive and normal prior distributions, respectively. In contrast, the structured temporal random effects are assigned random walk order 2 (RW2) prior distribution in order to account dynamic nature of the disease incidence. e random walk order 2 model is specified as follows: Moreover, hyperparameters are assigned noninformative flat priors.

Parameter Estimation and Model Comparison.
e Poisson spatiotemporal model is widely utilized for analyzing spatiotemporal count data. However, count data often display overdispersion, and using Poisson regression may underestimate the standard errors and overstate the Journal of Tropical Medicine significance of regression parameters [22]. In addition, a higher number of zeros in the data imposes problems, which might be occurred due to structural and sampling issues would be modeled using zero-inflated models [23,24]. Districts located in the highlands areas, altitudes with 2500 and above meters, have a lower malaria transmission intensity, are malaria free, and might have zero monthly malaria cases, especially in the lower transmission seasons [9,[25][26][27]. e spatiotemporal models were compared using the deviance information criterion (DIC) and Watanabe--Akaike information criterion (WAIC). A spatiotemporal model with the smallest DIC value is used to fit and interpret results. e model parameters and hyperparameters were estimated using Integrated Nested Laplace Approximations (INLA), which is a computationally efficient numerical approximation method for fitting complex spatiotemporal models and faster than Markov Chain Monte Carlo (MCMC) [28,29].

Temporal and Spatial Variation of Malaria Morbidity.
e malaria surveillance data showed that 4,565,506 malaria cases occurred in the Amhara region between July 2012 and June 2020, with a decline from July 2012 to June 2018 and an increasing trend since 2019. Monthly total malaria cases have fluctuated during months of the year, and there were higher variabilities at the beginning of the study periods between 2012 and 2015 (Figure 2(a)). e intensity of malaria infection varied among different age groups, and more than three-quarters (78%) of malaria patients were aged above 14 years. e remaining 1% and 21% of malaria patients were under 5 years and 5-14 years, respectively (Figure 2(b)). Plasmodium falciparum was a more prevalent malaria parasite in the Amhara region, accounting for 67% of the confirmed cases between 2012 and 2020 ( Figure 2(c)). e seasonal variation of yearly malaria cases is depicted in Figure 2(d), revealing that malaria cases were higher from September to November, following the primary rainy season between June and August, in 2012-2019. Monthly malaria cases were consistently higher in 2012 and 2013 than in other periods. e number of malaria cases reached the bottom level in 2018 and had a consistently lower number of patients throughout all months (Figure 2(d), red dot line). However, the number of malaria cases was higher in 2019 and 2020 compared to 2018, which has an increasing trend that might be due to various conditions limiting malaria prevention, control, and elimination endeavors of concerned bodies. e average annual malaria incidence rate was 24.8 per 1,000 persons between 2012 and 2020 in the region. e annual malaria incidence rate was declined between 2012 and 2018, but there was a higher annual incidence rate in 2019. e result also indicates that total malaria and confirmed cases had consistently similar annual incidence rates between 2012 and 2020. e average annual incidence rate of malaria in patients and malaria in pregnant women were 9.68 and 16.93 per 100,000 persons per year, respectively, and reached the bottom level in 2018/19 (Table 1). e annual malaria incidence rate was also varied on the types of Plasmodium species, and the result indicated that the average annual incidence rate of P. falciparum and P. vivax were 16.32 and 8.17 per 1,000 persons per year, respectively (Table 1). Besides, the result revealed that the annual malaria mortality rate was 2 per 1,000,000 persons per year. e annual incidence rates of total malaria and confirmed cases, P. falciparum cases, and malaria in pregnancy were lowest in the 2017/18 fiscal year. On the contrary, P. vivax cases and malaria in patients had the lowest incidence rate in 2018/19. However, the annual malaria mortality rate was stable throughout the study years (Table 1). e annual malaria cases of districts in the study region are presented in Figure 3. Maps in Figure 3 depicted the spatial and temporal variations of malaria cases among districts of the Amhara region between 2012 and 2020 with a similar legend. e map showed that districts in the western Amhara region had higher malaria cases than districts located in the eastern parts of the region. Between 2012 and 2015, there were significantly higher malaria cases in North Gondar, South Gondar, Awi, East Gojjam, and West Gojjam zones.
Districts located in the highland areas of North Shewa, South and North Wollo, and Wag-Himra zones had fewer malaria cases either imported or infected through travel history ( Figure 3). Further, the map showed that malaria incidence would have spatial and temporal autocorrelations due to the clustering of similar colors in the neighborhood districts throughout the study period. Hence, the spatiotemporal variation of malaria incidence would be considered to estimate and explore space-time trends and dynamics of malaria transmission in the region.

Parametric Spatiotemporal Trend.
e model comparison result is presented in Table 2, and the spatiotemporal negative binomial model has the lowest DIC and WAIC values and is used to estimate spatiotemporal linear trends and nonparametric dynamics of malaria incidence. e parametric spatiotemporal trends of malaria incidence results are presented in Table 3. e intercept and slope of time in a month were negative, indicating that the log scale of malaria incidence had decreasing linear trend between 2012 and 2020. e malaria incidence rate was reduced by a factor of 0.984 (95% credible interval: 0.983-0.986) for a unit increment in months in the study period. Moreover, the precision of spatial and space-time interaction effects was estimated, and their 95% credible intervals showed the presence of significant variations in spatial heterogeneity, spatial clustering, and differential trends across space over time. e space-time interaction effect is estimated using the differential trend representing the difference between district-specific trends and the average trend of the region and is depicted in Figure 4(a). e differential trends revealed that districts had a geographical variation on their districtspecific trends. e values within Figure 4(a) are posterior probabilities used to test the presence of significantly different district-specific trends compared to the mean trend of the region. e 45 districts had a steeper trend than the mean trend since they had 0.9 or more posterior probabilities (Figure 4(a)). Figure 4(a) also depicted that most districts in the old North Gondar Zone, now reorganized into three zones named Central, West, and North Gondar, had a steeper trend than the average regional trend. Even though the North Shewa zone had a smaller number of malaria cases than other zones, most districts had a steeper trend than the mean trend of the region. On the contrary, districts located in West Gojjam and South Wollo Zone had a less steep trend than the average trend of the region (Figure 4(a)).    Journal of Tropical Medicine e study region has included 23 districts, which are solely town administrations, and most of them had a less steep malaria incidence trend than the average regional level trend. On the contrary, Bahir Dar, Lalibela, Kobo, Woldia, Bati, Debre Birhan, and Shewa Robit were town administrations that had a steeper trend than the mean trend of the   region. Generally, districts shaded in yellow to light orange colors had a lower rate of decreasing trend than the average trend of the region (Figure 4(a)). e spatial variation of malaria incidence was significantly explained using spatial effects through its components, that is, spatial clustering (structured) and spatial heterogeneity (unstructured). e 99.96% geographic variation of malaria incidence was explained using the spatial clustering of districts. e Moran's I statistic of the structured spatial effects was 0.443 (P-value � 0.0017 < 0.01), revealing the presence of significant spatial clustering of districts based on their malaria incidence. e estimated district/spatial effects and its map are given in Figure 4(b), indicating that as one goes towards the west and spatial effects of malaria incidence get higher. Malaria incidence had higher estimated spatial effects in the northwestern and western Amhara region, and it decreased as one moved towards the central and eastern parts of the region. Districts, mainly located in the highland fringe areas, had relatively smaller spatial effects, that is, in North Shewa, South Wollo, and North Wollo zones (Figure 4(b)).

Nonparametric Dynamic Trend of Malaria Incidence.
e dynamic feature of the malaria incidence trend was estimated using negative binomial spatiotemporal dynamic model underlying separable space-time interaction effects, which has the smallest DIC value (Table 2). e nonparametric dynamic trends of malaria incidence were decomposed into structured and unstructured spatial and temporal effects, and their results are presented in Table 4. e result revealed that malaria incidence had significant spatial and temporal variation among districts of the study variations. e structured and unstructured temporal dynamic trends of malaria incidence are shown in Figure 5. e structured temporal trend has an overall decreasing nonlinear trend between 2012 and 2018 and increased between 2018 and 2020 with a noticeable seasonal variation in each year. However, the unstructured temporal effects had random oscillations. e series plot had many peaks and changes over months that indicate the seasonality of malaria morbidity and higher temporal effects occurred between September and November throughout the study year.

Discussion
e focus of this study was estimating linear trends and dynamic features of malaria morbidity in the Amhara region using monthly surveillance data. More than 4.6 million people were confronted with malaria problems in the region between June 2012 and July 2020. e average malaria incidence rate was 24.8 per 1,000 persons per year between Journal of Tropical Medicine 7 2012 and 2020, which is a bit greater than in 2004-2014 [7]. e malaria incidence rate was declined three to fourfold in 2018 compared to 2012 and increased starting from mid-2018 might be due to drug and insecticide resistance, social, demographic, cultural, and behavioral beliefs and practices, and unreformed health infrastructure [30]. e result suggests that malaria incidence of the region followed a similar pattern of changes as the national and global malaria transmission patterns [4,31]. Malaria is public health burden that affects individuals in all age intervals, and 78% of malaria patients were aged above 14 years. e most prevalent malaria parasite was P. falciparum, which accounted for 67% of confirmed cases in the region, a predominant malaria parasite in Ethiopia and the WHO African region [4,6,[31][32][33]. e number of malaria cases was higher between September and December following the main rainy season in the region, and there were also considerably higher malaria cases from May to July. Seasonal variability might be attributed to the suitability of environmental conditions for the reproduction of mosquito vectors [34] and influenced by crop cycle, crop weeding, and grass cover of lands that appeared on and after rainy seasons [35]. e suitability of environmental conditions determines the distribution of Plasmodium species in space and time, and the spatiotemporal distributions of malaria have been related to this [36][37][38]. e annual malaria incidence varied across districts of the study region, and space-time distribution of malaria declined between 2012 and 2018. However, compared to 2018, most districts in the Amhara region had an increased number of malaria patients from 2019 to 2020.
Districts in the northwestern and western parts of the region had higher annual malaria incidence than the eastern parts throughout the study period. Mainly, South Gondar, North Gondar, Awi, East Gojjam, and Wag-Himra zones had higher malaria incidence between 2012 and 2015, similar to a previous study finding in the west and near the border with Sudan and South Sudan [4]. e result is also supported by Alemu et al.'s [6] findings, which indicated that malaria transmission remained high in northwest Ethiopia between 2003 and 2012. On the contrary, districts located in the highlands of the North Shewa, North and South Wollo, and the Oromo special zone had lower malaria incidence between 2012 and 2020. e estimated linear spatiotemporal trend suggested that the Amhara region had a decreasing malaria trend between 2012 and 2020; that is, between July 2012 and June 2020, per a unit increment in months of the year, the rate of malaria incidence decreased by a factor of 0.984, which is consistent with findings from previous studies conducted in Ethiopia between 2004 and 2016, 2011 and 2016, and 2013 and 2018 [3,4,32]. e result contradicted Taye et al.'s [7] predictions of an increased malaria trend in the Amhara region in 2015-2020, which might be occurred due to differences in data aggregation time scale and spatial units. Among the 152 districts, 45 had a significantly steep trend than the average trend of the region. e majority of districts with a higher malaria incidence in the North Gondar zone had a higher declining rate than the mean trend and had spatial disparity that is in line with the findings of Yalew et al. [9]. Among 23 town districts, about 70% of them had a less steep trend than the average regional trend, supported by Doumbe-Belisse  [39]. However, Bahir Dar, Lal Yibela, Kobo, Woldia, Bati, Debre Birhan, and Shewa Robit were urban districts with a highly declining malaria incidence rate compared to the regional level trend. e nonparametric dynamics of malaria incidence decomposed into spatial and temporal effects, each with structured and unstructured heterogeneity [20]. e structured monthly temporal dynamics fluctuated in the study periods and had seasonal variations, and higher incidence occurred following the main rain seasons between July and December [9]. Overall, malaria incidence decreased between 2012 and 2018 and began rising since 2019 across all months of the year. Compared to districts located in the eastern Amhara, the structured spatial effect of districts in western Amhara was higher.
e results revealed that the North Gondar zone had a higher spatial effect due to peak malaria incidence. In contrast, North Shewa and South Wollo zones had lower spatial effects [6].
Spatiotemporal trend analysis is needed to provide insight and evidence supporting policy decision-making to prevent and control infectious diseases [40]. e parametric and nonparametric spatiotemporal trends are crucial for estimating area-level trends in considering spatial-time interaction effects. It is also used to detect areas that require emergency prevention and control interventions and can be reproduced anywhere with spatiotemporal areal data. e district-level trends would play a valuable role in evaluating the progress of malaria prevention and control performances of the districts' health offices. Furthermore, malaria prevention, control, and elimination could target districts with lower or higher malaria trends to achieve the malaria elimination target of the region. is study did not include climate, environmental, and other malaria transmission intervention impacts on the estimation of the nonparametric spatiotemporal dynamics of malaria, which is the limitation of the study.

Conclusions
e APHI malaria surveillance data have been used to explore and investigate the spatiotemporal linear and dynamic temporal trend of malaria incidence across districts of the Amhara region, Ethiopia. e Amhara region had a linearly declining malaria incidence trend between 2012 and 2020. e finding of this study revealed that districts in the western and southwestern parts of the Amhara region had a steeper trend than the average trend of the region. In contrast, districts in the eastern part of the region had a less steep trend than the average trend of the region. Further, twothirds of urban districts had a less steep trend incidence than the average regional trend. Generally, monthly malaria incidence had sinusoidal wave dynamics that varied across months, and where there was an overall decreasing trend between 2012 and 2018. However, the trend of malaria incidence was reversed and showed an increasing trend.
us, an intervention and controlling mechanism that considers malaria incidences and district-specific differential trends would be indispensable to mitigate malaria transmission in the region.

Abbreviations
WHO: World health organization EPHI: Ethiopian public health institute APHI: Amhara public health institute FMOH: Federal ministry of health DIC: Deviance information criterion WAIC: Watanabe-Akaike information criterion PHEM: Public health emergency management RDT: Rapid diagnostic test EPI: Epidemiology PF: Plasmodium falciparum PV: Plasmodium vivax.

Data Availability
e malaria data used to support the findings of the study are available from the corresponding author upon request.

Ethical Approval
Ethical clearance was approved by Bahir Dar University, College of Science Ethical Review Board, reference no. PGRCSVD/137/201.

Conflicts of Interest
e authors declare that they have no conflicts of interest.

Authors' Contributions
TZN analyzed and drafted the manuscript. TZ and EKM, supervisors for TZN, designed and supervised from the design to the write-up of the manuscript. All authors read and approved the final manuscript.