Forecasting Changes in Religiosity and Existential Security with an Agent-based Model

We employ existing data sets and agent-based modeling to forecast changes in religiosity and existential security among a collective of individuals over time. Existential security reflects the extent of economic, socioeconomic and humandevelopment provided by society. Ourmodel includes agents in social networks interactingwith one another based on the education level of the agents, the religious practices of the agents, and each agent’s existential security within their natural and social environments. The data used to inform the values and relationships among these variables is based on rigorous statistical analysis of the International Social Survey Programme Religion Module (ISSP) and the Human Development Report (HDR). We conduct an evaluation that demonstrates, for the countries and time periods studied, that our model provides a more accurate forecast of changes in existential security and religiosity than two alternative approaches. The improved accuracy is largely due to the inclusion of social networkswith educational homophilywhich alters theway inwhich religiosity and existential security change in the model. These dynamics grow societies where two individuals with the same initial religious practices (or belief In God, or supernatural beliefs) evolve di erently based on the educational backgrounds of the individuals with which they surround themselves. Finally, we discuss the limitations of our model and provide direction for future work.


Introduction
. Traditional approaches to the social-scientific study of religion have di iculty accounting for the varied ways in which education, existential security, and social networks shape the religiosity of an individual and a collective of individuals. Existential security reflects the extent of economic, socioeconomic and human development provided by society (Norris & Inglehart ). Furthermore, understanding the dynamics of the relationship among these factors over time adds complexity. Simulated environments have been used to test di erent theories of religious extremism and di erent proposals for reducing religious violence (Upal ; Bainbridge ; Iannaccone & Makowsky ). However, we are unaware of any previous work to model the relationships of the aforementioned factors with respect to a collective of individuals in agent-based model using established empirical data sources. In this paper we explore this issue. .
Our model can be parameterized for a given country and time period. Each agent has variables that reflect their education level, their existential insecurity, and di erent facets of their religiosity. In addition to these characteristics, each agent is connected to a social network of other agents. Based on their interactions with one another, the existential insecurity and religiosity variables of the agents change. We explore the behavior of the model to highlight those conditions that drive a collective of individuals over time to become more or less: ( ) religious and ( ) existentially secure.
. Next, we conduct an evaluation which demonstrates that, for the countries and time periods studied, our model provides a more accurate forecast of changes in existential security and religiosity than two alternative approaches. Both the data processing in our model that leverages existing data sets and the relative accuracy of the forecasts are novel contributions. We conclude the paper by discussing the context in which of our model is valid and our model's limitations. Despite the identified limitations, the model reflects an e ort to link empiricism and theory within an ABM as advocated by Flache et al. ( ).

Background
. In this section, we provide a brief overview of some of the major theoretical concerns related to the main variables in our model: religiosity, education, and existential security. In the sections that follow, we point to some of the previous simulation research related to these topics and present our own agent-based model.

.
Our first main variable is religiosity. Although the use of this term is highly contested, there is significant research from a wide variety of disciplines indicating that cohesion-enhancing ritual behaviors intended to engage supernatural agents that are believed to be interested in human a airs are virtually universal across cultures. These behaviors are the result of cognitive tendencies that are culturally evoked in a variety of ways (Atran ; Mc-Cauley & Lawson ; Wildman ; Shults ). In addition to this contextual variance, scholars of religion have also identified many ways in which these features vary across individuals in human populations. Despite this variance, the relevant features recur across cultures in recognizable patterns (Schuurmans-Stekhoven ; Caldwell-Harris ; Gebauer et al. ; Barber ). .
The second main variable is education. The relationship between education and religiosity is well researched. Most studies show that there is a significant negative correlation between education and religiosity (Blancke et al. ; Hill ; Hungerman ; Lewis ). It is well known that scientists, for example, are typically less religious than the general population; this is o en attributed to their educational background, which trains them to be analytical thinkers (Larson & Witham ; McCauley ; Norenzayan & Gervais ). There is some research that challenges the idea that education is always negatively related to religiosity. One study, for example, found that education was negatively associated with religiosity only for individuals with a weak religious background (Ganzach et al. ). However, the consensus in the literature is that analytical thinking skills taught in higher education are negatively correlated with religiosity even when controlling for the strength of one's religious background (Dutton & Van der Linden ; Lynn et al. ). For the purposes of our model, we use a specific measure of education from an international survey as a proxy for education variables (i.e. analytical thinking) that negatively correlate with religiosity. .
Existential security is our third main variable. The definition of existential security we adopt is related to the security axiom, which is based on the idea that "societies around the world di er greatly in their levels of economic and human development and socioeconomic equality -and consequently, in the extent to which they provide their people with a sense of existential security" (Norris & Inglehart ). This axiom is based on cross-cultural survey analyses that show that "the emergence of high levels of existential security tends to diminish anxiety and stress, promoting feelings of psychological well-being" -which, in turn, they argue, reduces the importance of religious values in people's lives (Norris & Inglehart , ). .
Other research has shown that as secular institutions increasingly provide existential security for a population, there is a reduced reliance on religious institutions, which appears to cause a decline in religious belief and participation (Norenzayan ; Stolz et al. ; Kay et al. ; Habel & Grant ). This is particularly clear in Scandinavian societies such as Norway, where people feel they are protected by a strong social safety net, and report higher levels of existential security (and lower levels of religiosity) than most other countries (Kvande et al. ). Quantifying existential security, as we do in our model, is similar to the quantification of other human perceptions such as happiness (Helliwell et al. ; Gore et al. ; Dodds et al. ; Mitchell et al. ) or well-being (Deaton ; McMahan & Estes ; Decancq & Lugo ). .
The agent-based model we describe in Section explores the relationship between religiosity, education, and existential security. In future work, we plan to explore other factors that influence religiosity, such as pluralism and freedom of expression. While related research e orts have included the use of statistics to support their analyses, none have used agent-based models to forecast the way in which these factors a ect changes in the religiosity of a population over time. Our work begins to fill this gap.

