Influenza-associated excess mortality in the Philippines, 2006-2015.

Influenza-associated mortality has not been quantified in the Philippines. Here, we constructed multiple negative binomial regression models to estimate the overall and age-specific excess mortality rates (EMRs) associated with influenza in the Philippines from 2006 to 2015. The regression analyses used all-cause mortality as the dependent variable and meteorological controls, time, influenza A and B positivity rates (lagged for up to two time periods), and annual and semiannual cyclical seasonality controls as independent variables. The regression models closely matched observed all-cause mortality. Influenza was estimated to account for a mean of 5,347 excess deaths per year (1.1% of annual all-cause deaths) in the Philippines, most of which (67.1%) occurred in adults aged ≥60 years. Influenza A accounted for 85.7% of all estimated excess influenza deaths. The annual estimated influenza-attributable EMR was 5.09 (95% CI: 2.20-5.09) per 100,000 individuals. The EMR was highest for individuals aged ≥60 years (44.63 [95% CI: 4.51-44.69] per 100,000), second highest for children aged less than 5 years (2.14 [95% CI: 0.44-2.19] per 100,000), and lowest for individuals aged 10 to 19 years (0.48 [95% CI: 0.10-0.50] per 100,000). Estimated numbers of excess influenza-associated deaths were considerably higher than the numbers of influenza deaths registered nationally. Our results suggest that influenza causes considerable mortality in the Philippines-to an extent far greater than observed from national statistics-especially among older adults and young children.

Introduction

Influenza is a serious public health concern that causes 3−5 million cases of severe illness and about 290,000 to 650,000 deaths worldwide each year [1, 2]. Estimates of the influenza burden in individual countries are needed to formulate public health policies and strategies to control influenza. This is particularly important to protect those at greater risk of developing life-threatening influenza complications, such as young children, older adults, and people with chronic illnesses [2, 3]. However, numbers of influenza-attributable deaths are difficult to measure directly because influenza infections are not typically laboratory-confirmed and may not be diagnosed [4]. In addition, direct mortality measurements can miss deaths from secondary complications triggered by influenza infection (e.g., influenza-triggered exacerbation of pre-existing chronic illnesses).

Although the burden of influenza has been extensively evaluated in temperate regions of Europe and North America, it is less well characterized for many countries in Asia. The Philippines is located in Southeast Asia’s tropical climate region, which is generally considered an important source of new viruses and global influenza epidemics because of the large and highly interacting human and animal populations [5]. Influenza A and B viruses circulate throughout the year in the Philippines, and there are often multiple annual peaks in influenza activity [6-8]. Circulating influenza strains tend to match the Southern rather than the Northern Hemisphere vaccine strains, hence the Southern Hemisphere influenza vaccine is used in each year’s national vaccination program [8]. The burden of influenza in the Philippines is largely unknown. The mean annual influenza incidence rate has been estimated as 5.4 per 1,000 individuals in an urban region of the country, with particularly high incidence (22.6 per 1,000) in young children [6]. However, the rate of influenza-associated mortality has not been quantified. Here, we used negative binomial regression models to estimate the influenza-attributable excess mortality in the Philippines from 2006 to 2015, and compared these influenza mortality estimates with death registry data to quantify under-reporting of influenza deaths.

Methods

Study design

This was a retrospective analysis of influenza-associated deaths over the period Jan 1, 2006 to December 31, 2015 in the Philippines. The objective was to estimate the overall and age-specific excess mortality rates (EMRs) associated with influenza. Ethical approval was not required for this analysis of aggregated administrative data.

Data sources

Weekly all-cause deaths were obtained from the death registration dataset of the Philippine Statistics Authority [9] (data summarized in ). In the Philippines, deaths are certified using International Classification of Diseases Tenth Revision (ICD-10) codes [10]. However, official reporting systems rarely capture every death that occurred; for instance, one study found that in one Philippine province, only 77% of deaths were captured in government records [11]. But since the extent of this under-reporting problem has yet to be examined nationally, we made no adjustment for under-registration of deaths for this study.

Weekly percentages of laboratory-confirmed influenza A and B cases in the Philippines were obtained from the WHO’s Global Influenza Surveillance and Response System (GISRS) FluNet database [12] (). The GISRS data for the Philippines was collected through passive surveillance of influenza-like illness (ILI) and severe acute respiratory infection cases at sentinel sites located throughout the country. In the surveillance, ILI was defined as an acute respiratory infection with measured fever of ≥38C° and cough with onset within the last 10 days while severe acute respiratory infection (SARI) was defined as an acute respiratory infection with history of fever or measured fever of ≥38C° and cough with onset within the last 10 days and requires hospitalization [13]. ILI sentinel surveillance sites are health centers and hospital outpatient departments while SARI sentinel surveillance sites are hospital inpatient departments [8]. Laboratory confirmation of influenza virus from clinical samples was performed at the Research Institute for Tropical Medicine, Metropolitan Manila, Philippines (the WHO-designated National Influenza Center) by real-time reverse-transcription polymerase chain reaction. More details about the surveillance can be found in a previous study [8].

To account for climatic variation in influenza transmission and seasonality, meteorological data (rainfall, mean temperature, and relative humidity) were obtained from the 52 weather stations of the Philippine Atmospheric, Geophysical and Astronomical Services Administration [14]. Fifteen weather stations were excluded because they had more than 6 months of missing data. Weekly nationwide values for rainfall and relative humidity were obtained by averaging data from the 37 included weather stations. The average weekly temperature was calculated by taking the mean of the weekly maximum and weekly minimum, following the logic of the World Meteorological Organization’s recommended method of computing for average daily temperature [15].

Data analyses

Excess mortality associated with influenza was estimated for the overall population and for five age groups (0 to 4 y, 5 to 9 y, 10 to 19 y, 20 to 59 y, and ≥60 y) using negative binomial regression models (regression equations for each age group can be found in the, ). Negative binomial regression was used instead of the Poisson regression since deaths were over-dispersed. Because our study included the time when the Philippines was affected by typhoon Haiyan, a tropical storm resulting in >6,000 deaths, this event was formally included in the models which were based on the following equation: where t denotes time, E[Y] is the expected value of weekly number of all-cause deaths Y, and β values are the coefficients. β0 is the intercept; β1 to β6 account for the polynomial time trends; β7 and β8 are coefficients associated with the percentage of samples confirmed positive for influenza A and influenza B, respectively; β9 to β11 are coefficients for the meteorological data, rainfall (β9), mean temperature (β10), and relative humidity (β11); β12 and β13 pertain to the coefficients of the dichotomous variable for the Typhoon Haiyan week (week 45 of 2013) and the 2009 flu pandemic, respectively; β14 and β15 pertain to annual cyclical terms; β16 and β17 pertain to semiannual cyclical terms; and lastly, e is the error term that follows exp(e) ~ Gamma(1/α,α), and α is the overdispersion parameter. Annual and semiannual cyclical terms were used as seasonality controls since there are two seasons in the Philippines and influenza seasonality is known to be semi-annual, peaking from around June to November [8, 16]. The polynomial time trend and seasonality time trends were entered consecutively to determine the regression equation that best fitted the data. We also included various lags (no lag, one-week lag, and two-week lags) of the flu positivity rates to account for the possible delayed effect of flu on mortality. This iterative process () resulted in 84 regressions calculated for each age group. The regression that had the lowest Akaike Information Criterion [17], i.e. the model that provides the best fit was selected (). Because the number of tested samples was not available for week 16 of 2008, the value was imputed by taking the mean of the data from the week before and the week after. To estimate the excess mortality associated with influenza, we first calculated the annual predicted all-cause deaths for each age group using the chosen regression model, and then subtracted the annual predicted deaths without influenza A or without influenza B (coefficients for influenza A or B set to zero), as detailed elsewhere [1, 4, 18–20]. The resulting annual average mortality due to influenza A or B were divided by the 2015 total population and multiplied by 100,000 to get the EMR per 100,000 persons. The 95% confidence intervals (CI) for the EMRs were estimated through bootstrapping of residuals and re-estimating the excess mortality (1,500 iterations). The estimated number of influenza deaths was compared to the number of influenza deaths registered by the Philippine Statistics Authority with ICD-10 codes J10 (‘influenza due to other identified influenza virus’) and J11 (‘influenza due to unidentified influenza virus with other respiratory manifestations’). Sensitivity analyses were also performed to investigate the robustness of the estimates by running the same regressions with imputed values for all-cause deaths for week 45 of 2013 (coinciding with Typhoon Haiyan) and for the influenza positivity rate for weeks with ≤1 sample tested or with ≤10 samples tested. Imputed values were the average of the data from the week before and the week after the data point where possible; otherwise, they were the average of the two weeks nearest to the data point to be imputed. All analyses were performed using R version 3.5.1 [21].

