Evolutionary model discovery of causal factors behind the socio-agricultural behavior of the Ancestral Pueblo.

Agent-based modeling of artificial societies allows for the validation and analysis of human-interpretable, causal explanations of human behavior that generate society-scale phenomena. However, parameter calibration is insufficient to conduct data-driven explorations that are adequate in evaluating the importance of causal factors that constitute agent rules that match real-world individual-scale generative behaviors. We introduce evolutionary model discovery, a framework that combines genetic programming and random forest regression to evaluate the importance of a set of causal factors hypothesized to affect the individual's decision-making process. With evolutionary model discovery, we investigated the farm plot seeking behavior of the Ancestral Pueblo of the Long House Valley simulated in the Artificial Anasazi model. We evaluated the importance of causal factors unconsidered in the original model, which we hypothesized to have affected the decision-making process. Our findings, concur with other archaeological studies on the Ancestral Pueblo communities during the Pueblo II period, which indicate the existence of cross-village polities, hierarchical organization, and dependence on the viability of the agricultural niche. Contrary to the original Artificial Anasazi model, where closeness was the sole factor driving farm plot selection, selection of higher quality land, distancing from failed farm plots, and desire for social presence are found to be more important. Finally, models updated with farm selection strategies designed by incorporating these insights showed significant improvements in accuracy and robustness over the original Artificial Anasazi model.


