Improving rural health care reduces illegal logging and conserves carbon in a tropical forest.

Tropical forest loss currently exceeds forest gain, leading to a net greenhouse gas emission that exacerbates global climate change. This has sparked scientific debate on how to achieve natural climate solutions. Central to this debate is whether sustainably managing forests and protected areas will deliver global climate mitigation benefits, while ensuring local peoples' health and well-being. Here, we evaluate the 10-y impact of a human-centered solution to achieve natural climate mitigation through reductions in illegal logging in rural Borneo: an intervention aimed at expanding health care access and use for communities living near a national park, with clinic discounts offsetting costs historically met through illegal logging. Conservation, education, and alternative livelihood programs were also offered. We hypothesized that this would lead to improved health and well-being, while also alleviating illegal logging activity within the protected forest. We estimated that 27.4 km2 of deforestation was averted in the national park over a decade (∼70% reduction in deforestation compared to a synthetic control, permuted P = 0.038). Concurrently, the intervention provided health care access to more than 28,400 unique patients, with clinic usage and patient visitation frequency highest in communities participating in the intervention. Finally, we observed a dose-response in forest change rate to intervention engagement (person-contacts with intervention activities) across communities bordering the park: The greatest logging reductions were adjacent to the most highly engaged villages. Results suggest that this community-derived solution simultaneously improved health care access for local and indigenous communities and sustainably conserved carbon stocks in a protected tropical forest.

Tropical forests lose more than 100 trees every second, altering landscapes and impacting livelihoods, health, biodiversity, and climate change (1). Across the tropics, forest loss now exceeds forest gain, leading to a net carbon emission from some of the most important natural carbon stocks in the world (2). Averting further forest loss is an important natural climate solution and a high priority for science, management, and policy from local to global scales (3, 4).

In biodiverse, carbon-rich tropical forests, the establishment of protected areas benefits both conservation and climate mitigation goals, but often involves excluding, and thus disenfranchising, local communities that surround protected areas (5, 6). Failure to address the needs of local people can in turn lead to unsustainable forest use, when communities with few alternatives illegally extract resources and convert land in order to survive (6, 7). Another major hypothesized driver of poverty is lack of access to high-quality, affordable health care, which can lead to vicious cycles of poor health and expanding out-of-pocket costs, further incentivizing poor families to rely on unsustainable resource use, like illegal logging, in order to raise cash to meet critical health care needs (8). Within this context, this study examines whether providing rural health care incentivizes reductions in illegal logging by local and indigenous communities living around a national park in Indonesian Borneo (Fig. 1), thereby improving health and well-being and conserving biodiversity and globally important carbon storage.

Fig. 1.

Cross-sector global health and forest conservation needs. (A) Maps of global aboveground forest carbon density and Universal Health Coverage (9): Tropical areas, particularly Africa and Asia, have high forest cover and low health care coverage. (B, Inset) Forest loss (resulting from deforestation and forest degradation) accelerates over time across all 32 terrestrial IUCN Category II National Parks (10) established before 2001 in Indonesia [boxplots; forest change data: Hansen et al. (1)]. (C) Study site and approach: locations of IUCN Category II National Parks (10) in Indonesia, with the intervention park highlighted, and an outline of the problems and hypotheses addressed in this analysis, along with hypothesized outcomes that were tested empirically through objective earth observation and health clinic records.

Indonesia contains some of the most carbon-dense forests in the world () (11), with the island nation representing only 1.4% of the world’s land area, but 3.6% of natural forest cover (12). A 2011 moratorium on new logging concessions was implemented to reduce total emissions from deforestation in Indonesia (12), but, at the same time, illegal logging was estimated to represent as much as 61% of all logging activity (13). Protected areas, which cover 12% of Indonesia’s land mass (10), may be key to preserving remaining carbon stocks, biodiversity, and forest ecosystem services (14). However, top-down law enforcement has proven insufficient to prevent conservation threats to protected areas, like illegal encroachment and fire (15). Indeed, despite protected status, more than 60% of lowland forests within protected areas in Borneo’s West Kalimantan region were lost to illegal logging between 1985 and 2001 (16, 17), and this trend continues to accelerate across the region (Fig. 1).

To better understand local attitudes toward forests, conservation, and drivers of illegal logging in West Kalimantan, Indonesia, a nonprofit organization conducted more than 400 h of focus groups between 2005 and 2007 with nearly 500 community representatives (community leaders, farm group leaders, religious leaders, teachers, women’s group leaders, and interested community members). The open-ended conversations identified access to affordable, high-quality health care as a major basic need, and lack of access to health care a potential driver of illegal logging in 23 districts near Gunung Palung National Park (GPNP) (Fig. 1) (18). This corresponds to a broader concern in Indonesia, which ranks in the lowest one-third of countries in terms of coverage of essential health services (Fig. 1) (9), and where many populations contend with unmet water, sanitation, and hygiene needs, high maternal and infant mortality, and high burdens of infectious and noninfectious diseases (19). In response to local health care needs and conservation implications, the nonprofit established a local health clinic in 2007, in close partnership with the district government and the national park management. Clinic services and alternative payment options (e.g., barter options including seedlings and manure used in conservation activities) were available to anyone seeking care, complementing the limited health care available provided by the government. With the support of the National Park management, memorandum of understanding (MOU) agreements were signed by the nonprofit and 21 of 23 districts (“desa” administrative units), representing 73 villages (“dusun” administrative units), near GPNP to participate in the health care–conservation exchange intervention. Through the intervention, clinic discounts were given to villages based on community-wide reductions in illegal logging activity, as reported by community liaisons and through monitoring of logging trail activity. At the request of community members in intervention areas, conservation programs, educational programs, and alternative livelihood trainings were also facilitated periodically in partnership with government entities. In 2007, 2012, and 2017, household surveys were conducted in random households across villages that engaged with the intervention to assess self-reported changes in well-being, knowledge, attitudes, behaviors, and livelihoods, in order to adapt intervention activities accordingly to meet community needs.

