We used World Health Organization (WHO) Global Alert and Response (GAR) updates and performed systematic searches of three databases (PubMed, Scopus, and Google Scholar) to compile all confirmed and possible human cases of HPAI H5N1 virus infection. We included all articles published between January 1, 1997 and April 19, 2013 with keywords “H5N1,” “human,” and “humans.” We excluded articles that described non-human cases (animal or molecular studies), did not report individual case data, did not include data on laboratory-confirmed HPAI H5N1 cases, or described asymptomatic infections (e.g., seroprevalence studies).

We defined confirmed human H5N1 cases using the World Health Organization guidelines, requiring isolation of HPAI H5N1 virus, a positive result by reverse transcription polymerase chain reaction (RT-PCR) testing of clinical specimens using H5-specific primers and probes, an elevated H5-specific antibody titer of ≥1∶80 (or equivalent using the WHO protocol), or at least a fourfold rise in H5N1 virus neutralization antibody titer in paired sera [17]. We defined possible cases as those lacking laboratory confirmation but having symptoms and known contact with a confirmed human HPAI H5N1 case. This definition is more restrictive than the WHO definition of suspected and probable cases. The WHO definitions involve information on exposures, such as to raw poultry meat or environments contaminated with wild bird feces, that was not available to us in the literature. We then created Abstraction Form S1 based on clinical and demographic variables known to be relevant to influenza A (H1N1)pdm09 and useful in clinical practice.

We initially created a database of all cases published on the WHO Global Alert and Response (GAR) website, which includes only HPAI H5N1 cases reported from November 2003 to present. Although clinical laboratory data were not provided for each case, these were assumed to satisfy WHO reporting criteria. We then attempted to match all cases identified through literature sources to this database.

Two independent investigators (RP, MM) evaluated each article for inclusion; a third investigator (NK) resolved all disagreements. We included articles in all languages. A professional translator (YX) evaluated the numerous Chinese language articles. For Japanese, Russian, French, and Spanish languages, we verified inclusion with native-language speakers. For all other languages, we used PDF OCR X Community Edition for file conversion into text format (version 1.9.32, Burnaby, British Columbia) and a web-based translator for translation into English [18].

We then extracted the predefined set of variables for each case (Abstraction Form S1), systematically comparing demographic variables to avoid duplication. One co-author (TU) reviewed discrepancies between the WHO GAR publications and the published medical literature; based on his familiarity with several of these cases, we resolved discrepancies using data from the literature.

Using our pre-defined Abstraction Form S1, we extracted our primary outcome variable (mortality), demographic predictor variables (age, sex, country, per capita government expenditure on health [PCGEH], season, body mass index [BMI], and comorbidities), infection-related predictor variables (contact with poultry, delay from symptom onset to hospitalization, and whether the case was part of a cluster of known cases), and hospitalization predictor variables (laboratory data, pneumonia, acute respiratory distress syndrome [ARDS], and mechanical ventilation). We defined a case cluster as at least two geographically and temporally proximal (epidemiologically-linked) confirmed human HPAI H5N1 cases.

We obtained PCGEH at the average exchange rate (USD) for each country and year through the World Health Organization Global Health Observatory Data Repository [19]. Since data beyond 2011 are unavailable, data from 2011 were carried forward for later cases. Finally, we created a season variable based on month (Summer: June–August; Fall: September–November; Winter: December–February; Spring: March–May).

We performed all statistical analyses using R software (Version 3.0.0, Vienna, Austria) and defined statistical significance by an alpha level of 0.05. Our primary analytic goal was to develop a parsimonious decision tree model with optimal predictive ability for mortality following HPAI H5N1 virus infection. We first assessed bivariate associations between each predictor variable and mortality using logistic regression models. For continuous predictors (age, delay from hospitalization to symptom onset, and PCGEH), we visually assessed for the linearity assumption (Figure S1).

All continuous predictor variables had potentially nonlinear relationships with death, so we analyzed them both in continuous and in categorical form. We divided age into four categories similar to those associated with mortality from influenza A (H1N1)pdm09: 0–4 years, >4–18 years, >18–25 years, and >25 years [20]. PCGEH was split in quartiles. We created three groups of countries because most countries had very few cases; including these countries as individual predictors would cause over-fitting and statistical instability [21]. We therefore defined three categories of countries: Indonesia (n = 171), Egypt (n = 169), and all others combined (“Other”). Fewer than 100 cases were documented in each of the countries in the “Other” category.

As we performed initial analyses, we found that several parameters were missing data. We therefore developed a decision tree using CART methods [16]. CART models generally yield comparable results to logistic regression models, but unlike logistic regression, CART methods do not require missing data to be deleted or imputed (minimizing possible bias), capture higher-order interactions more easily, do not assume an underlying linear model, and generate a graphical prediction tool that is easy for practitioners to use and interpret [22], [23].

