1. Epidemiology and Global Health
Download icon

Spatio-temporal associations between deforestation and malaria incidence in Lao PDR

  1. Francois Rerolle  Is a corresponding author
  2. Emily Dantzer
  3. Andrew A Lover
  4. John M Marshall
  5. Bouasy Hongvanthong
  6. Hugh JW Sturrock
  7. Adam Bennett
  1. Malaria Elimination Initiative, The Global Health Group, University of California, San Francisco, United States
  2. Department of Epidemiology and Biostatistics, University of California, San Francisco, United States
  3. Department of Biostatistics and Epidemiology, School of Public Health and Health Sciences, University of Massachusetts Amherst, United States
  4. Division of Epidemiology and Biostatistics, School of Public Health, University of California, Berkeley, United States
  5. Center for Malariology, Parasitology and Entomology, Ministry of Health, Lao People's Democratic Republic
Research Article
  • Cited 1
  • Views 971
  • Annotations
Cite this article as: eLife 2021;10:e56974 doi: 10.7554/eLife.56974

Abstract

As countries in the Greater Mekong Sub-region (GMS) increasingly focus their malaria control and elimination efforts on reducing forest-related transmission, greater understanding of the relationship between deforestation and malaria incidence will be essential for programs to assess and meet their 2030 elimination goals. Leveraging village-level health facility surveillance data and forest cover data in a spatio-temporal modeling framework, we found evidence that deforestation is associated with short-term increases, but long-term decreases confirmed malaria case incidence in Lao People’s Democratic Republic (Lao PDR). We identified strong associations with deforestation measured within 30 km of villages but not with deforestation in the near (10 km) and immediate (1 km) vicinity. Results appear driven by deforestation in densely forested areas and were more pronounced for infections with Plasmodium falciparum (P. falciparum) than for Plasmodium vivax (P. vivax). These findings highlight the influence of forest activities on malaria transmission in the GMS.

eLife digest

Biting mosquitos spread the malaria parasite to humans. Along the Mekong River in Southeast Asia, spending time in the surrounding forest increases a person's risk of malaria. This has led to a debate about whether deforestation in this area, which is called the Greater Mekong Sub-region (GMS), will increase or decrease malaria transmission. The answer to the debate is not clear because some malaria-transmitting mosquitos thrive in heavily forested areas, in particular in the GMS, while others prefer less forested areas.

Scientists studying malaria in the Amazon in South America suspect that malaria transmission increases shortly after deforestation but decreases six to eight years later. Some studies have tested this ‘frontier malaria’ theory but the results have been conflicting. Fewer studies have tested this theory in Southeast Asia. But deforestation has been blamed for recent malaria outbreaks in the GMS.

Using data on malaria testing and forest cover in the GMS, Rerolle et al. show that deforestation around villages increases malaria transmission in the first two years and decreases malaria rates later. This trend was driven mostly by a type of malaria called Plasmodium falciparum and was less strong for Plasmodium vivax. The location of deforested areas also mattered. Deforestation within one to 10 kilometer of villages did not affect malaria rates. Deforestation further away in about a 30 kilometer radius did affect malaria transmission. Rerolle et al. suggest this may be because villagers have to spend longer times trekking through forests to hunt or harvest wood when the wider area is deforested.

Currently, National Malaria Control Programs in the GMS focus their efforts on reducing forest-related transmission. This study strengthens the evidence supporting this approach. The results also suggest that different malaria elimination strategies may be necessary for different types of malaria parasite. Using this new information could help malaria control programs better target resources or educate villagers on how to protect themselves. The innovative methods used by Rerolle et al. reveal a more complex role of deforestation in malaria transmission and may inspire other scientists to think more carefully about environmental drivers of malaria.

Introduction

Engaging in forest activities, such as logging, hunting, or spending the night in the forest, is considered a primary risk factor for malaria infection in the Greater Mekong Sub-region (GMS) (Chaveepojnkamjorn and Pichainarong, 2004; Das et al., 2004; Lansang et al., 1997; Erhart et al., 2005; Trung et al., 2004) with recent outbreaks attributed to deforestation activities (Lyttleton, 2016). This is most likely the result of increased human exposure to the forest dwelling Anopheles dirus and Anopheles minimus, the most efficient and widespread malaria vectors in the GMS (Obsomer et al., 2007; Obsomer et al., 2013). However, deforestation may also alter this ‘forest malaria’ (Sharma and Kondrashin, 1991) ecosystem and eventually reduce malaria incidence, as is generally accepted to be the case in Southeast Asia (Guerra et al., 2006). Several previous studies have assessed the relationship between deforestation and malaria, but the majority focused outside of the GMS, most notably in the Amazonian forest (Wayant et al., 2010; Olson et al., 2010; Hahn et al., 2014a; Valle and Clark, 2013; Terrazas et al., 2015) where the evidence has been met with conflicting interpretations (Tucker Lima et al., 2017). As national malaria programs across the GMS target forest-going populations for prevention and treatment efforts (Guyant et al., 2015; WHO, 2016), improved understanding of the relationship between deforestation and malaria is critical for programs to assess and meet national 2030 malaria elimination goals (World Health Organization, 2016; WHO, 2018).

In the Amazon, the ‘frontier malaria’ hypothesis (Sawyer, 1988) posits that malaria temporarily increases following deforestation efforts to open a human settlement area in the forest. Subsequently, after approximately 6–8 years, settlements become more urbanized and isolated from the surrounding forest, and less suitable for malaria vectors, resulting in reduced malaria transmission (de Castro et al., 2006). Recent work has challenged this hypothesis, however, and found that some older settlements were also likely to have high malaria incidence (Ilacqua et al., 2018), highlighting the importance of assessing the relationship between deforestation and malaria at different spatio-temporal scales (Singer and de Castro, 2006).

A recent review of the literature on deforestation and malaria in the Amazon (Tucker Lima et al., 2017) recommended the integration of multiple socio-economic, demographic and ecological mechanisms to disentangle the relationship between deforestation and malaria. The complexity of land-use changes driving deforestation such as urbanization, agriculture, irrigation or resource mining can alter the environment in different ways. For example, deforestation in the Amazon has been shown to increase mosquitoes’ larval habitat (Vittor et al., 2009) through the creation of areas with abundant sunlight and pooling water, resulting in increased human biting activity (Vittor et al., 2006). Alternatively, immigration and rapid population movements, stirring human-vector interactions, are other mechanisms affecting malaria transmission in frontier areas (Moreno et al., 2007). A modeling study (Baeza et al., 2017) showed that the temporal pattern of increased incidence followed by a decrease can vary depending on ecological and socio-economic parameters in frontiers areas.

The importance of addressing complex confounding structures influencing the relationship between deforestation and malaria was also highlighted by Bauhoff and Busch, 2020. Variables such as temperature (Beck-Johnson et al., 2013; Mordecai et al., 2013), precipitation (Parham and Michael, 2010; Parham et al., 2012), or seasonality Hay et al., 1998 are known environmental predictors of malaria, although the spatio-temporal scale of those effects often varies across different areas (Teklehaimanot et al., 2004). Furthermore, remote areas may experience higher malaria rates because of poor access to public health services, but also have denser forest cover or lower deforestation rates (Busch and Ferretti-Gallon, 2017). Finally, forest-going populations in the GMS are also at higher risk for malaria (Parker et al., 2017) due to poor adherence to protective measures against mosquitoes such as insecticide-treated bed nets (ITNs) or long-lasting insecticidal hammocks (LLIHs) (Peeters Grietens et al., 2010; WHO, 2016) and inadequate access to treatment (Guyant et al., 2015).

Bauhoff and Busch, 2020 identified only 10 empirical studies that assessed the relationship between deforestation and malaria with appropriate adjustments for confounding. Of these, seven reported a positive association (Austin et al., 2017; Wayant et al., 2010; Olson et al., 2010; Terrazas et al., 2015; Pattanayak et al., 2010; Garg, 2015; Fornace et al., 2016), two did not find any associations (Bauhoff and Busch, 2020; Hahn et al., 2014a), and one disputed study found a negative association (Valle and Clark, 2013; Hahn et al., 2014b; Valle, 2014). Most recently, a study found deforestation to increase malaria risk and malaria to decrease deforestation activities in the Amazon, using an instrumental variable analysis to disentangle any reverse causality loop (MacDonald and Mordecai, 2019). However, only half of the above-mentioned studies used high-resolution forest data, with most studies using spatially aggregated data and exploring only a limited range of spatial and temporal scales. Only three of these studies were conducted in Southeast Asia (Pattanayak et al., 2010; Garg, 2015; Fornace et al., 2016), and none in the GMS. Importantly, all three found that malaria increases after deforestation, but all had limitations. The two studies in Indonesia Pattanayak et al., 2010; Garg, 2015 used coarsely aggregated forest data and potentially biased self-reported malaria data. The third study, in Malaysia (Fornace et al., 2016), focused on a specific and geographically confined malaria parasite, Plasmodium knowlesi, whose primary host is the long-tailed macaque and whose presence in the GMS, where P. falciparum and P. vivax dominate, is limited.

In this analysis, we examined the relationship between deforestation and malaria incidence by combining high-resolution forest coverage data (Hansen et al., 2013) and monthly malaria incidence data from 2013 to 2016 from two separate regions in the GMS: northern Lao People’s Democratic Republic (PDR) with very low malaria transmission and southern Lao PDR where P. falciparum and P. vivax are seasonal. By conducting the analysis at the village level, we were able to explore a wide range of spatial scales (1, 10, and 30 km around villages) that might be relevant in characterizing the relationship between deforestation and malaria. In addition, we leveraged the longitudinal nature of both the incidence data collected and the forest data produced from annual remote sensing imagery (Hansen et al., 2013) to explore the most relevant temporal scales. Finally, we considered alternative definitions of deforestation, restricted to areas with at least certain levels of forest cover, to investigate the type of deforestation driving the relationship with malaria.

To date, no prior studies have quantified the relationship between deforestation and malaria incidence in the GMS. Understanding this relationship is especially important in the GMS, where forest-going activities are a main source of income generation (Lao PDR, 2016) and malaria clusters in forest-going populations (World Health Organization, 2016; WHO, 2016). To assess the relationship between deforestation and malaria incidence, we modeled the monthly village-level malaria incidence in two regions of Lao PDR using health facility surveillance data and evaluated the most relevant spatio-temporal scale.

Results

Forest and environmental data

Figure 1 shows the average tree crown cover within 10 km for the year 2016 and the percent area within 10 km that experienced forest loss between 2011 and 2016 in two regions of southern and northern Lao PDR. Overall, the forest cover was denser in the north than in the south and deforestation over this period was higher in the north than in the south. Appendix 1—figure 6 and Appendix 1—figure 7 show the distribution of forest and deforestation variables as the temporal scales and spatial scales around study villages were varied. For example, the cumulative percent area within 30 km of a village that experienced forest loss between 2011 and 2016 ranged from 0 to 10% in the north, whereas it rarely exceeded 2.5% in the south. Deforestation rate in 2015 within 10 km of a village was of about 1% in the south and 2.5% in the north. The average tree crown cover increased with increasing buffer radius around villages (1, 10, and 30 km). However, the relationship with the percent area that experienced forest loss was less clear because both the area that experienced forest loss (numerator) and the area around villages (denominator) increased. Appendix 1—figure 15 and Appendix 1—figure 16 show the raw time series of forest cover and percent area with forest loss.

Average tree crown cover (%) in 2016 (left) and percent area that experienced forest loss between 2011 and 2016 (right) within a 10 km radius in northern (top) and southern (bottom) Lao PDR.

See Materials and methods for details on forest and deforestation metrics. Upper right indent maps northern and southern (Champasak province) Lao PDR regions.

The 284 villages in the north were overall less populated (mean 2015 population size: 498, IQR: [241; 548]) than the 207 villages in the south (2015 mean: 1095, IQR: [584; 1384]), but with some highly populated outliers. As expected, altitude differed substantially between villages of the mountainous northern region (mean: 557 m, IQR: [378; 679]) and the lowlands of the south (mean: 120 m, IQR: [98; 130]). Although both regions exhibited similar seasonal trends in precipitation and temperature, with a rainy season spanning from April to October, villages in the south experienced higher monthly precipitation and temperature than in the north over the study period (Appendix 1—figure 8).

Treatment-seeking data

For villages with an estimated travel time of 0 hr to the closest health facility (same 300 m2 pixel), the predicted probability of seeking treatment for fever was 0.87 (95% CI: [0.79; 0.92]) in the north and 0.78 (95% CI: [0.63; 0.89]) in the south. A 1 hr increase in travel time to the closest health facility was associated with a similar 0.79 (95% CI: [0.55; 1.13]) reduction in the odds of seeking treatment in the north and 0.76 (95% CI: [0.43; 1.34]) in the south, almost reaching statistical significance when pooling data from both regions: 0.77 (95% CI: [0.56; 1.04]). See detailed results in Appendix 1 - S2.

