We retrieved OTC product sales for Allegheny County, Pennsylvania (270 stores) from the NRDM database, covering the period from June 1, 2016, to December 30, 2021. NRDM product codes were aggregated into 7 analytically defined OTC categories (Table 1). The full time series was partitioned into a training interval (June 1, 2016-May 31, 2019) and a held-out test interval (June 1, 2019-December 30, 2021). The training interval was used exclusively for model parameter estimation (refer to the Model Development section), whereas the test interval was reserved for out-of-sample model evaluation and prospective abnormal event detection (refer to the subsequent sections under Methods).

To support the method that we used for estimating parameters in the NRDM Tracker model, we derived the daily incidence rate of ILI in the population, which serves as the prior probabilities in the model. To do so, we obtained daily counts of ED encounters associated with ILI-related viruses in the region, including adenovirus, rhinovirus or enterovirus, influenza, human metapneumovirus (HMPV), parainfluenza (PIV), and respiratory syncytial virus (RSV). For each day, we calculated the proportion of encounters related to ILI among all ED encounters. ILI counts were identified based on either positive laboratory test results or relevant ICD (International Classification of Diseases) diagnosis codes (Multimedia Appendix 1). The participating EDs included (1) University of Pittsburgh Medical Center (UPMC) Presbyterian University Hospital, (2) UPMC Children’s Hospital of Pittsburgh, (3) UPMC Shadyside, (4) McKeesport Hospital, and (5) UPMC Mercy Hospital. These 5 sites collectively represent a significant fraction (about 60%) of ED visits in Allegheny County and 30%‐40% of western Pennsylvania. The research data were sourced from the University of Pittsburgh Neptune Research Data Warehouse [26].

The NRDM and ED datasets share strong geographic and population overlap, as both focus on Allegheny County, Pennsylvania. We retrieved OTC sales data from 270 retail stores across the county, while the 5 EDs used for estimating prior disease probabilities also serve this region. The coverage of the NRDM for the Pittsburgh metropolitan statistical area, including Allegheny County, is 73% of all retail stores in the area. Although the datasets are not individually linked, their shared regional focus supports the use of ED-based ILI daily incidence estimates in helping model NRDM purchasing patterns.

Since many OTC products are used to relieve symptoms rather than treat specific diseases, we grouped conditions into two broad categories: (1) ILI diseases and (2) all other diseases, which we label as OTHER. The main abbreviations we used to describe our methods are summarized in Table 2.

As described further in this study, we applied the expectation-maximization (EM) algorithm [27] (refer to Algorithms 1 and 2 in Multimedia Appendix 2) to estimate the conditional probabilities Pd(D(p)|dxi), where Pd(D(p)=productj|dxi=ILI) represents the probability a customer purchases product j on day d due to ILI either in that customer or in someone they are making the purchase for (eg, a family member). Pd(D(p)=productj|dxi=OTHER) denotes the corresponding probability for patients with conditions other than ILI.

In applying EM to estimate the conditional probabilities above, we made the following assumptions:

The EM algorithm iteratively refines probability estimates through two main steps: the expectation step (E-step) and the maximization step (M-step).

In the E-step, the algorithm computes the posterior probabilities Pd(ILI|productj) and Pd(OTHER|productj) for each product j on day d, using the current estimates of P(productj|ILI) and P(productj|OTHER), along with the prior probabilities of ILI and OTHER on day d estimated from ED data (refer to details in the Research Data section).

In the M-step, the algorithm updates the conditional probabilities, P(productj|ILI) and P(productj|OTHER) by maximizing the expected likelihood. This is achieved by leveraging the weighted counts of product j purchases on day d, denoted as, Ndj, using the posterior probabilities obtained from the E-step.

The E-step and M-step are repeated iteratively until convergence is achieved—that is, until the maximum change in estimates of the 2 conditional probabilities falls below a predefined threshold, tol.

We used ED data in the Allegheny County region to estimate the priorprobabilities, Pd(dxi), used by the EM algorithm. To smooth the ED time series, we applied a centered 7-day moving average to the daily encounter counts. At the beginning and end of the series, where fewer neighboring values were available, a smaller window size (eg, 4 or 5 days) was used to compute the average. This adaptive window ensures that smoothing can still be applied at the edges without discarding data. Based on the resulting smoothed counts, we calculated daily prior probabilities for ILI and OTHER that were used by the EM algorithm.

