County-Level Impact of the COVID-19 Pandemic on Excess Mortality Among U.S. Veterans: A Population-Based Study.

BACKGROUND: As the novel coronavirus (COVID-19) continues to impact the world at large, Veterans of the US Armed Forces are experiencing increases in both COVID-19 and non-COVID-19 mortality. Veterans may be more susceptible to the pandemic than the general population due to their higher comorbidity burdens and older age, but no research has examined if trends in excess mortality differ between these groups. Additionally, individual-level data on demographics, comorbidities, and deaths are provided in near-real time for all enrolees of the Veterans Health Administration (VHA). These data provide a unique opportunity to identify excess mortality throughout 2020 at a subnational level, and to validate these estimates against local COVID-19 burden.
METHODS: We queried VHA administrative data on demographics and comorbidities for 11.4 million enrolees during 2016-2020. Pre-pandemic data was used to develop and cross-validate eight mortality prediction models at the county-level including Poisson, Poisson quasi-likelihood, negative binomial, and generalized estimating equations. We then estimated county-level excess Veteran mortality during 2020 and correlated these estimates with local rates of COVID-19 confirmed cases and deaths.
FINDINGS: All models demonstrated excellent agreement between observed and predicted mortality during 2016-2019; a Poisson quasi-likelihood with county fixed effects minimized median squared error with a calibration slope of 1.00. Veterans of the U.S. Armed Forces faced an excess mortality rate of 13% in 2020, which corresponds to 50,299 excess deaths. County-level estimates of excess mortality were correlated with both COVID-19 cases (R2=0.77) and deaths per 1,000 population (R2=0.59).
INTERPRETATION: We developed sub-national estimates of excess mortality associated with the pandemic and shared our data as a resource for researchers and data journalists. Despite Veterans' greater likelihood of risk factors associated with severe COVID-19 illness, their excess mortality rate was slightly lower than the general population. Consistent access to health care and the rapid expansion of VHA telemedicine during the pandemic may explain this divergence. FUNDING: This work was supported by grants from the Department of Veterans Affairs Quality Enhancement Research Initiative [PEC 16-001]. Dr. Griffith's effort was supported in part by the Agency for Healthcare Research & Quality [K12 HS026395].

Evidence before this study

Prior work was identified through a literature search in PubMed and Google Scholar for analyses of excess mortality due to COVID-19 in the United States. The prior literature identified was high-quality and risk of bias was not a concern as these analyses are predictive rather than causal evaluations. Most of the existing evidence on excess mortality in 2020 relies either on aggregated national-level data with limited information on individual demographics or uses site-specific data with low sample sizes.

Added value of this study

Our study builds on this prior work by using individual-level data on over 11 million Veterans to develop sub-national estimates of excess mortality attributable to the pandemic. We build and validate a predictive model that uses a rich set of demographics and comorbidities to estimate projected deaths among Veterans. This provides a more precise estimate of excess mortality than prior work and offers suggestive evidence of the importance of consistent access to health care and the ability to rapidly respond to a pandemic.

Implications of all the available evidence

This study can help researchers further examine excess mortality attributable to the pandemic and, using our approach, help evaluate what interventions may have successfully moderated the effects of the pandemic. Moreover, this work provides a framework to model excess mortality in future waves of the COVID-19 pandemic or pandemics of other novel infectious viruses.

Alt-text: Unlabelled box

Introduction

The novel coronavirus (COVID-19) has directly resulted in millions of deaths worldwide. In response, countries have implemented policies to mitigate disease transmission including travel restrictions, school closures, and mask mandates. , Hospitals and ambulatory care centres postponed most elective and routine procedures, and fear of COVID-19 also led households to voluntarily reduce their mobility and demand for health services.3, 4, 5 In the United States, an estimated 40% of adults avoided care during the COVID-19 pandemic, including 12% who avoided emergency care and 31% who avoided routine care. Early evidence suggests that while some of the forgone care was low value, some was high value. Timely access to care is an important determinant of both short and long-term health, and disruptions in healthcare access may increase the risk of avoidable morbidity and mortality for other conditions. Additionally, in some regions, hospitals’ inpatient and critical care facilities faced overcrowding during the pandemic's peak, which may have led to increased mortality and worsened health outcomes. Nationally, studies estimate that all-cause mortality among the general population increased by approximately 23% from March 2020 through the end of year. However, there is a dearth of research using individual-level data or on subnational excess mortality estimates.