Introduction
Exposing the mechanics of the human decision-making process that cause complex, societyscale phenomena is a difficult endeavor. Decision-making processes are often driven by multiple causal factors [1] with researchers having no direct means of measuring how these factors contribute to society-scale phenomena. Abductive reasoning via data-driven modeling and simulation techniques can overcome these issues by 'growing' artificial societies [2] and adjusting their configurations until adequate matches between simulation results and real world data are achieved. Agent-based modeling (ABM) in particular, offers the benefit of representing behaviors as human-interpretable rules. These rules are driven by the agent's autonomous evaluation of a variety of factors that are hypothesized by the modeler to be important in the decision-making process being modeled, indicated by the ability of ABM to simulate real world observations. However, a particular behavior rule only represents a single hypothetical decision-making process contained within a large space of possible, alternate decision-making processes. Exploring this vast space of rules requires the repeated re-implementation of multiple versions of the same ABM with different embedded decision-making processes [3]; this is a tedious task, involving the comparison of a massive number of combinations of causal factors. Thus, researchers often resort to modeling the most intuitive decision-making processes, if not process, which risks a subjective and inaccurate representation of the actual individual behavior [4]. The current standard of ABM exploration, parameter calibration, is a black-box technique and does not perform white-box rule exploration. Parameter calibration works on the assumption of the correctness of a predefined rule and fine tunes the coefficients of its constituent factors, but cannot easily experiment with different structures and operators through which these factors combine. Unless the importance of factors and how they are structured in the behavior rule are established, the underlying behavior rule merely remains an untested hypothesis of the actual individual behavior [2][3][4]. Parameter calibration tools are readily available for ABM frameworks, such as BehaviorSearch [5] for NetLogo [6], and OptQuest [7] for AnyLogic [8]. Inductive games have been used to infer the decision-making of societies via game theory [9], yet no established methodology exists for ABMs. As ABM rules are implemented as program instructions, genetic programming [10] is a highly suitable technique for model discovery. However, research into using genetic programming with ABMs for exploration of causality has been limited [11][12][13].
To meet this need, we introduce evolutionary model discovery, a technique for agent rule exploration and causal factor importance measurement, which combines the automated program generation capability of genetic programming [10] with the factor importance evaluation capability of random forest regression [14][15][16][17][18][19]. Unlike current standard techniques like pattern-oriented modeling [3] and model selection [20], evolutionary model discovery has the advantage of avoiding manual and repetitive re-implementation of models through automated program generation, resulting in a greatly reduced risk of implementation errors. Agent rules generated through genetic programming consist of functions of primitives that are easily comparable, as they follow a common representation. Using this representation, differences between candidate models are isolated to the code implementing the decision under scrutiny, to facilitate factor analysis and to avoid the need to compare two completely different implementations. The comparability of candidate models is important in drawing insights into the causes of the society-level phenomena being simulated. The stochasticity of genetic programming allows for the exploration of a vast space of possible agent rules, while selection of fitter models for breeding the next generation of rules ensures the exploitation of stronger factor interactions. Assumptions on agent behavior can be relaxed and rules with deeper factor interactions evaluated. Genetic programming, random forest training, and factor importance evaluation are all easily parallelizable techniques, which is important considering the large search space that can result even from a simple factor set.
In this study, we employ evolutionary model discovery to find causal factors that drive the farm plot selection decisions of the Ancestral Pueblo community modeled in the Artificial Anasazi model [21,22]. The ABM simulates the population dynamics of the Long House Valley between the years 800 AD to 1400 AD during which there was a sudden population collapse around 1350 AD. The original model demonstrated that this collapse was not caused by environmental factors alone. The model is data driven and simulations attempt to match the annual population time-series measuring households in the valley, which was estimated through data gathered from archaeological digs [21]. The agents in the model represent households, and are dependent on the agricultural success of their farm plot for sustenance and reproduction. The farm plot selection strategy originally implemented dictates that upon depletion of a household's current farm plot, the agent moves to the next closest available plot of land. In other words, the sole factor influencing this decision is the minimization of distance over the complete set of available plots of land in the valley.
We argue that this behavior does not capture the socio-agricultural decision-making of the Ancestral Pueblo, which would have most likely been influenced by factors other than distance. Archaeological and paleoecological findings have suggested that during the Pueblo II period (A.D. 890-1145) the Ancestral Pueblo organized into communities held together through cross-village polities, governed in a hierarchical and non-egalitarian manner [23]. This is evidenced through power-law and log-normal distributions of settlement sizes, intermittent power-law distribution in the Kiva (communal meeting sites, primarily used for ceremonial and ritualistic practices) sizes, and inequality in the decoration and grandeur of burial sites across the American Southwest at the time.
Evidence has shown that the periods of the Pecos classification [24], which defined the rise and fall of successful Ancestral Puebloan communities, were characterized by phases of exploration and exploitation [25]. During the exploration phases, locations suitable for farming were scouted and different organizational forms were experimented with. A successful exploration period would then result in a period of exploitation during which the discovered agricultural niches were exploited and communities held together through ritualistic practices and political force, particularly seen in the Great House system of Chaco Canyon during Pueblo II [25].
Accordingly, we hypothesize that nine different factors and four different social structures governing information flow may have driven the farm plot seeking behavior of the modeled Pueblo society. Specifically, we hypothesize that the following factors could have had significant importance: distance (F Dist ), dryness of the farm-land (F Dry ), quality of farm-land (F Qual ), yield of the land in the previous year (F Yield ), water availability (F Water ), social presence near the potential farm land (F Soc ), homophily by age (F Age ), homophily by agricultural success (F Agri ), and inter-zone migration (F Mig ), under the following possible social connectivity configurations: full information of the valley (S All ), information provided by family immediate family members (S Fam ), information provided by the most productive households (S Perf ), or information from the nearest neighbors (S Nhbr ). We consider the coefficients of these factors in the evolved agent behavior rules as the factors' 'presence' in that particular behavior rule. Each factor's presence is then analyzed for its importance at predicting the ABM's fitness through feature importance analysis on a random forest trained on data generated by the genetic program. Utilizing a random forest for this purpose allowed us to measure both main effects of the factors' presence and the joint contributions of factors towards the ABM's fitness. After identifying the most important factors, we determined the optimal presence for them. With these insights we were able to construct causally accurate and robust farm selection procedures.
Our results falsify the original assumption [21,22] that closeness was the sole causal factor governing farm plot selection of the Ancestral Pueblo society. Instead, evolutionary model discovery reports the most important factors as quality, social presence, migration from zone, distance, and dryness in order of decreasing importance. In particular, the selection of higher quality land that either had a higher social presence or was located in a different zone was shown to be more likely behavior and versions of the Artificial Anasazi with these farm selection strategies were significantly more robust against random initialization of parameters. Our results indicate that the farm selection strategy was likely more human-like than that implemented in original version of the model [21,22].

