Forests, Coca, and Conflict: Grass Frontier Dynamics and Deforestation in the Amazon-Andes

Population growth with weak economic development can promote tropical deforestation, but government infrastructure investment can also open new frontiers and thus increase deforestation. In the Andean region of South America, population growth has been a leading explanation for both deforestation and coca cultivation, but coca generates armed conflict and attracts counter-drug measures, obscuring the differences between population-driven and frontier-opening models of deforestation. Using a 15-year panel from Colombia, we model deforestation, coca cultivation, and conflict victims as interrelated responses with a suite of covariates encompassing land cover, land cover changes, population, population changes, counter-drug measures, and government infrastructure spending. Infrastructure spending suppresses coca, coca and eradication by aerial fumigation both increase conflict, and conflict promotes deforestation and is associated with depopulation. But the strongest predictor of deforestation is pasture growth, which covaries with coca. While these models show that infrastructure spending can help reduce coca, and coca’s influence on deforestation is indirect and mediated by conflict, the models also reveal the most important challenge to forest conservation is neither coca nor conflict, but an insatiable appetite for land that expresses itself through pasture growth.


Introduction
Tropical forests hold vast stores of carbon and the greatest biodiversity in the world, which makes their conservation essential to stabilize climate and prevent extinction (Myers et al. 2000;Mitchard 2018). Although the specifics of carbon dynamics in tropical forests are debated (Baccini et al. 2017;Hansen et al. 2019;Hubau et al. 2020), agriculture, forestry, and other land use contribute ~2 3% of annual global emissions, mainly through deforestation (Arneth et al. 2019). For example, while transportation is responsible for ~1 4% of annual emissions (Edenhofer et al. 2014), policy and technical advances promise to reduce those emissions (Santos 2017). In contrast, there are no obvious technological solutions for land use emissions and mitigating climate change within sustainable limits requires the maintenance or restoration of tropical forests, particularly in the Amazon (Soares-Filho et al. 2010, Chazdon et al. 2016, IPCC 2019). Yet, the Amazon-along with all tropical forests of Latin America-is threatened by infrastructure development in general (Nobre et al. 2016;Bebbington et al. 2020), and road construction in particular (Perz et al. 2008;Gallice et al. 2017). Justified by the need to grow their national economies through mining, ranching, and commodity production (Richards & VanWey 2015, Sonter et al. 2017, Maeda et al. 2021, infrastructure expansion requires massive investment, transforming both physical and sociodemographic landscapes (Nobre et al. 2016). At the same time, maintaining ecosystem services such as carbon sequestration and water cycling (Lovejoy & Nobre 2019;Staal et al. 2020), and averting catastrophic biodiversity loss (Gomes et al. 2019) requires curbing tropical deforestation. Thus, there is tension between forest conservation for local and global ecosystem services and economic development, and this makes discovering the factors influencing tropical forest transformation essential to avert global catastrophe.
Two models, immiseration and frontier, have been proposed to explain tropical forest loss (Rudel & Roper 1997). While the immiseration model refers to poverty-driven deforestation from growing populations seeking to expand subsistence agriculture, the frontier model involves powerful actors opening new frontiers by promoting economic development in forests (Fearnside & de Alencastro Graça 2006;Barber et al. 2014;Ferrante & Fearnside 2020). Growing populations are central to both models of deforestation, but the drivers of deforestation differ in each case. While poverty and rural population growth are fundamental in the immiseration model, investment and urban populations are crucial in the frontier model. With high levels of inequality (Armenteras et al. 2019a), large populations living in poverty, and vast forests, immiseration has been a traditional explanation for forest loss in Latin America (Barbier 1997;Carr 2006). But in the Amazon, at least, most deforestation accrues from large landholdings incompatible with immiseration (Fearnside 1993). Over the last two decades, globalization has made demand for commodities and infrastructure development-the latter more compatible with the frontier model-salient drivers of deforestation throughout Latin America (Grau & Aide 2008;DeFries et al. 2010). Both models, then, can potentially explain tropical deforestation in the region.
A particular version of the immiseration model has shaped deforestation analyses in Colombia and other Amazon-Andes countries. Coca cultivation, grown by migrant farmers for the global cocaine market, has featured prominently as a driver of land use change in Colombia, Peru, and Bolivia (Álvarez 2002;Bradley and Millington 2008b;Young 1996). Since coca growers are, on average, poorer than other farmers in subnational analyses (Davalos 2016;Dávalos & Dávalos 2020), poverty and immiseration seem to be logical explanations for coca deforestation. Yet, even studies identifying coca cultivation as an important driver of landscape change fail to fully support the immiseration model because legal land uses such as pastures dominate (e.g., a ratio of >80:1 pasture: coca or 495,917 ha of pasture to 6,013 ha for 2010 in San Lucas [Chadid et al. 2015]). At the municipality scale with relevant demographic data, the unique effects of coca cultivation vanish (Dávalos et al. 2011), and frontier dynamics from mechanized agriculture and infrastructure construction often dwarf coca-driven deforestation, for example, in Bolivia (Killeen et al. 2007). Nevertheless, because violent armed conflict is a central feature of illegal coca cultivation (Díaz & Sánchez 2004), conflicting estimates of demographic and cultivation effects may result from omitting armed conflict in deforestation analyses at the municipality level (Negret et al. 2019).
With this paper, our goal is to simultaneously model the rate of forest change (positive for gain, negative for loss or deforestation), coca cultivation, and conflict as likely interrelated responses to changes in human activities. We argue that the frontier model of deforestation better explains forest loss in Colombia and in particular its Amazon-Andes. In addition to testing for the role of state-sponsored economic development on deforestation, a long history of incomplete, abandoned, or otherwise failed development projects have left behind a legacy of coca and conflict that requires modeling the latter two interrelated factors. Here we test predictions from the immiseration and frontier hypotheses by analyzing deforestation, coca cultivation, and armed conflict simultaneously with data from Colombia. In support of the immiseration model, population growth and increases in coca cultivation (and other crops) should be important covariates of deforestation, while infrastructure development should be more important in the frontier model. Although multivariate causal or quasi-causal methods have recently been implemented to analyze deforestation at the municipality scale (Fergusson et al. 2014, Christiansen et al. 2020Mendoza 2020), none tested infrastructure spending as a proxy for state-led frontier-opening actions. Since our data, municipalities across multiple years, correlate in space and through time, we must also factor in spatiotemporal trends for testing these hypotheses. While frontier dynamics have been proposed to explain the ubiquity of coca cultivation at deforestation fronts (Dávalos et al. 2011;Dávalos et al. 2016), our analyses are the first to integrate forests, coca, conflict, and development into a set of interdependent spatiotemporal models at a subnational spatial resolution.

Tropical forests, coca cultivation, and conflict in Colombia
In Colombia, tropical forests are concentrated in the Amazon-Andes and the Chocó biogeographic region. While these two regions encompass just 11% of continental municipalities, they comprise 54% of the total area in the country (see Table 1). These are also two of the most biodiverse places in the world (Myers et al. 2000;Dávalos et al. 2011), holding a complement of endemic species disproportionate to their size. Despite their outsized conservation value, these forests, along with highly fragmented Andean and Caribbean remnants are currently threatened by conversion to agriculture in general and pastures for cattle ranching in particular (Dávalos et al. 2014), illegal mining (Anaya et al. 2020), coca cultivation (Chadid et al. 2015, and infrastructure development-albeit apparently unsanctioned by the state (Silva Numa 2016; Paz Cardona 2021).
Although Colombia's illicit crop production includes marijuana, opium poppy, and coca, the latter is by far the most extensive and disperse of the three crops (Dávalos & Dávalos 2020). Coca cultivation in Colombia is also concentrated in the Amazon-Andes and Chocó biogeographic region and has often been considered a poverty-driven cause of deforestation (ecogeographic regions are shown in Figure 1). About 95% of the municipalities in these environmentally rich regions have some coca cultivation, in contrast with only 37% in the rest of the country. Together, the two regions concentrate 23% of all municipalities with coca cultivation in the country, or about double what would be expected based only on the number of political units they hold. While coca has become a national challenge since the 2000s (Rincón-Ruiz & Kallis 2013), the Amazon-Andes have been a coca cultivation stronghold for much longer than anywhere else in Colombia (UNODC 2010;Dávalos et al. 2016). Close spatial correspondence between deforestation and coca cultivation at the landscape (pixel) scale suggests this is an important factor in deforestation (Negret et al. 2019). But this relationship disappears when examined at the municipality scale in which sociodemographic factors can be modeled (Dávalos et al. 2011;Armenteras et al. 2013a;Sánchez-Cuervo & Aide 2013). Nevertheless, across South America, coca cultivation is a key factor to understand deforestation, if not directly, then indirectly through either associations with other land uses such as pastures (UNODC 2006;Bradley and Millington 2008a;Armenteras et al. 2013b;Chadid et al. 2015), links to armed conflict and displacement (Dion & Russler 2008;Ballvé 2012), or because the counter-coca response of forced eradication can displace crops to generate new deforestation frontiers (Dávalos et al. 2009;Rincón-Ruiz & Kallis 2013). Therefore, it is also important to include counter-drug measures, even though eradicating coca by spraying herbicide from airplanes, or aerial fumigation, is not as important as urbanization or road improvement in explaining coca declines in parts of the Amazon-Andes (Dávalos et al. 2014).
Across South America's Andean region, coca cultivation has been explained as the result of state weakness (Shifter & Jawahar 2004), lack of territorial control (Thoumi 2005), and poor market integration (Dion & Russler 2008). While acknowledging enduring challenges to integrate economies across the coca frontier (Dávalos & Dávalos 2020), recent scholarship upends the weak state model by instead highlighting the long history of state-led development interventions across Amazon-Andes, and the role of those projects in opening new agricultural frontiers that attracted colonists, and later became centers of coca production (Dávalos et al. 2016). Instead of state absence, in this new framework, coca cultivation results from a series of state-led policy and investment decisions with the goal of economically integrating tropical lowlands into national economies, particularly in the Amazon (Gootenberg & Dávalos 2018). Then, when national governments all but withdrew their support for these projects in the late 1970s and 1980s, an illegal agricultural economy took hold (Dávalos 2018). In Colombia, centers of coca cultivation in the 1990s, before the massive expansion of the aerial fumigation program and subsequent dispersal of coca throughout the country (Dávalos et al. 2009, Rincón-Ruiz & Kallis 2013, overlap almost perfectly with former Amazon-Andes colonization projects in the Ariari region of Meta (Torres 2018), Guaviare, Caquetá, and Putumayo (Dávalos et al. 2016). In the wake of failed state-led efforts to instigate economic development in peripheral areas (Gootenberg 2020), eradication and alternative development became just another set of state-building tools that often involved great violence (van Dun 2012; Rincón-Ruiz & Kallis 2013). Although our time series is too short to test this alternative explanation for the emergence of coca, this history is relevant because state interventions, aerial fumigation and manual eradication (Davalos 2016), are part of the general counter-drug strategy that may influence deforestation. In Colombia, armed conflict overlaps with both tropical forests and coca cultivation. On average, municipalities with coca crops have almost four times more human right violations in comparison to municipalities without coca cultivation, and the number of human right violations in municipalities with coca in the Amazon-Andes and Chocó biogeographic region surpasses the national average. Conflict is also disproportionately associated with extensive forest cover (Holmes et al. 2018;Reilly & Parra-Peñas 2019), as well as particular economic activities associated with certain land uses, such as cattle ranching with pastures (Holmes et al. 2019). The effects of armed conflict on forest conservation and biodiversity has been studied for decades (Dudley et al. 2002;Fjeldså et al. 2005;Hanson et al. 2009;Gaynor et al. 2016), but quantitative estimates focused on Colombia have found conflicting scale-and region-specific effects (Sánchez-Cuervo & Aide 2013). For example, while paramilitary presence boosted growth in woody cover in the Cauca Valley montane region, this effect was reversed in northern Andean montane forests, with no trend for guerrilla presence (Sánchez-Cuervo & Aide 2013). In Urabá, part of the Chocó biogeographic region, land cover/land use analyses revealed violence as a key factor in forest loss, replacing forests (Fergusson et al. 2014). More recent analyses found a detectable but small contribution of conflict to deforestation (Negret et al. 2019). To summarize, tropical deforestation, coca cultivation, and armed conflict are interrelated and should be analyzed simultaneously while factoring in suites of land covers, changes, state interventions against coca production and for development, and sociodemographic variables.

Data
Our analyses include five types of data: 1) new spatial data on land cover from Google Earth Engine (GEE) (Gorelick et al. 2017), 2) data on illicit crops from the Illicit Crops Monitoring Global Program of the United Nations Office of Drugs and Crime (UNODC), 3) population data from the National Administrative Department of Statistics (DANE for its Spanish-language acronym) (DANE 2012), 4) conflict data from the Banco de Datos de Derechos Humanos, DIH y Violencia Política (CINEP for its Spanish-language acronym) (CINEP 2020), and 5) infrastructure spending data from the Colombian National Planning Department (DNP for its Spanish language acronym) (DNP 2019). While the first of these data sets was available at a 25-hectare (ha) resolution, the remaining four were available at the scale of municipalities. In Colombia, there are 1,122 municipalities grouped into 32 departments and one capital district, however, we excluded one Caribbean Island department (San Andrés y Providencia) and filtered the data. Therefore, the new land use spatial data was summarized at the municipality scale, and analyses were conducted on 916 municipalities filtered as outlined below. Table 2 summarizes the panel data (i.e., observations of the same municipalities over time) included in analyses.

