A forecast model for prevention of foodborne outbreaks of non-typhoidal salmonellosis

GLM, GLARMA and GSARIMA framework

Let Y be an RV related to the counting reports of infections by a pathogen of interest such as Salmonella enterica. We consider that statistical distribution of Y it belongs to the exponential family and that its PDF is (1)${\displaystyle f_Y\left(y;\vartheta ,\phi \right)=exp\left(\frac{y\vartheta -b\left(\vartheta \right)}{\phi }+c\left(y,\phi \right)\right),y\in R_Y,}$ where ϑ, φ are canonical and scale parameters, respectively, R_Y is the support of Y and b, c are specific functions whose characterize a member of the exponential family of statistical distributions.

This parametrization has the following properties for the description of mean and variance as the first and second derivatives of the canonical and scale parameters, where E(Y) = b′(ϑ) and Var(Y) = φb″(ϑ), respectively.

In GLM, a function denominated link function of the mean of Y is equal to a systemic component, then g(μ) = η. In turn, the mean of Y corresponds to the inverse link function of η, which is related to know values of r covariates x = (x₀, x₁, …, x_r)^⊤, with x₀ = 1, where (2)${\displaystyle \mu =\text{E}\left(Y\right)=g^{-1}\left(\eta \right)=g^{-1}\left({\mathit{x}}^{\top }\beta \right),}$ with β = (β₀, β₁, …, β_r)^⊤ being the regressors associated with x.

Now we will derive the expressions for the case of an RV as Y but ordered in a temporal sequence t, this RV indexed over time t, with t = 1, …, n is denominated Y_t. The conditional distribution of Y_t given the past data set (3)${\displaystyle {\mathit{H}}_t=\left\{{\mathit{x}}_1,\dots ,{\mathit{x}}_t,y_1,\dots ,y_{t-1}\right\},}$ which is also assumed to belong to the exponential family of statistical distributions. In this context, the past data conditional PDF similar to Eq. (1) is expressed as ${\displaystyle f_{Y_t|{\mathbf{H}}_t}\left(y_t;{\vartheta }_t,\phi \right)=exp\left(\frac{y_t{\vartheta }_t-b\left({\vartheta }_t\right)}{\phi }+c\left(y_t,\phi \right)\right),y_t\in R_{Y_t}.}$ Note that here the canonical parameter ϑ_t and values of covariates x_t depend on a temporal sequence, while the parameter φ,  remains constant and independent of the time sequence. In this conditions, we denote the conditional mean and variance of Y_t given H_t by μ_t = E(Y_t|H_t) = b′(ϑ_t) and Var(Y_t|H_t) = φb^′′(ϑ_t), for t = 1, …, n, respectively. g(μ_t) can be expressed as a GLARMA model of p and q orders. This model is denoted by GLARMA(p, q), where (4)${\displaystyle {\eta }_t=g\left({\mu }_t\right)={\mathit{x}}_t^{\top }\beta +\sum _{h=1}^p{\varphi }_h\left(g\left(y_{t-h}\right)-{\mathit{x}}_{t-h}^{\top }\beta \right)+\sum _{j=1}^q{\lambda }_j\left(g\left(y_{t-j}\right)-{\eta }_{t-j}\right).}$ In this model ϕ_h and λ_j correspond to the hth and jth components of an ARMA(p, q) model, related to the autoregressive and moving average components, respectively. β are the regressors, related with known values of r covariates, depending over time, denoted by x_t = (x₀, x_1t, …, x_rt)^⊤, with x₀ = 1. The link function η_t = g(μ_t) of the GLARMA model given in Eq. (4) can any as the identity function inverse or logarithmic (log) functions (allowing to consider non-linear associations). In this model the variance is assumed to be constant over time. In the case of the identity link function, we have that (5)${\displaystyle {\mu }_t={\mathit{x}}_t^{\top }\beta +\sum _{h=1}^p{\varphi }_j\left(y_{t-h}-{\mathit{x}}_{t-h}^{\top }\beta \right)+\sum _{j=1}^q{\lambda }_j\left(y_{t-j}-{\mu }_{t-j}\right).}$

