A Bootstrap Prediction Confidence Band for QMRA Beta-Poisson Dose-Response Models

Gang Xie

Easy Links

Research Article Open Access

A Bootstrap Prediction Confidence Band for QMRA Beta-Poisson Dose-Response Models

Gang Xie^1,2*

¹School of Medicine, Griffith University, Queensland, Australia
²Faculty of Science, Health, Education and Engineering, University of the Sunshine Coast, Queensland, Australia

*Corresponding authors address: Gang Xie, Adjunct Research Fellow, School of Medicine, Griffith University, Queensland 4111, Australia, Tel.: +61-7-37357485; Fax: +61-7-37355318; E-mail: g.xie@griffith.edu.au

Received: May 05, 2017; Accepted: May 23, 2017; Published: May 26, 2017

Citation: Gang Xie (2017) A Bootstrap Prediction Confidence Band for QMRA Beta-Poisson Dose-Response Models. Int Struct Comput Biol 1(1): 1-11.

Abstract

Let

P_{I} (d)

denote the probability of infection at a given mean dose d . In the quantitative microbial risk assessment (QMRA) framework, the beta-Poisson dose-response model

P_{I} (d) = 1 -_{1} F_{1} (α, α + β, - d)

where

_{1} F_{1} (., ., .)

denotes the Kummer confluent hypergeometric function and α,β are the model parameters, remains the most popular plausible dose response model in practice. One commonly accepted way of constructing the confidence band about the dose-response curve in QMRA literature is to follow a bootstrap procedure based on the maximum likelihood estimates of the model parameters α and β. Here, it is shown that this bootstrap confidence bands reported in the literature represent the confidence intervals for the mean value of the probability of adverse effect, (i.e.,

P_{I} (d)

which may represent the probability of infection at mean dose levels), not the confidence intervals for prediction. Therefore, the existing literature bootstrap (95%) confidence bands normally contain far less than 95% data points. In this study, a sample-size-dependent bootstrap (95%) confidence band for prediction is proposed which can normally contain 95% of the data points. Comparisons between the existing and the newly proposed confidence band algorithms were made through a comprehensive Monte Carlo simulation study by applying both algorithms to four well studied experimental dose-response datasets. The comparison results showed that the confidence band for prediction is a better representation of the degree of uncertainty of an estimated dose-response relation. The upper limit of the estimated prediction confidence band could be used as a better (i.e., more conservative but sensible) estimate for the worst case scenario in risk assessment.

Key words: QMRA; beta-Poisson dose-response model; bootstrap; confidence intervals for prediction; Monte Carol simulation.

Introduction

QMRA [1, 2] has provided a valuable alternative quantitative framework for analysis of adverse health outcomes associated with the environmental exposures to pathogenic organisms [3, 4] and is widely used by researchers to characterise microbial risks associated with food, water, and wastewater use in agriculture (e.g., [1, 3, 5]). The core part of the QMRA framework is the dose-response analysis which models the mathematical characterization of the relationship between the dose administered and the probability of adverse effect (typically, the probability of infection) in the exposed population [1, 2, 6]. Among different microbial dose-response models proposed in the literature, the beta-Poisson models remain the most popular plausible models in practice [1, 2, 7].The beta-Poisson dose-response model and its special case, the exponential doseresponse model have been well studied and widely employed to characterize infectivity of various viral, bacterial, and protozoan pathogens since the 1980s [1, 2, 7, 8]. Because a beta-Poisson model has two parameters, a confidence band about the doseresponse curve cannot be calculated directly from the confidence intervals of the parameter estimates. It is well accepted in QMRA practice to employ a parametric bootstrap algorithm [9] to construct a confidence band for the dose-response curve[1, 2, 8, 10].The upper limit of the estimated confidence band is used as the worst case scenario in risk assessment.