Farm plot selection in the Artificial Anasazi
The Artificial Anasazi is an agent-based model of the Kayenta Anasazi during the years of 800 AD to 1350 AD [21,22]. This model was initially developed as part of a larger effort to study the Ancestral Pueblo civilization that occupied the Long House Valley region. The ABM is implemented in NetLogo [6,22]. Archaeological excavations provide annual population time series data as estimated counts of households that existed in the valley during the period of study. The model attempts to match the historical population timeseries with its simulated household count. Annual data on water sources and estimated soil dryness (Adjusted Palmer Drought Severity Index) for each grid location on the map are provided. The model used a normal distribution to map relative quality of soil over the map. The agent-based model simulates the rise and fall of households over a geographic map of the valley over time and produces a time series of annual household count. The original purpose of the Artificial Anasazi was to test if environmental factors could have triggered the sudden disappearance of the Anasazi from the Long House Valley around 1350 AD.
Critics of the Artificial Anasazi have argued that the agent-based model itself is but a single candidate explanation of the social phenomenon at hand, the rise and fall of the Anasazi population over time [4]. However, we view this as an advantage as the Artificial Anasazi can be used as a test-bed to discover multiple plausible explanations of the population dynamics of the Long Valley at the time. Testing combinations of hypothesized factors that may have influenced actual decision-making processes of the individuals results in a vast search space of plausible Artificial Anasazi behavior results. In particular, the individual-scale behaviors modeled lack social factors and agricultural awareness, which could have generated the complex, hierarchical societies that have been shown to have existed in during the simulated period in the American Southwest, by the Village Ecodynamics Project and studies that followed [23,[25][26][27]. Societal complexity and hierarchical organization peaked during the Pueblo II (A.D. 890-1145) period and gradually faded through the Pueblo III period (A.D. 1145-1285) [23], which coincides with the period simulated by the Artificial Anasazi (A.D. 800-1350), making it crucial to consider social factors and agricultural awareness when modeling the individual-scale behaviors.
We concentrated on a particular sub-model of the Artificial Anasazi: the farm plot selection strategy. The households perform farm plot selection under two conditions: 1) when a new child household is hatched by a household that has enough resources to increase its family size, or 2) when the current farm plot is unable to produce enough yield to satisfy the nutrition needs of the household anymore. The original model, hypothesizes that the households simply selected the next closest available farm plot to the household's current farm plot during farm plot selection, i.e., minimizing over distance. A patch must be free of farms or households and not be located inside a water body to be available. Consequently, the original farm selection strategy ignores other sensory data available to the households regarding the land and the state of other households in the valley.