The above models can be extended to GSARIMA (p, d, q) × (P, D, Q)_s analogues by including seasonality (S) and differentiation (D) components as follows: ${\displaystyle g\left({\mu }_t\right)=\Phi \left(L\right){\left(1-L\right)}^d{\left(1-L^s\right)}^D{\Phi }^P\left(L\right)\left({\mathit{x}}_t^{\top }\beta -g\left(y_t\right)\right)+g\left(y_t\right)-{\Lambda }^Q\left(L\right)\Lambda \left(L^s\right)\left(g\left(y_t\right)-g\left({\mu }_t\right)\right)+g\left(y_t\right)-g\left({\mu }_t\right),}$

where s is the length of the period (s =12 for monthly data with an annual cycle), ${\Phi }^P\left(L^s\right)=1-{\varphi }_1^P{L}^s-\cdots -{\varphi }_p^P{L}^{s^P},\Lambda {\left(L^s\right)}^Q=1-{\lambda }_1^Q{L}^s-\cdots -{\varphi }_q^Q{L}^{s^Q}.$

Density forecast

The density forecast (DF) technique aims to assess predictive performance outside the sample. This technique consists of dividing the original sample of data into a training set (2/3 of the first data) that is used to estimate the parameters, and then evaluate the performance outside the sample with the third rest of the data. That is, an out-of-sample test. These out-of-sample tests should be considered in the model validation process, see details in Rojas (2016). A probability integral transform (PIT) is the cumulative probability evaluated at the actual, realized value of the target variable. It measures the likelihood of observing a value less than the actual realized value, where the probability is measured by the DF. The PIT is uniform, independent and identically distributed if the density forecast is correctly specified.

GSARIMA Forecast model for prevention of foodborne outbreaks by non-typhoidal salmonellosis

We adapted our methodological framework from Maëlle, Dirk & Michael (2014). In this context, we considered an a-head time series based on fitted values of GSARIMA to alert infectious outbreaks, whose can be represented for multivariate forecast time series of counts. Denote the counts as y_it;  i = 1; ⋯; m; t = n + 1; ⋯; n + h, where n + h is the length of the forecast time series, whose begin at time n + 1 and m is the number of entities, e.g., geographical regions, hospitals or age groups, being monitored. In the context of a forecast model for future disease outbreaks, it is essential to detect future changes in the process occurring at an unknown time τ. As noted by Sonesson & Bock (2003), this change can be a step increase of the counts of future cases or a more gradual change. Based on the possibility of such change, for each future time t we want to differentiate between the two states in-control and out-of-control. At any timepoint t₀ ≥ n + 1 the available information -i.e., fitted values counts - is defined as ${\hat{\mathbf{y}}}_{t_0}=\left\{{\hat{y}}_t:t\le t_0\right\}.$ Detection is based on a statistic m(⋅) with a resulting alarm time $T_A=min\left\{t_0\ge 1:m\left({\hat{\mathbf{y}}}_{t_0}\right)>a\right\}$ where a is a known threshold. The functions for outbreaks detection use fitted values to estimate y_t0,  and compare it to the threshold a, above which, the current count can be considered as suspicious and, therefore, doomed as out-of-control. Then, based on Farrington et al. (1996) we designed a forecast model of outbreaks that uses the function of an algorithm called algo.farrington of survillance a package, in R software, see Höhle & Riebler (2005). R is a non-commercial and open source software for statistical analyses and plotting, which can be obtained from http://www.r-project.org. We modified this function that summarize the Farrington et al. (1996) algorithm, to take a range of h-step ahead forecast values of the number of counts of infection reports. This h-step ahead forecast values are obtained from the multivariate time series GSARIMA forecast model show in ‘Forecast with GSARIMA or GLARMA models’. For each time ahead, we use a GSARIMA simulated time series to set the number of counts in the same periods as a comparison. To obtain GSARIMA simulated time series we take command garmsim of gsarima package in R software, see Briët, Amerasinghe & Vounatsou (2013). This command requires an autoregressive representation obtained by arrep function of same package. The estimated parameters from the time series of the observed data, necessary to define the autoregressive representation, can be estimated by glarma function of the same name package in R software, see Dunsmuir (2015). Then, this is compared to the forecast values number of counts. If the forecast value is above a specific quantile of the prediction interval given by simulated time series, then an alarm is raised, see Algorithm 1.