Malaria case data

Malaria infections

A total of 63,040 patient records were abstracted from the malaria registries of all public health facilities in four southern districts between October 2013 and October 2016 and 1754 from all health facilities in four northern districts between January 2013 and December 2016.

In the south, 91.2% of the patients in the registries were tested for malaria, of which 78.1% were tested by RDT and 26.2% by microscopy. Overall test positivity was 33.2% for any infection, 16.4% for P. falciparum and 18.2% for P. vivax. Monthly incidence peaked to about six cases per 1000 people in the 2014 rainy season, eventually decreasing to below one case per 1000 in 2016. Incidence and test positivity were similar between P. falciparum and P. vivax in the south (Figure 2).

Malaria incidence (per 1000) and test positivity (%) over time.

Upper left boxed indent zooms in malaria incidence in the North to better show the temporal variation (see y axis for scale).

In the north, 92.1% of the patients in the registries were tested for malaria, of which 96.3% were tested by RDT and 9.6% by microscopy. Overall test positivity was 23.8% for any infection, 2.8% for P. falciparum and 22.5% for P. vivax. Monthly malaria incidence in the north was very low, never exceeding 0.3 per 1000 people. Most infections in the north were P. vivax cases with only a few seasonal P. falciparum cases (Figure 2).

Appendix 1—figure 9 shows the number of patients and cases recorded per month in health facility malaria registries as well as how the smoothed test positivity rates varied over time.

Socio-demographics

Age, gender, and occupation of patients seeking treatment at health facilities were also recorded in the malaria registries. On average, patients in the south were older than patients in the north with mean age of 28 years and 23 years, respectively. In the north, about half of the patients were male (53.1%), while most patients in the south were male (71.1%). Finally, the vast majority (68.2%) of patients in the south were farmers, whereas only 8% of patients in the north were farmers. Most patients in the north reported being unemployed (41.7%) or a student (31.2%) (Appendix 1—figure 10).

Geo-referencing

Overall, 88.1% of malaria records were matched to one of the 491 villages in study districts. The remaining (11.7% in the south and 17.3% in the north) were removed from the analysis because of ambiguous village names, local nicknames for small villages and dissolving and grouping of villages over time. Test positivity in the south was similar in matched (33.1%) and unmatched (34.2%) records but higher in matched (26.5%) than unmatched (10.5%) records in the north. No substantial difference was found in the distribution of socio-demographic variables available in malaria registries between matched and unmatched records (Appendix 1—figure 11). Fewer than 0.3% of matched malaria records were missing dates and also removed from the analysis.

Spatio-temporal analysis

Deforestation

Table 1 and Figure 3 show the adjusted incidence rate ratio (IRR) associated with deforestation, measured by a 0.1% increase in the percent area that experienced forest loss, in the previous 1–5 years within 1, 10, and 30 km of villages. Models controlled for various environmental factors and accounted for the probability of seeking treatment and the spatio-temporal structure of the data. 

Table 1
IRR between malaria incidence and a 0.1% increase in the area that experienced deforestation within 1, 10, or 30 km (left-right) of a village in the previous 1– 5 years (top-down) in northern and southern Lao PDR.

Adjusted for the probability of seeking treatment, the spatio-temporal structure of the data, the environmental covariates selected in the model and forest cover within 30 km in the year before the deforestation temporal scale considered as well as for malaria incidence in the previous 1 and 2 years. See Materials and methods for details.

SouthNorth
Time lagBuffer radiusBuffer radius
1 km10 km30 km1 km10 km30 km
Previous11.011.1611.031.01
1 year[0.99; 1.01][0.99; 1.04][1.10; 1.22][1; 1.01][0.99; 1.06][0.94; 1.08]
Previous111.0811.010.99
2 years[0.99; 1.01][0.98; 1.01][1.04; 1.13][1; 1.01][0.99; 1.04][0.95; 1.03]
Previous0.990.980.9311.010.96
3 years[0.99; 1][0.97; 1][0.90; 0.97][1; 1.01][0.99; 1.02][0.94; 0.99]
Previous0.990.980.94110.97
4 years[0.99; 1][0.97; 0.99][0.92; 0.97][1; 1.01][0.99; 1.02][0.94; 0.99]
Previous10.970.9411.010.96
5 years[0.99; 1][0.96; 0.99][0.91; 0.97][1; 1.01][0.99; 1.02][0.93; 0.98]
Associations between malaria incidence and a 0.1% increase in the area that experienced deforestation within 1, 10, or 30 km (left-right) of a village in the previous 1–5 years (top-down) in Lao PDR.

Adjusted for the probability of seeking treatment, the spatio-temporal structure of the data, the environmental covariates selected in the model and forest cover within 30 km in the year before the deforestation temporal scale considered as well as for malaria incidence in the previous 1 and 2 years. See Materials and methods for details.

Deforestation within 1 or 10 km of a village was not associated with malaria incidence rate in either the south or the north, regardless of the temporal lag. However, in the south, deforestation within 30 km of a village in the previous 1 and 2 years was associated with higher malaria incidence rates (e.g. 1 year lag, IRR = 1.16, 95% CI: [1.10; 1.22]). In the north, where incidence was much lower, the results were not as clear, but a similar trend was observed with wide confidence intervals compatible with a short term increased risk. On the other hand, deforestation within 30 km of a village in the previous 3, 4, or 5 years was associated with approximately a 5% lower malaria incidence rate both in the south (e.g. 5 year lag, IRR = 0.94, 95% CI: [0.91; 0.97]) and in the north (e.g. 5 year lag, IRR = 0.96, 95% CI: [0.93; 0.98]).

These results suggest deforestation around villages, but not in the near vicinity (1 or 10 km), is associated with higher risk of malaria in the first 2 years and lower risk of malaria beyond. There was stronger evidence of associations with deforestation in the south than in the north.

The IRR effect estimates in Table 1 assume a linear relationship between deforestation and malaria. Appendix 1—figure 5 shows a few of these relationships when such linearity isn’t assumed in the models. The functional forms reveal that they can be reasonably well summarized linearly, especially in the south. In the north, the functional forms highlight potential non-linearities for long-term temporal lags but come with wide confidence intervals at extreme levels of deforestation.

P. falciparum and P. vivax

In addition to different overall levels of transmission in the north and south, the relative species composition also differs by region. In the north, P. vivax is more prevalent with only a few sporadic and seasonal P. falciparum infections, whereas P. falciparum and P. vivax are co-endemic in the south (Figure 2). We used the co-endemicity and the larger amount of malaria case data collected in the south to assess the relationship between deforestation and malaria for both species separately.

Table 2 and Figure 4 show that the pattern of adjusted spatio-temporal associations identified in Table 1 is primarily driven by P. falciparum, with no associations for deforestation in the near vicinity of villages (1 or 10 km) but a short-term increase (e.g. 1 year lag, IRR = 1.27, 95% CI: [1.18; 1.36]) and long-term decrease (e.g. 5-year lag, IRR = 0.83, 95% CI: [0.80; 0.87]) in P. falciparum malaria incidence for deforestation within 30 km of villages.

Associations between malaria incidence and a 0.1% increase in the area that experienced deforestation within 1, 10, or 30 km (left-right) of a village in the previous 1–5 years (top-down) in southern Lao PDR, differentiated by malaria species.

Adjusted for the probability of seeking treatment, the spatio-temporal structure of the data, the environmental covariates selected in the model and forest cover within 30 km in the year before the deforestation temporal scale considered as well as for malaria incidence in the previous 1 and 2 years. See Materials and methods for details.

Table 2
IRR between malaria incidence and a 0.1% increase in the area that experienced deforestation within 1, 10, or 30 km (left-right) of a village in the previous 1–5 years (top-down) in southern Lao PDR, differentiated by malaria species.

Adjusted for the probability of seeking treatment, the spatio-temporal structure of the data, the environmental covariates selected in the model and forest cover within 30 km in the year before the deforestation temporal scale considered as well as for malaria incidence in the previous 1 and 2 years. See Materials and methods for details.

P. falciparumP. vivax
Time lagBuffer radiusBuffer radius
1 km10 km30 km1 km10 km30 km
Previous11.041.27111.07
1 year[0.99; 1.02][1.01; 1.07][1.18; 1.36][0.99; 1.01][0.97; 1.02][1.01; 1.13]
Previous11.011.15111.06
2 years[0.99; 1.01][0.99; 1.03][1.08; 1.22][0.99; 1.01][0.98; 1.01][1.01; 1.11]
Previous0.990.990.8510.991.02
3 years[0.98; 1][0.97; 1.01][0.80; 0.90][0.99; 1.01][0.98; 1.01][0.97; 1.06]
Previous0.990.980.8510.991.01
4 years[0.98; 1][0.96; 0.99][0.81; 0.88][0.99; 1][0.98; 1.01][0.98; 1.04]
Previous10.970.8310.991.01
5 years[0.99; 1][0.95; 0.98][0.80; 0.87][1; 1.01][0.98; 1][0.98; 1.04]

On the other hand, all the associations were attenuated for P. vivax infections. In the previous 2 years and within 30 km of villages, deforestation is still associated with a higher incidence of P. vivax (e.g. 1 year lag, IRR = 1.07, 95% CI: [1.01; 1.13]) but less so than for P. falciparum. However, regardless of the temporal lag or spatial scale, deforestation was no longer associated with lower P. vivax malaria risks.

Appendix 1—figure 17 plots the species-specific relationships when not assuming linearity in the models.

Alternative definitions of deforestation and interaction with forest cover

In previous models, our definition of deforestation did not distinguish between forest losses in densely forested areas and less forested areas. To explore potential interactions between deforestation and baseline forest cover, Table 3 and Figure 5 show how the adjusted IRR estimates vary as we consider deforestation in more densely forested pixels only (tree crown cover over 68% and 87% – see Materials and methods for rationale on thresholds). We conducted this secondary analysis only for the non-null relationships previously identified, that is, when considering a 30 km buffer radius around villages.

Associations between malaria incidence and a 0.1% increase in the area that experienced deforestation within 30 km of a village in the previous 1–5 years (top-down) and within areas with tree crown cover density above 0%, 68%, and 87% (left-right) in Lao PDR.

Adjusted for the probability of seeking treatment, the spatio-temporal structure of the data, the environmental covariates selected in the model and forest cover within 30 km in the year before the deforestation temporal scale considered as well as for malaria incidence in the previous 1 and 2 years. See Materials and methods for details.

Table 3
IRR between malaria incidence and a 0.1% increase in the area that experienced deforestation within 30 km of a village in the previous 1–5 years (top-down) and within areas with tree crown cover density above 0%, 68%, and 87% (left-right) in Lao PDR.

Adjusted for the probability of seeking treatment, the spatio-temporal structure of the data, the environmental covariates selected in the model and forest cover within 30 km in the year before the deforestation temporal scale considered as well as for malaria incidence in the previous 1 and 2 years. See Materials and methods for details.

SouthNorth
Deforestation within areasDeforestation within areas
with tree crown cover density abovewith tree crown cover density above
Time lag0%68%87%0%68%87%
Previous1.161.321.281.011.041.34
1 year[1.10; 1.22][1.14; 1.53][1; 1.64][0.94; 1.08][0.96; 1.14][0.99; 1.81]
Previous1.081.181.350.990.980.94
2 years[1.04; 1.13][1.08; 1.28][1.15; 1.59][0.95; 1.09][0.93; 1.03][0.80; 1.11]
Previous0.930.930.890.960.960.86
3 years[0.90; 0.97][0.89; 0.97][0.81; 0.99][0.94; 0.99][0.92; 0.99][0.76; 0.98]
Previous0.940.940.870.970.960.87
4 years[0.92; 0.97][0.91; 0.97][0.80; 0.94][0.94; 0.99][0.93; 0.99][0.78; 0.96]
Previous0.940.930.830.960.950.83
5 years[0.91; 0.97][0.90; 0.96][0.76; 0.90][0.93; 0.98][0.92; 0.98][0.75; 0.92]

The associations with deforestation became more pronounced as we restricted forest losses to more forested areas: the adjusted IRR for deforestation in the previous 1 year, within 30 km of southern villages, increased from 1.16 (95% CI: [1.10; 1.22]) to 1.28 (95% CI: [1; 1.64]) when considering deforestation in areas with more than 0% and 87% tree crown cover respectively. On the other hand, the adjusted IRR for deforestation in the previous 5 years, within 30 km of southern villages, decreased from 0.94 (95% CI: [0.91; 0.97]) to 0.83 (95% CI: [0.76; 0.90]) when considering deforestation in areas with more than 0% and 87% tree crown cover, respectively. A similar trend was observed in the north, although statistical significance was not reached as frequently as in the south.

These evidence strengthen our previous results and suggest that deforestation in deep and dense forests is more closely associated with malaria incidence in villages than deforestation in less forested areas.

