2.1. Cultural evolution of commitment

Formal evolutionary game theoretical models show that increased threats to individual well-being contribute to the evolution of cooperative norm strength and policing; sudden reductions in within-group, individual payoffs induce stronger norm adherence and favour more costly punishment of norm violators (Roos et al., Reference Roos, Gelfand, Nau and Lun2015). Empirical evidence is consistent with these models; cross-culturally, such ‘cultural tightness’ corresponds to various threats to well-being (Gelfand et al., Reference Gelfand, Raver, Nishii, Leslie, Lun, Lim and Yamaguchi2011; Jackson et al., Reference Jackson, Gelfand and Ember2020). In general, then, it follows that under threatening conditions, individuals will flock towards and adhere more strongly to cultural systems that are especially adept at alleviating such stress.

As a social system that can (a) soothe anxiety (Lang et al., Reference Lang, Krátkỳ, Shaver, Jerotijević and Xygalatas2015; Pargament et al., Reference Pargament, Smith, Koenig and Perez1998; Sosis, Reference Sosis2007; Sosis & Handwerker, Reference Sosis and Handwerker2011), (b) bolster within-group cooperation and solidarity (Lang et al., Reference Lang, Purzycki, Apicella, Atkinson, Bolyanatz, Cohen and Henrich2019; Sosis & Bressler, Reference Sosis and Bressler2003; Sosis & Ruffle, Reference Sosis and Ruffle2003), (c) flexibly attend to novel perturbations in and threats to cooperation (Purzycki et al., Reference Purzycki, Stagnaro and Sasaki2020, Reference Purzycki, Bendixen, Lightner and Sosis2022a; cf. McNamara et al., Reference McNamara, Norenzayan and Henrich2016) and (d) include beliefs for the punishment of violations of moral and other social norms (Johnson, Reference Johnson2005, Reference Johnson2016; Norenzayan, Reference Norenzayan2013; Purzycki et al., Reference Purzycki, Willard, Klocová, Apicella, Atkinson, Bolyanatz and Ross2022c; Stark, Reference Stark2001), religion should be especially attractive for individuals enduring insecurity and stress. In other words, under duress, individuals are likely to intensify their commitment to the kinds of traditional norms and behaviours that have served prior generations and facilitated the kinds of social support that bolsters one's survival and reproduction, and religion has been shown to be especially supportive of these traits (Henrich et al., Reference Henrich, Bauer, Cassar, Chytilová and Purzycki2019; Inglehart, Reference Inglehart2020; Purzycki & Sosis, Reference Purzycki and Sosis2022; Sibley & Bulbulia, Reference Sibley and Bulbulia2012). As stressed by Inglehart (Reference Inglehart2020), the kinds of existential security that religion might uphold are both physical and intellectual (i.e. religious worldviews provide answers to life's bigger questions). We focus here on the former, hence, the ‘material insecurity hypothesis’.

2.2. Material insecurity hypothesis

While not without complications, a wide-ranging sociological literature is consistent with the material insecurity hypothesis (cf. Stark, Reference Stark1999). Yet a host of socio-demographic factors associated with material insecurity also appear to increase religiosity at the individual and group levels (see Storm, Reference Storm2017). For example, food insecurity, poverty and wealth inequality have been found to co-vary with religiosity (Baimel et al., Reference Baimel, Apicella, Atkinson, Bolyanatz, Cohen, Handley and Purzycki2022; Hekmatpour, Reference Hekmatpour2020; Höllinger & Muckenhuber, Reference Höllinger and Muckenhuber2019; Inglehart, Reference Inglehart2020; Norris & Inglehart, Reference Norris and Inglehart2011; Purzycki et al., Reference Purzycki, Ross, Apicella, Atkinson, Cohen, McNamara and Henrich2018b; Ruiter & Van Tubergen, Reference Ruiter and Van Tubergen2009; Solt et al., Reference Solt, Habel and Grant2011), but as we are presently interested in assessing whether or not food security causes changes in religiosity, we must declare a particular causal direction. Furthermore, in addition to the predicted rise in religiosity owing to insecurity associated with age and proximity to death (Jong & Halberstadt, Reference Jong and Halberstadt2018), some cross-cultural evidence suggests that important indices of religious commitment can change across the lifespan in important ways (Bengtson et al., Reference Bengtson, Silverstein, Putney and Harris2015; Purzycki & Bendixen, Reference Purzycki and Bendixen2020; Shaver & Sosis, Reference Shaver and Sosis2014).