Veterans of the U.S. Armed Forces may be at higher risk than the general population due to their greater comorbidity burdens and older age which are associated with severe COVID-19 illness. Prior research has demonstrated that, similar to the general population, COVID-19 caseloads and mortality disproportionately affects Veterans from marginalized population groups.13, 14, 15 The U.S. Veterans Health Administration (VHA) implemented an early, nationwide response to limit the spread of the pandemic. While VHA inpatient utilization fell by 42% during the pandemic's early months, telehealth encounters grew by 556% from March to April 2020. , VHA electronic health records also become available to researchers in near real-time. As such, the VHA is an ideal system within which to examine the impact of the pandemic on excess all-cause mortality.

This study examined the extent to which Veterans experienced excess mortality attributable to the COVID-19 pandemic at the sub-national level throughout 2020. We validated these estimates against measures of COVID-19 burden and provide our deidentified dataset as a public service to support the research community and data journalists. Our work provides valuable insights into the experience and needs of Veterans and may be used to help VHA leadership and other integrated health systems forecast changes in healthcare access and utilization demand for patients at continued high risk.

Methods

This cross-sectional study was approved by the VA Boston Healthcare System's institutional review board and adheres to the Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) reporting guideline for cross-sectional studies (see Appendix A1 for more information). Informed consent was waived per institutional policy because the research, which includes millions of Veterans, could not be practicably carried out without the waiver or an alteration of the study.

Study Data & Population

We obtained administrative data from the VHA's Corporate Data Warehouse (CDW) for 11.4 million unique Veterans who sought care during 2016-2020. These data include a variety of individual-level characteristics that have previously been associated with mortality risk, including demographics (i.e., age, gender, race, ethnicity), VHA priority group (an eligibility determination which reflects disability related to military service or economic hardships), , and major comorbidities. , Dates of death were obtained from the VHA CDW which assesses death status using notifications from numerous sources including the Social Security Administration's (SSA) Death Master File and government death certificates. County-level measures of COVID-19 deaths and confirmed cases were obtained from the Johns Hopkins Coronavirus Resource Center.

Our unit of analysis was the county-year. We created annual county-level snapshots of the VHA patient for each year 2016-2020. Demographics, socioeconomic characteristics, and comorbidities were identified for Veterans alive as of January 1st of each year. Veterans were assigned to county based on their last reported address of residence, which is available in the VHA CDW. County-level data on comorbidities and demographics were totalled for each county, and proportions were generated as a share of total veterans living in that county-year. Patients’ ages were calculated as of these dates, and the most recent information for other demographics and priority group were included. The VHA classifies each veteran into one of eight priority groups based on military service history, household income, presence of service-connected disabilities, and other factors. Following previous work, we used a one-year lookback period to identify patients’ Quan-Elixhauser comorbidities using International Classification of Disease (ICD) codes.

We used three VHA data holdings to identify deaths among enrolled Veterans; these included the SPatient table in the VHA CDW, the Master Veteran Index, and Vital Status Mini. While reporting delays may be of concern in some cases, these are unlikely to be a major issue in our analysis – the average death was reported within 30 days, and 96% were reported within 90 days. For those who died during a given calendar year, we used a two-year lookback to ensure complete coding of comorbidities. The final analytic aggregated data included 16,070 county-years. We restricted our analysis to counties with five or more Veterans to ensure stability of our estimates, and restricted counties to those with a full five years of data to ensure a balanced panel.

Analytic Approach

Our analysis proceeded in five steps. We first generated cross-validated 5-fold test training pairs using pre-pandemic data (2016-2019). Second, the training sets were each used to estimate eight county-level regression models:

Poisson, pooled

Poisson with county-level fixed effects

Poisson with county-level random effects

Generalized estimating equation (GEE) Poisson

Quasi-likelihood Poisson with county-level fixed effects

Negative binomial, pooled

Negative binomial with county-level fixed effects

Negative binomial with county-level random effects