Here, we use more than 10 y of de-identified patient records from the health clinic coupled with remotely sensed earth observation data to test the hypothesis that a multisector health care–conservation intervention can simultaneously improve human health care access and use, and disincentivize illegal deforestation in a carbon-rich tropical forest. First, we use a synthetic controls approach to compare park-level forest loss rates in GPNP before vs. after the intervention began in 2007, compared to all terrestrial Indonesian International Union for Conservation of Nature (IUCN) Category II National Parks as potential controls. MOU-signing villages bordered the national park and were thus nonrandomly assigned, precluding analysis as a randomized controlled trial. However, we were able to use clinic patient records collected between 2008 and 2018 to compare health care access, usage, and diagnosis trends between MOU-signing and non–MOU-signing patient groups in a quasi-experimental study design. We tested whether, near GPNP, 1) the clinic increased health care usage for patients from villages with signed MOUs compared to villages without signed MOUs, and 2) trends in disease diagnoses changed over the same time period for patients from villages with vs. without signed MOUs. To establish the plausibility of a causal relationship between the intervention and conservation outcomes, forest change data, clinic usage data, and records on village-level engagement with the intervention programs (i.e., health clinic use and periodic education and livelihood programs) were then used to test for a dose–response of village-level forest loss within the national park to village-level intensity of engagement with the intervention and its associated programs. Last, select responses from self-reported household survey data collected by the intervention team in 2007, 2012, and 2017 were used to assess changes in household livelihoods and income over the intervention period, and to gain further insight into potential mechanisms driving conservation–health linkages.

Results

Intervention Impact on Forest Change.

Remotely sensed forest loss rates (1) in a synthetic controls analysis (20, 21) were significantly lower within districts intersecting the focal national park, GPNP, compared to a synthetic control assembled from districts in 32 control parks across Indonesia during the postintervention period from 2008 to 2018, compared to a preintervention period from 2001 to 2007 (estimate, 69.8% reduction of forest loss; 90% CI, 50.8–81.4; P = 0.003; Table 1 and Fig. 2). This translates to an estimated 27.4 km2 of forest loss averted postintervention to 2018 (90% CI, 19.9–32.0 km2) (Table 1 and Fig. 2). This finding was robust to a number of alternative data subsets defining the “donor pool” of possible control parks to be included in the synthetic control (Table 1), as well as to comparison of the intervention effect with 500 permuted “placebo” treatment regions made up of randomly assigned units from the pool of control districts (estimate, 69.8% reduction; permuted P = 0.038; 90% CI, 26.3–83.7; Fig. 2) (21). Correspondingly, household responses collected through monitoring and evaluation activities during 2007 (baseline), 2012, and 2017 demonstrated a strong and significant reduction in self-reported illegal logging activity: The number of adult males who reported logging inside GPNP during the intervention period compared to baseline declined (generalized linear regression with logit link; estimate, 68.8% reduction; 95% CI, 60.8–75.5; P < 0.001; Fig. 2), as did the number of households reported to rely on logging as a primary income source (generalized linear regression with logit link; estimate, 90.6% reduction; 95% CI, 83.4–95.2; P < 0.001; Fig. 3).

Table 1.

Results from the synthetic controls analyses on park-level forest loss in GPNP compared to a counterfactual derived from three subsets of Indonesian IUCN Category II National Park controls: All nonmarine parks established prior to 2001, all nonmarine parks, and all parks

Model	Forest loss, treated, km2	Forest loss, control, km2	% Change	P value [90% CI]	Permuted P value [90% CI]	No. obs.	No. district	No. parks
Nonmarine parks, est. before 2001	11.891	39.30	−69.75	0.003 [−81.4, −50.8]	0.038 [−83.7, −26.3]	27,702	1,539	32
Nonmarine parks	11.891	28.36	−58.1	0.013 [−74.0–32.4]	0.062 [−78.3, −1.6]	36,738	2,041	44
All parks	11.891	28.36	−58.1	0.013 [−74.0, −32.4]	0.080 [−80.6, 0]	40,320	2,240	52

The first two columns provide estimates of forest loss (in square kilometers) in the treated region following the intervention and loss in the synthetic control region. P values and confidence intervals are calculated from a standard normal sampling distribution and Taylor series linearization. A permuted P value and CI were calculated using 500 permuted “placebo” treatment groups to satisfy a more robust set of assumptions and generate a more conservative estimate of the sampling distribution [Robbins et al. (21)]. In both cases, the CIs do not contain 0, and based on a lower-tailed, one-sided hypothesis test, the null hypothesis that there is no intervention effect is rejected [Robbins et al. (21)].

Fig. 2.

