Exploratory and Confirmatory Factor Analyses of Religiosity. A Four-Factor Conceptual Model

We describe an exploratory and confirmatory factor analysis of the International Social Survey Programme Religion Cumulation (1991-1998-2008) data set, to identify the factors of individual religiosity and their interrelations in quantitative terms. The exploratory factor analysis was performed using data from the first two waves (1991 and 1998), and led to the identification of four strongly correlated and reliable factors which we labeled Religious formation, Supernatural beliefs, Belief in God, and Religious practice. The confirmatory factor analysis was run using data from 2008, and led to the confirmation of this four-factor structure with very good fit measures. We also ran a set of structural equation models in an attempt to determine the causality links between these four factors. It was found that for the models which provide the best fit Belief in God does not cause Supernatural beliefs, Religious practice can cause Belief in God and that there are multiple paths leading to Belief in God, most of which include Religious formation as a source. The exploratory factor analysis also led to the identification of other factors related to traditional values, confidence in institutions and influence of religious leaders on politics, but these were found to have lower reliability or insufficient number of items to meet the acceptance criteria, and thus were not included in the confirmatory factor analysis and the investigation of causal links. The results obtained in this work have important material implications for the conceptualization of"religiosity,"and important methodological implications for the scientific study of religion.


Introduction
Religion has traditionally played a central role in human societies, shaping both personal development and cultural change. The complexity and diversity of religious phenomena have led to long-standing debates among scholars about how to conceptualize and measure "religiosity." Definitional and methodological proposals vary significantly across disciplines. Some focus on the formation of beliefs about supernatural agents within individuals, for example, while others focus on the function of ritual behaviours in group cohesion. We believe that the scientific study of these highly contested and extremely complicated phenomena calls for a multidisciplinary approach encompassing fields as diverse as cognitive science, evolutionary biology, moral psychology, anthropology, sociology, and political science [1]. In this article, we describe 1/22 arXiv:1704.06112v1 [stat.AP] 20 Apr 2017 a way of conceptualizing religion in light of quantitative factor analyses of a wide range of variables culled from large scale survey data. These analyses led to the discovery of four highly correlated factors, providing a new way of measuring "religiosity" that contributes to our understanding of the interrelationships (and possible causal links) among the core components of the relevant phenomena.
Several researchers have utilized statistical analyses of large data sets based on large scale surveys such as the World Values Survey (WVS), the European Values Survey (EVS) or the International Social Survey Programme (ISSP). For example, Gervais and Najle [2] used three variables in the WVS related to belief in God, religious upbringing and attendance to services to show, via a signal detection approach, that religious upbringing has a strong influence on later belief in God. Norris and Inglehart [3] presented a conceptual scheme connecting a number of hypotheses for explaining secularization ( [3], Fig. 1.1, page 15), which they define as the erosion of religious participation, values and beliefs. They classify the latter as "indicators of religiosity" ( [3], Table 2.1, page 41). Their factor analysis in Table 7.1 of [3] describes three factors of "Work Ethic," but no attempt was made to carry out a corresponding analysis on the religiosity "indicators." Inglehart and Welzel [4,5] used sets of selected variables in the WVS, EVS and ISSP to provide supporting evidence for their theories on secularization and the human development path (both of which are related to religion). More recently, Doebler has identified three "dimensions" of religion (believing, belonging, and attendance) based on a multilevel analysis of the EVS [6].
However, none of the previous works was specifically directed towards the quantitative description of the core factors of individual religiosity, using a sufficiently representative set of variables related to attitudes and opinions about religion. Moreover, none of these studies applied the methods of factor analysis and structural equation modeling (SEM) in a systematic way. For example, the three indicators of religiosity considered by Norris and Inglehart [3] were not proven to be factors of religiosity derived from statistical methods. Another problem is that even when factor analyses are used in the scientific study of religion, scholars rarely present detailed descriptions of the criteria used for acceptance or rejection of items and factors, the sequence of the elimination process, or estimates of reliability and goodness of fit for the solutions obtained.
In this article, we describe a four-factor conceptual model of religiosity based on statistical analyses of 34 selected variables in the ISSP Religion Cumulation data set [7], which combines the results from three waves of questionnaires on religion for years 1991, 1998 and 2008. The work was performed in three stages. First, we ran an exploratory factor analysis (EFA) using the first two waves (1991 and 1998). In this way we identified four reliable and strongly correlated core factors of religiosity, each with a clear meaning. These are "Supernatural beliefs," "Religious formation," "Belief in God," and "Religious practice." Then, we used the last wave (2008) to run a confirmatory factor analysis (CFA), as EFA and CFA should be based on different data [8]. This led to the confirmation of the results obtained in the first stage. Finally, we ran 588 different structural equation models (SEMs) in an attempt to determine the causal relationships between the four factors found in the previous stages, using the data from the 2008 wave (the same used for CFA).
The main contributions of the present work for the scientific study of religion are i ) the identification of four core factors of religiosity and their quantitative measurement based on a large set of variables and number of respondents (i.e. using the individual records instead of aggregating results by country or religious affiliation); and ii ) the systematic description of the statistical procedures used and the presentation of quantitative estimates of factor reliability and models' goodness of fit (GOF).

