Replanting unproductive palm oil with smallholder plantations can help achieve Sustainable Development Goals in Sumatra, Indonesia

Oil palm ( Elaeis guinensis ) is a controversial crop. To assess its sustainability, we analysed the contribution of different types of plantations (smallholder, industrial and unproductive) towards meeting six Sustainable Development Goals. Using spatial econometric methods and data from 25,067 villages in Sumatra, Indonesia, we revealed that unproductive plantations are associated with more cases of malnutrition, worsened school access, more air pollution and increased criminality. We also proposed a strategy for sustainable palm oil expansion based on replanting unproductive plantations with either industrial or smallholder palm oil. Smallholder replanting was bene ﬁ cial for ﬁ ve Goals (Zero poverty, Good health, Quality Education, Environmental preservation and Crime reduction), while the same intervention only improved two Goals in the industrial case (Zero poverty and Quality Education). Our appraisal is relevant to policymakers aiming towards the 2030 Agenda, organisations planning oil palm expansion, and retailers or consumers concerned about the sustainability of oil consumption.


I
n 2015, all United Nations Member States adopted the 2030 Agenda for Sustainable Development, an urgent call for action to all countries to achieve sustainable development by 2030 1 .The Sustainable Development Goals (SDGs) are at the core of the initiative -17 guiding principles encompassing economic growth, environmental protection, global peace, access to healthcare, basic education and others that should be followed by all countries to aim towards a more sustainable future.Since its adoption, many public and private organisations have also relied on the SDGs to evaluate the sustainability of their activities.Nonetheless, this analysis can be challenging, especially when the activities evaluated can have complex outcomes for many different development dimensions.Palm oil production in tropical countries is a typical example of such complex interactions between different SDGs 2 .
The palm oil industry is frequently criticised due to environmental and social impacts [3][4][5] .Specifically, for contributing to air and water pollution 6 , biodiversity loss 7 and a decline in social well-being in neighbouring communities due to land use conflict and labour exploitation 8 .However, palm oil production is also an important vehicle for economic growth and poverty alleviation with recent large-scale impact evaluations emphasising the tradeoffs associated with the industry's influence on the SDGs 7,9 .As global demand for vegetable oils increases from 201 to 307 million tones (Mt) by 2050 4,10 , identifying the role of palm oil production in achieving the multifaceted elements of the SDGs will become crucial to guide policymakers towards more sustainable production systems, given its relatively low land requirements compared to other oil crops 4 .
As the largest worldwide palm oil producer, Indonesia provides the ideal case study to evaluate the contribution of palm oil to the SDGs 11,12 .The Indonesian economy has thrived through palm oil production, with 17% of Indonesia's agricultural gross domestic product in 2014 depending on this crop 13,14 .However, its expansion has been associated for many years with deforestation, ecological damage and a mix of positive and negative social outcomes [3][4][5]9,11,15 . The Indnesian scenario is even more complex when considering the heterogeneity between the industrial and smallholder palm oil production profiles.Industrial plantations tend to be large, with thousands of hectares of land dedicated to the monoculture of palm oil.Smallholder plantations, on the other hand, are typically smaller -<25 hectares according to government definitionalthough wealthy individuals can sometimes own several hundred hectares.They also tend to exhibit a more heterogeneous spatial pattern, where palm oil is often mixed with other crops or land dedicated to other uses 16,17 .The differences between both production profiles could also entail discrepancies in their ecological footprint 18 and the socioeconomic characteristics of the villages where they are located.
In this study, we evaluate the differences between smallholder and large-holder industrial plantations through the perspective of the SDGs.The analysis relies on spatial econometric models to reveal correlations between the socio-economic characteristics of the villages in the region of Sumatra (Indonesia) and their land dedicated to different palm oil production profiles (active industrial plantations, active smallholder plantations and two types of unproductive plantations).To characterise the contribution of each production system to the SDGs, we explore the differences in seven SDG indicators if the unproductive (damaged) plantations in Sumatra were replanted through a smallholder or an industrial profile.Our analysis reveals how successfully replanting damaged plantations with smallholder palm oil exhibits more beneficial outcomes toward most SDGs than its industrial equivalent.However, smallholder palm oil production is also responsible for most unproductive (damaged) plantations (72% of the total damaged land has smallholder tenure), which indicates the need for careful planning and additional policies to ensure the viability of the replanting strategy.Our appraisal could be relevant for policymakers and stakeholders aiming to complete the 2030 Agenda in the context of palm oil production in Indonesia.