Here, it is shown that the (95%) bootstrap confidence bands reported in microbial risk literature normally contain far less than 95% of the data points because they are not confidence bands for prediction. This implies that the existing literature confidence band may not be a good representation of the degree of uncertainty of an estimated dose-response relation. To fill this knowledge gap in QMRA analysis, a new, sample-size-dependent bootstrap (95%) confidence band construction algorithm is proposed which can normally contain 95% of the data points and the results are validated by an extensive Monte Carlo simulation study using four well studied experimental dose-response datasets.

The rest of the paper is organized as follows. Section 2 begins with a description of the general setting for this study which includes: the notation and definitions; datasets; model specification and parameter estimation. The existing and the new bootstrap confidence band algorithms are specified in details in Section .The comparisons between the existing and the new algorithms were made and the results are presented and discussed in the last part of Section 2. The paper is completed with a conclusion section. For a better reading flow and repeatability of our research findings, an appendix section is added at the end which contains: table 4, the details of the selected real doseresponse datasets; table 5, the core parts of R code programs [11] for random sample generation from a beta-Poisson doseresponse model.

Data, model specification and parameter estimation

Dataset, Notation and Definition

Well-studied dose-response experiment datasets provide a reliable baseline for validation of our research findings. In this study, four experimental datasets from the literature are analysed for comparison between the single-hit beta-Poisson model

P_{I} (d | α, β)

and the generalized beta-Poisson model

P_{I} (d | α, β, r *)

with respect to parameter estimation and model performance. The selected experiment datasets are: rotavirus (CJN) and infection in healthy volunteers (Dataset 1) [8, 10]; Campylobacter and infection in healthy volunteers (Dataset 2) [8, 10]; Salmonella (Nontyphoid Strains) and infection in human volunteers (Dataset 3) (p399,[12]); Listeria monocytogenes and infection in mice (Dataset 4)[2]. These datasets have been used for dose-response analyses in [2, 8, 10, 12], although these references are not the original data sources. These four datasets are chosen because of their variety of coverage, e.g., virus and bacteria, few data points (e.g., ≤ 10) and many data points (> 40), human volunteers and mice, and old (back to 1996) and recent research sources. For the purpose of a clear and consistent model specification and easy comparison with the literature results, we have adopted a notation scheme as detailed in table I throughout this paper. The term ‘organism’ is used as a short name for pathogenic microorganism.

Model Specification

Beta-Poisson dose-response models are derived based on a plausible conceptual process with the following assumptions [14]: (1) only one viable organism is required to initiate the infection process in vivo; (2) the exact number of organisms being ingested in each dose representing a random sample from a Poisson distribution (with a known mean number d); (3) the survival of any organism within a single host is independent survival of any other organism within that host.

A detailed description and derivation of the exponential and the beta-Poisson dose-response models are given by Haas et al..[1] For the purpose of this paper, the following model specifications are needed.

If we assume that the actual number of organisms ingested, D, follows a Poisson distribution, i.e.,D~ Poisson

(D_{i}); r

a constant (a homogeneous susceptibility assumption); and K_min = 1 (a single-hit assumption), the simplest QMRA dose-response

P_{I} (d) = 1 - exp(- r d), (1)

where is the mean/effective dose.

In Equation (1), if we take into account the individual susceptibility/host sensitivity by allowing rtofollow a beta distribution with the probability density function (pdf)

f (r | α, β) = \frac{Γ (α + β)}{Γ (α) Γ (β)} r^{α - 1} {(1 - r)}^{β - 1}, (2)

where

0 < r < 1,

and parameters

α > 0

and

β > 0

we obtain the well-known beta-Poisson (single-hit) dose-response model[1]:

P_{I} (d) = 1 -_{1} F_{1} (α, α + β, - d), (3)

where

_{1} F_{1} (α, α + β, - d)

is the Kummer confluent hypergeometric function [15] defined by

_{1} F_{1} (α, α + β, - d) = 1 + \frac{Γ (α + β)}{Γ (α)} \sum_{j = 1}^{\infty} [\frac{Γ (α + j)}{Γ (α + β = j)} \frac{{(- 1)}^{j - 1} {(d)}^{j}}{j!}] . (4)

Since there are no analytic solutions to

_{1} F_{1} (α, α + β, - d),