Discussion

Based on a large dataset of health facility surveillance records in two regions of Lao PDR, we found evidence that deforestation around villages is associated with higher malaria incidence over the short-term but lower incidence over the long-term (e.g, in the south, within 30 km of villages: IRR = 1.16 [1.10; 1.22] for deforestation in the previous year and IRR = 0.93 [0.90; 0.97] for deforestation in the previous 3 years). Our evaluation of alternative spatial scales identified strong associations for deforestation within a 30 km radius around villages but not for deforestation in the near (10 km) and immediate (1 km) vicinity. Our results incorporated correction for the probability of seeking treatment, modeled as a function of distance to the closest health facility, as well as adjustment for several environmental covariates. Results appear driven by deforestation in densely forested areas and the patterns exhibited are clearer for infections with P. falciparum than for P. vivax.

The wide availability and longitudinal nature of malaria surveillance records collected routinely by the national program enabled exploration of the relationship between deforestation and malaria incidence over multiple spatio-temporal scales and across different levels of forest density. The spatio-temporal variability highlighted here provides insights into the causal mechanisms driving local-scale malaria incidence in the GMS. This approach not only quantified the deforestation-malaria incidence association in the GMS, but also strengthened the evidence for the key influence of forest-going populations on malaria transmission in the GMS.

This study’s results echo the frontier malaria hypothesis from the Amazon region, which posits an increase in malaria incidence in the first few years following deforestation and a decrease over the long term. However, we found an earlier inflexion point, 1–3 years after deforestation compared to 6–8 years in the Amazon (de Castro et al., 2006), most likely because of very different underlying human processes. Indeed, the frontier malaria hypothesis considers non-indigenous human settlements sprouting deeper and deeper in the forest whereas forest-going populations in the GMS are primarily members of established forest-fringe communities who regularly tour the forest overnight to hunt and collect wood (Dysoley et al., 2008). Industrial and agricultural projects or lucrative forest-based activities also attract mobile and migrant populations (MMPs) Guyant et al., 2015 in remote forested areas of the GMS but not on the same scale as the politically and economically driven unique colonization of the Amazon (de Castro et al., 2006).

Our results are also consistent with the three previous multivariable empirical studies (Garg, 2015; Pattanayak et al., 2010; Fornace et al., 2016) that assessed the effect of deforestation on malaria in Southeast Asia. Our study builds on these findings by using higher resolution forest data and exploring additional spatio-temporal scales. Using biennial village census data from Indonesia between 2003 and 2008 and district-aggregated remote sensing forest data, Garg, 2015 reported a 2–10.4% increase in the probability of a malaria outbreak in each village of districts that lost 1000 hectares of their forest cover in the same year. Using data from a 1996 cross-sectional household survey conducted in a quasi-experimental setting around a protected area in Indonesia, Pattanayak et al., 2010 found a positive association between disturbed forest (vs undisturbed) and malaria in children under 5, again using no temporal lag. Our analysis plan was largely inspired by Fornace et al., 2016, which used similar high-resolution forest data (Hansen et al., 2013) and 2008–2012 incidence data from Sabah, Malaysia. They reported a 2.22 (95% CI: [1.53; 2.93]) increase in the P. knowlesi incidence rate for villages where more than 14% (<8%, being the reference) of the surrounding area (within 2 km) experienced forest loss in the previous 5 years. On the other hand, our analysis explored wider spatial scales, bypassed any coarse categorization of forest and deforestation variables, corrected incidence for treatment-seeking probability, and most importantly focused on P. falciparum and P. vivax, the dominant malaria parasites in the GMS.

Engaging in forest activities, such as logging, hunting or spending the night in the forest, has been reported as a major risk factor by many studies in the region (Chaveepojnkamjorn and Pichainarong, 2004; Das et al., 2004; Lansang et al., 1997; Erhart et al., 2005; Trung et al., 2004). As countries of the GMS work toward malaria elimination, the literature stresses the key role of forest-going populations (Guyant et al., 2015; Nofal et al., 2019; Bannister-Tyrrell et al., 2019; Wen et al., 2016; Smith and Whittaker, 2014), although research programs highlight the challenges of accessing them (Bennett et al., 2021 ; Lover et al., 2019) as well as their diversity (Nofal et al., 2019; Bannister-Tyrrell et al., 2019). To our knowledge, no previous study has leveraged geo-spatial statistical analyses to characterize the importance of forest-going populations in the GMS. Our results suggest that deforestation in dense forests (Table 3) around villages, particularly areas further from the village (Table 1), is a driver of malaria in Lao PDR. We argue that this is indicative of the existence of a key high-risk group linking the deforestation patterns identified to malaria in the villages, namely a forest-going population. Deforestation captured by remote sensing in this setting likely reflects locations and times of heightened activity in the forest areas near villages, and therefore greater human-vector contact. We suspect longer and deeper trips into the forest result in increased exposure to mosquitoes, putting forest-goers at higher risk.

We conducted this study in northern and southern Lao PDR, where the malaria species composition differs, and assessed species-specific relationships in the south where P. falciparum and P. vivax are co-endemic. Our results highlight the challenges ahead of national programs with P. vivax elimination after successful P. falciparum elimination, as increasingly mentioned in the literature (Cotter et al., 2013; Kaehler et al., 2019). This study identified a clear pattern of spatio-temporal associations between P. falciparum and deforestation, but these were not apparent for P. vivax (Table 2). The increase in P. vivax incidence in the first 2 years following deforestation was identified as well but the associations were smaller than for P. falciparum. Importantly, deforestation was never associated with lower risks of P. vivax. A recent study in the Amazon MacDonald and Mordecai, 2019 reported a similar attenuation of the effects of deforestation on P. vivax compared to P. falciparum, most likely because of P. vivax parasites’ ability to relapse months or even years after infection, which decouples the association between transmission and incidence data. These species-specific differences may also explain why the pattern of spatio-temporal associations between malaria and deforestation were markedly clearer in the south than in the north where P. vivax dominates.

Our results did have some inherent limitations based upon routine health facility surveillance data. First, reliability of such records varies across and within countries of the GMS and may depend on malaria incidence level. This could lead to unmeasured residual confounding, further exacerbated by the lack of available data on malaria control activities in the region. Another challenge with these data is obtaining an accurate denominator for incidence, as not everyone attends a public health center when febrile. We addressed this issue by modeling the probability of seeking treatment as a function of travel time to the closest health facility using data from two cross-sectional surveys. Last, the village-level geo-referencing of malaria registries ignores the possibility that patients may become infected elsewhere. Unfortunately, these surveillance records did not include information about patients’ forest-going trips. Research to track and analyse micro-scale movements of forest-goers is needed to understand how they interact with the forest and where are the foci of infection.

The forest data we used has also been criticized, in particular for not distinguishing tropical forests from agroforestry (Tropek et al., 2014; Hansen et al., 2014) or man-made from natural causes of deforestation. The lack of temporal resolution for the forest gain variable (2000–2017 aggregate) as well as the assumption that forest loss happens all in 1 year are additional limitations of these data. Finally, our relative measure of deforestation, key to consistently compare the effects across different spatial scales, also implies that a 0.1% of the area that experienced forest loss within 30 km of a village is a much larger area (∼280 hectares) than within 1 km (∼0.3 hectare) and should be interpreted cautiously.

In conclusion, this study assessed the relationship between deforestation and malaria in Lao PDR. Our approach leveraged surveillance records collected by the national malaria program and high-resolution forest data and rigorously explored the spatio-temporal pattern of associations. As countries of the GMS work toward malaria elimination, our results highlight the challenges to transition from P. falciparum to P. vivax elimination, confirm and characterize the importance of high-risk populations engaging in forest activities and suggest malaria programs may benefit from monitoring areas of on-going deforestation using remotely sensed data.

Materials and methods

Study site and population

Request a detailed protocol

Lao PDR has seen a 92% reduction in cases between 2000 (280,000) and 2010 (23,000) (Lao PDR, 2016). Much of this progress has been attributed to heightened funding and better testing and treatments (Okayas, 2018).

This study was conducted in eight districts (Figure 6) to leverage the ecological and epidemiological diversity of Lao PDR. Four districts (Moonlapamok, Pathoomphone, Sanasomboon, and Sukhuma) are situated in the southern province of Champasak where both P. falciparum and P. vivax are endemic. The four other districts (Et, Paktha, Nambak, and Khua) each come from one of four northern provinces (Bokeo, Huaphanh, Phongsaly, Luang-Prabhang) where P. vivax is endemic but P. falciparum has reached historical lows (Lao national malaria database (dhis2), 2018).

Map of study’s districts.

The four districts in the north were chosen in consultation with district and provincial level malaria staff to represent the epidemiology of malaria in the region. They were selected as part of a cross-sectional survey designed to assess the prevalence and risk factors for malaria in northern Lao PDR (Lover et al., 2018). This region is very mountainous and characterized by a diverse climate, low-population density and limited road access (UNFPA, 2016). Land clearing using fires for agriculture is customary.

The four districts in the south were selected within a larger cluster randomized controlled trial (RCT) study designed to assess the effectiveness of high-risk group targeted active case detection in southern Lao PDR (Lover et al., 2019), where more than 95% of the country malaria burden is concentrated (Lao national malaria database (dhis2), 2018). This region is characterized by a moderately hilly and forested terrain and a workforce primarily engaged in forest-based and agricultural activities (Bennett et al., 2021).

When designing the study, in collaboration with the national control program, we purposefully excluded regions where we knew large programmatic activities where being implemented.

Malaria data

Malaria case data

Request a detailed protocol

We conducted a retrospective review of malaria registries recorded at all health centers in the study districts between January 2013 and December 2016 in the north and between October 2013 and October 2016 in the south. The registries included information on every patient that was tested (RDT and/or microscopy) for malaria at the health center. Date, species-specific test results, demographic variables (age, gender, and occupation) and the village of residence of the patient were recorded in the registries. With help from local Lao experts, village names were matched to a geo-registry of all villages in Lao PDR compiled from the 2005 and 2015 national census Lao national census, 2020 and provided by the Center for Malariology, Parasitology and Entomology (CMPE). The geo-registry contains GPS coordinates and population of Lao PDR’s villages. Unmatched records and records with missing date were removed from the analysis. Finally, these data were aggregated to extract the monthly village-level malaria incidence.

Treatment-seeking data

Request a detailed protocol

One issue with using passive surveillance data is that not everyone will seek treatment at a public health facility for a febrile illness, which can lead to an underestimate of the true incidence, if not accounted for. To correct for that, we modeled the probability that an individual in a given village of the study’s district would seek treatment at a public health facility when febrile. We assumed that such probability is essentially driven by the travel time to the closest health facility. See Appendix 1 - S1 for methods used to calculate travel times to closest health facilities.

To model the probability of seeking treatment, we used data from two cross-sectional household surveys conducted in the eight districts where registries were collected. In the north, 1480 households across 100 villages were surveyed in September-October 2016 (Lover et al., 2018). In the south, 1230 households across 56 villages were surveyed in the baseline assessment of the RCT (Lover et al., 2019) in December 2017. In particular, survey respondents were asked whether or not they would seek treatment at the closest health facility for a febrile illness and GPS coordinates of their household were recorded.

We then used the cross-sectional surveys to model the probability of seeking treatment (at a public health facility, implicit from now on), θ, as a function of travel time to the closest health facility, τ (Equation 1). To account for the correlation structure induced by the stratified sampling approach used in the surveys, we modeled the number of successes (febrile patients seeking treatment), Sh,v, at the household level and included a random intercept for village in the logistic regression.

(1) Sh,vBin(θh,v,Nh,v)logit(θh,v)=α0+α1τh,v+αv

where Nh,v is the number of febrile individuals in household h of village v and αv𝒩(0,σα).

We fit the models separately in the north and in the south and used the region-specific model to predict the probability of seeking treatment at all villages of the study districts based on their distance to the closest health facility. The population who seek treatment was then calculated by multiplying the village population by the probability of seeking treatment. See Appendix 1 - S2 for travel times and treatment-seeking probabilities results.

Forest data

For every 30 m pixel in Lao PDR, tree crown cover density for the year 2000 and year of forest loss between 2000 and 2017, were obtained from Hansen et al., 2013. These layers were produced using decision tree classifiers on Landsat remote sensing imagery (Hansen et al., 2013). Trees are defined as ‘all vegetation taller than 5 m in height’ (Hansen et al., 2013) and forest loss as ‘the removal or mortality of all tree cover in a Landsat pixel’ (Hansen et al., 2014). For example, as depicted in Figure 7, the Hansen data indicates that the tree crown cover in 2000 in pixel 1 is 54%, meaning that 54% of the 30 m pixel is covered by vegetation taller than 5 m. The Hansen data also indicates that forest loss occurred in pixel in 2013, meaning that all of the tree canopy disappeared in 2013.