Results

Descriptive analyses

An average of 485,412 all-cause deaths were registered per year in the Philippines between 2006 and 2015 (). Around 55% of all-cause deaths were of individuals aged ≥60 years (). A mean of 8,418 samples per year were tested for influenza during the same period, and a mean of 1,453 (17.3%) were positive for influenza virus (). Influenza A accounted for most (78.2%) of the confirmed influenza cases, although influenza B was detected more frequently than influenza A in 2008 and 2013. The proportion of samples testing positive for influenza A or B varied by year. One-third of samples (33.3%) were positive for influenza in the A/H1N1 2009 pandemic season, almost all of which were confirmed as influenza A. Excluding 2009, the proportion of samples positive for influenza A or B varied between 6.9% and 19.2%.

All-cause and influenza-associated mortality using the negative binomial regression method

The all-cause mortality predicted in the negative binomial regression replicated the registered weekly all-cause mortality for each age group (). In addition, the models fit the data well as they predicted a rise in the mortality rate from 2006 to 2015 among individuals aged 20 to 59 years and ≥60 years. There was a notable peak in registered deaths across age groups in week 45 of 2013 coinciding with Typhoon Haiyan.

Influenza was estimated to account for a mean of 5,347 excess deaths per year (about 1.1% of the average annual all-cause deaths for the study period) in the Philippines over the study period, most of which (67.1%) occurred in adults aged ≥60 years (). Influenza A accounted for 4,584 (85.7%) of all estimated excess influenza deaths. Although influenza B was estimated to cause fewer influenza deaths overall, it was responsible for about one-third (31.2%) of influenza deaths among 5−9-year-olds and one-fourth (24.1%) of influenza deaths among 0–4 year-olds.

The highest number of influenza-associated deaths was estimated to have occurred in 2009 (n = 8,784), almost all of which were caused by influenza A consistent with the 2009 A/H1N1 pandemic (). Deaths associated with influenza A or B were higher than average in the years 2013−2015 (n = 7,046−8,666) and were least frequent in 2008 (n = 1,807).

Overall, the annual EMR for influenza A and B-associated deaths was estimated as 5.09 (95% CI: 2.20–5.09) per 100,000 individuals (). The EMR was highest for individuals aged ≥60 years (44.63 [95% CI: 44.51–44.69] per 100,000), and second highest for children aged <5 years (2.14 [95% CI: 0.44–2.19] per 100,000). The lowest EMR was among individuals aged 10 to 19 years (0.48 [95% CI: 0.10–0.50] per 100,000). The overall annual EMR was 4.37 (95% CI: 1.49–4.37) per 100,000 for influenza A-associated deaths, and 0.73 (95% CI: –2.15–0.73) per 100,000 for influenza B-associated deaths.

A peak in influenza deaths also occurred in the week of Typhoon Haiyan (week 45 of 2013, figure not shown). However, in a sensitivity analysis, the EMRs were not considerably changed by using an imputed value for all-cause deaths at week 45 of 2013 (). The EMRs were also similar in other sensitivity analyses using imputed values for weeks with ≤1 sample tested (n = 1 week) and for weeks with ≤10 samples tested (n = 21 weeks).

Comparison of estimated influenza-associated deaths with national registry data

The estimated annual influenza-associated deaths were compared with influenza deaths registered by the Philippine Statistics Authority (). For all seasons, and across all age groups, the numbers of estimated excess influenza-associated deaths were considerably greater than those registered in the Philippines. Overall, compared to our estimated mean of 5,347 excess influenza-associated deaths per year, a mean of 96 deaths per year were registered.

Discussion

Our study provides the first estimate of influenza-associated mortality in the Philippines. The country’s influenza-associated EMR was considerable between 2006 and 2015: an estimated 5.09 influenza-attributable deaths occurred per 100,000 persons each year, and influenza was the cause of approximately one in every 100 deaths. A disproportionate percentage of influenza-attributable excess deaths (67.1%) occurred among individuals aged 60 years or older, considering that this age group represented just 7% of the total population [22].

The age group-specific EMRs estimated in our study are consistent with previous Filipino influenza mortality estimates by Iuliano et al., who estimated mortality rates per 100,000 of 4.0 (95% CI: 0.6−8.5) for those aged less than 65 years and 50.8 (95% CI: 12.9−96.1) for 65−74 year-olds [1]. However, these earlier estimates for the Philippines were extrapolated using EMRs of neighboring countries, and may not fully reflect the local parameters used in our study. On the other hand, our EMR for children aged <5 years of 2.1 per 100,000 is close to a recent meta-analyses by Wang et al. who found that EMR for the said age group for low-middle income countries is 1.7 [23]. The age-specific EMRs in our study are also aligned with those reported in other tropical countries where Southern Hemisphere influenza strains usually dominate [20, 24–26]. For instance, in Thailand, annual influenza-associated EMRs were estimated to be highest for those aged ≥65 years (42 per 100,000) [25] and in Western Kenya, the influenza EMR was estimated highest among those aged 50 years and older (74.0 per 100,000) and second highest in children aged less than 5 years (22.2 per 100,000) [26].

Our negative binomial regression was predictive of all-cause mortality in the Philippines, providing confidence in the modeling approach used. For example, we correctly predicted a considerable increase in influenza A deaths during 2009, coinciding with the A/H1N1pdm09 influenza pandemic that spread in most of Asia [8, 16]. Our regression models also projected a peak in influenza deaths in the week of typhoon Haiyan. Although this peak could be related to struggling health services and increased transmission of infections following the typhoon [27], it might also be partly influenced by the increased all-cause deaths at this time. By nature of regression analyses, trends of both dependent and independent variables are summarized by a fixed coefficient. Therefore, the peak in all-cause deaths–the dependent variable–could have translated to upward trends in estimated influenza deaths, even if the influenza positivity rates of tested samples were not elevated during the same period. Nonetheless, the overall EMRs were largely unaffected by this exceptional event since we controlled for the Typhoon Haiyan week in the main analyses and in the sensitivity analyses, where we used an imputed value to replace the peak in all-cause deaths for the Haiyan week. A similar effect from increasing all-cause deaths in some age groups may explain the upward trend in influenza deaths from 2013−2015.

