Orphanhood and caregiver death among children in the United States due to all-cause mortality 2000-2021: A Modeling Study

S0 Supplementary Figures and Tables

• Download figure
• Open in new tab

Supplementary Fig. S1: Definitions of orphanhood and grandparent caregiver loss used in this study (orange and purple), with comparison to definitions in previous studies [2021a] (blue).

• Download figure
• Open in new tab

Supplementary Fig. S2: Magnitude of children experiencing the death of a caregiver in the U.S. (A) Estimated number of U.S. children newly experiencing caregiver death as the sum of orphanhood and grandparent caregiver death from all causes (incidence). (B) Estimated number of U.S. children experiencing caregiver death from all causes in their lifetime (prevalence). (C) Incidence rates of caregiver death from all causes among U.S. children. (D) Prevalence rates of caregiver death from all causes among U.S. children.

• Download figure
• Open in new tab

Supplementary Fig. S3: Orphanhood prevalence rates by standardized race & ethnicity in each age group among U.S. children.

S1 Extended methods

This Supplementary Text describes the data sources and statistical methods used to estimate cause-specific incidence and prevalence of orphanhood and grandparent caregiver loss among U.S. residents between 2000 and 2021. Our methods are overall similar to Hillis et al. [2021b,a], and extended here to estimate the magnitude of children aged 0-17 years experiencing the death of one or both parents from all causes or the death of a grandparent caregiver from all causes.

The United Nations International Children’s Emergency Fund (UNICEF) defines orphanhood as children experiencing the death of one or both parents [UNAIDS, 2009]. As in Hillis et al. [2021a], we considered mothers of age 15-66 years and fathers of age 15-94 years so the maximum ages of parents at birth of a child were respectively 49 and 77 years. Mortality data were recorded among U.S. residents, and for this reason, orphanhood estimates are restricted to children of U.S. residents. Grandparents play indispensable roles as caregivers for children [United States Census bureau, 2021, Pew Research Center, 2013, US Census Bureau, 2023b, Sowden et al., 2021, Ellis and Simmons, 2014]. We further considered to provide minimum numbers of children who lost a grandparent caregiver, defined as the death of a grandparent caregiver aged older than 30 years living with children in the same household. Mortality data were recorded among U.S. residents, and so grandparent caregiver death estimates are also restricted to children of U.S. resident grandparent caregivers. In Hillis et al. [2021a], we previously stratified grandparent caregivers based on the care they provided into primary grandparent caregivers who were either the guardians of children or co-resided and provided for basic needs of children, as well as secondary grandparent caregivers who co-resided and provided for some needs of children, following the characterizations of Ellis and Simmons [2014]. Here, we report throughout the sum of primary and secondary grandparent caregiver loss (Supplementary Fig. S1).

S1.1 National-level NCHS mortality data by rankable causes of death, 1983-2021

We obtained line-list mortality data on U.S. residents from the NCHS Vital Statistics portal for each year from 1983 to 2021 [National Center for Health Statistics, 2023]. Data were collected from 1983 onwards because the corresponding children who lost a caregiver in 1983 at age 0 were of age 17 in 2000, and so entered our estimation of orphanhood prevalence in 2000. For each mortality record we retained year of death, the corresponding codes of the underlying cause-of-death (ICD-9 code and 282 cause recode before 1999; ICD-10 code and 113 cause recode after 1999), and demographic data of the decedent including sex, age at death, and information on race and Hispanic origin [National Center for Health Statistics, 2022].

Information about the race and Hispanic origin of decedents in death certificates is typically reported by the surviving next of kin, or on the basis of observation in the absence of an informant [Rosenberg, 1999]. Race and Hispanic origin were reported in different formats across the study period [Ingram et al., 2003], which we mapped to standardized race and Hispanic origin categories as described in Supplementary Table S7. Specifically, we grouped individuals of Hispanic origin and all individuals of non-Hispanic origin by their race, i.e. ‘Hispanic’, ‘Non-Hispanic American Indian or Alaska Native’, ‘Non-Hispanic Asian or Pacific Islander’, ‘Non-Hispanic Black’, ‘Non-Hispanic White’, and we refer to the resulting categories as “standardized race & ethnicity” for simplicity. Data on Hispanic origin were not available for 1983 (see below). Individuals of other races in 1984-1991 and individuals of more than one race in 2021 were not coded consistently and not included in this study. We note this approach to harmonising race reporting over the study period did not necessarily make the primary data fully comparable. Earlier research shows inaccuracies are limited [Heron, 2021], which indicates that the incremental implementation of multiple race reporting in the U.S. is unlikely to have introduced notable bias in our orphanhood estimates.

The underlying cause-of-death is defined by the World Health Organization (WHO) as ‘the disease or injury which initiated the train of events leading directly to death, or the circumstances of the accident or violence which produced the fatal injury’ [World Health Organization, 2016]. For 1983-1998, underlying causes-of-death were in the data classified with the Ninth Revision of the International Classification of Diseases (ICD-9) [World Health Organization, 1977], and grouped further into 282 selected causes of death, termed ‘282 cause recode’ [National Center for Health Statistics, 1979]. For 1999-2021, underlying causes-of-death were in the data classified with the Tenth Revision of the International Classification of Diseases (ICD-10), and grouped further into ‘113 Selected Causes of Death’ [National Center for Health Statistics, 2020].

Due to our interest in orphanhood and caregiver loss, we defined 53 non-overlapping rankable underlying causes-of-death that we termed ‘caregiver loss causes-of-death’ and that apart from small modifications are identical to the 52 rankable underlying causes-of-death from the NCHS 113 Selected Causes of Death list [National Center for Health Statistics, 2020, Anderson et al., 2001]. The modifications are listed in Supplementary Table S8. Specifically, we re-categorized drug-induced causes-of-death into ‘Drug overdose’, joining the ICD-9/ICD-10 causes-of-death ‘Intentional self-poisoning by and exposure to nonopioid analgesics, antipyretics and antirheumatics’, ‘Intentional self-poisoning by and exposure to antiepileptic, sedative-hypnotic, antiparkinsonism and psychotropic drugs, not elsewhere classified’, ‘Intentional self-poisoning by and exposure to narcotics and psychodysleptics [hallucinogens], not elsewhere classified’, ‘Intentional self-poisoning by and exposure to other drugs acting on the autonomic nervous system’, ‘Intentional self-poisoning by and exposure to other and unspecified drugs, medicaments and biological substances’ (E950.0-E950.5; X60-X64); ‘Assault by drugs, medicaments and biological substances’ (E962.0; X85); ‘Accidental poisoning by and exposure to nonopioid analgesics, antipyretics and antirheumatics’, ‘Accidental poisoning by and exposure to antiepileptic, sedative-hypnotic, antiparkinsonism and psychotropic drugs, not elsewhere classified’, ‘Accidental poisoning by and exposure to narcotics and psychodysleptics [hallucinogens], not elsewhere classified’, ‘Accidental poisoning by and exposure to other drugs acting on the autonomic nervous system’, ‘Accidental poisoning by and exposure to other and unspecified drugs, medicaments and biological substances’ (E850-E858; X40-X44); and ‘Poisoning by and exposure to nonopioid analgesics, antipyretics and antirheumatics, undetermined intent’, ‘Poisoning by and exposure to antiepileptic, sedative-hypnotic, antiparkinsonism and psychotropic drugs, not elsewhere classified, undetermined intent’, ‘Poisoning by and exposure to narcotics and psychodysleptics [hallucinogens], not elsewhere classified, undetermined intent’, ‘Poisoning by and exposure to other drugs acting on the autonomic nervous system, undetermined intent’, ‘Poisoning by and exposure to other and unspecified drugs, medicaments and biological substances, undetermined intent’ (E980.0-E980.5; Y10-Y14).