After estimating the conditional probabilities P(productj|dxi), we derive the expected number of customers with disease group dxi on day d, denoted Ed(dxi), using equation 1. Specifically, Ed(dxi) is calculated as the sum of the posterior probabilities Pd(dxi|D(p)) across all purchases on day d. For each individual purchase p, the posterior probability is computed via Bayes' rule (equation 2), where the numerator is the joint probability, calculated as the product of the conditional probability of the purchase given the disease (as estimated via the EM algorithm) and the prior probability of the disease. The denominator normalizes this quantity by summing over the 2 disease group categories that we model.

In contrast to the EM training phase, where disease priors were learned from ED data, we fixed the priors at Pd(ILI)=0.05 and Pd(OTHER)=0.95 for each day during the calculation with equation 3 to focus the contribution of NRDM data in estimating daily disease patterns. These values are within the bounds of the probabilities for ILI and OTHER estimated from ED data during the training phase, during which the daily ILI daily incidence ranged from 0 to 0.06 and the daily OTHER daily incidence ranged from 0.94 to 1 (refer to the Descriptive Statistics section; Figure 1).

For each day d, we used a probabilistic metric, Pd(D), as an indicator of the presence of unusual purchases on that day. This metric represents the probability of the observed data on day d (ie, the number and types of OTC purchases on that day). We refer to this probability as the daily likelihood. Using equation 3, we estimate this quantity by multiplying the probabilities of all individual purchase events Pd(D(p)), assuming that these purd purchase events are independent. A purchase behavior typically reflects a customer’s (or their household member’s) illness. However, the specific underlying disease is not recorded in the NRDM database. Therefore, for each purchase behavior, we sum the joint probabilities, Pd(D(p),dxi), across the 2 disease groups that we modeled.

Each joint probability is computed as the product of the prior probability of disease dxi on day d, denoted as Pd(dxi), and the conditional probability Pd(D(p)|dxi), which is the probability of observing the purchase behavior given the presence of disease group dxi.

To estimate these conditional probabilities in the model training period, we applied the EM algorithm, as discussed in the Model Development section. As mentioned, when training the NRDM Tracker model, we used the daily incidence of diseases in EDs in the same region on day d to estimate the prior probability, Pd(dxi). When applying the model to monitor for disease outbreaks, we used a uniform prior for Pd(dxi), Pd(ILI)=0.05 , Pd(OTHER)=0.95, which does not require having access to ED data.

We also derive the following normalized likelihood (Equation 4), which is the geometric mean of the daily likelihood across all purchases on day d:

where purd is the total number of the products purchased in the 7 OTC categories on day d. Gd(D) reflects the (geometric) average likelihood of individual product purchases. Changes in the unnormalized likelihood, Pd(D), may be due to changes in the relative distribution of the products being sold, the number of products being sold, or both. In contrast, changes in Gd(D) are only due to changes in the relative distribution of the products being sold.

This section describes the methods that NRDM Tracker applies to monitor a test set of NRDM data for disease outbreaks. It uses the log of the likelihood given by equation 3 to look for low values that may indicate an outbreak. When values are low, it can analyze the data at the product level to look for patterns that may indicate which product categories are being influenced by a putative outbreak; such a characterization may suggest the type of outbreak that is occurring. In the remainder of this section, we discuss how we derive P values at the product level.

For each product j, we obtain the daily sales count ndj for day d and computed a right-sided P value based on the empirical distribution of daily sales during a control period from day d-38 to day d-8, assuming a normal distribution. Approximately 69% (4326/6284) of Shapiro-Wilk normality tests [28] passed the normality check, indicating that most control period count distributions were consistent with a normal distribution.