For every 30 m Landsat pixel within a buffer radius r (1, 10, and 30 km) of study’s villages, the tree crown cover density in 2000 and the year of forest loss were combined to derive the deforestation and forest cover variables.

The two upper plots highlight the raw data at two example pixels from the lower plot.

Deforestation variable

Request a detailed protocol

To define our primary exposure variable, for all villages in the study districts and year of the study period, we calculated the percent area within a buffer radius of 1, 10, and 30 km that experienced forest loss in the previous 1, 2, 3, 4, and 5 years (Figure 7). These distances were chosen to explore a range of spatial scales at which the forest environment may be differentially relevant for village-based populations and forest-goers. To explore potential interactions between deforestation and forest cover, we computed an alternate exposure variable, restricting to areas that both experienced forest loss and had a tree crown cover density above 68% and 87%. Those thresholds are limits of the inter-quartile range (IQR) of the distribution of tree crown cover density in any 30 m pixels within 10 km of study’s villages that experienced forest loss between 2000 and 2017. This alternate definition captures deforestation activities occurring in areas with denser forest cover.

Forest cover variable

Request a detailed protocol

We also combined the two Hansen layers to produce annual tree crown cover maps of the study districts, assuming no changes prior to the year of forest loss but setting to 0 the pixel tree crown cover density afterwards (Figure 7). For all villages in the study districts and year of the study period, we calculated the average tree crown cover density within a buffer radius of 1, 10, and 30 km and for 0, 1, 2, and 3 year lags. This is a secondary exposure, adjusted for in the primary analysis.

Environmental covariates

Request a detailed protocol

Village population sizes were needed to estimate monthly malaria incidence. 2005 and 2015 population estimates for the 491 villages of study districts were obtained from the national census Lao national census, 2020. The annual population growth rate (3.7%) was used to impute population values for two villages missing 2005 estimates and for two villages missing 2015 estimates. Then, village-level population growth rates were used to estimate villages’ population per year between 2008 and 2016, assuming linear annual growth rate (median = 1.7%, IQR = [0%; 4.5%]).

Altitude, temperature, rainfall, and access to health care were considered as potential village-level confounders of the relationship between malaria and forest cover factors. Travel time to closest health facility, computed for the treatment-seeking model, was used as a proxy for health care access and villages’ remoteness. Altitude was extracted from SRTM (Jarvis et al., 2008) 1 km resolution layers. Monthly average day and night temperature were extracted from MODIS 1 km resolution product (MOD11C3 Wan et al., 2015). Finally, monthly total rainfall was extracted from CHIRPS (Funk et al., 2014) 1 km resolution publicly available data. The average and standard deviation of the annual total precipitation and the average monthly temperature from the monthly time series was computed over the 2008–2012 period, which corresponds to the 5 year time period directly before our malaria data (2013–2016). This ‘long-term’ aggregation of the climatic variables is included in the model to capture the spatial differences in overall climate between the villages of our study area. To account for the seasonal effect of these climatic variables, monthly temperatures and precipitation in the previous 1, 2, and 3 months were also extracted, as well as the average temperatures and total precipitation over the previous 1, 2, and 3 months (seven ‘short-term’ variations: in current month, in previous 1, 2, or 3 months and aggregated over current and previous 1, 2, or 3 months). See ‘Details on covariates’ below.

Altitude was missing for one village and we used an online elevation finder tool (FreeMapTools) for imputation. Temperature was missing for 2.4% of the village-months over the study period, most likely because of cloud coverage of the MODIS imagery. Monthly temperature was never missing more than 2 years in a row at villages of the study’s districts and we imputed the temperature of the same month of the following year (or prior year when needed), adjusting for average district-level monthly temperature differences between the 2 consecutive years. Monthly rainfall was not missing at any of the villages.

Statistical analysis

Statistical model

Request a detailed protocol

To model malaria incidence (Equation 2), the number of positive cases Yv,t at village v over month t was modeled using a generalized additive model (GAM) (Wood, 2017). To account for overdispersion, a negative binomial distribution was used, including an additional variance parameter ν. The probability of seeking treatment θv, estimated from the treatment-seeking model, was multiplied by the village population Popv,t to derive the population seeking treatment, Popv,tseek. This was included as an offset term in the incidence model. Spatial autocorrelation was accounted for by the bivariate thin plate spline smoothing function on coordinates, f(Lat,Long) and village random intercepts were included. A non-linear temporal trend was also included with the smoothing function on month, f(t). Finally, the primary exposure, deforestation, and potential environmental confounders, including forest cover, were modeled with splines in f(Xv,ti). Splines add up polynomial basis functions in between knots and allow to control for very flexible relationships with covariates and spatio-temporal trends. Regularization was used to integrate model selection into the model fitting step by adding an extra penalty to each term so that the coefficients for covariates can be penalized to zero, also meaning that splines can be kept minimal if the data does not support more flexibility. See Figure 8 for a graphical visualization of our conceptual model for this analysis.

(2) Yv,tNegBin(E[Yv,t],ν)log(E[Yv,t])=log(μv,tPopv,tseek)=log(μv,t)+log(Popv,tθv)log(μv,t)=iβif(Xv,ti)+f(t)+f(Lat,Long)+βv

with βv𝒩(0,σβ).

Conceptual model for our analysis showing how the raw input data (blue boxes) were combined via intermediate data (white boxes) and models (white diamonds) to produce our estimated outputs (red circle).

We ran 15 models separately in the north and the south, each varying the buffer radius (1, 10, and 30 km) and temporal scale for deforestation (previous 1, 2, 3, 4, and 5 years). The coefficients of the linear effect for deforestation were extracted and exponentiated to get the incidence rate ratio (IRR) associated with a 0.1% increase in the percent area that experienced forest loss around villages.

Secondary analyses

Request a detailed protocol

As secondary analyses, we ran the same models, separately by malaria species (P. falciparum and P. vivax), leveraging the large amount of data and co-endemicity in the south. We also used our alternative definitions for deforestation, restricted to areas that both experienced forest loss and had a tree crown cover density above 68% and 87%, to explore the interaction between deforestation and the amount of forest cover.

To further strengthen the robustness of our analysis, we conducted a sensitivity analysis where villages’ populations in the surveillance system registries were not adjusted for the probability of seeking treatment. See Appendix 1 - S3.

Details on covariates

Request a detailed protocol

To prevent collinearity in the final model, for each of the three monthly climatic variables (precipitation, day and night temperature), we first selected the one of its seven ‘short-term’ variations (in current month, in previous 1, 2, or 3 months and aggregated over current and previous 1, 2, or 3 months) that provided the best AIC fit in an univariate model, solely adjusted for the spatio-temporal structure of the data (f(t), f(Lat,Long) and village random intercepts). See Appendix 1 - S6.

Malaria incidence in the previous 1 and 2 months were included in the model. Results in Appendix 1 - S4 show this was necessary to fully address temporal autocorrelation and led to a better AIC fit. Different shape of the temporal trend f(t) were also explored (up to 25 spline knots, auto-regressive, cyclic cubic spline) but none accounted for temporal autocorrelation better.

In a preliminary analysis, before including the deforestation variable, we ran our model in Equation 2 with forest cover as the primary exposure. We ran 12 models separately in the north and the south, each varying the buffer radius (1, 10, and 30 km) and temporal scale (0, 1, 2, and 3 year lag) for the forest cover variable. The coefficients of the linear effect for forest cover were extracted and exponentiated to get the incidence rate ratio (IRR) associated with a 1% increase in the average tree crown cover density around villages (Appendix 1—table 3). The model including average tree crown cover density within 30 km of villages with no temporal lag provided the smallest AIC value. In the final models with the deforestation variables we therefore included the average tree crown cover density within 30 km of villages in the starting year of the temporal scale for the deforestation variable considered (e.g. 3-year lag in the model with percent area that experienced forest loss in previous 3 years as the deforestation variable) to adjust for baseline forest cover.

Appendix 1

S1: Travel times methods

To calculate the travel time along a path linking any two points of the map, we defined a transition matrix that gives the speed at which one may travel between two adjacents pixels. We followed the parameterization suggested by Alegana et al., 2012 and demonstrated by Sturrock et al., 2014 (see Appendix 1—table 1), which first uses Toblers’ hiking function to specify the travel speed between two points of different altitudes. Intuitively, it is faster to travel downhill than uphill. Second, the speed is adjusted based on the type of landcover traveled through: a forested or a flooded area for instance slows you down. Last, the network of roads and major rivers may be used to catch a bus or a boat and therefore increases the travel speed. Altitude (SRTM 90 m, Jarvis et al., 2008) was aggregated and resampled at the land cover (ESA GlobCover 2009 Project European space agency, 2009) 300m-resolution and roads and waterways from Open Street Map (Open street map, 2020) were rasterized to calculate the transition matrix all across Lao PDR. The ‘raster’ package in R (R Development Core Team, 2020) was used.

Appendix 1—table 1
Data used to parameterized the transition matrix with the travel speed between any two adjacent pixels of the map.
Data layerCategorySpeed (km/h)
Digital elevation (slope)0°(flat)5
5°(uphill)3.71
−5°(downhill)5.27
Land coverCroplandNo adjustment
Artificial and bare areasNo adjustment
Open deciduous forest0.8 * Hiking speed
Sparse herbaceous0.8 * Hiking speed
Closed deciduous forest0.6 * Hiking speed
Herbaceous0.6 * Hiking speed
Flooded0.5 * Hiking speed
Other forest cover0.4 * Hiking speed
Water0.2 * Hiking speed
Roads and riversMotorway/trunk80
Primary/secondary60
Tertiary/unclassified10
Major rivers5

We then used the Djisktra’s algorithm from the R-package igraph (Csardi and Nepusz, 2006) and the gdistance package (van Etten, 2017) to find the fastest route between every village (or every household in the cross-sectional surveys) and its closest health facility. Coordinates of health facilities across Lao PDR came from the 2017 stratification exercise and were provided by CMPE. We authorized travel through non-study districts but not across international borders.

S2: Travel times and treatment-seeking results

Appendix 1—figure 1a shows how the travel time to closest health facility varies across Champasak province in southern Lao PDR, influenced by both distance and road connectivity. Appendix 1—figure 1c presents a right-skewed distribution of travel time from study villages to the closest health facility. Most villages are within 2 hr of the closest health facility but some are as far as 6 hr away. The distribution is similar for villages in the northern and southern study districts.

Appendix 1—figure 1
Treatment-seeking modeling plots. Note that treatment-seeking at public health facilities is implied all along the manuscript.

In the southern household survey, 243 individuals reported fever in the past 2 weeks. 225 (92.6%) of them, from 156 households, answered whether or not they sought treatment and were included in the treatment-seeking model. 219 (97.3%) reported seeking treatment and they all reported where they did so: 154 (70.3%) of them sought treatment at a public health facility (Village malaria worker (VMW), health center, district hospital or provincial hospital) and would therefore appear in the malaria registries collected.

In the northern household survey, 378 individuals reported fever in the past 2 weeks. 360 (95.2%) of them, from 297 households, answered whether or not they sought treatment and were included in the treatment-seeking model. A total of 283 (78.6%) reported seeking treatment. Only 40 (14.1%) of them reported where they did so but all of them sought treatment at a public health facility and we therefore upweighted the population that sought treatment at a public health facility accordingly.

Most surveyed households included in the treatment-seeking model were within 2 hr of travel time to the closest health facility but some were almost 5 hr away (Appendix 1—figure 12). Appendix 1—figure 1b shows the modeled relationship between the probability of seeking treatment (at a public health facility, implied from now on) and distance to the closest health facility. For villages within the same 300 m2 pixel as a health facility (estimated travel time of 0 hr), the predicted probability of seeking treatment was 0.87 (95% CI: [0.79; 0.92]) in the north and 0.78 (95% CI: [0.63; 0.89]) in the south. A 1 hr increase in travel time to the closest health facility was associated with a similar 0.79 (95% CI: [0.55; 1.13]) reduction in the odds of seeking treatment in the north and 0.76 (95% CI: [0.43; 1.34]) in the south, almost reaching statistical significance when pooling data from both regions: 0.77 (95% CI: [0.56; 1.04]). Appendix 1—figure 1d shows the resulting distribution for the probability of seeking treatment for all villages in study’s districts. Monthly village-level malaria incidence was adjusted accordingly.

S3: Sensitivity analysis

We conducted a sensitivity analysis where village population at risk of appearing in the surveillance system registries were not adjusted for the probability of seeking treatment. The effect estimates and confidence intervals were virtually unchanged, strengthening the robustness of our primary analysis (Appendix 1—table 2 below).

Appendix 1—table 2
IRR associated with a 0.1% increase in forest loss.