Climate impacts. (A) As the number of adult men (age 19 or older) who report logging as a primary livelihood declined between 2007 and 2017 based on survey data (Inset), a synthetic controls analysis using remotely sensed earth observation data on forest change verified that forest loss rates in Gunung Palung National Park (green solid line, GPNP, the intervention park) were significantly lower than a synthetic control from which the counterfactual (black dashed line) for forest loss after the onset and ramp up of the intervention in GPNP (red vertical dashed line) was estimated. (B) The difference between forest loss in GPNP (green) and the 500 “placebo” synthetic control treatments (gray) made up of random permutations of the sampling units—the dotted black line on the x axis represents no difference between treatment and placebo groups. (C) Forest loss rates were converted to estimates of aboveground carbon biomass preserved, using average tree height calculated from a publicly available LiDAR-derived dataset and locally calibrated wood density equations (see for details). (D) Quantitative outcomes showing a dose–response in forest loss rates in GPNP to intervention effort: Engagement (see for details) was binned into low-, medium-, and high-engagement categories based on person-contacts across many intervention activities in each village, for 36 villages bordering GPNP with signed MOUs; changes in average forest loss rates (±1 SE) from the 5-y interval before the intervention (2002 to 2006) to the last 5 y of the intervention (2013 to 2017). •P < 0.10; ***P < 0.001; N.S., not significant.

Fig. 3.

Health impacts. (A) Individual visitation frequency (Left, average visits/patient to the health clinic during the study period) and health care use (Right, the percentage of the district population that were recorded at least once during the study period as patients at the clinic), among patients from districts that signed an MOU and thus received discounts on care, and those that did not; partial responses to MOU status are shown after controlling for distance effects (travel time to the clinic). (B) Change in odds of disease diagnoses from clinic patient records (presented as odds ratios for MOU and non-MOU patient populations [controlling for distance effects], comparing odds of diagnosis in 2008 to 2009 vs. 2017 to 2018 with 95% CIs; ). (C) Change in primary livelihoods including self-reported logging (proportion of households, 95% CIs) from 2007 to 2017. (D) Change in reported perceptions of neighborhood wealth (Left, where most responses are “average” in medium pink, versus "poor" in light pink, and "wealthy" in dark pink) and mean purchasing power parity (PPP)-adjusted household monthly incomes (Right), as reported from household surveys at 5- and 10-y follow-up periods (2012 vs. 2017). N.S., not significant; ***P < 0.001.

Intervention Impact on Forest Carbon.

Using published carbon equations parameterized specifically for Borneo (22) along with LiDAR (light detection and ranging)-estimated canopy heights in the focal national park, GPNP, in 2014 (Fig. 2 and ), the effect size of 69.8% reduction in annual forest loss was estimated to equate to a cumulative 0.59 Tg of aboveground carbon loss averted (90% CI, 0.27–1.13 Tg). Based on the maximum trade value of $30 per ton of CO2 realized on the European Emissions Trading System (23), the gross value of the total carbon loss averted in GPNP on the European carbon market would have been approximately $65.3 million USD in 2019. The estimate of aboveground carbon loss averted is conservative, because 1) the LiDAR flight data in GPNP covered mixed and some previously burned forest types, and the derived average vegetation height (27 m) is much lower than the tallest canopy height recorded in GPNP (71 m) (Fig. 2), and 2) the relationship between canopy height and aboveground carbon is convex and nonlinear, suggesting that averaging across 30 × 30-m pixels consistently underestimates the true carbon density (22). The carbon value of the intervention impact demonstrates a theoretical monetary return that would more than offset intervention costs if a carbon market were accessible to interventions aiming to couple rural health programs with forest conservation in a similar way.

Intervention Impact on Health Clinic Usage and Diagnoses.

Overall, 28,462 unique patients visited the clinic at least once over the study period from 2008 to 2018. Most patients came from districts located on the periphery of GPNP that signed MOUs to participate in the intervention, but a substantial fraction of patients (42%) came from districts without MOUs, who sometimes traveled many hours or days to use clinic services (). Clinic affordability (MOU status and associated discounts) and accessibility (estimated travel time to the clinic) jointly influenced two metrics of clinic usage: probability of clinic use, measured as the proportion of a district’s population that used the clinic at least once, and individual patient visitation frequency. Patients with shorter travel times to the clinic were more likely to use the clinic (Poisson generalized linear mixed model [GLMM] with population size offset and district random effect; estimate= −1.14; SE = 0.17; P < 0.0001) and visited the clinic more frequently (negative binomial GLMM with district random effect; estimate = −0.180; SE = 0.035; P < 0.0001) (). At the same time, controlling for distance, signing of an MOU (and receiving clinic discounts) increased clinic use: A larger proportion of MOU-signing district populations used the clinic (on average 27.8% vs. 2.76%; estimate = 1.93; SE = 0.36; P < 0.0001; Fig. 3 and ), and individual patients from MOU-signing districts visited the clinic 33% more often, on average (2.4 visits over 10 y vs. 1.8 visits; negative binomial GLMM with district random effect; estimate = 0.284; SE = 0.073; P < 0.0001; Fig. 3 and ). Patients that visited more than two to three times were usually returning repeatedly for health care related to a chronic health condition, such as epilepsy, emphysema, or hypertension. Overall, the clinic usage statistics confirm that, controlling for distance effects on clinic usage, signing an MOU to participate in the intervention incentivized increased use of health care services at the clinic. Even so, patients without MOUs represented more than 40% of all patient visits (), likely because noncash payment options like exchange of tree seedlings, manure, handicrafts, or labor made service affordable.

Time Trends in Disease Outcomes Based on Diagnoses at the Clinic.