only numeric approximation solutions or asymptotic solutions can be obtained[15, 16]. In their original paper, Furumoto and Mickey [17] derived the simple, attractive approximation beta-Poisson dose-response formula:

P_{I} (d) = 1 - {(1 + \frac{d}{β})}^{- α}

When

α ≪ β

and

β ≫ 1, P_{I} (d)

obtained from Equation (5) is a very good approximation to the true values of

P_{I} (d)

defined by Equation (3) [1, 10]

A single-hit model reaches its maximum risk limit when r=1 in Equation (1). This is an important property of the beta- Poisson model (Equation (3)) which provides an upper bound of the confidence level of the dose-response relation. It is worth noting that this maximum risk limit property is not retained in the approximate beta-Poisson model (Equation (5)) [10]

Model Parameter Estimation and Statistical analysis

Let denote the value of a ─2 log-likelihood ratio and we term this the deviance. According to [1](p285, Equation 8.33), for dose-response data analysed in this study, Y can be calculated as

Y = - 2 \sum_{i = 1}^{m} [y_{i} log \frac{π_{i}}{π_{i}^{o}} + (N_{i} - y_{i}) log \frac{1 - π_{i}}{1 - π_{i}^{o}}]

where

π_{i}

is the predicted response, typically the predicated probability of infection;m is the number of population groups at risk in terms of the mean dose levels as defined in Table 1. The true response (i.e., the true probability of infection) is estimated by

π_{i}^{o} = y_{i} / N_{i},

the ratio of the observed number of infected individuals to the number of all individuals at risk at each mean/ effective dose level

D_{i} .

Table 1: Notation and definition

Notation	Definition
Ds= ${D_{1,} D_{2}, ..., D_{m}}$	mean dose levels (i.e., the average number of organisms per dose or average concentrations); where the meaning is clear, denote $d \equiv D_{i}$ for $i = 1, 2, ..., m$ .
N = ${N_{1}, N_{2}, ..., N_{m}}$	$N_{i}$ is the number of individuals participated in the i th exposure group ( $m$ groups, each corresponding to one mean dose level $D_{i}$ ).
y = ${y_{1}, y_{2}, ..., y_{m}}$	$y_{i}$ is the number of infected individuals in the th exposure group
r	survival probability of each single organism ingested
K	number of organisms surviving to the target site inside a host body
K_min	the threshold value of the surviving organisms, i.e., the minimum number of organisms required for causing an infection event
$P_{I}$	probability of infection estimated by a dose-response model
$π_{i}^{o}$	observed proportion of infected individuals in the i th exposure group, i.e., $π_{i}^{o} = y_{i} / N_{i} for i = 1, 2, ..., m .$
$π_{i}$	predicted binary response by a model, e.g., if response is defined by probability of infection, $π_{i} = P_{I} (D_{i}) \equiv P_{I} (d)$ .
exp(.) and log	exp(1) = $e \approx 2.71828$ is the natural logarithm base such that log(e) = 1.
$Γ (.)$	the gamma function [13]

Equation (6) is a measure of deviance and can be used as an objective function for minimisation in model parameter estimation [18].The results are the maximum likelihood estimates (MLE) [1, 19]. Note that

π_{i} = P_{I} (d)

defined in Table 1.Then, the MLE of the model parameters, α ̂and β ̂, can be obtained using Equation (3) for an exact beta-Poisson model or Equation (5) for an approximate beta-Poisson model. The evaluation of the Kummer confluent hypergeometric function (Equation (4)) could be non-trivial. The statistical data analyses in this study were conducted using the open source statistical software R [11]. The R function ‘hyperg_1F1’ in the R package ‘gsl’ [20] is employed for evaluation of Equation (4). The built-in optimization function ‘optim’ in R is employed to obtain α ̂ and β ̂ by minimizing the model deviance Y defined by Equation (6).

A new Bootstrap confidence band Algorithm