Related Work
. Our model draws on social influence modeling and is related to a number of existing agent-based models. We review social influence modeling first and then describe the ABMs which are most closely related to our work.

Social influence modeling .
Social influence is an important factor in human interaction. In encounters, individuals can modify their opinions, attitudes, beliefs, or behavior to mimic or to oppose those they interact with. These modifications can be the result of persuasion, uncertainty or peer pressure (Flache et al. ).
. Initial social influence models included actors whose opinion on a position was socially influenced on a continuous spectrum (Abelson ; Berger ; DeGroot ; Lehrer ). Examples of these models include agents determining the appropriate speed on an interstate. Later, modelers assumed that opinions did not vary on a continuous scale but instead reflected a choice between mutually exclusive options (Axelrod ; Latané ; Liggett ; Sznajd-Weron & Sznajd ).
Here, examples include political parties, music and movies. In some cases it is possible to model mutually exclusive options as nominal traits on a continuous spectrum (Nowak et al. ; Flache et al. ). .
Recently, researchers have identified three classes of social-influence models that are the most popular in agent-based modeling. These classes of social influence models are: ( ) assimilative social influence ( . In assimilative social influence models, individuals are connected by a structural relationship and always influence each other to reduce opinion di erences. Here, if the network is connected the influence dynamics eventually create consensus (Flache et al. ). .
For models with similarity biased influence, only su iciently similar individuals can influence each other to reduce opinion di erences. How much similarity is su icient depends on other mechanisms included in the model (e.g. social identity, confidence in others, etc). With similarity based influence consensus can be avoided. However, if the similarity bias is su iciently strong, then multiple homogenous but distinct clusters of individuals emerge. Opinions, however, never leave the initial range (De uant et al. ; Hegselmann et al. ; Flache et al. ).
. In models with repulsive influence, when individuals are too dissimilar they can influence each other to increase opinion di erences. The amount of dissimilarity needed to trigger repulsive influence depends on other mechanisms included in the model (e.g. social identity, ego-involvement). Here, consensus can be avoided. Also, clusters can form and adopt maximally opposing views (bi-polarization). These dynamics allow opinions to leave the initial range (Flache et al. ). .
It is important to note that most social influence models represent opinions as a one-dimensional variable. However, recent research has shown the importance of using a multidimensional model to study opinion polarization (Li & Xiao ). We employ this multi-dimensional representation in our representation of religiosity.

Related simulations .
Other researchers have employed modeling and simulation methodologies to simulate culturally relevant change within individuals and collectives of individuals. Axelrod's dissemination of culture model simulated a variety of mechanisms showing how interaction between di erent cultural features challenges intuitive assumptions about individuals' beliefs and interpersonal behavior (Axelrod ). Chattoe and Hamill used an agent-based model to demonstrate that social networks situate beliefs. Specifically, the model shows that what one knows depends largely on who one knows (Chattoe & Hamill ). Chattoe also demonstrated how church survival, in terms of membership, is a function of rational choice theory (Chattoe ).
. More recently, Upal and Gibbon developed an agent based system for simulating the dynamics of social identity beliefs that aimed at isolating factors that contribute to intergroup conflict (Upal & Gibbon ). Shults et al. developed two di erent computer simulations to demonstrate the impact of mortality salience on religiosity (Shults et al. ). At the same time, other researchers have used simulations of artificial societies composed of agents of two psychological types to create preferential self-organization into groups of ideological a inity ). Finally, agent-based modeling has also been used to study role of the emergence of a priestly class in solving large-scale social network coordination (Dávid-Barrett & Carney ), as well as the mechanisms that shape di erent modes of ritual interaction (Whitehouse et al. ).