Land cover data excluding coca
We used Google Earth Engine (GEE) (Gorelick et al. 2017) to access and collect base information from the MCD12Q1 MODIS yearly land cover type (Sulla-Menashe & Friedl 2019) for the 2001-2015 period.
The data set has annual information and a spatial resolution of 500-m. We used the 'Cover Type 1' band from this product, which refers to the classification of land cover types and uses the IGBP (International Geosphere-Biosphere Program) methodology, reclassed it into four types: forests (Type 1), encompasses all forest covers with tree canopies >2 m and tree cover <60%; grasslands (Type 2), areas dominated by herbaceous annuals (<2 m); croplands (Type 3), encompasses areas where at least 60% of it is cultivated and also mosaics of small-scale cultivation 40-60%; and others (Type 4), includes mostly all shrublands, savannas, urban, water bodies, snow and barren. Natural grasslands dominated by herbaceous annuals (<2 m) in this region of the world correspond to páramos (natural savannas such as those of the Llanos are under category 4). Páramos are neither ubiquitous nor expanding rapidly. Therefore, we interpret the expansion of 'grasslands' as evidence of pasture-based land uses and hereafter use 'pasture' to designate change in this land cover over time. The spatial resolution of these data (25-ha pixel) is coarser than that of coca cultivation, which is collected by identifying 1-ha pixel (and then aggregated for municipalities), meaning croplands of any type under 25-ha extension could not be detected using these new spatial data.
Although the licit or illicit use of land in the cropland variable is unknown, most coca cultivation is not expected to be detected using the cropland data because over 60% of coca is grown in fields under 1 ha (UNODC and Gobierno de Colombia 2014). This is reflected in our data, which yields a negligible fraction (<4%) from dividing coca cultivation over croplands for 80% of observations. Thus, coca cultivation and croplands were not collinear (R = −0.01, t 13727 = −1.64, P-value = 0.101) and therefore were unlikely to bias model results. To match the data sets, we calculated the area in hectares of each cover type per year per municipality. We obtained land cover data for every year (2001-2018) and every municipality ( Table 2). While we focused on variables available for most years for 2000-2017 (18,972 observations, 1,116 municipalities), using lags to calculate growth rates reduced the data by one year (to 17,856 observations). We then ignored 184 municipalities completely lacking forest and grasslands at any point in time, resulting in data for 15 years from 2001-2015 (14,912 observations, 932 municipalities). Finally, infrastructure spending data was missing from 2017 and from some municipalities (13,729 observations, 916 municipalities).