Since a P value is calculated for each day, over many monitored days, we applied the Benjamini-Hochberg procedure [29] to control the false discovery rate (FDR) at a significance level of 0.05. Specifically, we used the multiple tests function from the Python statsmodels.stats.multitest module (statsmodels 0.15.0, 2024) with the method=‘fdr_bh’ option to correct for multiple hypothesis testing. This method adjusts raw P values by scaling each value according to its rank and the total number of tests m. Refer to equation 5, where m is 6284 in our experiment, pj are raw P values, and the term j is the rank of the P value in ascending order. A reverse cumulative minimum is then applied to ensure monotonicity across the adjusted P values (equation 6), as described by Dudoit et al [30] and implemented in the Statsmodels package. These steps ensure that any day with an adjusted P value at or below .05 is considered a statistically significant anomaly, while the expected proportion of false discoveries is controlled. Importantly, the adjusted P values, qk, are guaranteed to be nondecreasing with respect to the ordered raw P values.

The research protocol of the study described in this paper was approved by the institutional review board (IRB) of the University of Pittsburgh (study 20030193). As no patients were directly enrolled in this study and no compensation was provided, the IRB granted a waiver of informed consent. The NRDM dataset contains only aggregated purchase counts and does not include any individual-level purchase information. For the ED encounter data, a UPMC-approved honest broker deidentified the records to satisfy the requirements for a limited data set before distribution to our research team.

This section presents the performance and behavior of NRDM Tracker. We begin by summarizing descriptive statistics of daily OTC product sales across 7 symptom-related categories (Descriptive Statistics section). Next, we describe the estimation of model parameters using the EM algorithm (Model Parameters section). We then evaluate the model’s ability to detect population-level anomalies during the COVID-19 pandemic using NRDM data (Detection of COVID-19 section). Finally, we assess the statistical significance of daily deviations in product sales across individual categories to identify which product categories are influenced by a putative outbreak (refer to the Detection of COVID-19 section).

Table 3 lists descriptive statistics of the daily counts of the 7 categories from 274 stores in Allegheny County, PA from June 1, 2016, to December 30, 2021. The median of the daily OTC sales count over the 7 categories is 4062, with a wide range from 259 to 15,979, indicating substantial daily variation. The cough or cold category has the highest average daily sales (3168), accounting for most of the total OTC purchases recorded in NRDM. Moderate-sale categories include antifever for adult, baby or child electrolytes, and antidiarrheal. Antifever for pediatric, hydrocortisones, and thermometers have relatively low average daily sales. Most categories show right-skewed distributions (mean>median), suggesting occasional sales surges. The large variances (especially in total and cough or cold categories) indicate irregular and potentially seasonal demand.

Figure 2 illustrates seasonal patterns in OTC product sales over 5 years. Both total and “cough or cold” categories exhibit strong seasonal patterns, with recurring peaks during the winter months. A prominent spike in early 2020 likely corresponds to the onset of COVID-19, followed by a dip and gradual return to seasonal trends.

Sales trends among other product categories reveal varying degrees of seasonality and pandemic impact. Among antifever products, the adult antifever category shows relatively stable daily sales with modest fluctuations until a sharp surge in 2020, a pattern likely explained by the dual role of acetaminophen, ibuprofen, and related medications as both antipyretics and analgesics, which contributes to their steady baseline demand. The pediatric antifever category follows a similar seasonal trend but at a lower scale. Cough or cold products dominate in both volume and seasonality. Baby or child electrolytes show periodic increases with peaks in the winter. Antidiarrheal sales exhibit moderate daily variation with peaks around December. The peaks in hydrocortisone sales are less pronounced but show a consistent seasonal pattern, with increased use in warmer months, plausibly due to higher incidence of skin irritation, insect bites, and sun-related conditions. Thermometer sales remain low at baseline but exhibit dramatic, event-driven spikes—particularly in early 2020—likely due to increased public concern during the onset of the COVID-19 pandemic.

Figure 1 illustrates the temporal evolution of the estimated daily incidence rate of ILI (used as the prior probability in the model), derived from ED case counts for adenovirus, rhinovirus or enterovirus, influenza, HMPV, PIV, and RSV. The series exhibits clear seasonal peaks occurring annually—typically during the winter months (eg, early 2017, early 2018, spring 2019, and early 2020)—consistent with typical seasonal trends of respiratory viruses in this region. A pronounced decline is observed in 2020 and 2021, reflecting reduced ILI-related ED visits during the initial phase of the COVID-19 pandemic in the United States. A modest resurgence is seen in late 2021, coinciding with the return of seasonal respiratory virus activity.

