Food storage facilitates professional religious specialization in hunter–gatherer societies

Professional religious specialists centralised religious authority in early human societies and represented some of the earliest instances of formalised social leadership. These individuals played a central role in the emergence of organised religion and transitions to more stratified human societies. Evolutionary theories highlight a range of environmental, economic and social factors that are potentially causally related to the emergence of professional religious specialists in human history. There remains little consensus over the relative importance of these factors and whether professional religious specialists were the outcome or driver of increased socio-cultural complexity. We built a global dataset of hunter-gatherer societies and developed a novel method of exploratory phylogenetic path analysis. This enabled us to systematically identify the factors associated with the emergence of professional religious specialists and infer the directionality of causal dependencies. We find that environmental predictability, environmental richness, pathogen load, social leadership and food storage systems are all correlated with professional religious specialists. However, only food storage is directly related to the emergence of professional religious specialists. Our findings are most consistent with the claim that the early stages of organised religion were the outcome rather than driver of increased socio-economic complexity.


Hypotheses
We undertook a literature search to identify hypotheses that attempt to explain the presence of religious specialists in historical small-scale societies. The literature we identified focuses on a wide range of different types of religious specialists and practitioners. For example, Singh (Singh, 2017) focuses on shamans, who are defined as "practitioners who enter trance to provide services" and are described to be a restricted professional class, and who are viewed as having greater abilities or powers to interact with invisible forces and diagnose or treat problems in life. By contrast, Winkleman (M. Winkelman, 1989) focuses on "magicoreligious practitioners," who includes a broad range of healers, shamans and priests/priestesses that "had access to supernatural power". This means that the hypotheses we have identified in the literature may not always have been intended to apply specifically to professional religious specialists, as defined in this project. As such, we treat these hypotheses as a starting point for identifying the kinds of factors that could plausibly be causally related to the emergence of professional religious specialists in human history.
For previous research to contribute a hypothesis to our study, it needed to: (a) clearly identify one or more variables related to religious specialists in early human societies; (b) the proposed relationship(s) needed to plausibly apply across different world regions, and; (c) the variables identified needed to be operationalizable with available ethnographic source materials and databases. When multiple sources identified similar factors, and made similar predictions, this literature was collapsed into a single hypothesis. Table S1 provides a list of the 15 hypotheses we identified, along with selected quotes of text related to each hypothesis.

