A new method to attribute differences in total deaths between groups to population size, age structure and age-specific mortality rate

Background Two decomposition methods have been widely used to attribute death differences between two populations to population size, age structure of the population, and age-specific mortality rate (ASMR), but their properties remain uninvestigated. Methods We assess how the two established decomposition methods yield varying results with three-factor factorial experimental designs, illustrating that they are sensitive to the choice of the reference group. We then propose a novel decomposition method to obtain robust decomposition results and use three cases to demonstrate its advantage. Results The three decomposition methods differ fundamentally in their allocation of interactions to the contributions of the three factors. In comparison with the existing methods, the new method is robust to the choice of the reference group. Three case studies showed inconsistent attribution results for the two existing methods but robust results for the new method when the choice of the reference population changes. Conclusions The proposed method offers robust and more justifiable attribution results compared to the two existing methods. This method could be generalized to attribution of group differences of other health indicators.


Background
As population aging accelerates worldwide, it is valuable for researchers to quantify the impact of aging on the burden of diseases at global, regional, national and local levels. One important aspect of quantifying the effect of population aging is to attribute death differences between two populations to their differences in population size, age structure, and age-specific mortality rates (ASMR). Such attribution (or decomposition) analysis is key to interpreting death differences across time and place and to developing long-term strategies to improve public health at different levels.
Two approaches have been established and widely used to estimate the contribution of population aging to changes in the total number of population deaths over time [1,2]. We refer to the approach developed by Bashir et al [1]. as method I and the approach that the Global Burden of Disease (GBD) study group developed by extending the method of Das Gupta for decomposing rate differences as method II [2,3]. Both approaches ascribe differences in total deaths to the contributions of three factors: (a) age structure of the population (a proxy for population aging), (b) population size, and (c) age-specific mortality rates (ASMR) due to all other reasons (e.g. motorization, urbanization, poverty, environmental pollution, individual behaviors, policy interventions, genetic disorders) [1,2,4].
The two approaches provide decompositions with intuitive interpretations in many applications, but some inherent issues with the approaches have not been addressed in the literature. First, as we demonstrate below, the two approaches do not consider how to allocate two-way and three-way interactions between the three factors when decomposing their contributions. Previous studies indicate the presence of an interaction between age structure and ASMR when the mortality rate difference between two populations is attributed to the two factors [5][6][7]. Second, both methods fail to provide users with useful guidance on the selection of the reference population for a decomposition analysis, but their results heavily depend on such selection. These limitations have not been discussed before despite the fact that both methods are widely used to support health-related decision-making [4,8,9].
In this study, we demonstrate that methods I and II generate inconsistent results for the same data when the reference population changes and then propose an alternative decomposition method that overcomes these shortcomings. We apply all three methods to the decomposition of differences in total deaths from unintentional falls (a) between the United States and

Methods and results
We use formula derivation to elaborate the shortcomings of methods I and II and the development of new method. Suppose we are to decompose the difference in total number of deaths between two populations (j = 1, 2). The two populations could be defined by time periods, geographic places, or both. Each population has p age groups (i = 1, 2, . . ., p). Let d ij , n ij , and m ij denote the number of deaths, population size, and ASMR for age group j of population i, respectively. Let s ij represents the proportion of the i th age group among the j th population, (i = 1, 2, . . ., p, j = 1, 2) ( Table 1).