Description of the ISSP Religion Cumulation Data Set. Selected Variables
The ISSP Religion Cumulation data set [7] contains the cumulated variables of the ISSP "Religion" surveys of 1991, 1998 and 2008 and comes in two separate files: a main file (ZA5070) with topic-related and background variables that appear in at least two surveys, and an add-on file (ZA5071) with variables that could not be cumulated for various reasons. The analysis in this article is based on the information in main file, which includes 122 variables (fields) for 102454 respondents from 28 countries. Only the countries that participated in at least two surveys (1998 and 2008) were included in the cumulation data file. For Germany and the United Kingdom two subsamples were included (East and West German, and Great Britain and Northern Ireland, respectively). Details on the contents, structure and coding of the ZA5070 cumulation file can be found in [9].
Because of its size, the variables included and the proportion of valid respondents' answers, the ISSP Religion Cumulation data set is very valuable for extracting information about religious beliefs, practices and values. However, the ISSP "Religion" surveys have the important limitation that they mostly cover Christian Catholic, Protestant and Orthodox countries. Although Israel and Japan also contributed, there is no information from Muslim countries, India or China. Table 1 shows the initial set of questions selected for the exploratory factor analysis as well as the variables' labels, number of valid levels and indication of polarity inversion. First, we identified the questions about religious beliefs, practices, and values that were most clearly related to the scientific study of religion. Variables V11-16 tap attitudes and values that are often explored in relation to religiosity in the literature. Variables V20-24 were included to study the correlation between religiosity and confidence in institutions, and identify eventual differences between confidence in religious and secular institutions. This issue is important for the study of secularization. Variables V25-27 are related to the relationship between religion and politics, which we also took as possible proxies for secularization. Variables V28-33, V35, V37 and V51 are related to beliefs in the "supernatural," particularly in God (V28-29 and V35). Variables V46-48 are related to the attendance of the respondent's mother and father, and self-attendance during a formative period (ages [11][12], which is often considered to be particularly influential on religious belief and practice later in life [2,10,11]. Variables 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 current religious practice. Because science (which is methodically naturalistic) is often in competition with supernatural beliefs, we also included variables V64-65. Finally, we selected demographic variables that could help exploring other aspects often explored in the scientific study of religion such as age, gender and educational level. For instance, sex was expected to be important since previous studies reported that women tend to be more religious than men [12]. Variables such as religion of origin or employment status were not included because they have no clear ordinal meaning.

Software Tools
The data processing was done using R [13]. Functions were written to import the ZA5070 v1-0-0.dta Stata file and generate a R data frame, and to perform frequent queries and operations (retrieving a set of fields corresponding to selected questions, the records for a particular year and country or set of countries, etc.).
Exploratory factor analysis was done using the psych package [14], which includes a Highest.education.level.degree 6 No * The uppermost level "No mother/mother not present" was merged with the lowest level "Never." ** The uppermost level "No father/father not present" was merged with the lowest level "Never." large number of functions for doing principal component, factor and cluster analysis of items, determining the number of factors to extract, estimating factor reliability and graphically representing the solutions. Confirmatory factor analysis of the four reliable and strongly correlated factors of religiosity found in the exploratory factor analysis was done using lavaan, which a free open-source, but commercial-quality package for latent variable modeling [15]. The structural equation modeling (SEM) analysis was done using lavaan as well.