For each age group, we estimated significantly more influenza deaths than the numbers registered in national statistics. Although this might partly follow unreliable use of ICD-10 during death registration and low level of diagnostic confirmation [28], we suspect that our greater estimates largely result from the additional deaths from underlying health conditions exacerbated by influenza (e.g., cardiovascular conditions) besides those caused by influenza directly (e.g., influenza that leads to pneumonia). Our findings suggest that influenza mortality in the Philippines is greater than previously thought, and this information may help encourage improvements in the national influenza surveillance and public health programs. Increasing influenza vaccination coverage among risk groups could be an effective way to reduce influenza-attributable mortality in the country [3]. Indeed, vaccination coverage was only 2.3% for adults aged 60 years and older during the Philippines’ last public influenza vaccination program [29].

Our regression analyses had several limitations. First, because of missing FluNet data and to minimize selection bias, we were unable to estimate deaths associated with individual influenza A subtypes or influenza B lineages. Second, few ILI cases in the FluNet database had laboratory confirmation of influenza virus, likely because virologic testing is expensive and of limited clinical value. This may have introduced some selection bias that we could not adjust for, and suggests that the influenza positivity rate data is not representative. Nonetheless, in our sensitivity analyses, the EMRs were largely unaffected by imputed influenza positivity rates for weeks with ≤1 or ≤10 specimens tested, implying that our results are robust despite the limited influenza activity data. Third, FluNet data may not be nationally-representative. Sentinel surveillance sites were only present in 13 out of 17 regions in the country [8] and there is no indication that sites were distributed evenly across the country. Despite its limitations, FluNet data is the only source of influenza activity data of its scope in the country to date. Fourth, because of missingness of the data, we had to drop 15 out of 52 weather stations for the meteorological variables. We believe that this is not an issue since the meteorological variables were merely controls, not the independent variable of interest. Finally, unlike other studies [4, 30], our analyses did not control for infections caused by respiratory syncytial virus because this data was not obtainable for the Philippines. Since respiratory syncytial viruses co-circulate with influenza viruses [31], our results might have overestimated influenza mortality, particularly in the younger age groups [32].

Conclusions

Our results suggest that the numbers of excess deaths attributable to influenza in the Philippines are considerably greater than those recorded in the national death registry, especially among older adults and young children. These findings underscore the importance of prioritizing older adults and children less than 5 years of age for influenza vaccination, in line with recommendations by the World Health Organization [3].

Negative binomial regression equations for each age group.

(DOCX)

Click here for additional data file.

Average weekly all-cause deaths and influenza-positive samples, 2006−2015.

(DOCX)

Click here for additional data file.

Model fitting algorithm.

(DOCX)

Click here for additional data file.

Selected negative binomial regression models.

(DOCX)

Click here for additional data file.

Sensitivity analyses.

(DOCX)

Click here for additional data file.

% Influenza A and influenza B to total samples tested.

(TIF)

Click here for additional data file.

Meteorological controls.

(TIF)

Click here for additional data file.

Actual all-cause mortality and predicted all-cause mortality with influenza set to zero.

(TIF)

Click here for additional data file.

16 Oct 2019

PONE-D-19-24478

Influenza-associated excess mortality in the Philippines, 2006-2015

PLOS ONE

Dear Mr. Cheng,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

Two reviewers agree that your manuscript needs a major revision, so please address all of their comments on methodological and statistical modelling issues before resubmitting.

We would appreciate receiving your revised manuscript by Nov 30 2019 11:59PM. When you are ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter.

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). This letter should be uploaded as separate file and labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. This file should be uploaded as separate file and labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. This file should be uploaded as separate file and labeled 'Manuscript'.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

We look forward to receiving your revised manuscript.

Kind regards,

Joël Mossong

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

http://www.journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and http://www.journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. Thank you for providing the following Funding Statement:"Funding: Sanofi Pasteur."

We note that one or more of the authors are employed by a commercial company: name of commercial company.

a) Please provide an amended Funding Statement declaring this commercial affiliation, as well as a statement regarding the Role of Funders in your study. If the funding organization did not play a role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript and only provided financial support in the form of authors' salaries and/or research materials, please review your statements relating to the author contributions, and ensure you have specifically and accurately indicated the role(s) that these authors had in your study. You can update author roles in the Author Contributions section of the online submission for

Please also include the following statement within your amended Funding Statement.

“The funder provided support in the form of salaries for authors [insert relevant initials], but did not have any additional role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. The specific roles of these authors are articulated in the ‘author contributions’ section.”

If your commercial affiliation did play a role in your study, please state and explain this role within your updated Funding Statement.

b). Please also provide an updated Competing Interests Statement declaring this commercial affiliation along with any other relevant declarations relating to employment, consultancy, patents, products in development, or marketed products, etc.

Within your Competing Interests Statement, please confirm that this commercial affiliation does not alter your adherence to all PLOS ONE policies on sharing data and materials by including the following statement: "This does not alter our adherence to  PLOS ONE policies on sharing data and materials.” (as detailed online in our guide for authors http://journals.plos.org/plosone/s/competing-interests) . If this adherence statement is not accurate and  there are restrictions on sharing of data and/or materials, please state these. Please note that we cannot proceed with consideration of your article until this information has been declared.

Please include both an updated Funding Statement and Competing Interests Statement in your cover letter. We will change the online submission form on your behalf.

Please know it is PLOS ONE policy for corresponding authors to declare, on behalf of all authors, all potential competing interests for the purposes of transparency. PLOS defines a competing interest as anything that interferes with, or could reasonably be perceived as interfering with, the full and objective presentation, peer review, editorial decision-making, or publication of research or non-research articles submitted to one of the journals. Competing interests can be financial or non-financial, professional, or personal. Competing interests can arise in relationship to an organization or another person. Please follow this link to our website for more details on competing interests: http://journals.plos.org/plosone/s/competing-interests

3. We note that you have indicated that data from this study are available upon request. PLOS only allows data to be available upon request if there are legal or ethical restrictions on sharing data publicly. For information on unacceptable data access restrictions, please see http://journals.plos.org/plosone/s/data-availability#loc-unacceptable-data-access-restrictions.

In your revised cover letter, please address the following prompts:

a) If there are ethical or legal restrictions on sharing a de-identified data set, please explain them in detail (e.g., data contain potentially identifying or sensitive patient information) and who has imposed them (e.g., an ethics committee). Please also provide contact information for a data access committee, ethics committee, or other institutional body to which data requests may be sent.

b) If there are no restrictions, please upload the minimal anonymized data set necessary to replicate your study findings as either Supporting Information files or to a stable, public repository and provide us with the relevant URLs, DOIs, or accession numbers. Please see http://www.bmj.com/content/340/bmj.c181.long for guidelines on how to de-identify and prepare clinical data for publication. For a list of acceptable repositories, please see http://journals.plos.org/plosone/s/data-availability#loc-recommended-repositories.