The Simulation Model
. The goal of our model is for a given time period and country to predict changes in the existential security and the religiosity of a collective of individuals using an agent-based model. Here we describe the architecture of our model including: ( ) the entities within the model, ( ) the data sources used to initialize the entities, and ( ) the rules that dictate the interactions among them.

Model entities .
Our model is made up of agents interacting through social networks in an environment defined by its existential security level. Recall, existential security reflects the extent of economic, socioeconomic and human development provided by a country (Norris & Inglehart ). Each agent has an education level, an existential insecurity level, and four variables that reflect their religiosity. Each agent is also connected to a subset of the other agents in the model through a social network. An overview of these entities, with their attributes and representation, is shown in Table . Religiosity variables .
The four variables that reflect religiosity in each agent are: ( ) religious formation, ( ) religious practice, ( ) supernatural beliefs, and ( ) belief in God. The identification of these variables is based on a rigorous statistical analysis of questions and responses from the International Social Survey Programme Religion Module (ISSP) (Davis & Jowell ). This data contains the cumulated variables of the ISSP "Religion" surveys of , and and includes , respondents from more than countries.
. Table shows the questions selected for statistical analysis as well as the questions' labels, number of valid Likert levels, and an indication of if we performed polarity inversion on the question. Questions V11-16 tap attitudes and values that are o en explored in relation to religiosity. Questions V20-24 were included to study the correlation between religiosity and confidence in institutions. Questions V25-27 are related to the relationship between religion and politics. Questions V28-33, V35, V37 and V51 are related to beliefs in the "supernatural," particularly in God (V28-29 and V35). Questions V46-48 are related to the attendance of the respondent's mother and father, and self-attendance during a formative period (ages -), which is o en considered to be particularly influential on religious belief and practice later in life (Henrich ; Gervais & Najle ; Lanman ). Questions V49 (frequency of praying), V50 (participation in church activities other than regular religious services) and ATTEND (frequency of participation in regular religious services) are related to the respondent's How o en do you take part in the activities or organizations of a church or place of worship other than attending services?

V51
Would you describe yourself as religious? Yes V64 Agree/Disagree: Overall, modern science does more harm than good.

V65
Agree/Disagree: We trust too much in science and not enough in religious faith.

Yes
ATTEND How o en do you attend religious services? Yes Table : Initial set of selected questions in the ISSP data set: variables' names, labels, number of valid levels and indication of polarity inversion.
current religious practice. Since science is o en in competition with supernatural beliefs, we also included questions V64-65. .
The goal of the analysis of this data was to filter and aggregate the questions in Table into a lower number of unobserved variables called factors. The analysis results in a factor score for each respondent, for each factor, based on the respondent's answers to the questions included in the factor. Factor scores are continuous numbers, which reflect the extent to which each respondent manifests each factor. For each factor, scores are distributed normally with a mean of and a standard deviation of . However, we normalize scores for each factor on a to scale. Thus, a value of reflects a respondent where the factor is heavily present and a value of reflects a respondent where the factor is not present at all (DiStefano et al. ). .
Ultimately, we identified questions that form the aforementioned four religiosity factors. These questions and the factors they form are shown in Figure . Additional details about the choices associated with each question, iterative steps within the factor analysis, and the resulting measures of fit and statistical significance are provided in (Lemos et al. ).
Figure : Religiosity factors and associated questions from analysis of ISSP data. Analysis details, fit statistics, and statistical significance are provided in (Lemos et al. ).

Existential security
. Each agent is connected to the existential security level of the environment. The existential security level of the environment reflects the percentage of the agents that feel the level of economic, socioeconomic and human development support provided to them is su icient. An agent determines if they feel existentially secure by checking if their value for existential insecurity is below the existential security value of the environment. The initialization of these entities and the variables that make up each agent are described next. Then the description of the interactions among agents and the environment are presented.