Adjusted for the spatio-temporal structure of the data, the environmental covariates selected in the model and forest cover within 30 km in the year before the deforestation temporal scale considered and malaria incidence in the previous 1 and 2 months. See Materials and methods for details. Sensitivity analysis: village population unadjusted for probability of seeking treatment.

SouthNorth
Time lagBuffer radiusBuffer radius
1 km10 km30 km1 km10 km30 km
Previous11.011.1611.031.01
1 year[0.99; 1.01][0.99; 1.04][1.10; 1.22][1; 1.01][1; 1.07][0.94; 1.08]
Previous111.0911.010.98
2 years[0.99; 1.01][0.98; 1.01][1.04; 1.13][1; 1.01][0.99; 1.04][0.94; 1.01]
Previous0.990.980.9311.010.96
3 years[0.99; 1][0.97; 1][0.90; 0.97][1; 1.01][0.99; 1.02][0.93; 0.99]
Previous0.990.980.94110.97
4 years[0.99; 1][0.97; 0.99][0.92; 0.97][1; 1.01][0.99; 1.02][0.94; 0.99]
Previous10.970.9411.010.95
5 years[0.99; 1][0.96; 0.99][0.91; 0.97][1; 1.01][0.99; 1.02][0.93; 0.98]

S4: Inclusion of malaria cases in previous months

Appendix 1—figure 2 shows residual temporal auto-correlation plots in models from equation 2, when malaria incidence in previous 1 and 2 months are included or not. These plots show that including covariates for malaria incidence in the previous 1 and 2 months is necessary to address residual temporal autocorrelation and keep each lag-wise individual autocorrelation estimate below 5%. The plots presented here are for the model in the south with a 30 km buffer radius and a 1-year temporal lag but similar results were observed across all 15 models both in the north and in the south.

Appendix 1—figure 2
Residual temporal autocorrelation when malaria incidence in previous 1 and 2 months are included or not.

The AIC fit also substantially improved from 18,667 when no malaria incidence is included to 18,152 when malaria incidence in previous month is included and to 17687 when both malaria in the previous 1 and 2 months are included.

S5: Additional results

Environmental covariates

Appendix 1—figure 3 shows the relationship – via their individual contribution β*f(X) in equation 2 – between malaria incidence and the environmental covariates included in the model (30 km radius and 1-year temporal lag). Thses plots show that relationships differ slightly by region although the range covered by the environmental variables also differs by region. We also see the effect of regularization, that penalized some covariates to zero, like our long-term precipitation covariates. This penalization happened more frequently in the north, where we had much less data. Note that 95% confidence intervals (see Appendix 1—figure 13) have been hidden for better visualization. The larger amount of data in the south also allowed the identification of more precise relationships than in the north.

Appendix 1—figure 3
Relationships between malaria incidence and the environmental covariates in the multivariable model described in equation 2 (30 km radius and 1-year temporal lag), additionally adjusted for the probability of seeking treatment, the spatio-temporal structure of the data (f(t), f(Lat,Long) and village random intercepts) and malaria incidence in the previous 1 and 2 months.

See Materials and methods for details. Note that 95% confidence intervals (see Appendix 1—figure 13) have been hidden for better visualization.

Temporal trend

Appendix 1—figure 4 shows the relationship – via its individual contribution β*f(t) in equation 2 – between malaria incidence and the temporal trend included in the model (30 km radius and 1-year temporal lag). These plots show that relationships are quite similar in both regions with an increase in 2014, followed by a plateau in 2015 and a decrease in 2016. The larger amount of data in the south also allowed the identification of a more precise relationship than in the north.

Appendix 1—figure 4
Relationships between malaria incidence and the temporal trend in the multivariable model described in equation 2 (30 km radius and - year temporal lag), additionally adjusted for the probability of seeking treatment, the spatial structure of the data (f(Lat,Long) and village random intercepts) and malaria incidence in the previous 1 and 2 months.

See Materials and methods for details.

Forest cover

Appendix 1—table 3 shows the incidence rate ratio (IRR) associated with forest cover, measured by a 1% increase in the average tree crown density, in current and previous 3 years within 1, 10, and 30 km of villages. Forest cover within 1 km of a village was not associated with malaria incidence rate in either the south or the north, regardless of the temporal lag. However within 10 and 30 km of a village, increased forest cover tended to be associated with higher malaria incidence rates both in the north and the south (e.g. 30 km buffer, 1-year lag, IRR = 1.09, 95% CI: [1.03; 1.15] in the south; IRR = 1.12, 95% CI: [0.99; 1.26] in the north). The associations were higher when considering a larger spatial scale (30 km) but were already statistically significant for a 10 km buffer radius in the south. None of the associations reached statistical significance in the north, where the sample size is small. The temporal scale considered did not affect the associations much. Statistical significance wasn’t necessarily reached for all the associations highlighted, but the trends observed suggest forest cover around villages but not in the immediate vicinity (1 km) leads to higher risk of malaria both in the north and in the south. The model including average tree crown cover density within 30 km of villages with no temporal lag provided the best AIC. In the final models with the deforestation variables we therefore included the average tree crown cover density within 30 km of villages in the starting year of the temporal scale for the deforestation variable considered (e.g. 3-year lag in the model with percent area that experienced forest loss in previous 3 years as the deforestation variable).

Appendix 1—table 3
IRR [95% CI] associated with a 1% increase in average tree crown density.

Adjusted for the probability of seeking treatment, the spatio-temporal structure of the data, the environmental covariates selected in the model and malaria incidence in the previous 1 and 2 months. See Materials and methods for details.

SouthNorth
Time lagBuffer radiusBuffer radius
1 km10 km30 km1 km10 km30 km
Current11.071.060.991.011.10
year[0.99; 1.01][1.04; 1.10][1; 1.12][0.97; 1.02][0.96; 1.05][0.99; 1.23]
Previous11.071.0911.011.12
1 year[0.99; 1.02][1.05; 1.10][1.03; 1.15][0.97; 1.02][0.97; 1.06][0.99; 1.26]
Previous11.071.0911.021.10
2 years[0.99; 1.02][1.05; 1.10][1.03; 1.16][0.98; 1.03][0.97; 1.06][0.98; 1.25]
Previous11.071.101.011.021.10
3 years[0.99; 1.02][1.04; 1.10][1.04; 1.16][0.98; 1.03][0.97; 1.07][0.98; 1.24]

Deforestation - non-linearities

The IRR effect estimates in Table 1 and Figure 3 in the main manuscript assume a linear relationship between deforestation and malaria. Appendix 1—figure 5 shows a few of these relationships – via their individual contribution β*f(X) in equation 2 – when such linearity is not imposed in the GAM models. Although the AIC fit is slightly better when modeling non-linearities, these plots show that the linearity assumption is mostly warranted.

Appendix 1—figure 5
Adjusted relationship between deforestation and malaria incidence.

All models were adjusted for environmental covariates and forest cover on top of the probability of seeking treatment, the spatio-temporal structure of the data (f(t), f(Lat,Long) and village random intercepts) and malaria incidence in the previous 1 and 2 months. See Materials and methods for details. Note that scales are different between buffer radius for better visualization. Appendix 1—figure 14 shows the raw scatterplot between monthly village malaria incidence rate and deforestation. Appendix 1—figure 15 and Appendix 1—figure 16 show the raw time series of malaria incidence, forest cover and percent area that experienced forest loss.

S6: AIC fit of the seven monthly climatic variables variations

To avoid collinearity, we have selected (based on the best AIC fit – see Appendix 1—table 4) the 1 of the 7 variations (In current month, in previous 1, 2, or 3 months and aggregated over previous 1, 2, or 3 months) of the three monthly climatic variables (Precipitation, Day temperature and Night temperature) to be included in the final model. This was done independently for each of the four outcome models (South, North, South Pf, and South Pv). For South P. falciparum, the second best fitting AIC day temperature (in current month) was selected (rather than in previous month) because of very similar AIC fits and to ensure better comparability with the overall South and South P. vivax models.

Appendix 1—table 4
AIC fit of univariate models when including each of the seven monthly climatic variation one at a time as unique covariate in equation 2, solely adjusted for the probability of seeking treatment, the spatio-temporal structure of the data (f(t), f(Lat,Long) and village random intercepts).

AIC selected are in bold.

Outcome model
SouthNorthSouth P. falciparumSouth P. vivax
Day temperature
Current month18,546167113,22614,575
Previous month18,556170213,22414,590
2 months ago18,578166913,24914,594
3 months ago18,559167213,23214,593
Over current and previous month18,556167013,23114,583
Over current and previous 2 months18,570167013,24814,588
Over current and previous 3 months18,573168013,24914,592
Night temperature
Current month18,413166913,12014,474
Previous month18,453167013,15514,520
2 months ago18,547167313,23114,576
3 months ago18,581167213,25114,596
Over current and previous month18,296166413,04414,397
Over current and previous 2 months18,263166913,01414,385
Over current and previous 3 months18,262166313,00714,385
Precipitation
Current month18,532169313,19814,593
Previous month18,520166913,18114,575
2 months ago18,538165813,20714,594
3 months ago18,579166413,24314,596
Over current and previous month18,570167213,23914,596
Over current and previous 2 months18,543167013,18714,590
Over current and previous 3 months18,555167413,21214,591

S7: Additional figures

This section presents additional figures mentioned in the text and in the additional results section of appendix 1.

Forest and environmental variables

Appendix 1—figures 6, Appendix 1—figure 7 and Appendix 1—figure 8.

Appendix 1—figure 6
Distribution of average tree crown cover density within 1, 10, and 30 km of villages.
Appendix 1—figure 7
Distribution of percent area within 1, 10, and 30 km of villages that experienced forest loss between 2011 and 2016.

Note that the scales are different for every panel for better visualization of the distributions.

Appendix 1—figure 8
Distribution and time series of environmental covariates (population, altitude, monthly day temperature and monthly total precipitation) at study’s villages.

Malaria registries – malaria infections

Appendix 1—figure 9.

Appendix 1—figure 9
Additional figures from malaria registries: malaria infections.

Malaria registries - SES

Appendix 1—figure 10.

Appendix 1—figure 10
Distributions of socio-economomical variables of all patients recorded in the malaria registries.

Malaria registries – matched vs unmatched

Appendix 1—figure 11.

Appendix 1—figure 11
Additional figures from malaria registries: matched vs unmatched SES variables.

Treatment-seeking

Appendix 1—figure 12.

Appendix 1—figure 12
Distribution of travel time (in hours) from surveyed households to closest health facilities.

Statistical analysis – environmental covariates

Appendix 1—figure 13.

Appendix 1—figure 13
Relationships between malaria incidence and the environmental covariates in the multivariable model described in equation 2 (30 km radius and 1-year temporal lag), additionally adjusted for the probability of seeking treatment, the spatio-temporal structure of the data (f(t), f(Lat,Long) and village random intercepts) and malaria incidence in the previous 1 and 2 months.

Dashed lines are for 95% confidence intervals. Note that the y scale has been trimmed a bit for better visualization.

Raw association between malaria incidence and deforestation

Appendix 1—figure 14.

Appendix 1—figure 14
Raw scatterplot between monthly village malaria incidence rate and the percent area within 30 km of villages that experienced forest loss in the previous 1, 3, and 5 years.

Note that scales are different between regions for better visualization.

Raw time series of malaria incidence, forest cover, and deforestation

Appendix 1—figure 15 and Appendix 1—figure 16.

Appendix 1—figure 15
Time series of deforestation (percent area that experienced forest loss around villages), forest cover (average tree crow cover around villages) and malaria incidence, averaged over study’s villages and for varying buffer radius around villages (1, 10, and 30 km).
Appendix 1—figure 16
Time series of deforestation (percent area that experienced forest loss within 30 km of villages) and forest cover (average tree crow cover within 30 km of villages), for a few randomly sampled study’s villages.

Each color represents one village.

Statistical analysis – P. falciparum and P. vivax

Appendix 1—figure 17.

Appendix 1—figure 17
Adjusted relationship between deforestation and species-specific malaria incidence in southern Lao PDR.

All models were adjusted for environmental covariates and forest cover on top of the probability of seeking treatment, the spatio-temporal structure of the data (f(t), f(Lat,Long) and village random intercepts) and malaria incidence in the previous 1 and 2 months.

Data availability

All data generated or analysed during this study are included in the manuscript and supporting files. Source data files have been provided for Figures 2, 3, 4 and 5 and for Tables 1, 2 and 3.

