Abstract
We investigate the causal effects of religious service attendance on prosocial behaviours using longitudinal data from a nationally representative sample of 33,198 New Zealanders collected between 2018 and 2021. Our study innovates in three ways: (1) we use longitudinal rather than cross-sectional data; (2) we incorporate measures of help received alongside self-reported giving; and (3) our statistical models are designed to address causal questions, rather than simply to describe change over time. We model causal contrasts for three hypothetical interventions – increasing, decreasing, or maintaining religious service attendance – and assess their effects on eight distinct prosocial domains. Study 1 focuses on self-reported charitable donations and volunteering. Studies 2 and 3 examine receiving help – both personal and financial – from family, friends, and the wider community. Across all analyses, we find that the causal effects of religious service attendance are notably smaller than cross-sectional correlations suggest. However, even modest increases in regular attendance would result in charitable donations equivalent to approximately 4% of the New Zealand Government’s annual spending – a considerable public benefit. By applying robust methods for causal inference to national-scale panel data, our study provides insights that can inform public policy about the social functions of religious participation and advances a methodological framework for investigating the social consequences of cultural practices.
Keywords
Introduction
A central question in the scientific study of religion is whether religion fosters prosociality (De Coulanges, 1903; D. D. Johnson, 2005; Norenzayan et al., 2016; Schloss & Murray, 2011; Sosis & Bressler, 2003; Swanson, 1967; Watts et al., 2015, 2016; Wheatley, 1971; Whitehouse et al., 2023). However, quantifying causal effects for religion presents significant challenges (Major-Smith, 2023). Investigators have limited scope to randomise supernatural beliefs, community worship, and personal prayer. On the other hand, valid causal inferences from non-experimental, ‘real-world’ data require the integration of high-resolution, repeated-measures time-series data with appropriate causal inference methodologies. Few studies succeed in this integration. A recent survey of the religion and prosociality literature reveals that nearly all non-experimental studies examining the links between prosociality and religion, including several longitudinal studies, remain correlational (Kelly et al., 2024). Currently, studies using longitudinal (panel) data have not leveraged their repeated measures to derive reliable causal inferences.
Here, to obtain causal inferences from time-series data, we utilise comprehensive panel data from 33,198 participants in the New Zealand Attitudes and Values Study (NZAVS) spanning 2018 to 2021. We quantify the effects of clearly defined interventions in religious attendance across the population of New Zealand, assessing eight outcome domains related to charitable financial donations and volunteering. We obtain causal inferences by contrasting expected population averages under different modified treatment policies (Díaz et al., 2023; Haneuse & Rotnitzky, 2013; Hoffman et al., 2023).
Our initial causal contrast investigates: ‘What would be the average difference in each of the eight prosocial outcomes across the New Zealand population if everyone attended religious services regularly (at least four times per month) compared to if no one attended?’ This theoretical inquiry simulates a hypothetical experiment with random assignment to either regular attendance or complete non-attendance, mirroring the all-or-none contrast that is commonly used in experimental designs (Hernán et al., 2016).
A second causal contrast investigates: ‘What would be the average difference in each of the eight prosocial outcomes across the New Zealand population if everyone attended religious services regularly (at least four times per month) compared to maintaining current attendance levels?’ In this analysis, we compare the expected prosocial outcomes under regular religious service attendance (at least four times per month) with outcomes observed under existing attendance patterns. This causal contrast may inform practical policies aimed at influencing non-regular attendees who might stop attending religious services.
Our third causal contrast examines: ‘What would be the average difference in each of the eight prosocial outcomes across the New Zealand population if no one attended religious services compared to maintaining current attendance levels?’ In this analysis, we compare the expected prosocial outcomes in a scenario where all religious service attendance ceases with the outcomes observed under existing attendance patterns. This causal contrast may inform practical policies aimed at influencing regular attendees who might discontinue their participation in religious services.
Although the set of causal contrasts investigators might consider is theoretically limitless, the contrasts that we have selected address our targeted scientific and policy interests.
Importantly, our approach does not centre on testing specific hypotheses; rather it is designed to accurately and consistently compute our stated causal contrasts (Hernán & Greenland, 2024).
Method
Sample
Data for this study were collected as part of the New Zealand Attitudes and Values Study (NZAVS), an annual longitudinal national probability panel that assesses New Zealand residents’ social attitudes, personality, ideology, and health outcomes. The study has obtained responses from 72,910 participants since it started in 2009. The study operates independently of political or corporate funding and is based at a university. Data summaries for all measures used in this study are provided in Supplemental Appendices B–D: https://osf.io/vxd6m/. For more information about the NZAVS, visit: OSF.IO/75SNB.
The data analysed in this study were obtained from the NZAVS waves 10–12 cohort, covering the years 2018–2021. This is the largest cohort in the NZAVS. Although this cohort was obtained through national probability sampling and closely approximates the New Zealand population, we applied weights based on the 2018 New Zealand Census estimates for age, gender, and ethnicity to further enhance representativeness. Detailed information on these survey weights is available in the NZAVS documentation at OSF.IO/75SNB.
Treatment indicator
We assessed religious service attendance using the following question:
Do you identify with a religion and/or spiritual group? If yes, how many times did you attend a church or place of worship during the last month?
Responses were rounded to the nearest whole number. Because few participants reported attending more than eight times, we capped responses above eight at eight (refer to Online Supplement Appendix B: https://osf.io/vxd6m/).
Measures of prosociality
We assessed prosocial behaviour using the following measures:
Study 1: self-reported charity
Participants reported:
Volunteering: Please estimate how many hours you spent doing each of the following things last week . . . Volunteer/charitable work.
Annual Charitable Financial Donations: How much money have you donated to charity in the last year?
Study 2: help received from others in the last week: time
Participants estimated the amount of help they received in the past week in hours from:
Family . . . TIME (hours)
Friends . . . TIME (hours)
Community . . . TIME (hours)
Owing to high variability, we transformed responses into binary indicators: 0 = none; 1 = any.
Study 3: help received from others in the last week: money
Similarly, participants were asked:
Participants estimated the monetary help they received in the past week from:
Family . . . MONEY (dollars)
Friends . . . MONEY (dollars)
Community . . . MONEY (dollars)
Again, owing to high variability, responses were converted to binary indicators: 0 = none; 1 = any.
Studies 2 and 3 employ revealed measures of prosocial exposure to minimise self-presentation bias. We assume that if religious institutions foster prosociality, initiating regular attendance would increase exposure to prosocial behaviours. This approach allows us to infer prosociality using help-recieved responses as a measure. The help-received measure is less susceptible to self-presentation biases that might otherwise distort the association between religious attendance and genuine prosociality. Furthermore, including baseline outcomes and treatments in our analyses mitigates the threat of undirected correlated errors.
Comprehensive details of all measures are provided in Online Supplemental Materials: Appendix A: https://osf.io/vxd6m/.
Causal interventions
We defined three targeted causal contrasts based on pre-specified modified treatment policies:
Regular Religious Service Treatment: assumes that everyone attends religious services regularly, defined as at least four times per month. In this scenario, individuals who currently attend less than four times per month have their attendance increased to four times, while those already attending four or more times per month remain unchanged.
Zero Religious Service Treatment: assumes that no one attends religious services. Here, individuals who currently attend more than zero times per month have their attendance reduced to zero, whereas non-attendees remain unchanged.
The Status Quo: no treatment is applied. We retain each individual’s observed level of religious service attendance without modification.
Causal contrasts
Based on these policies, we computed three causal contrasts:
Regular vs Zero Religious Service Attendance: this contrast compares average prosocial outcomes in a society where everyone attends religious services regularly to one where no one attends. It simulates a hypothetical experiment where individuals are randomised to either regular attendance or no attendance, allowing us to assess differences in prosociality one year after the intervention.
Regular Religious Service Attendance vs The Status Quo: this contrast compares average prosocial outcomes in a society where everyone attends religious services regularly to the current state. It contrasts the causal effect of transitioning non-regular attendees to regular attendance.
Zero Religious Service Attendance vs The Status Quo: this contrast compares average prosocial outcomes where no one attends religious services to the current state. It examines the causal effect of entirely eliminating religious service attendance from society.
Identification assumptions
To consistently estimate a causal effect, investigators must satisfy three assumptions (refer to Bulbulia, 2024b):
Causal consistency: potential outcomes correspond to the observed outcomes under the treatments in our data. We assume that conditional on measured covariates, potential outcomes do not depend on how the treatment was administered (VanderWeele, 2009; VanderWeele & Hernan, 2013).
Conditional exchangeability: conditional on observed covariates, treatment assignment is independent of the potential outcomes being compared (i.e. there is no unmeasured confounding) (Chatton et al., 2020; Hernan & Robins, 2024).
Positivity: every individual has a non-zero probability of receiving each treatment level, regardless of their covariate values. We evaluated this by examining changes in religious service attendance from baseline to the treatment wave (Westreich & Cole, 2010).
Target population
Our target population comprises New Zealand residents represented in the baseline wave of the NZAVS during 2018–2019, weighted by 2018 New Zealand Census data for age, gender, and ethnicity (Sibley, 2021). Although the NZAVS is a national probability study designed to reflect the broader New Zealand population, it tends to under-sample males and individuals of Asian descent and over-sample females and Māori (the indigenous people of New Zealand). To address these disparities and enhance the accuracy of our findings, we applied survey weights adjusting for age, gender, and ethnicity. Survey weights were integrated into our statistical models using the weights option in the ‘lmtp’ package (Williams & Díaz, 2021), following protocols described by Bulbulia (2024e). Note that the sample in this study is a single cohort, enrolled in the NZAVS in NZAVS wave 10 and NZAVS wave 11, and which may have been lost to follow-up in NZAVS wave 12.
Eligibility criteria
Participants were included in the analysis if they:
Were enrolled in the 2018 wave of the NZAVS (Time 10).
Provided responses to the religious service attendance question at both Time 10 (baseline) and Time 11 (treatment wave).
Participants with missing covariate data at baseline were included, with missing data imputed using information available at baseline. Participants may have been lost to follow-up by the end of the study (Time 12); we adjusted for attrition and non-response using censoring weights, as described below.
A total of 33,198 individuals met these criteria and were included in the study.
Missing data
We adopted the following strategies for handling missing data:
Baseline missingness
We used the predictive mean matching algorithm from the mice package in R (Van Buuren, 2018) to impute missing baseline data (<2% of covariate values). We performed single imputation, using only baseline data for imputation (Zhang et al., 2023).
Outcome missingness
To account for confounding and selection bias from missing responses and panel attrition, we applied censoring weights obtained using nonparametric machine-learning ensembles via the ‘lmtp’ package in R (Williams & Díaz, 2021).
Confounding control
To address confounding, we employ a modified disjunctive cause criterion (VanderWeele, 2019), which involves:
Identifying all common causes of both the treatment and outcomes.
Excluding instrumental variables that affect the exposure but not the outcome.
Including proxies for unmeasured confounders affecting both exposure and outcome.
Controlling for baseline exposure and baseline outcome, serving as proxies for unmeasured common causes (VanderWeele et al., 2020).
The covariates included for confounding control are detailed in Supplement: Appendix B. These methods adhere to the guidelines provided by Bulbulia (2024e) and were pre-specified in our study protocol https://osf.io/ce4t9/.
Figure 1 presents a causal diagram of our identification strategy. The graphs are labelled

Ten causal-directed acyclic graphs clarify distinct threats to valid causal inference in a three-wave panel study.
Figure 1
Figure 1
In contrast, Figure 1
Figure 1
Figure 1
Figure 1
Statistical estimation
We used targeted minimum loss-based estimation (TMLE) to estimate causal effects from the observed panel responses (Van der Laan, 2014; Van der Laan & Gruber, 2012). Specifically, we used an implementation of TMLE that employs machine-learning algorithms to produce efficient statistical estimates for features of the data relevant to estimating our causal estimands from data. These algorithms do not make assumptions about the distributions of the data and are able to learn efficiently even in the presence of high-dimensional data such as ours.
We conducted our estimations using the ‘lmtp’ package in R (Williams & Díaz, 2021), utilising the SuperLearner library with predefined algorithms: ‘SL.ranger’, ‘SL.glmnet’, and ‘SL.xgboost’ (Chen et al., 2023; Polley et al., 2023; Wright & Ziegler, 2017). To ensure robust model performance, we employed 10-fold cross-validation, which guarantees that the data used for training the models are distinct from those used for testing. We produced graphs, tables, and output reports using the ‘margot’ package (Bulbulia, 2024a).
Sensitivity analysis using the E-value
To assess the robustness of our results to unmeasured confounding, we report VanderWeele and Ding’s ‘E-value’ in all analyses (VanderWeele & Ding, 2017). The E-value quantifies the minimum strength of association, on the risk ratio scale, that an unmeasured confounder would need to have with both the exposure and the outcome – after considering all measured covariates – to explain away the observed exposure–outcome association (Linden et al., 2020; VanderWeele et al., 2020). We used the bound of the E-value 95% confidence interval closest to 1 to evaluate the strength of evidence. We note that the E-value presents an approximation of confounding risk.
Scope of interventions
To illustrate the magnitude of the interventions, we present histograms (Figure 2) showing the distribution of religious service attendance during the treatment wave. In Figure 2(a), the intervention for regular religious service affects a more significant portion of the sample than the zero religious service intervention depicted in Figure 2(b). The ‘Regular vs Zero’ comparison addresses the question: what is the difference in effects between a society where religious service is universal compared to one in which religious service is completely absent?

This figure shows a histogram of responses to religious service frequency in the baseline + 1 (i.e. the treatment) wave. Responses above eight were assigned to eight, and values were rounded to the nearest whole number. The red dashed line shows the population average. (a) Responses in the gold bars are shifted to four on the regular religious service intervention. All those responses in grey (four and above) remain unchanged. (b) On the zero-intervention, responses in the blue bars denote those shifted under the zero-intervention treatment.
Changes in religious service attendance
Table 1 shows the transitions in religious service attendance from baseline to the treatment wave. Assessing changes in the treatment variable is essential for evaluating the positivity assumption of causal inference, which, as stated above, requires that every individual has a non-zero probability of receiving each treatment level, regardless of their covariate values (Danaei et al., 2012; Hernan & Robins, 2024; VanderWeele et al., 2020). Although causal inference does not rely on the assumption that every combination of covariate and exposure is realised, where such observations are absent, stronger modelling assumptions are required. We observe that attendance levels at 0 (no attendance) and 4 (weekly attendance) were the most common responses and, furthermore, that there was indeed a change in attendance over time. That is, at the population level, we find changes in reported attendance between the baseline wave and the exposure wave occurred in the data.
This transition matrix captures stability and change in religious service between the baseline and treatment wave.
Each cell in the matrix represents the count of individuals transitioning from one state to another. The rows correspond to the state at baseline (from), and the columns correspond to the state at the treatment wave (to).
Results
Study 1: causal effects of regular church attendance on self-reported volunteering and self-reported volunteering and donations
Regular religious service versus zero treatment contrast for donations and volunteering
Results for the treatment contrasts between Regular Religious Service and Zero Religious Service, focusing on self-reported volunteering and charitable donations, are displayed in Figure 3(a) and Table 2. These results are estimated on the causal difference scale.

This figure graphs the results of model estimates for the three causal contrasts of interest on reported charitable behaviours at the study’s end. The causal contrasts are: (a) Regular versus Zero Religious Service, (b) Regular Religious Service versus Status Quo, and (c) Zero Religious Service versus Status Quo. Contrasts are expressed in standard deviation units.
This table reports the results of model estimates for the causal effects of a universal gain of weekly religious service vs a universal loss of weekly religious service on reported charitable behaviours at the end of the study. Contrasts are expressed in standard deviation units.
For donations, the effect estimate is 0.132 [0.102, 0.161]. The E-value for this estimate is 1.507, with a lower bound of 1.426. At this lower bound, unmeasured confounders would need a minimum association strength with both the intervention sequence and outcome of 1.426 to negate the observed effect. Weaker associations would not overturn it. We infer evidence for causality. On the data scale, this intervention represents a difference of NZD 656.58 per adult per year in charitable giving compared with the zero attendance intervention.
The effect estimate for hours volunteered is 0.123 [0.09, 0.156]. The E-value for this estimate is 1.482, with a lower bound of 1.389. We infer evidence for causality. On the data scale, this intervention represents a difference of NZD 30.21 minutes per adult per week in volunteering compared with the zero attendance intervention.
Regular religious service versus status quo treatment contrast for donations and volunteering
Figure 3(b) and Table 3 present results for the treatment contrasts between regular religious service and status quo, focusing on self-reported volunteering and charitable donations. These results are estimated on the difference scale.
This table reports results of model estimates for the causal effects of a universal gain of weekly religious service attendance vs the status quo on reported charitable behaviours at the end of the study. Contrasts are expressed in standard deviation units.
For donations, the effect estimate is 0.121 [0.102, 0.14]. The E-value for this estimate is 1.477, with a lower bound of 1.422. At this lower bound, unmeasured confounders would need a minimum association strength with both the intervention sequence and outcome of 1.422 to negate the observed effect. Weaker confounding would not overturn it. We infer evidence for causality. On the data scale, this intervention represents an increase of NZD 601.87 per adult per year in expected charitable giving over the status quo.
For hours volunteered, the effect estimate is 0.095 [0.066, 0.123]. The E-value for this estimate is 1.404, with a lower bound of 1.317. We infer evidence for causality. On the data scale, this intervention represents an increase of 23.33 minutes per adult per week in hours volunteering over the status quo.
Zero religious service versus the status quo treatment contrast for donations and volunteering
Figure 3(c) and Table 4 present results for the treatment contrasts between zero religious service and status quo, focusing on self-reported volunteering and charitable donations. These results are estimated on the causal difference scale.
This table reports the results of model estimates for the causal effects of a universal loss of weekly religious service attendance vs the status quo on reported charitable behaviours at the end of the study. Contrasts are expressed in standard deviation units.
For donations, the effect estimate is −0.011 [−0.029, 0.008]. The E-value for this estimate is 1.111, with a lower bound of 1. We infer that the evidence for causality is not reliable. On the data scale, this intervention represents a difference of NZD −54.72 per adult per year in charitable giving compared to the status quo. Still, again, this effect is not reliable.
The effect estimate for hours volunteered is −0.028 [−0.042, −0.014]. The E-value for this estimate is 1.189, with a lower bound of 1.128. We infer evidence for causality. On the data scale, this intervention represents a difference of −6.88 in volunteering minutes compared with the status quo.
Study 2: Causal effects of regular Church attendance on support received from others – time
Regular vs zero causal treatment contrast for time received from others
Figure 4(a) and Table 5 present results for the treatment contrasts between regular religious attendance and zero, focusing on voluntary help received from others during the past week (yes/no). These results are estimated on the risk ratio scale.

This figure reports the results of model estimates for the three causal contrasts of interest on help received from others during the past week (yes/no). The causal contrasts are (a) Regular vs Zero Religious Service Attendance, (b) Regular Religious Service Attendance vs Status Quo, and (c) Zero Religious Service Attendance vs Status Quo. Contrasts are expressed on the risk ratio scale.
This table reports the results of model estimates for the causal effects of a universal gain of weekly religious attendance vs a universal loss of weekly religious attendance on voluntary help received from others during the past week (yes/no) at the end of the study. Contrasts are expressed on the risk ratio scale.
For community gives time, the effect estimate is 1.378 [1.231, 1.541]. The E-value for this estimate is 2.1, with a lower bound of 1.764. At this lower bound, unmeasured confounders would need a minimum association strength with both the intervention sequence and outcome of 1.764 to negate the observed effect. Weaker confounding would not overturn it. We infer evidence for causality.
For friends give time, the effect estimate is 1.187 [1.108, 1.271]. The E-value for this estimate is 1.658, with a lower bound of 1.454. We infer evidence for causality.
For family gives time, the effect estimate is 0.95 [0.901, 1.003]. The E-value for this estimate is 1.288, with a lower bound of 1. We infer that evidence for causality is not reliable.
Regular religious service vs status quo treatment contrast for time received from others
Figure 4(b) and Table 6 present results for the treatment contrasts between regular religious service and status quo, focusing on voluntary help received from others during the past week (yes/no). These results are again estimated on the risk ratio scale.
This table reports the results of model estimates for the causal effects of a universal gain of weekly religious attendance vs the status quo on voluntary help received from others during the past week (yes/no) at the end of the study. Contrasts are expressed on the risk ratio scale.
For community gives time, the effect estimate is 1.289 [1.174, 1.415]. The E-value for this estimate is 1.899, with a lower bound of 1.626. At this lower bound, unmeasured confounders would need a minimum association strength with both the intervention sequence and outcome of 1.626 to negate the observed effect. Weaker confounding would not overturn it. We infer evidence for causality.
For friends give time, the effect estimate is 1.128 [1.061, 1.199]. The E-value for this estimate is 1.508, with a lower bound of 1.315. We infer evidence for causality.
For family gives time, the effect estimate is 0.958 [0.913, 1.006]. The E-value for this estimate is 1.258, with a lower bound of 1. We infer that evidence for causality is not reliable.
Zero religious service vs status quo treatment contrast for time received from others
Figure 4(c) and Table 7 present results for the treatment contrasts between zero religious service and status quo, focusing on voluntary help received from others during the past week (yes/no). These causal effect estimates are again expressed on the risk ratio scale.
This table reports results of model estimates for the causal effects of a universal loss of weekly religious service attendance vs the status quo on voluntary help received from others during the past week (yes/no) at the end of the study. Contrasts are expressed on the risk ratio scale.
For family gives time, the effect estimate is 1.008 [0.991, 1.026]. The E-value for this estimate is 1.098, with a lower bound of 1. We infer that evidence for causality is not reliable.
For friends give time, the effect estimate is 0.95 [0.928, 0.973]. The E-value for this estimate is 1.288, with a lower bound of 1.197. At this lower bound, unmeasured confounders would need a minimum association strength with both the intervention sequence and outcome of 1.197 to negate the observed effect. We infer evidence for causality.
For community gives time, the effect estimate is 0.936 [0.889, 0.985]. The E-value for this estimate is 1.339, with a lower bound of 1.14. We infer evidence for causality.
Study 3: causal effects of regular Church attendance on support received from others – money
Regular vs zero causal contrast on money received from others
Figure 5(a) and Table 8 present results for the treatment contrasts between regular religious service and zero, focusing on money received from others during the past week (yes/no). These results are again presented on the risk ratio scale.

This figure reports the results of model estimates for the three causal contrasts of interest on help received from others during the past week (yes/no). The causal contrasts are: (a) Regular vs Zero Religious Service Attendance (b) Regular Religious Service Attendance vs Status Quo; (b) Zero Religious Service Attendance vs Status Quo. Contrasts are expressed on the risk ratio scale.
This table reports the results of model estimates for the causal effects of a universal gain of weekly religious service vs a universal loss of weekly religious service on financial help received from others during the past week (yes/no) at the end of the study. Contrasts are expressed on the risk ratio scale.
For community gives money, the effect estimate is 1.376 [1.112, 1.703]. The E-value for this estimate is 2.095, with a lower bound of 1.465. At this lower bound, unmeasured confounders would need a minimum association strength with both the intervention sequence and outcome of 1.465 to negate the observed effect. We infer evidence for causality.
For family gives money, the effect estimate is 1.137 [1.028, 1.258]. The E-value for this estimate is 1.532, with a lower bound of 1.198. We infer evidence for causality.
For friends gives money, the effect estimate is 1.137 [0.964, 1.342]. The E-value for this estimate is 1.532, with a lower bound of 1. We infer that evidence for causality is not reliable.
Regular vs status quo causal contrast on money received from others
Figure 5(b) and Table 9 present results for the treatment contrasts between regular religious service and status quo, focusing on money received from others during the past week (yes/no). These results are expressed on the risk ratio scale.
This table reports the results of model estimates for the causal effects of a universal gain of weekly religious service vs the status quo on financial help received from others during the past week (yes/no) at the end of the study. Contrasts are expressed on the risk ratio scale.
For community gives money, the effect estimate is 1.254 [1.098, 1.432]. The E-value for this estimate is 1.818, with a lower bound of 1.426. At this lower bound, unmeasured confounders would need a minimum association strength with both the intervention sequence and outcome of 1.426 to negate the observed effect. Weaker confounding would not overturn it. We infer evidence for causality.
For family gives money, the effect estimate is 1.13 [1.037, 1.232]. The E-value for this estimate is 1.513, with a lower bound of 1.233. At this lower bound, unmeasured confounders would need a minimum association strength with both the intervention sequence and outcome of 1.233 to negate the observed effect. We infer evidence for causality.
For friends gives money, the effect estimate is 1.041 [0.951, 1.139]. The E-value for this estimate is 1.248, with a lower bound of 1. We infer that evidence for causality is not reliable.
Zero vs status quo causal contrast on money received from others
Figure 5(c) and Table 10 present results for the treatment contrasts between zero religious service and status quo, focusing on money received from others during the past week (yes/no). These results are expressed on the risk ratio scale.
Table reports results of model estimates for the causal effects of a universal loss of weekly religious service vs the status quo on financial help received from others during the past week (yes/no) at the end of study. Contrasts are expressed on the risk ratio scale.
For family gives money, the effect estimate on the risk ratio scale is 0.993 [0.953, 1.035]. The E-value for this estimate is 1.091, with a lower bound of 1. We infer that evidence for causality is not reliable.
For friends gives money, the effect estimate on the risk ratio scale is 0.915 [0.809, 1.036]. The E-value for this estimate is 1.412, with a lower bound of 1. We infer that evidence for causality is not reliable.
For community gives money, the effect estimate on the risk ratio scale is 0.911 [0.796, 1.042]. The E-value for this estimate is 1.425, with a lower bound of 1. We infer that evidence for causality is not reliable.
Additional study: comparison of causal inference results with cross-sectional regressions
To clarify how our causal inferences compare with estimates from commonly used observational research methods, we quantified the statistical associations between religious service attendance and all prosocial outcomes using cross-sectional methods frequently employed in observational psychology. For each analysis, we included all regression covariates from the causal models, including sample weights, while omitting the baseline measurement of the outcome variable.
Cross-sectional volunteering result
The change in expected hours of volunteer work for a one-unit increase in religious service attendance is b = 0.31; (95% confidence interval (CI): 0.28, 0.34). Multiplying this by 4.2 gives a monthly estimate of 77.95 minutes. This result is 2.58% greater than the effect estimated from the ‘regular vs zero’ causal contrast, revealing an overstatement in the cross-sectional regression model.
Cross-sectional charitable donations result
The coefficient for religious service on annual charitable donations suggests a change in expected donation amount per unit increase in attendance is b = 451; (95% CI: 408, 494). When adjusted to a monthly rate by multiplying by 4.2, this value equals NZ Dollars 1894.45. It is 2.89% greater than our causal contrast estimate, again revealing an overstatement in the cross-sectional regression model.
For Studies 2 and 3, which focus on community help received, we handled the non-collapsibility of odds ratios by assuming a Poisson distribution for the outcome variables, and thus obtaining a rate ratio that approximates a risk ratio (Huitfeldt et al., 2019; VanderWeele et al., 2020).
Cross-sectional community assistance received result: Time
The exponentiated change in expectation for a one-unit change in religious service attendance is b = 1.17; (95% CI: 1.14, 1.19) approximate risk rate ratio. The monthly rate ratio derived by multiplying this coefficient by 4.2 is 1.921. This estimate is 1.39% greater than the ‘regular vs zero’ causal estimate, revealing an overstatement of effect-estimate from that the cross-sectional regression model suggests.
Cross-sectional community assistance received result: Money
Similarly, the exponentiated change for money received yields an approximate risk ratio of b = 1.18; (95% CI: 1.08, 1.27). The monthly risk ratio, after adjustment, is 1.996. This rate ratio is 1.45% greater than the causal estimate, again revealing an overstatement in the cross-sectional regression model.
These findings indicate that while cross-sectional regression results may be suggestive, they can differ substantially from those obtained through the causal analysis of panel data.
Economic effects of religious service attendance on charitable donations
We next use our results to estimate the approximate economic value of religious service attendance under different scenarios—‘Regular Religious Service Attendance’, ‘Zero Religious Service Attendance’, and the ‘status quo’, focusing on charitable donations. We find the expected total donations under the interventions we modelled are as follows:
Regular Religious Service: an increase in religious service attendance yields an individual average donation sum of
Zero Religious Service: Reducing religious service attendance to zero yields an average donation sum of
Status Quo: The expected individual average donation sum is currently
With 3,989,000 adult residents in New Zealand in 2021. 1 We find that:
First, multiplying the adult population by the average donation sum gives a status quo national estimate for charitable giving of
Second, the net gain to charity from country-wide regular attendance at religious services, compared to the status quo, is
Third, although the net cost to charity from a complete cessation of regular religious service attendance is
To provide context, we next consider the magnitude of these economic consequences by comparing them to the New Zealand government’s annual budget. In the year the outcomes were measured (2021–2022), this budget amounted to NZD 57,976,000,000.
First, the expected gain from a nationwide adoption of regular religious service represents 4.1% of New Zealand’s annual government budget 2021.
Second, we do not find a reliable one-year effect on charitable giving from the loss intervention.
Hence, the scenario in which all New Zealand adults were to regularly attend religious services implies a substantial increase in society-wide charitable support one year after the intervention compared to the status quo. On the contrary, the scenario in which New Zealand was to experience a complete cessation of religious service attendance is not reliably distinguishable from the status quo condition.
We emphasise that these population-wide estimates for aggregated one-year effects following interventions reflect short-term behavioural changes. They do not account for the longer-term implications of gaining or losing religious institutions on charity and volunteering.
Discussion
Our study provides causal evidence that adopting regular religious service attendance across New Zealand would increase levels of charity and volunteering. On the other hand, results also suggest that completely eliminating religious services would produce relatively minor changes in the year following cessation – a finding attributable to low baseline levels of religious attendance.
This study underscores the importance of careful causal inquiry for examining the social consequences of religion. The broad question, ‘Does religion cause prosociality?’ is too vague to yield meaningful insights. Effective investigation requires specifying causal contrasts with well-defined treatments, selecting relevant measures of religious belief and behaviour, defining the target population, and collecting appropriate repeated-measures data in sufficiently large samples over time. Furthermore, only after satisfying fundamental causal assumptions and identification criteria can we derive reliable statistical estimates that address our research questions. Employing flexible statistical estimators further mitigates the risk of model misspecification (Hoffman et al., 2023; Van der Laan & Rose, 2018; Wager & Athey, 2018). Following statistical estimation, sensitivity analyses are essential to assess the robustness of our quantitative estimates to confounding bias (Hernan & Robins, 2024; VanderWeele et al., 2020).
In this study, we defined clear causal contrasts based on interventions that either increased or decreased religious service attendance, or maintained attendance at the status quo. Our statistical models accounted for baseline confounders, including measures of baseline treatment and baseline outcomes. To minimise reliance on modelling assumptions, we employed flexible, doubly robust machine learning ensembles with cross-validation. We then compared expected average outcomes under different treatments one year after the interventions, ensuring that estimation adhered to a temporal order in which causes precede effects.
We introduced two novel measures of prosociality – help received and financial support from one’s community. Although these indicators may be subject to measurement error, they help mitigate concerns related to self-presentation bias, as few would boast of dependency. Findings from these indicators cooberate the findings we obtain from self-reported charitable donations and volunteering, lending support to evolutionary theories of religious prosociality, which propose that a fundamental, evolved function of religion is to enhance community-making (Sosis & Bressler, 2003; Watts et al., 2015, 2018; Whitehouse et al., 2023).
Importantly, our findings suggest that traditional effect size measures, such as Cohen’s d or
An intriguing question remains: who benefits from the charity of religious individuals? Although our data indicate that religious individuals often receive help from their communities, it is unclear whether those providing the help are members of the same religious group. For instance, religious service attendance might increase requests for assistance – begging – which non-religious community members could then fulfill. Similarly, religious giving and volunteering might disproportionately benefit religious elites, potentially at the expense of broader public goods. Although these scenarios may seem extreme, our data do not exclude them because our data do not describe the network structure of charitable giving. Nevertheless, we may look to other data sources for insights.
In New Zealand, religious organisations operate as public charities with transparent financial records. Public data show that religious institutions account for 40% of the charitable sector in New Zealand (McLeod, 2020 p. 17), with similar or even higher proportions observed in other countries (Brooks, 2004; Monsma, 2007; Woodyard & Grable, 2014). Furthermore, evidence suggests that religious organisations are efficient charities with low administrative costs and high volunteer engagement (Bekkers & Wiepking, 2011b; Khanna et al., 1995; McLeod, 2020 p. 26). Of course, we cannot dismiss the possibility that some religious leaders pocket a large share of benefits from religious charities. However, given that all charities in New Zealand require open books, with published institutional salaries and legal punishments for theft, such instances are presumably rare. Secular charities, or course, are not immune to corruption. Measurement error poses a threat. However, whether this threat is upwardly biased to favour religion-based charity and volunteering remains uncertain. Based on contextual information, it seems credible to speculate that most religious charity goes to genuinely charitable causes. Nevertheless, we reiterate that our data do not address this question.
Despite the strengths of our study, several limitations persist. It is important not to confuse the precision of our causal workflow with the precision of the estimates that we derive from it. Direct and correlated measurement errors can skew findings by inflating or reducing estimated magnitudes of true effects (Bulbulia, 2024d; VanderWeele & Hernán, 2012). As noted by Bekkers and Wiepking (2011a), even uncorrelated errors can distort estimates of charitable giving, leading to downward biases in religious charity estimates. Although we employed multiple measures and adjusted for baseline giving to reduce the effects of systematic errors, measurement errors may nevertheless influence our results. For instance, our models might underestimate the true costs of religious decline in the near term. Since religious institutions comprise 40% of New Zealand’s charitable sector, a significant drop in religious attendance could threaten the viability of these charities. Put differently, our one-year estimates of charitable giving and volunteering may present a false picture of security. Conversely, a sharp rise in religious attendance might yield more considerable benefits than our models suggest. We do not claim to possess a magical crystal ball; our results offer signals, not certainties.
The generalisability of our findings beyond New Zealand is yet to be examined. Although our results pertain to the New Zealand population, further research is necessary to determine whether they hold in other cultural or national contexts. Initiatives such as the Global Flourish Study (B. R. Johnson & VanderWeele, 2022), and similar efforts will eventually provide opportunities to apply causal methods across diverse cultural settings. Also remaining to be investigated are the causal effects on prosociality of prolonged exposure to religious services over many years, another fascinating horizon for future investigations.
In conclusion, by integrating robust causal methods with longitudinal data, we provide clear quantitative evidence that religious service attendance influences prosocial behaviours. Beyond the significance of our findings, we hope this study encourages psychological scientists to adopt causal methods in their research. Psychological questions are inherently causal; however, effectively addressing these questions requires careful causal methods that differ from the associational approaches that currently dominate teaching and research (Bulbulia, 2024e; Rohrer et al., 2022; VanderWeele, 2021). We hope this study inspires others to build upon this early effort to ask, and answer, clearly defined causal questions about the causal effects of religious service attendance on charity and volunteering.
Supplemental Material
sj-docx-1-prj-10.1177_00846724241302810 – Supplemental material for The causal effects of religious service attendance on prosocial behaviours in New Zealand: A national longitudinal study
Supplemental material, sj-docx-1-prj-10.1177_00846724241302810 for The causal effects of religious service attendance on prosocial behaviours in New Zealand: A national longitudinal study by Joseph A Bulbulia, Don E Davis, Kenneth G Rice, Chris G Sibley and Geoffrey Troughton in Archive for the Psychology of Religion/Archiv Für Religionpsychologie
Footnotes
Author contribution
J.B. conceived the study and approach and wrote the draft manuscript. C.S. led NZAVS data collection. All authors contributed to the manuscript.
Data availability
The data described in the paper are part of the New Zealand Attitudes and Values Study. Members of the NZAVS management team and research group hold full copies of the NZAVS data. A de-identified dataset containing only the variables analysed in this manuscript is available upon request from the corresponding author or any member of the NZAVS advisory board for replication or checking of any published study using NZAVS data. The code for the analysis can be found at:
.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The New Zealand Attitudes and Values Study is supported by a grant from the Templeton Religious Trust (grant nos. TRT0196 and TRT0418). J.B. received support from the Max Planck Institute for the Science of Human History. The funders had no role in preparing the manuscript or deciding to publish it.
Ethical approval
The University of Auckland Human Participants Ethics Committee reviews the NZAVS every 3 years. Our most recent ethics approval statement is as follows: The New Zealand Attitudes and Values Study was approved by the University of Auckland Human Participants Ethics Committee on 26 May 2021 for 6 years until 26 May 2027 (ref. no. UAHPEC22576).
Notes
References
Supplementary Material
Please find the following supplemental material available below.
For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.
For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.