Simulation study

The new methodology proposed in this paper is studied using a Montecarlo (MC) simulation study.

The GSARIMA forecast is simulated by using the autoregressive representation method, which is implemented in R, using a package named gsarima. This package contains methods for the generation of random numbers with univariate structure of time series, which is expandable to a multivariate response. This package has functions that allow the generation of generalized time series for a set of statistical counting distributions that belong to the exponential family. We focused on the negative binomial distribution (NBI), given its common use in this type of modeling, due to its properties of accumulation of counts.

We generated 5,000 scenarios of simulation, establishing different conditions given by a set of parameters in each scenario, in order to verify if the amount of alarms generated varies according to the configuration of the scenarios. For example, scenarios with a positive or negative trend must have values of d > 0 and / or D > 0, and positive or negative autoregressive coefficients, respectively. In turn, scenarios with coefficients of positive moving average reveal upward trends, and vice versa. To generate these 5,000 scenarios, we used the following indicators with the mentioned distribution, selected following Briët, Amerasinghe & Vounatsou (2013):

Statistical parameters of NBI statistical distribution

• dispersion parameter: {2,4}

• mean: {7,10}

GSARIMA parameters

• autoregressive parameter: p ∼ U(−1, 1)

• moving average parameter: q ∼ U(−1, 1)

• seasonal autoregressive parameter: P ∼ U(−1, 1)

• seasonal moving average parameter: Q ∼ U(−1, 1)

• integration of time series parameter: d = {2, 1, 0}

• integration of seasonal time series parameter: D = {2, 1, 0}

• covariate: x ∼ U(0, 1)

• regressor coefficient: β = 0.7

• intercept: β₀ = 1

Using these parameters, we simulated 5,000 predicted GSARIMA time series of 156 weeks, each with a frequency of 52 weeks/years (total of 3 years). In each scenario, the generated outbreak alarms were recorded following the Farrington algorithm adapted to our methodology, see Algorithm 1. To verify the differences in the quantities of alarms, in the different generated scenarios, we constructed the classification shown in Table 1, according to the parameters of our Monte Carlo simulation study.

In order to show the different options that we can explore in our simulation study, we showed only three examples of simulated GSARIMA time series (ts) of length 156 week (3 years), using the approach mentioned in ‘GSARIMA Forecast model for prevention of foodborne outbreaks by non-typhoidal salmonellosis’. The parameters of the simulated GSARIMA ts were: ts1 = {p > 0, q > 0, P < 0, Q < 0, d = 0, D = 2 }, ts2 = {p > 0, q > 0, P < 0, Q > 0, d = 0, D = 0 }, ts3 = {p < 0, q > 0, P > 0, Q > 0, d = 1, D = 0}. All these simulated GSARIMA ts were made considering high dispersion and high mean parameters, and can be examined in Fig. 1. On the abscissa axis, the time is shown in years divided weekly, where for example 2018.5 means the 26th week of the year 2018 (which marks the middle of the year 2018). In Fig. 2 the upperbound is shown as a dashed line (this band shows four-week aggregated data), while the alarms -timepoints where the upperbound has been exceeded- is shown as triangles. This date was obtained by the application of Algorithm 1 to generation of alarms for prevention of outbreaks in forecast GSARIMA ts1,ts2 and ts3.

We used the Kruskal-Wallis (KW) test to compare medians of indicators related to the grouped scenario denominations according to the denominations shown in Table 1. Results of comparative of mean/median, value of statistic KW and its p-value by grouped scenario denominations are showed in Table 2 .

Case study