A bootstrap approach is commonly accepted for construction of confidence bands on dose-response curves in microbial risk literature [1, 2, 8, 10]. Following the descriptions given in Haas et al. [1]and in Teunis et al. [8], the bootstrap algorithm employed in the literature (the existing algorithm) is summarized in the left column in Table 2. The new bootstrap algorithm, which produces sample-size-dependent confidence bands for prediction, is summarized in the right column of the same table. The comparison of both algorithms is depicted through figures1 to 7 for the four selected datasets.

Table 2: Parametric bootstrap algorithms for construction of 95% confidence bands about the mean dose-response curve

Existing bootstrap algorithm (confidence band for means)		New bootstrap algorithm (confidence band for prediction)
Step 1	Generate a large number of simulation samples, e.g., N=2000, based on Ds, N, $\hat{α}, \hat{β}$ ; Determine the plotting mean dose levels (D_sim*).	Step 1	Determine at which mean dose levels (Dsim) and the number of people (N_sim={max(N)}) exposed to hazard for generating simulation samples;
Step 2	Fit a beta-Poisson model to the generated simulation sample so that 2000 sets of ${{\hat{α}}^{}, {\hat{β}}^{}}$ are obtained;	Step 2	Generate a large number of simulation samples, e.g., N*=2000, based on D_sim, N_sim, $\hat{α}, \hat{β}$ ;
step 3	Calculate the probability of infection based on Equation (3) or (5) using parameter estimates obtained in Step 2 (2000 values at each D_sim level );	Step 3	Calculate the probability of infection using the ratio y_sim/N_sim, where y_simare the simulation sample values from Step 2 (2000 values at each D_sim level);
Step 4	Obtain a 95 percentile (2.5% to 97.5%) confidence interval for each D_sim level using the results from Step 3; connect all lower/upper bound points to form a confidence band about the infection probability curve.	Step 4	Obtain a 95 percentile (2.5% to 97.5%) confidence interval for each D_sim level using the results from Step 3; connect all lower/upper bound points to form a confidence band about the infection probability curve.

Note that the generation of random samples of ysim

= {y_{s i m} (1), y_{s i m} (2), ..., y_{s i m} (m)},

given Ds =

{D_{1}, D_{2}, ..., D_{m}}

(or D_sim) and N=

{N_{1}, N_{2}, ..., N_{m}},

is essential for employing a bootstrap algorithm to construct confidence bands. For a selected beta- Poisson dose-response model, the parametric bootstrap [9] random samples, ysim, is generated using R code programs provided in table 5.

The two different approaches to estimate the probabilities of infection,

P_{I} (d),

have a substantial impact on the resulting confidence bands. With respect to each bootstrap sample (i.e., the simulated data) ysim, the parameter estimates

{\hat{α}}^{*}

and

{\hat{β}}^{*}

are obtained first and the existing algorithm then estimates

P_{I} (d)

using Equation (3) (for an exact beta-Poisson model) or Equation (5) (for an approximate beta-Poisson model), whereas the new algorithm directly calculates the ratio ysim/N_sim as the model estimates for

P_{I} (d)

.

Although both algorithms construct the parametric bootstrap confidence bands, it is important to realize that the band based on the existing algorithm counts for the variation of the probability of infection given a mean dose level while the band based on the new algorithm counts for the variation for the probability of infection given a specific dose.

Therefore, the difference between these two confidence bands is an analogy of the difference between confidence bands for a regression line and for prediction of a new data point in a linear regression model. Furthermore, these two algorithms are fundamentally different in evaluating the discrepancy between Equation (3) and Equation (5): the existing algorithm treats the discrepancy as a difference caused by two different dose-response models (i.e., the exact versus the approximate models) while the new algorithm treats the discrepancy as a difference caused by two different parameter sets under the exact beta-Poisson doseresponse model.

Because the new algorithm employs the proportion data as the estimates for

P_{I} (d)

significant data discretization effects are expected in the resulting confidence bands when the maximum number of

N_{i} f s

(denoted by max(N)) is relatively small (e.g., max(N) < 15). Because the observed proportion of infection

\hat{π_{i}^{o}} = y_{i} / N_{i},

where