The two existing decomposition methods
Method I uses six steps to decompose the difference in total deaths between the two populations (D 2 -D 1 ) (note: the first population is selected as the reference for this example) [1]: i. Scale the size of each population to 100 000 persons while retaining the original age structure.
ii. Based on the scaled population size (100 000 persons), use original ASMRs of each population to calculate the expected total number of deaths for population 1 (D 1e ) and population 2 (D 2e ), respectively.
iii. Apply ASMRs from population 1 to the scaled 100 000 persons for population 2 to estimate the expected total number of deaths (DASMR 1 2e) for population 2.
With AS I , ASMR I and PS I representing deaths attributed to age structure, ASMR and population size defined by method I, respectively, we derive the following formulas according to the literature (details are provided in the S1 File) [1]: Method II differs slightly from method I by decomposing (D 2 -D 1 ) according to the following five steps (note: the first population is again selected as the reference) [2]: i. Apply the ASMRs and age structure from population 1 to population 2 to calculate the expected total number of deaths, D(AS 1 , ASMR 1 ) 2e, which denotes the expected deaths of population 2 when it has the same age structure and ASMRs as population 1.
ii. Apply age-specific mortality rates from population 1 to population 2 to calculate the expected total number of deaths, DASMR 1 2e, which represents the expected deaths of population 2 when it has the same ASMRs as population 1.
Note: d ij , n ij , m ij , and s ij represent the number of deaths, population size, age-specific mortality rate and proportion of group population among total population for the ij th subgroup. D 1 and D 2 , N 1 and N 2 , M 1 and M 2 represent the total number of deaths, population and crude mortality rate for groups 1 and 2, respectively. https://doi.org/10.1371/journal.pone.0216613.t001 iii. (D(AS 1 , ASMR 1 ) 2e -D 1 ) reflects the contribution of population size.
Using AS II , ASMR II and PS II to represent death differences attributed to age structure, ASMR and population size defined by method II, respectively, we have the following formulas according to the literature (note: calculation details are included in the S1 File) [2]: Factorial representation of decomposition methods I and II. The decomposition process can be represented in analogy to the analysis of a three-factor factorial experiment design. We can attribute the difference in total deaths between two populations to the main effect, three two-way interactions and one three-way interaction of population size, age structure, and ASMR. Using population 1 as the reference, the main effects (M p , M s and M m ), the twoway interactions (I ps , I pm , I sm ), and the three-way interaction (I psm ), where the subscripts p, s and m stand for population size, age structure and ASMR respectively, can be represented as follows: Applying the above expressions, the contributions of the three factors defined by methods I and II can be written as: Clearly, the allocation schemes of interaction terms to the three factors differ between methods I and II, but are asymmetric for both methods, which leads to the non-robustness of attribution results. A closer look at the terms in (7)-(13) reveals that, when the reference group is switched, all main effects will change both their signs and magnitudes, the two-way interactions will change their magnitudes but not their signs, and the three-way interaction will change their signs but not their magnitudes. For example, the formula for the main effect of population size will change from formula (7) to formula (16) when the reference is altered from population 1 to population 2. As a result of the asymmetric allocation of the interactions, the underlying contribution expressions of population size, age structure and ASMR, as shown in formulas (14) and (15), will change accordingly when the interactions exist and the reference population is changed from population 1 to population 2.
Evaluating the robustness of methods I and II. A simple criterion to evaluate the robustness of methods I and II is to check whether the absolute contributions remain the same but with opposite sign when the reference population changes. Following the principles of the two methods [1,2], we summarize the formulas for calculating the contributions of the three factors under the two choices of reference populations in Table 2. Detailed derivations of these formulas are provided in the S1 File. These formulas indicate that, for both methods, decomposition results vary with the choice of reference population.
Three illustrative case studies. We present three examples to illustrate how decomposition results depend on the choice of reference population for methods I and II. Specifically, we decomposed the difference in total number of deaths from unintentional falls: (a) between the United States and China in 2017, (b) between 1990 and 2017 in the United States, and (c) between 1990 and 2017 in China. All data were derived from the GBD Study 2017 update (GBD 2017) [10]. Age was divided into ten groups: <5, 5-14, 15-24, 25-34, 35-44, 45-54, 55-64, 65-74, 75-84, and �85 years ( Table 3). Data analysis was performed using R 3.4.0.
Case study one. In 2017, approximately 134,773 unintentional fall deaths occurred in China and 38,368 unintentional fall deaths occurred in the United States. The decomposition results of the three factors between China and United States in 2017 changed greatly when the reference population was changed from China to the United States. For method I, the contributions of ASMR, age structure and population size changed from -28% to 6% (deaths: from -27,190 to 6,187), 61% to -14% (deaths: from 59,253 to -13,561), and -133% to 108% (deaths: from -128,467 to 103,778), respectively (note: the signs of contributions reverse when the reference population is switched) (Fig 1A1). We observed changes not only in the signs but also in the magnitudes of contributions. Similarly, the contributions of the three factors changed respectively from -6% to 28% (deaths: from -6,253 to 26,903), 14% to -61% (deaths: from 13,627 to -58,966), and -108% to 133% (deaths: from -103,778 to 128,467) for method II (Fig 1A2).
Case study two. Between 1990 and 2017, the number of unintentional fall deaths rose from 64,102 to 134,773 in China. The estimated contributions of the three factors to the difference  between 1990 and 2017 in China changed greatly when we switched the reference year from 1990 to 2017. For method I, the contributions of ASMR, age structure and population size changed respectively from 26% to 3% (deaths: from 18,079 to 2,096), 45% to -87% (deaths: from 32,030 to -61,226), and 29% to -16% (deaths: from 20,562 to -11,541) ( Fig 1B1). For Method II, the contributions changed from 30% to 3% (deaths: from 21,333 to 1,776) for ASMR, 53% to -73% (deaths: from 37,796 to -51,885) for age structure, and 16% to -29% (deaths: from 11,541 to -20,562) for population size (Fig 1B2).
Case study three. For the United States, the contributions of the three factors to the differential numbers of unintentional fall deaths between 1990 and 2017 estimated by methods I and II are also sensitive to the choice of the reference year. For method I, the contributions changed from 52% to -48% (deaths: from 13,003 to -12010) for ASMR, 15% to -38% (deaths: from 3,771 to -9,491) for age structure, and 33% to -15% (deaths: from 8,436 to -3,704) for population size (Fig 1C1). For method II, the contributions of the three factors changed respectively from 66% to -37% (deaths: from 16,668 to -9,370), 19% to -29% (deaths: from 4,833 to -7,404), and 15% to -33% (deaths: from 3,709 to -8,436) ( Fig 1C2).