In order to validate our model, we performed a case study based on data reports of Salmonella enterica serovar Enteriditis cases in Sydney, Australia (2014–2016). These data were weekly collected during three years (2014–2016) by the National Notifiable Diseases Surveillance System (NNDSS) from the Health Department of the Australian Government.

Figure 3 shows weekly time series of number of reports of Salmonella enterica serovar Enteriditis cases, whereas Fig. 4 presents (a) Autocorrelation function (ACF), (b) Partial Autocorrelation function (PACF) and (c) Negative binomial (NBII) Quantil-Quantil (QQ)-plots of the variable under study (weekly time series of number of reports of Salmonella enterica serovar Enteriditis cases), respectively. ACF and PACF show measures of the correlation between the observations of a time series separated by k lag time units, in this case, k = 1 week. In our case, we have autocorrelation in several weeks of lag, whose are shown on the lines that exceed the confidence intervals shown by the scored lines. QQ plot indicates the comparison between theoretical quantiles (given by the proposed statistical distribution, NBII) and empirical data of the variables under study, whose all appear within the confidence intervals shown by the scored lines. Note that these figures show that the probability distribution of these variables are NBII but not independent and identically distributed. The adjusted forecast to deal with a zero inflation is covered by the proposed statistical distribution, of a negative binomial type, which is capable of admitting data at zero, assigning to this value a probability distribution according to its historical and temporal frequency. To confirm that the original time series is stationary, we applied Augmented Dickey-Fuller Test, obtaining a statistics Dickey-Fuller = −4.0722, Lag order = 3, p-value = 0.0144 for alternative hypothesis of time series with stationarity.

In this study we explored the relationship between the number of reports of Salmonella enterica serovar Enteriditis cases and different covariates obtained for the same time periods and city mentioned above. The consumption of undercooked eggs is known to be a risk factor for non-typhoid salmonellosis (Martelli & Davies, 2012) and infectious outbreaks caused by Salmonella enterica serovar Enteritidis have been widely reported worldwide (Pijnacker et al., 2019; Jiang et al., 2020; Muvhali et al., 2017; Rizzo, 2006). Additionally, environmental factors like temperature and humidity were included in this study. High temperatures have been previously correlated to an increased incidence of Salmonella infections (Akil, Ahmad & Reddy, 2014; Kovats et al., 2004), while humidity has been positively associated to Salmonella infections (Kim et al., 2015).

In this study, the covariates corresponded to the climatic variables of mean maximum temperature per week (°C/week), mean 3pm relative humidity (%/week), as well as the weekly demand for eggs (units/person*week). The climatological data were collected from the Bureau of Meteorology of the Australian Government, while egg consumption data were collected from the Euromonitor Agency. Figure 5 shows scatterplots between number of reports of Salmonella enterica serovar Enteriditis as response variable and the mentioned covariates. Note that, in general, the relationships between reports of Salmonella enterica serovar Enteriditis as response variable and the mentioned covariates are linear and positives. The covariates seemed to have a symmetric distribution.

We assumed a GSARIMA(p, d, q) model, with p = {0, 1, 2}, d = 0 and q = {0, 1} using a NBI statistical distribution for the number of reports de Salmonella enterica serovar Enteriditis cases and considered the covariates described with different link functions for the mean response. In order to select the best GSARIMA model, we used AIC, BIC and GD, whose are reported on Table 3.