Results
Socio-economic differences between palm oil production systems.Palm oil production can play a major role in the socioeconomic development of the villages in which is produced.However, it is still unclear how these impacts could vary depending on the production profile used for its production.To this end, we combined two datasets describing, respectively, the socio-economic characteristics of the villages of Sumatra (Indonesia) and their total area dedicated to palm oil plantations, classified by their type of ownership (industrial/smallholder) and production status (young, mature and unproductive damaged or unproductive but possibly replanted) 16,17 .See Methods section for a more detailed description.
We aggregated the palm oil data into four production profiles: active industrial palm oil, active smallholder palm oil, unproductive damaged plantations and unproductive but possibly replanted plantations (see Methods for more details).Then, we estimated the socio-economic characteristics significantly associated with each profile through two types of spatial regression analysis.The Spatial Durbin Model (SDM) 19 , a global spatial dependency model and the Spatial Durbin Error Model (SDEM) 19 , a local-dependency variant.Both models produced broadly consistent results (see Supplementary Tables 5-8).
Unproductive (damaged) plantations are significantly associated with more cases of malnutrition, crime, air pollution and a higher distance to primary education centres (see Fig. 1).These socio-economic characteristics are detrimental to SDG 2 (Zero hunger), SDG 4 (Quality education), SDG 15 (Life on land) and SDG 16 (World peace).Damaged plantations are also positively correlated with more healthcare facilities and more types of credit facilities, two features that are beneficial respectively for SDG 3 (Good Health and well-being) and SDG 8 (Economic growth).Nevertheless, these positive associations are not unique to unproductive (damaged) land.They were also observed in the case of active smallholder plantations, which could indicate that such infrastructures are built during the early stages of plantation development and remain after their failure.These characteristics differ from the ones observed for unproductive (but possibly replanted) areas.In that case, we observe fewer negative effects overall and slight beneficial outcomes for poverty reduction.
Villages with industrial or smallholder plantations share some common characteristics compared to non-oil palm producing villages: they are less associated with natural disasters and tend to have fewer households receiving Surat Keterangan Tidak Mampu (SKTM) poverty letters.These letters are issued to families below the national poverty level to facilitate access to essential resources such as fuel, healthcare or education 20 .Prior studies have already reported a positive influence of palm oil production on poverty alleviation 9,11 .In our case, these beneficial economic outcomes are only observed for villages with active plantations from both profiles and for possibly replanted areas, which highlights the negative role of damaged plantations in sustainable development.
Additionally, we found that several characteristics distinguish villages with smallholder and industrial plantations.Those with a smallholder profile are associated with more credit facilities and cases of malnutrition, while in the case of industrial land, the reduced distance to primary schools is one of the main differential traits.The diversity observed in both production profiles evidences the need to evaluate their sustainability separately.This is especially relevant for planning interventions over land use, as the same intervention could entail different sustainability outcomes depending on the village's predominant production system.

Counterfactual analysis.
With the global increase in vegetable oil demand and related environmental impacts of agricultural expansion, it is crucial to develop strategies to increase palm oil production while maintaining Indonesia's efforts toward net zero deforestation 21 .With damaged palm oil plantations having multiple detrimental outcomes on sustainable development, replanting them could be a key strategy for sustainable palm oil expansion.Nevertheless, it is still unclear whether smallholder or industrial palm oil is better at maximising the beneficial outcomes for the SDGs.We, therefore, implemented a counterfactual analysis to estimate the change in seven SDG indicators if all unproductive (damaged) plantations were replanted through an industrial or smallholder profile.The seven indicators used for the counterfactual analysis span six different development goals (see Methods for the selection criteria).These are SDG 1 (Zero poverty): log(num.SKTM poverty letters), SDG 3 (Healthcare and Well-being): Num.healthcare facilities, SDG 4 (Quality Education): log(dist.primary school) and log(dist.junior-high school), SDG 9 (Industry, innovation and infrastructure): log(dist.nearest market), SDG 15 (Life on land): percentage of village area covered by forest, SDG 16 (Peace and Justice): number of different types of crime.
The counterfactual analysis followed the pipeline described in Fig. 2. The SDM model was trained seven times, using each SDG indicator as the dependent variable y and all the village characteristics (socio-economic and palm oil production features) as independent features X.The trained models were used to predict the expected change in the SDG indicators if all unproductive (damaged) land in Sumatra was replanted with either industrial or smallholder plantations.
The predictions obtained for each SDG indicator can be visualised as the counterfactual maps shown in Supplementary Fig. 2.These maps provide a detailed description of the location and magnitude of the changes produced in the SDGs, but they do not indicate how the intervention affects the well-being of local communities.For example, removing one hospital from a village without other healthcare facilities might have a more drastic impact on the villager's well-being than the same change in a village with other healthcare centres.
To contextualise the results in terms of the village's well-being, we defined a deprivation threshold for each SDG indicator using the reference values proposed in the Indonesian Village Development Goals (Indeks Pembangunan Desa) 22 (see Methods section for more details).Interventions that induce a change from a deprived state to a non-deprived state or vice-versa will be considered highly positive or highly negative contributors to the SDGs, as they drastically affect the well-being of the villages.If the intervention does not induce a change in the deprivation state, it will only be considered a moderate influence, regardless of the magnitude of the change observed.Finally, interventions where no changes in the SDG indicators are induced, are considered neutral interventions.Using these criteria, the counterfactual maps were translated into Fig. 3, which summarises the effects of both interventions across all villages.
Replanting unproductive (damaged) plantations is associated with a positive contribution towards SDG 1 for both production profiles.The smallholder strategy resulted in a slightly more beneficial outcome, as it is associated with highly positive effects for more than 27% of the villages where the counterfactual was applied.In the case of SDG 3, replanting unproductive land with industrial plantations was associated with a reduction in the number of hospitals, a negative influence on the SDG, while smallholder expansion is positively impacting the objective.This outcome must be further validated as the apparent lack of medical coverage in industrial plantations could be because some company-owned health facilities may not be included in the socio-economic data.
For SDG 4 both production profiles have mainly negative outcomes towards reducing the distance to junior-high schools.However, the industrial strategy seems to result in more beneficial effects towards improving the distance to primary school centres.Thus, in this case, promoting industrial plantations could be more beneficial overall for SDG 4, also because similar to health care, many companies provide schooling for children of their workforce 23 .
Increasing industrial plantations is also associated with a small reduction in forest coverage, leading to a moderately negative  Overall, replanting damaged plantations was more beneficial when implemented through a smallholder profile than its industrial equivalent.Smallholder replanting was beneficial for SDGs 1, 3 and 15 and it exhibited mixed effects for SDGs 16 and 4 (positive effect for primary schools and negative for junior-high schools).Meanwhile, industrial replanting only presented synergistic associations with SDG 1 and one facet of SDG 4 (primary school).This last association was the only one more beneficial in the industrial replanting scenario than in the smallholder case.Interestingly, both counterfactual scenarios were not compatible with the progress of SDG 9 and SDG 4, in the case of the juniorhigh school indicator.This revealed the presence of a trade-off effect between the intervention proposed and the achievement of some development goals.The findings demonstrate the intricate nature of the connections between an intervention and the SDGs.An identical intervention can yield positive outcomes for certain indicators while having negative effects on others, leading to various synergistic or trade-off relationships with the development goals.