Data sources used to initialize entities .
The initialization of each entity is based on the country and time period being modeled. For a given country and a given year, the existential security level of the environment is initialized using data from the Human Development Report (HDR). The HDR is an annual multi-facetted analysis of wellbeing focused on key dimensions of human development including a long life, a healthy life, and a decent standard of living. The Human Development Index (HDI) is the summary measure used in the HDR for a country's achievement across these dimensions (Anand ). .
Similarly, based on a specific country and a specific year, an agent is initialized by randomly sampling an ISSP respondent. Using the respondent's data, the characteristics of an agent are parameterized. Specifically, the factor scores for the four religiosity factors of the chosen ISSP respondent initialize the agent's religiosity variables and the education level of the respondent parameterizes the agent's education level. .
The education level of the respondent is based on the answer to the question: "What is the highest level of education you've achieved?" Responses to the question are coded on a -point scale from: ( ) no formal education to ( ) university level with degree. We do not employ factor analysis for the education variable because this question measures a response that is very close to the variable we want to capture and there are no additional questions related to education within the ISSP. . Finally, the existential insecurity level of an agent is parameterized by sampling a uniform distribution between and . While this parameterization is not based on data, the interactions of the model tie each individual's existential insecurity to the existential security level of the environment. This relationship is described further in paragraph . . .
Each agent is also assigned a social network. The social network reflects an agent's set of active social relationships. The number of and weight of the links within the network is constructed by an algorithm (Conti et al. ) that mirrors human social networks observed in the wild (Hill & Dunbar ; Dunbar & Shultz ). This research indicates that human relationships have a hierarchical structure and an individual has ∼ active social relationships (Dunbar ). .
This limit is due to cognitive and time constraints of the individual (Hill & Dunbar ). The hierarchical structure of the network consists of a series of concentric layers of acquaintance with increasing sizes. The layers are hierarchically inclusive, so that each layer includes all inner levels. The innermost layer is the support group of the individual (size ∼ ) followed by the sympathy group (size ∼ ) and then entire active network (size ∼ ). This structure is depicted in Figure   .
Within the network each relationship is also characterized by a level of emotional closeness. Strong relationships have a higher level of emotional closeness than weaker ones. The closeness reflects both the frequency of contact and time since last contact within the period of one week (Roberts & Dunbar ; Conti et al. ). In our implementation the emotional closeness is also correlated with the similarity in the education level of agents. This relationship has been observed in a several di erent social anthropology studies (McPherson et al. , ; Mark ; Salzarulo ; Chattoe & Hamill ). The degree to which the educational level of the agent is correlated to the emotional closeness is determined by the parameter, Education Homophily (EH).

Interactions among entities .
Given these initialized entities we define rules that govern the interactions among our entities. First, we review rules that dictate how changing religiosity variables influence one another. Next, we describe how agent's interactions within their social network influence their existential insecurity and religiosity. Finally, we describe the rules related to agent interactions with the existential security level of the environment.
We use structural equation modeling (SEM) to organize the relationships among the four religiosity factors. SEM enables us to hypothesize an architecture of the relationships among the factors and assess the extent to which the hypothesized architecture matches the observed data from the ISSP. Using our four factors it is possible to construct unique SEM architectures. Of those , four models have identical fit statistics that are superior to the other . Of the four models with the best fit statistics, one places the factors in an order that is consistent with theories of religiosity posited in the scientific study of religion ( ; Norenzayan et al. ). This model is shown in Figure . Additional details and fit statistics for this and other candidate models are provided in (Lemos et al. ).
Figure : Structural equation model that organizes the four religiosity variables within an agent. .
The model describes quantitatively how changes in religious practices and religious formation factors of an agent predict changes in the agent's belief in God and the agent's supernatural beliefs. Religious formation reflects questions related to the religious upbringing of the agent. Since this factor reflects events that have happened in the agent's past, we do not update the value of this variable in our model. .
Recent research has shown the importance of using a multidimensional representation of opinion to study the dynamics of social influence on polarization (Li & Xiao ). We employ this same idea in our multi-dimensional representation of religiosity shown in Figure . This capability can create a situation where two agents come to a consensus in one dimension of religiosity (i.e. supernatural beliefs) while being influenced to have polar opposite opinions on another (i.e belief In God).
. Next, we describe how the value of an agent's religious practice and existential insecurity variables can change based on influence from the social network.

Updating Agent Religiosity Values Based On Social Network Interactions .
At each time step agents interact with the agents in their social network. Each time step corresponds to one week. In this interaction the religious practice variable and existential insecurity within the agent are influenced by the respective values of these variables in other agents within their social network. The extent to which the variable is influenced is determined by a time-dependent weighted average. Given a matrix IN that includes an entry for the influence of each of the N agents on every other agent, the total influence exerted on agent i is show in Equation . T T otalInf luence i and a set A that includes all agents enables us to define, A SNv,t,i . Set A contains the value of each variable v, at each time step t, for each agent j, throughout the simulation (A v,t,j ). A SNv,t,i is the influence exerted on agent i by his/her social network (SN ) for a given variable v at time t. Formally it is shown in Equation . A Based on the influence exerted on the agent's religious practice variable from their social network at time t (A SN RP,t,i ), new values for the agent's supernatural belief (A SN SB,t,i ), and belief in God (A SN BIG,t,i ) variables are generated using the equations in the model shown in Figure . .

Finally each value, based on the influence of the agent's interaction with the social network (A
, is combined with the agent's existing value for the respective variable. This combination is computed using Cobb-Douglas function (Cobb & Douglas ). We employ the Cobb-Douglas function because it is an established, flexible, and widely used method to aggregate the influence of the environment with the existing value of a variable though the parameter β. Formally, given a variable v, this combination is shown in Equation . A