Data Preparation
The R data frame generated from the ZA5070 v1-0-0.dta Stata file was split in two data frames, one including the data from years 1991 and 1998 with 61928 records and another the data from year 2008 with 40526 records, so as to run EFA and CFA with distinct data.
The data preparation was done in the following steps (identical for EFA, CFA and SEMs).
First the selected variables were converted to ordered factors except AGE and SEX which were converted to numeric and unordered factor (categorical, nominal) respectively. Then, all ordinal variables were inspected to mark invalid/inconclusive answers as NA and remove spurious factor levels (e.g. levels for which no answer was recorded. After that, the ordinal variables were further inspected for the convenience of polarity inversion (reversal of the levels' order). More specifically, we wanted the top levels to represent the highest degree of beliefs in the "supernatural," religious practice, traditional values, God giving life a meaning, favorable to religion influencing politics, unfavorable to modern science, and confidence in institutions (both religious and secular). In this way, all item loadings and correlations between factors were expected to be positive. This facilitates the interpretation of factor solutions, computation of factors' reliability and scores, and construction of summed scales. Table 1 shows the questions for which the polarity was inverted. Variables V46-47 further required merging of top to bottom levels, so that all cases of "no attendance" were treated in the same way.
In the EFA stage we computed the correlation matrices using three different methods via the psych::mixed.cor function: Pearson (simplest), Spearman, and mixed Spearman/polychoric with polychoric correlations being used for variables with up to six levels (inclusive). It was found that the differences between the correlation matrices computed using the three methods were negligible and Pearson correlation produced marginally better factor solutions. 1 Therefore, all variables were converted to numeric in a final stage, which allowed faster and more efficient computation of EFA, CFA and SEMs.

Exploratory Factor Analysis Method
The EFA was done in two stages. In the first stage we computed factor solutions and eliminated variables (items) with low communality and low or cross-loadings [8], until all variables met the acceptance criteria described below, and a statistically significant and theoretically meaningful solution was found. This process was iterative, for removing variables also led to changes in the optimal number of factors to be extracted and consequently the variables' communality and loadings. In the second stage, the factors' reliability was assessed using the criteria described below, and a factor solution was computed using only the variables loading on the factors deemed reliable. This final factor solution was found to meet all acceptance criteria and to be theoretically meaningful, and the ensuing factor structure was then considered in the CFA and SEM analyses.
The first stage consisted of the following steps: 1. Inspection of the correlation structure and estimation of the optimal number of factors to be extracted; 2. Computation of the factor analysis solution for the recommended number of factors as well as for one more and one less factor than recommended (i.e. "bracketing" on the number of factors) [8]; 3. Inspection of the factor solution for variables with low communality and weak or cross loadings [16], [8]; 4. Elimination of variables that did not meet the acceptance criteria, if any, and re-doing the previous steps; otherwise, proceed to the second stage (inspection of factors' reliability).
To determine the recommended number of factors, we used the scree test, the Very Simple Structure (VSS) and Velicer's Minimum Average Partial (MAP) criteria, as well as parallel analysis using psych::fa.paralell [14]. Parallel analysis consistently led to better estimates for the number of factors than the other methods, and was therefore the method of choice in the present work (bracketing on the number of factors as described in the second step above also confirmed this finding).
The acceptance criteria for the variables were: i ) communality > 0.40 [16] (reference [8] recommends > 0.50); ii ) loadings with absolute value > 0.32 (accounts for at least 10% variance); and iii ) no cross loadings (variables with loadings with absolute value at 0.32 or higher on two or more factors [16]).
In the second stage, the factors obtained at the end of the first stage were inspected for reliability, using Cronbach's α [17] computed via the psych::alpha function. The acceptance criteria for factors were α > 0.70 and at least three variables (items) loading on them [8]. This last restriction is important for the subsequent CFA. Factors that failed to meet these criteria were removed from the analysis by eliminating the variables loading on them, and a final factor solution was computed and checked against all acceptance criteria (for variables and factors).
We now describe further details necessary for the replication of the EFA reported herein. The factor solutions were computed using psych::fa using data frames as inputs (instead of correlation matrices). Correlations were estimated using the Pearson coefficient (cor = ''cor''), as noted in the comment above on data preparation, with pairwise complete observations and no imputation of missing values. Oblique rotations (using the oblimin method, default in psych::fa) were used, since strong correlations between factors were to be expected. 2 Communality estimates were computed using the minimum residual method [18,19](default in psych::fa), which does note require matrix inversions and generally improves over the principal axis factor (PAF) method [14]. We used one hundred bootstrap iterations and specified a maximum of 2000 iterations for convergence. The α level for confidence intervals of loadings and factor correlations was set to 0.01. The remaining input parameters required by psych::fa were set to default values.