References

  1. Report
    1. Bennett A
    2. Lover A
    3. Dantzer EU
    (2021)
    Report on Formative Assessment to Identify and Characterize Mobile and Migrant Populations in Champasak Province, Southern Lao Pdr
    University of California, San Francisco (UCSF), Malaria Elimination Initiative, Lao PDR.
    1. Chaveepojnkamjorn W
    2. Pichainarong N
    (2004)
    Malaria infection among the migrant population along the thai-myanmar border area
    The Southeast Asian Journal of Tropical Medicine and Public Health 35:48–52.
    1. Csardi G
    2. Nepusz T
    (2006)
    The igraph software package for complex network research
    InterJournal, Complex Systems 1695:1–9.
    1. Das NG
    2. Talukdar PK
    3. Das SC
    (2004)
    Epidemiological and entomological aspects of malaria in forest-fringed villages of sonitpur district, Assam
    Journal of Vector Borne Diseases 41:5.
  2. Website
    1. Lao national census
    (2020) Lao decide
    Accessed January 15, 2020.
  3. Report
    1. Lao PDR
    (2016)
    Lao Pdr National Strategic Plan for Malaria Control and Elimination 2016-2020. Technical Report, Centre for Malariology, Parasitology, and Entomology (CMPE)
    Department of Communicable Disease Control, Lao PDR Ministry of Health. Vientiane, Lao PDR.
  4. Report
    1. Okayas H
    (2018)
    Mekong Malaria Elimination Programme
    World Health Organization, Global Malaria Programme, Geneva.
  5. Website
    1. Open street map
    (2020) Open street map
    Accessed January 15, 2020.
  6. Report
    1. Pattanayak SK
    2. Corey CG
    3. Lau YF
    4. Kramer RA
    (2010)
    Biodiversity Conservation and Child Malaria: Microeconomic Evidence From Flores, Indonesia
    Economic Research Initiatives at Duke Working Paper.
  7. Software
    1. R Development Core Team
    (2020) R: A Language and Environment for Statistical Computing
    R Foundation for Statistical Computing, Vienna, Austria.
  8. Report
    1. Sawyer DR
    (1988)
    Frontier Malaria in the Amazon Region of Brazil: Types of Malaria Situations and Some Implications for Control
    PHO/WHO/TDR.
  9. Conference
    1. Sharma VP
    2. Kondrashin AV
    (1991)
    Forest malaria in Southeast Asia
    Proceedings of an Informal Consultative Meeting WHO/MRC, WHO/Malaria Research Centre, New Delhi, 1991.
  10. Report
    1. UNFPA
    (2016)
    Results of Population and Housing Census 2015 (English Version)
    Technical report, Lao Statistics Bureau.
    1. Vittor AY
    2. Pan W
    3. Gilman RH
    4. Tielsch J
    5. Glass G
    6. Shields T
    7. Sánchez-Lozano W
    8. Pinedo VV
    9. Salas-Cobos E
    10. Flores S
    11. Patz JA
    (2009)
    Linking deforestation to malaria in the Amazon: characterization of the breeding habitat of the principal malaria vector, anopheles darlingi
    The American Journal of Tropical Medicine and Hygiene 81:5–12.
  11. Report
    1. WHO
    (2016)
    Approaches for Mobile and Migrant Populations in the Context of Malaria Multi-Drug Resistance and Malaria Elimination in the Greater Mekong Subregion
    Regional Office for South-East Asia.
  12. Report
    1. WHO
    (2018)
    Geneva: World Health Organization, Technical Report, WHO
    World Report Malaria.
  13. Report
    1. Wood SN
    (2017)
    Generalized Additive Models: An Introduction with R
    Chapman and Hall/CRC.
  14. Report
    1. World Health Organization
    (2016)
    Eliminating Malaria in the Greater Mekong Subregion: United to End a Deadly Disease
    Technical report, World Health Organization.

Decision letter

  1. Ben S Cooper
    Reviewing Editor; Mahidol University, Thailand
  2. Eduardo Franco
    Senior Editor; McGill University, Canada
  3. James A Watson
    Reviewer; Mahidol Oxford Tropical Medicine Research Unit, Thailand

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Thank you for submitting your article "Spatio-temporal associations between deforestation and malaria incidence in Lao PDR" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and a Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: James A Watson (Reviewer #1).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission. Please submit a revised version that addresses these concerns directly.

Summary:

Deforestation is widely thought to have played an important role in the declining incidence of malaria in the Greater Mekong Sub-region but detailed analyses have been lacking. In this paper the authors assemble a large dataset of geolocated malaria cases from Lao PDR together with forest cover data. Using a spatio-temporal modelling framework, they find no evidence of an association of deforestation within 1 or 10km of a village and malaria incidence. However, deforestation within 30km of a village was associated with malaria incidence, though with the direction of the association depended on the timescale considered. In the short term (1-2 years) deforestation was associated with increases in malaria incidence in the south, though no association was found in the north. In the long-term (3-5 years) deforestation was associated with a lower malaria incidence rate, with decreases of about 5% in both the south and the north.

Essential revisions:

1) All reviewers had concerns about variable selection and a number of suggestions were made, including regularisation, the need to adjust for climate anomalies and malaria control programme activities. The differences between the north and south were also highlighted and, following the consultation process, there was a consensus that lumping both parasites together and using one region did not seem appropriate. There were also questions about the choice of initial climate covariates. These issues need to be addressed either by conducting the analysis in a different way (we leave it to the authors to decide on the best way) or by convincing the reviewers that such changes are not needed. A better explanation and consideration of the climate covariates that were considered is also needed.

2) There needs to be a more thorough exploration of the effect of inclusion of malaria cases in the previous month. As reviewer 3 points out, it is surprising that the effect of malaria cases in the previous month is left as a sensitivity analysis as the autocorrelation effect should be important in a transmission system.

3) The concerns about the way forest data are used (as raised by reviewer 2) need to be addressed. The authors should justify how they have derived canopy cover from a dataset not designed to quantify annual forest canopy cover. This is quite a widely used dataset and the interpretation is unconventional. The authors could either re-analyse these data or provide further explanation of how they have modelled canopy cover levels.

4) The authors should include the nonlinear trend f(t) in the results as well as the overall temporal pattern in the regions' deforestation.

5) There needs to be greater clarity about what was done in the results – e.g. regarding the treatment seeking model and adjustment for confounders.

6) There is a need for clarification of differences in malaria regimes in north and south in the Introduction.

Reviewer #1:

This is a carefully written paper on the complex associations between deforestation and malaria incidence. For someone with no background in the relevant literature it provides an easy to read introduction to previous work and shows how these new results fit with previous reports. On the whole I enjoyed reading it. My main suggestions are cosmetic.

I feel that the authors have done themselves a disservice in the Results section when reporting the observed associations between deforestation and malaria incidence. When I first read it, I had the impression that no adjustment for confounders was made! This does not seem to be mentioned at all in the Results. But in fact the underlying model is quite complex (Figure 8) and the authors use elevation, precipitation and temperature as the possible confounders (if I understood correctly). I would suggest adding a few sentences in the Results to clarify this. Maybe the Discussion could have a sentence or two on any unmeasured/residual confounding? Another really important adjustment is done for treatment seeking behaviour. This is rarely done in similar stats analyses and the authors have put quite a lot of work put into constructing a sophisticated model of treatment seeking (Figure 7B is really informative! Relegating this figure to the supplementary materials is a shame.)

I don't think the Results section even mentions the treatment seeking model. It is an important output and when mentioned in the Results, it should be made clear that the treatment seeking model is built using a different dataset, so there is no doubling dipping when constructing the geospatial model.

I'm not very familiar with GAMs so feel free to ignore this comment: instead of doing forward/backward AIC-based variable selection, could you instead fit the full model with appropriate regularisation (ridge regression type but in a GAM framework)?

Final comment concerns P. vivax. It is to be expected that the short term associations are dampened as probably the majority of infections are relapses (as mentioned in the Discussion). So you could argue that seeing a dampened association provides some kind of reassurance that there isn't massive unmeasured confounding going on (I would worry if larger associations were seen for vivax). This argument doesn't really work for the long term associations though. Any reasons for why no long term associations are observed?

Reviewer #2:

This study addresses an important research gap in understanding the ecology of malaria within Lao PDR and assembles an impressive dataset of geolocated malaria cases. However, there are major limitations in this study which question the validity of results and the interpretation of data sources.

First, without adjusting for precipitation or temperature differences annually, it is impossible to interpret whether the time lags between deforestation and malaria incidence reflect ecological changes and forest loss or whether this is due to interannual variation in other environmental factors. This time period includes a major El Nino event with extensive forest fires reported across Southeast Asia, including in Lao PDR (e.g. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6520341/). Without adjusting for the climate anomalies during this time period, it is difficult to attribute changes in risk to deforestation. Additionally, this analysis would ideally include some estimates of malaria control program activities during these time periods and whether this also varied either temporally or spatially. The authors address bias in health system reporting but no other control activities.

Additionally, assuming that the authors are simply using the forest data obtained from Hansen et al., this is not analyzed in a method consistent with how it was produced. This dataset includes tree crown coverage for 2000 (as a percentage) and a forest loss layer which reflects the change in state of a pixel from forest (defined as 50% canopy cover) to non-forest (less than 50%). As the forest loss layer only reflects a change in state (this could be 100% canopy cover to 0% or 52% to 49%), this is can only be interpreted by classifying the initial forest layer using a 50% threshold and cannot be used to estimate canopy cover percentages. There is also a forest gain layer (reflecting the change of state from non-forest to forest) which is not discussed at all within this paper. Following discussions about the accuracy of these datasets, there have also been updates to this data and additional layers available on World Resources Institute which could significantly improve this classification. As well, as forest configuration and fragmentation has widely been associated with malaria incidence in forest settings, this work would be improved by including additional metrics reflecting patterns of deforestation rather than just total area within a circular buffer.

Reviewer #3:

This manuscript addresses an important question in the ecology of malaria, namely what is the impact of deforestation on malaria incidence. Although there have been multiple studies and reviews on the subject for other regions especially the Amazon, an overall picture remains elusive because of the complexity of an environmental impact that is non-stationary in time as proposed in the hypothesis of "frontier malaria". Statistical analyses most often lack the longitudinal data needed to disentangle effects of deforestation as a function of time, and to also consider the spatial scales at which such effects are manifested. The work of Castro and colleagues for the Amazon did consider the temporal axis and provided evidence for a transient increase followed by a decrease to low endemic malaria.

Here, the authors take advantage of a longitudinal and spatially-resolved malaria data set in Lao People's Democratic Republic, together with a forest data set from remote sensing, to investigate how deforestation influences malaria incidence by both P. falciparum and P. vivax at different spatial and temporal scales. They contrast results for two regions within Lao PDR, North and South, with distinct transmission regimes and populations. They provide evidence consistent with the "frontier malaria" hypothesis, albeit with shorter times for the turn-around, from an increasing to a decreasing trend, than reported for the Amazon. The positive and short-term positive effect on malaria incidence is stronger in the South than in the North, and for P. falciparum. These results are important in demonstrating the need to consider statistical analyses that carefully introduce time since deforestation, and in supporting a positive transient effect of deforestation on malaria that should inform both public health strategies and forest management.

I have some comments that should be addressed to make the analyses and their interpretation clearer:

1) The way climate covariates are incorporated may ultimately "work" to capture a complex set of effects in a statistical model, but it is not very convincing. In particular, model selection is conducted on the basis of non-specific cases and the South only. The South and North show very different transmission regimes not just because of different population sizes but because of different environments, especially altitude (and therefore temperatures). In the more endemic regime of the South, malaria cases are likely to reflect much less clear effects of climate covariates than in the low transmission, epidemic, North. In addition, aggregating cases for the different parasites is problematic as their transmission is affected by climate covariates in different ways, given the relapses of P. vivax which typically lead to a different seasonality and less sensitivity to climate drivers. P. falciparum shows a more epidemic behavior in the North as expected for higher altitudes/cooler temperatures. These differences make one wonder whether it is valid to conduct model selection first in terms of the South and for both parasites together, and then fix these variables in the rest of the analysis.

2) The choice of initial climate covariates is also somewhat confusing. Why use WorldClim at all? One can obtain the means, or totals, and the CVs, from the time series used for the monthly data. WordClim uses a particular window of time (decades) to obtain these "typical" variables. It is my understanding that this time period is earlier than that considered here. Why use altitude and temperature? Aren't these strongly correlated within a region? An alternative consideration is to separate interannual and seasonal effects, by considering monthly values (as done here) and also the mean temperatures (of that same data set) over a window of time that precedes the transmission season (for P. falciparum) and corresponds to the rainy season.

3) The effect of malaria cases in the previous month is left for the sensitivity analysis as one variation of the main results. This is surprising as this autocorrelation effect should be important in a transmission system. That is, it would be natural to start with that model and see what is the significance of that variable. Its inclusion seems to do more than just weaken slightly the results. The effect in the North for P. falciparum becomes non-significant and that in the South weakens considerably. I wonder what would have been the result of model selection to start with, if this variable had been included. (I understood that the "adjusted" model does not repeat model selection).

4) The distinction of the malaria regimes in the North and South should be described more clearly in the Introduction for a general reader not familiar with malaria dynamics. The North appears seasonal low transmission or epidemic (given the altitudes and the time series), and the South, seasonal endemic. The differences in transmission characteristics between parasites should also be included early in the text; relapses are mentioned late in the Discussion. These considerations are important to how one looks at the time series and results.