Interacting With The Environment
.
Once every agent has computed the influence of their social network on the applicable variables, every agent interacts with the environment. An agent interacts with the environment by checking if their value for existential insecurity is below the existential security value of the environment. Recall for a given country and a given year, the environment is initialized with the country's HDI value. Also, recall that the existential insecurity value of agents is initialized with a random number sampled from a uniform distribution between and .
. This check reflects whether or not the agent is provided su icient existential security by the environment. Formally, the process of applying this check to all N agents at time step t to determine the influence of the agents on the existential security level of the environment (ES agentInf,t ) is shown in Equation .
Agents with high EI values are less likely to be provided su icient existential security, while agents with low EI values are more likely to be provided su icient existential security.
Finally, the value for the existential security of the environment in the next time step (t + 1) is computed by combining the value of the existential security value of the environment at time t with ES agentInf,t at time step t. This combination is also done with Cobb-Douglas function as described in Equation .
The choice to reuse this function and its parameter, β, improves parsimony by reducing the total number of parameters in our model. Recall, it also an established and widely used method to combine the intrinsic value of a variable with an environmental influence on the variable. An overview of the model and the data sources used to initialize the entities within the model is shown in Figure . Appendix A contains a parameter table for the model indicating the value of each parameter and the source of each parameter's value.
Figure : Overview our of model and the data sources used to initialize the entities. .
Each of the interactions among the entities use assimilative social influence. Recall, this social influence model leads to general consensus amongst agents. However, employing social networks with educational homophily and a multi-dimensional representation of religiosity allows our model to avoid consensus. These dynamics can grow societies where two individuals with the same initial religious practices (or belief In God, or supernatural beliefs) evolve di erently based on the educational backgrounds of the individuals with which they surround themselves. In the next section we evaluate our model against two alternatives to determine its e ectiveness in forecasting changes in religiosity and existential security.

Evaluation
. The e ectiveness of our agent-based forecasts for changes in religiosity and existential security is elucidated through empirical evaluation against alternative modeling approaches. Specifically, we compare our agentbased approach to: ( ) a baseline approach based entirely on historical data and ( ) a statistical approach that uses linear regression modeling (LR).
. For a given period of time, each model predicts changes in the religious practice, supernatural beliefs, belief In God, and existential security of the population of a given country. The baseline approach assumes there will be no changes in these factors from the most recent previous data. This approach mirrors predicting that the weather tomorrow will be the same as the weather today. The statistical approach uses regression to predict future changes in a variable using a weighted linear combination of the current variables (religious practice, education level, religious formation, existential security, supernatural beliefs, and belief in God).
. We evaluate each of these candidate models using a three part process. First, we identify similar time periods of measurement between the ISSP and HDR data. ISSP data was collected in , , and while HDR data was collected in , , , Given that there is no intersection of common data collection we use time periods of HDR data that are closest to the ISSP time periods. These time periods are ISSP: -/ HDR: -and ISSP: -/ HDR: -.

.
Next, for these time periods we identify the countries where data was collected in the ISSP and the HDR. There are countries where data is collected in both time periods for both data sources. These countries are: Germany, Hungary, Ireland, Italy, Netherlands, New Zealand, Norway, Philippines, Poland, United Kingdom, and United States. The full availability of all country data within the ISSP and HDR is described in Appendix B.
. Finally, we fit the parameters of each model. The baseline approach does not require this final step because it does not have any parameters. However, to fit the parameters of the LR models and our agent-based model, we perform an automated search over all the combinations of possible values for the parameters using data for all countries in the previous time period for which there is both ISSP and HDR data.
. The automated search identifies the parameters for the model that minimize the root mean squared error (RMSE) of the absolute error of the mean forecast (ē a ) for each variable for each country. Theē a for a country is determined taking the absolute value of the di erence between the mean value of a forecasted variable and the mean of the actual data for the variable in the ISSP and HDR. .
To avoid overfitting the LR models, we only include LR models with statistically significant variables. Furthermore, since the ISSP does not collect longitudinal data, each LR model is only trained on country-level data (i.e. mean RP, mean SB, mean BIG) as opposed to individual-level data. These two factors (statistical significance and only country-level data) result in LR models which include only one variable, the current country-level value of the variable being predicted. .
Recall, our data set only has two di erent time periods (ISSP: -/ HDR: -and ISSP: -/ HDR: -). As a result, we fit the parameters of the regression models and the agent-based model using data from the first time period. Then we evaluate the accuracy of each model's forecasts using data from the second time period. Since one time step in our agent-based model corresponds to one week, we simulation time steps ( years) in training and evaluation. The models and their parameters are shown in Table . . Next, we use the identified parameters to forecast the ISSP and HDR values for the upcoming time period. Accuracy of a forecast is measured by the RMSE of theē a for each variable for each country from the actual data ( ISSP: -/ HDR: -). In addition, we evaluate the RMSE of the absolute error of the standard deviation of the forecast (σe a ) for each country for the variables religious practice, supernatural beliefs, and belief in God from the actual data. We cannot evaluate the RMSE of σe a for the existential security variable because the HDR data is recorded on a per country level, not a per person level. This limitation means that the existential security data for each country is a single number and does not have a distribution or variance. The overall accuracy of each model is reported in Table . The results for each individual country for all approaches are reported in Appendix C and the results for our ABM are shown in Figure .   Table shows that our parameterized ABM outperforms the baseline and LR approaches. Recall, each approach forecasts changes in the four variables for a given country over a given time period. For each of the four variables that each model predicts, the ABM has the lowest RMSE for the: ( )ē a -absolute error of the mean forecast and ( ) σe a -absolute error of the standard deviation of the forecast. This means that the forecasts for each of the four variables from the ABM better match the central tendency and variation of changes observed in the ISSP and HDR than the competing alternatives. Based on the factor and the evaluation measure, the ABM ranges from as accurate (RP σea ) to × more accurate (ESē a ) than the next best alternative.