We will update your Data Availability statement on your behalf to reflect the information you provide.

Additional Editor Comments (if provided):

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: No

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: No

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: The authors estimated influenza-associated mortality in the Philippines, which added one more component to understand the global burden of influenza. In general, this paper was well written, and I have some comments below.

1. What contributes to the sharp peak in death in 2008 for children aged 5-9 years?

2. The authors could consider to add a dummy variable in the regression model to adjust for the impact of Typhoon Haiyan, like what has been done to adjust for SARS outbreak.

3. Why would the effects of temperature and relative humidity be opposite in different age groups?

4. The mathematical formula is incorrect. In the negative binomial model, there should be a log link between dependent and independent variables.

5. In table 1, why there is a substantial decrease of samples tested in 2014 and 2015. Is there any change of influenza surveillance system in Philippines?

6. The axes of the figures could be revised to make the figures more readable.

Reviewer #2: This manuscript estimates Influenza-associated excess mortality in the Philippines during the period 2006 through 2015 using negative binomial regression models. As the authors indicate, mortality burden of influenza in the Philippines has not been quantified before thus findings from this manuscript will help to inform public health policies and strategies for the control of influenza.

Here are some comments that the authors could consider to help the readers to better understand their study and the findings that they present.

1. In line 78: The authors state, “No adjustment was made for under-registration of deaths”. Before that, the authors indicate that deaths “must be registered within 48 hours”. Isn’t this a contradiction? It would be of help to the readers if the authors could state why no adjustments were made for under-registration. Is it that no data are available to quantify deaths that are not registered, or perhaps that there is so much variation by site/hospital or region in the country to allow for meaningful adjustments?

2. The authors indicate that GISRS data were collected through ILI and SARI cases at sentinel sites. Can they include the case definitions used for ILI and SARI or provide appropriate references for the readers?

3. I note that the authors used a weekly time-series to model the overall and age-specific EMRs associated with influenza. Did the authors consider using time-lagged independent variables (particularly for meteorological variables and influenza activity)? This is not stated in the text and one would expect a time-lagged effect of these variables on mortality.

4. Also, is there a particular reason why the meteorological variables (rainfall, temperature, and humidity) were included in the model to estimate EMRS? Has there been data to suggest that these are important for the Philippines? Did inclusion of these variables result in better fitting models, across the age groups assessed?

5. In lines 118-119, the authors state the “Other missing data were not replaced in the principal analyses”. This is somehow ambiguous. Could the authors state what these data are?

6. Reading through the methods, I assume that the data from the 2009 pandemic period were included in the analyses, is this correct? If so, I would suggest that the authors rerun the analyses with the pandemic data excluded. This is particularly important if they are seeking to estimate the mean annual EMRs associated with influenza. A quick look at the data suggests an additional few hundred deaths (~300 deaths) annually if the data from the pandemic period are included.

7. Related to that, I suggest that the authors only model the estimates using the imputed deaths for week 45 in 2013 (when the typhoon Haiyan occurred). Regardless of the fact that they conducted sensitivity analyses, it is clear that the spike in deaths during that week was out of the ordinary and thus it would only make sense if it were excluded.

8. The authors mention that the 95% CIs were estimated using bootstrapping methods with 1,500 iterations. However, I note that the CIs are very narrow, particularly for data among young children <5 years. Could the authors look again at this, and perhaps comment about it in the discussion?

9. In lines 208-209, the authors state “…. and influenza was the cause of approximately one in every 100 deaths”. This is not explicit in the results section but I suppose that this is based on dividing the total mean annual all-cause deaths by the mean annual influenza-associated deaths. Could the authors try to make this more explicit in the text? This further highlights why you should not include the 2009 pandemic data in your EMRs models.

Reviewer #3: The authors intend to quantify mortality associated to influenza in the Philippines, which is very important for public health and prevention. Hence, this is a very important paper.

Some general comments:

The authors use three metrological measures, rain, temperature and humidity, which each also have seasonal variation. Therefore, the inclusion of the yearly and half-yearly sines must have been included to adjust for residual seasonality not covered by these metrological variables.

Argue why a negative binomial regression (compensate for over-dispersion – what about under-dispersion?)

Additive or multiplicative model? – link function

Selection of elements to be included in the model is based on having positive A and B coefficients – argue why - and secondly the lowest AIC. However, in the S2 Table I miss this information.

Figure 2 show huge peaks in number of deaths associated to influenza in week 45 2013 (The typhoon Haiyan), This I don’t understand.

- Was there a huge peak in positive influenza samples that week? – I do not believe so, probably none or very few samples were sampled in that week.

- As there are no peaks in the model (red line) in figure 1, and the number of deaths associated to influenza was calculated as the prediction from the full model minus the prediction by the same model, but with A and B set to zero. There should not be calculated peaks.

This indicate that the calculation of influenza-associated number of deaths is wrong!

Would be nice to show both the full models and the models with A and/or B as zero in graph 1.

Suggest including graphs showing the A and B positive percentages over calendar time used in the model and as supplementary the metrological parameters.

The authors intend to compare all-cause influenza-associated mortality with cause-specific influenza-associated mortality. The most commonly used cause-specific mortality is respiratory (ICD10 …) e.g. references 1 and 2. The authors only look at J10 and J11 (influenza the main cause of deaths). This is of cause interesting, but it would have been of more interest, if they (also) had looked at respiratory cause of deaths, and made serious comparisons, for example if there is a general relation e.g. has it been suggested that influenza-associated mortality estimated using all-cause is the double of respiratory.

Alternatively, the authors could leave out cause-specific influenza-associated mortality, and write another article comparing all-cause and cause-specific.

Suggest to include a typhoon Haiyan parameter in the model: 1 in week 45 2013, else 0.

The Philippines consist of many islands with varying population. The metrological stations are properly distributed more-or-less evenly over the whole area, why average metrological measures should be population weighted. Likewise for the influenza data.

Is this possible? – if not, this should be discussed as a limitation.

Minor comments:

Page 3, line 64: I believe ‘mortality has not been accurately quantified’ attribute too much faith in the model. Suggest to leave out ‘accurately’.

Page 4, line 66-67: ‘… and compare …’. Suggest ‘… and compare all-cause estimates with influence cause-specific estimates’.

Page 4, line 78: Is it correct that all deaths in the Philippines are registered with cause of deaths within 48 hours? – Faster than in any other country, I know of.

Page 4, line 81-86: What do you mean by and what is the difference between ILI sentinel sites and SARI (Severe acute respiratory infections) sites? – how many of each?

Page 5, line 93. I cannot find ‘average weekly temperature’ in reference 13, only average daily min and max temperatures.

Page 5, line 105: Ok to use a polynomial spline, but why 6?

Page 6, line 118: ‘Other missing data …’ - how many?

Page 13, line 212-215. I believe it is highly surprising that this studies all-cause influenza-associated mortality estimates are consistent with the respiratory-cause-specific estimates from Iuliano et al. study

Page 13, line 223: I would not use the word ‘accurately’, but something like our model fitted data well

Page 13, line 223-4: Sometimes ‘official deaths registry statistics’ stand for cause-specific and here for all-cause. Suggest using all-cause and cause-specific.