5) The nonlinear trend f(t) is mentioned in the model description but not in the results. I could not tell whether deforestation itself exhibits a trend and how consideration or not of f(t) influences the results. Including the resulting f(t) as well as the overall temporal pattern in the regions' deforestation would be informative.

6) Similarly, I would have liked to see in the supplement the main results underlying variable selection and model selection.

https://doi.org/10.7554/eLife.56974.sa1

Author response

Essential revisions:

1) All reviewers had concerns about variable selection and a number of suggestions were made, including regularisation, the need to adjust for climate anomalies and malaria control programme activities. The differences between the north and south were also highlighted and, following the consultation process, there was a consensus that lumping both parasites together and using one region did not seem appropriate. There were also questions about the choice of initial climate covariates. These issues need to be addressed either by conducting the analysis in a different way (we leave it to the authors to decide on the best way) or by convincing the reviewers that such changes are not needed. A better explanation and consideration of the climate covariates that were considered is also needed.

We realize that our model selection approach could be improved based on the reviewers’ suggestions. We have now conducted this analysis in a different way by integrating model selection into the model fitting step (using regularization) without first lumping parasites or regions together.

First, we have selected the best fitting (based on AIC) of the 7 variations (In current month, in previous 1, 2 or 3 months and aggregated over previous 1, 2 or 3 months) for the 3 monthly climatic variables (Precipitation, Day temperature and Night temperature) independently for each of the 4 outcome models (South, North, South Pf and South Pv) to be included for further model selection. This change is described in the Materials and methods section and in Appendix 1—table 4 in supplementary Materials S6 in the appendix.

Second, we have used regularization to integrate model selection into the model fitting step by adding in the GAM an extra penalty to each term so that the coefficients for covariates can be penalized to zero (using the select=TRUE argument in the R mgcv::gam function). We thank the reviewers for this great suggestion, as we agree that it makes our model selection more rigorous and transparent while still offering great comparability between the 4 outcome models (South, North, South Pf and South Pv). This change is described in the Materials and methods section.

Third, we revisited our choice of initial climate variables. As suggested by reviewers, Worldclim variables were removed. Instead, we computed the average and standard deviation of the annual total precipitation and the average monthly temperature for the 2008-2012 years from the monthly time series (CHIRPS and MODIS respectively) already used to extract monthly climatic variables. The 2008-2012 period corresponds to the 5-year time period directly before our malaria data (2013-2016). Including the effect of “long-term” climatic variables in our model, on top of the “short-term” seasonal effect modeled via the monthly climatic variables, is important to capture the spatial differences in overall climate between the villages of our study area. This change is described in the Materials and methods section.

Fourth, we modified the result section to clarify that all of our analyses were controlled for climate variables. Previously, that was only mentioned in the labels of the tables and plots in the result section. We thank the reviewers for pointing out this important caveat in our original manuscript. This change is described in the Results section. In addition, and as suggested by reviewer 1, we also moved the figure of our conceptual model from the appendix into the Materials and methods section (Figure 8) to stress that all of our models were adjusted for environmental variables.

Finally, the reviewers are right to mention the importance of adjusting for malaria control program activities. When designing the study, in collaboration with the national control program, we purposefully excluded regions where we knew large programmatic activities were being implemented. Unfortunately though, there are no good data on coverage of other control program activities in these villages over time and we couldn’t formally adjust for it. That said, to be of importance, this residual unmeasured confounding (spatio-temporal variations in malaria control program activities) would need to be associated with deforestation. We have no reason to believe so, except maybe because of an area’s remoteness. We tried to adjust for this remoteness by including travel time to the closest health facility in the model. We included this additional information in the Materials and methods section. Additionally, we included a sentence about unmeasured residual confounding in the Discussion.

2) There needs to be a more thorough exploration of the effect of inclusion of malaria cases in the previous month. As reviewer 3 points out, it is surprising that the effect of malaria cases in the previous month is left as a sensitivity analysis as the autocorrelation effect should be important in a transmission system.

The reviewers are right to raise that point. In our new analysis, we included malaria cases in the previous month in the model and let regularization penalize it just as any other covariates. In this new analysis, we actually realized that malaria cases in the previous 2 months was necessary as well to fully account for any residual temporal autocorrelation and was therefore also included. Several different forms of the non-linear temporal trend f(t) were also explored (up to 25 spline knots, auto-regressive, cyclic cubic spline) but none accounted for the temporal autocorrelation better. The exploration of the effect of inclusion of malaria cases in the previous 1 and 2 months on residual autocorrelation is now presented in the appendix S4. This change is described in the Materials and methods section.

3) The concerns about the way forest data are used (as raised by reviewer 2) need to be addressed. The authors should justify how they have derived canopy cover from a dataset not designed to quantify annual forest canopy cover. This is quite a widely used dataset and the interpretation is unconventional. The authors could either re-analyse these data or provide further explanation of how they have modelled canopy cover levels.

We thank the reviewers for raising that point and asking for clarifications. This is a formidable dataset and we are grateful to the research team that produced it. We have spent quite some time reviewing documentation in the original article, the supplementary materials of the original article, authors’ replies to published critics of their paper and a Google Earth Engine tutorial for using Hansen et al. global forest cover and change data, developed by the authors themselves.

Reviewer 2 is correct that the dataset includes tree crown coverage for 2000 (as a percentage) and a forest loss layer which reflects the change in state of a pixel from forest (defined as 50% canopy cover) to a non-forest state. On the other hand, our interpretation of a non-forest state differs from what the reviewer suggests as follows:

In the Materials and methods section of Hansen’s paper (see supplementary materials), we read:

“Forest loss was defined as a stand-replacement disturbance. Results were disaggregated by reference percent tree cover stratum (e.g. >50% crown cover to ~0% crown cover) and by year.”

In the authors’ replies to published critics of their original paper, we further read:

“Forest loss was defined as a “stand-replacement disturbance”, meaning the removal or mortality of all tree cover in a Landsat pixel.”

Last, in the tutorial developed by the authors of these data, we read:

“The treecover2000 band […] provides information about the tree cover around the globe in the year 2000. The lossyear band provides information about if and when each 30-meter pixel on the planet was deforested during the study period. Given these two bands, it is possible to create an image of the state of the world's forest at any year in the study period.” Followed by code setting to 0 the tree crown the years after forest loss occurred.

Hence our understanding of the data is that the forest loss layer reflects the change in state of a pixel from forest (defined as 50% canopy cover) to a non-forest state (defined as 0% canopy cover).

With that said, we want to stress that our primary exposure, deforestation, quantifies the occurrence of forest loss rather than its magnitude. As a result, the important distinction discussed above has no impact on our primary exposure. Its impact only pertains to one of our covariates, tree crown cover, the value of which in a given year relies on forest loss leading to a 0% tree crown cover, as depicted in Figure 7 of our manuscript.

We have further replied in detail to reviewer 2’s other comments regarding our use of these data below, when answering their concerns directly. In light of the points raised, we have modified the Discussion to mention these limitations.

4) The authors should include the nonlinear trend f(t) in the results as well as the overall temporal pattern in the regions' deforestation.

We thank the reviewers for these suggestions. We have added Appendix 1—figure 4 in the “Additional Results” supplementary material S5 in the appendix, showing the non-linear trend f(t). Appendix 1—figure 15 and Appendix 1—figure 16 in the “Additional figures” supplementary material S6 in the appendix show the overall temporal pattern in the regions’ deforestation.

5) There needs to be greater clarity about what was done in the results – e.g. regarding the treatment seeking model and adjustment for confounders.

We thank the reviewers for pointing out this important caveat in our original manuscript. We have added a few sentences about the treatment seeking results in the Results section. In addition, and as suggested by reviewer 1, we also moved the figure of our conceptual model from the appendix into the Materials and methods section (Figure 8) to stress that all of our models were adjusted for confounders and highlight the importance of our treatment-seeking model. Last, we included further details in labels of all tables and figures.

6) There is a need for clarification of differences in malaria regimes in north and south in the Introduction.

We thank the reviewers for this important comment. We have modified the Introduction to highlight the different malaria regimes in the North and South. We also included additional reference to the different regional transmission regimes in the Results section.

Reviewer #1:

This is a carefully written paper on the complex associations between deforestation and malaria incidence. For someone with no background in the relevant literature it provides an easy to read introduction to previous work and shows how these new results fit with previous reports. On the whole I enjoyed reading it. My main suggestions are cosmetic.

I feel that the authors have done themselves a disservice in the Results section when reporting the observed associations between deforestation and malaria incidence. When I first read it, I had the impression that no adjustment for confounders was made! This does not seem to be mentioned at all in the Results. But in fact the underlying model is quite complex (Figure 8) and the authors use elevation, precipitation and temperature as the possible confounders (if I understood correctly). I would suggest adding a few sentences in the Results to clarify this. Maybe the Discussion could have a sentence or two on any unmeasured/residual confounding? Another really important adjustment is done for treatment seeking behaviour. This is rarely done in similar stats analyses and the authors have put quite a lot of work put into constructing a sophisticated model of treatment seeking (Figure 7B is really informative! Relegating this figure to the supplementary materials is a shame.)

I don't think the Results section even mentions the treatment seeking model. It is an important output and when mentioned in the Results, it should be made clear that the treatment seeking model is built using a different dataset, so there is no doubling dipping when constructing the geospatial model.

We thank the reviewer for pointing out these important caveats in our original manuscript.

We modified the result section to clarify that all of our analyses were controlled for climate variables. Previously, that was only mentioned in the labels of the tables and plots in the result section. This change is described in the manuscript.

We have also added a few sentences in a new “Treatment seeking” section in the Results section. We appreciate the reviewer’s enthusiasm and appreciation for that piece of work. We present it mainly in the appendix not to distract from the core analysis.

In addition, and as suggested below, we also moved the figure of our conceptual model from the appendix into the Materials and methods section (now Figure 8) to stress that all of our models were adjusted for confounders and highlight the importance of our treatment-seeking model. Last, we included further details in labels of all tables and figures.

Finally, we included a sentence about unmeasured residual confounding in the Discussion.

I'm not very familiar with GAMs so feel free to ignore this comment: instead of doing forward/backward AIC-based variable selection, could you instead fit the full model with appropriate regularisation (ridge regression type but in a GAM framework)?

The response to this comment is addressed in the third paragraph of our response to the editor’s essential revision comment #1.

Final comment concerns P. vivax. It is to be expected that the short term associations are dampened as probably the majority of infections are relapses (as mentioned in the Discussion). So you could argue that seeing a dampened association provides some kind of reassurance that there isn't massive unmeasured confounding going on (I would worry if larger associations were seen for vivax). This argument doesn't really work for the long term associations though. Any reasons for why no long term associations are observed?

We thank the reviewer for this constructive feedback. As mentioned in the Discussion, another recent study in the Amazon found similar attenuation of effects for vivax and we agree with the reviewer’s rational.

We think the same argument works for the long-term associations. Vivax long term infections, that can remain for years, when relapsing attenuate the decrease in incidence associated with deforestation in the long-term and pull the estimates towards 1. A study looking at additional larger long-term time lags for vivax might help.

Reviewer #2:

This study addresses an important research gap in understanding the ecology of malaria within Lao PDR and assembles an impressive dataset of geolocated malaria cases. However, there are major limitations in this study which question the validity of results and the interpretation of data sources.

First, without adjusting for precipitation or temperature differences annually, it is impossible to interpret whether the time lags between deforestation and malaria incidence reflect ecological changes and forest loss or whether this is due to interannual variation in other environmental factors. This time period includes a major El Nino event with extensive forest fires reported across Southeast Asia, including in Lao PDR (e.g. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6520341/). Without adjusting for the climate anomalies during this time period, it is difficult to attribute changes in risk to deforestation.

We thank the reviewer for raising that point and asking for clarifications. As raised by the other reviewers, it appears clearly that our original manuscript failed to convey the key information that all of our models were adjusted for environmental covariates, including precipitation and temperature.

We modified the result section to clarify that all of our analyses were controlled for climate variables including altitude as well as short term (monthly) and long-term (5-year average) precipitation and day and night temperature. Previously, that was only mentioned in the labels of the tables and plots in the result section. This change is described in the Results section. Additionally, and as suggested by reviewer 1, we moved the figure of our conceptual model from the appendix into the Materials and methods section (now Figure 8) to stress that all of our models were adjusted for confounders. Last, we included further details in labels of all tables and figures.

The reviewer is right to point out that our study does not distinguish between different causes of deforestation. We thank the reviewer for sharing additional reference, bringing to light the importance of natural fires in the GMS. The paper actually shows no particular changes in the occurrence of forest fires in Lao PDR over the study period, but we have to acknowledge that the Hansen forest data used does not distinguish man-made from natural causes of deforestation. We included that limitation in the Discussion.