π_{i}^{o} = {log}_{N_{i} \to \infty} (\frac{y_{i}}{N_{i}}),

is used as the optimization target in the

P_{I} (d)

parameter estimation process for Equations (3) and (5), the (median) curve estimated using the new algorithm should converge to the analogous value using the existing algorithm for large max(N), e.g., max(N) > 50. Whereas the large sample assumption is implicitly assumed in the existing algorithm, the new algorithm is critically dependent on the choice of N_sim. Smaller (larger) N_sim would result in a wider (narrower) confidence band. With respect to each mean dose level

(d \equiv D_{i}),

the number of individuals exposed to a hazard,

N_{i}

(or N_sim) may be considered as a sample size in a binomial process. In this sense, the new algorithm confidence bands are, therefore, sample-size-dependent (i.e. N_sim =max(N)-dependent) These theoretical properties are examined in the Figures presented in next section.

Comparison Results and Discussions

Results and Discussions for dataset 1

Using the rotavirus dataset (Dataset 1) as an example, the data and parameters used in both bootstrap algorithms for confidence band construction are listed in Table 3.Figures1 and 2 are created based on Table II (bootstrap algorithm procedures) and table 3 (input information).

Table 3: Data and parameters used for bootstrap confidence band construction in Figures1 and 2

Ds	{ 0.009, 0.09, 0.9, 9, 90, 900, 9000, 90000}
N	{5, 7, 7, 11, 7, 8, 7, 3}
y	{0, 0, 1, 8, 6, 7, 5, 3}
Approximate B-P model:		$\hat{α}$ = 0.253; $\hat{β}$ = 0.422	Exact B-P model:	$\hat{α}$ = 0.167; $\hat{β}$ = 0.191
D_sim	{0.005, 0.009, 0.05, 0.09, 0.5, 0.9, 5, 9, 50, 90, 500, 900, 5000, 9000, 50000, 90000}
N_sim=max(N)	{11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11 }

In figure1, within the observed mean dose level range, the smooth confidence bands (dot-dashed lines in green colour, based on the existing algorithm) are consistently and noticeably narrower than the less smooth confidence bands (dotted lines in black colour, based on the new algorithm) while the median lines (solid lines in green or black colour) match well with each other.

Table 4: Selected dose- response experiment datasets from the literature

Dataset 1: Rotavirus (CJN) and infection in healthy volunteers [8, 10]
Ds (mean dose)	N (total)	y(infected)
9×10^-3	5	0
9×10^-2	7	0
9×10^-1	7	1
9	11	8
9×10¹	7	6
9×10^2>	8	7
9×10³	7	5
9×10⁴	3	3
Dataset2: Campylobacter and infection in healthy volunteers [8, 10]
Ds(mean dose)	N(total)	y(infected)
8×10²	10	5
8×10³	10	6
9×10⁴	13	11
8×10⁵	11	8
1×10⁶	19	15
1×10⁸	5	5
Dataset3: Salmonella (Nontyphoid Strains) and infection in human volunteers (p399,[12])
Ds(mean dose)	N(total)	y(infected)
1.52×10⁵	6	3
3.85×10⁵	8	6
1.35×10⁶	6	6
1.39×10⁵	6	3
7.05×10⁵	6	4
1.66×10⁶	6	4
1.5×10⁷	6	4
1.25×10⁵	6	5
6.95×10⁵	6	6
1.7×10⁶	6	5
1.2×10⁴	5	2
2.4×10⁴	6	3
6.6×10⁴	6	4
1.41×10⁵	6	3
2.56×10⁵	6	5
5.87×10⁵	6	4
8.6×10⁵	6	6
8.9×10⁴	6	5
4.48×10⁵	6	4
1.04×10⁶	6	6
3.9×10⁶	6	4
1×10⁷	6	6
2.39×10⁷	6	5
4.45×10⁷	6	6
6.73×10⁷	8	8
1.26×10⁶	6	6
4.68×10⁶	6	6
1.2×10⁴	6	3
2.4×10⁴	6	4
5.2×10⁴	6	3
9.6×10⁴	6	3
1.55×10⁵	6	5
3×10⁵	6	6
7.2×10⁵	5	4
1.15×10⁶	6	6
5.5×10⁶	6	5
2.4×10⁷	5	5
5×10⁷	6	6
1×10⁶	6	6
5.5×10⁶	6	6
1×10⁷	6	5
2×10⁷	6	6
4.1×107	6	6
1.5×106	6	5
7.68×106	6	6
1×107	6	5
1.58×105	6	1
Dataset4: (292 and 295 pooling) Listeria monocytogenes and infection in mice [2]
Ds(mean dose)	N(total)	y(infected)
2	6	0
5	6	1
1.1102	6	2
5.5×103	10	7
3.24×104	10	7
3.9×104	6	4
5.5×104	10	9
2.51×105	10	10
5.5×105	10	10
2.82×106	10	10

