Variable selection

A decision to use outbreak or illness counts in an estimation strategy depends on two factors, with the first being whether one of the variables is better suited to the modelling of trends across time. The second is to determine whether there are significant differences between the outbreak counts and illness counts that might lead to contradictory conclusions.

The estimation of temporal trends is dependent on the strength of the underlying temporal pattern compared to the degree of variability in the data, which can be defined by the signal-to-noise ratio given by

$$SNR_x = {{E( x) } \over {\sqrt {Var( x) } }}.$$

The variable with the larger signal-to-noise ratio will be better suited for the estimation of trends.

To investigate the SNR statistic for outbreak counts, assume the annual number of outbreaks can be modelled as O _t ~ Poisson(λ), where t is the indicator for year t. Using this model, the signal-to-noise ratio for outbreak counts is

$$SNR_O = \displaystyle{{E( O_t) } \over {\sqrt {Var( O_t) } }} = \displaystyle{\lambda \over {\sqrt \lambda }} = \sqrt \lambda .$$

For each outbreak j in year t, let S _tj be the number of illnesses and assume a general discrete distribution with mean and variance terms defined by S _tj ~ (μ, σ ²). Then the number of illnesses in year t is defined by the random sum $I_t = \sum\limits_{j = 1}^{O_t} {S_{tj}}$. The properties of a variable that is a random sum are such that

$$E( I_t) = E( O_t) E( S_t) = \lambda \mu $$

and

$$Var( I_t) = E( O_t) Var( S_t) + E( S_t) ^2Var( O_t) = \lambda \sigma ^2 + \mu ^2\lambda , \;$$

which yields

$$SNR_I = \displaystyle{{\mu \sqrt \lambda } \over {\sqrt {\mu ^2 + \sigma ^2} }}.$$

Given the signal-to-noise ratio values for outbreak vs. illness counts, we can determine under what conditions SNR _I ≥ SNR _O. Investigating this inequality yields

$$\eqalign{& \displaystyle{{\mu \sqrt \lambda } \over {\sqrt {\mu ^2 + \sigma ^2} }} \ge \sqrt \lambda \cr & \displaystyle{\mu \over {\sqrt {\mu ^2 + \sigma ^2} }} \ge 1, \;} $$

with the condition of equality satisfied only for σ ² = 0. It might be reasonable to use illness counts if σ ² were small relative to μ ², but a previous analysis of outbreak illness counts found that they are highly variable and remain heavily right-skewed even after applying a log transformation [Reference Ebel32].

We conclude that significant trends are more likely to be observed using the outbreak counts variable O _t because it necessarily has a larger signal-to-noise ratio. The effect of the higher signal-to-noise ratio of O _t can be quantified by comparing the widths of the confidence intervals relative to the estimated trend model, described below, that was fitted to O _t and I _t. For example, if $\hat{O}_t$ is the fitted value at time t, the average relative interval width statistic is