.
Next, we evaluate how e ectively each of the trained models forecast changes in the four factors for a set of countries they have not been trained on. In this evaluation we use the countries for which we have ISSP and HDR data from the most recent time period (ISSP -/ HDR -), but do not have data from the training time period (ISSP: -/ HDR: -). These countries are: Chile, Cyprus, Czech Republic, France, Portugal, Slovenia, Spain, Sweden, Switzerland, Denmark, and Japan.

.
This evaluation gives us insight into the robustness of the parameterizations identified in training. We use the term robust to reflect forecast accuracy for previously unseen countries and time periods. Again, we evaluate the RMSE ofē a and σe a of each of the three modeling approaches for the new countries. The overall accuracy of each model is reported in Table . The results for each individual country for all approaches are reported in Appendix D and the results for our ABM are shown in Figure .   The evaluation shows that the ABM continues to perform as well as the best approach even when it is applied to forecast countries where training data did not exist. Furthermore, for several evaluation measures it continues to be multiple times more accurate than the best alternative. These results provide some evidence that our model can be used to forecast current changes in the religiosity and existential security for other countries not included in the evaluation. However, it should be noted the performance of our model in Table (no training data) relative to the alternative approaches is not as strong as in Table (training data).

RMSE Religious Practice
. This could be a result of incorrectly fitting our model. Recall, the use of our model in the evaluation follows a logic: first parameters within our model are fit using a data set and then the accuracy of a forecast is evaluated using a di erent data set from the same data source. However, an ABM fit in this manner, like any other model, may o er predictions at a certain level of accuracy despite using inaccurate parameter values (Chattoe-Brown ; Gore et al. ). This can happen because relevant variables that should be included in the model (i.e. gender and ethnicity) are not included. .
We are aware of this limitation of our model and it may be manifested in the results of our evaluation. However, we are not aware of any other data sets, besides the ones we are currently employing, that could be used to evaluate our model and the two alternative models. Furthermore, our statistical analysis of the data sets has yet to reveal relevant significant di erences across other possibly relevant variables (i.e. gender and ethnicity) with respect to the religiosity factors identified in Figure . Next, we explore the features of the model that distinguish it from alternative approaches included in the evaluation. Then, we discuss this issue further along with other limitations of our model.

Model Exploration
. Three features distinguish our model from the alternative approaches included in the evaluation. These features are: ( ) the existence of social networks with educational homophily, ( ) the ability of agents to influence the religious practice and existential insecurity variables of one another via social networks, and ( ) the ability of agents to influence the existential security level of the environment. To highlight how these features create conditions that enable accurate forecasts we explore the dynamics of our ABM predictions.
. First, we consider how the existence of social networks with educational homophily a ect change in the agent's religious practice variable throughout the simulation. Existing research has shown that the level of education one has is inversely correlated with the extent of one's religious practice (Albrecht & Heaton ; Larson & Witham ; Glaeser & Sacerdote ; McCauley ; Norenzayan & Gervais ). This relationship also exists within the ISSP dataset where the two variables have a -. correlation with a p-value less than . .