De-identified diagnosis records from more than 61,000 unique doctor visits recorded during 2008 to 2018 showed improvements in many health outcomes for MOU and non-MOU patient populations. We found significant declines over time in diagnosed cases of malaria, tuberculosis, childhood-cluster diseases, neglected tropical diseases (NTDs), chronic obstructive pulmonary disease (COPD), and diabetes (Fig. 3). The only preventable and treatable diseases considered here that increased over time were lower and upper respiratory infections (Fig. 3 and ). The increase in diagnosed cases of respiratory diseases regionally may have been related to region-wide fire activity that spiked in 2015 (24). Increases in upper respiratory infections might also be due to increased care-seeking for more minor illnesses as trust was built between communities and the program. Time trends in diagnosed cases of disease were consistent whether or not district-level distance to the clinic was included as a covariate.

Regional diagnosis records from government clinics were not available to use as controls against which to compare temporal trends at the intervention clinic, but trends in diagnosis of several diseases departed from population-based prevalence estimates published by independent Global Burden of Disease studies for West Kalimantan during the same time period. For example, tuberculosis and COPD showed an upward trend regionally (25), whereas tuberculosis declined strongly in the clinic population in our study from 2008 to 2018, when the intervention’s health clinic oversaw TB-DOTS treatment (i.e., directly observed treatment short course for tuberculosis) for all regionally diagnosed patients (including in local government facilities).

Comparing Disease Outcomes in Patients from MOU-Signing vs. Non-MOU Districts.

We were unable to assess whether time trends in patient diagnoses were attributable to the increased health care access and use available through the clinic (beyond care available to all individuals through government-supported clinics in the region), because it would have been unethical to withhold access to any patient and comparable records for patients at other health facilities were not available to us for comparison over the same time period. As a result, we lacked a matched control (or “no clinic access”) group, but were able to statistically compare time trends in diagnoses among patients from MOU-signing districts vs. patients from districts without MOUs, to test whether community health outcomes benefited from clinic discounts associated with the intervention. Controlling for the distance between the clinic and patients’ home districts, we found few differences, indicating largely equitable health outcomes in terms of change in the proportion of patient diagnoses across all diseases (). The few exceptions included cases of lower respiratory infections (LRIs) and upper respiratory infections (URIs), which increased across all patient populations over the 10 y study period, but increased significantly less in MOU-signing patient populations (LRI estimate = −0.499; SE = 0.222; P = 0.0025; URI estimate = −0.650; SE = 0.222; P = 0.0036), as did cases of dental diseases (estimate = −0.877; SE = 0.167; P < 0.0001). In contrast, NTD diagnoses increased more in the MOU group than the non-MOU group in the 10-y intervention period (estimate = 0.675; SE = 0.253; P = 0.0076), a trend largely driven by an increase in leprosy diagnoses in the MOU group over time (estimate = 0.864; SE = 0.354; P = 0.015). This may signify true increases in leprosy rates in MOU-signing districts, or may signify increased health seeking behavior for rare and difficult-to-treat diseases, like leprosy, in MOU-signing populations compared to non-MOU populations.

Self-Reported Well-Being and Livelihood Impacts.

Household surveys were conducted by the intervention team in 2007, 2012, and 2017 (see for survey demographic information). Between 2007 and 2017, annual birth rates and infant death rates declined significantly, and although the measurement method used in this survey is not directly comparable to standardized US Agency for International Development (USAID) Demographic and Health Survey (DHS) methods, these declines are consistent with substantial regional declines apparent in DHS data for the same region () (26, 27). As illegal logging declined as a livelihood in the 10-y intervention period, the decline did not correspond with significant changes in unemployment, as employment increased in other sectors (Fig. 3). Monthly household income across all surveyed households in all districts was unchanged from 2012 to 2017 (t test, P = 0.28), but after adjusting income for change in purchasing power parity (PPP) (28), median household PPP-adjusted income decreased by 2.6% (t test, P = 0.001 Fig. 3). However, in rural and low-income settings, national-level PPP adjustment of income may not accurately represent wealth, which might be better estimated by asset-ownership (29). Additionally, household perceptions of neighborhood wealth were not significantly changed over time (Fig. 3).

Dose–Response of the Intervention’s Effect on Deforestation.

We found evidence of a dose–response across 36 villages (dusun administrative units) with an access area >0.30 km2 inside GPNP, whereby forest loss declined with increasing intervention engagement (engagement was defined as the sum of recorded person-contacts across all intervention activities, including clinic patient visits, forest liaisons meetings, conservation education activities, livelihoods training, and a number of other smaller programs; ). Comparing average forest loss rates over time (in three time periods: the preintervention period in 2002 to 2006, the first 5 y of the intervention in 2008 to 2012, and the most recent 5 y of the intervention in 2013 to 2017), forest loss rates near highly engaged villages decreased significantly (−0.15 ± 0.048%; P = 0.007), while forest loss rates near medium engaged villages did not change (0.06 ± 0.042%; P = 0.147), and forest loss rates near the least engaged villages showed an increasing trend (0.16 ± 0.085%; P = 0.067; Fig. 2). There was also a dose–response in the probability that any 30-m2 forested pixel was lost across the entire intervention period: controlling for slope; elevation; distance to nearest river, road, and park edge; and logging pressure (forest loss) outside the park, we found that highly engaged villages’ access areas inside GPNP lost significantly fewer forest pixels compared to that lost in low-engaged villages’ access areas (estimate = −0.85; SE = 0.013; P < 0.0001; Table 2), whereas medium-engaged villages’ access areas lost equivalent forest pixels to low-engaged villages’ access areas (estimate = −0.0087; SE = 0.012; P = 0.46). GPNP forest loss also decreased with average elevation of forest in GPNP access areas (estimate = −1.83; SE = 0.75; P = 0.015; Table 2) and increased with logging pressure outside of the park (estimate = 0.11; SE = 0.0073; P < 0.0001; Table 2). The dose–response of the intervention effect is consistent with a causal association between the intervention—consisting of expanded health care access and use, plus livelihood, education, and conservation programs—and ultimate reduced deforestation outcomes.