We retained these four drug overdose causes-of-death sub-categories as separate drug overdose sub-groups for the purpose of data harmonization (see below). Correspondingly, we removed the related three drug overdose causes-of-death sub-categories respectively from ‘Intentional self-harm’, ‘Assault’ and ‘Accidents’. We renamed these resulting three categories as ‘Suicide excluding drug overdose’, ‘Homicide excluding drug overdose’ and ‘Unintentional injuries excluding drug overdose’ (see Supplementary Table S8). Finally, to simplify presentation as in Figure 1 in the main text, we primarily considered COVID-19 and the leading five causes-of-death between 1983 and 2021 due to our interest in all caregiver loss. Additionally, we included any cause in the top five for adult ages 15-44 years, as these ages specifically include younger adults who are more likely to be parents of younger children. The only one addition to our causes when taking this approach was ‘Homicide excluding drug overdose’. Supplementary Table S9 summarises our aggregation of the 53 rankable causes-of-death into these seven leading caregiver loss cause-of-death groups and ‘Other’ caregiver loss causes-of-death. The ‘Other’ caregiver loss causes-of-death comprise the remaining 46 rankable caregiver loss causes-of-death and any other causes-of-death that are not included in the NCHS 52 rankable causes-of-death.

We then mapped and aggregated line-list mortality records to the 53 rankable caregiver loss causes-of-death. This was done by mapping the 1983-1998 line list data using the 282 recodes to the 52 NCHS rankable causes-of-death as described in [National Center for Health Statistics, 1988] on page 282ff and page 290ff. For drug-induced causes, we used Table 2 in [National Center for Injury Prevention and Control (U.S.). Division of Unintentional Injury Prevention. Health Systems and Trauma Systems Branch. Prescription Drug Overdose Team, 2013]. Supplementary Table S9 summarises our approach.

For 1999-2021, we mapped ICD-10 113 cause recodes to NCHS 52 rankable causes-of-death based on Table A in [Curtin SC, 2023], and subsequently mapped the NCHS 52 rankable causes-of-death to the 53 rankable caregiver loss causes-of-death based on Supplementary Table S9.

Next, we aggregated line-list death records to annualised death counts by the 53 rankable caregiver loss causes-of-death for each year in 1983-2021, as well as by sex, age band (15-19, 20-24, 25-29, 30-34, 35-39, 40-44, 45-49, 50-54, 55-59, 60-64, 65-69, 70-74, 75-79, 80-84 and 85+ years) and standardized race & ethnicity of the decedent. Data on race and ethnicity were not available for 1983, and for this year we attributed mortality counts to standardized race categories according to the age-, sex- and cause-of-death-specific standardized race compositions in 1984. To harmonise cause-of-death data from 1983-1999 to the ICD-10 cause-of-death classifications and avoid discontinuities in our orphanhood estimates from 1998 to 1999, we used where available the comparability ratios in Table 1 of [Anderson et al., 2001]. Thirteen comparability ratios of rankable causes of death were not provided in [Anderson et al., 2001] when underlying estimations were considered imprecise, and in these cases we set the comparability ratios to 1. Supplementary Fig. S4 illustrates the aggregated annual mortality data among U.S. residents by standardized race categories, and Supplementary Fig. S5 by leading causes-of-death as relevant for caregiver loss and orphanhood.

• Download figure
• Open in new tab

Supplementary Fig. S4: Aggregated age- and sex-specific U.S. mortality counts in 1983-2021, by standardized race & ethnicity.

• Download figure
• Open in new tab

Supplementary Fig. S5: Aggregated age- and sex-specific U.S. mortality counts in 1983-2021, by leading causes-of-death as relevant for caregiver loss and orphanhood.

S1.2 National-level live births data, 1969-2021

To calculate fertility rates and attribute children experiencing orphanhood to decendents between 1983 to 2021, we required for our approach to estimating orphanhood live births data from before 1983 because the children who lost a caregiver in 1983 at age 1-17 years were born between 1966-1982. However due to limitations in publicly available population size data, we only considered in the central analysis live births data from 1990 and assumed constant fertility rates before 1990. We investigated the sensitivity of our orphanhood estimates to this assumption using live birth records since 1980 together with corresponding population size estimates, and found that orphanhood prevalence estimates since 2008 were not affected by our assumptions on historic fertility, while in the sensitivity analysis orphanhood prevalence estimates for 2000-2007 were slightly lower due to overall lower fertility rates in the 1980s (Supplementary Text S2). Line-list live births with demographic data on both mothers and fathers were available and downloaded from the NCHS Vital Statistics portal for each year between 1969 to 2021 [National Center for Health Statistics, 2023]. We considered only live birth records to U.S. resident mothers for consistency with the mortality data available. Information on residency status of fathers was not available, and we assumed fathers were also U.S. residents in the retained birth records associated with U.S. resident mothers. For each live birth, we retained year of birth, the age of mothers and fathers, and information on race and Hispanic ethnicity [National Center for Health Statistics, 2023].

Age of mothers and fathers was reported by single year of age. Following Hillis et al. [2021a], we considered live births to women aged 15-49 years and men aged 15-77 years, to match the mortality data of women aged 15-66 years and men aged 15-94 years. In total, less than 0.01% of line-list records were removed from further analysis because parent were outside of these age ranges. Age information was mapped to 5-year age bands 15-19 years, …, 45-49 years for mothers and age bands 15-19 years, …, 50-54 years, and 55 years and older (55+ years) for fathers.

Information about the race and Hispanic origin of mothers and fathers was between 1969-2021 reported in different standards due to the revision of the U.S. certificates of live births on race in 1989 and 2003 [National Center for Health Statistics, 2006]. Information on the race of mothers and fathers is publicly available since 1969, and information on their ethnicity since 1978. For 1978-2021, we mapped available information on race and ethnicity to the same standardized race & ethnicity categories used to stratify the mortality data, according to the mapping described in Supplementary Table S10.

Next, we aggregated line-list live birth records to annualized live birth counts to mothers (in age bands 15-19, …, 45-49 years) and fathers (in age bands 15-19, …, 55+ years) for each year in 1969-1977, and stratified further by standardized race & ethnicity categories for each year in 1978-2021. Before 1985, data reporting was incomplete for some U.S. states and we used NCHS sampling weights reported in each year (see Appendix A of National Center for Health Statistics [1978]) to extrapolate reported live birth counts to state populations. Supplementary Fig. S6 illustrates the aggregated annual live birth counts by mothers and fathers and standardized race & ethnicity categories since 1978.

• Download figure
• Open in new tab

Supplementary Fig. S6: Aggregated age-specific live births counts in 1969-2021, by standardized race & ethnicity categories. Prior to 1978, information on ethnicity was not publicly available and live birth counts are shown without standardized race & ethnicity stratification.

S1.3 National-level population size data, 1990-2021