Using the EM Algorithm, we estimated the conditional probabilities P(productj|ILI) and P(productj|OTHER). The EM algorithm converged after 36 iterations. Table 4 presents the estimated conditional probabilities, indicating the likelihood that each product category is purchased given an ILI-related or OTHER encounter. These values reflect how specific each product is to ILI activity. Cough or cold medications show a notably higher probability under ILI (0.792) than under OTHER (0.695), suggesting strong relevance to ILI-related symptoms. In contrast, the other products have lower probabilities for ILI than for OTHER, indicating they are more likely to be associated with non-ILI conditions. These conditional probabilities are central to inferring ILI activity from purchase patterns and help distinguish symptom-specific products from those with more general therapeutic use.

To assess the robustness of the assumption of time-invariant conditional probabilities, we conducted a sensitivity analysis using time-varying conditional probability tables (CPTs) estimated with a rolling 12-month EM framework and weekly updates (Multimedia Appendix 3). As illustrated in Figure S1 and Table S2 in Multimedia Appendix 3, the estimated CPTs exhibit moderate temporal variation, reflecting changes in product use patterns and disease presentation over time. However, when these dynamic CPTs were incorporated into the NRDM-ILI estimation pipeline, the downstream impact was limited. The resulting NRDM-ILI time series closely matched those obtained using constant CPTs (Figure S2 in Multimedia Appendix 3), with nearly identical correlations with ED-ILI (0.6642 vs 0.6599). These findings suggest that, although conditional probabilities are not constant, the assumption of time-invariant CPTs is reasonably robust in this setting, and the primary conclusions of the model are not very sensitive to this assumption.

Although fever is a defining symptom of ILI, antifever medications show relatively lower conditional probabilities under ILI compared with cough or cold medications. This likely reflects differences in product use behavior rather than inconsistency with clinical definitions. Similar patterns are observed for other product categories, including those associated with ILI symptoms, where medications that are less specific to respiratory illness tend to show weaker associations. For example, thermometers, which are used for general fever monitoring rather than specifically for respiratory illness, show very low probabilities under both ILI and OTHER. These patterns are consistent across both the main analysis and the supplementary sensitivity analysis (Multimedia Appendix 3).

Using equations 1 and 2, the NRDM Tracker calculated the expected counts of ILI from the NRDM data for the period from June 1, 2019, to December 30, 2021. The NRDM model was trained using data from June 1, 2016, to May 31, 2019, to estimate the CPTs; therefore, this latter period represents out-of-sample model performance.

Figure 3 illustrates the moving averages of those expected ILI counts in blue. The yellow curve shows ILI counts derived from the ED cases for the following viruses: adenovirus, rhinovirus or enterovirus, influenza, HMPV, PIV, and RSV, and shows strong temporal alignment with the curve in Figure 1, which depicts daily ILI priors. A disease caused by 1 of these viruses was considered present if a patient had a positive laboratory test result for the virus (eg, a positive PCR [polymerase chain reaction] test for influenza) or an ICD code indicative of the infection (listed in Multimedia Appendix 1). These ED-based ILI tests and diagnoses excluded COVID-19 cases, which emulates the situation in which a novel disease may occur before it is explicitly identified and tested for. The Pearson correlation coefficient between 2 time series shown in blue and yellow is 0.66.

The peak day of the NRDM-estimated ILI counts is March 13, 2020 (based on the underlying data, although not perfectly aligned with the moving average curve in Figure 3), coinciding closely with the clinical onset of the COVID-19 in Allegheny County. Because NRDM-ILI reflects ILI activity inferred from consumer purchasing behavior, the peak on March 13 is likely due to precautionary purchasing triggered by widespread media coverage and growing public concern about COVID-19. The ILI counts from ED data have a peak on January 4, 2020 (not exactly aligned with the moving average curve), which corresponds to increased influenza and RSV activity observed in our ED data in late 2019 and in January and February 2020. A corresponding increase is seen in the NRDM signal during that period as well (Figure 3, blue curve).