These models were chosen because they are either indicated for count outcomes, address issues with overdispersion, or have statistical properties making them robust to misspecification of correlation structure. Our outcome was the total number of deaths from any cause during a given year. All models were adjusted for the proportion of enrolees who were in various age groups (40-49, 50-64, 65-79, and 80+); male; married; divorced, separated, or widowed; experiencing houselessness (identified using ICD-10 codes Z59.0 or Z59.1); Hispanic or Latino; in each of priority groups 2 through 8; were identified as having each of the Quan-Elixhauser comorbidities; or self-identified their primary race as one of two groups (Black, other non-white). Models were also adjusted for linear time (in years with 2016 set as 0) and used the natural log of enrolled Veterans as an offset.

Third, we assessed each model's performance using the test dataset. We used the mean regression coefficients from each of the training folds to generate mortality predictions for each test fold and calculated two measures: 1) median squared error (MSE), with error being the difference between predicted and observed county-level mortality, and 2) calibration slope which is the regression slope of the linear predictor for observed and expected mortality. Calibration refers to the agreement between predicted and observed outcomes, with slope values closer to 1 indicating better model agreement.

Fourth, we estimated excess mortality among the Veteran population during the COVID-19 pandemic in 2020. We used the top-performing model (in terms of MSE) to predict county-level mortality in 2020; observed mortality was divided by predicted mortality to generate observed-to-expected (O/E) ratios. We calculated total excess deaths by taking the difference between observed and predicted deaths.

Lastly, we validated our excess mortality estimates against measures of national and local disease burden. We used bivariate regressions to estimate the association between our O/E ratios and either COVID-19 confirmed cases or mortality rates per 1,000 county population. Regressions were weighted by the number of enrolled Veterans to obtain nationally representative estimates. We also stratified the top-performing mortality prediction model by month. To perform monthly analysis, we created monthly snapshots of the VHA patient population for each month over our study period. This step allowed us to explicate the relationship between monthly estimates for excess mortality and COVID-19 deaths at a national level, and to understand how excess mortality among Veterans varied throughout the course of 2020. A priori, we expected the O/E ratios to be greater in areas and months that were most affected by the pandemic. All analyses were performed using Microsoft R Open version 4.0.2 (Redmond, WA).

Role of the Funding Source

The funders had no role in the design and conduct of the study; collection, management, analysis, and interpretation of the data; preparation, review, or approval of the manuscript; and decision to submit the manuscript for publication.

Results

Descriptive Statistics

We observed 426,069 deaths among 9.4 million VHA-enrolled Veterans throughout 2020. There were no substantive differences in the distribution of demographic characteristics or frequency of comorbid conditions in the Veteran population between the 2016-2019 period compared to 2020 (Table 1 ). For rates of specific comorbidities, please see Appendix A2. The enrolee population in 2020 was predominantly male (91.2%), white (65.9%), non-Hispanic (93.9%), married (55.7%), with a mean age of 63.9 years (SD 16.1).

Table 1

Sociodemographic characteristics and comorbidities among Veterans Health Administration enrollees, 2016-2020

	Total Enrollees1		Deceased Veterans2
	2016-2019	2020	2016-2019	2020
Unique Veterans	10,382,735	9,408,093	1,511,101	426,069
Age, mean (std. dev.)	63.5 (16.0)	63.9 (16.1)	79.4 (13.0)	79.8 (13.1)
Age Groups	%	%	%	%
Less than 40	14.7	15.0	1.2	1.3
40-49	10.9	11.0	1.1	1.2
50-64	24.4	24.1	9.1	8.2
65-79	35.8	36.5	29.8	32.3
80 and over	14.2	13.4	58.9	57.0
Male	91.9	91.2	91.1	91.0
Race
White	66.3	65.9	61.1	62.4
Black	14.1	14.4	7.9	9.2
Other	2.7	2.9	1.6	1.7
Not Available	16.8	16.7	29.4	26.6
Hispanic or Latino	5.7	6.1	2.7	3.0
Marital status
Married	56.2	55.7	53.8	53.9
Single/Never Married	14.1	14.8	7.2	7.6
Divorced/Separated/Widowed	25.8	24.9	30.0	29.5
Homeless	1.4	1.4	1.8	2.0
Priority Status3 7 or 8	28.5	27.4	23.6	23.5
Comorbid Conditions4
None	57.9	51.6	52.3	49.1
1 to 2	24.5	26.1	13.3	12.6
3 to 5	14.3	17.6	17.4	18.9
6 or more	3.3	4.7	16.9	19.5