To calculate fertility rates, we further required population size data of U.S. residents in the age-, sex- and standardized race & ethnicity categories of the live births data. We obtained from the CDC WONDER Vintage bridged-race postcensal & ethnicity population size estimates for each year in 1990-1999 and 2000-2020 [CDC, 2021a], and single race population size estimates for 2021 [CDC, 2021b]. Data were extracted by 5-year age bands (15-19 years, …, 70-74 years) and single years of age from 75 to 77 years, and data for individuals aged 55-77 years were summed into a single age band as required for the purposes of our analyses. Race-specific population size estimates were aggregated to the standardized race & ethnicity categories in Supplementary Table S11.

Bridged-race & ethnicity population size data were not publicly available prior to 1990. For sensitivity analyses, we also obtained national-level population size estimates by 5-year age bands, sex, and U.S. states without race & ethnicity stratification for each year in 1969-1989 from the U.S. National Cancer Institute Surveillance, Epidemiology, and End Results (SEER) Program [National Cancer Institute, 2023].

S1.4 Estimates of national-level fertility rates, 1990-2021

We calculated age-, sex- and standardized race & ethnicity-specific fertility rates in each year in 1990-2021 for women in one of the age bands a ∈ {15 − 19, …, 45 − 49} years and men in one of the age bands a ∈ {15 − 19,, 55 − 77} years according to where the number of live births and population sizes in each strata are denoted by B_y,a,s,rand P_y,a,s,r respectively. Calculations were done for each of the five standardized race categories ‘Hispanic’, ‘Non-Hispanic American Indian or Alaska Native’, ‘Non-Hispanic Asian or Pacific Islander’, ‘Non-Hispanic Black’ and ‘Non-Hispanic White’. In the central analysis, we assumed the same fertility rates in each year in 1966-1989 as in 1990, and considered alternative assumptions in several sensitivity analyses (see below). Supplementary Fig. S7 shows the calculated U.S. fertility rates by standardized race & ethnicity.

Mortality rates were calculated in analogy to Eq. (1). We found correlations by standardized race & ethnicity between fertility and mortality rates that changed primarily by age of mothers and fathers, and less so over calendar years (Supplementary Fig. S8 for women); findings were similar in men. These correlations prompted us to estimate national-level incidence and prevalence of orphanhood by standardized race & ethnicity and then sum the standardized race & ethnicity-specific estimates to obtain national-level estimates.

• Download figure
• Open in new tab

Supplementary Fig. S7: National-level fertility rates per 1,000 population by standardized race & ethnicity in 1990-2021, and assuming constant rates in 1966-1989. For comparison, the black lines show national-level fertility rates, which we calculated directly in 1969-1989 and calculated as weighted average of the standardized race & ethnicity-specific fertility rates in 1990-2021. The dashed vertical line represents 1990, where before 1990 we assumed the same fertility rates in each year in 1966-1989 as in 1990

• Download figure
• Open in new tab

Supplementary Fig. S8: Correlations between standardized race & ethnicity-specific female fertility rates and standardized race & ethnicity-specific female mortality rates by age of mothers and calendar year. Empirical fertility rates are shown as dots, and linear regression fits with 95% prediction intervals are added to visualize correlations.

S1.5 Estimates of national-level orphanhood, 1983-2021

We estimated the number of children who newly experienced orphanhood in year y, y = 1983, …, 2021, from the population-level mortality records of U.S. residents in year y and the number of children each decedent was expected to leave behind. Our approach follows the methods described in [Hillis et al., 2021a, Unwin et al., 2022], with modifications to account for causes-of-death.

We obtained the expected number of children per one U.S. resident of age a years, sex s, and standardized race & ethnicity r in year y who are of age b = 0, 1, …, 17 years (denoted with C_y,a,g,r,b) by multiplying the corresponding fertility rates of Eq. (1) with pediatric survival probabilities of children born in year y − b and surviving until age b + 1 (denoted with p^survive). Specifically, where y ranges from 1983 to 2021, the single year of age of mothers ranges from a ∈ {15, …, 66} years, the single year of age of fathers ranges from a ∈ {15, …, 94} years, the standardized race & ethnicity categories are as described before, and b = 0, 1, …, 17 years. We obtained the pediatric survival probabilities from the child mortality data in [United Nations Population Division, 2020] and use in Eq. (2) the fertility rates of the age band that includes age a − b. We then estimated the number of children aged b who newly experienced in year y the death of a parent s of age specified in 5-year age bands a^′ ∈ A = {15 − 19, · · ·, 80 − 84, 85 +years}, and standardized race & ethnicity r who died of caregiver loss cause-of-death c by Due to the correlations between standardized race & ethnicity specific fertility and mortality rates, we calculated Eq. (3) for each standardized race & ethnicity. In Eq. (3), we assumed that population-level fertility rates are not correlated with population-level mortality rates, which may lead to upwards or downwards bias in orphanhood estimates that we explored in several sensitivity analyses (see below). We also assumed that parents and their children have the same standardized race & ethnicity.

Since orphanhood considers children who experienced the death of their mother, father or both, we are interested in the sum of Eq. (3) for both mothers and fathers, but need to subtract children who lost their other parent in the previous b − 1 years, or who lost their other parent in the same year y. Throughout, we assumed that the other parent 1 − s is in the same age band a^′and of the same standardized race & ethnicity as parent s. The probability that the other biological parent of the children in Eq. (3) died in the same year y is based on standard life table calculations [Shryock et al., 1975, Lee and Wang, 2003, Arias et al., 2022], specifically the probability that an individual of sex 1 − s and standardized race & ethnicity r died in year y between (continuous) age a and a + n conditional on survival up to age a, where n = 5 corresponds to the width of the 5-year age bands considered. This mortality hazard is approximated using midpoints x = (a + a + 5)/2 in each age interval, through where for ease of readability we have suppressed y, s and r, and _nP_a and _nD_aare respectively the estimated population sizes and observed death counts by the end of the corresponding calendar year in each 5-year age band a^′ = [a, a + 5). The intermediate quantities in (4) are the approximated (unknown) mortality probability density function _nf (x) at age midpoint x and (unknown) survival function S(a) up to age a, that can be expressed in terms of the age-specific mortality rate _nq_a, the proportion of individuals alive at age a and who die before reaching age a + n in the corresponding calendar year. It is standard to estimate _nq_a with (_nD_a)/(_nP_a + ¹ _nD_a), from which Eq. (4) follows [Shryock et al., 1975, Lee and Wang, 2003, Arias et al., 2022]. Using Eq. (4), we estimated the number of children aged b who newly experienced in year y the death of one parent due to cause-of-death c and the death of the other parent due to any cause with where x is the midpoint age in age band a^′. In Eq. (5), we assumed that deaths among parents occurred independently of each other and ignored correlations of deaths among parents by the same cause-of-death such as COVID-19 [Hillis et al., 2021a], as well as correlations of deaths among parents who died of different causes-of-death. Similarly, we estimated that the number of children aged b who newly experienced in year y the death of one parent s due to cause-of-death c and the death of their other parent 1 − s due to any cause in any of the previous i = 1, …, b − 1 years with With these considerations, we estimated the number of children aged b who newly experienced orphanhood in year y = 1983, …, 2021, by the death of one or both parents of age a^′ and standardized race & ethnicity r who died of cause-of-death c with Eq. (7b) subtracts the children who lost in year y a parent due to cause c and the other parent due to any cause, which are counted twice in Eq. (7a). Without line-list family data and working from individual-level live birth and death statistics, the two terms O^new-double, O^new-double are not identical and for this reason we subtracted the average of both. Eq. (7c) subtracts the children who already lost the other parent in previous years. In our previous studies on COVID-19 associated orphanhood, we did not consider the possibility of reinfection with COVID-19 and for this reason did not subtract children who already lost the other parent due to COVID-19 in previous years in these studies [Hillis et al., 2021a]. Further arguments show that across ages, standardized race & ethnic groups, and caregiver loss causes-of-death, the values in Eq. 7b-7c are approximately equal to orphanhood prevalence divided by four, and so do not exceed 1% of the average values in Eq. 7a.