.
Our ABM implements this correlation through the educational homophily (EH) parameter. Recall, the EH parameter controls the degree of educational uniformity within one's social network. This feature is important because it provides a means for the sustained existence of groups of individuals with high religious practices, belief in God, and supernatural beliefs in countries where the majority of the population has low religious practices, belief in God, and supernatural beliefs. In other words, this feature ensures that every agent in the population does not take on all the preferences of the majority. Figure : Change in the religious practices, belief in God, and supernatural beliefs of agents from all countries for agents with similar initial respective variable values, broken down by education level. .
Figure elucidates this feature in our ABM. It shows the amount of change for agents in the ABM with similar initial values for their religious practices, belief in God, and supernatural beliefs broken down by education level. For each variable (religious practice, belief in God, and supernatural beliefs) the population is formed by running the ABM for all countries and then matching each agent attached to a network where the majority of the members completed higher secondary education with a similar agent attached to a network where the majority of the members did not complete higher secondary education. Agents are considered similar if the initial values of the variable for the two agents di ers by less than . .
. . These dynamics can grow a society where two individuals with the same initial religious practices (or belief In God, or supernatural beliefs) evolve di erently based on the educational backgrounds of the individuals with which they surround themselves.
. Next, we explore: ( ) how the social networks in our ABM create changes in existential security for each of the countries and ( ) how those changes in the existential security level of the ABM environment correlate with changes in religious practice, supernatural beliefs and, belief In God. For each of the countries presented in our evaluation, the ABM predicts that the existential security of the environment will increase over time. The extent of that increase, how accurate the prediction is, and how the predictions correlates with predicted changes in the religiosity variables are shown in Figure . Figure : Predicted change from our ABM for the existential security, mean religious practices, mean belief in God, and mean supernatural beliefs of agents for each country from -.
. Figure shows that when the existential security level of the environment in our ABM is high (> . ) almost all agents feel existentially secure and the existential security level of the environment immediately increases. These types of predictions for our ABM are seen for countries including: Denmark, France, Germany, Ireland, Japan, Netherlands, New Zealand, Norway, Sweden, Switzerland, United Kingdom, and United States.
. However, if the initial existential security level of the environment is less than . , there is not any growth in the existential security of the environment during the first three years. Eventually, the social network interactions of the agents within the ABM result in fewer agents with extreme existential insecurity values, creating fewer agents that feel existentially insecure and ultimately creating an increase in the existential security of the environment starting a er year three. This increase continues through year ten. These types of predictions for our model are seen for countries including: Chile, Cyprus, Hungary, Italy, Poland, Portugal, Slovenia, and Spain.
. This pattern of predicting staggered existential security growth is even more pronounced for countries where the initial existential security value is less than . . This is the case with our ABM's prediction for the Philippines. Here there is no growth in the existential security of the environment for the first seven years of the prediction until finally there are enough social network interactions of the agents within the ABM to create a su icient number of agents that feel existentially secure resulting in an increase in the existential security of the environment from years seven to ten. .
It is important to note that our model does not predict existential security growth for countries where the initial existential security value of the environment is below . . Under these conditions the existential security of the environment decreases over time. Recall, the existential insecurity level of agents is uniformly distributed from [ , ]. As the agents interact, the existential insecurity level of each agent becomes less extreme, but since the existential security of the environment is below . most agents in the population still feel existentially insecure despite less extreme existential insecurity values. .
The behavior of our model under these conditions needs to be explored further and refined. The HDR data shows that the existential security value of the environment increases for many countries with an initial value less than . . As currently constructed our model cannot replicate this behavior. However, it is also important to note that all the countries in the HDR that show a decrease over time in the existential security value of the environment are countries where the initial value is less than . . In future work we will look to produce a more refined algorithm that takes into account each of these possible trajectories for countries with low existential security.
. Figure also shows that each increase in existential security predicted by our ABM is coupled with predictions for a decrease in belief in God. Furthermore, the magnitude of the predicted increase in existential security is almost exactly the same as the magnitude of the predicted decrease in belief in God. It is important to note that this relationship between existential security and belief in God is not encoded in any of the rules or interactions within our model. It emerges from the model interactions and the data that parameterizes the model for each of the countries. We will explore this result and gather more data as to why the same trend does not exist for supernatural beliefs and religious practice in future work.

Model Validity and Limitations
. Construct, internal, and external validity threats a ect our model. In addition to these validity threats, our model has a number of limitations. Here we review each of these areas and discuss how they relate to our model.