Variables
Below we outline the definitions used to code each variable used in this study. The full data and code used in this study is available through our OSF project page (https://osf.io/6jdqr/).
The phylopath package assumes that variables are either continuous, or binary coded. For this reason, we dummy coded all categorical and ordinal variables. We also scaled numeric variables to have a standard deviation of one, and centered values to have a mean of zero.

Professional Religious Specialists
For an individual to be defined as a professional religious specialists they must: (a) be considered to have a special ability, expertise or authority in relation to the supernatural, and (b) receive material payments for their services. When members of communities are perceived to differ in supernatural abilities, expertise or authority, but do not receive material payments for their services, they were not considered professional religious specialists. For example, among the Pumé of South America, there was no "material gain to flow from attaining the position of shaman," so professional religious specialists were coded as absent (Petrullo, 1939). In cases where offerings were ostensibly made to gods or spirits, but clearly provided material benefits to religious specialists, these offerings were considered to be material payments to religious specialists. When there were both professional and nonprofessional religious specialists within a given society, professional religious specialists were coded as present. This variable was coded as part of our Hunter-gatherer Religion Database (Watts et al., 2019).

Social Class
Social class refers to the presence of social class distinctions within a society, including both wealth based distinctions and the stratification of the population into heritable social classes.
Social class provides a general measure of the presence of social inequality within societies.
This variable is a binary recoding of Binford's (Binford, Lewis, 2001) "CLASS" variable where 0 indicates social classes are absent, and 1 indicates that social classes are present.
There are no missing data for this variable.

Slavery
Slavery refers to the state of slaves within a society, including debt slaves, other nonhereditary slaves, and hereditary slaves. This variable is a binary recoding of Binford's (Binford, Lewis, 2001) "SLAVE" variable where 0 indicates that slavery absent and 1 indicates that slavery is present. There are no missing data for this variable.

Internal Conflict
Internal conflict refers to the frequency of internal conflict within a society. This refers to conflicts that occur within the largest consumer groups that regularly camp together within a year. This variable is a binary coding of Binford's (Binford, Lewis, 2001) INTCON variable where 0 indicates that internal conflict is rare or absent and 1 indicates that internal conflict is moderate or frequent. There are two societies with missing data for this variable.

Intergroup Conflict
Intergroup conflict refers to the frequency of violence between groups that do not regularly or periodically camp together. Intergroup conflict includes raids conducted on groups of the same cultural unit that do no camp together. This variable is a binary coding of Binford's (Binford, Lewis, 2001) GPGPCON variables where 0 indicates no violence reported, and 1 indicates that inter-group conflict was present. There is one society with missing data for this variable.

Warfare
Warfare refers to how widespread and frequent warfare is for a society. War can include periodic attacks on intruders as well as warfare undertaken to gain new territories. This variable is a binary recoding of Binford's (Binford, Lewis, 2001) WAR1 variable where 0 indicates there is no organized warfare, and 1 indicates that organized warfare is present.
There are no missing data for this variable.

Food Storage
Food storage refers to the quantity of food stored, including food use for special events and during low productive phases of the year. This variable is a dummy coding of Binford's QTSTOR variable, which originally had four levels; 1 = minor quantities of food stored; 2 = moderate quantitates of food stored; 3 = major quantities of food stored; and 4= massive quantities of food stored. We binary coded this variable so that 0 indicates minor food storage, and 1 indicates moderate, major or massive food storage. There are no missing data for this variable.

Social Leadership
Social leadership refers to the presence of community leaders that are widely recognized, visible, and whose performance is evaluated by others. This does not include leadership based only on age-based distinctions, or domain specific differences in expertise. This variable is based on Binford's (Binford, Lewis, 2001) LEADER variable where 0 indicates that social leadership is absent, and 1 indicates that social leadership is present. There is no missing data for this variable.

Number of Annual Moves
Number of annual moves refers to the average number of times that household units shift locations per year. This variable is based on Binford's scalar NOMOV variable, which has value ranges from 0.1 to 45. For the purposes of our analyses we scale and center the NOMOV variable. There are 14 societies with missing data for this variable.

Distance of Annual Moves
Distance of annual moves refers to the average total distance that households move per year.
This variable is based on Binford's DISMOV variable, which was originally in mile units and had values ranging from 8 to 570 miles (mean = 183.9, median = 90.0). The distribution of values for this variable was left skewed so we log10 transformed, scaled and normalized these values for the purposes of analyses. There are 19 societies with missing data for this variable.

Consumer Group Size
Consumer group size refers to the mean number of people that regularly camp together during the most aggregated season of the year. This variable is based on Binford's GROUP2 variable, which originally had values ranging from 20.2 to 330.0 (mean = 70.47, median = 46.50). The distribution of values for this variable was left skewed so we log10 transformed, scaled and normalized these values for the purposes of analyses. There are four societies with missing data for this variable.

Max Camp Size
Max camp size refers to the maximum number of people that camp together for subsistence purposes. This is a mean value and does not necessarily occur annually. This variable is based on Binford's GROUP3 scalar variable, which originally had values ranging from 62 to skewed so we log10 transformed, scaled and normalized these values for the purposes of analyses. There are 18 societies with missing data for this variable.

Total Population Size
Total population size refers to the total number of people within a cultural unit. These people do not necessarily interact with one another and this group size is based upon the classification boundaries of ethnographers. This variable is based on Binford's TLPOP scalar variable, which originally had values ranging from 70 to 13500 (mean = 1919.20, median = 1200.00). The distribution of values for this variable was left skewed so we log10 transformed, scaled and normalized these values for the purposes of analyses. There are no missing data for this variable.

Pathogen Load
Our variable on pathogen load is based on an established variables that Fincher and Thornhill ) generated from the Global Infectious Disease and Epidemiology Network (GIDEAN) (Yu & Edberg, 2005). This variable, known as "pathogen" in the "Binford" R package, refers to how common seven different classes of pathogens (totaling 22 different parasites) are within the region that the society is located. These seven classes of parasites were: leishmanias; trypanosomes; malaria; schistosomes; filariae; spirochetes; and leprosy. Each of these seven classes was coded on a three point scale by GIDEN or Fincher and Thornhill , on which 3 means that the pathogen is endemic, 2 means it is sporadic, and 1 means that it is neither endemic nor sporadic. The sum of these seven categories was then taken as a measure of overall parasite prevalence. The data used to calculate these measures of pathogen load are based on observations gathered from April-June 2007, and have been shown to be highly correlated (r = 0.77) with measures of historical pathogen prevalence . We centered and scaled pathogen scale scores for the purposes of analyses. There are no missing data for this variable.

Environmental Richness
Environmental richness refers to the mean grams of carbon uptake per square meter of land per month (gC/m 2 /month). The rate of carbon uptake is a general measure of the ecological productivity of an environment, and is known as Net Primary Productivity (NPP). We use environmental data collected by NASA's Moderate Resolution Imaging Spectroradiometer (Running et al., 1999), as processes and stored by the D-Place database of world cultures (Kirby et al., 2016). Further information about how mean NPP was calculated for each society is available on the D-Place website and in Kirby et al. (Kirby et al., 2016). There are no missing data for this variable.