To estimate the prevalence of orphanhood in year y = 2000, …, 2021, we accrued the number of children who newly experienced orphanhood in the previous 17 years and current year y, while accounting for ageing. Specifically, for each calendar year since 2000, we estimated the total number of children aged b = 0, …, 17 years in calendar year y and of race & ethnicity r who experienced orphanhood in their lifetime by cause-of-death c in one or both parents with which sums over children who newly experienced orphanhood at younger ages in previous years.

Following Eqs. (7-8), we derived additional key quantities such as the number of children aged b in calendar year y = 1983, …, 2021 and standardized race & ethnicity r who newly experienced orphanhood in year y by parental sex s and cause-of-death c,

We chose the specific form of Eq. (9) over simpler approaches to ensure that the sum of O^new and equals Eq. 7. Further, we derived the number of children aged b in calendar year y = 2000, …, 2021 and standardized race & ethnicity r who experienced orphanhood in their lifetime by parent sex s and cause-of-death c by and the number of children of standardized race & ethnicity r who experienced orphanhood in calendar year y = 2000, …, 2021 in their lifetime by cause-of-death c in one or both parents by All other total numbers reported in this manuscript are aggregations of Eqs. (7-11).

S1.6 Estimates of national-level grandparent caregiver loss, 1983-2021

We estimated the number of children who newly experienced grandparent caregiver death in year y, y = 1983, …, 2021 from the population-level mortality records in year y of U.S. residents aged 30 years and above (30+), and assuming that each decedent who is estimated to have been a grandparent caregiver leaves a minimum of one grandchild behind. Our approach follows the methods described in [Hillis et al., 2021a], with modifications to estimate the age of children who experience the death of a grandparent caregiver.

Starting from 2010, the American Community Survey (ACS) [Ellis and Simmons, 2014] collects data on the total number of adults aged 30+ years living with grandchildren of age 17 or under in the U.S., and proportions of these by sex, and separately by bridged-race and Hispanic origin [United States Census bureau, 2021].

For each year y = 2010, …, 2021, we estimated the number of adults aged 30+ years living with grandchildren of age 17 or under for each sex and standardized race & ethnicity by multiplying totals and proportions made available by the ACS. We then divided these numbers with the respective population size estimates to estimate sex- and standardized race & ethnicity-specific proportions of adults aged 30+ years who live with their grandchildren of age 17 or under, which we denoted with γ_y,s,r. For y = 1983, …, 2009, we assumed that γ_y,s,r is the same as in 2010. We further assumed that each grandparent caregiver leaves upon death a minimum of one child behind (corresponding to Eq. (2)) and estimated the minimum number of grandchildren who newly experienced grandparent caregiver death in year y = 1983, …, 2021, with a U.S. resident grandparent caregiver aged 30+ years, sex s and standardized race & ethnicity category r who died of leading cause c by The ACS surveys did not collect data on whether both grandparents are alive and live with their grandchildren of age 17 or under, and for this reason we did not adjust Eq. (12) further for loss of other grandparents. We investigated in Sensitivity analyses our assumption that γ_y,s,r was approximately constant from 1983 to 2010 using longitudinal United Nations Population Division data on Households and Living Arrangments of Older Persons for the U.S., which suggested that the proportion of older persons who live with children or who are the primary caregivers of children has in the U.S. remained fairly constant since 1990 (Supplementary Fig. S23).

To estimate the minimum number of children who experienced grandparent caregiver death in their lifetime, we additionally need disaggregations of Eq. (12) by single year of age b = 0, …, 17 years. For the central analysis, we assumed that the age composition of G^new is the same as the age composition of children who lost parents older than 30 years, are obtained from Eq. (7) by summing over the years y = 2000, …, 2021, and c is one of the leading caregiver loss causes-of-death. Supplementary Fig. S9 shows that the age compositions of children who lost parents older than 30 years differ across causes-of-death while they are relatively stable across standardized race & ethnicity. It is plausible that the age composition of G^new may differ from theage composition of children who lost parents older than 30 years, and we explored alternative approaches to Eq. (13) in sensitivity analyses; none of these had a considerable impact on our overall estimates.

To obtain minimum estimates of caregiver loss, we then summed Eq. (9) and Eq. (13).

S1.7 Uncertainty quantification in national estimates

This section describes how we derived uncertainty intervals in the national-level estimates of orphanhood incidence and prevalence (Eqs. (7-11)) and grandparent caregiver loss (Eq. (13)). To capture uncertainty, we followed the same phenomenological approach as in [Hillis et al., 2021a], and added Poisson noise around mortality, natality, and population size data for all year, sex, age, standardized race & ethnicity, and cause-of-death strata. We generated 10,000 Poisson noise data sets, and in terms of the uncertainty in mortality data, additionally considered uncertainty related to the harmonization of causes-of-death counts from 1983-1999 to the ICD-10 classification to avoid discontinuities in orphanhood estimates between 1998 and 1999. Specifically, we resampled the comparability ratios 10,000 times, assuming these were normally distributed with the standard deviations reported in Table 1 of [Anderson et al., 2001] where available, and multiplied these ratios with the Poisson noise mortality data. Then, we repeated orphanhood incidence and prevalence estimates on each data set, and calculated 2.5% and 97.5% quantiles to generate 95% uncertainty intervals around central estimates.

In the ACS data, we accounted for uncertainty in the estimated number of adults aged 30+ years living with grandchildren and in the attribution to sex, and standardized race & ethnicity groups using published uncertainty ranges. The ACS provides for each of their published estimates margins of error that correspond to 90% confidence intervals. We converted the 90% margins of error into standard deviations and bootstrap re-sampled the available totals and separate proportions 1,000 times assuming normal distributions around their central estimates and the corresponding standard deviations, multiplied the bootstrap resampled totals and proportions, and then divided by the corresponding population sizes. Across years, we found these sources of uncertainty amounted on average to 95% bootstrap intervals on the order of ±0.79% of the central estimate of the number of Hispanic women aged 30+ years who live with grandchildren, ±8.23% among Non-Hispanic American Indian or Alaska Native women, ±1.93% among Non-Hispanic Asian women, ±0.97% among Non-Hispanic Black women and ±0.69% among Non-Hispanic White women; and similarly for men. We similarly generated 10,000 resampled data sets, repeated grandparent caregiver loss incidence and prevalence estimates on each data set, and then calculated 2.5% and 97.5% quantiles to generate 95% uncertainty intervals around central estimates.

• Download figure
• Open in new tab

