The objective of this case study is to assess the risk of invasive listeriosis from consumption
of cold-smoked salmon in France.
The process of interest lays from the end of the production line in the factory, when the
cold-smoked salmon is vacuum-packed, to the consumption.
The data and the model are adapted to illustrate the use of `mc2d`: the results will not and
*should not* be interpreted as an assessment of the actual risk of listeriosis from consumption
of cold-smoked salmon.
Interested readers could refer to \cite{POUILLOT-2007} and \cite{POUILLOT-2009} for a complete
risk assessment on that issue.
The model will be developed in a first section, without considering variability or uncertainty
(deterministic model). Variability will then be introduced in a second section, and a last
section will consider variability and a part of the data uncertainty.
# The Model
In this section, no variability nor uncertainty is considered.
We assess the final level of *L.\ monocytogenes* in the product, the exposure and the risk of
invasive listeriosis for an "average" individual of the "healthy" French
population.^[It makes little sense, but it will help us introducing smoothly the model.]
During the logistic, the retail and the home step, a bacterial growth is modeled considering
*i*) the fluctuating temperature during the various stages and; *ii*) the bacterial competition
with the food flora.
We use the models developed and/or used in \cite{POUILLOT-2007}. The data are adapted from
\cite{POUILLOT-2007} and \cite{DELIGNETTE-2006}:
- The DMS model predicts the bacterial growth during a stage of duration $d$,
when the temperature is fluctuating, with an intra-stage average temperature $m_{T}$
and an intra-stage standard deviation of the temperature $s_{T}$. It is written:
\begin{equation}
N_{1}=
\min\left(
N_{0} +
\frac{\mu_{ref}}{\ln(10)} \times d \times
\frac{\left(s_{T}^{2}+\left(m_{T}-T_{min}\right)^{2}\right)}{\left(T_{ref}-T_{min}\right)^{2}},
\: N_{max}\right) \label{DMS}
\end{equation}
if $m_{T}>T_{min}$, with $N_{1}$ the $\log_{10}$ concentration of bacteria ($\log_{10}$ CFU/g)
at the end of the stage, $N_{0}$ the $\log_{10}$ concentration at the beginning,
$\mu_{ref}$ the specific growth rate (day$^{-1}$) at a reference temperature $T_{ref}$ (°C),
$T_{min}$ the minimal growth temperature (°C), and $N_{max}$ the maximum achievable
concentration ($\log_{10}$ CFU/g). If $m_{T} \leq T_{min}$, $N_{1}=N_{0}$.
- We will use $T_{ref}=25$°C. In this section $N_{max}=7.27\ \log_{10}$ CFU/g.
- The model for *L.\ monocytogenes* uses $\mu_{ref,Lm}=6.2\ \text{day}^{-1}$ and
$T_{min,Lm}=-2.9$°C.
- The same model is used for the food flora, with $\mu_{ref,ff}=4.1\ \text{day}^{-1}$ and
$T_{min,ff}=-4.5$°C.
- The growth model for bacterial competition considers the Jameson effect, i.e., the growth of
*L.\ monocytogenes* and the food flora are stopped as soon as one population reaches $N_{max}$.
In practice, one evaluates $d_{Lm}$ and $d_{ff}$, the time needed for *L.\ monocytogenes*
or the food flora to reach $N_{max}$, and models growth during an effective duration of
$\min(d,d_{Lm},d_{ff})$.
The time needed to reach $N_{max}$ is evaluated by inverting \eqref{DMS}:
$$
d_{\left(N_{1}=N_{max}\right)}=
\left(N_{max}-N_{0}\right)\times
\frac{\ln(10)}{\mu_{ref}} \times
\frac{\left(T_{ref}-T_{min}\right)^{2}}{\left(s_{T}^{2}+\left(m_{T}-T_{min}\right)^{2}\right)}
$$
The other assumptions are:
- A cold-smoked salmon package is homogeneously contaminated with *L.\ monocytogenes*
at the end of production at a level of $0.1$ CFU/g.
- The food flora level at the end of production is $10^{2.78}$ CFU/g.
- The time-temperature profile is:
- 1.1 days at an average temperature of 3.2°C (logistic step), with an intra-stage
SD of 2.1°C.
- 4.7 days at 5.5°C at retail, with intra-stage SD of 1.0°C.
- 4.3 days at 8.2°C in the consumer's home, with intra-stage SD of 2.0°C.
- A healthy, non-elderly, non-pregnant individual eats 35 g of this product.
- The individual dose-response model is a one-hit model
$\Pr(\text{Illness} \mid D)=1-(1-r)^{D}$
with $r=4.7\times10^{-14}$ for this healthy sub-population. The mean-population
dose-response (Poisson-distributed dose with mean $D$) is
$\Pr(\text{Illness} \mid D)=1-\exp(-r\times D)$.
The question is: *What is the risk for this "average" individual?*
``` r
Nmax <- 7.3
murefLm <- 6.2; TminLm <- -2.9
murefFF <- 4.1; TminFF <- -4.5
Lm0 <- log10(1); FF0 <- 2.78
d1 <- 1.1; mT1 <- 3.2; sdT1 <- 2.1
d2 <- 4.7; mT2 <- 5.5; sdT2 <- 1.0
d3 <- 4.3; mT3 <- 8.2; sdT3 <- 2.0
conso <- 35
r <- 4.7e-14
modGrowth <- function(duration, mTemp, sdTemp,
N0Lm, murefLm, TminLm,
N0FF, murefFF, TminFF,
Nmax, Tref = 25) {
N0Lm <- pmin(N0Lm, Nmax)
N0FF <- pmin(N0FF, Nmax)
dLm <- (Nmax - N0Lm) * log(10) / murefLm * (Tref - TminLm)^2 /
(sdTemp^2 + (mTemp - TminLm)^2)
dLm <- ifelse(mTemp < TminLm & N0Lm != Nmax, Inf, dLm)
dFF <- (Nmax - N0FF) * log(10) / murefFF * (Tref - TminFF)^2 /
(sdTemp^2 + (mTemp - TminFF)^2)
dFF <- ifelse(mTemp < TminFF & N0FF != Nmax, Inf, dFF)
realDuration <- pmin(duration, dLm, dFF)
xLm <- N0Lm + (mTemp > TminLm) * murefLm / log(10) *
(sdTemp^2 + (mTemp - TminLm)^2) / (Tref - TminLm)^2 * realDuration
xFF <- N0FF + (mTemp > TminFF) * murefFF / log(10) *
(sdTemp^2 + (mTemp - TminFF)^2) / (Tref - TminFF)^2 * realDuration
return(list(xLm = xLm, xFF = xFF))
}
x1 <- modGrowth(d1, mT1, sdT1, Lm0, murefLm, TminLm, FF0, murefFF, TminFF, Nmax)
x2 <- modGrowth(d2, mT2, sdT2, x1$xLm, murefLm, TminLm, x1$xFF, murefFF, TminFF, Nmax)
x3 <- modGrowth(d3, mT3, sdT3, x2$xLm, murefLm, TminLm, x2$xFF, murefFF, TminFF, Nmax)
x3## $xLm
## [1] 3.21
##
## $xFF
## [1] 5.35## [1] 1637## [1] 57281## [1] 2.69e-09modGrowth is a convenient function for the growth model.
Within this function dLm is the time needed for
L. monocytogenes to reach Nmax, dFF
for the food flora, and realDuration is the effective
duration of growth during the stage. Note that:
pmin (parallel minimum) is
used instead of min (which returns the minimum of
all values), and ifelse instead of
if.x1, x2 and x3 are the
bacterial concentrations at the end of the logistic, retail and home
steps, respectively.
We now specify variability distributions for some inputs, following and .
## [1] 10001For the initial contamination levels in L. monocytogenes, we
have 62 enumeration data from a representative sample of positive
packages of cold-smoked salmon: 43 samples < \(0.2\) CFU/g, 7 at \(0.2\) CFU/g, 4 at \(0.4\) CFU/g, 2 at \(0.6\) CFU/g, and values of 0.3, 1.0, 1.6,
2.4, 5.4 and 7.0 CFU/g . We use fitdistrplus to fit a
normal distribution to the \(\log_{10}\) values, accounting for
censoring. The initial concentrations are then modelled through a
truncated normal1 on \([-2,\infty)\
\log_{10}\) CFU/g.
For the food flora, we use the distribution from : \(N_{0ff} \sim N(2.78, 1.14)\).
dataC <- data.frame(
left = c(rep(NA, 43), rep(.2, 7), .3, rep(.4, 4), 1, 1.6, .6, .6, 2.4, 5.4, 7),
right = c(rep(0.2, 43), rep(.2, 7), .3, rep(.4, 4), 1, 1.6, .6, .6, 2.4, 5.4, 7)
)
fit <- fitdistcens(log10(dataC), "norm")
fit## Fitting of the distribution ' norm ' on censored data by maximum likelihood
## Parameters:
## estimate
## mean -1.117
## sd 0.764Lm0V <- mcstoc(rnorm, mean = fit$est["mean"], sd = fit$est["sd"], rtrunc = TRUE, linf = -2)
FF0V <- mcstoc(rnorm, mean = 2.78, sd = 1.14)By default, the type of variability modelled is "V"
(variability).
Distributions derived from :
The time-temperature profiles are modelled using Table \(\ref{tab:4}\) (adapted from ).2 We assume a shelf life of 28 days with \(d_{1}+d_{2}+d_{3}\leq28\).3
d1V <- mcstoc(rexp, rate = 1/1.1)
mT1V <- mcstoc(rnorm, mean = 3.2, sd = 2.2, rtrunc = TRUE, linf = -3, lsup = 25)
sdT1V <- sqrt(mcstoc(rgamma, shape = 1.16, scale = 4.61))
d2V <- mcstoc(rexp, rate = 1/4.7, rtrunc = TRUE, lsup = 28 - d1V)
mT2V <- mcstoc(rnorm, mean = 5.5, sd = 2.2, rtrunc = TRUE, linf = -3, lsup = 25)
sdT2V <- sqrt(mcstoc(rgamma, shape = 0.65, scale = 2.09))
d3V <- mcstoc(rexp, rate = 1/4.3, rtrunc = TRUE, lsup = 28 - (d1V + d2V))
mT3V <- mcstoc(rnorm, mean = 8.2, sd = 3.8, rtrunc = TRUE, linf = -3, lsup = 25)
sdT3V <- sqrt(mcstoc(rgamma, shape = 0.35, scale = 19.7))From observed data, a discrete empirical distribution is used : values \(V = \{10, 12, 19, 20, 30, 34, 40, 50, 60, 67.5, 80, 100, 250\}\) g, observed \(F = \{11, 1, 1, 29, 12, 1, 41, 4, 4, 1, 4, 1, 1\}\) times.
r <- mcdata(4.7e-14, type = "0")
x1V <- modGrowth(d1V, mT1V, sdT1V, Lm0V, murefLmV, TminLmV, FF0V, murefFFV, TminFFV, NmaxV)
x2V <- modGrowth(d2V, mT2V, sdT2V,
x1V$xLm, murefLmV, TminLmV, x1V$xFF, murefFFV, TminFFV, NmaxV)
x3V <- modGrowth(d3V, mT3V, sdT3V,
x2V$xLm, murefLmV, TminLmV, x2V$xFF, murefFFV, TminFFV, NmaxV)
contaV <- 10^x3V$xLm
expoV <- consoV * contaV
riskV <- 1 - exp(-r * expoV)
Lm1 <- mc(Lm0V, FF0V, NmaxV, murefLmV, TminLmV, murefFFV, TminFFV,
d1V, mT1V, sdT1V, d2V, mT2V, sdT2V, d3V, mT3V, sdT3V,
consoV, r, contaV, expoV, riskV)
Lm1## node mode nsv nsu nva variate min mean median max Nas type outm
## 1 Lm0V numeric 10001 1 1 1 -2.00e+00 -9.30e-01 -9.88e-01 1.76e+00 0 V each
## 2 FF0V numeric 10001 1 1 1 -1.28e+00 2.78e+00 2.78e+00 6.85e+00 0 V each
## 3 NmaxV numeric 10001 1 1 1 3.97e+00 7.26e+00 7.27e+00 1.06e+01 0 V each
## 4 murefLmV numeric 10001 1 1 1 2.85e+00 6.24e+00 6.25e+00 9.25e+00 0 V each
## 5 TminLmV numeric 10001 1 1 1 -1.01e+01 -2.83e+00 -2.85e+00 3.69e+00 0 V each
## 6 murefFFV numeric 10001 1 1 1 1.31e-02 4.19e+00 4.17e+00 1.13e+01 0 V each
## 7 TminFFV numeric 10001 1 1 1 -3.51e+01 -4.52e+00 -4.46e+00 2.66e+01 0 V each
## 8 d1V numeric 10001 1 1 1 5.36e-05 1.10e+00 7.69e-01 9.69e+00 0 V each
## 9 mT1V numeric 10001 1 1 1 -2.98e+00 3.20e+00 3.15e+00 1.14e+01 0 V each
## 10 sdT1V numeric 10001 1 1 1 3.37e-02 2.08e+00 1.96e+00 6.57e+00 0 V each
## 11 d2V numeric 10001 1 1 1 1.70e-03 4.69e+00 3.29e+00 2.67e+01 0 V each
## 12 mT2V numeric 10001 1 1 1 -2.53e+00 5.53e+00 5.50e+00 1.35e+01 0 V each
## 13 sdT2V numeric 10001 1 1 1 1.60e-03 9.66e-01 8.63e-01 4.76e+00 0 V each
## 14 d3V numeric 10001 1 1 1 3.84e-04 4.09e+00 2.89e+00 2.53e+01 0 V each
## 15 mT3V numeric 10001 1 1 1 -2.98e+00 8.24e+00 8.21e+00 2.26e+01 0 V each
## 16 sdT3V numeric 10001 1 1 1 1.31e-06 1.91e+00 1.41e+00 1.16e+01 0 V each
## 17 consoV numeric 10001 1 1 1 1.00e+01 3.55e+01 4.00e+01 2.50e+02 0 V each
## 18 r numeric 1 1 1 1 4.70e-14 4.70e-14 4.70e-14 4.70e-14 0 0 each
## 19 contaV numeric 10001 1 1 1 1.16e-02 4.09e+06 2.22e+01 5.04e+09 0 V each
## 20 expoV numeric 10001 1 1 1 1.70e-01 1.18e+08 6.85e+02 1.01e+11 0 V each
## 21 riskV numeric 10001 1 1 1 7.99e-15 5.55e-06 3.22e-11 4.73e-03 0 V eachsLm1 <- mc(contaV = Lm1$contaV, expoV = Lm1$expoV, riskV = Lm1$riskV)
summary(sLm1, probs = c(0, 0.5, 0.75, 0.95, 1))## contaV :
## mean sd Min 50% 75% 95% Max nsv Na's
## NoUnc 4092015 76520808 0.0116 22.2 867 1962680 5.04e+09 10001 0
##
## expoV :
## mean sd Min 50% 75% 95% Max nsv Na's
## NoUnc 1.18e+08 1.88e+09 0.17 685 27251 60077024 1.01e+11 10001 0
##
## riskV :
## mean sd Min 50% 75% 95% Max nsv Na's
## NoUnc 5.55e-06 8.82e-05 7.99e-15 3.22e-11 1.28e-09 2.82e-06 0.00473 10001 0Lm1 contains all parameters and outputs.
sLm1 extracts a short summary.
If 6.5% of cold-smoked salmon packages are contaminated, 49,090,000 French people are in the “non-susceptible” population, and they consume smoked salmon 6.4 times per year on average, the expected number of listeriosis cases is:
meanRisk <- mcapply(riskV, "var", mean)
expectedN <- round(0.065 * unmc(meanRisk) * 6.4 * 49090000)
expectedN## [1] 113We include both variability and uncertainty, focusing on the uncertainty in initial contamination, growth parameters, and prevalence.
The uncertainty around initial L. monocytogenes
contamination is modelled using a bootstrap procedure with
bootdistcens from fitdistrplus.
## [1] 101bootLm0 <- bootdistcens(fit, niter = ndunc())
MLm0 <- mcdata(bootLm0$est$mean, type = "U")
SLm0 <- mcdata(bootLm0$est$sd, type = "U")
Lm0VU <- mcstoc(rnorm, type = "VU", mean = MLm0, sd = SLm0, rtrunc = TRUE, linf = -2)For the food flora, uncertain hyperparameters \(M_{N0ff}\) and \(\sigma_{N0ff}\) govern the variable parameter \(N_{0ff}\) : \[ N_{0ff} \sim N(M_{N0ff},\sigma_{N0ff}), \quad M_{N0ff} \sim N(2.78, 0.265), \quad \ln(\sigma_{N0ff}) \sim N(0.114, 0.172) \]
Uncertainty around \(\mu_{ref,Lm}\), \(T_{min,Lm}\), \(\mu_{ref,ff}\), \(T_{min,ff}\) and \(N_{max}\) is modelled via hyperparameters .4
\[\begin{eqnarray*} \mu_{ref,Lm} & \sim & N(M_{\mu ref,Lm},\sigma_{\mu ref,Lm}), \quad M_{\mu ref,Lm} \sim \Gamma(69.7,\ 0.0896), \quad \ln(\sigma_{\mu ref,Lm}) \sim N(1.03,\ 0.191)\\ T_{min,Lm} & \sim & N(M_{Tmin,Lm},\sigma_{Tmin,Lm}), \quad M_{Tmin,Lm} \sim N(-2.86,\ 0.459), \quad \ln(\sigma_{Tmin,Lm}) \sim N(0.638,\ 0.208)\\ \mu_{ref,ff} & \sim & N(M_{\mu ref,ff},\sigma_{\mu ref,ff}), \quad M_{\mu ref,ff} \sim \Gamma(32.5,\ 0.127), \quad \ln(\sigma_{\mu ref,ff}) \sim N(-0.656,\ 0.221)\\ T_{min,ff} & \sim & N(M_{Tmin,ff},\sigma_{Tmin,ff}), \quad M_{Tmin,ff} \sim N(-4.52,\ 1.23), \quad \ln(\sigma_{Tmin,ff}) \sim N(2.00,\ 0.257)\\ N_{max} & \sim & N(M_{Nmax},\sigma_{Nmax}), \quad M_{Nmax} \sim N(7.27,\ 0.276), \quad \ln(\sigma_{Nmax}) \sim N(-0.172,\ 0.218) \end{eqnarray*}\]
with \(\mu_{ref}>0\) and \(T_{min}<25\).
MmurefLm <- mcstoc(rgamma, type = "U", shape = 69.7, scale = 0.0896)
SmurefLm <- mcstoc(rlnorm, type = "U", meanlog = 1.03, sdlog = 0.191)
murefLmVU <- mcstoc(rnorm, type = "VU", mean = MmurefLm, sd = SmurefLm, rtrunc = TRUE, linf = 0)
MTminLm <- mcstoc(rnorm, type = "U", mean = -2.86, sd = 0.459)
STminLm <- mcstoc(rlnorm, type = "U", meanlog = 0.638, sdlog = 0.208)
TminLmVU <- mcstoc(rnorm, type = "VU", mean = MTminLm, sd = STminLm, rtrunc = TRUE, lsup = 25)
MmurefFF <- mcstoc(rgamma, type = "U", shape = 32.5, scale = 0.127)
SmurefFF <- mcstoc(rlnorm, type = "U", meanlog = -0.656, sdlog = 0.221)
murefFFVU <- mcstoc(rnorm, type = "VU", mean = MmurefFF, sd = SmurefFF, rtrunc = TRUE, linf = 0)
MTminFF <- mcstoc(rnorm, type = "U", mean = -4.52, sd = 1.23)
STminFF <- mcstoc(rlnorm, type = "U", meanlog = 2.00, sdlog = 0.257)
TminFFVU <- mcstoc(rnorm, type = "VU", mean = MTminFF, sd = STminFF, rtrunc = TRUE, lsup = 25)
MNmax <- mcstoc(rnorm, type = "U", mean = 7.27, sd = 0.276)
SNmax <- mcstoc(rlnorm, type = "U", meanlog = -0.172, sdlog = 0.218)
NmaxVU <- mcstoc(rnorm, type = "VU", mean = MNmax, sd = SNmax)x1VU <- modGrowth(d1V, mT1V, sdT1V,
Lm0VU, murefLmVU, TminLmVU, FF0VU, murefFFVU, TminFFVU, NmaxVU)
x2VU <- modGrowth(d2V, mT2V, sdT2V,
x1VU$xLm, murefLmVU, TminLmVU, x1VU$xFF, murefFFVU, TminFFVU, NmaxVU)
x3VU <- modGrowth(d3V, mT3V, sdT3V,
x2VU$xLm, murefLmVU, TminLmVU, x2VU$xFF, murefFFVU, TminFFVU, NmaxVU)
contaVU <- 10^x3VU$xLm
expoVU <- consoV * contaVU
riskVU <- 1 - exp(-r * expoVU)
Lm2 <- mc(Lm0VU, FF0VU, NmaxVU, murefLmVU, TminLmVU, murefFFVU, TminFFVU,
d1V, mT1V, sdT1V, d2V, mT2V, sdT2V, d3V, mT3V, sdT3V,
consoV, r, contaVU, expoVU, riskVU)
Lm2## node mode nsv nsu nva variate min mean median max Nas type outm
## 1 Lm0VU numeric 10001 101 1 1 -2.00e+00 -9.29e-01 -9.90e-01 4.71e+00 0 VU each
## 2 FF0VU numeric 10001 101 1 1 -3.90e+00 2.81e+00 2.81e+00 9.30e+00 0 VU each
## 3 NmaxVU numeric 10001 101 1 1 1.02e+00 7.30e+00 7.31e+00 1.32e+01 0 VU each
## 4 murefLmVU numeric 10001 101 1 1 6.55e-05 6.48e+00 6.40e+00 2.43e+01 0 VU each
## 5 TminLmVU numeric 10001 101 1 1 -1.58e+01 -2.86e+00 -2.86e+00 8.90e+00 0 VU each
## 6 murefFFVU numeric 10001 101 1 1 2.58e-01 4.15e+00 4.12e+00 9.05e+00 0 VU each
## 7 TminFFVU numeric 10001 101 1 1 -5.67e+01 -4.51e+00 -4.51e+00 2.50e+01 0 VU each
## 8 d1V numeric 10001 1 1 1 5.36e-05 1.10e+00 7.69e-01 9.69e+00 0 V each
## 9 mT1V numeric 10001 1 1 1 -2.98e+00 3.20e+00 3.15e+00 1.14e+01 0 V each
## 10 sdT1V numeric 10001 1 1 1 3.37e-02 2.08e+00 1.96e+00 6.57e+00 0 V each
## 11 d2V numeric 10001 1 1 1 1.70e-03 4.69e+00 3.29e+00 2.67e+01 0 V each
## 12 mT2V numeric 10001 1 1 1 -2.53e+00 5.53e+00 5.50e+00 1.35e+01 0 V each
## 13 sdT2V numeric 10001 1 1 1 1.60e-03 9.66e-01 8.63e-01 4.76e+00 0 V each
## 14 d3V numeric 10001 1 1 1 3.84e-04 4.09e+00 2.89e+00 2.53e+01 0 V each
## 15 mT3V numeric 10001 1 1 1 -2.98e+00 8.24e+00 8.21e+00 2.26e+01 0 V each
## 16 sdT3V numeric 10001 1 1 1 1.31e-06 1.91e+00 1.41e+00 1.16e+01 0 V each
## 17 consoV numeric 10001 1 1 1 1.00e+01 3.55e+01 4.00e+01 2.50e+02 0 V each
## 18 r numeric 1 1 1 1 4.70e-14 4.70e-14 4.70e-14 4.70e-14 0 0 each
## 19 contaVU numeric 10001 101 1 1 1.00e-02 1.00e+07 1.90e+01 4.20e+11 0 VU each
## 20 expoVU numeric 10001 101 1 1 1.00e-01 3.70e+08 5.79e+02 2.10e+13 0 VU each
## 21 riskVU numeric 10001 101 1 1 4.66e-15 1.67e-05 2.72e-11 6.27e-01 0 VU eachsLm2 <- mc(contaVU = Lm2$contaVU, expoVU = Lm2$expoVU, riskVU = Lm2$riskVU)
summary(sLm2, probs = c(0, 0.5, 0.75, 0.95, 1))## contaVU :
## mean sd Min 50% 75% 95% Max nsv Na's
## median 6299495 9.57e+07 0.0114 18.36 1244.5 5485146 6.12e+09 10001 0
## mean 10025362 2.54e+08 0.0120 26.79 2893.2 8564635 2.02e+10 10001 0
## 2.5% 1486775 2.06e+07 0.0101 3.86 97.4 202346 8.85e+08 10001 0
## 97.5% 53421962 1.94e+09 0.0174 82.67 14469.4 31267142 1.62e+11 10001 0
##
## expoVU :
## mean sd Min 50% 75% 95% Max nsv Na's
## median 2.35e+08 4.07e+09 0.153 566 37475 1.64e+08 2.59e+11 10001 0
## mean 3.70e+08 1.09e+10 0.162 814 90618 2.66e+08 9.07e+11 10001 0
## 2.5% 5.39e+07 7.65e+08 0.105 117 2946 6.30e+06 2.85e+10 10001 0
## 97.5% 1.96e+09 8.21e+10 0.233 2505 491756 1.00e+09 6.68e+12 10001 0
##
## riskVU :
## mean sd Min 50% 75% 95% Max nsv Na's
## median 1.09e-05 0.000191 7.22e-15 2.66e-11 1.76e-09 7.69e-06 0.01209 10001 0
## mean 1.67e-05 0.000452 7.60e-15 3.83e-11 4.26e-09 1.25e-05 0.03642 10001 0
## 2.5% 2.53e-06 0.000036 4.94e-15 5.48e-12 1.38e-10 2.96e-07 0.00134 10001 0
## 97.5% 8.22e-05 0.003279 1.09e-14 1.18e-10 2.31e-08 4.71e-05 0.26748 10001 0The summary provides the mean, SD, minimum, median … and a 95% credible interval. The point estimate is the median of the 101 values in the uncertainty dimension; the credible interval runs between the 2.5th and 97.5th percentiles of the uncertainty dimension.
meanRiskU <- mcapply(riskVU, "var", mean)
expectedNU <- round(prevU * meanRiskU * 6.4 * 49090000)
summary(expectedNU)## node :
## NoVar
## median 223.0
## mean 342.3
## 2.5% 48.5
## 97.5% 1586.5This estimate reflects uncertainty from initial contamination, bacterial growth parameters, and sampling of positive packages. Many other uncertainties (notably around the dose-response model) are not considered here; see for a complete analysis.
The tornado chart in the variability dimension suggests a large impact from the growth rate of L. monocytogenes, the storage duration during the consumer step, and the initial contamination level. The tornado in the uncertainty dimension highlights the impact of uncertainty around \(N_{max}\) on the mean risk.
## Spearman's rho statistic
## Output: riskVU
## $riskVU
## Lm0VU FF0VU NmaxVU murefLmVU TminLmVU murefFFVU TminFFVU d1V mT1V sdT1V d2V
## median 0.301 -0.0796 0.0690 0.452 -0.241 -0.03035 0.115 0.0437 0.0329 0.00669 0.281
## mean 0.303 -0.0845 0.0726 0.446 -0.240 -0.03230 0.123 0.0433 0.0332 0.00628 0.282
## 2.5% 0.191 -0.1568 0.0301 0.324 -0.337 -0.06209 0.033 0.0226 0.0181 -0.00653 0.207
## 97.5% 0.402 -0.0286 0.1377 0.557 -0.159 -0.00414 0.229 0.0613 0.0482 0.02273 0.347
## mT2V sdT2V d3V mT3V sdT3V consoV contaVU expoVU
## median 0.158 0.00511 0.409 0.260 0.0292 0.125 0.991 1
## mean 0.160 0.00539 0.411 0.261 0.0288 0.126 0.991 1
## 2.5% 0.129 -0.00727 0.319 0.212 0.0146 0.101 0.986 1
## 97.5% 0.185 0.01854 0.495 0.308 0.0446 0.155 0.994 1## Tornado on uncertainty
## Spearman's rho statistic
## Output: riskVU
## $riskVU
## mean Lm0VU sd Lm0VU 97.5% Lm0VU mean FF0VU sd FF0VU 97.5% FF0VU mean NmaxVU sd NmaxVU
## mean riskVU -0.0087 -0.0347 -0.0527 -0.0404 -0.0322 -0.08867 0.401 0.550
## sd riskVU -0.0572 -0.0808 -0.0985 -0.0366 -0.0909 -0.12235 0.162 0.684
## 97.5% riskVU 0.0507 -0.0159 -0.0191 -0.0331 0.0659 -0.00686 0.545 -0.075
## 97.5% NmaxVU mean murefLmVU sd murefLmVU 97.5% murefLmVU mean TminLmVU sd TminLmVU
## mean riskVU 0.823 0.415 0.348 0.493 0.0298 0.0923
## sd riskVU 0.755 0.235 0.303 0.353 0.1061 0.0785
## 97.5% riskVU 0.346 0.667 0.358 0.625 -0.1985 0.0646
## 97.5% TminLmVU mean murefFFVU sd murefFFVU 97.5% murefFFVU mean TminFFVU sd TminFFVU
## mean riskVU 0.0846 -0.0801 -0.0872 -0.102 0.1059 -0.0897
## sd riskVU 0.1139 -0.0135 -0.0120 -0.011 0.0843 -0.1136
## 97.5% riskVU -0.0470 -0.3539 -0.2595 -0.421 0.1261 -0.0200
## 97.5% TminFFVU mean contaVU sd contaVU 97.5% contaVU mean expoVU sd expoVU
## mean riskVU -0.0331 0.988 0.892 0.636 1.000 0.896
## sd riskVU -0.0644 0.857 0.939 0.325 0.897 1.000
## 97.5% riskVU 0.0278 0.649 0.317 0.996 0.637 0.330
## 97.5% expoVU
## mean riskVU 0.638
## sd riskVU 0.330
## 97.5% riskVU 1.000Tornado chart for the L. monocytogenes example (Variability).
Tornado chart for the L. monocytogenes example (Uncertainty).
As a conclusion, this example illustrates how predictive growth
models can be implemented within mc2d…
so that at least one CFU is included in one 100 g package↩︎
\(\Gamma\) denotes the Gamma distribution parameterized as \(\Gamma(shape,\ scale)\). Exponential(\(x\)) denotes the exponential distribution with mean \(x\).↩︎
See the code for how to model this using truncated distributions.↩︎
Note: there was a typo in that led to an error in : the standard-error for \(\ln(\sigma_{\mu ref,Lm})\) is \(1.03\), not \(-1.03\). We use the correct value here.↩︎