A new decomposition method
When there are no interactions between the three factors, the two established methods generate the same decomposition results. However, in practice, interactions exist in almost all cases. To overcome the deficiencies of methods I and II, we propose a new approach, referred to as method III. This method relies on the principle that a reliable decomposition method should generate stable and consistent results regardless of the choice of reference population. When the reference population is switched, the signs of the contributions of the three factors are expected to reverse, but their magnitudes should remain the same.
Applying this principle, we suggest the three two-way interactions should be equally divided between relevant factors. As there is minimal theoretical guidance about how to allocate the three-way interaction to the three factors, we recommend dividing them equally. The detailed derivations are included in the S1 File. Method III uses four steps to decompose the contributions of the three factors. These steps are illustrated with population 1 as the reference, but results are the same with population 2 as the reference: i. Calculate the main effects of population size (M p ), age structure (M s ), and ASMR (M m ) using formulas (7)-(9).
ii. Calculate two-way interactions between population size and age structure (I ps ), age structure and ASMR (I sm ), and population size and ASMR (I pm ) using formulas (10)- (12).
iii. Calculate three-way interaction of the three factors (I psm ) using formula (13).
iv. Using AS III , ASMR III and PS III to represent deaths attributed to age structure, population size and ASMR by method III, calculate the contribution of each factor based on the formulas below.  Fig 1A1 (Fig 1A2) represents the decomposition results between the United States and China in 2017 from method I (method II); Fig 1B1 (Fig 1B2) represents the decomposition results between China in 1990 and 2017 from method I (method II); Fig 1C1 (Fig 1C2) represents the decomposition results between the United States in 1990 and 2017 from method I (method II). https://doi.org/10.1371/journal.pone.0216613.g001 The key to the preservation of magnitude and reversed sign is the equal allocation of the two-way interactions. For example, when population 2 is designated as the reference group, (s i2 -s i1 ) in formula (20) is replaced by (s i1 -s i2 ), leading to a reversed sign but a result with the same magnitude. As the three-way interaction also changes its sign but not magnitude when the reference group is switched, the principle holds no matter how we allocate the three-way interaction, i.e., even for unequal allocations. However, equal allocation is a reasonable choice when there is no prior information on the relative importance of the three factors.
Evaluating the performance of method III. Based on the same evaluation criteria mentioned above, we assessed the robustness of method III using the same three case studies. The results suggest that method III generates the same decomposition results regardless of the choice of reference.
Case study one. The gap in the total number of unintentional fall deaths between the United States and China in 2017 (-96,404) is explained -17% (-16,615) by ASMR, 38% (36,370) by age structure, and -120% (-116,159) by population size, using US as the reference (Fig 2A) (note: the sum of the percentages of three factors deviates slightly from 100% due to rounded numbers). Using China as the reference, these contribution estimates merely reverse their signs.
Case study two. Using 1990 as the reference, the gap in the total number of unintentional fall deaths between 1990 and 2017 in China (70,671) is explained 13% (9,183) by ASMR, 65% (46,032) by age structure, and 22% (15,456) by population size (Fig 2B). Robust attributions of three factors were observed when the reference was changed.
Case study three. Finally, using 1990 as the reference, the gap in the total number of unintentional fall deaths between 1990 and 2017 in the United States (25,210) is explained 51% (12,848) by ASMR, 26% (6,460) by age structure, and 23% (5,902) by population size (Fig 2C). Similar and robust attributions were observed when the reference was changed.