Supplementary Fig. S9: Percent contributions of 1-year age bands to children newly experiencing the loss of parents older than 30 years, stratified by standardized race & ethnicity and leading caregiver loss causes-of-death.

S1.8 Estimates of state-level orphanhood, 2021

To guide policy at the state level, we further estimated orphanhood and grandparent caregiver loss in each of the 50 U.S. states and District of Columbia. We focussed on estimating the number of children who newly experienced orphanhood (incidence) and the number of children who experienced orphanhood in their lifetime (prevalence) in each state in 2021, by the leading caregiver loss causes-of-death (see Supplementary Table S9). To obtain state-level orphanhood prevalence estimates for 2021, as before we had to accrue and age state-level orphanhood incidence estimates from 2004, and to obtain state-level fertility rates from 1987. Overall, we proceeded analogously as for the national estimations, but obtained data from CDC WONDER [CDC, 2023] as state-level mortality data were not publicly available from NCHS from 2005 onwards [National Center for Health Statistics, 2023]. Counts below 10 were suppressed in CDC WONDER and for this reason we performed all estimations without stratification by standardized race & ethnicity. We adjusted for data suppression, and known biases in this approach that are due to correlations in mortality and fertility rates across standardized race & ethnicity. We note these are modelled adjustments and we cannot exclude that this resulted in bias in state-specific orphanhood estimates.

We extracted annual death counts by state, sex, age bands (15 − 19 years, …, 95 − 99 years, > 100 years), and cause-of-death (ICD-10 113 Selected Causes of Death) from the CDC WONDER mortality portal [CDC, 2023] from 2005 to 2021; and combined these data with the previously described, cause-specific NCHS mortality data sets from 2004 for which information on the state of residence of the decedent was publicly available. Throughout, we focused only on the eight leading caregiver-loss causes of death chronic liver disease and cirrhosis, COVID-19, diseases of the heart, drug overdose, homicide excluding drug overdose, malignant neoplasms, suicide excluding drug overdose, and unintentional injuries excluding drug overdose. For the leading caregiver-loss causes of death, the average discrepancy in the state-, sex-, age-, and cause-of-death-specific data from CDC WONDER was, when aggregated to national-level and compared to the corresponding NCHS mortality data, 22.9% among women and 13.5% among men. As a first step, we imputed suppressed entries with a value of 2, which resulted in average discrepancies of 6.80% among women and 2.53% among men. To account for these discrepancies, as a second step, we adjusted the state-level mortality counts according to where y in 2005 to 2021, a^′ denotes the age band of the parent, s the sex of the parent and c the leading caregiver loss cause-of-death of the parent, and D_y,a′_,s,c are the national-level deaths derived from the line-list NCHS data in Eq. (3), and D_y,a′_,s,l,c are the state-specific deaths derived from CDC WONDER for each of the 50 U.S. states and District of Columbia, indexed by l. The suppression adjusted, cause-specific, state-level mortality counts matched the cause-specific, national-level mortality counts well, except when all deaths D_y,a′_,s,c,l were entirely suppressed across all states, which was the case for deaths due to homicide in women aged 55 years or older (Supplementary Fig. S10, compare pink versus blue). Overall, these suppression adjustments did not noticeably change the contribution of causes-of-death to mortality, because the majority of death counts were not suppressed when the data were aggregated to leading caregiver loss causes-of-death (Supplementary Fig. S11).

We next extracted live birth counts for mothers by state and age bands (15 − 19, …, 44 − 49 years) from the CDC WONDER natality portal [CDC WONDER, 2022] from 2005 to 2021; and combined these data with the previously described NCHS live birth data sets from 1987 to 2004 for which information on the state of residence of mothers was available. We proceeded analogously for fathers using the age bands 15 − 19, …, 44 − 49, 50 − 54, 55+ years from 2016 to 2021, as natality records were not available by demographic characteristics of fathers from 2005 to 2015. Data suppression was not an issue for the stratifications required to estimate orphanhood: for the years 2005-2021, the average discrepancy in the CDC WONDER and NCHS natality data sets was 0.026% for mothers and 0.032% for fathers (as compared to male NCHS live birth data with ages up to 77 years). To compute fertility rates at the state level, we further extracted age- and sex-stratified annual population size estimates for each state from 1987 to 1989 from the NIH National Cancer Institute Surveillance, Epidemiology, and End Results (SEER) database [National Cancer Institute, 2023]; and CDC WONDER from 1990 to 2021 [CDC, 2021c, 2022]. State-level fertility rates were calculated as in Eq. (1) when live births counts were available, but without stratification by race & ethnicity. For mothers, live birth counts were not publicly available for women aged 44-49 years in several years and 8 states (Alaska, Delaware, Montana, North Dakota, South Dakota, Vermont, West Virginia and Wyoming), and in these cases we interpolated state-specific fertility rates with locally estimated scatterplot smoothing (LOESS) as implemented in the R stat package version 4.2.3 with span argument 0.85. For Wyoming, no data were publicly available after 2019 and so interpolation was not possible and we used 2019 values. For fathers, live births counts were not publicly available from 2005 to 2015, and for these years we estimated state-specific fertility rates with LOESS with span argument 0.85, and using NCHS data from 2000 to 2004 and CDC WONDER data from 2016 to 2020 (Supplementary Fig. S12).

We then estimated state-specific incidence and prevalence of orphanhood from the suppression-adjusted state-level, cause-specific mortality counts and partially imputed fertility rates as outlined in Eqs. (7-11), with one modification. As described above the state-level mortality and fertility data were not disaggregated by race & ethnicity. To account for potential bias arising from correlations in mortality and fertility rates across race & ethnicity, Overall, the correction factors ν_y,s,c tended to be close to one, except for Homicide excluding drug overdose in women and ‘Other’ caregiver loss causes-of-death (Supplementary Fig. S13). We then adjusted the state-specific estimates as follows, we compared the national-level estimates of the number of children who newly experienced orphanhood in (9) to the state-specific estimates with the correction factors and used Eq. (16) to calculate the state-specific analogues to Eqs. (8-11).

• Download figure
• Open in new tab

Supplementary Fig. S10: Completeness of state-level mortality counts. State-level mortality counts were extracted from CDC WONDER by year, sex, age band, and cause-of-death; and subject to suppression of counts below 10 individuals. Counts were aggregated to national level to assess their completeness as described in the Supplementary Text. Shown are the cause-specific, state-level mortality counts extracted CDC WONDER, the suppression-adjusted, cause-specific, state-level counts, and for comparison the cause-specific, national-level mortality counts obtained from NCHS.

We estimated the number of children who newly experienced grandparent caregiver death in year y, y = 2004, …, 2021 from the state-level mortality records in year y of U.S. residents aged 30 years and above (30+), and assuming that each decedent who is estimated to have been a grandparent caregiver in each state leaves a minimum of one grandchild behind. Overall, we proceeded as for the national-level estimation of grandparent caregiver loss as ACS also published state-specific data on adults aged 30+ years who live with grandchildren of age 17 or under [United States Census bureau, 2021]. For each year y = 2010, …, 2021, we calculated the expected number of adults aged 30+ years who live with their grandchildren of age 17 or under by multiplying the corresponding total co-resident numbers with the sex-specific proportions for each state; and then divided these with corresponding population sizes to obtain the corresponding proportions γ_y,s,l. Overall, these proportions were considerably more uncertain than in the national-level analysis due to smaller ACS sample sizes (see below). For y = 2004, …, 2009, we assumed that γ_y,s,lis the same as in 2010. We further assumed that each grandparent caregiver leaves upon death a minimum of one child behind (corresponding to Eq. (2)) and estimated the minimum number of grandchildren who newly experienced grandparent caregiver death in year y = 2004, …, 2021, with a U.S. resident grandparent caregiver aged 30+ years, sex s and state category l who died of leading caregiver loss cause-of-death c by

