A statistical analysis revealing the complex dynamics of offertory collection changes across U.S. Catholic parishes during the pandemic.
Jesus promised that the "gates of hell shall not prevail against" His Church, but students of history know that promise does not extend to all places and all times. The Catholic Church has often been faced with floods, earthquakes, famine, and plagues that have challenged the people of God. Sometimes the Church rises to these challenges — helping the afflicted, increasing in devotion, and sharing their faith. Other times, faith is shaken and souls are lost. Looking back to 2020, the global Catholic Church was faced with a trial of the strength of her institutions, devotion, and charity.
The essence of perfection is to embrace the will of God in all things, prosperous or adverse. ... Our conduct in such instances is the measure of our love for God.
How did the Catholic Church conduct itself during the COVID-19 pandemic? Bishops and pastors made difficult decisions to navigate health recommendations and restrictions while ministering to the spiritual and physical needs of the faithful. Looking back with the benefit of hindsight, this trial can serve as a wide-scale stress test. How did the Catholic Church persevere in adversity when faced with difficult decisions, critical communications, and a desperate need for spiritual care? As St. Alphonsus de Liguori notes, it is easy to unite ourselves to God's will in times of prosperity, but adversity is the measure of our love of God. Taking that principle to a broad application, measures of Church performance during the COVID-19 pandemic are informative of the overall resilience and devotion of the Church in America. Ideally, there would be readily available data on the practice of the Catholic faith in America collected at the parish-level, including estimates of attendance, levels of engagement, and conversions. However, these data were not readily available and only with significant effort was reporting by The Pillar able to extract data from parish bulletins, revealing an average 12% drop in offertory collections across 100 parishes. Considering other research by Hill and Strandberg that indicates that this shift coincided with increased average donation size and larger drops in numbers of donors, this 12% drop may indicate an even larger drop in attendance and engagement.
But descriptive statistics and averages tell only part of the story. Some parishes saw collections decline by 44%, while others actually increased by 35%. What explains this dramatic variation? The Pillar examined the link between demographic and COVID-19 impact metrics in Parish ZIP-codes and counties, but found it difficult to identify many significant explanatory factors. Using the drop in offertory collections as a proxy for parish resilience may offer insight into the critical factors that strengthen the Catholic Church in the face of future challenges. In this research, Acutis Computing used cutting-edge statistical methods to better understand this critical moment in the life of the Church in America.
The main part of the Pillar dataset included weekly offertory collections reported in a selection of parish bulletins, which were processed by hand. Due to significant differences in parish size, collection values were all scaled by the average weekly collection for each parish. As seen in the scaled weekly collection shown in red in the figure below, there was a significant shift in collections starting the second half of March 2020. Before that point, parish collections were about 1.06, or 106% of their average in the overall dataset, but after recovering from the catastrophic drop in March, settled on an average value of 0.9, or 90% of each parish's average.
This pattern appeared very clearly in the overall averages, but it was less evident in each individual parish's reported collections. As seen in this image, parish giving varies week-to-week with local maxima near Christmas and Easter. Averaging limited the size of these variations in the overall weekly average seen in red, but the unexplained variations week-to-week were much more pronounced the raw dataset seen in gray.
A natural first approach might be to examine how each demographic or parish factor relates to collection changes. However, this one-factor-at-a-time approach fails to capture the complex reality. As this example shows, a naïve statistical approach identifies a statistically-significant correlation between reported unemployment in August 2020 and the change in collections. This simple statistical model explains 15.8% of the variation in the observed collections — but is this a real relationship?
Although the positions of the markers and the best fit line (in dark red) suggest a significant relationship between unemployment and the change in collections, a look at the color of the markers show that this relationship is dominated by clusters within the same diocese. For example, the lowest unemployment occurred in the diocese of Fargo, shown in green in the above chart, which also observed increases in parish collections. Furthermore, Chicago, shown in blue, had some of the highest unemployment and these coincided with the largest drops in parish giving. Some clusters, like that of San Antonio, shown in orange, which had large decreases in collections, did not experience unemployment in August to the same degree as New Orleans, shown in pink. This suggests that some other common factor for these dioceses may be introducing a misleading relationship in the data. Strict statewide COVID-19 policy could be both affecting unemployment and collections, generating a spurious correlation. Furthermore, the data included many measures of unemployment with various levels of relationship observed with changes in parish giving. Ultimately, these simple, one-factor-at-a-time approaches often overstate the relationship of individual responses that may have other factors that better contribute to the statistical model as well as obscure relationships that may only exist after controlling for other significant factors.
Observing this diocesan clustering suggested the use of a hierarchical modeling approach to capture the relationship of parishes contained within diocese. This type of statistical model can be complex to implement, without always having a direct mathematical formulation. Therefore, advanced Bayesian methods were used to implement these hierarchical models. Although others, including the Villanova School of Business and Hill and Strandberg have noticed the magnitude of the effect of COVID-19 on parish collections, this is a novel attempt to review this data with hierarchical statistical modeling and an eye towards overall resilience.
Building on The Pillar's work (and generously shared data), presented in three articles: link, link, and link, Acutis Computing designed a sophisticated statistical analysis to improve understanding of the drivers of these changes. Rather than examining simple correlations, this study leveraged complex feature engineering, followed by regularized hierarchical regression, including attempts to extract latent variables and alternative model structures. This approach accounted for the nested structure of parishes within dioceses, as well as the structure of diocese within states and provinces. With the improvements of Large Language Models since this work in 2021, Acutis Computing also used Artificial Intelligence to review the COVID policy responses of each diocese and provide a score for the relative severity and innovation. Further data on diocesan vitality in 2024 were also joined from this article: link.
This analysis used hierarchical mixed effects models in a Markov Chain Monte Carlo solver with horseshoe prior regularization for variable selection, explicitly modeling diocese, state, and province-level clustering while accounting for 28 potential predictors spanning demographics, economics, COVID impacts, and parish vitality metrics. Diagnostics were performed on the resultant statistical model to ensure compliance with the model assumptions. Other approaches, including CatBoost, Random Forests, and other machine learning techniques were also applied to the data. Due to the limited number of observations, these methods did not provide additional insight.
*Diocese vitality data was included from a recent 2025 article. Ideally, 2019 data would have been used to demonstrate causality. Use of 2024 data introduces the possibility of a correlative relationship (i.e. the COVID-19 response caused the diocese vitality measures, rather than vice versa, as presented here). However, U.S. Catholicism has been largely affected by consistent long-term trends, so these observations should be strongly correlated with those from 2019.
†Diocesan policy was analyzed through a series of prompts first using OpenAI ChatGPT 5, followed by a review of those results with Claude Opus 4.1. These reviews were done only by providing a list of diocese names to the Artificial Intelligence -- no information on the other predictors or response was provided in the prompts, which should prevent any hallucinated relationships.
Although these data were collected by individual review of weekly bulletins, the questions are about the parishes and dioceses. For this reason, it is necessary to transform these many observations taken over a two-year time interval into a singular measurement for each parish. A fairly intuitive measurement was used in the work by the Pillar applying this data: percent change from 2019 to 2020 total collections. However, the drop in collections was not evident in the data until March — and so the signal of the before and after effect of COVID could be obscured by behavior in January and February of 2020. The first task was designing a data feature that emphasized this signal while removing this and other sources of noise.
Two example parishes are included in the following section to illustrate the new feature.
Parish #47 had consistent collections before the COVID-19 lockdowns with a very consistent downturn in April, with some recovery throughout the year. Although there was a local maximum in collections in December 2020, collections in December 2019 were similar to other months that year. This parish was in the San Bernardino diocese, California, where the diocese adopted strict COVID-19 precautions. The bishop also encouraged outdoor worship, which signified some innovation in attempting to bring the Blessed Sacrament to the faithful. The median age was 60 years, 23% of the region surrounding the parish reported Hispanic or other heritage, and 804 COVID deaths were recorded for every 100,000 people in the county. Furthermore, as of 2024, they had 1 seminarian for 91,424 Catholics.
This parish had a 43% decrease in giving between 2019 and 2020, which was the largest observed. By taking comparing the difference of the mean of the red points in the blue region (shown as a solid, darker blue line in the right plot) from the mean of the red points in the tan region (shown as a solid, lighter brown line in the right plot), then expressing that difference in the means in terms of the standard deviation of the difference (each mean's standard deviation is shown by the dashed lines of their respective color), a t-statistic of -7.21 was calculated and used as the COVID metric for further prediction. Although this was a very large drop due to the significant decrease in collections and the consistency of the pattern, it was not the largest observed in the dataset — 5 other parishes had a greater COVID impact with the new metric.
Parish #8 was selected to provide some contrast with parish #47. Its week-to-week collections, seen in gray, were considerably more variable. Furthermore, there is a small difference in the after-lockdown period, but the average after-lockdown month is approximately equivalent to a 'bad' before-lockdown month. In this parish, collections in 2020 were elevated in November 2019, December 2019, and January 2020, as well as a sharp recovery in December 2020. As a parish in the diocese of Columbus, Ohio, parish #8 was subject to very typical COVID-19 precautions, with fairly typical innovation to address those limitations as well. The median age was 35 years, only 4% of the surrounding area was Hispanic or other heritage, and 780 COVID-19 deaths were reported for every 100,000 people in the county. Finally, as of 2024, they had 1 seminarian for every 15,317 Catholics.
This parish only had a 4% decrease in giving between 2019 and 2020, which was above the average decrease and less than a tenth of the drop seen in parish #47. Following the same process of collecting the monthly (red) observations in the blue and tan regions, whose respective histograms can be seen in the right part of the figure, and computing the mean difference and its standard error, a COVID metric of -1.43 was computed. Although this value is more abstract than a direct percent decrease (i.e. the mean difference is 1.43 standard deviations below 0), it enabled the development of statistical models that were better able to identify the relevant factors that contributed to changes in collections following the COVID lockdowns. The value observed for this parish was close to the median — 48% of parishes had a COVID impact larger than this one.
Between the original 2019-2020 percent change metric, 79% of parishes experienced a decrease in parish collections. Using the new metric, 81% experienced this decrease. The median percent decrease was 7% and the median COVID metric was -1.40. The standard deviation of the percent change was 14%, but the same for the COVID metric was 2.7.
Other time-series and highly correlated variables were included in the initial dataset, including cumulative COVID-19 deaths, month-by-month unemployment, and demographic indicators of socio-economic status. Each of these variable sets were also treated with detailed feature engineering to remove highly collinear terms from the model, enabling better model convergence and improved interpretability.
Even before model results, the enhanced COVID metric highlights the differences between dioceses and parishes. Each black dot in the following chart represents the impact of COVID lockdowns on a parish's collections, as defined in the section on feature engineering. Note that most of these markers are below zero, indicating a drop in collections following the lockdown. The larger the negative value, the greater the drop in collections. Other markers are close to or above zero, which indicates that there was either no change or an increase in collections during the COVID period. The parishes are ordered so that parishes within the same diocese are adjacent to each other. This makes it clear that parishes from the same diocese are more likely to have similar changes in their COVID impact.
The red horizontal line is the average COVID impact across each diocese in the Pillar dataset. The red shaded area surrounding each horizontal line indicates ±1 standard deviation around the mean within each diocese. Some dioceses, like Sioux Falls, had small variation between parishes and a COVID metric close to 0, while others, most notably Fort Worth, had large variation in parish impacts, with some very large impacts and other noticeable increases in collections.
The two lines in the separate panel on the right side of the plot indicate ±1 standard deviation of the diocese-level data around its mean, and the standard deviation of the within-diocese variation of parishes around the same mean. The larger shaded area for diocese indicates that more of the observed variation in the COVID metric is due to each parish's diocese than other parish-specific impacts. In fact, the Intraclass Correlation Coefficient (ICC) for diocese in these data is ~50%, indicating that about half of the variation observed between parishes is actually variation between dioceses. Based on this observation, the first statistical model was a multiple linear regression model on the diocese-level aggregations. Several predictors, including diocesan policy and reported sacrament ratios, were only observed at the diocesan level and the remaining variables were averaged.
Only four variables of 25 remained as meaningfully related with the COVID metric after regularization. Several other models were fit, including those that included all variables with non-zero effect sizes, others that modeled interactions and latent socio-economic status, etc. However, the model that only included these four terms most accurately represented the data.
In the final model, only Hispanic/Other had a 95% credibility interval that excluded zero, which is indicated by the red bar in the above figure. This suggested that Catholic ethnic sub-groups responded to, or were more vulnerable to, disruptions from the COVID-19 pandemic or lockdowns differently than other Catholic ethnicities. These other variables also had a large statistical effect; however, there is a greater chance that this same procedure conducted on a larger sample of parishes or dioceses could change these observations.
To further contextualize the results, the diocesan predictions were plotted on the same chart as the overall dataset. The red lines indicate the actual mean for the diocese observed in the data and the blue line represents the value that the model would predict based on the log-ratio hispanic/other, infant baptism ratio, etc. for that specific diocese. The blue shaded area indicates the error that the model expected for its prediction. The black arrow designates the residual error of the model — an arrow pointing up indicates that the diocesan average was higher than that anticipated by the model and an arrow pointing down indicates a diocese that had a larger drop than expected by the model. For example, the Diocese of Columbus should have had a drop of -1.7 according to the model, but the actual value was -1.0. Therefore, the Diocese of Columbus outperformed by model by 0.7 and the arrow points up, indicating the amount of overperformance. This falls within the blue shaded credibility interval, indicating that this was an expected variation. The largest residuals show where the model may have had some trouble fitting the data or significant factors were missing from the dataset. The Diocese of San Antonio, Texas, the Diocese of San Bernadino, California, and the Diocese of Yakima, Washington all significantly underperformed the model. On the other hand, the Diocese of New Orleans, Louisiana, the Diocese of Orange, California, and the Diocese of Seattle, Washington had smaller drops in giving than expected by the model.
The plot on the right shows the diocesan mean ±1 standard deviation for all of the parishes (in all of the dioceses), the standard deviation between the dioceses, and the standard deviation of the residual error. The smaller the variation in the residual error, the better the model was able to predict the data. However, because the model was fit on the diocese-level dataset, it performed worse when used to predict the effects on individual parishes.
The collected data was computed from parish bulletins, half of the observed variation occurred at the parish-level, and this analysis set out to identify which parishes demonstrated the most resilience. However, we also identified that comparing parishes without accounting for differences in the diocese could be very misleading and generate spurious conclusions. Therefore, Acutis Computing applied several hierarchical models to properly capture nested relationships in the data: parishes belong to a diocese, dioceses belong to a province, and dioceses are in a U.S. state. Having controlled for the proper nested structure, there can be greater confidence in the results identified by the model. Using cross-validation, a final hierarchical structure that viewed parishes within diocese within states was selected. This selection was interesting, because it suggested that the diocese reaction to state policy and relationship with state policy-makers had a greater effect on resilience than internal church structure.
As before, four variables were identified by the model after regularization. Three were the same as those selected by the diocese model; however, the effect of Hispanic/Other and Infant Baptism Ratio decreased with the inclusion of diocesan structure. This may be due to the demographic information being drawn from the area surrounding each parish, rather than registered parishioners or some other metric of the composition of the parish. Infant Baptism Ratio decreased less than the ethnic group effect, leaving it as the most significant predictor in the dataset. The seminarian ratio effect increased in this statistical model, and a new variable, median age, replaced the 2019 parish collection variable. This suggests that parishes with a strong intergenerational link — indicated by high infant baptisms, vocations to the priesthood, and participation from parents and grandparents — were most resilient to the effects of the pandemic.
Unlike the diocesan model, none of the 95% credibility intervals for these parameters excluded zero. Despite that, many of these effects were large and selected under multiple selection regimes, increasing confidence that they have a meaningful relationship with the COVID resilience metric.
The red dots in the following chart represent the each parish's observed COVID metric and the blue marker indicates the statistical model's estimate for that parish based on its other measures. Recall that smaller numbers, toward the bottom of the chart, indicate less resilience and larger numbers, toward the top of the chart, indicate more resilience. The bars represent the 95% credibility interval for each estimate from the model. The overall diocesan averages are not presented because they are now included in the statistical model. Note that, while many of the estimates were near the observed data, there were also many observations that fell far outside of the estimated error expected by the model. This was more likely to occur in diocese like Fort Worth or Chicago, where the variation between parishes was large. Due to the limited amount of data, the model could not account for the complexity of differences in variance between dioceses, reducing the accuracy of the error bars. Furthermore, while the model was able to capture much of the variation between dioceses, the limitations of the data collection methodology often prevented good predictions for individual parishes. For example, consider the two parishes in Orange County. Both parishes may have been in the same ZIP-code or county and received similar or identical demographic measures and, considering that the seminarian and baptism data was only available at the diocese-level, received the same measures for those as well. This resulted in identical model predictions, despite one parish having a positive COVID metric and the other negative. These data support some insights, but there are considerable limitations when attempting to measure meaningful differences between parishes in the same diocese.
Overall, the hierarchical model explained 53% of the variance between parishes as seen in the right panel of the figure above. This was the same R2 value as if the diocesan average was used as a prediction, suggesting that additional data is needed to fully predict individual parish performance. However, this was also an improvement over the 31% found by the diocesan model when applied to individual parishes.
Parishes in strong dioceses showed the most resilience, especially those with a high median age in dioceses with many infant baptisms and seminarians per capita.
Although socio-economic status and COVID-19 impact metrics did not significantly effect outcomes, minority ethnic groups were more likely to be less resilient. Direct parish vitality was not well-represented in the data, but parishes in vital dioceses were more resilient.
About 50% of the variation was driven by the diocese and mostly the same factors that influenced the diocesan outcomes also influenced parish outcomes. However, limited measures directly linked to individual parishes suggests that opportunities to understand parish-level performance may have been missed or may require further data mining for further insights.
Neither COVID-19 severity in an area nor the severity of the diocesan lockdown had a measurable effect on parish resilience.
1. Diocesan Leadership and Culture Matter
Differences between dioceses accounted for about half of the variation in observed resilience as indicated by the Acutis Computing COVID metric. These differences were more significant than those of the U.S. state and even adjacent dioceses had considerably different mean responses. Even accounting for other demographic and policy-based variables, only half of the diocesan variation was able to be predicted from the model. This suggests that the overall strength of the bishop, chancery, and culture of the diocesan leadership has an undeniable impact on resilience.
Furthermore, ethnic sub-groups within the diocese may be more susceptible to shock or stress. In the past, Catholic ethnic groups that immigrated to the U.S. were subject to increased risk of secularization. This stress test indicates that more recent ethnic groups may be subject to the same secularizing influences and further work is necessary to minister directly to these communities.
2. Vitality Creates Resilience
These results indicate that building resilience is, in many respects, straightforward. A diocese that prioritizes robust sacramental participation and cultivates priestly vocations simultaneously strengthens its ability to withstand future shocks to the faith. By contrast, declines in new seminarians or in the faithful actively living a sacramental life should alert a bishop to potential vulnerabilities. Such weaknesses threaten not only long-term vitality but can yield significant and immediate consequences when the next crisis emerges.
3. Missed Opportunities from Weak Data
As a data company entering into the work of the New Evangelization, it is somewhat shocking to see the limited amount of data available. That a news company found it profitable to go through online bulletins by hand in order to capture a small portion of the finances of the Catholic Church in America, speaks volumes to the limited transparency and technological weakness of Catholic institutions. With the onset of some online platforms, including online giving, communication, and prayer services, some of these attitudes toward technology are shifting. The weakness of overall technical competency, however, is allowing non-Catholic, private aggregation of key information that Catholic leadership and faithful need to promote the health of the Church in America and fuel the New Evangelization. The COVID-19 pandemic provided an opportunity for Church leadership to gain key insights into what's working and what isn't across the United States and the world -- however, from a lack of interest in data and technology, the opportunity was largely missed.
To ensure that future opportunities are not missed, parishes should be surveyed to establish baseline demographics for current parishioners. Registration processes should be modified to include key demographic data, including ethnic groups, age, socio-economic indicators, proximity to the parish, etc. This will enable identification of the current demography of the parish and allow Church leadership to observe changes over time. Furthermore, information about the devotional life of parishes, including belief and practice of its parishioners, as well as community involvement would provide key indications of parish vitality and engagement. Finally, publication of statistics of seminarians and baptisms at the parish-level would immediately improve these statistical models.
Modeling Approach: We employed hierarchical linear mixed effects models with random intercepts for state/diocese, using brms with horseshoe priors for regularized variable selection followed by gaussian priors for final parameter estimation. Direct modeling with STAN was used for latent variable extraction and other prediction methods as needed. All continuous predictors were standardized prior to modeling.
Sample: 99 parishes with complete data from 20 dioceses across 10 ecclesiastical provinces. Missing data (1 parish) handled through complete-case analysis due to small sample size. The data was provided in an anonymized format and the sampling method was unclear to the author.
Model Comparison: We initially tested standard machine learning approaches (Ridge, ElasticNet, tree ensembles, CatBoost) and maximum likelihood based approaches, but found Bayesian hierarchical models substantially outperformed them and better identified misidentified model structures.
Statistical Significance: Given the modest sample size and high dimensionality, many effects showed larger coefficients without reaching significance at a 95% level. In keeping with ASA recommendations, p-values were not highly emphasized in lieu of result contextualization.
Variance Structure: Although there was some exploration of various variance/covariance structures, the limited data size prevented significant exploration. Future work would benefit from a much larger dataset.
Limitations: This analysis uses ZIP-code level demographics as proxies for actual parish demographics, which introduced significant measurement error. Sacramental ratios come from 2024 diocesan reports and may not perfectly reflect 2020 parish vitality. The relatively modest variance explained suggests important unmeasured factors (e.g., pastor communication style, digital infrastructure, pre-existing giving systems) that future research should address. Other research of larger, non-public datasets has indicated that Mass streaming and online giving were highly significant. (Hill and Strandberg)
Code Availability: Analysis conducted in R (brms, loo, rstan packages) and Python (pandas, statsmodels, scikit-learn, plotly). Reproducible code available on GitHub.
Data-driven insights can reveal the factors that drive resilience in your context
Start a Conversation