From Table 3, note that the smallest AIC, BIC and GD correspond to the GSARIMA(2.0,0) model with an identity link function, which is equivalent to a GLARMA model (2.0). In order to model the response variable of number of reports of Salmonella enterica serovar Enteriditis (Y), we propose the GLARMA model given by (8)${\displaystyle \text{E}\left(Y_j\right)={\mu }_j={\beta }_0+{\beta }_1{x}_j+{\varphi }_1\left(y_{j-1}-{\beta }_0-{\beta }_1{x}_{j-1}\right)+{\varphi }_2\left(y_{j-2}-{\beta }_0-{\beta }_1{x}_{j-2}\right),}$ where β₀ and β₁ are the regression coefficients, x_j is the value of the covariate vector X = x₁, x₂, x₃, where x₁ = Mean Maximal Temperature (° C/week), x₂ = Mean 3pm Relative Humidity (%/week), x₃ = weekly demand for eggs (units/person*week) and ϕ₁, ϕ₂ are the autoregressive coefficients. We fit the GLARMA model by using the command garmaFit. The maximum likelihood estimates of the parameters of the model given in Eq. (8), with approximate estimated standard errors in parenthesis, are: ${\widehat{\beta }}_0=-15.28\left(0.86\right)$, ${\widehat{\beta }}_1=\left\{{\widehat{\beta }}_{x_1}=0.25\left(0.019\right),{\widehat{\beta }}_{x_2}=0.19\left(0.02\right),{\widehat{\beta }}_{x_3}=0.95\left(0.58\right)\right\}$, ${\widehat{\varphi }}_1=-0.021\left(0.045\right)$, ${\widehat{\varphi }}_2=-0.014\left(0.045\right)$ and ${\left(\widehat{\text{Var}}\left(Y_t|{\mathbf{H}}_t\right)\right)}^{1∕2}={\left(\widehat{\phi }\left({\mu }_t\right)\right)}^{1∕2}=0.041\left(0.016\right)$. All coefficients are significant at 10%. The coefficient of determination (R²) of this model is 0.88. This indicates that the variables considered in the model (y_t−1, y_t−2, x₁, x₂, x₃) explain 88% of the variations of the response variable (y). In this model not all explanatory variables contribute the same in determining the variation of the response variable. For example, extracting the variable x₃ from the model R² remains almost unchanged (0.87), while extracting x₃ and x₂, R² decreases to 0.73, leaving only the lags y_t−1, y_t−2 explaining the response variable y, R² decreases to just 0.09. Then, prediction model can be expressed as ${\displaystyle {\widehat{\mu }}_{j+1}=-15.28+0.25x_{1_{j+1}}+0.19x_{2_{j+1}}+0.95x_{3_{j+1}}-0.021\left({\widehat{y}}_j-\left(-15.28+0.25x_{1_j}+0.19x_{2_j}+0.95x_{3_j}\right)\right)+0.041\left({\widehat{y}}_{j-1}-\left(-15.28+0.25x_{1_{j-1}}+0.19x_{2_{j-1}}+0.95x_{3_{j-1}}\right)\right).}$

To confirm the correct fit of this proposed GLARMA model, six plots were examined in Fig. 6: (a) Observed time series related to fixed effect of GLM estimation or GLARMA estimation;, (b) Disposition of Pearson Residuals in the time; (c) Histogram of Uniform PIT, (d) Histogram of Randomized Residuals normalized (randomized for a discrete response distribution); (e) Quantile-Quantile (QQ) Plots for randomized residuals of a fitted GLARMA object; and (f) Plot of the ACF of the residuals. These results indicated that is possible to apply a GLARMA model to original data, considering the covariates and logarithmic link function to the mean of response variable in each time of its realization. The selection of the best GLARMA model order is according Akaike criteria. Addiotonally, we also checked Box–Pierce test type Ljung to corroborate aleatory disposition to 2 lags with statistic χ² = 1.18, degree of freedom (df) = 2, and p-value = 0.582, and normal distributions of the randomized residuals of the model, with Shapiro–Wilk normality test obtained statistic W = 0.9874, p-value = 0.9321. The forecast PDF serves to simulate future scenarios of number of reports of Salmonella enterica serovar Enteriditis cases.

Following Algorithm 1 we can apply our model for prevention of future foodborne non-typhoidal salmonellosis outbreaks considered the above mentioned forecast covariates projected for the next year. These weeks results are shown in Fig. 7. Our model makes it possible to predict 3 alarms in the third quarter of next year, given the expected weather and egg consumption conditions. Remember that upperbound showed as a dashed line shows four-week aggregated data. In this case, the model does not forecast outbreaks.