The 2009 H1N1 influenza outbreak highlighted the need for methods that not only detect an outbreak but also estimate incidence so that public-health decision makers can allocate appropriate resources in response to such an outbreak. Several studies [1–3] have attempted to estimate incidence during the 2009 H1N1 outbreak, but they were performed months after the outbreak and were based on data that are either not easy to collect (e.g., confirmed laboratory cases) or not available in a timely fashion (e.g., surveys). Some influenza surveillance systems were also in operation during the H1N1 outbreak, but they did not provide daily estimates of influenza activity [4,5].

Previous research indicates that data on over-the-counter (OTC¹) sales of retail healthcare products, such as cough and cold medications, can support the early detection of disease outbreaks [6,7]. In this paper, we provide support that thermometer sales (TS) increase significantly and early during an influenza outbreak.

We explore the hypothesis that TS data can be used to estimate accurately the number of patient cases that present to EDs and are symptomatically like influenza (SLI). We developed and evaluated a model to estimate the number of SLI cases in EDs in Allegheny County (AC), Pennsylvania based on TS. To train this model, we used ED visit data from hospitals in AC and TS data from AC. We then generalized the model to predict daily ED SLI cases from the daily TS in any region in which we are able to collect TS data. Unlike counts of ED SLI cases, which are not readily available for most counties in the US, there exists TS data for a significant number of locations throughout the US, with at most a one day delay. Given this availability, the method and model we present is intended to help public health epidemiologists throughout the US to estimate the extent of an influenza outbreak as it presents to EDs with little delay; the evaluation of this method beyond AC is proposed as a next step in the future pursuit of this line of research.

The remainder of this paper is organized as follows. In the next section, we provide an overview of previous research that aimed to detect and characterize influenza outbreaks using OTC sales. We summarize a project that facilitates obtaining daily OTC sales data across the US with little delay. We then discuss a method we developed to estimate SLI cases in AC. We describe the experimental methodologies we used to evaluate the performance of the models in AC, and we present and discuss the results of those experiments. We present methods to adjust TS data so that these models can produce estimates for other regions, not just AC. Finally, we present conclusions and close the paper with a summary.

This section provides background information about the use of OTC healthcare-product sales to detect disease outbreaks. We then provide an overview of the 2009 H1N1 influenza epidemic in the US and of several published studies that estimate its level of activity. The section then describes a system we have developed to estimate from ED reports the expected number of daily SLI cases that visit EDs in a region. We close the section by describing two influenza surveillance systems that can use the predictions produced by the models we present in this paper.

In this section we present a novel method for estimating SLI cases from OTC TS in EDs. Due to data availability, we created this model for estimating ED SLI cases in AC, but the methodology is flexible and can be extended to make estimates for other regions (Section 6).

To create a model for AC, we first gathered training data as follows. We obtained non-negative daily sales of pediatric and adult thermometers from May 1 through December 31, 2009 (n = 245 days) in AC from the NRDM. We will refer to the data from this period as the 2009 dataset. We compared these TS data with the expected number of daily SLI cases that presented to EDs monitored in AC during the same period according to BCD. Fig. 1 shows a plot of TS during that eight-month period (solid line) and the number of ED SLI cases estimated by BCD during the same period (dash-dotted line). We found that the strongest cross-correlation of the two curves is 0.9 when the expected cases lag behind the TS data by one day; a zero-day lag had almost the same cross-correlation (0.89). Fig. 2 shows a plot of the estimated cases as a function of TS, which supports a linear relationship between those two quantities (black line).

The rationale for the selection of the 2009 dataset is twofold. First, the BCD system officially started operating in AC in May 2009, and hence, we only had data starting at that time. Second, we wanted to devise a model that could estimate from thermometer sales not only SLI activity during outbreaks (high activity periods) but also during non-outbreak periods (low activity periods). The 2009 dataset contains sub-periods of both high (during the 2009 H1N1 outbreak) and low SLI activity and thus was the suitable for creating the aforementioned model. A model that is capable of estimating SLI activity only during outbreaks is more limited in its usefulness.

Eq. (1) describes a linear model of SLI as a function of TS for AC, where S and I denote the slope and intercept parameters of the model, respectively, ED_SLI_BCD represents SLI cases in the EDs in AC that have BCD deployed, and TS_NRDM denotes thermometer sales reported by NRDM for AC. (1)ED_SLI_BCD(AC,date)=I+S⋅TS_NRDM(AC,date)

Eq. (1) only predicts the number of SLI cases that present to EDs in AC that have BCD installed. To predict the total number of ED cases of those diseases in AC, we need to adjust that equation to incorporate the ED coverage of BCD in AC. Let ED_SLI_T(AC, date) be the total number of ED SLI cases (the T stands for total) during date in AC. We can estimate ED_SLI_T using Eq. (2), where EDcov is the ED coverage of BCD, which is discussed in more detail later in this section. (2)ED_SLI_T(Ac,date)=ED_SLI_BCD(AC)∕EDcov(AC)

