In order to test the performance of the algorithms, we use four main datasets, containing information on a total of 26,766 deaths: three physician COD diagnosed VA datasets, namely the Indian Million Death Study (MDS) ²⁰ , South African Agincourt Demographic and Health Survey (DHS) dataset ²¹ , and Bangladeshi Matlab DHS dataset ²² , and one health facility diagnosed COD dataset, namely the PHMRC VA data collected from six sites in four countries (India, Mexico, the Philippines and Tanzania) ^{23, 24} . We use four combinations of the PHMRC data by age group (adult and child) and by site (all versus India-only); this filtering was done to determine the effect on results when deaths are collected from the same geographical setting. A total of seven datasets were used and are summarized in Table 1. These datasets are publicly available, except for the MDS, and have been used in other studies ^{11, 15, 24} .

The MDS VA dataset used in this study contains information on 12,225 child deaths from ages one to 59 months. For each death, two trained physicians independently and anonymously assigned a WHO ICD version 10 code ²⁵ . In the cases where the two physicians did not initially agree or reconcile on a COD, a third senior physician adjudicated ²⁰ . Similarly, the Agincourt dataset ²¹ underwent dual physician COD assignment on its 5,823 deaths for ages 15 to 64 years. COD assignment was slightly different for the Matlab dataset which has 2,000 deaths for ages 20 to 64 years; a single physician assigned a COD, followed by review and verification by a second physician or an experienced paramedic ²² . In contrast, the PHMRC dataset is comprised of 6,718 hospital deaths that were assigned a COD based on certain clinical diagnostic criteria, including laboratory, pathology, and medical imaging findings ^{23, 24} . For each VA datasets, we grouped the physician assigned CODs into 17 broad categories, refer to Table 2. We also show the distribution of records for each COD for each of the seven datasets used in our study.

An overview of our approach is shown in Figure 1. We transformed each VA dataset into binary format with VA survey questions being the attributes (columns), answers being the values of cells in rows (re-coded into binary format with ‘Yes’ coded as 1 and ‘No’ as 0) and CODs (group number identifier listed in Table 2) being the last (or the first) column. For all VA datasets, a death is represented as a row (record).

We divided each VA dataset into training and testing datasets. We trained multiple NBC models ¹⁷ on the transformed training datasets using the one-against-all approach ^{18, 19} . We choose NBC because it has shown better results on VA surveys in the past ¹¹ . The one-against-all approach was used because it improves the algorithm’s classification accuracy on datasets with several categories of dependent variables as demonstrated by past literature ^{18, 19} . This will be explained in detail in the next section. During testing, the trained NBC models assign CODs to each death in the testing dataset. The assigned causes are ordered by their probabilities with the assumption that top cause is most likely the real cause.

Training Naïve Bayes using one-against-all approach. NBC uses a training dataset to learn the probabilities of symptoms and their CODs ^{11, 17} . NBC first measures the probability of each COD, P(COD), in the training dataset. Secondly, it determines the conditional probabilities of each symptom given a particular COD, P(Sym|COD). Thirdly, NBC determines the probability of every COD given a VA record in the test set, i.e., P(COD|VA).

P(COD| VA) = P(COD) ∏ S y m ∈ V A P ( S y m | C O D )

Equation 1. Conditional probability of COD given a VA record.

P(COD|VA) is determined by taking the product of all P(Sym|COD) (i.e., all symptoms in the VA record) and P(COD). The highest P(COD|VA) value determines that COD as the correct COD. In particular, we chose the Naïve Bayes Multinomial classification algorithm that estimates probabilities by using a maximum likelihood estimate which is readily available in data mining software applications like Weka ^{17, 18} .

COD NBC = a r g m a x C O D ∈ C O D s P ( COD | VA )

Equation 2. Select the class with the maximum probability.

In the one-against-all approach, we built an NBC model for each COD instead of one model for all CODs. In this approach, a dataset with M categories of CODs (dependent variables) is decomposed into M datasets with binary categories (CODs). Each binary dataset Di has a COD Ci (where i = 1 to M) labelled as positive and all other CODs labelled as negative with no two datasets having the same CODs labelled as positive. Finally, NBC is trained on each dataset Di resulting in M Naïve Bayes models, as shown in Figure 2. Each model is then used to classify the CODs for records in the test dataset producing a probability of classification. The cause Ci (where i=1 to M) with the highest probability is considered as the correct classification.

Testing OAA-NBC on new surveys. During testing, each Naïve Bayes model predicts a COD for each VA record in the test dataset, resulting in a list of CODs for each VA record in the test dataset. The list of assigned CODs is sorted by the COD probabilities. We made a minor modification in the one-against-all approach; instead of selecting a COD with the highest probability, we ranked the CODs in descending order of their probabilities for each VA record. We kept the ranked probabilities to generate cumulative performance measures, which are described in detail in the next section.