Industrial expansion in Sumatra
Province analysis.To assess the scale dependency of the findings, we repeated the counterfactual analysis for each of the ten provinces in Sumatra individually (see Supplementary Figs.3-12).We then estimated the production profile that yielded the most favourable and detrimental outcomes for each province and compared it to the results observed for the entire island.
The results presented in Table 1 confirmed the trends observed in the study for all of Sumatra.Replanting by smallholders was associated with the most beneficial outcomes for the SDGs across all provinces.The only exception was the case of SDG 4 (for the primary school indicator), where industrial replanting consistently outperformed smallholder replanting across the regions evaluated.This analysis also revealed consistent results in their trade-off effects, as replanting damaged plantations also entailed detrimental outcomes for SDG 9 across all provinces.
In general, the proposed counterfactual approach enabled us to differentiate between two seemingly similar interventions concerning palm oil usage.This tool could allow policy-makers to make more informed decisions to achieve sustainable palm oil development.However, it is crucial to recognise that the interplays described by the counterfactual evaluation are conditional on successfully transforming unproductive terrains into active plantations.Thus, guaranteeing the viability of the replanted plantations is of utmost importance for the sustainability of the strategy.
Sumatran villages tend to have more unproductive (damaged) plantations associated with smallholder tenure than with industrial palm oil (see Fig. 4).In terms of total area, 72.5% of the total unproductive damaged land is associated with the smallholder production profile (calculation in the Supplementary Methods).It is crucial to identify the mechanism behind the increased failure rate of smallholder plantations.Whether it is caused by sub-optimal soil and planting materials 24 ; or if indeed smallholder producers sometimes lack the capital or knowledge to make plantations thrive 24 .This could reveal if replanting damaged plantations with a smallholder profile is indeed sustainable by itself, or if other supplementary programs (e.g., educational, economic, etc.) should also be introduced to increase its effectiveness.