To estimate the minimum number of children who experienced grandparent caregiver death in their lifetime in each state, we additionally need disaggregations of Eq. (17) by single year of age b = 0, …, 17 years. For the central analysis, we assumed that the age composition of G^new is the same as the age composition of children who lost parents older than 30 years in the same state, are obtained from Eq. (16) by summing over the years y = 2004, …, 2021, and c is one of the leading caregiver loss causes-of-death. Supplementary Fig. S14 shows that the age compositions of children who lost parents older than 30 years differed across causes-of-death while they are relatively stable across states.

• Download figure
• Open in new tab

Supplementary Fig. S11: Contribution of leading causes-of-death to U.S. deaths by states. State-level mortality counts were extracted from CDC WONDER by year, sex, age band, and cause-of-death and adjusted for suppression. Shown are the contributions of leading causes-of-death to the unadjusted state-specific deaths in 2021, and the adjusted state-specific deaths in 2021.

• Download figure
• Open in new tab

Supplementary Fig. S12: Interpolation of missing state-specific male fertility rates. State-level live birth counts were extracted from CDC WONDER by year, sex, and age band; and combined with line-list NCHS live birth data up to 2004 when information on state residence was made available. For fathers, demographic information was only available since 2016 in CDC WONDER and state-specific male fertility rates were interpolated with LOESS estimators as described in the Supplementary Text. Shown are empirical male fertility rates (dots) for 5 states from 1987 to 2004 and from 2016 to 2021, along with the LOESS estimates (lines) for 2000 to 2020 and 95% prediction intervals (ribbon).

S1.9 Uncertainty quantification in state-level estimates

This section describes how we derived uncertainty intervals in the state-level estimates of orphanhood incidence and prevalence (Eqs. (7-11)) and grandparent caregiver loss (Eq. (13)). Overall, the approach here is similar to the national-level uncertainty analysis.

We added Poisson noise around state-specific mortality, natality, and population size data for all year, sex, age, and cause-of-death strata. We generated 10,000 Poisson noise data sets, repeated orphanhood incidence and prevalence estimates on each data set, and then calculated 2.5% and 97.5% quantiles to generate 95% uncertainty intervals around central estimates.

To estimate grandparent caregiver loss, we accounted for uncertainty in the estimated number of adults aged 30+ years living with grandchildren by U.S. state using published uncertainty ranges. The ACS provides for each of their state-specific estimates margins of error that correspond to 90% confidence intervals. We converted the 90% margins of error into standard deviations and bootstrap re-sampled the available totals and separate proportions 10,000 times assuming normal distributions around their central estimates and the corresponding standard deviations, and then multiplied the sampled totals and proportions. As expected given available ACS sample sizes, across years we found that uncertainty in the proportions γ_y,s,l was in some states considerably larger than for the national race & ethnicity-specific analysis, with 95% bootstrap intervals of up to ±13.5% of the central estimates (Supplementary Fig. S15). We similarly generated 10,000 resampled data sets, repeated grandparent caregiver loss incidence and prevalence estimates on each data set, and then calculated 2.5% and 97.5% quantiles to generate 95% uncertainty intervals around central estimates.

S1.10 Estimates of racial & ethnic groups impacted by leading causes of caregiver loss in U.S. states, 2021

Characterizing the racial & ethnic groups of children that are impacted by orphanhood and grandparent caregiver loss at state-level is challenging due to limits in publicly available data. We focused on estimating and comparing 2021 orphanhood and grandparent caregiver loss prevalence rates by standardized race & ethnicity in the 10 U.S. states with highest orphanhood prevalence rates and for the primary (first-ranked) caregiver-loss cause-of-death only, West Virginia (‘Drug overdose’), New Mexico (‘Drug overdose’), Mississippi (‘Unintentional injuries’), Louisiana (‘Drug overdose’), Kentucky (‘Drug overdose’), Tennessee (‘Drug overdose’), Alabama (‘Diseases of heart’), Oklahoma (‘Unintentional injuries excluding drug overdose’), Ohio (‘Drug overdose’) and Florida (‘Drug overdose’). Throughout, we obtained publicly available data from CDC WONDER [CDC, 2023] from 2005 onwards [National Center for Health Statistics, 2023] and adjusted for data suppression, but note these are modelled adjustments and we cannot exclude that this resulted in bias in state-specific orphanhood estimates.

Annualized live birth counts were extracted as described for the state-level analysis for 2005 to 2019, but now also stratified by bridged-race and Hispanic origins, and then combined with NCHS live births data from 1987 to 2004. For 2020-2021, live birth counts to women were not available by bridged-race at state-level. Analogously, annual death counts for the leading caregiver loss causes-of-death were compiled as described above for 2005 to 2021, and stratified into the standardized race & ethnicity categories in Supplementary Table S7, and then combined with the NCHS mortality data sets from 2004. Suppressed values were imputed by 1 and adjusted further for suppression as in Eq. (14). We followed previous criteria on data reliability and excluded from consideration for each state those race & ethnicity groups that had more than two age bands with fewer than 20 live birth counts Mathews and Hamilton [2019]. Supplementary Table S6 lists the corresponding strata as ‘small populations’, and also describes the average completeness of live birth records for the years 1995-2004 when directly comparable line-list live birth records were also available from NCHS. For the remaining year-state-sex-race & ethnicity strata, we considered the corresponding mortality data for the leading caregiver-loss cause-of-death. Again, we excluded from consideration those standardized race & ethnicity groups that had more than two age bands with fewer than 20 death counts Mathews and Hamilton [2019]. Supplementary Table S6 lists the corresponding strata as ‘small death counts’, and also describes the average completeness of death records for the years 1999-2004 when directly comparable line-list death records were also available from NCHS. Supplementary Fig. S16 illustrates the state and standardized race & ethnicity specific live birth data from CDC WONDER, the line-list live birth data from NCHS and the state and standardized race & ethnicity specific live birth data for the three states with highest orphanhood prevalence rates in 2021. Similarly, Supplementary Fig. S17 illustrates the state and standardized race & ethnicity specific mortality data from CDC WONDER for the primary caregiver loss cause-of-death, the line-list live death from NCHS and the adjusted state and standardized race & ethnicity specific mortality data for the three states with highest orphanhood prevalence rates in 2021. Population-size estimates for each state were obtained from the CDC WONDER population size portal for 1990 to 2021 [CDC, 2021c, 2022] by bridged-race, and converted to the standardized race & ethnicity categories described in Supplementary Table S11.

• Download figure
• Open in new tab

