1. Introduction
Forests provide a myriad of ecosystem services (Millennium Ecosystem Assessment, 2005), including, in many settings, the mitigation of diseases like malaria by reducing suitable mosquito habitats and by housing species that feed on mosquitos. Malaria causes the loss of approximately 46.5 million disability adjusted life-years (DALYs) per year worldwide (Millennium Ecosystem Assessment, 2005). In developing countries with widespread malaria prevalence, protecting forests may be an effective way to reduce the burden of the disease (e.g., Berazneva and Byker, Reference Berazneva and Byker2017, Reference Berazneva and Byker2024; Garg, Reference Garg2019).
While forests can help mitigate the malaria burden, the impact and effectiveness depend on both biophysical and socioeconomic contexts (Yasuoka and Levins, Reference Yasuoka and Levins2007; Pattanayak and Yasuoka, Reference Pattanayak, Yasuoka and Colfer2008; Pattanayak and Pfaff, Reference Pattanayak and Pfaff2009). For example, though some mosquito species prefer more shade, in aggregate more sunlight is associated with increases in mosquito densities (Yasuoka and Levins, Reference Yasuoka and Levins2007). This is either because of the increased suitability for malaria vectors or due to the decrease in species (e.g., bats, dragonflies) that feed on mosquitoes and help control their populations (Pattanayak and Pfaff, Reference Pattanayak and Pfaff2009). Undisturbed old-growth forests with dense canopies are more likely to mitigate malaria than disturbed forests as the latter tend to have changed microclimatic conditions favoring mosquitoes (Yasuoka and Levins, Reference Yasuoka and Levins2007; Pattanayak and Yasuoka, Reference Pattanayak, Yasuoka and Colfer2008). Households who work in agricultural fields close to disturbed forests are more likely to be exposed to malaria-carrying mosquitoes (Pattanayak and Pfaff, Reference Pattanayak and Pfaff2009). Further, socioeconomic conditions mediate the ability of households to mitigate their own susceptibility to the disease. For example, poorer households are less able to engage in averting behavior or seek treatment upon infection (Wangdi et al., Reference Wangdi, Gatton, Kelly, Banwell, Dev and Clements2016; Busch and Ferretti-Gallon, Reference Busch and Ferretti-Gallon2017).
Empirical evidence on the role of forests in mitigating malaria is mixed (Colfer et al., Reference Colfer, Sheil and Kishi2006): while a large number of studies find that forests reduce the incidence and spread of malaria (e.g., Berazneva and Byker, Reference Berazneva and Byker2017, Reference Berazneva and Byker2024; Garg, Reference Garg2019), some find no statistically significant effect (e.g., Bauhoff and Busch, Reference Bauhoff and Busch2020) or a negative relationship (e.g., Valle and Clark, Reference Valle and Clark2013). The prevalence of malaria is associated with higher levels of forest cover, but also higher levels of forest loss in Malaysia (Fornace et al., Reference Fornace, Abidin, Alexander, Brock, Grigg, Murphy, William, Menon, Drakeley and Cox2016). In fact, previous studies have hypothesized an inverted U-shaped relationship between deforestation and malaria incidence, where small clearings can create beneficial conditions for mosquito larvae; larger clear cuts and development for ranching, agriculture or urbanization result in reduced larvae and therefore reduced malaria incidence (de Castro et al., Reference de Castro, Monte-Mór, Sawyer and Singer2006), a hypothesis supported empirically in the context of Brazil (Santos and Almeida, Reference Santos and Almeida2018). For these reasons, the magnitude of the impact of forests in mitigating malaria may depend on the context and remains an empirical question. Further, the complicated biophysical and socioeconomic nature of the determinants of malaria incidence underscores the need for careful analysis of the relationship between forests, deforestation, and malaria in each setting of interest.
Efforts to quantify the monetary value of forests in mitigating malaria are still limited (Ferraro et al., Reference Ferraro, Lawlor, Mullan and Pattanayak2012; Pattanayak et al., Reference Pattanayak, Kramer and Vincent2017): to our knowledge, Garg (Reference Garg2019) offers the only estimate of the malaria mitigation service forests provide relevant to this study; the study reports the morbidity-related malaria-reducing benefits of primary forests of, at least, $1–2 per hectare in Indonesia. Without the monetary value of the forest ecosystem services as input into decision-making, cost effective conservation and development policies may not be implemented well (Atkinson et al., Reference Atkinson, Bateman and Mourato2012).
Despite the fact that India has one of the highest burdens of malarial disease globally (World Health Organization, 2021), little empirical work has been conducted to determine the relationship between forests, deforestation, and malaria within the country. For example, using longitudinal data for 2000–2019, Ranjha and Sharma (Reference Ranjha and Sharma2021) find a positive correlation at the district level between forest cover and levels of malaria. Using pairwise comparisons and ${\chi ^2}$
tests, Saxena et al. (Reference Saxena, Nagpal, Singh, Srivastava, Dev, Sharma, Gupta, Tomar, Sharma and Gupta2014) find that deforestation between 2000 and 2009 led to the dispacement of one mosquito species with another, potentially increasing malaria rates in the Assam region. Using ${\chi ^2}$
and Student's t tests, Sharma et al. (Reference Sharma, Tyagi, Padhan, Upadhyay, Haque, Nanda, Joshi, Biswas, Adak, Das, Chauhan, Chitnis and Subbarao2006) find higher incidence of malaria inside forest villages relative to those on plains. While the studies noted above provide correlations, none of them use causal approaches or attempt to quantify the value of forests in preventing malaria in the country.
Combining socioeconomic data from the national Demographic and Health Survey (DHS) with spatially explicit remote sensing data, we quantify the benefits from forests in terms of reducing the malaria mortality in India between 2013 and 2015. We find that one additional square kilometer of forest lost increases the deaths from Plasmodium falciparum by 0.16 per 100,000 people. This translates into a marginal value of forests of, at least, $1.26–85.9/ha/year in 2015 US$. We make two contributions to the existing literature. First, we use rigorous econometric methods (e.g., heterogeneity robust difference-in-difference estimators) to causally identify the role of forests in mitigating the malarial burden in India. Second, we use non-market valuation methods, specifically estimating damage costs, to quantify the malaria-reduction benefits from forests to inform future policies. Our results are a first step in the design of effective conservation and development interventions to decrease the burden of malaria in the country.
2. Study area
Our study spans rural areas in 628 districts in India between 2013–2015 (figure 1). During the study period, India lost 335,900 hectares of primary forest, comprising about 1 per cent of its forest cover in 2010 (Global Forest Watch, n.d.). On average, 91 per cent of the annual forest cover loss during the study period is attributed to commercial forestry, followed by shifting agriculture (∼5 per cent), commodity driven deforestation (∼1.8 per cent), urbanization (0.9 per cent), and wildfires (0.13 per cent) (Global Forest Watch, n.d.).
Distribution of the forest area in 2013 (left panel), changes in malaria mortality due to P. falciparum (middle panel) for 2013–2015, and forest loss (in km2) at the district level for 2013–2015 (right panel). In the left panel, non-forest areas are given in white; in the right and middle panels white indicates no change in the forest cover. We focus on conterminous India and do not include data from islands
Malaria is still prevalent in the country: in 2021, it was estimated that 83 per cent of the cases and 82 per cent of the deaths in Southeast Asia were in India (World Health Organization, 2021). About 95 per cent of the country spans areas suitable for the transmission of malaria (Gething et al., Reference Gething, Patil, Smith, Guerra, Elyazar, Johnston, Tatem and Hay2011). Most of the deaths are attributed to P. falciparum; approximately one-half of all cases and one-third of the deaths in the country are attributed to P. vivax (World Health Organization, 2021).
The prevailing mosquito species carrying malaria include Anopheles baimaii, Anopheles culicifacies, Anopheles fluviatilis, Anopheles minimus, Anopheles stephensi, and Anopheles sundaicus; their distribution varies based on the ecosystem type (Subbarao et al., Reference Subbarao, Nanda, Rahi and Raghavendra2019). Of these, Anopheles (An.) culicifacies and An. fluviatilis jointly contribute to over 75–80 per cent of malaria cases in India; by itself, An. culicifacies accounts for about 60–70 per cent of the malaria cases in the country (Subbarao et al., Reference Subbarao, Nanda, Rahi and Raghavendra2019). The densities of An. stephensi, An. fluvialitis and An. culicifacies have been shown to increase following deforestation in India (Yasuoka and Levins, Reference Yasuoka and Levins2007). Because of the biology of the mosquito species, we expect that forest loss in India is likely to exacerbate the malaria incidence, ceteris paribus. Empirically, this finding is supported by the correlations between forest loss and malaria mortality from P. falciparum in figure 2.
Correlations between forest loss (left panel) and malaria mortality (right panel) at the district level. The values represent averages across all districts in our sample
The Indian government initiated the large-scale distribution of insecticide-treated mosquito nets (ITNs) in 2016, with the goal of eradicating malaria by 2027 (Narain and Nath, Reference Narain and Nath2018; Indian Ministry of Health and Family Welfare, 2020).Footnote 1 To avoid any confounding effects of this program, we limit the analysis to 2015. Indoor Residual Spraying (IRS) is also used to prevent malaria; while we do not have data on the prevalence of this approach, previous studies report that less than 20 per cent of the Indian households use ITNs or IRS (Wangdi et al., Reference Wangdi, Gatton, Kelly, Banwell, Dev and Clements2016). We assume that household behaviors do not change during the three years spanned by our analysis. That is, we assume the panel data estimators we use also control for the use of IRS.
3. Methods
3.1 Non-market valuation
We use a damage cost approach following Dickie (Reference Dickie, Champ, Boyle and Brown2017) and adopt the notation therein to estimate the value of malaria-reduction benefits from forests. Specifically, we model each person in our study area to consume two goods –a market good, x, and a health-related good, h. In our case, h is malaria morbidity and is the output of a (household) production function h = f(I, q), where I is a private good input such as an ITN and q-captures environmental quality such as forest loss. Thus, if the environmental quality decreases, I can increase to compensate for that; therefore, we refer to I as ‘averting behavior’. Assuming a quasi-concave utility, U = U(x, h), with V(p, q, y) as the corresponding indirect utility, and a budget constraint y = x + pI, where y represents the household income as a sum of expenditure on a numeraire consumption good x with price normalized to 1 and the averting expenditure, I at price p, we derive the following first-order conditions (FOC) (same as in Dickie, Reference Dickie, Champ, Boyle and Brown2017):
where MWTP indicates marginal willingness to pay, $\lambda$
is the Lagrange multiplier and equals the marginal utility of income at the maximum (Dickie, Reference Dickie, Champ, Boyle and Brown2017), and I* indicates the Marshallian demand for the averting good. The first term $(\partial U/\partial h)/\lambda \;$
captures the value of statistical life (VSL) and can be obtained from previous studies (Thaler and Rosen, Reference Thaler, Rosen and Terleckyj1976; Deschênes and Greenstone, Reference Deschênes and Greenstone2011).
The second term $dh/dq\;$
in equation (1) is the total effect of the change in malaria morbidity due to forest loss and can be estimated using our data. Because price data for insecticide nets are not available for our study area and the distribution of nets is oftentimes free in India (e.g., Raghavendra et al., Reference Raghavendra, Chourasia, Swain, Bhatt, Uragayala, Dutta and Kleinschmidt2017), we ignore the $p(\partial I\ast{/}\partial q)$
term. We therefore estimate a damage cost, which is a lower bound of the true benefit from forests as $\partial I\ast{/}\partial q$
is negative (Dickie, Reference Dickie, Champ, Boyle and Brown2017).
3.2 Estimation
To obtain an estimate for $dh/dq$
in equation (1), we assume a linear relationship between malaria morbidity, forest cover loss, and socioeconomic and biophysical covariates:
where h indicates malaria morbidity, c indicates a sampled DHS cluster, d is district, and t is year, with $t = 1$
being the first year of the treatment and T the end year; $t = 0$
indicates a baseline and ${\alpha _t}$
is year fixed effects. z represents time-varying cluster-level precipitation. All other relevant covariates like poverty levels and the presence of ITNs are time-invariant by assumption and are captured by the cluster-level fixed effects, ${\alpha _c}$
. Forest loss is captured by q, calculated at the district level. Because the impact of forest loss may vary through time, we estimate $\beta$
as a function of when an observation is treated and by how much. That is, $\beta$
captures the event-study effect of forest loss in each period after the baseline. Since a DHS cluster is smaller than a district, we treat the district-level forest cover as exogenous. ${\varepsilon _{cdt}}$
is an error term assumed to be independently and identically distributed (i.i.d.).
Forest loss occurs at different times in the different districts. Further, within a district, forest loss may occur multiple times and with different magnitudes. For these reasons, we apply a dynamic treatment heterogeneity-robust panel data estimator to recover the effect of continuous forest cover loss (de Chaisemartin and D'Haultfœuille, Reference de Chaisemartin and D'Haultfœuille2024). The estimator allows us to calculate the total effect on morbidity per unit forest loss, holding everything else constant: first, for each t, the approach estimates the expected difference in the outcome of the sample units treated at the same time to a counterfactual outcome of a group with the same baseline treatment that remains unchanged. Under the assumptions of parallel trends and no anticipation, the estimator tests for impacts of lagged treatments on the outcome and compares the impacts of current and lagged treatments. It then aggregates the impacts across time periods and groups comprised of observations treated at the same time and normalizes them by the number of times a sample unit is treated, to generate the average total effect per unit treatment. In contrast to the traditional two-way fixed effects (TWFE) model, the heterogeneity-robust estimator does not include observations that were treated in previous periods as part of the control group at time t (de Chaisemartin and D'Haultfœuille, Reference de Chaisemartin and D'Haultfœuille2020, Reference de Chaisemartin and D'Haultfœuille2024). Further, by using a dynamic treatment estimator, we avoid any issues related to arbitrary and, potentially, negative weights associated with traditional TWFE estimators when treatments take place at different times (de Chaisemartin and D'Haultfœuille, Reference de Chaisemartin and D'Haultfœuille2024). We perform the estimation using the did_multiplegt_dyn statistical package in Stata (de Chaisemartin et al., Reference de Chaisemartin, Ciccia, D'Haultfœuille, Knau, Malézieux and Sow2023).
Identification is determined by the number of clusters in districts that experienced forest cover loss during the study period. Because the estimator requires a finite number of values for the treatment variable, we discretize the extent of forest loss using 10 km2 increments. That is, we recode the forest loss variable to be equal to 1 for forest cover loss greater than 0 but smaller than 10 km2 in a given year, 2 for forest loss between 10 and 20 km2, 3 for forest loss between 20 and 30 km2, 4 for forest loss between 30 and 40 km2, 5 for forest loss between 40 and 50 km2, 6 for forest loss between 50 and 60 km2, 7 for forest loss between 60 and 70 km2, 8 for forest loss between 70 and 80 km2, and 9 for forest loss between 80 and 90 km2. Because of the few observations in the tail of the forest loss variable distribution, treatment for observations with more than 90 km2 forest loss in a given year is given a value of 10. Untreated observations are given a value of 0. For the main specification, all observations are given a value of 0 for the treatment (forest loss) at the baseline.
Since we do not have data on forest regrowth, we observe only forest cover loss. For this reason, forest cover is monotonically non-increasing over time; similarly, the lost forest area is monotonically non-decreasing over time: a cluster may have lower forest at $t = 1$
than at $t = 0$
(baseline) or it may have the same forest area. An observation with forest cover loss at $t = 1$
and no forest loss at $t = 2$
is still considered treated at $t = 2$
. If it lost forest at both $t = 1$
and $t = 2$
, we recorded the cumulative effect at $t = 2$
. For example, the West Godavari district in the state of Andhra Pradesh lost 1.69 km2 of forest between 2013–2014 and 1.43 km2 between 2014–2015, respectively. Thus, the treatment for clusters in that district is coded as 1 for 2014 and 1 + 1 = 2 in 2015.
The estimator allows for a very small number of exogenous time-varying covariates; we include precipitation as it is exogenous and relevant for mosquito abundance. Because the treatment is at the district level, we also cluster the standard errors at this level (Abadie et al., Reference Abadie, Athey, Imbens and Wooldridge2022).
3.2.1 Robustness checks
We perform a number of robustness checks: (a) alternative definitions of the treatment variable, (b) sensitivity of the results to initial conditions, (c) comparison with the traditional TWFE and first difference estimators, (d) contemporaneous heterogeneity robust estimators with a discretized multi-value forest loss area variable; (e) a panel data instrumental variable estimator with a control function approach to account for endogenous change in forest cover, and (f) a comparison with different ways the malaria mortality variable is calculated. Finally, we estimate the model using districts as the unit of analysis.
3.2.1.1. Treatment definitions
We repeat the main estimation using a binary treatment equal to 1 if there is any forest loss in 2014 or 2015 and 0 if no forest loss: we assume that once an observation is treated, it remains so. In additional specifications, we define the forest treatment variable as independent from previous years. That is, if an observation lost forest at $t = 2$
, we recorded the treatment at $t = 2$
. For example, the West Godavari district in the state of Andhra Pradesh lost 1.69 km2 of forest between 2013–2014 and 1.43 km2 between 2014–2015, respectively. Thus, in contrast to the main specification, here the treatment for clusters in that district is coded as 1 for 2014 and 1 in 2015.
Similarly, we repeat the main estimation by replacing the treatment levels by the means of each treatment bin in each year and adjusting for cumulative impacts. For example, instead of the treatment taking a value of 1 in 2014 if the forest loss is between 0 and 10 km2 as in the main estimation, we use the mean value of the forest loss (= 0.995 km2) for the treatment in 2014. If the observation lost forest between 0 and 10 km2 between 2014 and 2015, we used the mean forest loss for that bin in 2015 (= 0.88644046) and added that to the mean forest loss per bin from the previous year. That is, for this observation, the treatment takes a value of 0.995 in 2014 and 1.8814747 in 2015.
3.2.1.2. Sensitivity to initial conditions
The main specification has treatment = 0 for all observations in 2013. To address concerns that prior treatments could affect the outcomes, we repeat the dynamic estimation using the discretized forest loss between 2012 and 2013. We use the same bins as for the main analysis.
All districts lost some forest between 2012 and 2013. For this reason, we consider a district as ‘untreated’ at the baseline if it lost less than 10 km2 of forest between 2012 and 2013. In additional specifications, we split the sample by the value of the forest loss between 2012 and 2013: specifically, we estimate the dynamic model in the sub-samples the smallest (<10 km2) and largest forest loss (>90 km2) between 2012 and 2013.
3.2.1.3. First difference estimator
The main estimation relies on changes in the forest cover each year. If those are too small, the impact on malaria mortality may be insignificant. To address potential concerns that the small area of forest loss each year may render the impact on morbidity insignificant, we estimate a two-period first difference equation:
where $\Delta h$
indicates the change in malaria mortality in any given cluster c in district d from 2013 to 2015; t indicates a time period, with $t = 0$
being the baseline. Because the outcome variable is calculated as the malaria mortality in 2015 less the values in 2013, negative values indicate a decrease in malaria, whereas positive indicate an increase. $\alpha$
captures the temporal trend in malaria mortality; $\beta$
captures the effect of the district-level forest loss variable calculated as the difference in forest area between 2013 and 2015. Our hypothesis is that $\beta$
is positive. $\gamma$
captures the impacts of exogenous time-varying covariates $\Delta {z_{cdt}}$
. In our models, $\Delta {z_{cdt}}$
captures the change in precipitation. ${W_{cd0}}$
contains baseline characteristics associated with malaria incidence and mortality: malaria mortality, forest area, population, altitude, wealth, presence of nets and toilets, caste, and tribe affiliations. Because of the short time span of the study period, we assume the socioeconomic drivers of malaria remain unchanged; therefore, they cancel out in a first difference equation. Since the treatment variable, the area of forest cover, is at the district level, we cluster the standard errors at that level (Abadie et al., Reference Abadie, Athey, Imbens and Wooldridge2022).
We recognize that some of the covariates in (3) may be endogenous – for example, the presence of mosquito nets. We control for factors that may be driving the adoption of nets (e.g., baseline malaria mortality, elevation). As a robustness check, we also estimate the model without these potentially endogenous covariates.
3.2.1.4. Traditional two-way fixed effects model
We estimate
as a traditional TWFE model. As before, ${h_{cdt}}$
captures the malaria mortality in cluster c in district d at time t. ${\alpha _c}$
and ${\alpha _t}$
capture cluster and year fixed effects, respectively. ${z_{cdt}}$
contains the exogenous time-varying precipitation. The forest loss variable is given in ${q_{dt}}$
. As before, we hypothesize $\beta$
is positive.
Note that a traditional TWFE model estimation with a continuous treatment that takes place at different times is likely to produce arbitrary and potentially negative weights used in the aggregation of the per period effects and, therefore, may introduce bias that cannot be signed a priori (de Chaisemartin and D'Haultfœuille, Reference de Chaisemartin and D'Haultfœuille2020). For this reason, we interpret the estimates with caution and use this specification as a baseline comparison with our preferred model.
3.2.1.5. Contemporaneous heterogeneity robust estimators
We compare the results from the dynamic heterogeneity-robust estimator with an estimator that calculates the instantaneous treatment effects. Specifically, we use the did_multiplegt command in Stata 16 (de Chaisemartin et al., Reference de Chaisemartin, D'Haultfœuille and Guyonvarch2019). To run the latter model, we use a discretized forest loss variable as well as an estimation with the stable_treatment option and ‘0’ to indicate no change in the treatment, i.e., to define the non-switchers. We use 1,000 bootstrap replications for the standard errors.
3.2.1.6. Instrumental variables
Because the unit of observation is a circle with a radius of 5 km, it is unlikely a single cluster can affect the forest cover loss within a district. However, as a robustness check, we use annual data on particulate matter (pm 2.5) to instrument for forest loss. The data are available at 50 km resolution (Inness et al., Reference Inness, Ades, Agustí-Panareda, Barré, Benedictow, Blechschmidt, Dominguez, Engelen, Eskes, Flemming, Huijnen, Jones, Kipling, Massart, Parrington, Peuch, Razinger, Remy, Schulz and Suttie2019). Our reasoning is that particulate matter should be negatively correlated with forests but should not affect malaria mortality directly. The use of air quality to proxy for forest cover is similar to MacDonald and Mordecai (Reference MacDonald and Mordecai2019) who use aerosols in the dry season as an instrument.
We estimate a panel data instrumental variable model with a control function approach (Wooldridge, Reference Wooldridge2015) using cluster and year fixed effects as well as precipitation as an exogenous control. While the first stage of the approach of the control function is identical to that of an instrumental variable technique, the second stage includes a constructed error term. Because of the statistical significance of the latter, we use a wild t bootstrap and 1,000 replications to generate the standard errors. In all specifications, we cluster the standard errors at the district level (Abadie et al., Reference Abadie, Athey, Imbens and Wooldridge2022).
3.2.1.7. District-level analysis
Finally, we repeat the heterogeneity-robust estimators and the instrumental variable approach using data aggregated at the district level. The new dataset contains 583 observations; we excluded districts with less than 0.07 km2 forest in 2013 as well as those where malaria is not found.
4. Data
4.1 Unit of analysis
We use the spatial coverage of the 2015 wave of the DHS as a sampling frame (International Institute for Population Sciences – IIPS/India and ICF, 2017) the dataset provides the centroids of surveyed villages and distinguishes between rural and urban areas. We retain only the rural villages in our analysis. Most coordinates for the surveyed rural clusters are displaced by a random number within a 5-km radius; in addition, 1 per cent of those may be randomly displaced within additional 10 km (DHS, n.d.). However, even with the displacement, all clusters remain in the original district in which they were located (DHS, n.d.) We therefore consider a 5-km buffer around the DHS centroid as our unit of analysis; we refer to these as ‘clusters’.
We focus on the rural clusters with non-missing coordinates (n = 19,920), which constitute about 70 per cent of all DHS clusters. Because of the often-significant overlap between the 5-km buffers, we randomly selected 10 per cent of the sample (n = 2,053). In choosing a 10 per cent sample rather than all clusters, we try to minimize the overlap between clusters while simultaneously retaining a large enough sample for the statistical analysis. The average distance between cluster centroids in the sample is about 27 km (range 15–163 km). This means that the average distance between the 5-km buffers around each cluster point is 17 km. Because of the distance between clusters, we assume away any spatial dependence except for clusters in the same district having the same treatment. Without sampling, the average distance between the DHS rural cluster points is 5,519 m (range: 31–53,973.25 m).
We select the sample by first creating a fishnet of 200 rows by 200 columns over the extent of the country, with each cell a square with a side of roughly about 14 km. We randomly selected a cell (id = 22,256) and, using that as a starting point, dropped all immediately adjacent fishnet cells. This resulted in 9,404 fishnet cells being retained. Using these, we randomly selected 1 DHS point from each grid cell. A map of the retained DHS clusters and their 5-km buffers is available in the online appendix (figure A1).
We drop the clusters within districts with no forest as well as those in districts with less than 0.07 km2 of forest (n = 29). We also exclude clusters for which no malaria caused by P. falciparum is possible using the P. falciparum spatial limits from the Malaria Atlas (‘Plasmodium falciparum Spatial Limits for 2010’).Footnote 2 The random offsets of the sampled DHS clusters ensure that all clusters fall within the district where the original surveyed villages are located. Therefore, in case of mismatches between the geospatial district boundary data and the DHS survey, we retained the DHS district designation. The final dataset consists of a panel of 1,985 clusters for each of the three years (2013–2015).
Because of the different algorithm used to generate forest cover loss data after 2015 (Weisse and Potapov, Reference Weisse and Potapov2021) as well as the increased distribution of mosquito nets due to government efforts after 2016, we limit the analysis to 2013–2015 only. The short time frame also lends support to the assumption that behaviors with regards to IRS and ITNs remained unchanged and are controlled for by time fixed effects.
4.2 Outcome variable
Because data on malaria were not available in the India DHS data, we obtain the annual number of deaths per 100,000 for Plasmodium falciparum at 5 km resolution from the Malaria Atlas. The data are spatially variable and available for the whole country; because of the common methodology employed in creating the dataset, it is possible to compare different regions of the country. We calculate the total number of deaths per 100,000 people at the cluster level. Mortality data for P. vivax were not available.
We compare different measures of P. falciparum mortality from the Malaria Atlas and the World Health Organization (see online appendix, table A1a). We find that the Malaria Atlas data overestimate the malaria mortality by at least two times (see online appendix, table A1b). However, that the mean mortality increases monotonically between 2013 and 2015 is consistent across data sources (tables A1a and b). Because it is most conservative, we use median mortality layer from the Malaria Atlas. As a robustness check, we repeat the estimation using the lower and upper confidence interval values for malaria mortality as well as the modeled values for P. falciparum malaria: for the latter, we calculate the mean, max, and minimum values per cluster.
4.3 Forest cover
We use the Hansen data available at a 30 m resolution (Hansen et al., Reference Hansen, Potapov, Moore, Hancher, Turubanova and Tyukavina2013) to quantify changes in the forest area. We chose this dataset because it spans the whole country and provides annual changes in forest cover. It has known limitations: for example, the forest loss data we use are conservative and less likely to detect impacts from tree logging prior to 2016 (Weisse and Patapov, Reference Weisse and Potapov2021). Further, the Hansen dataset does not differentiate between the type of forest being lost (e.g., planation vs. native forest, moist broadleaf vs. dry broadleaf, etc.). However, previous studies in India have suggested most of the forest cover in 2000 was comprised of natural forests, with only 0.11 per cent being tree plantations (see ‘Location of forest in India’ in Global Forest Watch (n.d.)). The Hansen dataset also does not allow us to distinguish between temporary forest loss and permanent forest conversion. Finally, our data exclude reforestation as the available reforestation data from the Hansen et al. (Reference Hansen, Potapov, Moore, Hancher, Turubanova and Tyukavina2013) dataset are only cumulative, spanning a period of 2000 to 2012 only. For this reason, our analysis excludes forest regrowth as well as any reforestation and afforestation efforts and instead focuses on the loss of existing forest.
The Hansen data provide the annual percentage tree cover in 2000. Using that layer as the baseline, we create a binary forest/no forest layer using a 25 per cent tree cover per pixel as a cutoff to define forests (Sexton et al., Reference Sexton, Noojipady, Song, Feng, Song, Kim, Anand, Huang, Channan, Pimm and Townshend2015); we use the binary forest-no forest layer for 2000 to filter any forest loss events taking place on cells that were not forested initially. Note that ‘forest loss’ in our data refers to any event that results in a pixel losing tree cover, so that tree cover falls under 50 per cent of the pixel area in a given year.
To obtain the forest cover for each year, we subtract forest loss events from the baseline binary forest layer. Then, using district boundaries for 2015, we calculate the total area of forest annually for each district. We exclude districts with less than 0.07 km2 of forest in 2013.
4.4 Biophysical covariates
To obtain data on precipitation, we use the Climate Hazards group Infrared Precipitation with Stations (CHIRPS) dataset (Funk et al., Reference Funk, Peterson, Landsfeld, Pedreros, Verdin, Rowland, Romero, Husak, Michaelsen and Verdin2014), available at 5 km resolution. We calculate the mean precipitation within a DHS cluster. We use the altitude data from the DHS surveys for each cluster.
4.5 Socioeconomic covariates
We include socioeconomic covariates in the first difference models as a robustness check only. Assuming that the socioeconomic characteristics do not change within the short time period of the analysis, we use the 2015 wave of the DHS as it is the first that is geo-located; these data give us a proxy for the baseline socioeconomic characteristics in the area. The data contain spatially explicit information on 1,315,617 individuals (601,509 households), the majority of whom (∼75 per cent) are located in rural areas. The survey uses responses from men aged 15–54 and women aged 15–49 and provides detailed information on households and individuals. At the cluster level, we calculated the average share of the households with nets, the average shares of households considered to be poor and rich, the average share of households with a toilet, and the average share of households that belong to a tribe or a caste. In the calculation of these covariates, we ignored any survey weights. To obtain the population density within a cluster, we use data from Landscan, available annually at 1 km resolution (Bright et al., Reference Bright, Rose and Urban2016). We calculate the total number of people within a cluster. Descriptive statistics of the variables used in the estimation are given in table 1.
Descriptive statistics for the cluster-level sample
Notes: ‘Cluster-level’ pertains to a 5-km buffer around the DHS 2015 geospatial points. ‘HH’ indicates households. For precipitation and population density, negative values indicate an increase between 2013 and 2015; for malaria mortality, positive values indicate an increase between 2013 and 2015.
5. Results
5.1 Impact of forests on malaria mortality
The heterogeneity-adjusted panel data estimators indicate that a reduction in the forest area resulted in more deaths from malaria (table 2). Using the main model estimates in column [1] and dividing by 10 due to the treatment definition, we find that on average, the loss of 1 km2 of forest resulted in 0.16 additional deaths per 100,000 people. The estimate is qualitatively the same as in the dynamic model where the treatment reflects changes in forest cover for each period (column [2]) rather than cumulative changes.
Results from the panel data estimators for the 2013–2015 data
Notes: The outcome variable is the number of deaths per 100,000 people. The standard errors (given in parentheses) are clustered at the district level. The reported sample sizes reflect the number of observations per period multiplied by the number of time periods. Columns 1–4 contain the results from the heterogeneity robust estimators. The TWFE results are in column [5]. The results from an instrumented variable regression, using particulate matter as an instrument, are given in column [6]. Because the control function approach yields the same coefficients but slightly smaller standard errors, we do not present the estimates here. We consider [1] to be the main model.
The sign and statistical significance of the estimate is consistent across the heterogeneity-robust specifications (table 2) as well as the traditional models (tables A6, A8 and A9, online appendix). However, when we account for the impact of treatment lags, the magnitude is larger: in the models with continuous treatment, the total average treatment effect is about twice as large as the models that calculate the instantaneous effect (columns [3]–[5]). The event-study per period effects for the main estimator are plotted in figure 3: the magnitude of the average total impact is driven by the impact at $t = 1$
; the event study effect for $t = 2$
(coefficient = 0.83, st. error = 0.15) is consistent with the heterogeneity-robust instantaneous effects in table 2. After accounting for differences in the treatment levels in table 2, the average treatment effect from the main model is smaller than the coefficient from the instrumental variable two-stage estimation (column [7]) and the model using mean forest area loss levels to define the treatment (column [8]) in table 2.
Graph of the results from the dynamic panel data estimation. The horizontal axis shows the period when the treatment changes, with −1 indicating a placebo
The results are consistent and robust across different definitions of the outcome variable and treatment but vary in magnitude. The average treatment effect is about twice as small as in the ones from the models with binary treatment (table A4, online appendix). Except for the model using the average lower bound of the malaria morbidity layer per cluster, the main estimate is the most conservative (online appendix table A10). The results from the cluster-level sample are very consistent in magnitude with the estimation using the district-level sample (online appendix table A11).
Repeating the dynamic estimation with the number of people in a cluster as the outcome, we find an average effect of −140.619 (st. error = 46.38,262) in clusters with forest loss. The negative impact suggests that the increases in the malaria mortality in the districts with forest loss are not driven by an influx of people.
The analysis on the role of the baseline forest loss indicates positive coefficients that are consistent with the main specification (tables A3a and 3b in the online appendix). However, in some specifications, the number of observations drops significantly, raising concerns about statistical power (de Chaisemartin and D'Haultfœuille, Reference de Chaisemartin and D'Haultfœuille2024). The coefficient on forest loss is statistically significant for the sub-sample with large values (>90 km2) of forest loss between 2012 and 2013 (table A3b). The robustness check with districts as the unit of analysis and forest loss less than 10 km2 between 2012 and 2013 results in consistent and statistically significant estimates comparable with the main estimation (table A11, column [2]). For these reasons, we are not concerned about initial conditions potentially biasing our results significantly.
5.2 Value of forests in reducing malaria mortality
We use the estimates from table 2, column [1] to calculate the MWTP for avoided forest loss per equation (1). The reported VSL estimates for India vary significantly. For example, Majumder and Madheswaran (Reference Majumder and Madheswaran2016) report VSL for avoided mortality ranging from $153,000–358,000 (in 1999 US$) taken from a study by Simon et al. (Reference Simon, Cropper, Alberini and Arora1999), to $3.74 million (in 1990 US$) taken from a study by Shanmugam and Madheswaran (Reference Shanmugam, Madheswaran, Haque, Enamul, Murty and Shyamsundar2011). Using data on India's Consumer Price Index data (World Bank, n.d.) to ‘age’ the estimates as in Das and Vincent (Reference Das and Vincent2009), these values translate to 146,755–343,388 to 6.97 million, respectively in 2015 US$. Sweis (Reference Sweis2022) report a VSL for India of 0.2–1.2 million US$ in 2019. After using the same formula from Sweis (Reference Sweis2022) and the gross national product for India for 2015, we find a range of 0.28–2.80 million US$ for 2015. The lowest value is from Ozawa et al. (Reference Ozawa, Stack, Bishai, Mirelman, Friberg, Niessen, Walker and Levine2011) who report a VSL of $41,100 per capita for malaria in 2015. Because of the wide dispersion of the VSL values, we use the range US$41,100–2.8 million for the VSL in 2015. Plugging into equation (1) the VSL values and the estimated total impact of forest loss on malaria morbidity from table 2, we obtain a WTP for each avoided square kilometer of forest loss of between $0.07–4.47/km2 in 2015 US$.
To aggregate and convert to average values per hectare, we multiply the estimates by the rural population in a district in 2015 and divide by the forest area, converting the square kilometers to hectares. Because the DHS survey weights are aggregated and adjusted for non-responses, making it difficult to generate the population density per district, we use the Landscan population density data for 2015 (Bright et al., Reference Bright, Rose and Urban2016) to obtain the population within rural areas. For the latter, we use the GRUMP dataset (CIESIN et al., 2011) to define the rural zones within districts. We then calculate the total number of people in the rural parts of districts with forest cover in 2015 using the zonal statistics option in Arc Map 10.1. We find a total of 678,000,000 people living in rural areas in 2015. We rescale the estimates based on the ratio of World Health Organization-to-Malaria Atlas data for 2015 (table A1a). The approach results in annual average marginal forest value of $1.26–85.9/ha/year in 2015 US$.
6. Discussion
Using panel data for 2013–2015, we find robust evidence that forest loss increased malaria mortality in rural India. The annual value of forests in avoiding malaria mortality is between $1.26 and $85.9/ha/year in 2015 US$. Our estimates are a lower bound for two main reasons. First, we exclude any defensive expenditures and other behaviors households may undertake to minimize the risk of malaria. Second, we exclude the impact of forests on malaria morbidity. For example, using data on foregone wages in Indonesia, Garg (Reference Garg2019) finds that forests provide $1–2/ha of value in reducing malaria morbidity.
Our results are subject to three caveats. First, we focus on a relatively short time span. We limit our analysis to 2013–2015 because a nationwide government program aimed at distributing insecticide-treated nets was introduced in 2016 (Indian Ministry of Health and Family Welfare, 2020; World Health Organization, 2021). The wide distribution of nets raises econometric concerns related to endogeneity – understanding which households receive a net, when, and who within the household receiving nets uses them and how much a net costs; because of data limitations we cannot model for the household decision variables driving the use of nets and cannot address the endogeneity adequately if we were to use data after 2016. Without sufficient data on behaviors and the prices of averting expenditures, we may not be able to separate the effects of forests from those of mosquito nets. The data on malaria incidence and mortality indicate a rapid decrease prior to COVID (World Health Organization, 2021), suggesting the program may be effective. However, concerns have been raised about the potential of increasing insecticide resistance for ITNs (e.g., Faizi and Kaur, Reference Faizi and Kaur2021; World Health Organization, 2021) and IRS (Sahu et al., Reference Sahu, Tediosi, Noor, Aponte and Fink2020).
Second, our estimates capture the short-term impact of forest loss on malaria mortality. The conversion of forests to urban areas may alter exposures and habitats for malaria; for this reason, the longer term impact of forest loss on morbidity and mortality is unclear (MacDonald and Mordecai, Reference MacDonald and Mordecai2019).
Third, our work does not differentiate between the different types of forests. An emerging body of work has suggested heterogeneous impacts of forests by distinguishing between old-growth versus disturbed (Pattanayak and Yasuoka, Reference Pattanayak, Yasuoka and Colfer2008), core versus edge and perforated (Blackman and Leguízamo, Reference Blackman and Leguízamo2023; Cheng and Miteva, Reference Cheng and Miteva2024) and primary versus secondary forests (Garg, Reference Garg2019). Owing to data limitations, we do not examine heterogeneity in the impacts of forests based on whether they are secondary versus primary forests, native versus plantation forests, or old-growth versus new forests; exploring the differences between core and disturbed areas is beyond the scope of the current analysis. Previous studies have suggested that native forests and plantations may differ in their role in malaria transmission, with plantations offering more opportunities for prolonged exposure to mosquitoes and, therefore, higher morbidity (Kar et al., Reference Kar, Kumar, Singh, Carlton and Nanda2014). If the latter holds, we expect our estimates to be downward biased. However, based on the data from Du et al. (Reference Du, Yu, Yang, Xu, Chen, Peng, Zhang, Fu, Harris and Gong2022), most of the forest areas in India are either natural or planted (sal and alders), with only a small area on the east coast occupied by cashew plantations. For this reason, we expect the bias to be small.
Globally, multiple interventions like ITNs, IRS, rapid malaria diagnostic and treatment, and vaccines against malaria are currently being implemented to reduce the burden of the disease (World Health Organization, 2021). However, these interventions are implemented in isolation without any landscape – and, specifically, forest conservation – considerations, which could complement the traditional approaches. The latter rely on funding and the goodwill of multiple global and local partners and may be disrupted, for example, because of the COVID pandemic (World Health Organization, 2021). We show that forest conservation can complement existing efforts to reduce the cost of malaria mitigation and generate additional benefits like climate change mitigation, support for biodiversity, and the provision of fuelwood. Our work can be used as a first step in the design of suitable interventions, like payments for ecosystem services, to incentivize local actors to protect forests.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S1355770X25000075.
Competing interest
The authors declare none.
Open access