Hypothesized alternate factors influencing farm plot selection
Human social behavior is rarely entirely rational. Accordingly, our hypothesis proposed that the farm selection decisions of the Ancestral Pueblo were complex, and took into account the state of the potential farm plots available to them and the social influences of other households around them. Agent_Zero [2] models the human decision making process into three dimensions: social, emotional and rational. Similarly, we defined factors that we hypothesized to influence the farm plot selection process within these dimensions. The social component is expressed through four mutually exclusive social connectivity configurations through which the agent could receive information on a subset of potential farm plots, s, out of the entire set of potential farm plots in the valley, S All . The received information is then processed through a utility function f(x) defined as a combination of factors and operators, F, which consider both the internal state of the household and the conditions of the farm plot and its surroundings in order to determine the next farm plot x 0 2 s � S All as in Eq (1).
Households in the original Artificial Anasazi model consider a single factor, distance, which we will refer to as F Dist , and choose the potential farm plot with minimal distance to their current farm location. No further factors are considered in the decision making process. Furthermore, the original model assumes that the households have complete information of the valley, and every potential farm plot is compared. Therefore, the farm selection process of the original Artificial Anasazi can be represented as in Eq (1).
Arguing that the farm selection decision may have been more complex, considering a variety of other factors, we proposed an extended factor set consisting of four social and five rational factors, namely: homophily by age (F HAge ), homophily by agricultural productivity (F HAgri ), social presence (F Soc ), migration from current zone (F Mig ), comparison of quality (F Qual ), comparison of dryness (F Dry ), comparison of yield (F Yeild ), comparison of water availability (F Water ), and comparison of distance (F Dist ). Additionally, the numerical operators + and − are included in F, for the aggregation of sub-scores reported by the social/emotional and rational factors.
Four hypothesized configurations of social connectivity were included F. These configurations determined the subset of all viable farm plots that were to be considered by the households for comparison. 1) Full information (S All ): Households had complete knowledge of all potential farm plots in the valley. Full information was used by agents in the original version of the model, assuming that each household knew and compared every potential farm plot in the Long House Valley. 2) Family inherited information (S Fam ): Households solely depended on information available through their 'family'. Families are defined as a household's parent household, sibling households, any surviving grandparents, and the household itself. 3) Nearest-neighbor information (S Neigh ): agents only consider the farm plots known to their neighboring households within a fixed radius of their current location. 4) Best performers S Perf : Households only consider potential farm plots known to the best performing households, demonstrating a leadership dynamic.
Four social/emotional factors were included in F: two types of homophily (the tendency for social entities to congregate among those with similar traits), need for social presence, and one of fleeing/migration. Each social/emotional factor returned a sub-score representing the desirability of each evaluated farm plot. Sub-scores were normalized within the factors, to lie in the range of 0 to 1, for fair comparison. 1) Homophily by age (F HAge ): Households prefer to select farm plots near other households that are of similar age, where age is measured as the number of simulation steps the household has survived since splitting from its parent. 2) Homophily by agricultural productivity (F HAgri ): Households tend to select farm plots near other households with a similar corn stock to itself. 3) Social presence (F Soc ): Agents score potential farm plots with many nearby households higher than those in isolation. 4) Fleeing/ migration (F Mig ): Agents score potential farm plots that are in a completely different zone than the current one with a full sub-score, while patches in the same zone receive a sub-score of zero.
Five Rational factors considered for the farm selection process were logical comparisons of sensory data on the potential farm plots already available to the households in the original model. Similar to the social/emotional factors, rational factors also returned a normalized subscore of farm plot desirability between 0 and 1. 1) Comparison of quality (F Qual ): Higher subscores were reported for potential farm plots with higher quality of land. 2) Comparison of dryness (F Dry ): Higher sub-scores were reported for potential farm plots with higher dryness of land. 3) Comparison of yield (F Yeild ): Higher sub-scores were reported for potential farm plots that were known to have higher yield in the previous year. 4) Water availability (F Water ): Higher sub-scores were reported for potential farm plots with more nearby water sources. 5) Comparison of distance (F Dist ): Higher sub-scores were reported for potential farm plots that were closer to the current farm plot location.

Evolutionary model discovery
Evolutionary model discovery allows agent-based modelers to explore the importance of a hypothesized set of factors affecting individual-level decision making towards a macro, society-level outcome. Accordingly, evolutionary model discovery requires the modeler to identify the particular agent behavior rule being evaluated within the original agent-based model. The modeler must also provide a set of hypothesized factors and combining operators that the modeler hypothesizes to affect the decision-making process represented by the agent behavior rule.
A factor F i 2 F, where F is the modeler's set of hypothetical factors and operators, is defined as in Eq (3). Where C is the set of commands defined within F that are applied on the n number of input parameters P to produce an output return value R, where the type of each parameter t P j and the type of the return value t R are each an element of the set T of all possible parameter and return types defined by the modeler. A factor is considered an operator if C resembles an operation on one or more factors, which it accepts as parameters, rather than resembling a decision-making step. In order for a factor or operator F i to accept another F j as an input, the condition Eq (4) must be met.
An agent behavior rule b 2 B is represented as a tree of factors combined under this condition. Depending on T and the factor definitions, the space of behavior rules B can be infinitely large. To prevent the construction of such undesirably large trees, we specify a maximum depth for all b. There must be at least one F i of which t R F i is the return type expected by the entire agent behavior rule.
Given the ABM and F, evolutionary model discovery performs two stages of analysis. First, models driven by alternate decision making processes consisting of combinations of elements of F are evolved through genetic programming [10,28,29]. Genetic programming performs automated program implementation and is a suitable approach towards automating the rule discovery process [11][12][13]30]. Genetic programming evolves generations of programs through crossover and mutation operators performed on a representation consisting of primitives and terminals that combine to define program statements. Primitives are defined as a set of functions that encode program statements and may be strongly typed to only accept child and parent primitives that are compatible with the arguments and return statements accepted by its program statement. Primitives with no arguments are considered terminals. The syntax tree representation is perhaps the most common representation used in genetic programming, and arranges the primitives and terminals into a tree structure, a representation compatible with b. Programs in a generation that have a closer fit to data are more likely to be selected for reproduction through crossover and mutation to populate the next generation of programs.
Second, factor and factor interaction importance was assessed by random forest feature importance measurement. A random forest regressor was trained on the factor presence to fitness data produced by the genetic program. Random forests are an ensemble learning algorithm consisting of a forest of randomized decision trees [14][15][16]. The two most common factor importance measurement techniques for random forests are gini importance (or mean decrease in impurity), and permutation importance (or mean decrease in accuracy) [14][15][16]. However, both gini importance and permutation importance are unable to quantify the importance of factor interactions, as they consider the global importance each factor has for the random forest. Functional analysis of variance [19,31] is able to quantify the importance of factor interactions, yet lacked precision considering the inherent heteroskedasticity of the data produced by the genetic program, caused by its tendency to explore and test models of higher fitness. Instead, joint contribution [17,18] was used for this purpose as it has been successfully used to assess the importance of variable interactions in a large number of recent studies [32][33][34][35][36][37][38].
Twenty genetic programming runs were executed with the objective of minimizing the (RMSE) between the simulated household count to the actual household count over 550 simulation ticks of the Artificial Anasazi. Details on the RMSE calculation can be found in [13]. In order to ensure robustness of the evolved rules, the parameters of the ABM were randomly initialized with values ±5% about the optimal parameter values found through Stonedahl's calibration of the Artificial Anazasi through a genetic algorithm [5] [39] and parallelized by SCOOP [40]. Each genetic program run was executed for 100 generations over populations of 50 individuals. Syntax trees of minimum depth 4 and maximum depth 10 were used to avoid trees exhibiting bloat. The Half-and-Half tree builder was used for initialization [10]. To accommodate the high computational cost, the genetic program runs were distributed across a 48 vcpu Amazon Web Services EC2 instance. The random forest and gini importance algorithm of Scikit-learn [41] were used, while ELI5 [42] was used for permutation accuracy importance, and tree interpreter [18] for joint contribution measurement.
Finally, new farm selection strategies were designed taking into account the insights gained through evolutionary model discovery. The robustness of the Artificial Anasazi with these new strategies were tested against the original model by comparing the RMSE of 100 runs of each model under randomized initialization of parameters within the ranges above.