Discussion
The sustainability of palm oil production is frequently called into question, especially by commentators in high-income countries.This criticism focuses strongly on the environmental impacts of oil palm but it is less clear the extent to which negative impacts on the environment are balanced by positive impacts elsewhere.This is at the core of sustainable development considerations.The SDGs are intended to be targeted simultaneously, resulting in positive outcomes for all these objectives.In reality, trade-offs are, however, likely to occur.Here we have analysed the extent to which oil palm development in Indonesia, the world's largest palm oil producer, results in synergies and trade-offs between different SDGs, and how this varies in different production profiles.This analysis provides a basis for a better understanding of the overall impact of oil palm on sustainable development and also a way to steer future developments that maximise positive outcomes and minimise negative ones.
Oil palm expansion, if indeed sought by the Indonesian government, should preferably occur on existing agricultural land, or in areas where previous attempts to cultivate oil palm have been ineffective 25 .Replanting failed or poorly performing plantations with smallholder palm oil was associated with beneficial effects for poverty, healthcare provision, terrestrial environments, primary education access and in some villages crime variability reduction.This indicates a synergistic association between the proposed intervention and SDGs 1, 3, 4, 15 and 16.Meanwhile, industrial replanting only entailed beneficial outcomes for education and poverty (SDGs 4 and 1).These outcomes are quite consistent across all provinces and would indicate that smallholder production may be more beneficial for the SDGs than industrial plantations.
However, we also observed that 72,5% of unproductive (damaged) land in Sumatra was associated with smallholder tenure.The beneficial outcomes of the intervention proposed are conditional on transforming unproductive land into viable plantations.Thus, it is crucial to ensure the replanted areas do not fail again after the intervention.To this end, policymakers should invest in comprehending the mechanism behind the increased Table 1 Summary of the palm oil production system with the most beneficial and detrimental effects in all provinces of Sumatra.
"S" indicates that smallholder production has the largest (positive or negative) influence on the SDG indicators."I" indicates that industrial production is the one with the greatest (positive or negative) influence on the SDG indicators."I/S" indicates that both production systems exhibit the same positive or adverse effects.Missing values ("−") indicate cases where no positive or negative associations are found.The geographical location of each province can be found in Supplementary Fig. 1.The SDG icons are property of the United Nations, reprinted for informational purposes in accordance with ref. 67 (all rights reserved).
failure rate of smallholder production to determine if it is linked to land quality or management issues.
Both mechanisms are often difficult to differentiate.Peatlands and hilly areas are sub-optimal for oil palm, the former because of low pH and low productivity (and high fertiliser needs), the latter because of potential lack of water and soil erosion 26,27 .These are often marginal lands where aspiring farmers might try to start a plantation, but if the costs for generating good yields are too high they abandon it 5 .Complicating this is the fact that many smallholders opt for low-input low-output systems 28 that can generate similar profits to more intensively managed plantations but may have higher failure rates (i.e., the cost of walking away from a plantation is lower when investments have been low).Palm oil can also fail due to natural factors, like the presence of diseases, fire 29 or floods 30 .In that case, it is necessary to comprehend why some landowners may not be rapidly re-purposing plantations after their failure.This could be linked to land speculation practices 31,32 , contested land ownership, insufficient funding, or the lack of technical skills to re-purpose the plantation.Overall, comprehending the causes behind plantation failure could reveal the locations most suitable for the replanting strategy proposed and whether supplementary policies (educational programs or financial/technical support) are required to maximise its success chances.
Even if replanting damaged plantations with smallholder palm oil could aid sustainable development in Sumatra, these conclusions should not be directly extrapolated to other areas.The region of Sumatra is in a more advanced stage of oil palm development compared to other parts of Indonesia, such as Kalimantan and Papua 16 .Additionally, it is also the region with the largest amount of unproductive plantations, making it an ideal scenario to test the replanting strategy.If there is interest in assessing the sustainability of similar interventions in other Indonesian regions, or other countries, the whole counterfactual analysis should be repeated using adequate data.
Despite its caveats, our study provides valuable insights for government authorities and policy-makers involved in planning oil palm and other monocultural plantation crops.We offer the possibility of evaluating simultaneously the impact of an intervention over seven different SDG indicators, revealing some of the existing synergies and trade-offs among them.With adequate data, this approach could be expanded to evaluate all 17 SDGs, providing a handy tool to evaluate the sustainability of palm oil production as a whole.Moreover, bringing together land cover, production scale and socio-economic indicators is an important step towards improving understanding of how agricultural development can have the greatest positive SDG impacts and avoid negatively affecting people, especially those in forest environments.
Our study also opens new lines of research in different fields.For example, it could be used as a hypothesis-generating tool to design on-site sociological studies.A similar approach combining data-driven models and experimental studies is starting to be used in SDG assessment 33 .However, to the best of our knowledge, it has not yet been applied in the context of palm oil research.Complementing our study with household surveys could reveal the causal linkages behind the correlations identified in this study.It could also shed some light on the impacts of palm oil expansion on the villages' minority groups.
The model proposed could also be used to test other types of interventions in palm oil use.For example, there would be value in studying the social and environmental benefits of reforesting damaged plantations.This could help clarify if less intensive land uses (agroforestry cultivation) would be preferred in unproductive plantations with poor-quality soils.Another interesting counterfactual scenario could be assessing the impact of climate change on sustainable palm oil production.Tropical production areas are predicted to become drier and hotter 34 , causing alterations in rainfall patterns and water stress.This makes oil palm more prone to pest and disease outbreaks, potentially reducing its yields.Other risks include sea-level rise in coastal production areas, especially on peat soils 35 , and increased wildfires 36 .These changes will likely influence interplays between palm oil production and the SDGs and should be modelled specifically to guide sustainable policies.
Other lines of research could be directed towards improving the quality of the data used in this analysis.The socio-economic data for 2018 was obtained from the Village Potential Survey (PODES) survey, an Indonesian survey of village well-being with a much higher spatial resolution than the one usually reported in the SDG Global Database (the United Nation's official countrywide SDG reports) 37 .However, it only contains a small subset of variables that can be used as SDG indicators, only providing partial insight into each goal.For example, the use of the distance to education centres as a proxy indicator for SDG 4, does not capture other facets of this goal like the quality of the education provided.Improving the resolution of the socio-economic data could allow us to better characterise the SDGs from a multidimensional perspective.
Regarding the palm oil dataset, our study would significantly benefit from an improved ownership categorisation.The binary differentiation between smallholder and industrial plantations could be masking different sustainability outcomes between smallholders with different characteristics.This categorisation also disregards whether the owner resides in the same village where the plantation is located.However, we worked under the assumption that many of the costs and benefits of palm oil accrue irrespective of the owner's residence (e.g.ecological impacts).We also chose not to incorporate palm oil trade as a predictor in our analysis, despite its role in driving palm oil production through international and national demand.Estimating the influence of demand on individual villages would require extensive surveys to gather specific data, and we believe that palm oil demand does not have a distinct impact on the SDG indicators at the village level beyond its influence on palm oil production itself.
Our study also provides some methodological improvements to ease the introduction to the field of spatial econometrics for researchers in the palm oil community.The use of spatial data in palm oil research has thrived in recent years 3,16,17,26,38 .However, few studies have decided to follow spatial econometric approaches to evaluate the sustainability of palm oil plantations 26 .Instead, the main focus has been on assessing the impacts of interventions retrospectively.This approach is very useful for describing the causal effects of an intervention but it does not take advantage of the spatial nature of the data 3,7,39 .The use of spatial regression models offers the possibility for predictive evaluation, which could be a very valuable tool for policy-makers planning new interventions.To facilitate the use of spatial regression models in sustainability research we have improved the model selection pipeline described in ref. 40 , which was one of the most complete approaches towards model selection in spatial econometrics.
Overall, we hope this study can open the door to more research exploring the environmental and social interests of oil production.The current analysis and methodological approach also allow for future comparison between the SDG impacts of different vegetable oil crops, such as soybean, rapeseed/canola, coconut, peanut or sunflower.All these crops have different oil yields and thus land requirements and come with different sustainability profiles because of their different geographical location.The trade-offs and synergies between different crops and production profiles have never been analysed, making it impossible to compare their sustainability meaningfully.The current study thus provides a first methodological step towards such a global assessment of the sustainability of vegetable oil production.

Methods
Data.We relied on high-spatial-resolution data from 26,505 villages in Sumatra describing four different palm oil production situations, the area of forest and the socio-economic characteristics of each village.The land use data was generated by combining two satellite-derived datasets from refs. 17,41and refs. 16,42, which also specifies their type of ownership (smallholder or industrial).The approach used to classify the plantation's ownership relies on their differences in size, land homogeneity and infrastructure development and it is further described in the Supplementary Methods.
References 17,41 used radar and optical images (Sentinel-1 and Sentinel-2) to detect the area covered by closed-canopy oil palm plantations in 2019.These plantations are already productive, generally when oil palm is older than 3 years.On the other hand, refs. 16,42used historical optical images (Landsat) to identify the year when the land was cleared to plant oil palm (between 2001 and 2019).This analysis considers all land cleared for palm oil production, including immature, damaged and failed plantations; and thus surpasses the estimated extent of the closed-canopy plantations.By combining the two datasets above with an additional analysis of recent replanting (from Sentinel-2 images) we generated the following four plantation classes: 1) mature oil palm corresponds to land cleared for oil palm before 2017 and detected as closed-canopy oil palm in 2019; 2) young oil palm is the land cleared for oil palm after 2017 (inclusive) and not detected as closed-canopy oil palm in 2019; unproductive oil palm corresponds to land cleared for oil palm before 2017 and is not closed canopy in 2019; 4) possibly replanted plantations is the land cleared for oil palm before 2017, is not closed-canopy in 2019, but is detected cleared between 2017 and 2019 if the normalised difference vegetation index from Sentinel-2 images presented a sharp decrease during 2017-2019."Possibly replanted plantations" can be interpreted as failed or damaged plantations established before 2017, cleared and replanted during 2017-2019.We specified separately the four production classes for each production profile (industrial and smallholder), resulting in the eight different types of palm oil plantations detailed in Supplementary Table 2.
To measure the forest area in 2019 at the village level we combined a map of the "primary" forest extent in 2000 43 and an annual forest loss map ref. 44 .The forest cover for the year 2000 was produced by including intact and degraded primary forest types described in ref. 43 .This map was updated to 2019 using the forest change maps in ref. 44 .Areas of forest in 2000 cleared between 2001 and 2018 were marked "cleared previously" and removed in the updated 2019 forest map.
For our regression analysis, the village area occupied by the eight plantation classes and the forest coverage was reported relative to the total village area.Absolute land units were only used for the estimation of the total percentage of unproductive land.
The socio-economic characteristics of the Indonesian villages were obtained from the PODES for 2018 11 .PODES is a national survey carried out by the Indonesian Central Bureau of Statistics every 2-4 years and includes more than 100 questions about the economic, social and environmental characteristics of Indonesian villages.We selected questions of most relevance to the SDGs and synthesised response variables into 27 features, 6 describing general village characteristics and 21 reflecting different indicators directly related to the SDGs.While the PODES survey was conducted in 2018, the oil palm datasets were created for 2019.The datasets have a one-year mismatch, but this mismatch is negligible because the oil palm plantation area does not change substantially from year to year.SDG indicators.SDG indicators are quantifiable characteristics that reflect the progress towards achieving each development goal.Using the "Global indicator framework for the SDGs" 45 , we identified seven features in our dataset associated with some objectives of the SDGs: • Zero poverty: SDG 1, is dedicated to poverty eradication "End poverty in all its forms everywhere" 45 .More precisely, objective 1.2 in ref. 45 states, "By 2030, reduce at least by half the proportion of men, women and children of all ages living in poverty in all its dimensions according to national definitions".The closest indicator in the PODES dataset is the number of SKTM poverty letters.SKTM poverty letters have already been used as a proxy for poverty in other studies 20,46,47 , and they reflect a poverty alleviation system created by the Indonesian government to fight against national poverty.SKTM letters are issued to families below the national poverty level to facilitate access to essential resources such as fuel, healthcare, or education 20 .As such, they should provide a good proxy to determine the amount of poverty present in each village (see Supplementary Methods for a correlation analysis between poverty letters and income).The number of poverty letters is negatively associated with SDG 1 as the objective of this SDG is to reduce the number of poor families; • Good health and well-being: objective 3.8 states the need to "Achieve universal health coverage" 45 .Therefore, we have relied upon the number of healthcare facilities as an indicator for SDG 3.This indicator is positively associated with SDG 3 since increasing the number of hospitals is beneficial to the development goal; • Quality education: SDG 4 considers access to primary and secondary education 45 .In this case, two variables were considered as SDG indicators, the distance to primary schools (km) and the distance to junior-high schools.Improving access to education facilities is crucial to ensure that all youth achieve literacy and numeracy, as proposed in objective 4.6 of the SDGs 45 .Both indicators exhibit a negative association with the development goal, as the distance to education centres needs to be reduced to improve this objective; • Industry, innovation and infrastructure: SDG 9 has been represented using the distance to the nearest market.Objective 9.1 states the need to promote the creation of resilient infrastructures that help to develop local communities and their economy 45 .The most related indicator in the PODES dataset was the distance to the nearest market (km), which exhibits a negative association with the SDG it represents; • Life on land: for SDG 15, which addresses environmental protection, we have relied on the percentage of village area covered by forest as the SDG indicator.This feature is directly proposed in the Global indicator framework of the SDGs as an indicator for objective 15.1 "Ensure the conservation, restoration and sustainable use of terrestrial and inland freshwater ecosystems and their services, in particular, forests" 45 .The indicator is positively associated with the SDG; • Peace and Justice: for Goal 16 "Promote peaceful and inclusive societies for sustainable development" 45 , we have relied on the number of different types of crimes as the SDG indicator.This variable indicates the variety of crimes occurring in a village and was used as a proxy for criminality.Achieving world peace requires a reduction of criminality; thus our indicator is negatively associated with the SDG.
Note that the signs of the associations between the indicators and the SDGs play a major role in interpreting the positive and negative effects.Increasing one indicator does not always translate into beneficial outcomes for the SDG they represent.
The choice of SDG indicators was influenced by the methodological prerequisites of the spatial models utilised.The spatial econometric models used are designed to use continuous variables as their dependent feature.Consequently, we exclusively selected indicators of a similar nature as our SDG indicators.
Preprocessing.The dataset contains continuous, discrete, categorical and binary variables.Thus, the preprocessing differed for each type of feature.Categorical variables were transformed into binary "dummy" variables through a one-hot-encoding approach 48 .In the process, one dummy variable for each categorical feature was removed to avoid perfect colinearity between the rest of the variables.Some other features were preprocessed to remove potential colinearities.The number of core natural disasters was subtracted from the total number of natural disasters to create an independent variable describing only low-incidence disasters.Also, the "village's main source of income" was combined with the feature describing the "village's main agricultural product" to avoid the colinearity between them.Given the nonnormal distribution of some variables, we applied a log transformation to bring them closer to normality.For each feature x, the transformation logðx þ 1Þ was applied.The discarded features are detailed in the Supplementary Methods.
Finally, 1,438 villages were removed from the analysis either because of missing data or due to inconsistencies between the number of families and the number of poverty letters issued.The final dataset (after filtering) contains data for 25,067 villages in Sumatra and is described in Supplementary Table 1.
Model selection.Through the use of the Global Moran's I (MI) 49 we unveiled the presence of significant spatial clustering in all the features of our dataset and the residuals of the OLS regression, see Supplementary Table 3. Thereby, we relied on spatial regression models for the analysis.This methodology has been widely used in spatial econometrics, where it is common to observe that the geographical location of the samples can affect their characteristics [50][51][52] .
To determine the preferred spatial regression model for the analysis, we followed the pipeline described in Fig. 5.This pipeline is based on the approach described in ref. 40 , but it was expanded to evaluate the Spatial Lagged X model (SLX) and the Manski model.The pipeline starts from the OLS model and evaluates the need to include spatial dependency in the dependent variable y or in the error term ε through the use of the Lagrange multiplier (LM) 53 and the Robust Lagrange Multiplier (RLM) 54 .Both tests are based on the residuals of the OLS model and follow a chi-squared distribution with one degree of freedom 40 .Suppose the hypothesis is rejected in favour of one type of dependency (pvalue ≤ 0.05), the pipeline is continued only using spatial regression models with the preferred type of dependency, respectively the Spatial Autoregressive model (SAR) and the Spatial Error model (SEM), see ref. 19 for a formal description of these models.If both models are identified as significant, the pipeline is continued for both types of models.
The SAR and SEM models can be expanded by including spatial dependency in the independent features X, which results in the SDM and the SDEM.We relied on the Likelihood ratio (LR) to evaluate if this extension is needed.Using the LR it is possible to test if the SDM and SDEM can be simplified to their nested versions SAR and SEM, respectively.The LR was also used to evaluate if the SDM and SDEM models could be simplified into the Spatial Lagged X model (SLX) 19 , a model only considering spatial dependency over the independent features.A significant LR (p-value ≤ 0.05) indicates that the H 0 assuming that the complex model can be simplified to its nested version is rejected, thus the complex model is preferred.Contrarily if the test is not significant, the simplified version is preferred.If both models cannot be simplified they are both recommended for the sustainability analysis.
Finally, we evaluate if the SDM or SDEM models capture the whole spatial dependency in the data through the MI on the residuals of both regressions.If the residuals of the regressions are still spatially correlated, the preferred model is a model considering all types of spatial dependency, the Manski model 55 .More details regarding the application of the pipeline to our dataset are available in the Supplementary Methods.Spatial Durbin Model.The primary analysis was performed using the Spatial Durbin Model, which is defined as for a model with N spatial units and K features.The term y is an N × 1 vector describing the dependent variable for all the spatial units in the dataset.The αI N term is a constant vector of dimensions N × 1 that reflects the model's intercept.The term Xβ describes the dependency on the exogenous explanatory variables, where X is the N × K matrix describing independent features, and β is its K × 1 coefficient vector.The spatial auto-regressive term ρWy introduces spatial dependency over y by assuming that W is the N × N weight matrix describing the connectivity pattern between spatial units and ρ is the scalar parameter reflecting the strength of spatial dependence.Similarly, the term WXθ describes the spatial dependency over the independent features, and the K × 1 vector θ reflects its associated parameters.Finally, the N × 1 disturbance vector ε contains independent, normally distributed terms with a vector mean zero (N × 1), and constant variance, (σ 2 ) 19 .The SDM model is regarded as a global spatial dependency model, as the presence of the spatial lag term (ρWy) causes that changes applied in a specific location to propagate through the whole W matrix affecting the N-order neighbours recursively 19,53 .
The spatial influence between samples (villages in this case) was defined through an adjacency matrix (A ij ) created using a queen contiguity criteria 53 , in which villages are only considered neighbours if they share a common point along their borders.A ij was row-normalised to build the weight matrix (W), defining the strength of the association between villages 19 .The spatial weights matrix was decomposed into 30 power of traces using Monte Carlo Markov Chain simulations to reduce the computational burden of the simulations 19 .Both models were implemented in R 4.2.1 using the package spatialreg (version 1.2-1) 56 .Sparse Cholesky decomposition was used to compute the Jacobian of the W matrix (method 'Matrix' in the same package) 56 .
Spatial impacts.The marginal contribution of each independent feature (x k ) in spatial regression models is described by an N × N matrix (where N is the number of samples), describing the expected change in the y of all locations when x k is modified.To ease the interpretation of the results, the marginal contribution matrices can be transformed into the direct, indirect and total spatial impacts 19 .We report the average of the diagonal elements in the matrix as the direct impact, and it can be interpreted as the average change in y in a location when their own x k is increased by 1 unit.The average for the off-diagonal effects can be interpreted as the influence of a change on their neighbouring locations, and it is reported as the indirect impact (or spatial spillover effect).Finally, summing both effects we obtain the total spatial impact and it describes the average influence of modifying x k across all spatial units 19,40 .We estimated the significance of the spatial impacts numerically using the regression coefficients and the covariance matrix to perform 10,000 estimations of the impact measures 19 .For our analysis, we only considered the features with statistically significant direct, indirect and total impacts (p-value ≤ 0.05) as significant associations.
We would like to remark that the associations identified by the spatial models should not be interpreted as causal outcomes but as the spatially averaged contribution of each feature x to estimating the dependent variable y.
Prediction.In-sample prediction in spatial regression models allows estimating the expected value of y (ŷ) when the independent variables are modified in the same locations where the model was trained.There are two possible predictors for spatial regression models the "trend" and the "trend-corrected" predictors, refer to ref. 57 for more detailed definitions.In models with spatial auto-correlation (ρWy), the trend-corrected predictor provides more accurate estimations than the "trend" model.Thereby, the counterfactual analysis for the SDM model was performed using the "trend-corrected" predictor while for the SDEM model the "trend" predictor was used instead.
Quality metrics.To validate the goodness of fit of the models, we evaluated several statistical metrics.First, we report two different pseudo-R 2 metrics, a pseudo-R 2 adjusted using Ezekiel's correction 58 and the Nagelkerke's (NK) pseudo-R 259 .Pseudo-R 2 cannot be interpreted as an OLS R 2 , they can only be used as a model selection metric to compare models trained with the same dependent and independent variables.Other metrics for model selection are the Akaike and Bayesian Information Criteria (AIC and BIC respectively) 60,61 .We also evaluated the model's accuracy in describing the spatial dependency in the data.To this end, we estimated the MI 49 on the residuals (ε) of the regressions.If the spatial regression model accurately characterises the spatial dependency between samples, we expect that the residuals are spatially uncorrelated and have a non-significant (p-value > 0.05) MI 49 .Finally, we relied on the Normalised Root Mean Squared Error (NRMSE) as a goodness-of-fit metric.The prediction error was estimated using the difference between the real value for the dependent variable (y) and the in-sample predicted variable (ŷ), using the "trend-corrected" predictor for the SDM model and the "trend" estimator for the SDEM 57 .This measure was transformed Fig. 5 Schematic representation of the model selection pipeline.The pipeline is based on the approach proposed by ref. 40 , but it extends it by also evaluating the Spatially Lagged X and the Manski models.The pipeline relies on the LM, the LR and the MI to identify an adequate spatial regression model for each regression.The pipeline should be applied by starting from the OLS model and only progressing to the next model if the test meets the significance condition indicated on the arrow.The preferred model (or models) is identified when it is not possible to advance further in the pipeline.
into the Root Mean Squared Error by where N is the number of villages and ðy i À ŷi Þ is the prediction error described before 62 .The RMSE was normalised by dividing it by the range of the dependent variable y 63 , as shown in Sustainability of different production profiles.For this analysis, we aggregated the palm oil production data into four independent production profiles: active smallholder plantations, active industrial plantations, unproductive (damaged) plantations and unproductive (but possibly replanted) plantations.We considered an "active plantation" the sum of young and mature plantations for that production profile.Contrarily, the two "unproductive land" profiles were defined respectively as the sum of industrial and smallholder land for each class.Note that the four production profiles are considered independent because their sum does not reflect the totality of the village area, only the area dedicated to palm oil production.We estimated the socio-economic characteristics significantly associated with each production profile by training four SDM models with the data for each production profile as the dependent variable y and all the socio-economic variables from PODES as independent features X.The complete list of exogenous variables used can be found in the Supplementary Methods.For this analysis, the dummy variables describing the "village's main source of income" was excluded to avoid colinearity between villages reporting palm oil as primary income and villages with a high palm oil coverage.
The significant associations identified for the four regression modes were reported using the total spatial impacts.Additionally, several statistical quality metrics were reported for each model trained (refer to the Quality metrics section).
Counterfactual analysis.The counterfactual analysis was implemented following the structure described in Fig. 2. Seven regressions (SDM model) were trained using each SDG indicator as the dependent variable and the socio-economic and palm oil production features as the dependent variables (see Supplementary Methods for more details).To avoid perfect colinearity, when estimating SDG 15, variables describing village location, forest use and forest dependency were excluded from the analysis.
The counterfactual scenario proposed evaluates the influence on the SDGs of replanting all damaged plantations through a smallholder or an industrial profile.Thus, the percentage of land dedicated to unproductive (damaged) plantations (industrial + smallholder) was transferred to the young smallholder or young industrial categories, respectively.
The effect of the counterfactual was quantified by computing the difference between the dependent variable ŷ predicted after the intervention with the value predicted when the intervention is not applied.Before estimating the difference between both scenarios, the predictions obtained for each SDG indicator were bounded to reflect the characteristics of each type of variable.For all SDG indicators, the predictions were lower-bounded to zero, and only the percentage of forest coverage (indicator for SDG 15) was upper-bounded to 100%.
The effects produced by the counterfactual analysis were contextualised by determining their influence on the well-being of the villages affected.For each SDG indicator, a deprivation threshold was defined according to the values indicated by Village Development goals (Indek Pembangunan Desa) 22 .These values define the minimum levels required to consider that the village is not deprived of that socio-economic characteristic.The following list describes the thresholds selected.
• SKTM poverty letters: villages with more than 10% of the population with SKTM poverty letters are considered deprived of this indicator; • Number of healthcare facilities: villages without healthcare facilities are considered deprived of this indicator.
• Distance to primary school: villages without a primary school (distance to nearest primary school > 0 km) are considered deprived of this indicator; • Distance to junior-high school: villages with more than 3 km to the closest junior-high school are considered deprived of this indicator; • Distance to market: villages with more than 10 km to the nearest market are considered deprived of this indicator; • Percentage forest coverage: villages without forest coverage (0% of the village covered by forest) are considered deprived of this indicator.Note that the Multidimensional Poverty Index does not consider Forest Coverage as an indicator; thus, the threshold was defined as strictly as possible; • Number of different types of crime: Villages with more than 3 different types of crime are considered deprived of this indicator.These thresholds were used to estimate the deprivation status for each village and SDG before and after the intervention of the counterfactual.Interventions causing a village to change from a deprived to a non-deprived state are classified as highly positive, while interventions causing the transition from a non-deprived to a deprived status are labelled as highly negative.Contrarily, interventions that have some impact on the SDG indicator but do not induce a transition between deprivation states are only regarded as moderately positive or negative interventions, as they modify the SDG indicator, but this is not translated into changes in village well-being.Finally, in villages where replanting damaged plantations did not change the SDG indicators, the intervention is considered to have a neutral effect.
Province analysis.For this analysis, Sumatra was divided into the ten provinces shown in Supplementary Fig. 1.The counterfactual analysis was repeated using the same trained models as for the previous analysis but assuming that only the unproductive plantations in the region evaluated are being replanted.The impact on village well-being of the regional interventions was evaluated using the thresholds described in the previous section.
The results for the counterfactual analysis in ten provinces are shown in Supplementary Figs.3-12.However, we summarised them in a table describing the most beneficial and detrimental profile for each region (Table 1).The beneficial production profile for each region was estimated by selecting the interventions with the highest amount of villages with positive outcomes (highly pos.+ moderate pos.) effects for each SDG indicator.The most negative outcome was selected accordingly.The production profile with the highest magnitude (positive/negative) is reported using the following abbreviations: "I" if industrial plantations or "S" for smallholder plantations, "I/S" if the percentage of positive and negative outcomes is equivalent for both production profiles, and "-" if no outcomes are observed.