Supplementary Fig. S13: Comparison of estimated national-level incidence of orphanhood and estimated state-level incidence of orphanhood. At the national level, incidence of orphanhood was estimated from line-list live birth and death data stratified by race & ethnicity; whereas as the state-level, incidence of orphanhood was estimated from aggregated and partially interpolated live birth data and partially suppressed cause-specific deaths that were not stratified by race & ethnicity. Shown are incidence estimates using either approach by leading caregiver loss causes-of-death. Overall, incidence estimates were very similar except for ‘Homicide excluding drug overdose’ in women, ‘Chronic liver disease and cirrhosis’ in women and ‘Others’ causes-of-death. State-level incidence estimates were adjusted to match the national-level totals as described in the Supplementary Text.

• Download figure
• Open in new tab

Supplementary Fig. S14: State-specific percent contributions of 1-year age bands to children newly experiencing the loss of parents older than 30 years. Percent contributions were stratified by sex and leading caregiver loss causes-of-death of parents (facets) and are shown for each U.S. state (colours).

• Download figure
• Open in new tab

Supplementary Fig. S15: Magnitude of uncertainty in the number of adults aged 30+ who live with grandchildren in each U.S. state. Estimates of uncertainty were derived using bootstrap resampling from the margins of error reported by ACS for each of their estimates. Uncertainty for women and men is shown separately in colour.

• Download figure
• Open in new tab

Supplementary Fig. S16: Data completeness evaluation for race & ethnicity- and state-specific live birth counts for three U.S. states, West Virginia, New Mexico and Mississippi. Standardized race & ethnicity live birth counts were extracted from CDC WONDER by state, year, sex, and age band; and subject to suppression of counts below 10 individuals. Counts from NCHS were aggregated to state-level by standardized race & ethnicity to assess their completeness as described in the Supplementary Text. Shown are the state-level live birth counts extracted from CDC WONDER, the suppression-adjusted, standardized race & ethnicity-specific counts, and for comparison the live birth counts obtained from NCHS by the same stratification. Points shown in opaque colours represent standardized race & ethnicity categories that were excluded from further analysis due to the extent of suppressed data.

• Download figure
• Open in new tab

Supplementary Fig. S17: Data completeness evaluation for race & ethnicity- and state-specific mortality counts for three U.S. states, West Virginia, New Mexico and Mississippi. Standardized race & ethnicity-specific mortality counts were extracted from CDC WONDER by state, year, sex, age band, and leading caregiver loss cause-of-death; and subject to suppression of counts below 10 individuals. Counts were aggregated to state-level by standardized race & ethnicity to assess their completeness as described in the Supplementary Text. Shown are the mortality counts extracted from CDC WONDER, the suppression-adjusted, race & ethnicity-specific counts, and for comparison the mortality counts obtained from NCHS by the same stratification. Points shown in opaque colours represent standardized race & ethnicity categories that were excluded from further analysis due to the extent of suppressed data.

We estimated race & ethnicity- and state-specific incidence and prevalence of orphanhood as outlined in Eqs. (7-11), except that female race & ethnicity- and state-specific fertility rates in 2020-2021 were as assumed to be as in 2019 due to limitations in publicly available data. To investigate estimation accuracy, we compared the resulting sum of the race & ethnicity- and state-specific orphanhood incidence estimates to the previous state-specific orphanhood incidence estimates (Supplementary Fig. S18). For Alaska and Oklahoma, the sum of the race & ethnicity- and state-specific orphanhood incidence estimates was more than 20% below the state-specific orphanhood incidence estimates attributable to the leading caregiver loss cause-of-death. Supplementary Table S6 lists these states as ‘large discrepancy in estimates’. The table then reports, for all remaining states, 2021 race & ethnicity- and state-specific orphanhood prevalence rate estimates that are attributable to the leading caregiver loss cause-of-death.

• Download figure
• Open in new tab

Supplementary Fig. S18: Reliability evaluation of orphanhood incidence estimates by state, year, standardized race & ethnicity and leading caregiver loss cause-of-death. Race & ethnicity- and state-specific orphanhood incidence rates for the leading caregiver loss cause-of-death were estimated from suppression-adjusted live births and deaths counts from CDC WONDER as described in the Supplementary Text. To assess the reliability of these estimates, we considered weighted averages of the orphanhood incidence estimates across race & ethnicity (yellow triangles), and then compared the weighted averages against state-level estimates (blue circles). Triangles and circles are shown in opaque colour when race & ethnicity-specific estimates disagreed by more than 20% and were considered unreliable.

Sensitivity in mortality data

In the central analysis, we derived mortality counts by year, sex, age band, standardized race & ethnicity and caregiver loss causes-of-death from NCHS line-list mortality records from 1983 to 2021. Aggregate mortality counts are also available from CDC WONDER CDC [2023] from 1999 to 2021. We extracted data from the NCHS Vital Statistics portal because this allowed us to use the same data source across all years, bypass data suppression, and incorporate uncertainties in standardized race & ethnicity reporting. The CDC WONDER mortality counts allowed us to check our in-house data aggregations, if we obtain data from CDC WONDER at coarser population strata. For this purpose, we extracted annual death counts by sex, and age bands (15-19 years, …, 95-99 years, >100 years) from the CDC WONDER mortality portal CDC [2023] from 2000 to 2021. The overall mortality counts that we aggregated from the NCHS line list data were identical with the CDC WONDER mortality data without cause-of-death stratification. Secondly, we extracted annual death counts by sex, age bands (15-19 years, …, 95-99 years, >100 years) and cause-of-death (ICD-10 113 Selected Causes of Death) from the CDC WONDER mortality portal CDC [2023] from 2000 to 2021, without stratification by race & ethnicity. We then compared the mortality counts from the two data sources for each year and the leading caregiver loss causes-of-death (Supplementary Fig. S19). We found that across years, the maximum discrepancy in the CDC WONDER counts relative to the NCHS mortality counts was 0.0068% in women and 0.0057% in men, which reflected data suppression especially due to homicide excluding drug overdose in women, and overall indicated consistency of our data with that from CDC WONDER.

Sensitivity in live births data

In the central analysis, we derived live birth counts by year, sex, age band, and standardized race & ethnicity from NCHS line-list natality records from 1969 to 2021. Aggregate live birth counts are also available from CDC WONDER CDC WONDER [2022] from 1995 to 2021. We chose to extract data from the NCHS Vital Statistics portal because this allowed us to use the same data source across all years. The CDC WONDER live birth counts allowed us to check our in-house data aggregations. For this purpose, we extracted annual live birth counts from the CDC WONDER natality portal for women by age bands 15-19 years, …, 45-49 years from 2000 to 2021, and separately for men by age bands 15-19 years, …, 55+ years from 2016 to 2021 without stratification by race & ethnicity, which avoids data suppression. We then compared the live birth counts from the two data sources for each year and found that both data sets matched across all years.

Sensitivity in national-level orphanhood estimates to assumptions on historic fertility rates

In the central analysis, we assumed that male and female fertility rates by age and standardized race & ethnicity were in 1966-1989 constant and as in 1990 (Supplementary Fig. S7). This assumption was made because population size data stratified by detailed race categories (more than three categories) were not publicly available before 1990. Considering that age- and sex-specific live births data were available by race & ethnicity since 1978, we estimated in this sensitivity analysis the historic composition of population sizes by race & ethnicity in 1980-1989, and then updated the corresponding fertility rates and national-level orphanhood incidence and prevalence estimates. More specifically, we considered time trends in the proportion of each race & ethnicity in each 5-year age band and both sexes in the U.S. intercensal population size estimates in 1990 to 1995 (Supplementary Fig. S20A), and extrapolated these with linear models backwards in time to 1980-1989. We do not think that this estimation approach resulted in accurate estimates of population sizes, but rather that this approach conveys possible sensitivities in our orphanhood estimates. We identified minor sensitivities to national orphanhood prevalence estimates up to 2007, with no impact on incidence and prevalence estimates for recent years (Supplementary Fig. S20).