Coca, fumigation, and manual eradication data
To quantify coca cultivation, we used the net area under coca cultivation at the cut-off date of the Annual Coca Survey at the municipal scale (December 31). Each year, the Illicit Crop Monitoring Global Program of the United Nations Office of Drugs and Crime (UNODC) captures satellite images that cover the entire continental Colombian territory (1,142,000 Km 2 ). The accuracy in identifying coca fields from the satellite images ranges between 87% and 90% (UNDCP, 2002;UNODC, 2003). After the images are captured, the UNODC conducts field verification to calculate the extension of the area under coca cultivation with gaps or covered by clouds. Although there is measurement error in coca quantification, estimates are still unbiased because poor weather and satellite location are external and random with respect to coca grower decision making (Davalos 2016). Importantly, the coca cultivation area identified in the images is also adjusted for aerial and manual eradication activities performed during the same period. After the corrections for gaps and clouds and adjustments for eradication activities, the result is the net area under coca cultivation annually for each municipality. Manual and aerial eradication are also reported in annual hectares by municipality, obtained from the Colombian Antinarcotics Police (DIRAN for its Spanish-language acronym) (DIRAN 2019). We modeled coca count as a response variable in our analyses.
Based on land cover data (forests, grasslands, croplands, and coca cultivation), we calculated dynamic growth variables (r land use ) per municipality using equation 6 of Puyravaud (2003) for rates of change. To avoid division by 0 for cases in which no land cover had been observed the previous year, units equivalent to a single pixel observation (i.e., 1 ha to coca, 25 ha to other land uses) were added to the corresponding land cover variables before calculating rates. All dynamic growth variables were negative to indicate loss and positive to indicate gain, with 0 indicating stability. An important focus of these analyses is to determine how growth in forest and grassland cover types relate to each other. To avoid including municipalities lacking either land cover over the study period, we filtered the data to include only those observations with both forest and grass for any year sampled (or >25 ha summed for both covers). This yielded 916 municipalities available for analyses. We modeled the rate of change in forest cover for each municipality as a response variable in our analyses.