Table 2.

Dose–response of forest change to the intervention: Results of a generalized linear mixed-effects regression of forest loss within GPNP over time and the effect of engagement level of each village with the intervention’s programs and activities (see , for details on engagement activities and quantification of engagement levels)

	Log-odds	CI	P
Intercept	−0.12	−8.40–8.16	0.977
Population	−0.47	−1.15–0.20	0.171
Forest lost outside	0.11	0.09–0.12	<0.001
Average elevation	−1.83	−3.31 – −0.36	0.015
Average slope	1.7	−0.68–4.08	0.162
Distance to nearest river	0.47	−0.48–1.42	0.335
Distance to nearest road	−0.12	−1.00–0.76	0.792
Distance to park edge	0.03	−0.60–0.65	0.936
Medium engagement	−0.02	−0.88–0.83	0.955
High engagement	0.80	−0.14–1.74	0.096
Year	0.34	0.32–0.35	<0.001
Interaction terms estimating engagement effect
Medium engagement*year	−0.01	−0.03–0.01	0.456
High engagement*year	−0.85	−0.88 – −0.83	<0.001
Random effects
σ2	3.29
𝜏00: Village	0.83
Marginal R2	0.134
Conditional R2	0.308
No. obs.	108 obs.
	36 villages

Log-odds are presented for centered and scaled predictors. The effect of interest is the interaction of engagement level with year, with log-odds estimates representing the outcome in villages with that engagement level compared to outcomes in low-engaged villages (as a comparison group). Coefficients can be backtransformed to the response scale using a logit link function. Bolded P values represent statistically significant factors.

Discussion

Our results offer objective evidence that increasing access to affordable, high-quality health care as part of a comprehensive conservation intervention—in this case, to rural communities with limited resources and income options living near a densely forested national park in Indonesia—benefits both conservation and human health. In addition, community members self-reported that the intervention was working: By 2012, more than 97% of surveyed households indicated that they believed the intervention was reducing illegal logging. Further insight into mechanisms by which the intervention was reducing illegal logging was gained in 2017 via a household survey question asking, “which programs are most helpful” to stop logging in GPNP. Among the subset of households that interacted with intervention programs, roughly half identified health care discounts alone or in combination with other intervention activities (representing a plurality of responses) as the most important incentive to reduce illegal logging in the park, roughly one-quarter identified livelihood programs alone or in combination with other activities (including health care) as most important, while only a few (6%) indicated that the intervention is not effective at reducing illegal logging. Further investigation is required to establish whether this approach may be effective in other tropical forest parks where high tree cover, high poverty, and lack of access to affordable, high-quality health care fuel illegal logging and forest loss, even within protected areas (30).

Early and continued collaboration with local communities, who identified mechanisms driving linked health–environment problems and potential regional solutions, was essential to the intervention’s multisector success. Globally, about 35% of protected areas are traditionally owned, managed, used, or occupied by indigenous and local communities, yet the perspective and guidance of indigenous peoples and local communities is rarely considered in the design of conservation and climate mitigation programs (31). As a result, many interventions have had negative consequences for local communities that rely on natural resources for subsistence (31). Incentive-based conservation approaches, developed to integrate community development and conservation, have had mixed success, as benefits are not always distributed equitably or do not reflect community needs (32). In contrast, we found that community leadership in the design and implementation of a conservation intervention focusing on pressing health and well-being needs resulted in strong positive benefits to local communities as well as to global conservation goals.

This work demonstrates an actionable framework for aligning cross-sectoral goals. Frameworks such as this are urgently needed to advance effective policy efforts aimed at achieving the United Nations’ Sustainable Development Goals (SDGs) (31, 33, 34). Here, we evaluated outcomes related to conservation (Life on land, SDG 15) and health (Good health and well-being, SDG 3) resulting from an intervention that actually addressed several additional goals, including Climate action (SDG 13), Decent work and economic growth (SDG 8), and Partnerships for the goals (SDG 17). Because the SDGs are deeply interconnected, there is both opportunity and urgency to address multiple targets at once. This intervention offers a case study of how programs can be designed, implemented, and evaluated to address health and conservation goals simultaneously.

The forest carbon results reported here do not include measures of belowground carbon conservation in mineral soils or peatland, the latter of which stores more carbon than aboveground forest biomass in Borneo (35) and is particularly vulnerable to carbon loss and subsidence following deforestation events (36). We also do not include measures of forest regrowth in preserved areas or previously degraded areas being restored through intervention activities (37), which undoubtedly amplified carbon storage and sequestration benefits of the intervention. Furthermore, over the long term, preserving and restoring forest-related ecosystem services might also benefit human health by reducing the risk of waterborne diarrheal disease (38), lowering heat stress (39), or reducing vectors of malaria and arboviruses (40). Measuring these longer-term effects of ecosystem integrity on human health remains an important goal for future linked conservation and public health interventions.