Page 13, line 224: It is not correct that your model predicted the typhoon Haiyan peak! – see figure 1

Page 14-15, line 252-255: You might have used, what is often called the Goldstein index: ILI-rate * positive-percentage, where the ILI rate reflect the population dynamics and the positive-percentage limit the ILI-rate to influenza i.e. exclude other circulating respiratory pathogens.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Xiling Wang

Reviewer #2: No

Reviewer #3: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files to be viewed.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. Please note that Supporting Information files do not need this step.

Submitted filename: PONE-D-1924478-Review Comments.pdf

Click here for additional data file.

Submitted filename: Phillippines review.docx

Click here for additional data file.

9 May 2020

Point-by-point response to Reviewers’ comments

Reviewer #1:

The authors estimated influenza-associated mortality in the Philippines, which added one more component to understand the global burden of influenza. In general, this paper was well written, and I have some comments below.

1. What contributes to the sharp peak in death in 2008 for children aged 5-9 years?

Upon reviewing our files, we found out that this sharp peak was due to a typographical error that arose when we transferred our estimates to Excel to generate the graphs. Apologies for this inconvenience and thank you for noticing. Please rest assured that all estimates were done appropriately – it was just in the graphing that the error arose.

2. The authors could consider to add a dummy variable in the regression model to adjust for the impact of Typhoon Haiyan, like what has been done to adjust for SARS outbreak.

Thank you very much for your suggestion. We revised the model to include a dummy for Haiyan and a dummy for the pandemic year. Our new modelling algorithm is documented in S1 Text and S2 Table.

3. Why would the effects of temperature and relative humidity be opposite in different age groups?

While these effects are interesting and could be investigated, we feel that this could be explored and further discussed in a different study. We are primarily interested in flu activity as predictors of deaths and these meteorologic controls were included as previous studies. We feel that the changes in direction are secondary to the goal of improving model fit and we have shown that the fits of our models that include these parameters fit the data well.

4. The mathematical formula is incorrect. In the negative binomial model, there should be a log link between dependent and independent variables.

Thank you for letting us know. We rewrote it as:

E[Y_t ]=exp{β_0+β_1 t+β_2 t^2+β_3 t^3+β_4 t^4+β_5 t^5+β_6 t^6+ β_7 [Influenza A]_t+β_8 [Influenza B]_t + β_9 [Rainfall]_t+β_10 [Mean Temperature]_t+β_11 [Relative Humidity]_t+β_12 [Haiyan]_t+β_13 [Pandemic]_t+β_14 [sin⁡(2πt/52) ]+β_15 [cos⁡(2πt/52) ]+β_16 [sin⁡(2πt/26) ]+β_17 [cos⁡(2πt/26) ]+e_t }

where e is the error term that follows exp(e_t )~Gamma(1⁄α,α), and α is the overdispersion parameter following past research.

5. In table 1, why there is a substantial decrease of samples tested in 2014 and 2015. Is there any change of influenza surveillance system in Philippines?

According to the Research Institute of Tropical Medicine, the National Influenza Center, the data collection was scaled down from 2012 onward. However, this should not be an issue since we used positivity rates instead of raw counts.

Source: http://ritm.gov.ph/reference-laboratories/national-reference-laboratories/influenza-and-other-respiratory-viruses/

6. The axes of the figures could be revised to make the figures more readable.

We have done so.

Reviewer #2:

This manuscript estimates Influenza-associated excess mortality in the Philippines during the period 2006 through 2015 using negative binomial regression models. As the authors indicate, mortality burden of influenza in the Philippines has not been quantified before thus findings from this manuscript will help to inform public health policies and strategies for the control of influenza.

Here are some comments that the authors could consider to help the readers to better understand their study and the findings that they present.

1. In line 78: The authors state, “No adjustment was made for under-registration of deaths”. Before that, the authors indicate that deaths “must be registered within 48 hours”. Isn’t this a contradiction? It would be of help to the readers if the authors could state why no adjustments were made for under-registration. Is it that no data are available to quantify deaths that are not registered, or perhaps that there is so much variation by site/hospital or region in the country to allow for meaningful adjustments?

Apologies for the confusion. We meant that the policy dictates that deaths should be recorded within 48 hours. But this does not always happen. In fact, there is one small study that quantifies the extent of underreporting. This was done in only one province in the Philippines so we cannot confidently say that the extent is true nationwide. Given the uncertainty of the degree of underreporting according to location and according to cause of disease, we have decided to not adjust for underreporting. Our findings would likely underestimate the burden of disease due to flu. We have included these points in the limitations section of our discussion. Moreover, since the analysis does not depend on 48 hour policy, we decided to remove that information to avoid confusing the readers.

Source: Carter, K. L., Williams, G., Tallo, V., Sanvictores, D., Madera, H., & Riley, I. (2011). Capture-recapture analysis of all-cause mortality data in Bohol, Philippines. Population health metrics, 9(1), 9.

2. The authors indicate that GISRS data were collected through ILI and SARI cases at sentinel sites. Can they include the case definitions used for ILI and SARI or provide appropriate references for the readers?

We have done as you suggested in lines 88-94: “In the surveillance, ILI was defined as an acute respiratory infection with measured fever of ≥38C° and cough with onset within the last 10 days while severe acute respiratory infection (SARI) was defined as an acute respiratory infection with history of fever or measured fever of ≥38C° and cough with onset within the last 10 days and requires hospitalization.” Our definitions are aligned with the definitions by FluNet and WHO.

3. I note that the authors used a weekly time-series to model the overall and age-specific EMRs associated with influenza. Did the authors consider using time-lagged independent variables (particularly for meteorological variables and influenza activity)? This is not stated in the text and one would expect a time-lagged effect of these variables on mortality.

Thank you for your suggestion. We tried introducing up to two-period lags for the flu variables and the model with two lags yielded the best model fit. We did not introduce lags for the meteorological controls since this is not commonly done in similar studies (e.g. Aungkalanon et al 2015; Chow et al 2006).

Source: Aungkulanon, S., Cheng, P. Y., Kusreesakul, K., Bundhamcharoen, K., Chittaganpitch, M., Margaret, M., & Olsen, S. (2015). Influenza‐associated mortality in T hailand, 2 006–2011. Influenza and other respiratory viruses, 9(6), 298-304.

Chow, A., Ma, S., Ling, A. E., & Chew, S. K. (2006). Influenza-associated deaths in tropical Singapore. Emerging infectious diseases, 12(1), 114.

4. Also, is there a particular reason why the meteorological variables (rainfall, temperature, and humidity) were included in the model to estimate EMRS? Has there been data to suggest that these are important for the Philippines? Did inclusion of these variables result in better fitting models, across the age groups assessed?

Many previous studies (e.g. Aungkulanon et al 2015; Chow et al 2006) included meteorological controls in their models. A local paper by Lucero et al (2016) points out that flu circulates year-round, but is more pronounced from June to November. Given that the extant literature suggests that meteorological controls are important in modelling EMR, we included meteorological controls in all our models by default.

Source: Aungkulanon, S., Cheng, P. Y., Kusreesakul, K., Bundhamcharoen, K., Chittaganpitch, M., Margaret, M., & Olsen, S. (2015). Influenza‐associated mortality in T hailand, 2 006–2011. Influenza and other respiratory viruses, 9(6), 298-304.