Notes:111,428,134 unique Veteran enrollees during the study period whose comorbidities were identified using data for each calendar year they were enrolled. 2Includes Veterans who died from any cause. Comorbidities of Veterans who died were identified using a 24-month lookback period from time of death and calculated from CDW. 3VHA priority groups are an eligibility determination which reflects disability related to military service, economic hardships, and other factors. 4Indicates the number of Quan-Elixhauser comorbidities observed during the lookback period (see Appendix A1 for a list).

However, there were some notable changes in underlying comorbidities and demographics among Veterans who died. The share of Veterans that died who were aged 65-79 increased by 2.5 percentage points (29.8% vs. 32.3%) in 2020, while the proportion who were Black increased by 1.3 percentage points (7.9% to 9.2%). Approximately 34.4% of Veterans who died during 2016-2019 had 3 or more comorbidities, compared to 38.4% during 2020. Compared to those who died in 2016-2019, those who died in 2020 were more likely to have cardiac arrhythmias, hypertension, diabetes, renal failure, and obesity. These demographics and comorbidities are known to be associated with higher likelihood of severe COVID-19 infection. ,

Mortality Prediction Models

Across our sample, there were more deaths in 2020 than in prior years. The median county experienced 39 deaths pre-2020, which increased to 44 deaths in 2020. Similarly, the average number of deaths in a county rose from 118 deaths pre-2020 and 133 deaths in 2020.

The results of our eight annual models are presented in Table 2 . All models demonstrated excellent agreement between observed and predicted mortality on the test data; calibration slopes ranged from 0.97 to 1.00 and MSE ranged from 23.9 to 25.5. The quasi-likelihood Poisson model and Poisson with county fixed effects outperformed all other models in terms of MSE and had a calibration slope of 1.00. Either model provided estimated an O/E ratio of 1.13 in 2020, suggesting there were 50,299 excess deaths among Veterans in 2020 (95% CI 19,230 to 81,370). The O/E ratios for other models ranged from 1.12 to 1.14. We selected the quasi-likelihood Poisson model as preferred due to its greater flexibility ability to account for overdispersion.

Table 2

Measures of model performance and excess mortality

Model Name	MSE1	Calibration Slope2	O/E Ratio3 (95% CI)
Quasi-Likelihood Poisson, FE	23.9	1.00	1.13 (1.05, 1.24)
Poisson, FE	23.9	1.00	1.13 (1.05, 1.24)
NB, FE	24.3	1.00	1.13 (1.04, 1.23)
NB, Pooled	24.3	0.97	1.13 (1.10, 1.16)
Poisson, GEE	24.5	0.98	1.14 (1.09, 1.19)
Poisson, Pooled	25.1	0.98	1.14 (1.12, 1.15)
Poisson, RE	25.3	1.00	1.12 (1.04, 1.23)
NB, RE	25.5	1.00	1.12 (1.02, 1.24)

Notes: FE: fixed effects; RE: random effects; GEE: generalized estimation equations; NB: negative binomial; MSE: median squared error; O/E: observed-to-expected mortality. 1Compares predicted versus observed mortality within test sets during 5-fold cross validation. 2Coefficients from a bivariate linear regression of observed and predicted mortality within test sets during k-fold cross validation. 3National O/E ratios were calculated by dividing observed mortality versus predicted mortality from any cause during calendar year 2020.

There was substantial geographic variation across counties in terms of observed versus expected mortality for 2020 (Figure 1 ). Counties in the bottom quartile had an average O/E ratio of 0.84, while counties in the upper quartile had an O/E ratio of 1.60. Counties with higher O/E ratios were generally located in Alaska and the Great Plains, South Atlantic, and West South Central regions of the United States.

Figure 1