A more nuanced assessment of how health care–conservation exchange programs influence disease occurrence is another important future direction for research. In the context of this study, clinic health records offered a rich dataset on more than 1,250 unique ICD-10 (10th revision of the International Statistical Classification of Diseases and Related Health Problems) codes detected in the patient population (). However, ethical and logistical constraints prevented the establishment of a control group for evaluating health care outcomes: Denying health care access to certain individuals was antithetical to the intervention’s goal to improve health and well-being, and measuring disease occurrence for hundreds of ICD-10 codes in a group of nonpatient individuals was unrealistic. Therefore, we cannot fully account for the contribution of 10 y of regional improvements in government health systems and infrastructure development. However, by comparing clinic usage and diagnosis trends in MOU-signing patients receiving clinic discounts vs. non-MOU signing patients, we established that the intervention incentivized increased health-seeking behavior (with higher clinic usage in MOU-signing districts) and led to potential benefits for common diseases including respiratory infections and dental diseases. In the future, statistical inference on health outcomes of improved clinic access and affordability could be achieved through 1) close collaboration with existing health care facilities that serve distinct patient populations not exposed to the intervention in order to generate comparable datasets, and/or 2) adherence to pragmatic experimental designs, such as stepped wedge cluster randomized trials, that randomly expose groups to an intervention incrementally, thus generating a varying number of control groups at each time point while eventually exposing all groups to the intervention (41).

Deforestation in tropical rainforests has doubled since 2008 (1). Tropical Asia contains some of the most carbon-dense forests in the world (), and Indonesia has ranked consistently among the top countries for forest loss worldwide, not far behind much larger and wealthier countries including Brazil, Russia, and the United States (1). Meanwhile, human health and livelihoods are intimately linked to environmental change. Here, we show that amid this challenging context, local community stewards are both critical actors in, and beneficiaries of, integrated conservation and health solutions. Where health care access is limited and the conservation value of tropical forests is high, reducing rural health care gaps through conservation–health interventions may offer a synergistic means to enhance health and well-being benefits to local communities, while simultaneously conserving critical forest carbon and biodiversity resources.

Materials and Methods

We examined the effect of the intervention on rates of forest loss in GPNP compared to other national parks across Indonesia, using an ex post facto research design (nonrandomized control groups designated after the fact) and a synthetic controls analysis, with nearly two decades of earth observation data quantifying forest change. “Annual forest loss” and “total forest cover” were downloaded from the Hansen Global Forest Change dataset (1) [version 1.6, 2000 to 2018; accessible through the Google Earth Engine data repository (42), 30-m pixel resolution]. Area of forest lost by year was extracted by district-level administrative unit (desa) in Google Earth Engine for all districts whose boundaries intersected GPNP and all other Indonesian national park boundaries.

We used time-varying and time-invariant characteristics to match treated units (districts in GPNP) to untreated units (districts in other parks) to assemble the weighted synthetic control group. Time-varying characteristics included forest lost inside and outside park boundaries, forest fires inside and outside park boundaries, and human population density inside and outside park boundaries. Time-invariant characteristics included area of the district inside and outside park boundaries, area of the focal park, marine area of the focal park (to capture information indicating a coast adjacent or a primarily marine park), year established as a National Park, and average slope within the park (to capture ease of logging access). Total forest cover estimates inside national park boundaries were included to represent total forest available to log, and total forest cover estimates outside national park boundaries were included as a proxy for potential logging pressure outside the park.

To estimate changes in human population density during the evaluation period, we extracted population density estimates from WorldPop (www.worldpop.org/) by district and year from 2000 to 2018, both within and adjacent to park boundaries, using the Google Earth Engine platform (42). The effect of forest fires was controlled for using the MODIS Burned Area Monthly Global data product (500 m), which provides the burn status of each 500-m pixel at a monthly resolution. The park characteristics were acquired from the World Database on Protected Areas (10).

We ran synthetic controls models using three different data subsets defining the donor pool of possible control units. These subsets included 1) all districts in terrestrial (nonmarine) parks established before 2001 (i.e., dropping any parks that are designated as marine only parks and those that were established after the start of the deforestation dataset in 2001); 2) all districts in all terrestrial parks (i.e., dropping only entirely marine parks, but ignoring year of establishment); and 3) all districts in all National Parks in Indonesia (i.e., the most inclusive group of possible control districts in National Parks). Models were run using annual data using 2001 to 2007 as the preintervention period (since the intervention was not expected to lead to immediate changes in deforestation rates in mid-2007, when the intervention started), and 2008 to 2018 as the postintervention period.

In each model, P values and 90% confidence intervals (for a one-tailed lower test) were calculated using a standard normal sampling distribution and Taylor series linearization to estimate the variance and produce CIs (21). In addition, P values and CIs were also calculated using 500 permuted placebo treatment groups for comparison with the estimated effect for the actual treatment group to satisfy a more robust set of assumptions and to generate a more robust and conservative estimate of the sampling distribution (21). These “permutations” are placebo tests in that they randomly assign districts in the “control” group to the placebo treatment group, the synthetic controls model is rerun, and the magnitude of the placebo treatment result is compared to the actual treatment group result (21). All models were run using the “microsynth” package in R (43, 44) following established methods outlined in Robbins et al. (21).