Population data
To include population, population types (rural, urban), and population change as variables, we used data from DANE (DANE 2012). To quantify changes in population we computed the rate of change-population growth or r population -from one year to the next for total, urban, and total rural population. These rates are positive for population increases and negative for loss. We also included the rural fraction of the population, hereafter 'rural proportion' as a covariate.

Conflict data
CINEP collects municipal data on the number of victims of all kinds of human rights violations perpetrated by different armed groups, such as guerrilla, paramilitary, or the national army. 1 Although these data are usually collapsed into rates (e.g., homicides per 100,000 inhabitants), even when rare, violence can affect land use (Murillo-Sandoval et al. 2021). In addition, most municipalities most of the time lack conflict victims, with a few developing temporary clusters of conflict. We therefore modeled rates as a response in initial maximum likelihood analyses but used counts in Bayesian analyses. In the latter, the count can be modeled as a negative binomial or zero-enriched negative binomial while also accounting for population size by including population variables as covariates.

Infrastructure data
Infrastructure spending data refers to expenditures on fixed capital, such as: land, roads, buildings, and equipment. Expenditures are measured in nominal thousand pesos from the DNP (DNP 2019). Annually, each municipality reports to the DNP their income and expenditures desegregated by activity. To standardize these measures and make them comparable across municipalities and over time, we used per capita values and the consumer price index (CPI) to adjust for inflation.

Empirical model
Our analyses can be summarized in three steps. First, we used maximum likelihood regressions to determine which covariates were significantly associated with each of three responses: change in forest cover, coca cultivation, and conflict. Second, since these methods accounted for neither spatial and temporal correlations, nor the non-normal frequency distribution of the responses, we used Bayesian methods to estimate the magnitude and direction of relationships with covariates. These models included smooths to capture variation in the responses arising from a municipality's location in space and time so that the covariates explained the variance remaining in the observations. Lastly, we used posterior predictive checks-simulating data using the modeled equations-to assess whether response frequency distributions and spatiotemporal trends adequately captured the variance observed in the data.