$$W_O = \displaystyle{1 \over T}\sum\limits_{t = 1}^T {( {\hat{O}}_t + t_\ast{\times} s.e.( {\hat{O}}_t) -( {\hat{O}}_t-t_\ast{\times} s.e.( {\hat{O}}_t) ) /{\hat{O}}_t} , \;$$

where $s.e.( \hat{O}_t)$ is the estimated standard error and $t_\ast$ is the t-statistic derived from the fitted model described below.

The second step in the variable selection process is to determine what, if any, differences can be detected between the annual occurrence of outbreaks and illnesses for each commodity. To start, we assessed the correlation between the number of outbreaks and illnesses [Reference Hollander, Wolfe and Chicken33]. Next, an analysis of variance model is used to determine if significant differences exist in the number of illnesses per outbreaks for any of the commodity groups.

Trends

The goal of this study is to determine if significant trends exist in the proportion of outbreaks associated with any of the commodities. If the total number of outbreaks or illnesses is constant across time, then changes in the annual number of outbreaks or illness for a commodity group, rather than the ratio of outbreaks in a commodity group relative to the total number of outbreaks, can be modelled using a single methodological approach.

The simplest test for possible trends is a Mann–Kendall test [Reference Mann34] applied to each of the commodity groups to assess if a significant monotonic trend was present in the data. To test for trends that are not necessarily monotonic, a general additive model was fitted to the data. For the outbreak counts, the model used for outbreaks is of the form

$${\rm log}( O_t) = \beta _0 + f( Year_t) {\rm} + \varepsilon _t, \;$$

with t = 1998, 1999, …, 2016, 2017 and f(Year _t) being a penalised B-spline regression function with a Poisson link function, as described in previous food safety applications [Reference Powell35–Reference Wood37]. The Bayesian information criterion is used to select the smoothing parameters of the best fitting model [Reference Neath and Cavanaugh38]. This model employs a visual interpretation for significant trends using a horizontal line test. A significant trend is indicated in any situation where a horizontal line, drawn across the graph, intersects both the lower and upper bounds of the interval. For this study, the lower and upper bounds are defined as the 2.5th and 97.5th percentiles. Therefore, a horizontal line that fails to intersect these boundaries results in a conclusion that we cannot reject the null hypothesis of no trend at an α = 0.05 level of significance. Estimation of this model's parameters was performed using the R software package Mixed GAM Computation Vehicle (mgcv) [39, Reference Wood40]. This model was fitted to summarise the annual totals O _t and I _t. The confidence intervals determine the average relative confidence interval width statistics (W _O, W _I).

Interpretation of trends in the proportion of outbreaks or illnesses associated with a specific commodity is not straightforward because a change in one commodity necessarily affects the proportions for all remaining commodities. The term spurious correlation was used to describe this phenomenon at the end of the 19th century [Reference Pearson31]. This effect can be demonstrated by considering a simplified example. Suppose that outbreaks in year t consist of only two groups, with the groups being outbreaks associated with commodities derived from animals (o _ta) and the other group being outbreaks from plant-derived foods (o _tp). Furthermore, assume the number of animal-associated outbreaks is increasing linearly (o _ta = α _a + φ_at, with α _a, φ_a > 0) with time, while the number of plant-associated outbreaks remains constant (o _tp = α _p). The proportion of plant-associated outbreaks is

$$p_p( t) = \displaystyle{{o_{tp}} \over {o_{tp} + o_{ta}}}.$$

The effect of the increase in animal-associated outbreaks on the proportion of plant-associated outbreaks can be determined by considering the first derivative, given by

(1)$$\displaystyle{{\partial p_p( t) } \over {\partial t}} = \displaystyle{{-\alpha _p\varphi _a} \over {{( \alpha _p + \alpha _a + \varphi _at) }^2}}.$$

This example demonstrates that a significant change in the annual proportion of outbreaks associated with a specific commodity category is necessarily offset by changes in the opposite direction in the remaining categories. Furthermore, this relationship demonstrates why changes in the proportion of outbreaks for any single commodity cannot be considered in the absence of the remaining commodities.

The trend model uses the ratio of the number of outbreaks in year t attributed to food categoryd, denoted o _t,d, d = 1, …D, with respect to the total number of outbreaks from the year, denoted O _t. The term

$$\hat{\,p}_{t, d} = \displaystyle{{o_{t, d}} \over {O_t}}$$

is the proportion of outbreaks in category d at time t. A similar model can be constructed for outbreak illnesses. The constraint that the d proportions add to one implies that they are not independent of one another (i.e. the proportions are negatively correlated because if the proportion for one category increases, this increase must be offset by decreases in other categories). As a result, these proportions cannot be analysed using statistical methods that are intended for unconstrained data [Reference Pearson31, Reference Aitchison41]. The solution employed is to use the additive log-ratio transformation technique used in compositional data analyses [Reference Faes42, Reference Van den Boogaart and Tolosana-Delgado43], which reparametrises the outbreak count variable as

$$\log \left({\displaystyle{{o_{t, d}} \over {o_{t, D}}}} \right), \;\quad d = 1, \;\ldots D-1.$$

The data were fitted to a multivariate linear model of the form

$$\log \left({\displaystyle{{o_{t, d}} \over {o_{t, D}}}} \right) = \alpha _{1d} + \alpha _{2d}t + \varepsilon _{td}$$

and tested for significance of the α _1d and α _2d, where the vector of error terms is

$$\varepsilon _t\sim MVN( 0, \;\;\Sigma ) , \;$$

with Σ being the D − 1 × D − 1 covariance matrix. A quadratic model was also tested for each commodity, though none were significant. A similar model could be constructed using illness counts, but these results are not presented because the signal-to-noise results from above are still applicable.

In cases where no outbreaks occurred for commodity group d in year t, the fraction o _td/o _tD was set to 1% of the fraction of outbreaks for that commodity. Parameter estimates are obtained using multivariate linear regression and the trend model for each category d is obtained by back-transforming the model predictions. Commodity group D is a reference group, and it consists of all outbreaks not being considered as part of the specific analysis. The trend for commodity group D is determined from the constraint that the sum of the proportions is equal to one. The results of the multivariate linear model are used to determine which commodities demonstrated significant changes in the proportion of outbreaks using a framework similar to one developed previously for assessing trends in antimicrobial resistance [Reference Faes42].

Consumption analysis

The goal of the meat consumption analysis is to provide basic estimates of the change in the number of servings between 1998 and 2017 for the beef, chicken, pork and turkey commodity groups. The change is estimated by fitting a simple linear model of the form

$$\eqalign{& {\hat{N}}_{d.t} = {\hat{C}}_{t, d} \times K \cr & {\hat{C}}_{d, t} = \gamma _{1d} + \gamma _{2d}t, \;\quad t = 1, \;\ldots , \;20, \;} $$

where $\hat{N}_{d, t}$ is the number of servings per year (or capita) for commodity group d, $\hat{C}_{d, t}$ is the annual consumption (kg) estimate per capita for commodity group d and K is a constant that converts consumption data to servings per year (or capita). For example, if it is assumed a typical serving size of 100 g applies to all commodities, then K = 10 servings per kilogram consumed per capita per year.

The proportional change in annual number of servings is estimated as

$$\Delta N_d = \displaystyle{{( {{\hat{N}}_{d, 20}-{\hat{N}}_{d, 1}} ) } \over {{\hat{N}}_{d, 1}}}.$$

Assuming a constant serving size, an improved understanding of the factors affecting changes in the estimated proportion of outbreaks for each commodity can be derived by considering the relationship for any year t

(2)$$I_{d, t} = I_{total, t} \times p_{d, t} = N_{d, t} \times P_{d, t}( ill) , \;$$

where I _d,t is the number of illnesses from the consumption of commodity d, I _total,t is the total number of foodborne illnesses from all commodities, p _d.t is the proportion of outbreaks from commodity d (assumed to be equivalent to the true attribution fraction of illness cases) and P _d.t(ill) is the probability of illness per serving. The total number of domestically acquired illnesses per year can be estimated from the reported case rate per 100 000 [44] multiplied by the total population of the United States [45], adjusted for foreign travel and non-foodborne cases [Reference Scallan46]. The P _d.t(ill) term can be decomposed to describe multiple risk factors [Reference Williams, Ebel and Vose47], such as dose concentration and difference in virulence between serotypes. Efforts undertaken by government or industry to reduce contamination would also be captured by changes in P _d.t(ill).

By assuming serving size is constant across time and commodities, an empirical measure of the change in the probability of illness per serving for commodity group d is estimated as

$$\Delta P_d( ill) = \displaystyle{{( {\Delta {\hat{\,p}}_d + 1} ) ( {\Delta {\hat{I}}_{total, t} + 1} ) } \over {( {\Delta {\hat{N}}_d + 1} ) }}-1.$$

As described for $\Delta \hat{N}_d$, the other changes are measured based on the linear difference between 2017 and 1998. This estimate of the change in the probability of illness per serving is intuitive because its numerator estimates the change in annual illness counts (as the product of the attribution fraction and total illnesses) and its denominator estimates the change in total servings. Nevertheless, a linear model fitted to the case rates of salmonellosis from 1998 through 2017, based on culture-confirmed test results [44] has no significant slope term (p = 0.156).