Results
The resulting best farm selection strategies evolved by the genetic program by run are provided in Table 1 along with their respective RMSE values. 15 of the runs produced RMSE values lower than the current best RMSE in the literature obtained through parameter calibration of the Artificial Anasazi model with the original farm plot selection by closeness (733.6) [43]. All best scoring rules for each run utilized S All , i.e., the model produced best results when the agents had full information regarding available farm plots as shown in Fig 1, comparing S All , S Fam , S Neigh , and S Perf over the complete factor presence to fitness data. One-tailed Mann-Whitney U tests comparing the fitness of all rules by their social connectivity configurations confirmed that rules with S All had significantly (α = 0.05) lower RMSE than the other three argmax x2S Perf f ðxÞ (p = 0.012). Accordingly, the rest of the analyses detailed in this paper were performed on rules where the social connectivity configuration was S All . Fig 2 displays the distribution of RMSE against factor presence, for presence values that were recorded in at least 200 rules across the 20 genetic program runs. Negative correlations to RMSE (higher fitness) were seen between F Dist , F Qual , F Water , F Yield , F Mig , F Soc , and F Age , and in general the genetic program favored the positive presence of these factors, and evolved more rules with these factors having a positive effect on farm selection. F Dry on the other hand had a negative correlation to RMSE for presence less than 2.
The random forest fit the factor presence to fitness data best for a forest of 520 regression trees, testing from 10 to 1000 trees with a train/test split 90%-10%. Accordingly, a forest of 520 trees was used for factor importance determination. Factor importance under S All obtained