Additionally, this analysis would ideally include some estimates of malaria control program activities during these time periods and whether this also varied either temporally or spatially. The authors address bias in health system reporting but no other control activities.

The response to this comment is addressed in the sixth paragraph of our response to the editor’s essential revision comment #1.

Additionally, assuming that the authors are simply using the forest data obtained from Hansen et al., this is not analyzed in a method consistent with how it was produced. This dataset includes tree crown coverage for 2000 (as a percentage) and a forest loss layer which reflects the change in state of a pixel from forest (defined as 50% canopy cover) to non-forest (less than 50%). As the forest loss layer only reflects a change in state (this could be 100% canopy cover to 0% or 52% to 49%), this is can only be interpreted by classifying the initial forest layer using a 50% threshold and cannot be used to estimate canopy cover percentages.

The response to this comment is addressed in our response to the editor’s essential revision comment #3.

There is also a forest gain layer (reflecting the change of state from non-forest to forest) which is not discussed at all within this paper.

The reviewer is right to mention the forest gain layer, present in the data, but not used nor mentioned in our paper.

Unlike the forest loss layer, which indicates the year in which forest loss occurred, the forest gain layer is simply a binary variable that takes the value 1 if forest gain (“defined as the inverse of loss or the establishment of tree cover from a non-treed state within a Landsat pixel”, as described in the authors’ replies) occurred between 2000 and 2017 and 0 if not. That is, we don’t know when it happens. In addition of lacking temporal resolution, the forest gain layer does not specify the percent tree crown cover reached after forest gain occurred. One of our key covariates, tree crown cover, would have benefited from information on forest gain but the layer present in the data is not informative enough.

On the other hand, our primary exposure, deforestation, is about the occurrence of forest loss which is fully quantified in the forest loss layer and doesn’t need the forest gain layer.

With that said, the assumptions behind our forest variables, presented in Figure 7, and related to both our understanding of the forest loss layer and the non-use of the forest gain layer but also the raw data’s limitation are key, and we thank the reviewer for asking for clarifications. We have modified the Discussion to mention this limitation.

Following discussions about the accuracy of these datasets, there have also been updates to this data and additional layers available on World Resources Institute which could significantly improve this classification. As well, as forest configuration and fragmentation has widely been associated with malaria incidence in forest settings, this work would be improved by including additional metrics reflecting patterns of deforestation rather than just total area within a circular buffer.

We have used the latest dataset available at the time (Summer 2019) which was an updated version of the dataset, first published in 2013. The reviewer is correct to mention alternative definitions of deforestation such as forest configuration and fragmentation, which have been used in some of the references cited in our paper. We understand that adding other definitions of deforestation would strengthen our work, but we feel this goes beyond the scope of our analysis, which is primarily focused on the exploration of associations over varying spatial and temporal scales. Also, please note that in Table 3 and Figure 5, we have explored one alternative definition to our primary deforestation exposure, restricting deforestation to the occurrence of forest loss in densely forested area.

Reviewer #3:

This manuscript addresses an important question in the ecology of malaria, namely what is the impact of deforestation on malaria incidence. Although there have been multiple studies and reviews on the subject for other regions especially the Amazon, an overall picture remains elusive because of the complexity of an environmental impact that is non-stationary in time as proposed in the hypothesis of "frontier malaria". Statistical analyses most often lack the longitudinal data needed to disentangle effects of deforestation as a function of time, and to also consider the spatial scales at which such effects are manifested. The work of Castro and colleagues for the Amazon did consider the temporal axis and provided evidence for a transient increase followed by a decrease to low endemic malaria.

Here, the authors take advantage of a longitudinal and spatially-resolved malaria data set in Lao People's Democratic Republic, together with a forest data set from remote sensing, to investigate how deforestation influences malaria incidence by both P. falciparum and P. vivax at different spatial and temporal scales. They contrast results for two regions within Lao PDR, North and South, with distinct transmission regimes and populations. They provide evidence consistent with the "frontier malaria" hypothesis, albeit with shorter times for the turn-around, from an increasing to a decreasing trend, than reported for the Amazon. The positive and short-term positive effect on malaria incidence is stronger in the South than in the North, and for P. falciparum. These results are important in demonstrating the need to consider statistical analyses that carefully introduce time since deforestation, and in supporting a positive transient effect of deforestation on malaria that should inform both public health strategies and forest management.

I have some comments that should be addressed to make the analyses and their interpretation clearer:

1) The way climate covariates are incorporated may ultimately "work" to capture a complex set of effects in a statistical model, but it is not very convincing. In particular, model selection is conducted on the basis of non-specific cases and the South only. The South and North show very different transmission regimes not just because of different population sizes but because of different environments, especially altitude (and therefore temperatures). In the more endemic regime of the South, malaria cases are likely to reflect much less clear effects of climate covariates than in the low transmission, epidemic, North. In addition, aggregating cases for the different parasites is problematic as their transmission is affected by climate covariates in different ways, given the relapses of P. vivax which typically lead to a different seasonality and less sensitivity to climate drivers. P. falciparum shows a more epidemic behavior in the North as expected for higher altitudes/cooler temperatures. These differences make one wonder whether it is valid to conduct model selection first in terms of the South and for both parasites together, and then fix these variables in the rest of the analysis.

The response to this comment is addressed in the first three paragraphs of our response to the editor’s essential revision comment #1.

(2) The choice of initial climate covariates is also somewhat confusing. Why use WorldClim at all? One can obtain the means, or totals, and the CVs, from the time series used for the monthly data. WordClim uses a particular window of time (decades) to obtain these "typical" variables. It is my understanding that this time period is earlier than that considered here. Why use altitude and temperature? Aren't these strongly correlated within a region? An alternative consideration is to separate interannual and seasonal effects, by considering monthly values (as done here) and also the mean temperatures (of that same data set) over a window of time that precedes the transmission season (for P. falciparum) and corresponds to the rainy season.

The response to this comment is addressed in the fourth paragraph of our response to the editor’s essential revision comment #1.

In addition, the reviewer is right in that temperature and altitude are not independent of each other but the landscape in northern Lao PDR is very craggy and much more variable than temperature. Both of them are included as initial covariates and regularization may penalized them to zero if needed.

(3) The effect of malaria cases in the previous month is left for the sensitivity analysis as one variation of the main results. This is surprising as this autocorrelation effect should be important in a transmission system. That is, it would be natural to start with that model and see what is the significance of that variable. Its inclusion seems to do more than just weaken slightly the results. The effect in the North for P. falciparum becomes non-significant and that in the South weakens considerably. I wonder what would have been the result of model selection to start with, if this variable had been included. (I understood that the "adjusted" model does not repeat model selection).

The response to this comment is addressed in our response to the editor’s essential revision comment #2.

(4) The distinction of the malaria regimes in the North and South should be described more clearly in the Introduction for a general reader not familiar with malaria dynamics. The North appears seasonal low transmission or epidemic (given the altitudes and the time series), and the South, seasonal endemic. The differences in transmission characteristics between parasites should also be included early in the text; relapses are mentioned late in the Discussion. These considerations are important to how one looks at the time series and results.

The response to this comment is addressed in our response to the editor’s essential revision comment #6.

5) The nonlinear trend f(t) is mentioned in the model description but not in the results. I could not tell whether deforestation itself exhibits a trend and how consideration or not of f(t) influences the results. Including the resulting f(t) as well as the overall temporal pattern in the regions' deforestation would be informative.

The response to this comment is addressed in our response to the editor’s essential revision comment #4.

6) Similarly, I would have liked to see in the supplement the main results underlying variable selection and model selection.

Regularization has now integrated model selection and model fitting into one step, simplifying our methodology and making it more rigorous and transparent. See Appendix 1—table 4 in supplementary Materials S6 in the appendix for results on the monthly climatic variable selection.

https://doi.org/10.7554/eLife.56974.sa2

Article and author information

Author details

  1. Francois Rerolle

    1. Malaria Elimination Initiative, The Global Health Group, University of California, San Francisco, San Francisco, United States
    2. Department of Epidemiology and Biostatistics, University of California, San Francisco, San Francisco, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing
    For correspondence
    francois.rerolle@ucsf.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3837-5700
  2. Emily Dantzer

    Malaria Elimination Initiative, The Global Health Group, University of California, San Francisco, San Francisco, United States
    Contribution
    Data curation, Project administration, Writing - review and editing
    Competing interests
    No competing interests declared
  3. Andrew A Lover

    Department of Biostatistics and Epidemiology, School of Public Health and Health Sciences, University of Massachusetts Amherst, Amherst, United States
    Contribution
    Data curation, Funding acquisition, Project administration, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2181-3559
  4. John M Marshall

    Division of Epidemiology and Biostatistics, School of Public Health, University of California, Berkeley, Berkeley, United States
    Contribution
    Supervision, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0603-7341
  5. Bouasy Hongvanthong

    Center for Malariology, Parasitology and Entomology, Ministry of Health, Vientiane, Lao People's Democratic Republic
    Contribution
    Project administration, Writing - review and editing
    Competing interests
    No competing interests declared
  6. Hugh JW Sturrock

    1. Malaria Elimination Initiative, The Global Health Group, University of California, San Francisco, San Francisco, United States
    2. Department of Epidemiology and Biostatistics, University of California, San Francisco, San Francisco, United States
    Contribution
    Supervision, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
  7. Adam Bennett

    1. Malaria Elimination Initiative, The Global Health Group, University of California, San Francisco, San Francisco, United States
    2. Department of Epidemiology and Biostatistics, University of California, San Francisco, San Francisco, United States
    Contribution
    Conceptualization, Supervision, Funding acquisition, Methodology, Project administration, Writing - review and editing
    Competing interests
    No competing interests declared

Funding

Bill and Melinda Gates Foundation (OPP1116450)

  • Francois Rerolle
  • Emily Dantzer
  • Andrew A Lover
  • Bouasy Hongvanthong
  • Hugh JW Sturrock
  • Adam Bennett

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

For their expertise and assistance, we thank Michelle Roh (Study design and interpretation), Ricardo Andrade Pacheco (GAM and spatial modeling), Alemayehu Midekisa (Remote sensing), Stephen Shiboski (Statistical analyses) and Maria Glymour (Causal inference and manuscript writing). FR was funded was the Bill and Melinda Gates Foundation (Grant ID OPP1116450).

Ethics

Human subjects: This study was approved by the National Ethics Committee for Health Research at the Lao Ministry of Health (Approval #2016-014; 8/22/2016) and by the UCSF ethical review board (Approvals #16-19649 and #17-22577). The informed consent process was consistent with local norms, and all study areas had a consultation meeting with, and approvals from, village elders. All participants provided informed written consent; caregivers provided consent for all children under 18, and all children aged 10 and above also provided consent directly. The study was conducted according to the ethical principles of the Declaration of Helsinki of October 2002.

Senior Editor

  1. Eduardo Franco, McGill University, Canada

Reviewing Editor

  1. Ben S Cooper, Mahidol University, Thailand

Reviewer

  1. James A Watson, Mahidol Oxford Tropical Medicine Research Unit, Thailand

Publication history

  1. Received: March 16, 2020
  2. Accepted: February 19, 2021
  3. Accepted Manuscript published: March 9, 2021 (version 1)
  4. Version of Record published: April 6, 2021 (version 2)

Copyright

© 2021, Rerolle et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 971
    Page views
  • 113
    Downloads
  • 1
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

  1. Further reading

Further reading

  1. How does cutting down forests influence the spread of malaria?

    1. Epidemiology and Global Health
    Benjamin Amoah et al.
    Research Article

    Background:

    Monitoring malaria transmission is a critical component of efforts to achieve targets for elimination and eradication. Two commonly monitored metrics of transmission intensity are parasite prevalence (PR) and the entomological inoculation rate (EIR). Comparing the spatial and temporal variations in the PR and EIR of a given geographical region and modelling the relationship between the two metrics may provide a fuller picture of the malaria epidemiology of the region to inform control activities.

    Methods:

    Using geostatistical methods, we compare the spatial and temporal patterns of Plasmodium falciparum EIR and PR using data collected over 38 months in a rural area of Malawi. We then quantify the relationship between EIR and PR by using empirical and mechanistic statistical models.

    Results:

    Hotspots identified through the EIR and PR partly overlapped during high transmission seasons but not during low transmission seasons. The estimated relationship showed a 1-month delayed effect of EIR on PR such that at lower levels of EIR, increases in EIR are associated with rapid rise in PR, whereas at higher levels of EIR, changes in EIR do not translate into notable changes in PR.

    Conclusions:

    Our study emphasises the need for integrated malaria control strategies that combine vector and human host managements monitored by both entomological and parasitaemia indices.

    Funding:

    This work was supported by Stichting Dioraphte grant number 13050800.