Additionally, sex is variously associated with religiosity (Schnabel et al., Reference Schnabel, Hackett and McClendon2018; Vardy et al., Reference Vardy, Moya, Placek, Apicella, Bolyanatz, Cohen and Atkinson2022; Walter & Davie, Reference Walter and Davie1998) as are fertility and family size (Blume, Reference Blume, Voland and Schiefenhövel2009; Glavatskaya et al., Reference Glavatskaya, Borovik and Thorvaldsen2018; Inglehart, Reference Inglehart2020; Shaver et al., Reference Shaver, Power, Purzycki, Watts, Sear, Shenk and Bulbulia2020). On evolutionary theoretical grounds, family size in particular has a complicated relationship with material insecurity (Lawson et al., Reference Lawson, Alvergne and Gibson2012; Purzycki et al., Reference Purzycki, Ross, Apicella, Atkinson, Cohen, McNamara and Henrich2018b; Strassmann & Gillespie, Reference Strassmann and Gillespie2002), but evidence points to a relationship nonetheless. Similarly, education level is directly implicated in the number of children one has (Becker et al., Reference Becker, Cinnirella and Woessmann2010; Goodman et al., Reference Goodman, Koupil and Lawson2012) and shows a widely touted – but no less complicated – relationship with security and well-being (see Boarini & Strauss, Reference Boarini and Strauss2010; Desjardins, Reference Desjardins2008; Giambona et al., Reference Giambona, Porcu and Sulis2022).

In addition to its role in material security, educational attainment is also of particular interest in the study of religion. Exposure to formal, secular education has long been asserted to weaken religious faith owing to the cultivation of critical thinking skills (Dawkins, Reference Dawkins2016; Weber, Reference Weber1993 [1920]). In other words, educational attainment might both directly and indirectly decrease religiosity. However, while some cross-national assessments do indeed show negative associations between years of formal education and religiosity (e.g. Inglehart, Reference Inglehart2020; Norris & Inglehart, Reference Norris and Inglehart2011, pp. 271–274), others reveal that both the magnitude and even the direction of the estimated association exhibits substantial cross-cultural variability (e.g. Albrecht & Heaton, Reference Albrecht and Heaton1984; Pew Research Center, 2016; Ruiter & Van Tubergen, Reference Ruiter and Van Tubergen2009; Schwadel, Reference Schwadel2011, Reference Schwadel2015).

As insightful as this literature is, it has some important limitations worth addressing. Theoretically, despite the explicit causal language across the theories that motivate this work, there is a distinct lack of causal modelling in the literature. This makes it difficult to synthesise various works into a single causal framework. Indeed, the empirical efforts are primarily correlational and do not formally model the manifold causal pathways by which education and food security could affect individual religious commitments. So, while theory is causal, the models are often informal and analyses typically focus on significant correlations rather than counterfactual states.

In terms of data, the bulk of this work often focuses on group-level relationships. Yet, group-level correlations can mask or muddle individual-level effects (Robinson, Reference Robinson2009; Simpson, Reference Simpson1951; Storm, Reference Storm2017; see Supporting Information Section 2 of the present report for further discussion and demonstration of this point). As such, individual-level, high-resolution data are important for adequately assessing individual-level hypotheses. In terms of sampling, most of the data that populates these studies are collected among individuals in fully market-integrated industrialised societies. Given the variation in food security and access to education the world exhibits, we should be able to see the anticipated effects of these indices of material security on religiosity across the range of communities. Here, we offer a modest step towards rectifying these issues by (a) synthesising a host of hypotheses into a single causal model, (b) being explicit about and simulating the assumptions of our model, (c) applying this model to data sampled from a diverse range of societies and (d) attending to individual-level effects using a general and principled causal inference method.

3.1. Causal inference

Across academic fields, there is increasing recognition that (a) causal explanatory accounts of natural phenomena require explicitly derived models and (b) estimating causation is best approached from an interventionist perspective (variously known as the manipulationist, counterfactual or potential outcomes frameworks; Hernan & Robins, Reference Hernan and Robins2020; Morgan & Winship, Reference Morgan and Winship2015; Pearl et al., Reference Pearl, Glymour and Jewell2016; Woodward, Reference Woodward2005). In this view, establishing causation requires that we not only examine what our data tell us about worlds that exist, but also form deliberately-calculated inferences of causal relationships derived from worlds that would exist, given particular conditions.

Quite often, the standard procedure of including a statistical model with all variables under consideration simply assumes that all variables are predictive of the dependent variable and ignores any internal causal structure. Yet, ‘controlling’ for everything can actually create more problems than it solves because otherwise ignored causal relationships that are internal to those predictors can introduce further confounding effects (e.g. conditioning on a common descendant of the exposure and outcome – collider bias) or suppress actual causal effects (e.g. conditioning on a mediator or a post-treatment variable). Put simply, under some conditions, holding some variables constant (i.e. controlling) can be inconsequential, while under others it can be disastrous to inference (Achen, Reference Achen2005; Cinelli et al., Reference Cinelli, Forney and Pearl2020; Westreich & Greenland, Reference Westreich and Greenland2013; Wysocki et al., Reference Wysocki, Lawson and Rhemtulla2022). The conditions that matter here are the causal structure of variables germane to testing hypotheses. Making these conditions explicit is the task of the researcher.