A VA algorithm’s performance is measured by quantifying the similarity between the algorithm’s COD assignments to physician review (or clinical diagnoses in PHMRC) assignments. Since the community VA datasets included in this study come from countries that have weak civil and death registration, physician review is the most practical and relatively accurate (and only) option to use for assessing algorithm performance. Moreover, given that these deaths are unattended, it follows that there is no ‘gold standard’ for such community VA datasets. Nevertheless, we are confident in the robustness of dual physician review as initial physician agreement (i.e. where two physicians agreed right at the onset of COD coding) was relatively high; e.g., 79% for MDS and 77% for Agincourt.

We measured and compared the individual and population-level performance of all of the algorithms using the following metrics: sensitivity, partially chance corrected concordance (PCCC) and cause-specific mortality fraction (CSMF) accuracy. These measures are commonly used in VA studies ^{11, 15, 26} . They are shown in Equation 3 – Equation 5. They are helpful in objectively assessing the performance of VA algorithms, as they provide a robust strategy to assess an algorithm’s classification ability for test datasets with widely varying COD distributions ^{12, 26} .

S e n s i t i v i t y = T r u e p o s i t i v e T r u e p o s i t i v e + F a l s e ​ N e g a t i v e

Equation 3. Sensitivity of classification

P C C C ( k ) = S − k n 1 − k n

Where S = T r u e p o s i t i v e T r u e p o s i t i v e + F a l s e ​ N e g a t i v e

Equation 4. Partially chance corrected concordance (PCCC) of classification: S is the fraction of positively (correctly) assigned causes when the correct cause is in the top k assigned causes out of total n causes.

Sensitivity and PCCC are metrics that assess the performance of an algorithm for correctly classifying the CODs at the individual level. Sensitivity measures the proportion of death records that are correctly assigned for each COD ¹² . Similarly, PCCC computes how well a VA classification algorithm classifies the CODs at the individual-level while also taking chance (likelihood that it was randomly assigned a COD) into consideration ^{8, 11, 12, 15} .

C S M F A c c u r a c y = 1 − ∑ j = 1 n | C S M F j T r u e − C S M F j P r e d | 2 ( 1 − M i n i m u m ( C S M F j T r u e ) )

Where C S M F P r e d = ( T P + F P ) N a n d C S M F T r u e = ( T P + F N ) N

Equation 5. Cause-specific mortality fraction (CSMF) Accuracy of classification: n is the total COD and N is the total records.

In contrast, CSMF accuracy (hereafter referred to as ‘agreement’) is a measure for assessing how closely the algorithms classified the overall COD distribution at the population level ¹² . It can be observed from Equation 5 that CSMF Accuracy computes the absolute error between predicted COD distributions by an algorithm (pred) and the observed (true) COD distributions.

We measure the cumulative values of sensitivity, PCCC and agreement on each rank and for each algorithm; e.g., sensitivity at rank two represents the sensitivity of both rank one and rank two classifications, which facilitates measuring the overall performance of the algorithms for classifications at the top two or more ranks. If 60% of the individual classifications are correct at rank one and 15% more are correct at rank two then the overall accuracy is 75% at both ranks. Finally, we also perform a statistical test of significance on the results of all the datasets to ascertain that the difference in results is not by chance. This type of the statistical test depends on the data distribution and association between experiments. We use Wilcoxon signed rank test as we are unsure about normal data distribution of our results. Our null hypothesis is that there is no significant difference between OAA-NBC and another algorithm. This is further discussed in the results section.

In order to compare the performance between OAA-NBC, InterVA-4 ⁷ , Tariff ⁶ , InSilicoVA ⁸ and NBC ¹¹ , we follow a seven step procedure. In Step one, we partitioned each VA dataset using the commonly used evaluation criteria in data mining: 10-fold cross validation ¹⁸ . In 10-fold cross validation, a dataset is divided into 10 parts. Each part is created by using stratified sampling method—i.e., each part contains the same proportion of standardized CODs as the original dataset. In Step two, we selected one part for testing and nine parts for training from each VA dataset. In Step three, we trained OAA-NBC, InterVA-4, Tariff, InSilicoVA and NBC on the designated training data subsets from each partitioned VA dataset. In Step four, we generated classifications with ranks for each algorithm on the test part per VA dataset. In Step five, we calculated the cumulative sensitivity, PCCC and agreement for each rank per each VA dataset. In Step six, we repeated the process from Step two to Step five up to 10 repetitions with a different part for testing in each turn and for each VA dataset. In Step seven, we computed the average sensitivity, PCCC and agreement for each rank per VA dataset and algorithm.

We implemented OAA-NBC in Java and with Weka API ¹⁸ . Weka provides APIs for one-against-all approach and Naïve Bayes Multinomial classifier ¹⁸ . We used the OpenVA package version 1.0.2 in R to implement InterVA-4, Tariff, InSilicoVA and NBC algorithms. The data format also was transformed into InterVA-4 input format (Y for 1 and empty for 0 values). It is important to note that the Tariff version provided in the OpenVA package is computationally different from the IHME’s SmartVA-Analyze application tool. We used custom training option for InterVA-4 and InSilicoVA as present in OpenVA package in R. In custom training, the names of symptoms do not need to be in the WHO standardised format, and the rankings of the conditional probability P(symptom|cause) are determined by marching the same quantile distributions in the default InterVA P(symptom|cause). The reason for choosing customized training instead of using pre-trained global models is that different datasets have different proportions of symptoms and causes of deaths, and custom training allows algorithms to generate models customized for the dataset. It also allows for fair evaluation across algorithms, especially for the ones that only work by using customized training on datasets and acquire more knowledge of the dataset during testing.