PLOS ONE
Evolutionary model discovery of the Ancestral Pueblo through both the gini importance and permutation accuracy importance techniques can be seen in Fig 3. Gini importance generally had less precise estimations than permutation accuracy importance. Yet both techniques indicated F Qual as the factor of highest importance towards RMSE prediction. F Soc , F Mig , and F Dist also scored higher importance values than the other factors hypothesized. Fig 4 displays the p-values of one-tailed Mann Whitney U tests (alpha = 0.05), comparing the permutation importance of each factor A against every other factor B, testing the alternate hypothesis: importance of A > importance of B. According to the results, 7 of the 9 factors showed significant difference and could be ordered in terms of permutation accuracy importance as F Qual , F Soc , F Dist , F Mig , F Water , F Yield , F HAgri , F HAge , and F Dry . Fig 5 compares the top ten joint contributions towards RMSE prediction of the random forest by individual factors, and joint contributions of factors considered in pairs and triples. Again, F Qual demonstrated far higher importance than any other factor or factor interaction. The factor pairs (F Qual , F Mig ) and (F Qual , F Soc ) also demonstrated high importance, followed by (F Qual , F Mig , F Soc ), (F Dry , F Qual , F Mig ), and (F Dist , F Qual , F Soc ). Overall, F Qual was present in all highest scoring joint contributions. Despite F Dry having very low individual importance, F Dry showed higher importance when considered in combination with F Qual and F Mig .
Considering the evidence of F Qual , F Soc , F Mig , F Dist , and F Dry as important factors, Fig 6 demonstrates Mann Whitney U tests conducted for each factor F i , for the alternate hypothesis that RMSE when presence of F i was A, is less than the RMSE when presence of F i was B in rules with S All . Models with positive presence of F Qual , F Soc , F Dist , and F Mig showed significantly higher fitness (with the exception of when presence of F Mig = -2). Models with strong positive or negative presence of F Dry showed lower RMSE overall, most likely a result of F Dry 's interaction with F Qual , F Soc , or F Mig . The lowest median RMSE for (F Qual , F Soc ) was 985 at presence of F Soc at 5 and presence of F Qual at; the lowest median RMSE for (F Qual , F Mig ) was 997 at presence of F Mig at 3 and presence of F Qual at 5.
Finally, rules following the three highest joint contributions were constructed using the best values for each factor concerned: argmax x2S All ðF Qual ðxÞÞ, argmax x2S All ð5F Soc ðxÞ þ 6F Qual ðxÞÞ, and argmax x2S All ð3F Mig ðxÞ þ 5F Qual ðxÞÞ, and RMSE was compared against the original farm selection strategy argmax x2S All ðÀ F Dist ðxÞÞ for 100 runs each under random initialization of parameters within the ranges specified in section. Fig 7 shows that all three of these rules derived through evolutionary model discovery had significantly lower RMSE than that of the original farm selection strategy under randomized parameter initialization.

Discussion and conclusion
ABMs are an excellent tool for simulating and analyzing individual-scale, human-interpretable explanations of complex social phenomena. However, the design of agent rules is at the modeler's discretion and may not accurately represent the decision-making strategies of the individuals being modeled. Modelers may lack sufficient individual-scale data or observations required to identify individual-scale causal factors, and fail to provide a complete and representative design of the decision-making strategies at work. Treating the ABM as a blackbox and performing parameter calibration alone, cannot adequately explore the vast space of possible interactions of causal factors and operators that may combine to form more realistic representations of the actual decision-making processes of interest. Instead, a whitebox exploration of the agents' behavior rules must be performed, which treats the ABM as a sandbox upon which different combinations of hypothesized causal factors and operators are tested in order to identify important factors and their interactions. We address this issue with the introduction of evolutionary model discovery, which is able to distinguish, out of a hypothesized set of causal factors, those causal factors that are important towards the generation of the behavior of interest and their role in the decision-making process. By combining automated program generation via genetic programming with feature importance evaluation via random forests, evolutionary model discovery is able to quantify the importance and optimal presence of these factors in the decision-making process that result in society-level phenomena simulated by the ABM. This allows for the construction of agent rules that more accurately represent the real-world decision-making process of individuals and result in more robust models. In addition to discovering the Ancestral Puebloan socio-agricultural behaviors in this study, evolutionary model discovery has been successful in exposing causal factors of several agent-based models of complex social phenomena; these include factors driving residential