In addition, Figure 3 includes a red curve representing ED-derived counts of ILI and COVID-19 beginning in March 2020, corresponding to the initiation of COVID-19 testing in Allegheny County. As COVID-19 subsequently spread throughout the county, it became the dominant respiratory infectious agent, while influenza and RSV cases declined sharply. Consequently, the NRDM-ILI trend (blue curve) began to align more closely with the ED-derived counts that include COVID-19 (red curve), rather than the ED counts that do not (yellow curve).

Figure 4 contains 2 plots displaying the daily log-likelihood of the OTC data computed by NRDM Tracker to monitor for unusual data patterns that may indicate the outbreak of a novel disease. The top panel shows the daily log-likelihood. A prominent drop is visible on March 13, 2020. Interestingly, the top panel also shows sharp declines on other dates, such as July 2020, suggesting potential anomalies in purchasing patterns. The bottom panel depicts the normalized daily log-likelihood, computed as the geometric mean of category-level probabilities. This measure indicates whether there are changes in the relative proportions of product category purchases. The normalized likelihood curve exhibits seasonal patterns and a sustained decline in mid-2020, followed by gradual recovery through 2021. The start of the COVID-19 outbreak does not appear to be associated with a strong change in the relative purchase of product categories, and thus, the prominent drop in the red curve in Figure 4 was driven by spikes or drops in total purchase counts rather than a shift in product mix.

Using the Benjamini-Hochberg FDR procedure, we computed the adjusted P values for the daily counts of product sales for each of the 7 OTC categories. Figure 5 shows the time series of FDR-adjusted P values on a logarithmic scale, with the red dashed line indicating the significance threshold at P=.01. Sharp drops in P values—particularly around March 2020—are visible across all categories, reflecting significant deviations in product purchasing behavior during the early phase of the COVID-19 pandemic. Several other spikes and dips are also observed in late 2020 and late 2021, some of which may correspond to the resurgence of the COVID-19 pandemic due to evolving variants of SAR-CoV-2.

This paper investigates the potential for NRDM Tracker, using OTC data, as a tool for public health surveillance and early outbreak detection. By leveraging OTC product sales and applying a probabilistic modeling approach, we were able to estimate daily ILI activity and identify abnormal population-level purchasing patterns. The EM-based algorithm estimated conditional probabilities that distinguish ILI-related purchasing behavior from that associated with other health conditions. The resulting ILI estimates from NRDM data exhibited moderately strong correlation (r=0.66) with ILI-related ED visits, which provides support for the model’s alignment with measurable clinical rates of ILI disease in the region.

Our findings indicate that the NRDM-based ILI estimates, derived from OTC purchasing patterns, can capture public behavioral responses to disease outbreaks. These responses likely include purchases of cough and cold medications, antipyretics, and thermometers, even in the absence of corresponding clinical ILI encounters. Such dynamics could explain the divergence between the NRDM and ED curves in late 2019 and very early 2020, when NRDM estimates rose sharply while ED diagnoses for ILI remained comparatively stable. In general, while there is a close correspondence between the NRDM and ED curves in Figure 3, it is apparent that they are at times different, with an NRDM-derived estimate of increased ILI counts sometimes preceding the signal from ED data.

It is important to note that ED-based ILI surveillance, while widely used, may underestimate the true burden of ILI in the population [31], as many individuals with mild symptoms do not seek care in EDs. As a result, ED-derived ILI counts may primarily reflect more severe cases and may lag behind broader community-level disease activity. This limitation may help explain some of the observed differences between NRDM-based estimates and ED-based ILI signals. In particular, OTC purchasing behavior may capture earlier or milder symptom responses in the population, such as self-treatment before seeking professional medical care. A previous study [32] has noted temporal lags between community-level infection trends and subsequent increases in health care use, including ED visits and hospitalizations. In this context, the NRDM-based signals observed in this study may reflect earlier-stage or less severe disease activity, complementing ED-based surveillance. Together, these findings suggest that NRDM-based and ED-based estimates capture complementary aspects of disease activity, and that combining these signals may provide a more comprehensive assessment of population-level ILI than either source alone. Although fever is a defining symptom of ILI, antifever medications show a relatively weaker association with ILI activity compared with cough or cold medications in our analysis. This likely reflects differences in product use behavior rather than inconsistency with clinical definitions. Antifever medications are commonly used across a wide range of conditions and may already be available in households, reducing the need for purchase during illness episodes. In addition, OTC purchasing behavior may capture less severe cases of ILI, where individuals may not require or choose to take antipyretic medications, whereas ED-based surveillance is more likely to reflect more severe cases, potentially with higher fever. Similar patterns are observed for other general-purpose products, which tend to show weaker associations with ILI compared with more symptom-specific medications.