Data availability
The palm oil plantation data was obtained by combining the datasets refs. 17and 16 , which can be freely accessed via the following databases refs. 41and 42 .The geometries used to generate the maps were extracted from the following sources: Indonesia (obtained with the preprocessed socio-economic data), Malaysia 64 , Thailand 65 and Singapore 66 .A data usage license for the PODES data can be requested and purchased from the Indonesian Bureau of Statistics via https://pst.bps.go.id.Because of this, it is not possible to upload a minimal dataset on a public repository.However, we have uploaded the output from the trained regression models to the following public GitHub repository AFosch/Figs_PalmOil_SDGs.This data can be used to obtain some of the figures shown in the manuscript.

Fig. 2
Fig.2Structure of the counterfactual analysis.A complete description of the pipeline can be found in the Methods section.The SDG icons are property of the United Nations, reprinted for informational purposes in accordance with ref.67 (all rights reserved).

Fig. 3
Fig.3Contribution towards the SDGs of replanting unproductive land with smallholder or industrial plantations.The values reported indicate the percentage of villages where the implementation of the counterfactual caused high (positive/negative), moderate (positive/negative), or neutral outcomes.Note that the percentage is computed relative to the number of villages affected by the counterfactual analysis.Interventions causing the transition from a deprived to a non-deprived village are regarded as highly positive interventions.In contrast, interventions with a beneficial outcome for the SDG indicator that do not entail a change in the village's deprivation state are only regarded as moderately positive.Similarly, interventions inducing the transition from the non-deprived to the deprived state are considered highly negative, while interventions with negative influences that do not cause major changes in well-being are only considered moderately negative.Interventions where no changes in the SDG indicator are produced are regarded as neutral.
ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01037-4 10 COMMUNICATIONS EARTH & ENVIRONMENT | (2023) 4:378 | https://doi.org/10.1038/s43247-023-01037-4| www.nature.com/commsenv Significant total spatial impacts for different palm oil production profiles.The profiles evaluated are: unproductive (damaged) plantations, unproductive (but possibly replanted) plantations, active industrial palm oil and active smallholder palm oil.Positive spatial impacts indicate a positive correlation between village land and the socio-economic characteristic involved, while negative spatial impacts indicate negative correlations (see Methods for a formal description of the spatial impacts metric).Grey panels indicate not statistically significant associations (p-value > 0.05).The spatial impacts were estimated using the SDM model trained as described in the Methods section, which resulted in a Nagelkerke's pseudo − R 259 between 0.38 − 0.63 for the four models.Other quality metrics are available in Supplementary Tables5-8.
and trade of goods usually supplied to local markets, such as fruits.Finally, for SDG 16, we identified that the interventions are associated with an increase in the number of different types of crime or a decrease in crime variability depending on the village where it is implemented.