CART procedures build a decision tree by selecting locally optimal splits that minimize “impurity” on the outcome measure of the two child nodes. Low impurity on the outcome measure indicates that the classifier performs well at separating observations with one outcome (e.g., death) from observations with another outcome (e.g., survival). For example, if sex is a strong risk factor for mortality, then mortality will be similar within each sex and different between sexes. All possible binary splits are considered for both continuous and categorical variables. The initial split is chosen as the single best classifier on the outcome measure; then, within each child node, the splitting procedure is recursively repeated until no further splits are possible. All observations, including those with missing data, are included in model-building: at each split, the impurity index is simply calculated over only those observations not missing the relevant predictor variable. To avoid over-fitting, the initial large tree is pruned based on a cost-complexity index, which captures the tradeoff between better fit and added complexity due to each additional node in the tree [21]. The optimal final tree, defined as the tree with the lowest expected misclassification, is selected using cross-validation [16].

We assessed model performance using a receiver operating characteristic (ROC) curve and corresponding area under the curve (AUC) [23]. We constructed bootstrapped confidence intervals for the sensitivity and specificity thresholds of the ROC curve and for the AUC. To replicate the results of the CART model, we developed corresponding prognostic models using complete-case logistic regression (Table S1) and multiply-imputed logistic regression (Table S2).

Our search identified 3,227 potentially relevant articles published since 1997 (Figure 1). After removing 1,540 duplicate articles and 7 unavailable articles, two investigators (RP, MM) independently reviewed the remaining 1,680 articles. The review process yielded 163 articles meeting inclusion criteria, comprising 617 unique cases. Nearly all studies meeting inclusion criteria were case reports, so we did not assess methodological quality criteria.

The quality of data reporting on HPAI H5N1 cases was inconsistent. Several variables had missing data, most with homogeneity among the non-missing values; these were not included in the analysis (Figure 2). In other words, BMI, presence of comorbidities, whether the case was part of a cluster, laboratory data, presence of pneumonia, ARDS, and use of mechanical ventilation were mentioned almost exclusively for cases in which the parameter was present rather than absent. For example, only 10 case reports noted that a case had not received mechanical ventilation, 63 noted that a case had received ventilation, and the remaining 544 offered no information. In light of this limitation of the case-reporting process, we narrowed the scope of our prognostic model to include only demographic and infection variables, as post-admission hospital variables were usually unavailable. We removed these variables from analysis because we would have limited statistical power to assess their effects, and infrequently recorded homogenous variables may have little clinical utility (see Figure 2 for details on excluded variables). The eight variables used in analysis had 9% missing observations.

Demographic and clinical characteristics of the 617 cases are presented in Table 1. Overall mortality was 59%. Cases occurred in 15 countries, with 55% reported from Indonesia and Egypt combined (n = 171 and n = 169, respectively). 140 cases occurred within epidemiologically- and geographically-linked clusters, 13 occurred sporadically, and cluster status was unknown for the remaining 464. Cases tended to be young adults (median age 18 years; IQR 6–29). There was a marginal predominance of females (54% female, chi-square = 3.75, df = 1, p = 0.05). For all demographic parameters reported in the WHO GAR case summary (year, country, and mortality) [8], distributions for our literature-based dataset and the WHO summary were nearly identical (Table 1). We find little indication, therefore, of reporting or selection bias in our sample.

In bivariate logistic regression models, risk factors for mortality were longer delay to hospitalization, infection not in Egypt, older age, lower PCGEH, likely contact with poultry, female sex, and illness onset during summer months (Table 2).

The decision tree, trained on all 607 cases with observed mortality, evaluated seven candidate predictor variables: age, PCGEH, country group, delay from symptom onset to hospitalization, sex, contact with poultry, and season. The variables are listed here in descending order of variable importance in the decision tree, a measure based on split quality. The first four were used as splitting variables in the final, pruned tree (Figure 3).

The first node splits on age, with higher mortality in patients at least 4.5 years of age. In the second level of the tree, young children (<4.5 years) in high-PCGEH settings (≥32.65 USD) are predicted to survive, with the lowest mortality (4%) of all groups. Young children in low-PCGEH settings (<32.65 USD) are predicted to die (57% mortality).

Older patients (≥4.5 years) are further partitioned by country. Older cases in Indonesia are predicted to die, with the highest mortality (84%) of all groups. Older cases not in Indonesia are classified by one final split based on delay to hospitalization: cases with a short delay to hospitalization (<2.5 days) are predicted to survive (33% mortality), while cases hospitalized later after illness onset (≥2.5 days) are predicted to die (65% mortality).

We assessed the decision tree's performance with an ROC curve and corresponding AUC (Figure 4). Tested on all cases without missing data on mortality or any of the seven candidate model variables (n = 301), the AUC was 0.81 (95% CI: 0.76–0.86), indicating very good discrimination. Performance remained strong (AUC = 0.75; 95% CI: 0.71–0.78) when the model was instead tested on all cases with observed outcome but potentially missing any subset of the predictors (n = 607). The latter is a more difficult prediction task, requiring surrogate splitting based on the missing predictors. In surrogate splitting, the value of the missing predictor is estimated using the other predictors. This estimated value is then used as usual to classify the observation.

We also performed complete-case logistic regression and multiply-imputed logistic regression including all predictor variables that were candidates for inclusion in the CART model. All models yielded similar results (Tables S1 and S2).