The NRDM Tracker system also derived daily likelihoods (in log form) of observed sales patterns as an overall measure of how well its model fit the OTC data. If the log-likelihood value on a given day is very low, relative to the values on previous days, it suggests the possibility of an outbreak. If so, those product categories with low P values are likely contributing to the low log-likelihood. Our results suggest that likelihood-based signals may provide an early indication of abnormal health-related activity. The sharp drop in log-likelihood observed in March 2020 aligns well with the initial COVID-19 outbreak in the United States. Importantly, our right-tailed P value framework—combined with multiple comparison correction—allowed for detection of statistically significant anomalies across multiple product categories. This analysis provides insights into the product sales that drive the detection of an anomaly.

Compared with previous studies that primarily examined associations between OTC sales and health care use [19,33] or evaluated NRDM for specific outbreak investigations [20,21,23,24], our study extends this body of work by developing a probabilistic surveillance framework capable of both estimating population-level ILI activity and identifying statistically significant anomalies that may indicate emerging or unexpected disease activity. Although previous studies have demonstrated that OTC purchasing data collected through NRDM can provide complementary syndromic signals of infectious disease activity [19,23,24], these studies have primarily focused on evaluating relationships with known clinical outcomes or detecting predefined outbreak events. In contrast, the NRDM Tracker presented in this paper uses probabilistic modeling of product category–specific purchasing patterns to continuously estimate ILI activity and detect statistically significant deviations from expected behavior, providing a framework that may support both seasonal disease monitoring and early identification of abnormal health-related signals potentially associated with emerging infectious disease activity.

There are several limitations of our approach. First, NRDM sales data do not directly record patient diagnoses, making it challenging to validate findings against gold-standard clinical data without assumptions. Second, behavioral factors may influence OTC purchasing patterns independently of underlying disease activity. For example, during the early phase of the COVID-19 pandemic, widespread public concern and media coverage likely contributed to “panic purchasing” of symptom-related products such as fever-reducing medications and thermometers, leading to transient spikes in sales that may not directly reflect true disease incidence. In addition, changes in supply chains, product availability, and public health advisories—both historically and prospectively—may alter purchasing behavior over time. Public statements or widely disseminated health communications, as for example, those related to the safety of specific medications, may also influence consumer behavior and health care–seeking patterns. These factors may introduce deviations in observed OTC signals that are not directly attributable to disease dynamics, potentially leading to false-positive signals of outbreaks. Finally, our modeling assumes stable conditional probabilities over time, which may not hold true under significant shifts in disease presentation or consumer behavior. This limitation could be attenuated by periodically updating the conditional probabilities.

Taken together, these limitations point to challenges that remain to be addressed in future work. Future research may also enhance the current work by incorporating additional data sources (eg, participatory surveillance systems [9,11,12], social media signals [34,35], and other digital epidemiologic data streams [36]) and by applying adaptive models that account for evolving behaviors and changing disease dynamics [37]. These directions are consistent with recent advances in AI [38] for integrating heterogeneous data sources to improve infectious disease surveillance and epidemic modeling [39].

In summary, our results provide support that NRDM-based data, when paired with probabilistic modeling and anomaly detection techniques, may provide an effective supplement to traditional disease monitoring systems. In particular, the observed alignment between NRDM-derived ILI estimates and ED-reported ILI cases, along with the detection of known public health events (eg, COVID-19 onset), support the system’s potential for early detection of seasonal disease activity and for identifying anomalies that may be associated with novel disease outbreaks. By identifying significant deviations in symptom-related OTC product sales, the system may offer timely public health insights that could support outbreak monitoring and response. Continued refinement and integration with other data sources and modeling frameworks may further strengthen the role of OTC data in the national biosurveillance ecosystem.