We performed data partitioning, as discussed in Step 1, using Java and Weka’s ¹⁸ stratified sampling API. Each algorithm was executed on that partitioned data. We used a separate Java program to compute the cumulative measures of sensitivity, PCCC and agreement (CSMF accuracy) on the COD assignments of each algorithm for each VA dataset. This process ascertained that our evaluation measures are calculated in the exact same manner. Our source code for all the experiments is available on GitHub and is archived at Zenodo ²⁷ .

Figure 3 shows the averaged agreements (CSMF accuracy) by algorithms across all VA datasets using rank one (most likely) cause (COD) assignments and the fifth most likely cause assignments (rank five). OAA-NBC produces the highest agreements for most of the VA datasets, ranging from 86% to 90% for rank one; it comes second or identical to NBC for the PHMRC Child datasets (global and India). Furthermore, OAA-NBC agreements were relatively consistent across the VA datasets compared to some of the other algorithms that varied considerably, such as Tariff, InterVA-4 and InSilicoVA. As expected, the cumulative agreements increased the overall agreements for each algorithm when including the top five ranked classifications for every VA dataset.

Individual-level cumulative sensitivity results for classification ranks one and five are shown in Table 3; cumulative PCCC values were not shown as the values were very close to the cumulative sensitivity values. It can be observed from Table 3 that OAA-NBC has the highest sensitivity values for the first ranked (most likely) COD assignments compared to the other algorithms, ranging between 53-63%. When considering all top five ranked classifications, OAA-NBC has improved sensitivity values by 31–38%, with cumulative values ranging from 91–95%. In the case of Tariff, InterVA-4 and InSilicoVA, the sensitivity values are significantly lower (10–40%) in comparison to OAA-NBC; NBC does not differ substantially from OAA-NBC, as differences only range from 3–7%. These results show that OAA-NBC consistently yields closer agreement with the physician review or clinical diagnoses at the individual-level than the other algorithms on most of the VA datasets.

We also performed a Wilcoxon signed rank statistical test on the reported sensitivity in Table 3, generated from the five algorithms (we also included rank two to rank three values which are not shown in the table to minimize space). For 35 observations, the Wilcoxon signed ranked test yielded Z-score=5.194 and two tailed p-value=2.47 x 10 ^-7 between OAA-NBC and NBC. It yielded the same Z-scores and two tailed p-values against InSilicoVA, InterVA-4, and Tariff. Thus, there is a statistically significant difference between the sensitivity values generated by OAA-NBC and the four other algorithms (p < 0.05). Similarly, we conducted the Wilcoxon signed rank test on 35 observations of agreements for the five different algorithms, finding a statistically significant difference between OAA-NBC and the other algorithms (Z=4.248, p < 0.05).

Thus, the use of one-against-all approach with NBC (OAA-NBC) improves the performance of COD classification for VA records, and yields better COD assignments at the population- and individual-level, which are statistically different and not attributed to chance compared to the four other algorithms. This also conforms to the machine learning literature that the one-against-all approach improves the performance of classification algorithms when there are more than two classes (CODs) ¹⁹ . However, this does not indicate that OAA-NBC does not require improvement, as the overall sensitivity for the top ranked CODs per VA record is still lower than 80%. We also made an additional assessment on the COD sensitivity; Table 4 shows the sensitivity per cause for first ranked predictions and VA dataset for each algorithm (PHMRC Indian datasets are excluded as their results are similar to PHMRC global datasets and this minimizes space too). The sensitivity values vary per VA dataset and cause for all of the algorithms; road and transport injuries and other injuries were the only causes that OAA-NBC predicted consistently well for four out of the five VA datasets. However, there are several causes where the sensitivity of the classifications by OAA-NBC were lower than 50%, and in some cases, 0% (four causes in MDS and two causes in PHMRC – Child global datasets). Sensitivity values are 0% for COD groups that have proportion of records near 1% per VA dataset (number of records for each COD in VA datasets are shown in Table 2). In general, the algorithms performances vary on different CODs for certain conditions in VA datasets. For example, classifications were equal to or under 10% across all algorithms for HIV/AIDs, cancers, cardiovascular disease, and chronic respiratory diseases in the MDS dataset. Algorithms like OAA-NBC and NBC mostly have better sensitivity for COD groups that have higher proportion of records in training dataset. However, this is not always the case, and better sensitivity values also depend on how distinguishable VA records of a COD group are from all other COD groups. In the next section, we discuss the problem and effects of imbalance within datasets on the algorithms’ classification accuracy.