Ds= {D_1,D_2,…,D_m } representsm mean dose levels (i.e., the average number of organisms per dose or average concentrations) applied tom exposure groups; N = {N_1,N_2,…,N_m }where N_i represents the number of individuals participated in the ith exposure group (corresponding to mean dose level D_i); y = {y_1,y_2,… ,y_m }where y_irepresents the number of infected individuals in the ith exposure group.

Table 5:

This is the most important feature to note in the comparison of these two algorithms. With the new algorithm, the resulting confidence band is critically dependent on the choice of N_sim. By choosing N_sim= max(N) as proposed in table II and exemplified in table 3, the confidence band based on the new algorithm provides a lower bound of the band-width in estimation of the prediction confidence levels (i.e., the confidence band for prediction should be at least as wide). As expected, in all datasets we examined, the existing algorithm band was narrower than the analogous band based on the new algorithm within the observed mean dose level range. This provides strong empirical evidence to confirm the theoretical property concluded in the last section that the confidence bands based on the algorithm currently proposed in the literature are not confidence bands for prediction.

Figure 1: Comparison of bootstrap confidence bands (dataset 1). Plotting scheme: circles represent the data points; solid lines depict the estimatedmedian

P_{I} (d)

Figure 2: Comparison of bootstrap confidence bands (dataset 1, approximate beta-Poisson model). Plotting scheme: circles represent the data points; solid lines depict the estimated median

P_{I} (d)

, and the 2.5 and 97.5 percentile lines are drawn with dot-dashed lines (the existing algorithm in green colour) or with dotted lines (the new algorithm in black colour); the bold red dashed line is the asymptotic maximum risk curve (i.e., when N_sim→ ∞); the thin red dashed line is the maximum risk curve in finite sample case (i.e., when N_sim is finite), a bootstrap 97.5 percentile line.

Another important feature to note in Figure1 is that, except for the confidence band of the exact beta-Poisson model based on the existing algorithm (the green dot-dashed lines in figure 1(b)), all three other confidence band upper limit lines exceed the maximum risk line (bold red dashed line) in the verylow dose region (the bottom left corner).This exceedance issue is further illustrated and discussed in figure2.

An approximate beta-Poisson model can perform poorly in that its predicted P_I (d) curve may go above the maximum risk line (the bold red dashed line); figure1(a) essentially reproduced the result reported by Teunis in 2000[10]. Figure 2(a) is a reproduction of figure 1(a) but with one extra curve added in – the thin red dashed line in the farthest left part of the plot.