Chow, A., Ma, S., Ling, A. E., & Chew, S. K. (2006). Influenza-associated deaths in tropical Singapore. Emerging infectious diseases, 12(1), 114.

Lucero, M. G., Inobaya, M. T., Nillos, L. T., Tan, A. G., Arguelles, V. L. F., Dureza, C. J. C., ... & Rodriguez, T. (2016). National Influenza Surveillance in the Philippines from 2006 to 2012: seasonality and circulating strains. BMC infectious diseases, 16(1), 762.

5. In lines 118-119, the authors state the “Other missing data were not replaced in the principal analyses”. This is somehow ambiguous. Could the authors state what these data are?

Apologies for the confusion. There is no other missing data apart from the positivity rate of 2008 week 16. Therefore, we deleted this sentence.

6. Reading through the methods, I assume that the data from the 2009 pandemic period were included in the analyses, is this correct? If so, I would suggest that the authors rerun the analyses with the pandemic data excluded. This is particularly important if they are seeking to estimate the mean annual EMRs associated with influenza. A quick look at the data suggests an additional few hundred deaths (~300 deaths) annually if the data from the pandemic period are included.

Thank you for this excellent point. We tried running the models without the pandemic year but the results seem to be close to the original EMR of 5.33 per 100,000. In fact, we got a slightly higher EMR of 5.81 when we excluded the pandemic year. Thus, we decided to just keep the whole dataset and run the models with the pandemic year controlled for. After controlling for the pandemic, our new runs yielded a more conservative estimate of EMR of 5.03.

7. Related to that, I suggest that the authors only model the estimates using the imputed deaths for week 45 in 2013 (when the typhoon Haiyan occurred). Regardless of the fact that they conducted sensitivity analyses, it is clear that the spike in deaths during that week was out of the ordinary and thus it would only make sense if it were excluded.

Thank you for pointing this out. We agree that Haiyan caused an unusual spike in death for week 45 of 2013. However, we felt that it would be better to keep the data as it is for the base case, and then run the modelling on imputed values. We dealt with it by adding a dummy variable for Haiyan as a control variable instead. We conducted a sensitivity analyses that either excludes this week or used an imputed value instead and found that the results did not change to a large degree.

8. The authors mention that the 95% CIs were estimated using bootstrapping methods with 1,500 iterations. However, I note that the CIs are very narrow, particularly for data among young children <5 years. Could the authors look again at this, and perhaps comment about it in the discussion?

Save for the number of iterations, we have little control over the resulting CIs since bootstrapping is entirely data driven. We have double-checked our code and re-ran the analysis to ensure correctness of analysis procedure.

9. In lines 208-209, the authors state “…. and influenza was the cause of approximately one in every 100 deaths”. This is not explicit in the results section but I suppose that this is based on dividing the total mean annual all-cause deaths by the mean annual influenza-associated deaths. Could the authors try to make this more explicit in the text? This further highlights why you should not include the 2009 pandemic data in your EMRs models.

This was calculated by dividing the average excess mortality per year by the annual average all cause mortality. We made it more explicit on lines 185-186: “Influenza was estimated to account for a mean of 5,347 excess deaths per year (about 1.1% of the average annual all-cause deaths for the study period)…”

Reviewer #3:

The authors intend to quantify mortality associated to influenza in the Philippines, which is very important for public health and prevention. Hence, this is a very important paper.

Some general comments:

1. The authors use three metrological measures, rain, temperature and humidity, which each also have seasonal variation. Therefore, the inclusion of the yearly and half-yearly sines must have been included to adjust for residual seasonality not covered by these metrological variables.

We agree. Therefore, we added the sine and cosine terms.

2. Argue why a negative binomial regression (compensate for over-dispersion – what about under-dispersion?) Additive or multiplicative model? – link function

We also considered using Poisson models in the analysis but results from using negative binomial models showed evidence of over-dispersion (alpha is not equal to zero and the point estimate is greater than 1). We have added these details in the methods.

Line 112: “Negative binomial regression was used instead of the Poisson regression since deaths were over-dispersed.”

As is usual practice, we used a log-link in the models. We apologize for failing to include this detail. We have updated the formula to be more explicit about the link function used.

3. Selection of elements to be included in the model is based on having positive A and B coefficients – argue why - and secondly the lowest AIC. However, in the S2 Table I miss this information.

We agree with this comment and have decided to just use the model fit (AIC) as the sole criterion to select the model. We have repeated the model selection and have presented updated results.

In addition, based on the comments of the other reviewers, we decided to introduce up to two lags of the flu variables and control for the Haiyan week and the 2009 pandemic. To facilitate readability of the paper, we opted to show just the results for the selected model.

4. Figure 2 show huge peaks in number of deaths associated to influenza in week 45 2013 (The typhoon Haiyan), This I don’t understand. - Was there a huge peak in positive influenza samples that week? – I do not believe so, probably none or very few samples were sampled in that week.

We agree that this appears as a surprising result. We feel that this has been addressed by the results of the sensitivity analysis and have been mentioned in the Discussion:

Results, Lines 209-211:

“A peak in influenza deaths also occurred in the week of Typhoon Haiyan (week 45 of 2013, figure not shown). However, in a sensitivity analysis, the EMRs were not considerably changed by using an imputed value for all-cause deaths at week 45 of 2013 (S4 Table).”

Discussion, Lines 250-260:

“Our regression models also projected a peak in influenza deaths in the week of typhoon Haiyan. Although this peak could be related to struggling health services and increased transmission of infections following the typhoon [27], it might also be partly influenced by the increased all-cause deaths at this time. By nature of regression analyses, trends of both dependent and independent variables are summarized by a fixed coefficient. Therefore, the peak in all-cause deaths – the dependent variable – could have translated to upward trends in estimated influenza deaths, even if the influenza positivity rates of tested samples were not elevated during the same period. Nonetheless, the overall EMRs were largely unaffected by this exceptional event since we controlled for the Typhoon Haiyan week in the main analyses and in the sensitivity analyses, where we used an imputed value to replace the peak in all-cause deaths for the Haiyan week.”

5. As there are no peaks in the model (red line) in figure 1, and the number of deaths associated to influenza was calculated as the prediction from the full model minus the prediction by the same model, but with A and B set to zero. There should not be calculated peaks.

This indicate that the calculation of influenza-associated number of deaths is wrong!

While there was a peak in actual all-cause mortality in 2014 due to Haiyan, we do not think that models are made to predict all the peaks and troughs as exactly as the real data. We would run the risk of overfitting the model if we try to predict every single peak and trough. However, since Haiyan is a special event that needs to be accounted for, we have revised the modelling strategy to control for the Haiyan week. We found that our updated results did not lead to drastically different estimates. We also found that the original models without the Haiyan and pandemic controls and the revised models perform well as its predicted all-cause mortality is 90+% correlated with the actual all-cause mortality. For the final paper, we’ve decided to include the Haiyan variable in the main model.

6. Would be nice to show both the full models and the models with A and/or B as zero in graph 1.

We feel that adding this detail in the main text might be too much but we have included these graphs as part of the supplemental material.

7. Suggest including graphs showing the A and B positive percentages over calendar time used in the model and as supplementary the metrological parameters.