We evaluated the intervention impact in terms of forest loss rates, and then estimated the quantity of aboveground forest carbon conserved from the quantity of forest loss averted estimated by the synthetic controls analysis. Following Jucker et al. (22), we estimated forest carbon stocks in GPNP using canopy heights derived by LiDAR (light detection and ranging), accessed from National Aeronautics and Space Administration and Oak Ridge National Laboratory Distributed Active Archive Center ( and see , for conversion equations and details) (45). Then we used the effect size of 69.8% forest-loss-rate reduction to estimate the total aboveground carbon stock (in teragrams carbon per hectare) conserved in the period from 2008 to 2018. This method is conservative, as using the mean pixel height of the LiDAR flight is an underestimation of the average top-of-canopy height of trees targeted by illegal loggers, who target the largest and most valuable trees. Even so, carbon densities were high, in part, because this region of Southeast Asia contains forests with some of the largest and most carbon-dense trees in the world ().

The Alam Sehat Lestari (ASRI) medical center opened in July 2007 and remains open. Our analyses consider the period from 2008 to 2018, beginning with the first full year of data and ending with the last full year of data before comprehensive evaluation began. For all patients who visited the ASRI medical center, patient records included only a unique ID to maintain patient privacy, the date of the visit, the patient’s home village (dusun administrative unit) and district (desa administrative unit), the patient’s age, and the diagnosis (ICD-10 code) given during the visit and/or the reason for visiting (e.g., medical checkup, fill prescription). Acquisition and analytical use of fully de-identified health clinic records was submitted for review to the Institutional Review Board of Stanford University and was determined to contain no identifiable data, requiring no further review as human subjects research. The clinic data were gathered within routine operations of the ASRI medical clinic (Indonesian nonprofit Alam Sehat Lestari is registered by the Indonesian Ministry of Cultural Affairs, #AHU-08962.50.10.2014; clinic operations are permitted by the Kayong Utara Department of Health in Indonesia).

We were interested in understanding how increases in clinic affordability (discounts provided to patients from MOU-signing villages/districts) impacted clinic usage, while controlling for clinic accessibility (estimated mean travel time from the patients’ district to the clinic). We quantified clinic usage at the district-level in two ways: 1) the frequency of patient visits, defined as the number of visits per patient (i.e., the number of unique months that a patient occurred in the database; see for details); and 2) the proportion of the population in each district that used the clinic at least once during the evaluation period. Clinic access was estimated for each district as the mean travel time (in minutes) from 10 randomly distributed points in each district, weighted by the point’s population density (see , for details).

To understand how patient-level visit frequency was affected by affordability (MOU status) and access (travel time), we ran a negative binomial GLMM with a random effect for district that quantified how the number of visits per patient varied with MOU status and the estimated travel time (in minutes) from the patient’s district to the clinic. Next, we used a Poisson generalized linear model to quantify how the number of unique patients per district varied with MOU status and estimated travel time to the clinic, where the 2018 population size in each district was included as an offset. For this analysis, we excluded patients from unknown districts (recorded only as “far” in the patient records) because district-level population size for this group of patients could not be determined. For both clinic usage analyses, we were unable to estimate travel time for one island district (Pelapis), and therefore excluded it from analysis. Details on how travel time was estimated for the far () districts are available in .

Overall, 1,255 unique ICD-10 codes (46) were applied to at least one patient between 2008 and 2018. Before analyzing how the proportion of unique patients that received a diagnosis (ICD-10 code) of a particular disease changed over the intervention period, we classified 824 ICD-10 codes that contributed to the most common disease categories (see , for details) into the following groups to be tracked: childhood-cluster diseases, COPD, dental disease, diabetes, diarrheal diseases, heart disease, liver disease, lower respiratory infections, upper respiratory infections, malaria, malnutrition, NTDs, trauma, and tuberculosis. Other unspecified diseases were grouped and appear in the text and figures as “untracked” (see , for a full list of the 824 ICD-10 codes tracked), and the unique patients to which the untracked ICD-10 codes were assigned are included in the total patient population.

To track changes in disease occurrence in the patient population, we estimated the proportion of unique patients that received a diagnosis for each disease in each district, annually. Each patient only counts toward one instance of any particular disease per year, and the denominator is the total unique patients that received any diagnosis in a year. For each disease, we used a binomial GLMM with a logit link (which yields coefficients that can be exponentiated to derive odds ratios) and a random effect for district to quantify differences in the proportion of disease diagnoses over time (early, 2008 to 2009, vs. late, 2017 to 2018) for two populations: patients from districts with and without MOUs. The model was run with and without controlling for district distance to the clinic, with nearly identical outcomes. The full time series showing changes in the period prevalence of each disease in the patient population are shown in .

We also tested whether MOU status impacted disease outcomes among patients over time, while controlling for a patient’s average distance (in minutes) from district to clinic (see for details on calculating travel time). To do so, we used binomial GLMMs with a logit link, using probability of a diagnosis of each disease as the outcome, and including scaled travel time and an interaction term for time by MOU status as predictors, and a random effect for district to control for repeated-measures and unmeasured district-level effects.