It is noted that, in the bottom left corner region of figure 2(a), the upper dotted line exceeds the upper dot-dashed line which is already above the bold red dashed line. However, the dotted line remains under the thin red dashed line. Therefore, the confidence band based on the new algorithm (dotted lines) does not fail here, because it is a sample-size dependent confidence band and the thin red dashed line is the finite sample maximum risk line (choosing the same N_sim values, as shown in table 3, to take account of the sampling variation in estimating the maximum risk level). On the other hand, the existing algorithm for computing the approximate beta-Poisson confidence band (dot-dashed lines) does violate the maximum risk limit for MLEs are obtained under the large sample assumption and the bold red-dashed line (N_sim→ ∞) applies as the maximum risk limit. Figure 2(b) is drawn to verify how the new algorithm confidence band will change as N_sim changes (N_sim = {200}, compared with N_sim = {11} in figure 2(a)). Figure2(b) shows that the thin red dashed line moves much closer to the bold red dashed line and the new algorithm band (dotted lines) shrinks significantly, with the upper dotted line being always under the thin-red dashed line. It is noted that the (existing algorithm) upper dot-dashed line even exceeds the thin red dashed line in figure 2(b) which confirms the finding reported by Teunis [10]: the maximum risk limit property is not retained in the approximate beta-Poisson model. Noting the fact that N_i is seldom more than 20 in the real dose-response data, it is, therefore, necessary to take into account the sampling variation effects in determination of the maximum risk limits using Equation (1).

As shown by figure 1(b), by ignoring the sampling variation effects, the existing algorithm could underestimate the maximum risk in the very low dose level range (e.g., mean/ effective dose d< 0.05or log(mean/effective dose) = -3 ) which is of particular interest (and concern) in microbial risk assessment practice.

Figure 3: Comparison of bootstrap confidence bands (dataset 2). Plotting scheme: circles represent the data points; solid lines depict the estimated median

P_{I} (d)

Results and Discussions for datasets 2, 3, and 4

Figures 3 through to 7 are produced for the analysis of datasets 2, 3, and 4 following the same procedures as exemplified with dataset 1. Similar to figure1, figure 3 shows that the median PI (d) lines match very closely for both algorithms and the black dotted lines (the new algorithm) form a wider confidence band than the band formed by green dot-dashed lines (the existing algorithm) over the observed mean dose level range. In contrast toFigure1, however, data points in figure 3 only spread in the high dose level range (log(mean dose) > 6.5). In figure 3, the new algorithm band is much narrower than the existing algorithm band outside the observed mean dose range (roughly where log(mean dose) < 6.5). It is also noted that the new algorithm bands are well under the limiting maximum risk (bold red dashed) line. In contrast, the approximate beta-Poisson model confidence band based on the existing algorithm (figure 3(a)) badly exceeds the maximum risk line in the lower mean dose level range (roughly where log(mean dose) < -2). With a 95% confidence band, we expect to observe, on average, only one out of 20 data points falling outside the band. However, we notice that, with six data points, the existing algorithm confidence band clearly missed one point and another point is only marginally included. In contrast, all six data points are well included in the new algorithm band.

Figures 4 and 5 are drawn in order to show the impacts of the sampling variation on the new algorithm confidence band while the existing algorithm band should remain unchanged. The sampling variation is determined by the sample size N_sim. The plots in figure 4(a), (b), (c), and (d), show how the bands, constructed based on the new algorithm, change as the sample size N_sim changes.

Figure 4: Comparison of bootstrap confidence bands by setting N_sim to different values (dataset 2, approximate beta-Poisson model). Plotting scheme: circles represent the data points; solid lines depict the estimated median

P_{I} (d)

, and the 2.5 and 97.5 percentile lines are drawn with dot-dashed lines (the existing algorithm in green colour) or with dotted lines (the new algorithm in black colour); the bold red dashed line is the asymptotic maximum risk curve (i.e., when N_sim→ ∞); the thin red dashed line is the maximum risk curve in the finite sample case (i.e., when N_sim is finite), a bootstrap 97.5 percentile line.

Figure 5: Comparison of bootstrap confidence bands (dataset 2, exact beta-Poisson model, N_sim={2} ). Plotting scheme: circles represent the data points; solid lines depict the estimated median

P_{I} (d)

Note that the thin red dashed lines move away / diverge from the bold red dashed lines as the sample size decreases. Figure5 is with the Nsim = {2} setting, i.e., it is assumed that only two subjects are exposed to hazard at each mean dose level. Settings in Figures4 and 5 are supposed to cover most of the real life experimental scenarios. Figure4 examines the approximate beta-Poisson model situation and Figure5 examines the exact beta-Poisson model situation. All these cases have shown a consistent pattern – the existing algorithm confidence band is possibly too wide outside the observed mean dose level range. Since the parameter estimation results are