To obtain our results, we use a general causal inference method for eliciting and contrasting outcomes under hypothetical exposures in observational settings (e.g. Ahern et al., Reference Ahern, Hubbard and Galea2009; Naimi et al., Reference Naimi, Cole and Kennedy2017; Robins, Reference Robins1986; Snowden et al., Reference Snowden, Rose and Mortimer2011). This method – variously referred to as g-formula, g-computation and standardisation (see Bendixen, Reference Bendixen2023; Hernan & Robins Reference Hernan and Robins2020, chap. 13; VanderWeele et al. Reference VanderWeele, Jackson and Li2016) – involves three main steps: first, given a well-defined estimand and explicit causal assumptions, an appropriate statistical model is fitted that adjusts for relevant confounding variables. In the context of the present study, this first step is undertaken in Sections 3.2 and 4.3 below. Second, using the fitted model in the first step, we predict (or impute) outcome values under varying hypothetical levels of a focal predictor while holding all other variables as observed. This step follows directly from the interventionist approach to causal inference outlined above, in that it simulates changes in potential outcomes under manipulated, counter-factual exposures. Third, the imputed or predicted values from the previous step are contrasted and summarised by some appropriate inferential statistic, such as the mean and an interval. This step effectively amounts to averaging over the joint distribution of the observed covariates, yielding an average or marginal effect of the exposure in the sample. For the present analysis, for each individual and outcome response option, we compute the posterior mean under hypothetical exposure levels. For this effect to have a valid causal interpretation, a set of identifiability assumptions needs to be met (Hernan & Robins, Reference Hernan and Robins2020). We discuss these in Section 7 of the Supporting Information.

3.2. Causal model

Our guiding causal model (Figure 1) incorporates the hypotheses examined in the aforementioned literature (Section 2.2). The variables and causal relationships we consider are drawn from predicted directions of previous studies. The model is also informed and constrained by our knowledge of what variables we could use in the dataset.

Figure 1. Directed acyclic graph of the assumed causal structure. Bold variables are focal variables.

The focal causal hypotheses are that religiosity increases as a function of exposure to education and food insecurity. To the extent that years of formal education constitute exposure to a particular form of instruction that has some effect on the degree to which one is religious, and concern about access to food approximates material insecurity, we use these direct effects as our focal theoretical predictions and empirical estimands: Education (Amount) → Religiosity, and Food Insecurity → Religiosity. We also incorporate some of the other aforementioned individual-level predictors into our model, including sex, age, and number of children. The model wagers that they comprise a causal network of influence between amount of education and religiosity. Specifically, it posits an internal causal structure to these variables: Sex → Food Security ← Age. If Children ← Food Security, conditioning on children would not be necessary if Food Security were held constant (for an introduction to causal graph analysis, see Pearl et al., Reference Pearl, Glymour and Jewell2016). However, if Children → Food Security, this creates a collider bias with Food Security and Education (Amount). We therefore take this more conservative route and hold the number of children constant. We also posit that Education (Amount) → Children. While this creates an indirect path through Children, as we are conditioning on number of children already and there are no other posited confounds, we are avoiding this problem. Furthermore, the amount of education should contribute to one's food security which should, in turn, contribute to religiosity. Sex and age confound the focal causal path, and there is an indirect effect of education amount on religiosity that runs through food security.

It is important to note that the type of education might matter here. We recognise that individuals could have had religious education in some contexts whether by default or by choice, and that choice could be facilitated by food security, and so forth. In this model, we posit only that it has an effect on religiosity. This is a very important assumption for a couple of reasons. First, this is unmeasured in the dataset. According to this model structure, however, its absence is inconsequential to estimating the target effects; we have no need for conditioning on type of education. If we posit that it was caused by food security, not having measured it would be inconsequential to the target estimate as long as we hold food security constant (i.e. the indirect effect through education type would be blocked). However, if the type of education were posited to cause food security, this would both confound the effect through food security, but also open a confounding path between education amount through food security. With this important note in mind, we continue under the assumption that type of education only affects religiosity. In Section 3 in the Supporting Information, we walk through a simulation of this model that shows that we can recover target estimates with and without this particular variable.

4.1. Participants