We combine Eqs. (1) and (2) to obtain Eq. (3) that approximates the total number of daily SLI cases that present to all EDs in AC. (3)ED_SLI_T(AC,date)=I+S⋅TS_NRDM(AC,date)EDcov(AC)=1EDcov(AC)+SEDcov(AC)⋅TS_NRDM(AC,date)

At the time this research was performed, the ED coverage of BCD in AC was around 39%. We calculated it by dividing the number of visits to EDs in AC where BCD is deployed by the total number of visits to EDs in AC from May 1st, 2009 through December 31st, 2009. We obtained the latter count from a system that our laboratory has deployed in all EDs in AC, except for one ED in a federal hospital.

For the sake of example, suppose that in the model given by Eq. (3) the slope (S) is 0.5 and the intercept (I) is 20. If on a given day d, the net number of thermometers sold in AC is 10, then the model would predict that the total ED SLI cases in AC on that day is approximately 64 = 20/0.39 + (0.5/0.39)·10, where 0.39 is the estimated ED coverage of BCD in AC.

This section describes experiments that evaluate how well Eq. (1) performs when using AC data for training a model with data from one period of time of m days and testing with data from a subsequent period of n days. To do so, we create training and testing sets, perform least squares linear regression on the training set to obtain slope (S) and intercept (I) parameters, and then use Eq. (1) with those parameters to predict ED cases for each day in the test set. The next training period begin n days after the start of the former training period.

For example, for a training set of m = 42 days (6 weeks) and a test set of n = 14 days (2 weeks) we proceed as follows. We train a model with TS and BCD data from for example 5/6/2009 (a Sunday) through 6/16/2009; we then apply that model to predict BCD-estimated SLI cases from 6/17/2009 through 6/30/2009. The TS data and BCD estimates from 5/20/2009 (i.e., n = 14 days after 5/6/2009) through 7/02/2009 (i.e., 42 days after 5/20/2009) become the next training set.

To evaluate the performance of the models after the 2009 H1N1 outbreak had largely dissipated, we included the first three months of 2010. For the sake of uniformity, we start on 5/6/2009, which is a Sunday, instead of on 5/1/2009, and include data up to 4/3/2010, which is a Saturday. We will refer to this dataset of 333 days as the experimental dataset. A plot of that set is shown in Fig. B1 in Appendix A.

Based on the above arrangement, we performed three experiments with different training and testing set sizes. In all these experiments the smallest training set consists of three weeks of data, and larger training sets increase in periods of three weeks, with the largest set consisting of to 15 weeks. That is, m = 21, 42, 63, 84, and 105 days. We stopped at 105 days because that is approximately the number of days in the primary wave of the 2009 H1N1 outbreak in AC. The experiments differ in the number of days used for the testing set:

Experiment 1 (E1): Each testing set consists of 7 days (1 week) of data. The total number of training sets in this experiment is 190.

The rationale of picking the sizes of the testing sets in the first and second experiment is to see how well the model can estimate disease activity during periods of time close to the date when the estimate is produced, which may be of greater public-health interest. The reason for selecting the size of the testing set on the third experiment was to see how well the models performed on a testing set the same size as its training set.

We evaluated the goodness of the regression performed to create each of the linear models described above with the coefficient of determination, denoted as R². It measures the fraction of the variability of the modeled response variable (e.g., expected number of daily ED SLI cases computed by BCD) that is explained by a regression [18]. Relationships with a strong positive (or negative) slope that are close to being linear would be expected to have an R² value near 1. Associations with a slope near zero would be expected to have an R² value near zero. We calculate R² with Eq. (4) [18], where

m is the number of days in the training set used to generate a model (i.e., m can be 21, 42, 63, 84, and 105).

We evaluated the predictive performance of each of those models using a measure known as the mean absolute percentage error (MAPE) [19,20], which is computed using Eq. (5). A MAPE value of X means that, on average, a model's predictions are X% greater or smaller than the actual quantity. Thus, the closer the MAPE values are to zero (with zero indicating that the predictions and actual values are the same), the better is the forecasting performance of an estimation method. As in the case of R², the value of n in Eq. (5) for this experiment is the number of days in the testing set for which we generate a prediction (i.e., n can be 7, 14, 21, 42, 63, 84, and 105), the variable ActualValue_i in Eq. (5) is the number of SLI cases estimated by BCD on day i in that set, the variable PredictedValue_i represents a given model's prediction of daily ED SLI cases from thermometer sales on that same day, and the parallel bars (∥) denote the absolute value of the quantity they enclose. (5)MAPE=1n∑i=1n(∣ActualValuei−PredictedValuei∣ActualValuei⋅100)

