key: cord-1003635-bfi1sb80 authors: Fulcher, Isabel R; Boley, Emma Jean; Gopaluni, Anuraag; Varney, Prince F; Barnhart, Dale A; Kulikowski, Nichole; Mugunga, Jean-Claude; Murray, Megan; Law, Michael R; Hedt-Gauthier, Bethany title: Syndromic surveillance using monthly aggregate health systems information data: methods with application to COVID-19 in Liberia date: 2021-05-31 journal: Int J Epidemiol DOI: 10.1093/ije/dyab094 sha: 156551445e7c4e66ef880a35e99eb9ee2ac72529 doc_id: 1003635 cord_uid: bfi1sb80 BACKGROUND: Early detection of SARS-CoV-2 circulation is imperative to inform local public health response. However, it has been hindered by limited access to SARS-CoV-2 diagnostic tests and testing infrastructure. In regions with limited testing capacity, routinely collected health data might be leveraged to identify geographical locales experiencing higher than expected rates of COVID-19-associated symptoms for more specific testing activities. METHODS: We developed syndromic surveillance tools to analyse aggregated health facility data on COVID-19-related indicators in seven low- and middle-income countries (LMICs), including Liberia. We used time series models to estimate the expected monthly counts and 95% prediction intervals based on 4 years of previous data. Here, we detail and provide resources for our data preparation procedures, modelling approach and data visualisation tools with application to Liberia. RESULTS: To demonstrate the utility of these methods, we present syndromic surveillance results for acute respiratory infections (ARI) at health facilities in Liberia during the initial months of the COVID-19 pandemic (January through August 2020). For each month, we estimated the deviation between the expected and observed number of ARI cases for 325 health facilities and 15 counties to identify potential areas of SARS-CoV-2 circulation. CONCLUSIONS: Syndromic surveillance can be used to monitor health facility catchment areas for spikes in specific symptoms which may indicate SARS-CoV-2 circulation. The developed methods coupled with the existing infrastructure for routine health data systems can be leveraged to monitor a variety of indicators and other infectious diseases with epidemic potential. developed methods coupled with the existing infrastructure for routine health data systems can be leveraged to monitor a variety of indicators and other infectious diseases with epidemic potential. Key words: Syndromic surveillance, disease monitoring, COVID-19, infectious disease, time series modelling Background Limited infrastructure and testing capacity in many low-and middle-income countries (LMICs) requires the use of novel approaches to disease monitoring. [1] [2] [3] Over the past several years, many LMICs have invested in developing routine health information systems. These systems collect data on the presentation of particular symptoms, health service utilization and treatment outcomes at the patient level or aggregated at the facility level for a given time period. These data systems provide important operational support to LMIC health systems, nad recent reviews have found they have been underused in research and intervention monitoring. [4] [5] [6] Globally, the COVID-19 pandemic response has been challenged by limited information on the magnitude and spread of the virus. To that end, syndromic surveillance can be used to identify potential SARS-CoV-2 hotspots and thereby inform resource allocation, lockdown strategies and enhanced seroprevalence testing strategies. Many symptoms of COVID-19-including fever, cough, diarrhoea and difficulty breathing-are captured in existing routine health system data, making these data ideal for monitoring purposes. 7 Although syndromic surveillance cannot replace direct monitoring of the disease with specific SARS-CoV-2 diagnostics, this approach can be used as a rapid, cost-effective strategy that could help identify areas for more specific testing when resources are limited. In addition, such methods can be integrated within regular monitoring activities across a variety of indicators, potentially leading to early identification of future emerging diseases. In the early months of the COVID-19 pandemic, the Cross-Site COVID-19 Syndromic Surveillance Working Group, including researchers and representatives from sites in Liberia, Lesotho, Malawi, Haiti, Mexico, Rwanda and Sierra Leone, collaborated to develop syndromic surveillance tools appropriate for aggregated health facility data. In this paper, we detail: (i) the processes of data preparation; (ii) the statistical modelling approach, including the adaptation of the parametric bootstrap approach for constructing prediction intervals; and (iii) data visualization tools. These are supported with links to resources to support replication. We demonstrate these methods with an example of monitoring acute respiratory infections at health facilities in Liberia. Finally, we conclude with recommendations to support syndromic surveillance activities in other locations and for other outbreaks using these methods. Methods for COVID-19 syndromic surveillance using monthly aggregate data This research only contained data aggregated at the health facility level. Ethical approval was not required. Different types of routine health information systems are used in different countries, and each requires a different • Routine health information systems data can be leveraged for disease monitoring and detection of emerging infectious diseases with epidemic potential. • We used time series modelling to detect health facility catchment areas experiencing higher than expected rates of COVID-19-associated symptoms, potentially indicating SARS-CoV-2 circulation. • Although syndromic surveillance cannot replace the use of SARS-CoV-2 diagnostics for direct monitoring of the disease, it can be used as a rapid, cost-effective strategy to identify local areas for more specific testing when resources are limited. approach for syndromic surveillance. Here, we focused on monthly counts aggregated at the health facility level for a given indicator. The format of these data is typical of Health Management Information Systems (HMIS) in LMICs, such as the popular District Health Information Software 2 (DHIS2) which has been endorsed for disease surveillance by the World Health Organization. 8, 9 Five of the seven country sites (Liberia, Lesotho, Malawi, Rwanda and Sierra Leone) report monthly aggregated data for their syndromic surveillance indicators, and therefore were equipped to use the approach described herein. The availability of health indicators varied by country ( Table 1 ). The most common indicator selected for COVID-19 syndromic surveillance was acute respiratory infection (ARI). As persons with COVID-19 often present with symptoms typical of ARI, we hypothesized that COVID-19 infected individuals would possibly be classified as ARI cases in the health systems data-meaning that the presence of SARS-CoV-2 circulation in a facility catchment area could plausibly appear as an increase in cases of ARI at that health facility. 10 Fever and pneumonia were also common indicators tracked across countries. We developed an automated data processing pipeline to systematically clean, model and visualize routinely collected data on respiratory infection indicators for each country. Annotated R code was developed for each country indicator to streamline data cleaning for future monthly extractions. Cleaned indicators typically included raw counts, such as the number of ARI cases reported to a health facility in a given month. Since March 2020, sites have submitted new indicator data on a monthly basis to be run through the data processing pipeline. Potential outliers were reported to the in-country monitoring and evaluation (M&E) officers for review. If the outlier data were suspected to be miscounted and determined by the M&E officer to be unresolvable, the value for the specific health We considered 1 January 2016 to 31 December 2019 as the baseline period across all sites, and data from this window were used to establish an expected (predicted) baseline count. For facility-level assessments, we fit a generalized linear model with negative binomial distribution and log-link to estimate expected monthly counts: where Y indicates monthly indicator count, t indicates the cumulative month number, K indicates the number of harmonic functions to include (we take K ¼ 3). The year term captures long-term annual trend and the harmonic terms capture seasonality. This mean model was chosen to allow smoothing without imposing strong assumptions on the seasonal behaviour, and also aligns with other models used during the pandemic. 11 We could have alternatively chosen to model trend as a monthly (t) or quarterly linear term, but this specification was chosen as it performed well across indicators and facilities. Further, if information is available on external covariates, such as annual rainfall, it would be possible to incorporate this information into the modelling procedure. We chose to use a negative binomial distribution (instead of Poisson) to account for overdispersion. The baseline model provided a predicted count with corresponding 95% prediction intervals for each health facility in a given month. To calculate the prediction intervals, we used a parametric bootstrap procedure, drawing realizations for the model coefficients from a multivariate normal distribution and using this to then compute predicted counts from a negative binomial distribution. This is similar to the methodology employed by the ciTools R package but uses a negative binomial distribution instead of quasi-Poisson. We defined the evaluation period as beginning on 1 January 2020, based on the first investigation of atypical pneumonia cases in Wuhan, China by the World Health Organization in early January. 12 We used the baseline model, which was fit using only data from the baseline period, to predict indicator counts with 95% prediction intervals for months during the evaluation period. This prediction reflects what we would expect to observe at a specific facility in a given month in the absence of any impact from COVID-19. We defined a deviation as the difference between the predicted and observed count (or proportion) for a given month. To facilitate interpretation across months, indicators and facilities of different sizes, we standardized the deviation measure by dividing by the predicted count. A negative value indicated a count (or proportion) less than expected and a positive value indicated a count (or proportion) higher than expected based on the baseline model. In our data visualizations, we reported this scaled deviation measure and indicated if the observed count fell outside the 95% prediction interval. Facilities with 20% or more missing observations from the baseline period (equivalent to 10 out of 48 months) were excluded from this facility-level syndromic surveillance activity. The consensus of the M&E officers at our collaborating sites is that higher levels of missing data likely indicated broader data reporting issues that may (i) compromise the baseline model and resulting prediction intervals and (ii) potentially perpetuate in the reported counts during the COVID-19 evaluation period. For the facilities included in the analysis, we assumed residual missing monthly counts were at random and conducted a complete case analysis for the baseline model. This missing-at-random assumption would be violated if counts were missing for a reason related to the count value (e.g. values more likely to be missing during high patient load months). Importantly, health-seeking behaviour may change during the course of a disease outbreak, which can affect the expected number of cases of a particular syndrome presenting to a health facility. During the Ebola Virus Disease outbreak, health care utilization declined by 18% across Liberia, Guinea and Sierra Leone over 2013-2016. 13 To account for such changes in overall health service utilization, we also modelled the proportion of indicator counts at that health facility, with a generalized linear model with negative binomial distribution and log-link using an offset term, D t , the total number of monthly outpatient visits: Similar to the counts, we performed a parametric bootstrap procedure to calculate the prediction intervals. If there were missing values for the monthly outpatient visit count (denominator), we conducted an additional step to impute the missing outpatient visit value in the parametric bootstrap. We repeated the imputation procedure 500 times and took the median count for each month to be the predicted value and the 2.5th and 97.5th percentiles as the 95% bootstrap prediction interval. Further details on this procedure are provided in the Supplementary File 2, available as Supplementary data at IJE online. It may be of interest to report syndromic surveillance results on a wider geographical level, such as a district or county, that contain multiple health facility catchment areas. If there are no missing monthly counts at the facility level, one could simply sum the count indicator for each month across all facilities within a geographical region and fit the model in Equation 1 . In this case, although both the regional and health-facility models may be valid, the predicted counts from each set of models would not be equivalent: see Supplementary File 2 for explanation. However, in the presence of missing monthly data, regional models need to account for missingness at the facility level for a given month. We continued to exclude facilities that had more than 10 months (20%) or more of missing data during the baseline period. Even after this exclusion, there was still some missingness at the facility level which would impact on the summed counts at the wider geographical level. We used a parametric bootstrap to impute any missing values from the facility level models. We drew realisations of the predicted indicator counts from the facility-level model for each month and facility and then summed these values for a district-(or county) levelestimate. We repeated this procedure 500 times and took the 2.5th and 97.5 th percentiles to create 95% prediction intervals. For proportions, the numbers of outpatient visits were summed across facilities and a proportion was computed. If there were missing values in the outpatient visits, another step was included in the above parametric bootstrap procedure where missing outpatient visits were generated from equation (1) where Y indicates monthly outpatient visit count. The details on this procedure are given in the Supplementary File 2. We note that bootstrap procedures have been used to account for missing values in previous literature 14 but have not yet been used to aggregate time series results in multi-level data, such as aggregative facility-level data to the county-level. To ensure that the specified regression models do not have autocorrelation in the residuals, we performed several diagnostic procedures. For the residuals from facility-and county-level models, we generated plots of residuals by time, autocorrelation functions and partial autocorrelation functions. We also conducted the Breusch-Godfrey test for serial correlation for the facility-level models. Data visualizations were an essential tool for communicating results back to country sites. Time series plots provide information on all observed data points in the baseline and evaluation periods in addition to the predicted values and 95% prediction intervals. Although this provides granular information for individual health facilities and indicators, it is difficult to compare information across multiple health facilities or indicators. To facilitate comparison data, we report the difference in case counts between observed and predicted per 100 000 persons. In the Supplementary materials, we provide two additional visualizations: tiled heat maps (Supplementary File 3, available as Supplementary data at IJE online) and interactive geographical maps (Supplementary File 4, available as Supplementary data at IJE online). Tiled heat maps were coloured based on the deviation from expected during the evaluation period months and flagged when counts (or proportions) were outside the 95% prediction intervals. Geographical maps show the deviation from expected across geographical regions for a single time point. These maps were made interactive to allow for toggling between indicators and months using the leaflet R package. All analyses and figure generation were done in R V3.6.0. Our code with working examples is available on our public GitHub repository to facilitate adoption of this methodology in related contexts [https://github.com/isabelfulcher/global_covid19_ response]. We present an application of the syndromic surveillance in Liberia for January through August 2020. Specifically, we detail the monthly counts of acute respiratory infection (ARI) cases and proportion of outpatient visits with a focus on Maryland County, the seventh largest county in Liberia. We present results at the facility level for JJ Dossen Hospital, the largest health facility in Maryland County [median of 262 monthly ARI cases with an interquartile range of (228, 299) during baseline] and at the county level for Liberia's 15 counties. We use the time series models presented in Equations (1) and (2) with K ¼ 3 and a baseline period of January 2016 through December 2019 for the facility-and county-level models. We report the facility-level models for ARI counts and proportions (Figure 1 ). In January 2020, JJ Dossen Hospital had a higher than expected number of ARI cases with 470 cases observed [compared with an expected 258 cases with a 95% prediction interval of (145, 414)] which promptly dropped to below expected after February 2020. However, from January through April 2020 the proportion of observed outpatient visits that were ARI cases (13-14%) was consistently higher than the proportion that was expected (7-8%) . Similarly, at the Maryland county level, which includes JJ Dossen Hospital and 23 other facilities, the ARI case count and proportion remained higher than expected through March and then promptly fell within or slightly below the predicted range until July (Figure 2 ). Residual plots for these time series models are provided in Figure 3 and do not exhibit signs of residual autocorrelation. Breusch-Godfrey tests were performed for JJ Dossen Hospital with a p-value > 0.05 for both counts and proportions. For the 15 county-level models, 325 facilities met the criteria for inclusion, with at least 80% complete acute Table 2 shows the breakdown of included facilities by county. Maryland, Bong, Grand Kru and River Gee have at least 90% of their facilities included, whereas Montserrado, Margibi, Grand Bassa and Rivercess have less than 40% of facilities included. The percentage of facilities missing ARI counts in 2020 was similar to 2019 with the exception of Margibi, which had more facilities with incomplete data in 2020, and Bomi, which had less facilities with incomplete data in 2020 ( Figure 4) . In January 2020, all 15 counties reported higher than expected ARI case counts, with seven counties' counts outside the 95% prediction interval ( Figure 5 ). Of the seven counties, the largest deviations were in Rivercess and Maryland with 842 and 556 higher than expected cases per 100 000 persons, respectively. For the proportion measures, only four of these counties-Maryland, Grand Kru, Rivercess and Sinoe-had deviations larger than the 95% prediction interval, which persisted in Maryland and Grand Kru through March 2020 ( Figure 6 ). Interestingly, five counties-Bong, Montserrado, Nimba, River Gee and Lofa-had consistently lower than expected ARI proportions after April 2020. The first confirmed COVID-19 case was on 16 March 2020 in Liberia. 15 The higher than expected number of ARI cases in January 2020 across multiple counties is unlikely due to COVID-19, as global evidence indicates that there was not sustained transmission of SARS-CoV-2 outside China during this time. 16 However, in the case of Maryland County, which had a large spike in ARI cases from January through March, it is possible that the containment measures taken for COVID-19 following the national health emergency contributed to a reduction in ARI cases in the months of April through August. 17 The lack of significant increases detected in the number and proportion of ARI cases seems to align with the low number of confirmed COVID-19 cases in Liberia during March through August 2020. Cumulatively from 16 March to 31 August, there were only 38 confirmed cases per 100 000 persons (1305 total) in Liberia, with 24 cases per 100 000 persons (32 cases total) from Maryland County (see Supplementary File 5, available as Supplementary data at IJE online). To detect a significant deviation in ARI cases in Maryland County, there would need to have been at least 362 additional ARI cases per 100 000 persons on a monthly basis-more than 75 times the number of confirmed COVID-19 cases in the county. Although the comparison between confirmed COVID-19 cases and ARI cases helps contextualize the interpretation of these syndromic surveillance results, there are two important caveats: (i) the number of confirmed COVID-19 cases may be underestimated if testing capacity is low; and (ii) individuals with COVID-19 may not present as an ARI case at a health facility, either because they are not exhibiting ARI symptoms or they did not seek care. We developed methodology to automate data cleaning, time series modelling and data visualizations for identification of deviations in COVID-19-associated symptoms. Although the presented results were for an 8-month period, we continue the data processing and modelling procedure on a monthly basis when new data become available. The monthly updates are shared with the country's health management team to investigate facilities that have higher than expected ARI cases. The results will now inform which facilities will be chosen for a seroprevalence study among health care workers in Maryland and Montserrado counties to validate these methods. The use of data on COVID-19-associated symptoms for syndromic surveillance is subject to several limitations that should be considered before employing such methods. First and foremost, syndromic surveillance is not equivalent to directly monitoring a disease via widespread viral or antibody testing, and the choices of indicators by each country have not yet been validated for monitoring purposes. We are currently validating these methods and indicators in regions that had high levels of COVID-19 testing. Second, this exercise relies on the availability of both previous data to establish a baseline and future data to detect deviations-large amounts of missing data, variations in the quality of data over time and lack of information on important COVID-19-related symptoms will hinder this effort. Changes in indicator definitions during outbreaks may also affect the utility of these methods; for example, patients with suspected COVID-19 may no longer be included in the acute respiratory infection indicator. Importantly, this was not the case in Liberia. Third, substantial changes in health-seeking behaviour or in data entry during the pandemic will impact on the interpretation of any deviations (or lack thereof). In the case of Liberia, as patients became informed and/or educated on COVID-19 case definition, they may have under-reported ARI symptoms or chosen not to go to a facility for fear of being quarantined. Additionally, clinicians may not have formally diagnosed clinical symptoms as an ARI because of concern that would trigger COVID-19 investigation; this could potentially explain some of the much lower than expected ARI cases and proportions. Last, countries must also have the capacity to develop a streamlined process for data cleaning, analysis and visualization that can be readily updated when new data are available. To facilitate a streamlined model fitting every month, the same parametric model specification for counts and proportions was used across all facilities (i.e. Equations 1 and 2). In principle, one could fine-tune each model for each facility and indicator by investigating autocorrelation functions and comparing various model specifications. However, this was not feasible in our setting with seven countries and a multitude of facilities and indicators within each. We instead recommend using the Breusch-Godfrey test for serial correlation as a way to flag facility-level models that need additional investigation (based on a prdetermined significance cut-off). In addition, we found it beneficial to have the granular time series plots for all facilities and indicators accessible on an online tool hosted by Shiny App. This enabled site leads to visualize the results and flag potentially ill-fitting models to the study team. Further, we note that excluding facilities with high levels of missingness may underestimate the raw indicator counts at wider geographical levels (e.g. Montserrado County in Liberia). However, because facilities with high levels of missingness are excluded from both the baseline models and the observed count, comparing deviations from expected counts can still be used for the purposes of syndromic surveillance, as long as reporting rates are independent of ARI caseloads and time period. We attempted to alleviate this concern by Figure 4 which showed similar reporting rates across most counties during the COVID-19 pandemic and the preceding year. Finally, although these models were developed in response to the COVID-19 pandemic, this approach could be used and adapted for ongoing surveillance on a wide range of health indicators that are included in most health information systems. In general, deviations from expected should warrant further investigation by local public health officials, but do not necessarily indicate an outbreak or emerging infectious disease. When tracking symptoms related to a specific disease, such as COVID-19, deviations from expected can identify local areas for more specific testing when resources are limited, but need to be further validated before being used to inform resource allocation or lockdown strategies. In conclusion, syndromic surveillance can be used to monitor for potential COVID-19 or other disease outbreaks in health facility catchment areas. We leveraged data from existing data collection systems and developed data science tools in open access software to routinely monitor for COVID-19 symptoms on a monthly basis. The methods and accompanying resources can be applied to other regions or diseases with minimal adaptation. The data underlying this article will be shared on reasonable request to the corresponding author. We have provided an example dataset on our GitHub repository: [https:// github.com/isabelfulcher/global_covid19_response]. Supplementary data are available at IJE online. Intensive care unit capacity in low-income countries: a systematic review Sustainable clinical laboratory capacity for health in Africa Access to lifesaving medical resources for African countries: COVID-19 testing and response, ethics, and politics Using routine health information systems for well-designed health evaluations in low-and middle-income countries Using routine health information data for research in low-and middleincome countries: a systematic review Using district health information to monitor sustainable development Clinical characteristics of coronavirus disease 2019 in China DHIS2. DHIS2 In Action World Health Organization. DHIS2-Based Surveillance Support Tools Clinical features, diagnostics, and outcomes of patients presenting with acute respiratory illness: a retrospective cohort study of patients with and without COVID-19 Timeline of WHO's Response to COVID-19 Estimation of excess deaths associated with the COVID-19 pandemic in the United States Utilization of non-Ebola health care services during Ebola outbreaks: a systematic review and metaanalysis A unified approach to measurement error and missing data: overview and applications Health Alert: Liberia, Confirmed COVID-19 Cases Timing the SARS-CoV-2 index case in Hubei province Ministry of Health and National Public Health Institute of Liberia (NPHIL) We thank our team of junior analysts for their contributions to this ongoing effort: Donald Fejfar, Katherine Tashman None declared.