We used data from the publicly available Evolution of Religion and Morality Project dataset (for motivation, design details, and further details about each field site, see Lang et al., Reference Lang, Purzycki, Apicella, Atkinson, Bolyanatz, Cohen and Henrich2019; Purzycki et al., Reference Purzycki, Apicella, Atkinson, Cohen, McNamara, Willard and Henrich2016, Reference Purzycki, Henrich, Apicella, Atkinson, Baimel, Cohen and Norenzayan2018a, Reference Purzycki, Lang, Henrich and Norenzayan2022b). The entire dataset consists of 2228 participants from 15 diverse populations around the world (see Fig. 2), but owing to there being no data on the outcome, we dropped one site (Hadza). Participants were recruited in a variety of ways; some of them were recruited on the basis of self-identification with a locally relevant religious tradition. There are no indications that recruitment on this basis contributed to the selection of individuals who are more religiously committed. All participants were 17 years of age or older. All participants engaged in a variety of economic game experiments designed to examine the breadth of religiously motivated cooperation. Participants were informed that their participation would be anonymous, any identifying information would be kept under lock and key and they had the right to withdraw at any point. Table 1 details the key descriptive statistics of each field site.

Figure 2. Map of 14 field sites. Note that there are two field sites in Tanna (Inland and Coastal).

Table 1. Means (standard deviations) for some individual-level demographic variables and group-level properties of communities sampled. ‘Mjr trad.’ refers to the major colonial and/or overarching tradition of each site while ‘Mnr trad.’ refers to relatively minor – but locally salient – traditions in each site.

4.3. Statistical model

To model the level of commitment across the two variables, we used Bayesian ordered-logistic regressions (Bürkner & Vuorre, Reference Bürkner and Vuorre2019). In each population, j, we asked each individual, i, the commitment questions, outcomes of which, K _i, are ordered categorical responses. As such, we model them using a multilevel ordered categorical likelihood distribution that is parameterised using a linear model term, ξ _i, and a vector of random cut-points, κ:

(1)$$K_i\sim {\rm Ordered\;Categorical}( {\xi_i, \;\kappa } ) $$

The linear model ξ _i is then given by:

(2)$$\xi _i = \alpha _{{\rm SITE}[ i ] } + \beta _{{\rm SITE}[ i ] }^M M_i + \beta _{{\rm SITE}[ i ] }^A A_i + \beta _{{\rm SITE}[ i ] }^S S_i + \beta _{{\rm SITE}[ i ] }^E E_i + \beta _{{\rm SITE}[ i ] }^C C_i$$

where SITE[i] gives the cultural group of individual i, A _i is the age of individual i, S _i is a variable indicating if individual i is male, E _i is the years of formal education completed by individual i, C _i indicates an individual's number of children and M _i indicates an individual's food insecurity. To account for the fact that increments in the predictor (e.g. an additional year of education and year of age or an additional child) might have differential associations with the outcome, we model the effects of years of education, age and number of children monotonically, such that a covariate's β coefficient represents the expected average difference between two adjacent levels of the predictor (Bürkner & Charpentier, Reference Bürkner and Charpentier2020). Note that we also modelled age and years of education with Gaussian processes (see Section 4 in the Supporting Information) but found the monotonic model to yield a better fit. Food insecurity and sex are modelled as indicator variables. Since unobserved site-specific factors are likely to affect our predicted effects in different ways across sites, our statistical model fully varies effects across sites. It is plausible that, for example, one year of education in one place deeply affects one's worldview whereas in others, it might have no effect. By-site clustering partly handles site-specific confounds.

We define our priors, which we checked through prior predictive simulation (see Section 4 in the Supporting Information) as follows. The thresholds κ are modelled with Normal(0, 10), all β coefficients with Normal(0, 0.5) and all monotonic covariates with a Dirichlet distribution with α denoting a series of 2s equalling the number of levels of the predictor. The Dirichlet prior encodes the a priori assumption that any of the levels of the predictor could be more or less likely than the others (McElreath, Reference McElreath2020, pp. 392–394). Variance components for the varying effects were set at Exponential(1) and the correlation matrix of the variance components to LKJCorr(4) (Lewandowski et al., Reference Lewandowski, Kurowicka and Joe2009). Taken together, the specified set of priors is weakly regularising.

We performed all analyses using the brms package (Bürkner, Reference Bürkner2017, Reference Bürkner2018, Reference Bürkner2021) for R (R Core Team, 2021), an interface to the probabilistic programming language Stan (Carpenter et al., Reference Carpenter, Gelman, Hoffman, Lee, Goodrich, Betancourt and Riddell2017). Chain and sampling diagnostics ($\hat{R} < 1.02$, well-mixing trace and rank plots, few divergent transitions, and relatively large effective sample sizes) were acceptable. Section 8 in the Supporting Information includes the full list of R packages used for this project as well as their dependencies and version number.