County-level estimates of excess all-cause mortality during 2020, Source: Authors’ analysis of data from the Veterans Health Administration's Corporate Data Warehouse. Notes: The figure displays observed-to-expected (O/E) mortality ratios for 2020, defined as total observed deaths versus total deaths predicted by the regression model. O/E ratios were estimated via a covariate-adjusted quasi-likelihood Poisson regression with county fixed effects and 5-fold cross validation as described in the text.

Validation with COVID-19 Disease Burden

We observed strong associations between the resulting county-level O/E ratios and measures of both COVID-19 cases (Figure 2 ) and deaths per 1,000 in the general population (Appendix A3). Each one-unit increase in O/E ratio was associated with +49.9 cases (95% CI 48.2 to 50.0, R2=0.77) and +0.86 deaths (95% CI 0.84 to 0.88, R2=0.59) per 1,000 population. Additionally, our monthly estimates for nationwide O/E ratios and excess deaths for 2020 followed similar trends as COVID-19 death rates. Observed mortality during the pre-pandemic months (January and February) was similar to prior years. Observed versus expected mortality peaked in April during the first wave of the pandemic (1.24), fell during June (1.06), but returned to higher levels before peaking once more in December (1.40). (Table 3 ) Taken together, the bivariate correlations and monthly models provide additional validation that our estimates of excess mortality are associated with the pandemic's severity.

Figure 2

County-level burden of confirmed COVID-19 cases in the general population during 2020, Source: Authors’ analysis of data from the Johns Hopkins Coronavirus Resource Center. Notes: The figure displays data on COVID-19 cases for the general population and is not veteran-specific.

Table 3

Observed versus expected mortality and excess deaths in 2020, by month

Month	O/E Ratio1 (95% CI)	COVID-192		Excess Deaths3 (95% CI)
		Cases	Deaths
January	0.97 (0.89, 1.05)	7	0	-1,237 (-4,134, 1,660)
February	1.00 (0.93, 1.09)	17	1	100 (-2,477, 2,677)
March	1.04 (0.96, 1.13)	185,930	4,389	1,221 (-1,524, 3,966)
April	1.24 (1.15, 1.36)	886,086	55,506	7,623 (5,070, 10,176)
May	1.11 (1.03, 1.21)	709,679	45,049	3,493 (973, 6,013)
June	1.06 (0.98, 1.15)	826,336	21,992	1,734 (-671, 4,140)
July	1.16 (1.07, 1.26)	1,898,110	26,158	4,671 (2,218, 7,125)
August	1.14 (1.05, 1.24)	1,449,357	29,163	4,128 (1,655, 6,601)
September	1.13 (1.04, 1.23)	1,195,863	22,946	3,724 (1,315, 6,134)
October	1.11 (1.02, 1.21)	1,898,227	23,548	3,313 (722, 5,904)
November	1.26 (1.16, 1.38)	4,430,369	36,653	8,050 (5,439, 10,660)
December	1.40 (1.29, 1.50)	6,306,460	76,386	13,479 (10,644, 16,314)

Notes:1Observed versus expected (O/E) mortality ratios and excess deaths were estimated via a covariate-adjusted quasi-likelihood Poisson regression with county fixed effects and 5-fold cross validation as described in the text. National O/E ratios were calculated by dividing total observed deaths versus total deaths predicted by the regression model. 2County-level measures of COVID-19 deaths and confirmed cases were obtained from the Johns Hopkins Coronavirus Resource Center. These estimates refer to the general population and are not veteran-specific. 3Excess deaths were calculated as the absolute difference between predicted and actual deaths.

Comparability with National Estimates of Excess Mortality

To aid in comparability of our findings with previously published national estimates excess mortality during the pandemic, we re-estimated the national O/E ratio and excess deaths when restricting our data to the months of March through December. Excess mortality was estimated to be higher during this phase of the pandemic, with an excess mortality rate of 16.7%. This equates to 51,436 excess Veteran deaths out of 359,779 total deaths observed from March through December.

Discussion

While the U.S. Veteran population experienced substantial variation in county-level excess mortality during the COVID-19 pandemic, the VHA's overall excess mortality rate of 16.7% was markedly lower than the 20.8% reported for the general population. A 16.7% excess mortality rate suggests an additional 51,436 Veterans died in 2020 than would have died in previous years, after controlling for demographic characteristics and comorbidities.