Discussion
This study presents evidence toward the two proposed research questions. First, we assessed two existing decomposition methods for quantifying the contributions of ASMR, age structure, and population size to the difference in total number of deaths between two populations. We found that both methods are sensitive to choice of the reference population. By mathematically formulating the decomposition procedures in analogy to a three-factor factorial experiment design, we found that the asymmetric allocations of interactions between relevant factors by both methods are responsible for the inconsistent decomposition results when the reference population is altered. Second, we proposed a new method that overcomes the limitations of the two existing methods.
The three case studies show that the new method III we propose generates different and robust attribution results compared to methods I and II. It overcomes potentially ambiguous or even conflicting interpretations caused by the choice of reference, and thus offers clear and consistent evidence to assist health policy-making.
Our new method (method III) is motivated by previous work by Das Gupta in decomposing a rate difference between two populations into contributions of structural factors (age structure for this study) and structure-specific rate (ASMR in our situation) [3,5]. However, there are notable differences between Das Gupta's method and ours. First, our decomposition method applies to differences in the total numbers of events of interest, while Das Gupta's method was designed to attribute difference in total rates. Second, Das Gupta's strategy is to average the contributions over all possible decomposition orders of the structural factors (the contribution of structure-specific rates is always decomposed first) and takes the mean as the contribution of the factors. In contrast, our method is built on intuitive and rational allocation of interactions to relevant factors and is a one-step approach that offers more efficient computation. Although modern computers can quickly enumerate all decomposition orders, taking average seems to be an ad hoc way to even out the potential inconsistency among the different decomposition orders. Our new method is theory-driven and totally irrelevant to decomposition order and therefore has robust and ubiquitous interpretation. Last, perhaps because of the complexity of Das Gupta's method, some previous publications present results based on one decomposition order rather than averaging over all possible decomposition orders [2,4], leading to conclusions that could be unstable in the presence of interactions. In contrast, our method III is easy to implement and less likely to be misused in practice.
We note that the equal allocation of two-way interactions in our method may not be optimal in aspects other than the robustness of results to the choice of reference population. The equal allocation of the three-way interaction is rather arbitrary, as the allocation proportion does not affect the final attribution. However, given that all two-way interactions are equally allocated to reach robustness, it seems logical to allocate the three-way interaction equally as well. In fact, if Das Gupta's approach is applied to decomposing rate differences between two groups, we found that its core idea is also to equally divide the interaction of two factors. This idea is similar to the strategy we recommend for attributing differences in total number of deaths between two groups to three factors.
As interactions almost always exist, we recommend use of method III for decomposition analysis, or use of method III in combination with traditional methods. For previous publications that adopted method I or method II to decompose death differences between populations [2,4,8,9], it would be informative to re-analyze those data using method III, or at least to perform sensitivity analyses by changing the reference population to ensure the results are robust.
Although we did not provide uncertainty estimates for the attribution results in our case studies, it is possible to quantify the uncertainty using standard statistical approaches. One possibility is to use a parametric formulation. For example, one could assume the sizes of age groups in each population follow a multinomial distribution with probability parameters p ij 's, for which s ij 's are estimates. Given the size of each age group, one could assume the number of death follows either a binomial distribution or a Poisson distribution with the group size as an offset [11,12]. Direct calculation of variances of the attribution results can be complex due to the presence of products of s ij and m ij ; however, parametric bootstrap is a straightforward alternative [13]. In addition, nonparametric bootstrap is another possibility, i.e., sampling individuals from each population with replacement and performing the decomposition analysis repeatedly [13].
We also note that Method III is readily generalizable to attribute differences in a broad spectrum of health outcomes other than number of deaths between two populations to similar factors (structure-specific rates, structural factors, and population size). For example, we could use this method to attribute differences in total incident cases from a specific infectious disease to population size, vaccination status (structural factor), and vaccine-specific incidence rates. In addition, the decomposition method for death rate is essentially the same to that for the number of deaths, the only difference exists in the number of factors to be attributed (two factors for rate difference and three factors for count difference). Therefore, our method is readily generalizable to decomposing difference in mortality rates.

Conclusion
When interactions between ASMR, age structure, and population size exist, the two existing decomposition methods generate inconsistent results that vary with the choice of reference population. The new method that we propose overcomes this limitation and is a compelling option to attribute a wide range of differences between two populations to contributing factors.
Supporting information S1 File. The details of formula derivation of three methods. (DOCX)