Environmental Predictability
Environmental predictability refers to year-to-year variation in the cyclical patterns of NPP, as defined by Colwell's (Colwell, 1974) information theoretic index. An environment can be predictable if is highly stable (high constancy), or if there are highly regular seasonal changes (low contingency). For an environment to be unpredictable, it must have both low constancy and high contingency. We use D-Place's measure of NPP predictability (Kirby et al., 2016), which is calculated using data from NASA's Moderate Resolution Imaging Spectroradiometer (Running et al., 1999) and has been used in previous research on the cultural evolution of religion (Botero et al., 2014). Further information about how NPP predictability was calculated for each society is available on the D-Place website and in Kirby et al. (Kirby et al., 2016). There are no missing data for this variable.

Environmental Variation
Environmental variation refers to variation in NPP per month, as calculated and defined in the D-Place database (Kirby et al., 2016). Environmental variation is distinct from environmental predictability because environments can have large but predictable variation in NPP, such as those that occur over seasons, or large and irregular variation in NPP, such as intermittent droughts. The NPP values used here are calculated using data from NASA's Moderate Resolution Imaging Spectroradiometer (Running et al., 1999). There are no missing data for this variable.

Language Tree
We built a tree representing the historical relationships between the languages spoken in each society. This global tree was constructed using language classification data from the Glottolog catalogue of world languages (Hammarström et al., 2020). We assigned Glottocodes to societies in our sample using the codes assigned to Binford's sample in the D-Place database of world cultures (Kirby et al., 2016). Glottolog contains a comprehensive index of classification data for over 20,000 world languages and dialects. Glottolog classifications are provided in a nested structure based on linguists' consensus judgements of the genealogical relationships between languages. This language classification provides the underlying structure of language phylogenies, but does not contain information on the scale of the branches. We developed a method for simplifying the classification scheme of Glottolog and scaling branches according to different potential distances between taxonomic classifications.
Our method of tree building restricts language classifications to four taxonomic levels; language family, language sub-family, language and dialect. Language family refers to the broadest taxonomic grouping within Glottolog. Language isolates are each treated as being of their own distinct language family. There are a total of 37 language families included in our sample, including the Austronesian, Dravidian and Uto-Aztcan language families. We define a language sub-family as the highest order taxonomic classification within a language family that differentiate at least 10% of the documented languages within that language family. We follow Glottolog classification codes for languages and dialects.
Two societies in our sample, Cholanaikkan and Nayaka, are classified as speaking the Jennu Kurumba language. For the purpose of representation on our language phylogeny, we treat these two societies as speaking different dialects of the same language.
We built three trees, each representing a different approach to scaling the distances between societies within our phylogenetically adjusted models ( Figure S3). These three trees provide a proxy for cultural ancestry within phylogenetic based statistical methods and we replicate analyses across all three to test the sensitivity of our results to branch scaling assumptions. We emphasize that the scaling of branches on these three trees is not proportional to the timing of language diversifications in human history. Across all trees, the root of a tree is a polytomy with each language family diverging. All language isolates are treated as representing different language families. In the first tree, the phylogenetic distance between languages increases linearly with taxonomic distance. The length of branches representing the language family splits from the root of the tree have a length of one, branches representing the splitting of language subfamilies from families have a length of one, and branches representing the splitting of languages from subfamilies have a length of one. Branches representing the distance between dialects and languages was set to 0.5, and the branch representing the length from the parent language to the language sub-family was scaled to 0.5 in order to retain ultrametricity. In this tree, two dialects of the same language have a distance of one, two languages of the same language sub-family have a distance of two, two languages from different sub-families of the same language family have a distance of four, and two languages from different language families have a distance of six.
In the second tree (Tree 2 in Figure S3), the phylogenetic distance between languages progressively increases according to the taxonomic distance between languages. The length of branches representing the language family splits from the root of the tree have a length of three, branches representing the splitting of language subfamilies from families have a length of two, and branched representing the splitting of languages from subfamilies have a length of one. Branches representing the distance between dialects and languages was set to 0.5, and the branch representing the length from the parent language to the language sub-family was scaled to 0.5 in order to retain ultrametricity. In this tree, two dialects of the same language have a distance of one, two languages of the same language sub-family have a distance of two, two languages from different sub-families of the same language family have a distance of six, and two languages from different language families have a distance of twelve.
In the third tree (Tree 3 in Figure S3), the phylogenetic distance between languages increases sub-linearly with the taxonomic distance between languages. The length of branches representing the language family splits from the root of the tree have a length of one, branches representing the splitting of language subfamilies from families have a length of two, and branched representing the splitting of languages from subfamilies have a length of three. Branches representing the distance between dialects and languages was set to 1.5, and the branch representing the length from the parent language to the language sub-family was scaled to 1.5 in order to retain ultrametricity. In this tree, two dialects of the same language have a distance of three, two languages of the same language sub-family have a distance of six, two languages from different sub-families of the same language family have a distance of ten and two languages from different language families have a distance of twelve.
This tree building approach results in a series of three simplified, ultrametric phylogenies that provide a proxy for the ancestral relationship between languages. In the main manuscript, we report only the results using the first tree, discuss the results of our replicated analyses in the Supplementary Results section, and provide the outputs on the OSF project page (https://osf.io/6jdqr/). Our findings remain largely unchanged across these three trees, supporting the generality of our findings.

Exploratory phylogenetic path analysis
Standard path analyses are unable to account for the non-independence of societies, so are subject to Galton's Problem (Mace & Pagel, 1994). Phylogenetic comparative research methods have the ability to help address Galton's Problem by incorporating a family tree representing the historical relationships between cultures (Mace & Pagel, 1994). The method of confirmatory path analysis developed by Hardenberg and Gonzalez-Voyer (Hardenberg & Gonzalez-Voyer, 2013), and implemented in the R package phylopath (van der Bijl, 2018), provides a framework for comparing different causal path structures while adjusting for the phylogenetic non-independence of societies. This method provides a way to model and compare a defined set of path models representing the causal relationships while addressing Galton's Problem.
We build on existing the confirmatory phylogenetic path analysis method by developing a package called BEPA (Brute-force Exploratory Path Analysis). The BEPA package includes functions to builds an index of all possible causal path structures for a set of variables, and to convert models from formulae to directed acyclic graphs. The approach taken to indexing all possible causal paths is computationally intensive and the number of possible models increases exponentially with each additional variable included. In order to improve performance, this package includes an option argument to exclude models that contain certain kinds of path relationships. For example, users can choose to generate an index of all possible models, except those in which Social Leadership predicts environmental predictability. This, and the multi-threading argument, aid in efficiency but we warn users that the package remains computationally intensive. The code used to perform BEPA is publicly available through GitHub (https://github.com/Joseph-Watts/BEPA).
In the main article we summarize the best fitting models from our analyses. The best fitting models are defined as those models within 2 CICc of the best fitting model (36, 43).
Averaging across the best fitting models accounts for uncertainty in model selection and conditional averaging calculates average path coefficients only when they are included in a model (rather than treating absent paths as zero) (36, 43).

Phylogenetic signal
We tested for phylogenetic patterning of professional religious specialists on our language phylogeny by calculating the Fritz & Purvis (Fritz & Purvis, 2010) D statistic using the phylo.d function from the R package named caper (Orme et al., 2013). The number of permutations was set to 100,000.

Phylogenetic and geographic distances
We examined regional clustering in professional religious specialists and generated a matrix representing the geographic distances between societies. The geographic distribution of professional religious specialists suggests differences in the frequency of professional religious specialists across world regions (Figure 2). Of the 26 hunter-gatherer societies located in North America, 25 (96%) were found to have professional religious specialists, whereas 2 of the 6 (33%) hunter-gatherer societies in Australia had professional religious specialists.
Using a Mantel test, we found evidence of significant geographic clustering in the distribution of professional religious specialists (z = 446.75, p = 0.001). Mantel tests (z = 19.00, p = 0.001) and Fritz and Purvis D (D = -0.33, p = 0.001) also indicate that there is phylogenetic signal in the distribution of professional religious specialists. These tests suggest that there is geographic and spatial autocorrelation in the distribution of professional religious specialists that must be adjusted for in comparative analyses (Dow & Eff, 2008;Jordan, 2013). We note that the geographic and phylogenetic distances between societies are themselves significantly correlated (z = 112.15, p < 0.001). Here we use phylogenetic comparative methods to adjust for the non-independence of societies, treating phylogenetic distance between societies as a general proxy for the non-independence of societies in our sample.

Non-Phylogenetic Models
In addition to the phylogenetically adjusted models reported in the main article, we performed a series of standard Generalised Linear Models (GLMs). These models tested the bivariate relationship between predictor variables and our variable professional religious specialists, without adjusting for the common ancestry of societies using a phylogeny. GLM analyses were performed using the GLM function in R (R Core Team, 2015), with a binomial distribution.
Of the 16 variables analyzed using GLMs, seven were found to be significantly correlated with Professional Religious Specialists in hunter-gatherer societies (Table S2). By contrast, only five of these relationships remained statistically significant after adjusting for common ancestry within our phylogenetically adjusted GLMs (Table 1). The results of our GLMs are included to show the importance of adjusting for common ancestry when modelling the cross-cultural distribution of professional religious specialists. We do not discuss the results of the GLMs further because this statistical method assumes an independent sample, which the results of our phylogenetic analysis do not support (see the section below on Tree Construction).

Phylogenetically-adjusted Bivariate Relationships
As reported in our main analyses, Phylogenetic Generalized Linear Models (PGLMs) were used to identify the cultural and environmental factors associated with professional religious specialists (Table 1). Phylogenetically adjusted bivariate analyses were performed using the phyloglm function in the phylolm R package (Ho et al., 2016). For our PGLM analyses we used logistic MLPE as the method of estimation and set the tolerance of the fitted values to 30. These analyses do not account for potential co-variance between predictor variables, but do adjust for dependence due to common cultural ancestry. In order to understand how adjusting for the non-independence of societies affects our results, we also include standard non-phylogenetic bivariate analyses in our Supplementary Methods (Table S2).
The results of our bivariate PGLM analyses are reported in the main article. One limitation of these bivariate analyses is that they do not assess the directionality of the relationship between variables, so cannot always distinguish between competing hypotheses.
For example, the positive association between social leadership and professional religious specialists is compatible with the hypothesis that professional religious specialists facilitate social leadership (H.7), as well as the hypothesis that social leadership facilitates professional religious specialists (H.8). Another limitation of these analyses is that it is not possible to tell whether the effects of one variable on another is partially of fully mediated by another. To address these limitations, we performed phylogenetically adjusted path analyses (see the main article and the exploratory phylogenetic path analysis section of the Supplementary Methods).

Follow-up analyses
We performed a series of follow-up analyses to test how robust the results of our exploratory phylogenetic path analyses (main article) are to variations in tree structure and model constraints.
First, we replicated our analyses using the second and third trees described in the Language Tree section of the Supplementary Methods (Tree 2 and Tree 3 in the Figure S3).
In our first replication, using Tree 2 ( Figure S3), the results of our phylogenetically adjusted bivariate analyses find the same variables are significantly associated with professional religious specialists, and the direction of their relationship remains the same. For the exploratory phylogenetic path analyses using Tree 2, there are 19 best fitting models (within 2 CICc of the lowest CICc). Across all 19 of these best fitting causal models, the only variable that is significantly and directly causally linked to professional religious specialists is food storage. In all 19 best fitting models the direction of this relationship is the same as that found using Tree 1.
The result of our second replication, using Tree 3 ( Figure S3), were also highly consistent with the results obtained using Tree 1 and Tree 2. In our bivariate phylogenetic analyses, we find that that the same variables are significantly associated with professional religious specialists, and in the same direction. For the exploratory phylogenetic path analyses using Tree 3, there are 15 best fitting models (within 2 CICc of the lowest CICc). As is the case with Tree 1 and Tree 2, across all of these best fitting causal models, the only variable that is significantly and directly causally linked to professional religious specialists is food storage. In all 15 of the best fitting models the direction of this relationship is also the same as that found using Tree 1 and Tree 2.
We also replicated our exploratory phylogenetic path analyses using Tree 1 without any restrictions on the possible causal paths. This means that cultural variables, such as professional religious specialists, could be predictors of environmental variables, such as environmental predictability. This unrestricted analysis was highly computationally intensive and involved testing the fit of a total of 3,771,782 possible model structures. This includes a large number of theoretically implausible models, such as models in which professional religious specialists affect environmental predictability. The results of these models found 24 models within 2 CICc of the best fitting model. Of these best fitting causal models, the only variable that is significantly and directly causally linked to professional religious specialists is food storage. In 22 of the 24 best fitting models the direction of the causal path was the same as that reported in the main text (Food Storage → Professional Religious Specialists).
However, in two of the best fitting models the direction of the causal path was reversed (Professional Religious Specialists → Food Storage). These two models were the 16 th (m1852899) and 17 th (m1852829) best fitting models and were 1.67 and 1.77 CICc away from the best fitting model. This means that while these models meet the threshold for inclusion in the best fitting models, they represent only a small proportion of the total sample of best fitting models and are assigned relatively low weighting.
The results of our follow up analyses indicate that the central findings of our study hold when using alternative approaches to branch length scaling, and without placing restrictions on the path structures of models. Full details of these follow up analyses, including reproducible code, are available through the OSF project page (https://osf.io/6jdqr/).      Tables   Table S1. Index of 15 hypotheses about the cultural and ecological factors associated with professional religious specialists in small-scale human societies. Relationships indicate the direction of the predicted causal relationship between religious specialists and each other variable. "Soviet scholars regard shamanism as a form of religious practice in pre-class societies. Some posit an early collective stage in which shamanic practice was open to anyone, followed by the development of professional specialization. As specialists who demanded payments and sacrifices from clients, shamans served to promote and legitimate growing economic inequality and the rise of feudal society." (Atkinson, 1992) Social Classes "What we suggest, then, is that this emphasis on spiritual authority and power is one of several strategies for acquiring power that has a number of advantages for aggrandizers in transegalitarian and chiefdom communities …. by using spiritual power as a subterfuge to justify and control excess resources and labor, to blunt possible accusations of self-aggrandizement, and as a commodity that cannot be stolen, aggrandizers could circumvent a number of problems in trying to promote their self-interests." (Hayden & Villeneuve, 2010) Slavery 2 Social Inequality → RS (positive effect) "the central forces transforming shamanism are a sedentaryagriculture lifestyle, the process of political integration of local communities into politically stratified societies, and the presence of classes (with castes or hereditary slavery)" (M. J. Winkelman, 1990) Social Classes "The initial hypothesis that systematic differences in magicoreligious practitioners occur as a function of social complexity is refined to suggest that agriculture, political integration, and the presence of social classes are the specific conditions that are responsible for the creation of the different types and configurations of magico-religious practitioners….. A second major aspect of the magico-religious phenomenon is the control of social power in sedentary agricultural societies, particularly those with hierarchical political integration." (M. Winkelman, 1986) "The notion that religion has its origin in society or social expression is consistent with the findings here, with the Priest role arising as a consequence of social and economic conditions (agriculture, political integration). However, the activities of the Shaman and related trance practitioners cannot be seen as arising out of religion in any fundamental sense … although these experiences are apparently shaped by social conditions (e.g., soul flight versus possession; see Bourguignon and Evascu 1977;Winkelman 1984, chapter 6;1986), the basic altered state that gives rise to the role and capacity is not determined by culture or social experience, but by human psychophysiology." (M. Winkelman, 1986) Slavery "Becoming a shaman also provided a way for low-status individuals to attain prestige, such as in some hierarchical societies of the Pacific Northwest (Gunn 1966), while in other instances, shamans were regarded as attractive sexual partners." (Singh, 2017) 3 RS → Cooperation (positive effect) "Shamans are healers, ritual leaders, and influential members of society whose keen insight and success in solving social problems (Rossano 2007;Winkelman 1990) can lead to wealth, power, and access to mates. Shamanism acts as a mechanism to reinforce social norms, encouraging group cooperation through ritual and social bonding, and calming anxiety during times of resource stress (Hayden 1987;Rossano 2007;Winkelman 1990;Winkelman 2010)." (Peoples et al., 2016) Internal Conflict "By binding supernatural authority to the punishments incurred when taboos are broken, the shaman reinforces group norms. Whatever the specific task, the good shaman serves as an emissary to the supernatural world, intent on repairing or enhancing the community's relationship to the supernatural and to each other-or, as Vitebsky (2000, pp 63) puts it, 'The mystic is also a social worker." (Rossano, 2007) 4 Outgroup warfare → RS (positive effect) "According to this theory, and the writings of many other researchers (Malinowski 1948), magical practices result from individuals trying to control important, uncontrollable events. As Intergroup Conflict Warfare individuals are better able to deal with this uncertainty -or if they confront less of it -they will depend less on shamans. Although a world completely devoid of uncertain outcomes seems implausible, Hart and Pilling (1960) offered this hypothesis to explain the absence of shamanism among the Tiwi. Citing the adequate rainfall, the consistent food supply, a dearth of dangerous animals, the rarity of tropical diseases, climatic docility, and the absence of antagonistic neighbors, the ethnographers concluded, "[The Tiwi] never invented magic to control their environment because their physical environment was on the whole a satisfactory and not a hostile universe" (p. 88)." (Singh, 2017) 5 RS → Outgroup warfare (positive effect) "Religious authority and military prowess go hand in hand. Religious authority allows a particular tribal leader to solve the large-scale collective action problem of uniting a group of autonomous tribes. To a much larger degree than economic benefit, religious authority can explain why a free tribal people would be willing to make a permanent delegation of authority to a single individual and that individual's kin group. The leader can then use that authority to create a centralized military machine that can conquer recalcitrant tribes as well as ensure domestic peace and security, which then reinforces the leader's religious authority in a positive-feedback loop." (Fukuyama, 2011) Intergroup Conflict Warfare 6 Food storage → RS (positive effect) "The Neolithic Revolution is believed to have paved the way for civilization, class society, and the state. The basic assumption is Food Storage always that only agriculture was able to generate a regular economic surplus sufficient to maintain a nonproductive class, such as priests, warriors, bureaucrats, and the like….This view, however, can no longer be maintained, for quantitative studies (e.g. Lee 1968Lee , 1969 show that, contrary to what has generally been assumed, huntergatherers do not work hard to make a living…. The explanation in terms of surplus definitely has to be replaced by a new one… preliminary evidence suggests that the relevant factor for the development of inequalities is not the presence or absence of agriculture, but the presence or absence of a storing economy, whether it be hunting-gathering or agricultural." (Testart et al., 1982) 7 Sociopolitical leadership → RS (positive effect) "the central forces transforming shamanism are a sedentaryagriculture lifestyle, the process of political integration of local communities into politically stratified societies, and the presence of classes (with castes or hereditary slavery)" (M. J. Winkelman, 1990) Social Leadership "The initial hypothesis that systematic differences in magicoreligious practitioners occur as a function of social complexity is refined to suggest that agriculture, political integration, and the presence of social classes are the specific conditions that are responsible for the creation of the different types and configurations of magico-religious practitioners….. A second major aspect of the magico-religious phenomenon is the control of social power in sedentary agricultural societies, particularly those with hierarchical political integration." (M. Winkelman, 1986) "The notion that religion has its origin in society or social expression is consistent with the findings here, with the Priest role arising as a consequence of social and economic conditions (agriculture, political integration). However, the activities of the Shaman and related trance practitioners cannot be seen as arising out of religion in any fundamental sense … although these experiences are apparently shaped by social conditions (e.g., soul flight versus possession; see Bourguignon and Evascu 1977;Winkelman 1984, chapter 6;1986), the basic altered state that gives rise to the role and capacity is not determined by culture or social experience, but by human psychophysiology." (M. Winkelman, 1986) 8 RS → Sociopolitical leadership (positive effect) "certain quality of an individual personality, by virtue of which he is set apart from ordinary men and treated as endowed with supernatural, superhuman, or at least specifically exceptional powers or qualities. These are such as are not accessible to the ordinary person, but are regarded as of divine origin or as exemplary, and on the basis of them the individual concerned is treated as a leader" (Weber, 1958) Social Leadership "There is evidence that in shamanism lie the origins of formal politics.... though there have been societies with shamans with no acknowledged political leaders, there have been few if any societies with a political leaders but no religious experts. And in some societies the shaman and the political leader have been one and the same." (Wright, 2010) 9 Sedentism → RS (positive effect) "the central forces transforming shamanism are a sedentaryagriculture lifestyle, the process of political integration of local communities into politically stratified societies, and the presence of classes (with castes or hereditary slavery)" (M. J. Winkelman, 1990) Number of Annual Moves "The initial hypothesis that systematic differences in magicoreligious practitioners occur as a function of social complexity is refined to suggest that agriculture, political integration, and the presence of social classes are the specific conditions that are responsible for the creation of the different types and configurations of magico-religious practitioners….. A second major aspect of the magico-religious phenomenon is the control of social power in sedentary agricultural societies, particularly those with hierarchical political integration." (M. Winkelman, 1986) "The notion that religion has its origin in society or social expression is consistent with the findings here, with the Priest role arising as a consequence of social and economic conditions (agriculture, political integration). However, the activities of the Shaman and related trance practitioners cannot be seen as arising out of religion in any fundamental sense … although these experiences are apparently shaped by social conditions (e.g., soul flight versus possession; see Bourguignon and Evascu 1977; Distance of Annual Moves Winkelman 1984, chapter 6;1986), the basic altered state that gives rise to the role and capacity is not determined by culture or social experience, but by human psychophysiology." (M. Winkelman, 1986) 10 Population size → RS (positive effect) "Shamanism is a complex set of practices, often involving ceremonies, initiations, mythologies, and music. Thus, like other forms of complex culture, it can disappear following demographic fragmentation, isolation, and bottlenecks (Henrich 2004;Shennan 2001)." (Singh, 2017) Consumer Group Size Max Camp Size Total Population Size 11 Resource stress → RS (positive effect) "The shaman may also be asked to relieve community suffering from such things as a natural disaster, poor harvest, social discord, hunger, or illness. Individuals may seek out the shaman when suffering from sickness or misfortune. In either situation, the shaman's function is to communicate with the ancestors or spirits who may be responsible for the perils being suffered." (Rossano, 2007) Mean Environmental Productivity "According to this theory, and the writings of many other researchers (Malinowski 1948), magical practices result from individuals trying to control important, uncontrollable events. As individuals are better able to deal with this uncertainty -or if they confront less of it -they will depend less on shamans. Although a world completely devoid of uncertain outcomes seems implausible, Hart and Pilling (1960) offered this hypothesis to explain the absence of shamanism among the Tiwi. Citing the adequate rainfall, the consistent food supply, a dearth of dangerous animals, the rarity of tropical diseases, climatic docility, and the absence of antagonistic neighbors, the ethnographers concluded, "[The Tiwi] never invented magic to control their environment because their physical environment was on the whole a satisfactory and not a hostile universe" (p. 88)." (Singh, 2017) 12 Natural disasters → RS (positive effect) "The shaman may also be asked to relieve community suffering from such things as a natural disaster, poor harvest, social discord, hunger, or illness. Individuals may seek out the shaman when suffering from sickness or misfortune. In either situation, the shaman's function is to communicate with the ancestors or spirits who may be responsible for the perils being suffered." (Rossano, 2007) Environmental Predictability "According to this theory, and the writings of many other researchers (Malinowski 1948), magical practices result from individuals trying to control important, uncontrollable events. As individuals are better able to deal with this uncertainty -or if they confront less of it -they will depend less on shamans. Although a world completely devoid of uncertain outcomes seems implausible, Hart and Pilling (1960) offered this hypothesis to explain the absence of shamanism among the Tiwi. Citing the adequate rainfall, the consistent food supply, a dearth of dangerous animals, the rarity of tropical diseases, climatic docility, and the absence of antagonistic neighbors, the ethnographers concluded, "[The Tiwi] never invented magic to control their environment because their physical environment was on the whole a satisfactory and not a hostile universe" (p. 88)." (Singh, 2017) 13 Resource unpredictability → RS (positive effect) "The shaman may also be asked to relieve community suffering from such things as a natural disaster, poor harvest, social discord, hunger, or illness. Individuals may seek out the shaman when suffering from sickness or misfortune. In either situation, the shaman's function is to communicate with the ancestors or spirits who may be responsible for the perils being suffered." (Rossano, 2007) Environmental Predictability "According to this theory, and the writings of many other researchers (Malinowski 1948), magical practices result from individuals trying to control important, uncontrollable events. As individuals are better able to deal with this uncertainty -or if they confront less of it -they will depend less on shamans. Although a world completely devoid of uncertain outcomes seems implausible, Hart and Pilling (1960) offered this hypothesis to explain the absence of shamanism among the Tiwi. Citing the adequate rainfall, the consistent food supply, a dearth of dangerous animals, the rarity of tropical diseases, climatic docility, and the absence of antagonistic neighbors, the ethnographers concluded, "[The Tiwi] never invented magic to control their environment because their physical environment was on the whole a satisfactory and not a hostile universe" (p. 88)." (Singh, 2017) (positive effect) "The shaman may also be asked to relieve community suffering from such things as a natural disaster, poor harvest, social discord, hunger, or illness. Individuals may seek out the shaman when suffering from sickness or misfortune. In either situation, the shaman's function is to communicate with the ancestors or spirits who may be responsible for the perils being suffered." (Rossano, 2007) Pathogen Load "According to this theory, and the writings of many other researchers (Malinowski 1948), magical practices result from individuals trying to control important, uncontrollable events. As individuals are better able to deal with this uncertainty -or if they confront less of it -they will depend less on shamans. Although a world completely devoid of uncertain outcomes seems implausible, Hart and Pilling (1960) offered this hypothesis to explain the absence of shamanism among the Tiwi. Citing the adequate rainfall, the consistent food supply, a dearth of dangerous animals, the rarity of tropical diseases, climatic docility, and the absence of antagonistic neighbors, the ethnographers concluded, "[The Tiwi] never invented magic to control their environment because their physical environment was on the whole a satisfactory and not a hostile universe" (p. 88)." (Singh, 2017) 15 Seasonal resource variation → RS (positive effect) "According to this theory, and the writings of many other researchers (Malinowski 1948), magical practices result from individuals trying to control important, uncontrollable events. As

Environmental Variation
Environmental Predictability individuals are better able to deal with this uncertainty -or if they confront less of it -they will depend less on shamans. Although a world completely devoid of uncertain outcomes seems implausible, Hart and Pilling (1960) offered this hypothesis to explain the absence of shamanism among the Tiwi. Citing the adequate rainfall, the consistent food supply, a dearth of dangerous animals, the rarity of tropical diseases, climatic docility, and the absence of antagonistic neighbors, the ethnographers concluded, "[The Tiwi] never invented magic to control their environment because their physical environment was on the whole a satisfactory and not a hostile universe" (p. 88)." (Singh, 2017) "The shaman may also be asked to relieve community suffering from such things as a natural disaster, poor harvest, social discord, hunger, or illness. Individuals may seek out the shaman when suffering from sickness or misfortune. In either situation, the shaman's function is to communicate with the ancestors or spirits who may be responsible for the perils being suffered." (Rossano, 2007)  * Indicates a p value of less than 0.05, ** indicates a p value of less than 0.01, and *** indicates a p value of less than 0.001