To assess the forecasting performance of a method, Lewis [21] suggested four ranges of MAPE values, namely, [0%, 10%], (10%, 20%], (20%, 50%], and (50%, ∞). He labeled these ranges as having highly accurate, good, reasonable, and inaccurate forecasting performance, respectively. We use these as regions of reference to facilitate our analysis.

In the next section we summarize the performance for each experiment using the average MAPE, the minimum and maximum MAPE values, and the minimum and maximum R² values.

In this section we present and discuss the results of the experiments outlined in Section 4. These experiments evaluate models that estimate daily ED SLI cases in AC during different periods of time. The result of experiments E1, E2, and E3 are shown in Tables 2–4, respectively. Table 2 shows that experiment E1 achieves the lowest Min MAPE values, relative to the other experiments, and achieves average MAPE values of 16.7% or less. However, E1 also has the largest MAPE estimation errors, as given by the Max MAPE column in Table 2). Recall from Section 4 that a MAPE greater than 10% and less than or equal to 20% is considered to be good forecasting performance, whereas one greater than 50% is inaccurate performance [21]. Table 2 shows that the average MAPE is well within the good forecasting range. A more stringent criteria, however, is to consider the maximum MAPE over all 44 sets. When we do so, most (4 out of 5) of the maximum MAPE values in E1 are uncomfortably close to inaccurate performance. As a result, the models produced by the method presented in this paper to estimate SLI cases over a period as short as 7 days should be used with caution.

Table 3 shows that the average MAPE values are almost as good as those of E1 (the training sets are of the same sizes in both experiments but the testing sets are of different sizes, except for the first), but the maximum MAPE values are not as high. All of the former are in the “good” forecasting bracket, whereas all the latter are well within the “reasonable” forecasting bracket.

For experiment E3, Table 4 shows that the average MAPE values are not as good as those obtained in E1 and E2, but the maximum MAPE values are lower than in the other two experiments. It is also clear that in this experiment the errors for the last two configurations (m = n = 84, and m = n = 105) are higher than those of the other configurations for this same experiment. The latter seems to occur when the training set contains mostly data from a non-outbreak period but the majority of data of the testing set comes from an outbreak period (or vice versa).

We now briefly discuss individual results in experiments E1 and E2, rather than their aggregated statistics. Figs. 3 and 4 show, for experiment E2 and E3, respectively, the MAPE values for each testing set in the configuration where the training set consists of m = 42 days (6 weeks) of data and the testing set consists of i = 14 days of data. Appendix B includes similar figures for the other configurations in these experiments.² In these figures it can be seen that the forecasting performance of the models declines under two circumstances:

A decline occurs when there is a sharp drop in the number of SLI cases in a test set relative to the associated training set. For example, the testing set consisting of data from 10/29/2009 through 11/11/2009 (n = 14 days) shown in Fig. B2, is a period when the number of BCD-estimates SLI cases was falling sharply. However, its corresponding training set consisting of data from 10/7/2009 through 10/28/2009 (m = 21 days) coincides with the period of highest disease activity. A similar phenomenon happens for the latter (Fig. B2) and other configurations, for the testing sets that start at various points during December 2009 (e.g., Figs. 3 and B3), and during February 2010 (e.g., Figs. B4 and B6).

The performance decline in the two cases described above can be explained in terms of a mismatch between the training and testing periods, due to changes in disease activity between those periods. The forecast errors worsened both before periods of high disease activity and after periods in which the disease activity had peaked and was waning. Fortunately, such performance declines did not last for more than a single estimation period in the majority of cases. For example, this pattern is observed for the testing periods that start on 7/15/2009 and on 12/16/2009 in Fig. 3, and those that start on 9/9/2009 and on 12/02/2009 in Fig. 4. In both figures it can be seen that each of those periods is followed by an estimation period with a noticeably lower estimation error. The same transient decline is visible in the other figures as well. One possible approach to addressing this problem would be to combine the methods presented in the present study with outbreak detection techniques. The detection of the start of an outbreak and the beginning of its decay would then guide which data to use to train the models.

In Tables 2–4 the Max R² values were consistently associated with periods of high and dynamic disease activity, where we would expect strong positive relationships between TS and daily ED SLI cases. In contrast, Min R² values were generally associated with periods of no disease activity, where the slope of the relationship between TS and ED SLI cases is close to zero.

We close this section with a summary of the main findings from the experiments described here.

On the whole, the linear models achieved “good” forecasting (MAPE) performance under most circumstances, and with only a few exceptions they achieved “reasonable” forecasting performance otherwise.