Piecewise structural equation modeling (piecewise SEM)
Although Structural Equation Models (SEM) have been applied in the social sciences for decades (Tarka 2018), these approaches have only recently been deployed in attempts to understand the relationships between armed conflict and forest loss in Colombia (Christiansen et al. 2020), without controlling for other covariates known to influence forest loss. Here we use SEM, or the multivariate relationships modeled as the sum of direct and indirect relationships among variables, to test predictions from the immiseration and frontier hypotheses relating land cover and human activities. In a departure from prior work (e.g., Armenteras et al. 2011;Dávalos et al. 2011;Sánchez-Cuervo et al. 2012;Armenteras et al. 2013a), we seek to model the multivariate spatiotemporal behavior of forest cover change, coca cultivation, and conflict, each as a function of demographic and other land use covariates. Each response corresponds to one equation in a set of three, with responses and covariates all observed by municipality by year (except for municipality area, which is constant over time).
Combining demographic, land use, and anti-coca measures generated a suite of 15 potential covariates (listed in Table 2), not counting interrelationships among the three response variables: forest cover change (hereafter r forest ), coca cultivation, and conflict. Both the temporal and spatial nature of the observations pose substantial challenges to standard univariate analysis, more so in the case of traditional SEM, which assumes independence among observations (Lefcheck 2016). This assumption of independence is particularly relevant when evaluating the goodness of fit of the estimated to the observed covariance matrix (Lefcheck 2016) that enables rejecting or confirming the entire model (Grace et al. 2012). To overcome these assumptions especially with ecological data, directed acyclic or piecewise SEMs have been developed (Shipley 2000(Shipley , 2009. Instead of estimating directional relationships or paths globally, paths are estimated within individual models, then pieced together (Shipley 2009). In evaluating a piecewise SEM, the key criterion is no longer a global goodness of fit, but the test of d-separation (Shipley 2000). This statistic evaluates whether paths are missing from the model, and how the model would be improved by including missing paths. A significant d-separation test indicates the model is not a good fit, and several missing paths contain information useful to explain the observations.
To explore the universe of relationships between each of the three responses -r forest coca cultivation, and conflict -and potential covariates, we used piecewise SEM as implemented in the piecewise R package (Lefcheck 2016). Both fitting Poisson distributions to count data (coca and conflict) and modeling spatial autocorrelation proved computationally prohibitive. Hence, to generate responses that might approximate normal distributions, both coca and conflict variables were log 10 -transformed (+1 for coca and +3 for conflict rates per thousand inhabitants). To account for spatiotemporal patterns, we fitted a series of linear hierarchical models with municipality-and year-specific effects and with year as a covariate using the nlme R package (Pinheiro et al. 2012). The initial set of equations included all covariates for the r forest response and a single covariate (year) for the coca cultivation and conflict responses. Relationships among the three responses were evaluated by adding coca cultivation or conflict as covariates of r forest and of one another until no significant P-values in d-separation tests remained for responses as covariates. Then, covariates with significant P-values in d-separation tests were incorporated iteratively until there were no more significant P-values in d-separation tests at all. Although both the municipality-and year-specific effects are expected to account for some of the spatiotemporal structure of the data, they do not fully address it, and count data need to be modeled using distributions for counts such as Poisson or negative binomial. For these reasons, we used the results of piecewise as starting points to fit Bayesian models.

Bayesian structural equation modeling
To overcome the strong spatial autocorrelation and temporal patterns in the observations, we implemented piecewise SEMs in a hierarchical Bayesian framework using the R package brms (Bürkner 2017) which implements models in stan (Carpenter et al. 2017). While both traditional and piecewise SEM rely on maximum likelihood to estimate both model parameters and the Akaike information criterion (Akaike 1974) for model selection (Shipley 2013), stan estimates posterior parameter distributions resulting from the likelihood and (uninformative) priors for said parameters. By incorporating complex data structures as part of the hierarchy of observations, brms/stan provided the flexibility needed to overcome the spatial and temporal correlations among observations. Although spatial autoregression models are implemented in brms (e.g., Morris et al. (2019)), spatiotemporal models are not. Instead, we modeled spatiotemporal effects using Markov random fields (MRF) smoothing on municipalities, with repeated measures by year. Specifically, we took advantage of the brms implementation of thin-plate regression smooths via mgcv (Wood 2003). By using the neighborhood matrix of each discrete unit (i.e., municipality), this method models the spatial data as an intrinsic Gaussian Markov random field (GMRF) that is then passed onto the brms model as a municipality-specific effect. While the GMRF defines a spatial effect that varies smoothly over the unit's neighbors, by itself, it does not include smoothing by year, which were added with each year as a factor. Thus, the combination of GMRF spatial and factor year smooths aims to capture the spatiotemporal structure of the data through sets of group-specific effects defined by sample-wide smooth means ('sample-wide smoothing' or µ g ) and groupspecific terms described by the corresponding standard deviation ('group-specific smooth sigma' or σ g ). Such that the smooth for each municipality m, at year t, in group g was given by: In which three groups of parameters g defined at least six and no more than nine sets of municipalities with similar smooths for each response variable. These smooths were included in the models by using the t2 (or smoothing) argument with year and municipality as variables, and options bs = c("fs", "mrf") and xt = list(year = NULL, MPIOS = list(nb = nb)). While the first option describes how to build the smooth as factor interaction (fs) for year and as a Markov random field (mrf) for municipalities, the second passes information about the variable structure such that year has no structure beyond being sequential, but municipalities (MPIOS) are structured by geographic vicinity given in the neighbor list 'nb'. The smooth terms were included in each model formula for every response variable.
With both negative and positive values and a preponderance of values centered on zero (no change), r forest rates calculated following Puyravaud (2003) can be approximated by normal distributions. At the scale of analysis, however, r forest centered on zero but exhibited an excess of zero values, generating a narrowly kurtotic or leptokurtotic frequency distribution. This was coupled with "fat tails," or a higher prevalence of extreme values than expected for a normal distribution. To address these mismatches between the model and the data, we reformulated r forest as a Student's t distribution, in which ν ('nu') is a positive integer shape parameter that determines both the kurtosis and the frequency of extreme values, such that an observation y m,t at municipality m at year t was modeled thus: In which µ = 0 for ν >1, and 2 2      , with ν values between 1 and 2 yielding σ 2 = ∞.
In which smooth g,m,t are municipality-and year-specific effects estimated in Equation β 0 ; is a sample-wide intercept; β 1 is a vector of estimated coefficients; and X m,t comprises forest, grassland and cropland cover, r pasture , conflict victims, infrastructure spending (lagged), and municipality area (as in ' covariate' of Table 2, all but r pasture in log10 scale).
In contrast with r forest , coca cultivation and conflict victims are both counts and strongly non normal even when transformed by logarithm after adding a constant (as was done to expedite exploratory piecewise analyses). Despite the prevalence of coca and conflict throughout Colombia, most municipalities lack both, and those that do have either coca or conflict, generate a positive long-tailed distribution. To better approximate the frequency distribution of both responses, these were analyzed as negative binomial or zero-inflated negative binomial variables. In this case each of the two responses were modeled as either: In which r is the shape parameter that describes the mean of a mixture of Poisson (count) distributions (i.e., r = gamma-Poisson shape parameter), and 1 pr pr  defines a rate parameter for the gamma distribution that defines the mixture. This allows for relaxing the expectation of equality of mean and variance of the Poisson distribution and is generally a better fit to observational data (O'Hara & Kotze 2010). The zero inflated negative binomial formulation adds one more parameter corresponding to the zero-inflation probability in the data that is not already modeled by the negative binomial distribution. Thus, the count variables were modeled as: In which smooth g,m,t are municipality-and year-specific effects estimated in Equation 1; β 0 is a sample-wide intercept; β 1 is a vector of estimated coefficients; and X m,t comprises forest, grassland and cropland cover, r pasture , r cropland , infrastructure spending (lagged), rural proportion, r urban population , r coca , aerial fumigation (lagged), manual eradication (lagged) and area for coca cultivation, and forest and cropland cover, infrastructure spending (lagged), total rural population, rural proportion, r urban population , r coca , coca cultivation, aerial fumigation (lagged), and manual eradication (lagged) for conflict victims (as in ' covariate' of Table 2, all but r variables and proportion given in log10 scale).
Although brms implements models in the stan statistical language, which uses a Hamiltonian Monte Carlo sampler algorithm more efficient than traditional Bayesian Gibbs sampling (e.g., Plummer 2003), complex analyses still require sampling separate Bayesian chains-or samplers starting from separate starting points-for many iterations. These models ran for 2,000 iterations, across four chains, with 800 iterations as burn-in, and sampling posteriors every five iterations across four simultaneous chains. The potential scale reduction factor (PSRF) (Gelman & Rubin 1992), which approaches one at convergence among model posterior distributions across chains, was used to determine if chains had converged. All parameters were summarized only after posterior estimated sampling sizes exceeded 200.

Model comparisons and visualizations
We used approximate leave-one-out ('loo') cross-validation statistics (Vehtari et al. 2017) in posterior predictive checks to both compare models and evaluate model fit. To compare the two model sets, we calculated leave-one-out probability integral transformations (loo-pit) on the posterior likelihood of the models. If the model is well calibrated, resulting loo-pit values are asymptotically uniform (Gabry et al. 2019), and this comparison can be visualized in quantile-quantile plots comparing data sets simulated from a standard uniform distribution against the model's loo-pit values. We also examined the posterior estimates of the zero-inflation probability to determine if this more complex model was necessary. After determining which model was best calibrated, we used the R package bayesplot (Gabry 2017) to visualize posteriors and MCMCvis (Youngflesh 2018) to summarize model parameters. Relationships between selected covariates and responses were extracted using the conditional_effects routine in brms.
Although SEMs rely on interpreting sample-wide coefficients, group-specific parameters based on the rich spatial and temporal data analyzed here are also informative. Group-specific Gaussian Markov random fields smooths display spatiotemporal trends in the data. We extracted these using native routines in brms. Since these smooths varied in both space and time, we plotted four-time steps to illustrate changes in trends.

Results
Analyses using piecewise SEM found suites of variables explaining variance across responses (Table S1). Standardized coefficients for r forest as dependent variable identified strong positive associations for forest cover, and negative relationships with municipality area and r pasture . For coca cultivation, the strongest associations were positive for municipality area, r coca , and aerial fumigation. For conflict, the greatest associations were negative for year and proportion of the population that is rural (hereafter, 'rural proportion'), then positive for both total rural population and forest cover. After five iterations, final d-separation tests estimated by Fisher's C 28 = 23.993, P = 0.682 indicated no paths were missing from these models and the models could not be improved with the suite of variables analyzed. Since these analyses did not completely factor the spatiotemporal structure of the data, we recoded the final piecewise models into multivariate equations in brms.
We ran two sets of Bayesian analyses corresponding to the negative binomial and zero-inflated negative binomial while accounting for the spatiotemporal structure of the data. All Bayesian analyses converged with PSRF<1.05, and large effective sampling sizes. Despite convergent chains and large effective sampling size of the posterior, initial runs had a few divergent transitions in posterior sampling which corresponded to biased sampling of parameter space. To avoid these divergent transitions, we increased the target average proposal acceptance probability (adapt_delta parameter) from the default of 0.8 to 0.9. Although this slowed down sampling, the change eliminated all posterior divergent transitions.
Comparisons between models using loo-pit statistics revealed the negative binomial and zero-inflated negative binomial models to be similar ( Figure S1). Both models show an acceptable fit for r forest , and poor fit for extreme values for both count variables. Zero-inflation probability estimates, however, were low, estimated at 95% high probability density intervals of 0.002-0.01 for coca cultivation, and 0.001-0.007 for conflict victims. Therefore, we summarize and discuss results for the simpler negative binomial models. Figure  2 summarizes most of the sample-wide parameters, and group specific parameters are shown in Table S2. Relationships inferred using piecewise were confirmed for the most part, but Bayesian coefficient estimates for six predictor-response combinations overlapped with zero (i.e., significant piecewise coefficients were not corroborated once spatiotemporal dynamics were fully modeled), and there were a few sign changes. For r forest , the most influential variables were negative for r pasture , forest cover-a sign change compared to piecewise-and conflict, and positive for grassland cover. For coca cultivation, while-in order-coefficients for r urban population , municipality area, rural proportion, forest cover, manual eradication, aerial fumigation, r pasture , and r coca were positive, grassland cover and infrastructure spending had negative effects. Total rural population, coca cultivation, forest cover and aerial fumigation were positive influences on conflict victims, with r urban population and rural proportion having negative effects.
To illustrate relationships among responses, Figure 3 shows how conflict victims nonlinearly decrease r forest and coca cultivation nonlinearly increases conflict victims (i.e., as coca goes up, conflict increases, as conflict increases r forest decreases). Effects from two predictors-forest cover and year-with effects across the three responses are shown in Figure 4. Although r forest nonlinearly decreased with forest cover, the latter nonlinearly increased both coca cultivation and conflict victims. Year had complex effects as a predictor: while r forest increased over time but remained negative (i.e., forest loss) for the average municipality, coca cultivation increased through 2004 and with high variance in estimates, then decreased, and conflict victims declined with low variance, except for a secondary and low peak in 2011. Figure 5 summarizes three predictors of r forest and coca cultivation: grassland cover which has a positive nonlinear effect on r forest and a negative nonlinear effect on coca cultivation; r pasture , whose effect is linear and negative on r forest but nonlinear positive for coca cultivation; and infrastructure spending, surprisingly finding no effect of the latter on r forest and a negative nonlinear effect on coca cultivation. While there is more coca cultivation in municipalities with less grassland area, coca cultivation is also associated with a higher r pasture . Figure 6 summarizes predictors common to coca cultivation and conflict. With only one exception, covariates had opposing nonlinear effects on coca cultivation and conflict: while cropland cover had a negative effect on coca cultivation, rural proportion, r urban population , aerial fumigation and manual eradication (whose relationship with conflict overlaps with zero), all had positive nonlinear effects on coca cultivation. In contrast, aerial fumigation nonlinearly increased both coca and conflict. As each set of models defined relationships among three responses with multiple covariates, Figure 7 summarizes relationships among variables based on sample-wide coefficients. Modeled GMRF smooths showed distinctive patterns for each response, with the strongest trend in the decline of conflict over time ( Figure S2).

Discussion
Through a series of multivariate models applying both structural equations and spatiotemporal smooths on data for Colombia for the first time, we estimated complex interactions among r forest , coca cultivation, and conflict that fit uneasily in both the immiseration and frontier models. While the immiseration model predicts a negative relationship between r forest and population, no such direct link was detected and the strongest negative influence on r forest was instead r pasture . Diverging from the frontier model as well, infrastructure spending had no effect on r forest . Yet, as the only development-related factor, infrastructure spending had a strong negative effect on coca cultivation, illustrating how development investment helps suppress illicit coca production. Both static and dynamic population variables had opposing relationships with coca cultivation and conflict victims, even as the latter two responses were positively associated. Despite including only infrastructure spending and population as sociodemographic predictors, our results have implications for understanding land cover changes, particularly from forest to grassy vegetation. Throughout Latin America, the transformation of forest to grass is a ubiquitous land cover change that signals the human-induced conversion of forests to pastures, especially at the forest frontier (Graesser et al. 2015;Armenteras et al. 2017) and is often mediated by fire (Armenteras et al. 2021). Our results thus show quantitative evidence for dynamic replacement of forest with grass cover for pastureland use, and how this change relates directly or indirectly to coca production, counter-drug policy, and conflict.    Neither r urban population nor infrastructure spending related to r forest , rejecting predictions from both the immiseration and frontier theories of deforestation. Nevertheless, finding r pasture and conflict as strong covariates of deforestation points to frontier dynamics largely independent of government investment. Three lines of evidence support r pasture as a signature of the moving frontier. First, if immiseration were the main mechanism of deforestation, smallholder crops-including coca cultivation-would dominate forest transformation. Instead, past analyses for Colombia have shown pastures to be the primary land cover replacing forests from the northern Andes in San Lucas (Chadid et al. 2015), to the south in Guaviare (Dávalos et al. 2014), and the Amazon more broadly (Armenteras et al. 2013a), with our study corroborating prior findings for the whole country. Second, this effect holds despite controlling for forest and grassland cover, but larger forests have lower r forest (Figure 4), indicating advances into larger forest blocs of the kind that require some form of investment, and in line with the frontier model (Rudel & Roper 1997). Lastly, there is evidence that as pastures replace forests, private financial transactions also increase in the road-mediated frontier of Guaviare (Dávalos et al. 2014). Therefore, r pasture might also be an indicator of the incorporation of forests-which are nominally public-into private hands (as it has been in Brazil [Hecht 1993;Fearnside & de Alencastro Graça 2006]), and is by itself a form of spontaneous economic development (White et al. 2000).
Armed conflict is a tool of territorial control that instigates both forest conversion and forced depopulation. We found more conflict victims resulted in more forest loss (Figure 3), with conflict concentrating in municipalities with greater forest (Figure 4) and cropland cover (Figure 6), but seemingly unrelated to infrastructure spending (Figure 2). While linking conflict to forest loss corroborates prior studies relating conflict to deforestation (Castro-Nunez et al. 2017;Negret et al. 2019), it does so while modeling conflict as a response with covariates of its own. We therefore confirm there is more forest cover in municipalities experiencing more conflict (Holmes et al. 2018;Reilly & Parra-Peñas 2019), which also connects conflict to the more forested deforestation frontier (Figures 2 and 4). Yet, the strongest covariates of conflict victims were demographic, with steep declines in urban population indicated by r urban population (Figure 6) and larger total rural populations (Figure 2). Past analyses have shown armed conflict leads to forest loss (Murillo- Sandoval et al. 2021), as it likely displaces local populations (Fergusson et al. 2014) and correlates to both cattle ranching as a land use (Holmes et al. 2019) and land grabbing (Castro-Nunez et al. 2017). Here we show a nationwide effect of conflict victims on r forest (Figure 3) while demonstrating that armed conflict and displacement (as measured by r urban population ) go together (Figure 6) (Holmes & Gutiérrez De Piñeres 2011;Rincón-Ruiz & Kallis 2013). Lastly, this conflict frontier is not driven by government spending (Figure 2), indicating that infrastructure investment could be an effective component of anti-drug policy without instigating conflict.
In a break with previous scholarship, we consistently identify both coca cultivation (Figure 3) and aerial fumigation (Figure 6) as positive covariates of conflict. Previous analyses have shown that coca promotes violence (Angrist & Kugler 2008), specifically by paramilitary (Holmes et al. 2019) and guerrilla actors vying for territorial control (Holmes et al. 2018), and armed conflict increases coca cultivation (Díaz & Sánchez 2004). Similarly, coca eradication by aerial fumigation has been shown to covary with armed conflict (Holmes et al. 2006, Holmes et al. 2018, although time series for specific years show the opposite relationship (Fisher & Meitus 2017). Yet, past analyses linked conflict with either coca cultivation or aerial fumigation, but not both. Importantly, aerial fumigation (but not manual eradication) was a positive covariate of conflict ( Figure  6), making this tool of anti-drug policy both a direct driver of conflict and an indirect driver deforestation via conflict victims. Past analyses proposed aerial fumigation drives deforestation by displacing coca cultivation through the balloon effect (Rincón-Ruiz & Kallis 2013). Unlike Davalos and Morales (2019) or Rincón-Ruiz & Kallis (2013), we did not set out to test the balloon effect, but we did include a measure of change in coca which did not covary with conflict: r coca (Figure 2). After controlling for the spatiotemporal structure of the data, we find evidence that aerial fumigation exacerbates conflict, much more than r coca , likely signaling displacement-mediated deforestation.
Discovering that conflict promotes deforestation nationwide while confirming that aerial fumigation increases conflict poses a conundrum: how to reduce coca cultivation without exacerbating conflict, which instigates forest loss? Infrastructure spending, already known to reduce farmer propensity to grow coca (Dávalos & Dávalos 2020), was neither related to r forest nor to conflict (Figure 2), but was negatively related to coca cultivation (Figure 5), suggesting it is feasible to reduce coca cultivation through social and infrastructure investment (Davalos 2016) without unduly exacerbating deforestation. While finding no strong positive relationship between infrastructure spending and deforestation is surprising (e.g., these are tightly linked in Brazil [Fearnside & de Alencastro Graça 2006;Barber et al. 2014;Ferrante & Fearnside 2020]), this can be explained in two complementary ways. First, infrastructure of all kinds is precarious at the deforestation frontier, with decades-old old settlements still lacking basic services in both San Lucas (Dávalos & Dávalos 2020) and the Ariari in the Amazon-Andes (Torres 2018). As a result, official roads-one of the main manifestations of infrastructure spending-fail to relate to deforestation, as in San Lucas (Chadid et al. 2015), or along the border with Ecuador in Putumayo (Viña et al. 2004). Second, private and illegal infrastructure with no environmental oversight sometimes fills these voids, as in the Macarena-Guaviare region (Silva Numa 2016), or in San Lucas itself (Chadid et al. 2015). Thus, government infrastructure spending fails to represent access roads that promote deforestation on the ground, but it does reflect some level of social investment that reduces coca cultivation (Davalos 2016; Dávalos & Dávalos 2020). Based on our findings, more infrastructure spending can reduce coca cultivation, perhaps enabling greater control through manual eradication which, unlike aerial fumigation, does not promote conflict.
In conclusion, by integrating deforestation, coca cultivation, and conflict into sets of interrelated models, our analyses have, for the first time, found quantitative, nationwide support for frontier deforestation dynamics driven by conversion to pasture and armed conflict, both of which are associated with coca cultivation. As in other Amazonian countries, forest loss involves replacement with pastures-and seldom crops-on a vast scale, whether domestic or international demand for cattle warrant it or not (Hecht 1993;Laurance et al. 2002;Rudel et al. 2002;Fearnside 2005;McAlpine et al. 2009). But unlike neighboring countries' frontiers, Colombia's are not instigated by government development investment, or not to a level detectable with these data. Instead, both coca cultivation and its largest-scale countermeasure, aerial fumigation, promote armed conflict that, in turn, generates deforestation and is associated with displacement. Because government infrastructure spending reduces coca cultivation while promoting neither deforestation nor conflict, targeted development investment could be a powerful tool for countering coca. But such interventions must consider both the history of failed state-sponsored development projects, and that private, illegal infrastructure has already led to deforestation spikes in the same regions where coca and conflict concentrate (Silva Numa 2016;Paz Cardona 2021). Right now the status quo ante-both coca cultivation and aerial fumigation generating thousands of conflict victims each year and promoting depopulation and deforestation-has given way to massive environmental costs at newly opened post-peace accord frontiers (Armenteras et al. 2019b;Van Dexter & Visseren-Hamakers 2019;Clerici et al. 2020;Murillo Sandoval et al. 2020) where lingering armed conflict generates much greater deforestation than before (Murillo- Sandoval et al. 2021). While our models show that infrastructure spending can help reduce coca, and coca promotes deforestation via conflict, they also reveal the most important challenge to forest conservation in Colombia is neither coca nor conflict, but the insatiable appetite for land that expresses itself through r pasture . We therefore expect this transformation of forests into pastures to continue along an endless grass frontier unabated by decreases in conflict or fluctuations in coca cultivation.

Additional File
The additional file for this article can be found as follows: • Supplementary File. Supplementary Data and Results. DOI: https://doi.org/10.31389/jied.87.s1

Data Accessibility Statement
The data and all scripts necessary to generate our results are available through GitHub (https://github.com/ lmdavalos/endless).