PLOS ONE
Evolutionary model discovery of the Ancestral Pueblo segregation and integration, and factors governing the prioritization of response under information overload on online social media [45].
Evolutionary model discovery provided several advantages in the discovery of individualscale behaviors in this study. Automation of rule design and model implementation through genetic programming eliminated the likelihood of programmatic error under manual reimplementation and made it possible to search the vast space of possible behavior rules. Being an optimization algorithm, the genetic program identified and searched the space of relevant behaviors thoroughly. The grammatical representation used ensured that hypothesized factors were specified as multiple reusable, human-interpretable units of causality, which were then easily analyzed. Both the genetic program and random forest were easily parallelized and run on high-performance cloud computing environments. Finally, the random forest was computationally inexpensive, with the time complexity of training increasing linearly with the number of unique hypothesized factors evaluated [45]. There were some limitations to the use of evolutionary model discovery in this study. It was assumed that the entire population of agents embodied the same behavior rule. This made it unable to distinguish between archetypes of individuals that might have embodied competing strategies. Additionally, evolutionary model discovery did not perform a systematic selection of features by their ability to fit multiple patterns organized in a hierarchical ordering of complexity, as recommended in pattern-oriented modeling [3]. Finally, the ability to reproduce 2D spatial patterns of the households was not tested as is typical of pattern-oriented modeling, which may have provided more insight into formation of meso-scale aggregations and expose preferred locations of the simulated space.
Applying evolutionary model discovery on the Artificial Anasazi, we show that the socioagricultural behavior of the Ancestral Puebloans of the Long House Valley was more deliberative and informed than assumed in the original Artificial Anasazi model. Instead, these results concur with several other archaeological and paleoecological findings in the literature regarding the Ancenstral Puebloan societies of the American Southwest during the Pueblo II period. Models designed through evolutionary model discovery insights are significantly more robust. Comparison between the RMSE of 100 runs of three models with farm selection strategies designed taking into consideration the insights from evolutionary model discovery, 1) argmax x2S All ðF Qual ðxÞÞ, 2) argmax x2S All ð5F Soc ðxÞ þ 6F Qual ðxÞÞ, and 3) argmax x2S All ð3F Mig ðxÞ þ 5F Qual ðxÞÞ, against 100 runs of the original farm selection strategy argmax x2S All ðÀ F Dist ðxÞÞ in [21,22,43,44], under random initialization of parameters. The three farm selection strategies derived from evolutionary model discovery were far more robust under random parameter initialization and showed significantly better RMSE scores compared to the original model. Bocinsky et al., highlight the importance of the viability of the agricultural niche for the growth of rain-fed maize dependent Ancestral Pueblo communities [27]. Similarly, our results indicate that, contrary to the original farm selection behavior of the Artificial Anasazi model, where households would select the next closest possible plot of land once their present farm was depleted, the households most likely selected potential farming land with higher soil quality (F Qual ), an indicator of the viability of the agricultural niche in the area. Furthermore, instead of choosing farm plots closer to the current, failed plot (−F Dist ), choosing farm plots that were farther away from the household's current, failed farm plot (F Dist ) or moving to a completely different zone in the region (F Mig ) was found to be a more likely behavior. These results agree towards a significant degree of awareness and consideration towards the viability of the agricultural niche.
Bocinsky et al. discover that phases of exploration preceded the exploitation phase [25]. The exploration phase would thus ensure that there was considerable knowledge of the region and resource distribution among the communities when deciding their next farming location. Interestingly, our results agree that it was highly likely that the households had near-complete knowledge of the potential arable land throughout the valley, since S All was the best social connectivity configuration for information spread.
Archaeological studies show that the periods of exploitation were likely characterized by cross-village polities, social organization, ritualistic practice, and political force [23,25]. Similarly, our results indicate the desire to congregate into communities was strong, as positive desire for social presence (F Soc ) was the second most important factor, and acting on information on arable land known to neighboring households (S Neigh ) was the second most successful social connectivity configuration.
Overall, versions of the Artificial Anasazi where farm plot selection was driven by seeking either higher quality land, higher quality land with more social presence, or higher quality land in different zones, all proved to be significantly more robust than the decision to move to the next closest available plot of land (Fig 7). An interesting extension of this work would be to incorporate the decision-making strategies discovered in this paper into the agent-based model of the Village Ecodynamics Project's agent-based model, in order to test the generalizability of the behaviors found on the data of the Long House Valley used in the Artificial Anasazi model to other Puebloan communities at the time.
Supporting information S1 File. Evolutionary model discovery Artificial Anasazi. This archive contains the Evolutionary Model Discovery Python source code. This Python package is also being actively maintained at: https://github.com/chathika/evolutionarymodeldiscovery and documentation is available at https://evolutionarymodeldiscovery.readthedocs.io/en/latest/. (ZIP) S1 Table. Factor scores. Factor presence to model fitness data produced by the 20 genetic programming runs on the Artificial Anasazi.