At baseline in 2007 and at follow-ups in 2012 and 2017, detailed household surveys were conducted in districts surrounding the national park (; and see , for details on surveyor selection and training). Within each district, a list of all households in the villages was provided by village heads and from that list, ∼10% of households (and, correspondingly, ∼10% of the total population of ∼60,000 people) were randomly selected for participation. In total, 1,348, 1,498, and 1,379 households were surveyed in 2007, 2012, and 2017, respectively (see , for details on the surveyed household demographics). At each time point, the survey instrument contained modules for the following: demography (age, sex, births, etc.), health, wealth (income), perceptions of wealth (designation of the household or neighborhood as “wealthy,” “average,” or “poor”), livelihoods (including logging activity and other occupations), and perceptions of nature, natural resources, conservation, and the intervention. Acquisition and analytical use of fully de-identified household survey data was submitted for review to the Institutional Review Board of Stanford University and was determined to contain no identifiable data, requiring no further review as human subjects research. Surveys were administered as part of routine monitoring and evaluation of ASRI program activities in 2007, 2012, and 2017 (Indonesian nonprofit Alam Sehat Lestari is registered by the Indonesian Ministry of Cultural Affairs, #AHU-08962.50.10.2014), following approval by local heads of participating desa and “kecamatan” administrative units in Indonesia.

From demographic characteristics, including ages, gender, births, and deaths in surveyed households, we calculated average annual infant mortality rates, defined as average annual infant deaths per 1,000 live births among household women in the 3-y period preceding the survey, and average annual births, defined as average annual births per women ages 19 to 59 (). For reference, we also extracted infant mortality rates (IMR) and general fertility rates (GFR) in a 5-y period from USAID Demographic and Health Survey data for Indonesia in 2007 and 2017 (26, 27). We report IMR and GFR for the province of West Kalimantan and for “rural” West Kalimantan, which we expect to be a more accurate representation of the rural communities in and around GPNP (). From reported incomes and household perceptions of neighborhood wealth, we calculated the average monthly income at the household level and the proportion of households that felt neighborhood wealth had increased, decreased, or remained the same at each time point. We also calculated the proportions of households that reported members engaged in various livelihoods (logger, fisher, farmer, civil servant, company employee, unemployed, and “other”). Last, to infer mechanisms by which the intervention may have reduced illegal logging rates, we calculated household responses to a 2017 survey question that specifically asked what intervention programs provide the strongest incentive to help stop logging inside GPNP.

For reported livelihoods and perceptions of neighborhood wealth, we quantified change over time using generalized linear models with binomial error distributions and logit links. For livelihoods, we computed the percent change in the proportion of households reporting each livelihood as a primary income source derived from the log-odds that livelihood changed over time. We used t tests to test for a change in average monthly income and PPP-adjusted monthly income over time. Due to slight variations in survey question wording, only 5-y comparison data (2012 vs. 2017) were available for monthly incomes and perceptions of neighborhood wealth.

Dose–Response of the Intervention’s Effect on Deforestation within GPNP.

After demonstrating a significant correlation between the intervention and pre-intervention to postintervention forest loss trajectories in GPNP compared to a synthetic control, we tested whether there was any evidence of a dose–response relationship within GPNP, among villages (dusun) with varying levels of engagement with the intervention programs (including use of the health clinic, and other periodic programs; ) and forest loss rates. To answer this question, we quantified engagement effort as cumulative person-contacts (i.e., number of contacts with persons reached by all program activities associated with the intervention from 2007 to 2017, allowing for repeated contacts with the same individuals over time) achieved through the following: the health care intervention (ASRI clinic visits, mosquito net distribution), conservation programs (Community Conservation Liaisons or “Forest Guardian” Program, Chainsaw Buyback Program, and Reforestation Program), alternative livelihood trainings (Organic Agriculture Program, Goats for Widows Program, Green Kitchen Program), and education activities (ASRI Kids Program, Community Education Program) (). Engagement effort was not distributed evenly across all villages and was predominated by frequent engagement with community liaisons for the intervention as well as doctor–patient contacts at the clinic (). Variation in engagement across the participating villages intersecting GPNP allowed us to test for evidence of a dose–response of intervention effort on deforestation within different access areas nearest each village around the park.

We used a k-means clustering algorithm to bin engagement [cumulative person-contacts in each village across all of the intervention programs from 2007 to 2018 ()] into low, medium, and high categories [R package “classInt” (47)]. We examined the effect of cumulative engagement effort on the proportion of forest lost in each village’s access area in the national park (number of remotely sensed 30-m2 pixels lost out of total forested pixels remaining). Village-level access areas inside GPNP were determined by a local team that mapped the parts of each village bordering GPNP that extended into GPNP and represented that village’s typical access area for illegal logging. Ultimately, 36 villages bordering GPNP with logging access areas >0.30 km2 were included in the analysis. As in the synthetic controls analysis, forest change data were obtained from the Hansen Global Forest Change dataset [version 1.6, 2000 to 2018, 30-m pixel resolution (1)]. Total forest cover by village was estimated annually by subtracting forest loss in each subsequent year from total remaining forest cover in the previous year (e.g., forest cover in 2002 = forest cover in 2000 – forest loss in 2001 – forest loss in 2002).

Changes in forest loss rates associated with engagement were estimated in two ways. First, for each engagement category (high, medium, and low), we estimated change in average annual forest loss rates over time, from before the intervention (2002 to 2006) to forest loss rates during the first 5 y (2007 to 2012) to the most recent 5 y (2013 to 2017) using a mixed-effects linear model with time as a predictor and nested random effects of village within district repeated over time. Next, we fit a binomial GLMM to estimate the probability that, over time, any 30-m2 forested pixel in a medium- or highly engaged village’s access area within GPNP was lost, compared to the probability that any 30-m2 forested pixel was lost in a low-engaged village’s access area. We controlled for village population size (supplied by village leaders in 2017), proportion of forest lost within village boundaries outside of the park (as a proxy for outside logging pressure), average slope and elevation of pixels inside the park, average distance of pixels to the nearest river, road, and park boundary, and nested random effects of village within district. We did not include data on fires because fire activity was relatively low in GPNP during the time period under consideration.