Sensitivity in national-level orphanhood estimates to potentially correlated fertility rates

In the central analysis, we assumed that population-level fertility rates were not correlated with population-level mortality rates. We considered in sensitivity analyses possible deviations from this assumption that might lead to upwards bias in our central estimates. For example, a person who died in year y may have experienced poor health in the preceding years y − 1, y − 2, … and in this case may also have been less likely to mother or father a child in years y, y − 1, y − 2, We modelled this possibility through four scenarios of dampened fertility rates close to individual death events. Specifically, we considered cumulative logistic adjustment factors that dampened fertility rates to either zero or half of the corresponding population-level, year-, age-, sex-, and standardized race & ethnicity-specific fertility rates in the year of death y. We also considered sensitivity analyses so that the onset of lower fertility rates preceded the year of death y by 1 or 3 years as shown in Supplementary Fig. S21A. We found that orphanhood incidence estimates were up to 8.4% lower than in the central analysis, and orphanhood prevalence estimates were up to 15% lower than in the central analysis (Supplementary Fig. S21B). It is also possible that our central orphanhood estimates might be biased downward. For example, women experiencing premature death are more likely economically disadvantaged and in turn are more likely to have had reduced access to contraceptives and corresponding higher fertility rates [Price et al., 2013], or earlier sexual debut resulting in higher cumulative fertility until death. It is also possible that overall, at population-level, there is no strong correlation between fertility and mortality rates especially as the lag between birth and death events is typically relatively large, and in the absence of data we have opted for this middle approach in the central analysis.

• Download figure
• Open in new tab

Supplementary Fig. S19: Comparison of national-level U.S. mortality counts from NCHS and CDC WONDER. In the central analysis, we derived mortality counts by year, sex, age band, standardized race & ethnicity and caregiver loss causes-of-death from NCHS line-list mortality records from 1983 to 2021. To check consistency of our data aggregations, we extracted annual death counts by sex, age bands (15-19 years, …, 95-99 years, >100 years) and cause-of-death (ICD-10 113 Selected Causes of Death) from CDC WONDER (i.e. without further stratification by standardized race & ethnicity). Mortality counts from both data sources were identical except for deaths due to homicide excluding drug overdose in women, and overall indicated consistency of our data with that from CDC WONDER.

Sensitivity in national-level grandparent caregiver loss estimates to assumptions on the age of children experiencing loss of a grandparent caregiver

In the central analysis, we considered in Eq. (13) as a proxy to the age composition of children experiencing grandparent caregiver loss the age composition of children experiencing orphanhood. Calculations were done independently for each leading caregiver loss cause-of-death. We performed two analyses to characterise the sensitivity of this approach to national-level grandparent caregiver loss estimates. First, we repeated calculations using as proxy the age composition of children experiencing orphanhood across all causes of caregiver loss, i.e.

• Download figure
• Open in new tab

Supplementary Fig. S20: Sensitivity in national-level orphanhood estimates to assumptions on historic fertility rates. (A) Composition of U.S. population by standardized race & ethnicitiy (fill colours) from 1990 to 2021, for both men and women in the age bands relevant for orphanhood estimation. (B) Estimated national-level orphanhood prevalence rates by standardized race & ethnicity (colours) in the sensitivity analysis (dashed) and central analysis (solid line). We identified minor sensitivities to national orphanhood prevalence estimates up to 2007, with no impact on incidence and prevalence estimates for recent years.

Second, we used data from the ‘Topical survey’ of the National Survey of Children’s Health (NSCH) [US Census Bureau, 2023a], filtered to respondents reporting to be grandparents and living with at least one child in the household. Respondents were asked to report on demographic characteristics including the age and race & ethnicity of one randomly selected child. We pooled these data across six years from 2016 to 2021 to characterise the age composition of children who live with grandparents due to small sample sizes in each survey after filtering (approximately 2, 500 grandparents per survey round). Supplementary Fig. S22A compares the age compositions used in the central analysis to those in the two sensitivity analyses. The comparison of the grandparent caregivers prevalence estimates are shown between the central estimate and two sensitivity analyses by age of child in Supplementary Fig. S22B. We found that grandparent caregiver incidence estimates deviated by up to ±0.028% and ±0.032% of the central estimate in the two sensitivity analyses due to the rounding numbers, respectively. Additionally, in the two sensitivity analyses, grandparent caregiver prevalence estimates were respectively 12.0% fewer and 13.6% higher than in the central analysis.

• Download figure
• Open in new tab

Supplementary Fig. S21: Sensitivity in national-level orphanhood estimates to assumptions on correlations between fertility and mortality rates. (A) We considered four scenarios of lower fertility rates preceding a death event, which we implemented through time-dependent multiplication factors to fertility rates by the time prior to a death event. (B) Estimated national-level orphanhood prevalence rates by race and ethnicity (facets and colours) in the sensitivity analyses (dashed, dotted, dotdash and twodash lines) and central analysis (solid line).

• Download figure
• Open in new tab

Supplementary Fig. S22: Sensitivity in national-level grandparent caregiver loss estimates to assumptions on the age of children experiencing loss of a grandparent caregiver. (A) Proxies for the age distribution of children who lost grandparent caregivers in two sensitivity analyses. In the first sensitivity analysis (dashed line), we computed the age distribution of children experiencing the death of a parent aged 30+ years across all causes-of-death and by race & ethnicity of the decedent parent. In the second sensitivity analysis (dotted line), we computed the age distribution of randomly selected children living with a grandparent caregiver by race & ethnicity of the child from NSCH surveys in 2016-2021. For comparison, we show the age distribution of children experiencing the death of a parent aged 30+ years due to drug overdose (solid line). (B) Estimated national-level grandparent caregivers loss incidence by age of child (colours) in the two sensitivity analyses (dashed, dotted line) and central analysis (solid line). Labels in colour represent the age of child.

Sensitivity in national-level caregiver loss estimates to assumptions on historic numbers of grandparent caregivers

In the central analysis, we assumed that the proportion of adults aged 30+ years who live with their grandchildren of age 17 or under (denoted γ_y,s,r) was the same in the years y = 1983, …, 2009 as in 2010. This assumption was made because data were not available prior to 2010. To investigate the implications of this assumption, we considered longitudinal United Nations Population Division data on Households and Living Arrangments of Older Persons for the U.S. since 1960 United [2022]. These data indicate that the proportion of older persons who live with children or who are the primary caregivers of children has in the U.S. remained fairly constant since 1980 (Supplementary Fig. S23), and for this reason suggests that our assumptions on the historic number of grandparent caregivers is unlikely to have a substantive impact on caregiver loss estimates.

• Download figure
• Open in new tab

Supplementary Fig. S23: Time trends in the proportion of U.S. residents older than 60 years who are living or caring for children. Data were obtained from United Nations Population Division, data on Households and Living Arrangments of Older Persons between 1960 and 2015. The dashed line represents the year since when grandparent caregiver loss estimates were generated, 1983.