The excess mortality rate among Veterans in 2020 is lower than the rate in the general population, even after restricting our analysis to the months after COVID-19 was declared a national emergency in the United States (March-December). This result is perhaps counterintuitive, since Veterans are older and have a greater burden of comorbidities on average compared to non-Veterans. Rapid VHA expansion of pre-existing telehealth infrastructure and the fact that VHA healthcare access is decoupled from employment may explain these differences. COVID-19 disrupted the in-person delivery of many health services, and some private medical professionals and hospital systems struggled to transition to telemedicine during the pandemic's early phases. In contrast, the VHA has used telehealth to deliver behavioural healthcare and improve access for Veterans in rural areas for more than a decade. Lower rates of excess mortality among Veterans may be partially attributable to the pre-existing telehealth infrastructure in the VHA, which facilitated the transition to telemedicine and mitigated the disruption in care delivery relative to other health systems that did not have this infrastructure in place prior to the pandemic. Another possible explanation has to do with the structure of the U.S. health insurance system. The U.S. experiences the largest socioeconomic disparities in healthcare access compared to any wealthy nation, and most Americans rely on receiving health insurance through their employers. Veterans, by contrast, are less likely to lose coverage due to pandemic-related layoffs, compared to the general population.

In this analysis, we developed and validated an approach to estimating excess mortality at the sub-national level. While prior work has estimated excess mortality at a national level using aggregated data, our study is the first to use rich individual-level data to develop sub-national measures of excess mortality for a population that is at heightened risk of severe COVID-19 illness. Access to VHA data better allows researchers to examine the impact of policies on health outcomes in near-real time, and to develop new methods to do so. In addition, we have provided our data as a resource for the community, so that researchers and data journalists may examine the burdens of COVID-19 on their communities and assess interventions that may have moderated the pandemic's effects.

Given the importance of timely access to care for patient's health short- and long-term health status, future research should examine whether continuity in healthcare access or early telehealth adoption moderated the association between local COVID-19 burden and excess mortality.

Limitations

Study limitations include a lack of data on confirmed COVID-19 diagnoses due to the absence of widespread testing, especially in the pandemic's early stages. While the CDW incorporates death records from several federal sources, a small number of Veteran deaths may be excluded due to reporting delays. To estimate the extent of this limitation, we conducted a retrospective analysis and found that >96% of death records are entered into the CDW within 90 days from the date of death. Additionally, because cause-of-death data is imported into the CDW annually and is currently unavailable for 2020, we are not yet able to attribute excess deaths to COVID-19 or other causes. For sparsely-populated rural areas, estimates for changes in mortality or utilization may be more sensitive to annual fluctuations due to the limited number of VHA enrolees. The ecological nature of our study design does not allow us to identify causal mechanisms underlying county-level variations in mortality. Lastly, while our analytic approach may be replicated in a wide variety of settings and populations, our specific results may not be generalizable to non-Veteran populations. While we attempt to control for a variety of comorbidities and other factors, comparisons between VHA and non-VHA data may still reflect underlying population differences. In addition, reporting of deaths may vary between the VHA and non-VHA health systems, which may contribute to the observed differences we report.

Despite these limitations, our analysis and findings are consistent with the growing literature documenting excess mortality during the COVID-19 pandemic. Our finding that rates of excess mortality among Veterans compared favourably to the general population, despite a greater burden of risk factors for severe COVID-19 illness, warrants further investigation. Future research should explore whether specific VHA characteristics such as continuous healthcare access or rapid expansion of telehealth moderated the pandemic's effects on excess Veteran mortality.

Contributors

YF, AL, KS and KG conceived the idea for the study. YF, AL, KS and KG contributed to the study design and analysis plan. YF, AL, KT and KG conducted the statistical analysis. All authors contributed to and approved the final manuscript. KG takes responsibility for the contents of the article.

Data Sharing Statement

Our analytic dataset, statistical code, and county-level estimates of excess mortality are both publicly available as part of a Mendeley Data repository (https://data.mendeley.com/datasets/zmj5jr2r69).

Declaration of Interests

The authors have no conflicts of interest to disclose. The views expressed in this article are those of the authors and do not necessarily reflect the position or policy of the Department of Veterans Affairs, United States Government, Boston University, or Vanderbilt University.