Confirmatory Factor Analysis Method
We used CFA to test the measurement and structural model of the four-factor solution suggested in EFA, i.e. to determine how well the measured variables represent the factors and the correlation structure among the factors, respectively. We used the 2008 data to run the CFA, since EFA and CFA should be done using different data sets [8]. Admittedly, the variables' loadings and correlations between factors could have changed due to social, cultural, political or economic changes. But confirmation of the model suggested by EFA implies that the factor structure is "solid," in the sense that it both fulfilled the statistical requirements and withstood the "test of time." For these reasons, we preferred to use the 2008 part of the data set (the latest and largest of the three waves) for CFA instead of random samples extracted from the three years.
The CFA was done by first preparing the data as done for EFA. Then, the measurement and structural models of the final four-factor solution obtained in EFA were set up in and run in "lavaan" using lavaan::cfa. Following the standard procedures in CFA, no cross loadings were permitted (measured variables were forced to load on a single factor), and non null correlations between factors were allowed. The GOF of the confirmatory model was assessed using several different estimators (as described below) and the acceptance criteria were [8]: i ) standardized loadings with (absolute) value > 0.50; ii ) average variance extracted (AVE) > 0.50; iii ) reliability of constructs (factors) > 0.70; and iv ) standardized residuals of measured variables < 2.50.
The details necessary for the replication of our CFA are as follows. The lavaan::cfa function was ran using a data frame of the variables loading on the four factors obtained in EFA with the records for ISSP Religion 2008. The model was specified so that each variable (item) loaded on a single factor. The scale of the latent variables was set by fixing the factor loading of the first indicator to one (default in lavaan::cfa). The estimator used was diagonally weighted least squares (DWLS) [20]. One hundred bootstrap iterations were used, as was done in EFA. The remaining parameters required to run lavaan::cfa were set to their default values.

Structural Equation Modeling
Next, we used structural equation modeling (SEM) to organize the relationships among the four factors obtained in EFA. SEM enables us to specify path diagrams that hypothesize relationships among the four factors. Each diagram can be converted into a model with two components [15]: 1. a measurement model that describes the relationship of factors and their indicators (survey questions) 2. structural equations that depict regressions among factors as causal paths from one factor to another.
Given a model, the structural equation modeling algorithm (lavaan::sem) identifies the values for parameters within the measurement model and structural equation model that best match the observed data. The best parameters were estimated using Full Maximum Likelihood Estimation (FML). 3 The χ 2 is the primary fit index for assessing the GOF between the observed values in the data and those produced by a SEM model. A χ 2 that is not significant implies a model with good fit. However, achieving a χ 2 that is not significant for any model, regardless of fit quality, becomes increasingly difficult as the size of the sample increases. Most SEMs are based on samples of size200. However, the sample size in our data set is measured in tens of thousands. Given the size of our data set, identifying a χ 2 that is not significant is exceptionally unlikely. This is not uncommon. In situations where the sample size is sufficiently large, four fit indices are employed in combination to evaluate SEMs. These four fit indices are split into two categories: (1) incremental fit indices and (2) absolute fit indices [21][22][23][24].
An incremental fit index measures the fit of a given SEM with respect to the SEM that has the worst possible fit and the best possible fit. This places the fit of a SEM on a continuum from 0 (worst fit) to 1 (best fit). The two incremental fit measures employed in the evaluation of SEMs are the Comparative Fit Index (CFI) and Tucker-Lewis Index (TLI). Both apply penalties for including additional parameters in the model, but the TLI applies a more severe penalty than the CFI. Values of CFI and TLI > 0.95 reflect models with acceptable fit [21]. The TLI is also referred to as Non-Normed Fit Index or (NNFI).
An absolute fit index presumes that the best fitting model has a fit of zero, and the corresponding measure of fit expresses the model's deviation from perfect fit. We used two absolute fit indices, namely the Root Mean Squared Error of Approximation (RMSEA) and the Standardized Root Mean Square Residual (SRMR). RMSEA is the standardized difference between the observed data and the SEM predicted data, while SRMR is defined as the standardized difference between the correlation of the observed data and the predictions of the SEM. Values of SRMR and RMSEA < 0.08 reflect models with acceptable fit [21].
The details necessary for the replication of our SEM are as follows. The SEM was done by first preparing the data as done for EFA and CFA. Then 588 different candidate models were generated. Each model was fit using a data frame of the variables, loading on the four factors obtained in EFA, with the records for ISSP Religion 2008. The lavaan::sem function was run on the data with all function parameters set to their default values. The quality of fit was assessed using the four fit indices previously mentioned: CFI, TLI, RMSEA and SRMR.