We have done so as a supplementary file. See S1 Figure.

8. The authors intend to compare all-cause influenza-associated mortality with cause-specific influenza-associated mortality. The most commonly used cause-specific mortality is respiratory (ICD10 …) e.g. references 1 and 2. The authors only look at J10 and J11 (influenza the main cause of deaths). This is of cause interesting, but it would have been of more interest, if they (also) had looked at respiratory cause of deaths, and made serious comparisons, for example if there is a general relation e.g. has it been suggested that influenza-associated mortality estimated using all-cause is the double of respiratory.

Alternatively, the authors could leave out cause-specific influenza-associated mortality, and write another article comparing all-cause and cause-specific.

We agree with this point since other papers have used cause-specific causes of death in the modelling. There is however the problem of inaccurate coding of causes of death. Many doctors were not cognizant to the right way of coding cause of death on the death certificates (Hufana et al 2009). Due to this uncertainty in coding, we went with the conservative route of just using all-cause mortality instead of the cause-specific deaths because there is more certainty in the counts of total deaths than cause-specific deaths.

Source: Carter, K.L., Williams, G., Tallo, V. et al. Capture-recapture analysis of all-cause mortality data in Bohol, Philippines. Popul Health Metrics 9, 9 (2011) doi:10.1186/1478-7954-9-9

Hufana L, Cajita J, Morante L, Lopez J, Tan CL, Mikkelsen L, et al. Assessing the production, quality and use of national vital statistics: a case study of the Philippines. 2009. [cited 7 Mar 2019]. Available from: https://crvsgateway.info/file/7572/98

9. Suggest to include a typhoon Haiyan parameter in the model: 1 in week 45 2013, else 0.

We agree with this suggestion and have revised our model.

10. The Philippines consist of many islands with varying population. The metrological stations are properly distributed more-or-less evenly over the whole area, why average metrological measures should be population weighted. Likewise for the influenza data.

Is this possible? – if not, this should be discussed as a limitation.

We had to drop weather stations that had 6-months’ worth of missing data so in effect, our meteorological averages were computed from 37 instead of 52 weather stations. We therefore cannot conclusively say that the weather stations we ended up using are proportionate to the population. Furthermore, using a weighted mean calls for another layer of assumption, specifically, the weight to be used. Because of these reasons, we decided to use the simple mean instead.

Similarly, the FluNet does not provide the information on a sentinel-site basis. In addition, there is no indication that the sentinel sites are distributed evenly throughout the country. Lucero et al (2016) states that the sites were only present in 13 out of 17 regions in the country. Thus, we cannot weigh the flu positivity rates by population size.

Following the previous suggestions, we added the following sentences to specify that we had limitations with the weather station and flu surveillance data.

Lines 283-290: “Third, FluNet data may not be nationally-representative. Sentinel surveillance sites were only present in 13 out of 17 regions in the country [8] and there is no indication that sites were distributed evenly across the country. Despite its limitations, FluNet data is the only source of influenza activity data of its scope in the country to date. Fourth, because of missingness of the data, we had to drop 15 out of 52 weather stations for the meteorological variables. We believe that this is not an issue since the meteorological variables were merely controls, not the independent variable of interest.”

Lucero, M.G., Inobaya, M.T., Nillos, L.T. et al. National Influenza Surveillance in the Philippines from 2006 to 2012: seasonality and circulating strains. BMC Infect Dis 16, 762 (2016) doi:10.1186/s12879-016-2087-9

Minor comments:

11. Page 3, line 64: I believe ‘mortality has not been accurately quantified’ attribute too much faith in the model. Suggest to leave out ‘accurately’.

We deleted the word “accurately.”

12. Page 4, line 66-67: ‘… and compare …’. Suggest ‘… and compare all-cause estimates with influence cause-specific estimates’.

Apologies for being unclear. What we meant to say was that we wanted to compare the flu mortality estimates with the official flu mortality data from the death registry data. So that it will be clearer, we revised the sentence into (lines 66-69): “Here, we used negative binomial regression models to estimate the influenza-attributable mortality in the Philippines from 2006 to 2015, and compared the flu mortality estimates with death registry data to understand under-reporting of influenza deaths.”

13. Page 4, line 78: Is it correct that all deaths in the Philippines are registered with cause of deaths within 48 hours? – Faster than in any other country, I know of.

Death registration in the Philippines is not always implemented as intended by policy. A study on death registration in one province showed that only 77% of deaths are recorded. However, since we do not know much about underreporting of deaths, we chose not to adjust it in our current modelling exercise.

Source: Carter, K.L., Williams, G., Tallo, V. et al. Capture-recapture analysis of all-cause mortality data in Bohol, Philippines. Popul Health Metrics 9, 9 (2011) doi:10.1186/1478-7954-9-9

14. Page 4, line 81-86: What do you mean by and what is the difference between ILI sentinel sites and SARI (Severe acute respiratory infections) sites? – how many of each?

ILI sentinel surveillance sites are health centers and hospital outpatient departments while SARI sentinel surveillance sites are hospital inpatient departments. The FluNet data did not specify how many sites there were, but according to this study which documented flu surveillance in the Philippines, from 2006 to October 2008 there were 18 health centers and 18 outpatient departments. Then by 2009, 16 more health centers were added. As of 2009, there were 52 sentinel sites. Unfortunately, we cannot find a dossier that gives an update on the number of sentinel sites from 2009-2015, 2015 being the last data point of our study. We have added these details in the methodology.

Lucero, M.G., Inobaya, M.T., Nillos, L.T. et al. National Influenza Surveillance in the Philippines from 2006 to 2012: seasonality and circulating strains. BMC Infect Dis 16, 762 (2016) doi:10.1186/s12879-016-2087-9

15. Page 5, line 93. I cannot find ‘average weekly temperature’ in reference 13, only average daily min and max temperatures.

We apologize for the confusion. The World Meteorological Organization has no specific guideline on average weekly temperature. We just applied their recommendation for average daily temperature to come up with the weekly temperature. Specifically, on page 78 of their manual, they say that “All ordinary climatological stations observe a daily maximum and minimum temperature … Hence, the recommended methodology for calculating average daily temperature is to take the mean of the daily maximum and minimum temperatures.”

To make our point clearer, we revised the concerned sentence into (lines 104-107): “The average weekly temperature was calculated by taking the mean of the weekly maximum and weekly minimum, following the logic of the World Meteorological Organization’s recommended method of computing for average daily temperature.”

Source: World Meteorological Organization. Guide to Climatological Practices 2011. [cited 16 July 2019]. Available from: http://www.wmo.int/pages/prog/wcp/ccl/guide/documents/WMO_100_en.pdf

16. Page 5, line 105: Ok to use a polynomial spline, but why 6?

We only have 6 years’ worth of data so we thought of keeping our polynomials up to 6. We also added the polynomial terms one-by-one (up to 6 polynomial terms) until we reach the model that provides the best fit.

17. Page 6, line 118: ‘Other missing data …’ - how many?

We did not have other missing data, apart from the number of tested samples for week 16 of 2008. We therefore deleted the “Other missing data…” sentence.

18. Page 13, line 212-215. I believe it is highly surprising that this studies all-cause influenza-associated mortality estimates are consistent with the respiratory-cause-specific estimates from Iuliano et al. study