Construct validity .
Threats to construct validity concern the appropriateness of the measures used to represent the entities in our model. While the data processing we use to leverage existing data sets reflects a novel means to construct an ABM, it limits the construct validity of our model in four ways. .
First, recall that we employ factors scores to aggregate and convert the Likert scores in the ISSP survey into continuous variables in the construction of the religious architecture for our agents (Figure ). Constructing agents that reason with a continuous representation of discretized data (i.e. Likert scores) can bias the manner in which the survey data informs the model in some cases (Flache & Macy ). It is important to note though, that o en fundamental results for models that use a continuous representation of discrete data generalize to other models that with a similar formalization of social influence (Flache et al. ). .
Second, constructing a SEM that is based on aggregate data may reduce the means in which agents can store and change their individual opinions in the model (Epstein ; Chattoe-Brown ). However, the previous section shows that our model is capable of producing societies where two individuals with the same initial religious practices (or belief In God, or supernatural beliefs) evolve di erently based on the educational backgrounds of the individuals with which they surround themselves. .
Third, the social networks employed in our model are static. While the networks enable the variables associated with each agent to change over time, the structure of the network and the weight of each link does not change.
As a result, our networks do not reflect exogenous events (e.g., work change, death of a node, marriages, etc.). In future work we will explore the interplay of individual behavior and network dynamics by taking advantage of existing research in computational sociology. Specifically we will explore work on the co-evolution of behavior and social network structure (Fehl et al. ) as well as approaches to study the coevolution of networks in social dilemmas (Corten ; Bravo et al. ). We expect that work in this direction may enable us to highlight the self organization of religious groups related to education and existential security as shown in (Abrica-Jacinto et al. ).
. Finally, the largest concern related to construct validity involves equifinality. The construction of our model required us to make a series of reasonable assumptions to fill unanticipated specificity gaps. While we did our best to address each gap with the most reasonable assumption, many other reasonable assumptions could also have been made, and it is likely that some of the models resulting from those assumptions would produce models with similar, possibly even superior, levels of forecast accuracy (Poile & Safayeni ). The choices in our model serve as recommendations for other researchers tasked within similar problems but we are also open to a dialog on alternative selections within our model construction.

Internal and external validity
. Internal validity threats arise when factors a ect the dependent variables without the modelers knowledge. It is possible that some implementation flaws could have a ected the evaluation results. However, the algorithms we used to within our model passed several internal code reviews and the RMSE of theē a and σe a of the model for each of the countries reflect the mean value taken over replications.
. Threats to external validity occur when the results of the model cannot be generalized. Although the evaluation was performed for using ∼ years of data from two well known sources the results of our model cannot be generalized to: ( ) other countries, ( ) during di erent years or ( ) di erent data sets.

Limitations .
It is important to note that our model does not o er a more accurate prediction for religious practice, supernatural beliefs, belief In God and existential security than the two alternative models for every country. However, for the countries on which the models were trained our model o ers more accurate predictions for: / countries for Religious Practice, / countries for Supernatural Beliefs, / for Belief In God, / for Existential Security. This reflects a statistically significant superior level of accuracy for our model compared to the best of the two alternatives for each country with p < . (p = . ). This is evidence that the superior level of accuracy for our model on countries for which it was trained is more than just random noise in the output of the models. .
For the countries on which the models were not trained our model o ers more accurate predictions for: / countries for Religious Practice, / countries for Supernatural Beliefs, / for Belief In God, / for Existential Security. This does not reflect a statistically significant superior level of accuracy for our model compared to the best of the two alternatives for each country. As discussed in the Evaluation, this may be a result of over fitting the model for countries on which it was trained or under fitting the model by not including enough relevant variables. Additionally, the lack of a statistically significant superior level of accuracy may indicate that our assumption that religious practice, supernatural beliefs, belief In God and existential security are all times series that are instances of the underlying data generation process is incorrect. Each of these possible explanations will be explored in future work. .
Finally, we emphasize that we have only evaluated our ABM against two alternative approaches, neither of which is an ABM. While, the Baseline and LR models reflect accessible alternate strategies to forecasting religious practice, supernatural beliefs, belief In God and existential security, they are simple models. However, our evaluation serves as a platform for other researchers to try and construct alternative ABMs with superior accuracy in the same manner. Providing such a platform is a necessary step to creating models of complex systems that can inform policy decisions (Ahrweiler et al. ).

Conclusion
. Our work employs existing data sets and agent-based modeling to forecast changes in the religiosity and existential security among a collective of individuals over time. Our model includes agents in social networks interacting with one another based on the education level of the agents, the religious practices of the agents, and each agent's existential security within their natural and social environments. The data used to inform the values and relationships among these variables is based on rigorous statistical analysis of the International Social Survey Programme Religion Module (ISSP) and the Human Development Report (HDR).

.
Our results show that for a given country and a given time period, our model provides a more accurate forecast of changes in the existential security and the religiosity than two alternative approaches for a specific time period for specific countries. While the context in which of our model is valid is constrained and the model has a number of limitations it reflects an e ort to link empiricism and theory within an ABM. In future work, we will explore additional mechanisms that may help to further clarify the changes in religiosity observed in our model and the possible adaptive role of secularization.  E_K_SYM_LOWER_LIMIT lower limit on weight of kin relationship in sympathy group . low k,sym in (Conti et al. ) E_NK_SUP_LOWER_LIMIT lower limit on weight of non-kin relationship in support group . low nk,sup in (Conti et al. ) E_NK_SYM_LOWER_LIMIT lower limit on weight of kin relationship in sympathy group .