Exploratory Factor Analysis
As outlined above, EFA was done in two stages. The first stage started by computing factor solutions using all variables in table 1 for a recommended (and best) twelve factors. Examination of this solution according to the criteria outlined above led to the exclusion of variables "V12," "V16," "V21," "V24," "AGE," "SEX" and "DEGREE." In a second iteration, factor solutions were computed excluding these variables for a recommended (and best) ten factors. This led to a third iteration in which variables "V35," "V37," and "V64" were excluded and factor solutions were computed for a recommended (and best) eight factors. Examination of the best eight-factor solution led to further elimination of variables "V27" and "V65," the first due to low communality and the second because it did not load on any factor. The final iteration in stage one consisted of computing factor solutions for a recommended eight factors, which led to the best fit indexes. Fig 1 shows the factor solution diagram for the eight-factor solution obtained at the final iteration of the first stage of the EFA. Table 2 shows the reliability estimates for the eight factors for several different estimators, computed using psych::alpha.
All factors in Fig. 1 are seen to have a clear meaning. MR1 can be labeled "Supernatural beliefs;" MR2 "Religious formation;" MR3 "Religion and politics" (influence of religious leaders on vote and government); MR5 "Belief in God;" MR6 "Confidence in institutions;" MR7 "Religious practice;" and the two factors MR4 and MR8 are related to "traditional values" (or "Conservatism") with respect to opinions about abortion and sex, respectively. 4 It can be observed that factors related to religion and traditional values are significantly correlated with two or three other factors, whereas "Confidence in institutions" and "Religion and politics" are weakly correlated (e.g. nearly orthogonal) to all other factors. Inspection of  acceptance criteria stated above. Also, they are also the ones with top ratio of reliable to unreliable variance (S/N). Because of their clear meaning, strong interrelation and reliability, these four factors emerge from the analysis as representing the "core dimensions" or individual religiosity. Therefore, a final factor solution was computed using only the variables loading on these factors, to obtain final estimates for the loadings and factor correlations as well as fit measures, before stepping to CFA.
The fact that confidence in institutions and opinions about religious leaders influencing politics are weakly correlated with religiosity has interesting implications for e.g. the study of secularization. Likewise, it is interesting to question why opinions about abortion and sex (particularly homosexual sex, which Inglehart et al. [5] took as a sensitive indicator of tolerance to out-groups) came out as independent factors. These aspects will be briefly addressed in the discussion below, but not further pursued in the subsequent analyses. Fig. 2 shows the diagram for the four-factor model of religiosity. Table 3 shows the standardized factor loadings, communality, uniqueness and Hoffmann complexity index [25,26] for all variables, and table 4 shows the corresponding estimates and 99% confidence intervals for the factor correlations. Description: raw α is Cronbach's estimate based on covariances; std. α is Cronbach's estimate based on correlations; λ 6 is Guttman's λ 6 reliability coefficient; average r is the average inter item correlation; S/N is the signal-to-noise ratio (or ratio of reliable variance to unreliable variance); ase is α's standard error; mean is the mean of the summed scale constructed from the items; and sd is the standard deviation of the total factor scores [14]. Factor labels: MR1 ⇐⇒ "Supernatural beliefs;" MR2 ⇐⇒ "Religious formation;" MR3 ⇐⇒ "Belief in God;" MR4 ⇐⇒ "Religious practice." It can be observed that the model is strongly consistent with the specified acceptance criteria, and that the GOF estimates are very good. The factor correlations also provide interesting information. The strongest correlation is between "Supernatural beliefs" and "Belief in God," and its 99% confidence interval does not intersect those of the other factor correlations. Next come the correlations between "Supernatural beliefs" and "Religious practice," and between "Belief in God" and "Religious practice," whose 99% confidence intervals overlap. At lower and disjoint 99% confidence intervals come the correlations between "Religious formation" and "Religious practice," "Religious formation" and "Belief in God," and "Supernatural beliefs" and "Religious formation."      12/22 section above. The CFA was computed using lavaan:cfa, which converged after 76 iterations, using 23636 out of 40526 observations. From the additional information presented in Fig. 3 it can be concluded that the model has very good fit estimates, which are consistent with those obtained in the EFA (e.g. the RMSEA fit index). Comparing Figs. 2 and 3, it can be concluded that in the CFA model both items' loadings and factor correlations have higher values than for the final model obtained in the EFA. Also, the differences between the (absolute) values of the items loading on each of the factors are smaller in the CFA model (particularly for the "Belief in God" factor). Additional information: α is Cronbach's reliability coefficient. The theoretical description of the ω 1 , ω 2 and ω 3 coefficients can be found in [27], [28] and [29] respectively. Table 5 shows some reliability estimates for the CFA model. Combining the information presented in Fig. 3 and table 5, it can be observed that the factors' reliability and AVE meet the criteria stated the the "Methods" section. Therefore, the four-factor model was further studied using SEMs to investigate how the fit measures depend on other types of links (e.g. regression) between factors, in an attempt to identify statistically robust causal relationships among them.
Given our four factors it is possible to construct 588 unique path diagrams. Each diagram can be converted to a SEM. Recall, while the χ 2 is the primary fit index for assessing the difference between the observed values in the data and those produced by a SEM, it is not useful for our data set due to its size (> 40,000 samples). As a result four fit indices are used to assess our models: CFI, TLI, SRMR and RMSEA. Values of CFI and TLI > 0.95 and values of RMSEA and SRMR < 0.08 reflect models with acceptable fit [21][22][23][24].
We tested each of the plausible 588 SEMs that can be formed from the four strongly factors given our knowledge of the domain. Our goal was to determine the extent to which each plausible model fit our data. Fig 4 shows the distribution of the four fit indices for the candidate models.
Of the 588 candidate models, there are four that share the same fit indices and exhibit the lowest SRMR (0.0344) and RMSE (0.0735) and the highest CFI (0.9675) and TLI (0.9575). Each of these measures is below the specified threshold for the fit of the models to be considered acceptable. In addition, the models do not have invalid parameter values (negative residuals or variances or correlations larger than 1.0). Furthermore, the regression coefficients, factor loadings and covariances all have the expected sign (positive). Finally, the supplemental information associated with the paper shows that for these four models all free parameters are statistically significant and there are no excessively large standard errors [21][22][23][24]. These four models are shown in Figure 5. The number of parameters ( # Pars), degrees of freedom (DF) and four fit indices for the four models are provided in Table 6. Table 7 provides several additional fit indices for the models including the Bentler-Bonett Index or Normed Fit Index (NFI) and upper and lower confidence intervals for the RMSEA. All of these four models share three common properties. First, "Belief in God" (MR3) does not play a role in causing one's "Supernatural beliefs" (MR1). In two of the four models "Supernatural beliefs" (MR1) can cause one's "Belief in God" (MR3) and in two of the models there is only a covariance between the two factors, not a causal relationship. As far as we know we are the first researchers to find quantifiable evidence that a causal relationship from "Belief in God" (MR3) to "Supernatural beliefs" (MR1) does not exist. We discuss the implications of this finding further in the Discussion Section. Second, in each of the four models "Religious practice" (MR4) can directly cause "Belief in God" (MR3). The relationship between these factors has been a frequently discussed topic in the scientific study of religion. Several pieces of research into the causal ordering of these factors mirror our results, while several other pieces conflict with our findings. We discuss those findings with respect to our work in the Discussion Section.
Finally, Figure 5 demonstrates that there are multiple paths to "Belief in God" (MR3) in each of the four models. Every model includes a direct path from "Religious formation" (MR2) to "Belief in God" (MR3). However, in each model at least one more additional path is possible and in one of the four models there are five different paths, originating at three different factors, to "Belief in God" (MR3). We discuss this evidence for multiple paths to "Belief in God" (MR3) with respect to existing research in the Discussion Section.
Based on the evaluation measures shown in Table 6 and 7 for the SEMs, we cannot express a further preference for any one of these four models. Instead they reflect four equally acceptable theories of the organization of the four factors identified in our EFA. However, the existence of the three properties shared among the models serves as a template for the organization of these factors given the data.

Discussion
The four-factor model that emerged out of the EFA and CFA described in this work is significant in conceptual terms, because it bears on one of the most contentious arguments among the scholars of religion, namely how (or whether) to define "religion." Some authors who analyze this problem from the perspectives of psychology, anthropology, history, or sociology of religion prefer definitions that incorporate variables related to "values" and "culture," which are not per se related to the "supernatural." In contrast, the elimination process carried on in the EFA and confirmed via the CFA left us with only four reliable and strongly correlated factors -"Supernatural beliefs," "Religious formation," "Belief in God," and "Religious practice"all of which are related to supernatural beliefs and practices. This complements the findings of other factor analyses, e.g. of psychological measurement scales, which have identified "supernatural belief/practice" as the "only unique factor" associated with "religiosity" [40] and "spirituality" [41]. It is also more consistent than the indicators proposed in [3].
We also provided significant statistical evidence that the correlation between the two factors related to beliefs is stronger than the correlations between beliefs and religious practice, and between religions formation and supernatural beliefs (table 4). Also, it is not surprising that a strong correlation was found between factors such as belief in the supernatural (and God) and attendance at religious services. The relationship between the tendency to believe in disembodied intentional forces and the tendency to signal commitment to an in-group through ritual participation is extremely well documented in fields such as the cognitive science of religion and cultural anthropology [1,[30][31][32][33][34][35]. However, our results were novel and somewhat surprising in several ways. For example, variables related to traditional values (namely opinions about sexual behaviors and abortion) formed factors that, although correlated with the four core factors of religiosity, had to be discarded in the CFA and SEM analyses because of insufficient reliability or number of items. One possible explanation for this surprising result is that the traditional values indicated in MR4 and MR8 in Fig 1 are more superficially related to religion than those that survived the winnowing process. Sexual selection in earlier ancestral environments would have been more closely tied concerns about family cohesion and child-bearing [36,37], but the values reflected in V11, V13, V14, and V15 are arguably more reliant on later forms of social entrainment associated with modern cultures. Another possible explanation is that the social changes associated with the shift from survival values to self-expression values and increasing emancipation from authorities (religious or secular) [4] in the countries included in the ISSP Religion Cumulation data set, have led to a dissociation between religiosity and attitudes towards behaviors such as abortion, adultery, and homosexuality. Still another possible explanation is that the number of variables included in the ISSP Religion Cumulation data set is insufficient to extract reliable factors related to values. To clarify these issues, we plan to analyze other data sets such as the WVS or the EVS.
Another unexpected and interesting result was the weak correlation between factors related to attitudes toward the political influence of religious leaders and confidence in institutions, and the four core factors of religiosity. This has important implications for the study of secularization. For instance, if we were to define "secular" in terms of high confidence in public institutions to guide society, or the desire to separate church and state, then the overall population of all the countries represented in the ISSP Religion Cumulation data set is already secular! This suggests that "religiosity," as understood by most of the respondents, has less to do with political organization and trust in institutions than with beliefs in supernatural agents and participation in supernatural rituals. In other words, our findings suggest that neither confidence in public institutions nor lack of confidence in religious institutions is a good proxy for secularity.
In the SEM analysis, we tested different combinations of causality links between the four factors of religiosity. It was found that for the four models with the best fit indices considered in the analysis, "Belief in God" is the common sink into which the remaining factors flow, and thus never causes "Supernatural beliefs". Also, "Religious practice" always causes "Belief in God," although its effect on later beliefs is mainly indirect. However, the differences between the measures of fit for these four models were not statistically significant, and therefore we cannot decide, based on statistical considerations alone, which one is best. Nevertheless, all the four models with best fit are consistent with previous research. It has been shown by other authors that the tendency to perceive disembodied intentional forces (expressed by factors "Supernatural beliefs" and "Belief in God") and the tendency to participate in in-group rituals (expressed by factors "Religious formation" and "Religious practice") are reciprocally reinforcing. Exploring the precise connections between these tendencies (and the mechanisms that produce them) requires a multi-disciplinary approach that attends to a variety of methodologies [1,48].
The causal precedence of "Supernatural beliefs" with respect to "Belief in God" is also consistent with other findings in the evolutionary theories of religion. Belief in the (miraculous) intervention of supernatural agents who have the power to reward or punish (heaven or hell in an afterlife) appears to be a relatively ancient phylogenetically evolved default. The role of cognitive and coalitional biases about disembodied, potentially punitive intentional forces in strengthening the cohesion of religious in-groups is well-documented in the literature [49][50][51]. "Belief in "God," however, is far more recent, having emerged only in the wake of the axial age (c. 800-200 BCE). When humans lived in small-scale, hunter-gatherer societies, the threat that relatively localized and familial supernatural forces (animal-spirits or ancestor-ghosts) might be watching was enough to keep everyone in line and enhance cooperation and commitment. As groups grew, however, so did the size and power of their gods. Although scholars disagree on precisely what factors were most decisive in the emergence of cooperation within large-scale societies, there is general consensus that belief in one "God" (or an Ultimate Reality such as Dharma or the Dao) -who is watching over, and capable of punishing everyone, regardless of their in-group -is correlated with life in contexts strongly influenced by the so-called "axial" traditions [52][53][54][55]. This makes it easier to understand why "Supernatural beliefs" (MR1) would be conditionally independent of "Belief in God" (MR3) and not vice versa. This also helps to explain why the SEMs with best fit indices had multiple paths to MR3 all of which included MR4 (religious practice), since the latter would arguably be required to facilitate acceptance of a maximally counter-intuitive concept (God) that is 17/22 not as easily delivered by cognitive defaults as, e.g., belief in an afterlife [1].

Conclusion
The theoretical definition of individual religiosity and its underlying dimensions is an important yet open problem in the scientific study of religion. In this work, we presented quantitative evidence that individual religiosity has four strongly inter-correlated factors -"Supernatural beliefs," "Belief in God," "Religious formation" and "Religious practice" -based on EFA and CFA of the ISSP Religion Cumulation data set. The EFA began with 34 selected variables using the data for the 1991 and 1998 waves, and ended with 13 variables loading onto the four factors identified above. That final model was tested using CFA based on the data for the 2008 wave, confirming that the measurement and structural models suggested by the EFA provided a very good fit to the data. The fact that CFA was done using data that was collected ten years apart from the data used for EFA indicates that although religiosity may be affected by social, cultural and political changes, the four-factor model not only resists the sieve of statistical methods but also stands the "test of time." We plan to confirm this as soon as a new wave of the ISSP Religion questionnaire becomes available, using the methods described in the present article.
After doing the EFA and CFA of the ISSP Religion Cumulation data set, we also performed a systematic exploration of the four-factor model using SEMs. We investigated how different combinations and types of links between factors (e.g. correlation or regression) affected the goodness of fit. In the SEMs with the best fit to the data, "Belief in God" never causes"Supernatural beliefs," whereas "Religious practice" always does. Three different paths leading to "Belief in God" were identified, all of which include "Religious formation," but the direct link from the latter to "Belief in God" is the weakest in the corresponding SEMs. This suggests that the effect of "Religious formation" on later religious beliefs is mainly indirect.
To the best of our knowledge no similar studies to those presented here have been published, and thus no independent validation of the factor structure obtained in the current work could be conducted. However, to further validate the four-factor model of religiosity, and determine whether or not other reliable factors related to values can be extracted, we plan to extend the analysis reported here to the upcoming wave of the ISSP Religion questionnaire. This will also further clarify the relationships and possible causal links among the factors. In order to tap other factors that are not available in ISSP, such as existential security and self-expression values, we also plan to perform similar analyses using the WVS and EVS longitudinal data sets.
The work reported here is part of the "Modeling Religion in Norway" (MODRN) project, whose overall goal is to construct computer models that can simulate changes in religiosity and secularity in order to provide tools for clarifying and evaluating public policy proposals. We plan to utilize the data reported here to guide the development, parametrization and calibration of agent-based models that can illuminate the mechanisms at work in religious and secular change in contemporary societies.