In this section, we show how to extend the model developed in Section 4 to predict daily ED SLI cases in other US second-level administrative divisions (SLDs), such as counties and independent cities [22,23]. For the sake of simplicity, we will use the term “county” instead of the more precise term SLD. The main assumption behind this method is that the linear relationship we found for AC between daily TS and ED cases of SLI (Fig. 2) is also present in other counties. We believe this is a reasonable first-order approximation. However, due to the non-availability to us of ED reports from other counties, we do not attempt here to validate the accuracy of the AC model when applied to other counties, but leave it as future research.

Eq. (3) is only applicable to predict cases in AC's EDs based on TS in AC, as reported by NRDM. We extended the model by adjusting the AC model for differences between (1) the population size of AC and a given county, and (2) the TS in AC and that county. To do so, we follow a three-step procedure. First, we normalize thermometers sold to be per person (PP) in the total population (pop) of a county. We use NRDM's market share (MS) in that county to estimate the TS. NRDM's MS (which we denote with MS_NRDM) is the fraction of the total OTC sales by all retailers in a particular county that are reported by NRDM retailers. We explain how we obtain MS_NRDM information in Appendix A. In Eq. (6), TS_PP is thermometer sales per person, r denotes the county of interest, pop(r) represents its population, and d is a particular date. (6)TS_PP(r,d)=TS_NRDM(r,d)∕MS_NRDM(r)pop(r)

Second, we incorporate Eq. (6) into Eq. (3) and adjust it to produce an estimate for AC of daily ED SLI cases per person in the population. The result is shown in the following equation. (7)ED_SLI_PP(AC,d)=IEDcov(AC)+SEDcov(AC)⋅TS_PP(AC,d)⋅(pop(AC)⋅MS_NRDM(AC))pop(AC)=IEDcov(AC)⋅pop(AC)+SEDcov(AC)⋅TS_PP(AC,d)⋅MS_NRDM(AC)

Note that in Eq. (7), the term MS_NRDM(AC) is part of the calculation of the adjusted slope (Eq. (3)), whereas TS_PP(AC, d) serves as the independent variable. Hence, to calculate ED SLI cases in other regions, we replace TS_PP(AC, d), which is specialized for AC, with the more general term TS_PP(r, d), which is the thermometer sales per person in an arbitrary country r, as shown in the following equation. (8)ED_SLI_PP(r,d)=IEDcov(AC)⋅pop(AC)+SEDcov(AC)⋅TS_PP(r,d)⋅MS_NRDM(AC)

To obtain an estimate of the number of daily ED SLI cases in a county, denoted with ED_SLI, we simply adjust ED_SLI_PP by multiplying it by the population of the county. We obtain Eq. (9) as a result. Note that if we use r = AC, Eq. (9) reduces to Eq. (3), which is the model in the previous section to predict total ED cases for AC. (9)ED_SLI(r,date)=ED_SLI_PP(r,date)⋅pop(r)

Continuing with the example of the previous section (where for illustration purposes we assumed that the S and I parameters were 0.5 and 20, respectively), suppose that we wish to compute the total ED cases in some county SC in which a net of 30 thermometers were sold on day d by the retailers in SC that report to NRDM. Furthermore, suppose that the population of SC is 200,000 people, that MS_NRDM in SC is 0.6 (i.e., 60% of all TS in SC are made by NRDM retailers), and that MS_NRDM in AC is 0.7. Using the population of AC in 2009 (1,218,494 people [23]), Eq. (9) predicts that there were about 13 ED cases of SLI in SC on day d: (10)ED_SLI(SC,d)=(200.39⋅(1,218,494)+0.50.39⋅(30∕0.6200,000)⋅0.7)⋅200,000=12.9

Based on the results of the experiments presented here, we conclude that it is possible that linear models trained to estimate SLI ED cases based on thermometer sales can have good forecasting performance. We showed this in the case for Allegheny County (AC), Pennsylvania. In our experiments we found that under most testing conditions the methods achieved good forecasting performance, with only occasional, transient declines. Given these results and the timeliness with which the TS data can be obtained from NRDM, we think the methods presented in this paper can be helpful for epidemiologists and public health officials.

The results presented here raise the possibility that similar models could be applied in other regions of the country, where thermometer sales data are commonly available, but daily ED SLI counts are not. In this paper, we describe how to transform the AC model to produce such estimates. At the time that this study was performed, we did not have access to estimates of daily ED SLI cases in other counties in the US. Thus, we could not test how accurately the AC-derived model predicts ED SLI counts in other counties in the US, many of which do have TS data that are available through the NRDM. Even more broadly, it will be important to test how well the approach described here will work in countries other than the US. In addition, in the future, it will be useful to evaluate how well models developed using data from one outbreak are able to perform in making predictions in subsequent outbreaks. We view these tasks as important directions for future research.

We hypothesize that other signals of SLI activity, such as laboratory confirmed SLI cases, may exhibit a similar linear relationship with TS, as well as with the sales of other relevant OTC medications. Investigating such relationships is another promising direction for future research.