We disagree that our results are surprising in this aspect. While we modelled all-cause mortality, we expect that most of these deaths would be respiratory deaths and cardiovascular deaths since these are the top causes of mortality in our population and other studies have linked flu to respiratory and cardiovascular deaths. Since our outcome includes respiratory deaths, our estimates should be close to Iuliano et al’s estimates.

Source: http://www.healthdata.org/sites/default/files/files/country_profiles/GBD/ihme_gbd_country_report_philippines.pdf

19. Page 13, line 223: I would not use the word ‘accurately’, but something like our model fitted data well

We removed “accurately” from our discussions; see line 247.

20. Page 13, line 223-4: Sometimes ‘official deaths registry statistics’ stand for cause-specific and here for all-cause. Suggest using all-cause and cause-specific.

We agree with the suggestion and have revised the sentence to use “all-cause” mortality rather than ‘official deaths registry statistics’.

21. Page 13, line 224: It is not correct that your model predicted the typhoon Haiyan peak! – see figure 1

We apologize for this reporting error. We have deleted that sentence in the manuscript.

22. Page 14-15, line 252-255: You might have used, what is often called the Goldstein index: ILI-rate * positive-percentage, where the ILI rate reflect the population dynamics and the positive-percentage limit the ILI-rate to influenza i.e. exclude other circulating respiratory pathogens.

To clarify, we did not use the Goldstein Index. What we used was the flu positivity rate which was computed by (specimens tested positive)/(total specimens tested). We do not have good estimates of daily ILI rates in the population precluding use of the Goldstein index.

Submitted filename: Response to reviewers.docx

Click here for additional data file.

2 Jun 2020

Influenza-associated excess mortality in the Philippines, 2006-2015

PONE-D-19-24478R1

Dear Dr. Cheng,

We are pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it complies with all outstanding technical requirements.

Within one week, you will receive an e-mail containing information on the amendments required prior to publication. When all required modifications have been addressed, you will receive a formal acceptance letter and your manuscript will proceed to our production department and be scheduled for publication.

Shortly after the formal acceptance letter is sent, an invoice for payment will follow. To ensure an efficient production and billing process, please log into Editorial Manager at https://www.editorialmanager.com/pone/, click the "Update My Information" link at the top of the page, and update your user information. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, you must inform our press team as soon as possible and no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

With kind regards,

Joël Mossong

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: All comments have been addressed

Reviewer #2: All comments have been addressed

Reviewer #3: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: (No Response)

Reviewer #3: Yes

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: (No Response)

Reviewer #2: (No Response)

Reviewer #3: No further comments, the manuscript have been correctly and thoroughly revised and I believe it's ready for publication.

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

4 Jun 2020

PONE-D-19-24478R1

Influenza-associated excess mortality in the Philippines, 2006-2015

Dear Dr. Cheng:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Joël Mossong

Academic Editor

PLOS ONE

Table 1

All-cause deaths and laboratory-confirmed influenza cases recorded in the Philippines, 2006−2015.

Year	Number of weeks FluNet data was available	All-cause deaths, N	Samples tested, N	Influenza-positive		Influenza A		Influenza B
				n	%	n	%	n	%
2006	52	439,772	5,955	557	9.4%	381	6.4%	176	3.0%
2007	52	440,651	6,291	536	8.5%	504	8.0%	32	0.5%
2008a	52	458,793	11,676	807	6.9%	211	1.8%	596	5.1%
2009	50	460,462	23,169	7,706	33.3%	7,556	32.6%	150	0.6%
2010	51	476,754	10,101	1,724	17.1%	875	8.7%	849	8.4%
2011	52	497,166	9,689	894	9.2%	627	6.5%	267	2.8%
2012	50	491,110	7,346	650	8.8%	403	5.5%	247	3.4%
2013	51	517,767	5,052	972	19.2%	388	7.7%	584	11.6%
2014	51	535,336	3,765	522	13.9%	269	7.1%	253	6.7%
2015	50	536,305	1,140	163	14.3%	149	13.1%	14	1.2%
Mean	51.1	485,412	8,418	1,453	17.3%	1,136	13.5%	317	3.8%

a Number of samples tested not available for week 16 of 2008.

Table 2

Estimated average annual influenza excess mortality per age group and annual excess mortality rate per 100,000 individuals, 2006−2015.

Age group	Influenza A & B		Influenza A		Influenza B
	Mean excess deaths, n (%)	EMR (95% CI)	Mean excess deaths, n (%)	EMR (95% CI)	Mean excess deaths, n (%)	EMR (95% CI)
0 to 4 y	249 (4.7)	2.14 (0.44−2.19)	189 (4.1)	1.62 (−0.02−1.66)	60 (7.9)	0.52 (−1.14−0.55)
5 to 9 y	169 (3.2)	1.51 (1.20−1.56)	116 (2.5)	1.04 (0.75−1.08)	53 (6.9)	0.48 (0.17−0.51)
10 to 19 y	100 (1.9)	0.48 (0.10−0.50)	84 (1.8)	0.40 (0.04−0.42)	16 (2.1)	0.08 (−0.29−0.09)
20 to 59 y	1,076 (20.1)	2.02 (2.00−2.03)	1,049 (22.9)	1.97 (1.96−1.97)	27 (3.5)	0.05 (0.04−0.06)
≥60 y	3,587 (67.1)	44.63 (44.51−44.69)	3,062 (66.8)	38.15 (38.07−38.18)	525 (68.8)	7.08 (6.45−6.57)
All agesa	5,347 (−)	5.09 (2.20−5.09)	4,584 (−)	4.37 (1.49−4.37)	763 (−)	0.73 (−2.15−0.73)

Abbreviations: CI, confidence interval; EMR, excess mortality rate per 100,000.

a Mean deaths for all ages were derived using the same estimation strategy as the other age groups, and therefore do not equal the sum of all ages.

Table 3

Estimated excess influenza-associated deaths versus nationally registered influenza deaths per age group.

Age group	2006	2007	2008	2009	2010	2011	2012	2013	2014	2015	Mean (2006−2015)
0 to 4 y
Estimated	165	176	121	400	277	195	183	367	296	309	249
Registered	11	1	4	3	5	1	3	5	0	1	4
5 to 9 y
Estimated	116	109	95	254	195	127	115	281	209	191	169
Registered	2	2	3	2	0	1	1	2	0	0	2
10 to 19 y
Estimated	66	74	38	177	99	80	66	129	133	139	100
Registered	8	7	5	2	6	3	3	2	2	1	4
20 to 59 y
Estimated	657	824	232	1,986	901	866	711	1,151	1,496	1,936	1,076
Registered	43	36	27	26	21	7	17	29	14	11	24
≥60 y
Estimated	1,979	2,317	1,192	5,740	3,291	2,823	2,492	4,822	5,233	5,978	3,587
Registered	99	81	85	58	52	26	47	121	32	35	64
All ages
Estimated	3,105	3,602	1,807	8,784	4,984	4,222	3,692	7,046	7,563	8,666	5,347
Registered	163	127	124	91	84	38	71	159	48	48	96

The estimated number of influenza deaths was compared to influenza deaths registered by the Philippine Statistics Authority with ICD-10 code J11.1 (‘influenza due to unidentified influenza virus with other respiratory manifestations’).