\hat{α} = 0.145; \hat{β} = 7.59;

(approximate beta- Poisson model [10]) for dataset 2. This implies that the confidence band for dataset 2, no matter whether the parameter estimates are obtained from Equation (3) or Equation (5), should almost always be valid. While this numeric evidence is consistent with the new algorithm bands as shown in figures3, 4, and 5, the existing approximate beta-Poisson model bands clearly contradict this expectation. The result of the deviance difference analysis (not reported in this article) also suggested that, with dataset 2, the approximate beta-Poisson model is a very good approximation. All these assessments strongly suggest that the confidence band based on the existing algorithm may not be a proper estimation of the degree of uncertainty of the dose-response relation in these circumstances. Figures 6 and 7, which present the confidence band comparison for datasets 3 and 4, provide more evidence to show the comparatively poor performance of the confidence bands based on the existing algorithm.

Figure 6: Comparison of bootstrap confidence bands (dataset 3, exact beta-Poisson model). Plotting scheme: circles represent the data points; solid lines depict the estimated median

P_{I} (d)

Figure 7: Comparison of bootstrap confidence bands (dataset 4, exact beta-Poisson model). Plotting scheme: circles represent the data points; olid lines depict the estimated median

P_{I} (d)

Figure 6 (a) is essentially a reproduction of figure 9-2 in [12] (p402) and figure 6(b) shows the corresponding confidence band based on the new algorithm. It is obvious that only a minority (roughly 10 out 47) data points are included in the existing algorithm band while the new algorithm band contains 45 out of 47 data points as expected. Similar comparisons can be made for dataset 4 as shown in figure 7(a) and (b) -- 5 out of 10 data points are outside the existing algorithm band while all data points are contained in the new algorithm band.

The literature dose-response analysis results of dataset 4 can be found in [2].

http://qmrawiki.msu.edu/index.php?title=Listeria_monocytoge nes_%28Infection%29%3A_Dose_Response_Models

(accessed on 23/04/2017). In the analysis, the author(s) compared the fit of the beta-Poisson and exponential models to the data and concluded that the beta-Poisson model was a adequate model. Figure 7 depicts the results of the doseresponse analysis with dataset 4. Clearly, the new algorithm has provided a much better fit.

Conclusions

A new parametric bootstrap confidence band algorithm is proposed for the popular QMRA beta-Poisson dose-response model in this study. The resulting confidence band is sample-size dependent (i.e., N_sim-dependent)and it provides the lower bound of the band-width for estimating the true confidence level for prediction. Through an intensive Monte Carlo simulation study, the performances of the confidence bands based on the existing literature algorithm and the new algorithm are compared by fitting the beta-Poisson dose-response models to four well studied datasets. It has been shown that, asymptotically (i.e., as N_sim→ ∞) the estimated median dose-response curve based on the new algorithm converges to the existing algorithm median curve; in finite N_sim case, these two curves differ only in smoothness (the existing algorithm curve is always smoother). The new algorithm confidence band substantially improves the existing algorithm band in data point coverage. Due to the big sampling variation in proportion data under small sample conditions, the existing algorithm is likely to underestimate the maximum risk limit of dose-response data. On the other hand, outside the observed mean dose level range, the existing algorithm tends to overestimate the upper bound of the confidence band. This study has shown that the new algorithm is able to facilitate both types of bias. Therefore, the confidence band for prediction based on the new algorithm is a better representation of the degree of uncertainty of an estimated dose-response relation. The upper limit of the estimated prediction confidence band could be used as a better estimate for the worst case scenario in risk assessment.

Acknowledgements

The early version of this article was finished in 2015 when the author worked for the pond project which was funded by Department of Science, Information Technology and Innovation (DSITI) of Queensland government, Australia. The author would like to thank Professor Kerrie Mengersen (Queensland University of Technology) for her valuable feedback on the early version of the draft manuscript.

References