FINDING SIMULTANEOUS TOLERABLE DOSAGE COMBINATIONS FOR MULTIPLE OUTCOMES AND MULTIPLE STRESSORS

Introduction

Unavoidable hazardous chemicals are used in our daily lives. In agriculture, for example, pesticides and fertilizers are essential for efficient food production, but they may have adverse effects on consumers and the environment. In our homes, cleaning products contain a mixture of hazardous chemicals that could affect our health. It is therefore of interest to researchers to find the chemical dosage that maximizes productivity while minimizing the side effects. Understanding a chemical’s acceptable threshold helps in decreasing any adverse effects. Numerous toxicological studies have been conducted to model the adverse effects of chemicals that associated with dose-response data, such as Mumtaz et al. ^[1998], Groten et al. ^[2001], Ivanova et al.

^[2009], and Chen and Wang ^[2018]. Doseresponse models explain the relationship between chemical dosage and the response of the experi mental subject. Figure 1 (a) represents the dose-response model as a solid curve fitted to the data; it shows a decreasing relationship between chemical dosages on the x-axis (x = Dose) and the response on the y-axis (y = Response). In cases where there is more than one variable defines a response, each subscale is considered as an outcome. For example in a pesiticide study by Moser et al. ^[2005] conducted on rodents, neurotoxicity is the desired response to be measured. However, there is no clear definition of neurotoxicity as it can be measured by several outcome variables such as Motor Activity, Stimulus Respose (such as tail pinch in rats), Blood Cholesterinase and Brain Cholesterinase.

As there are multiple outcome variables in the response one needs to consider each subscale as each outcome variable may be more or less sensitive to each of the multiple chemicals (stressors). At tempts to compress the multiple outcomes into a single composite response can result in considerable loss of information. The tolerable dose is the chemical dose that can be administered without causing adverse effects. Popular chemical thresholding techniques use an ANOVA approach, such as the No Observable Adverse Effect Level (NOAEL) and the Lowest Observable Adverse Effect Level (LOAEL); in experiments, these thresholds determine the dose that causes no adverse effect and the dose that causes the lowest adverse effect, respectively, as shown in Figure 1 (b).

For more about chemical thresholds, see Cox ^[1987], EPA ^[2012], Jakubowski and Czerczak ^[2014], and Chen and Wang ^[2018]. Instead of using NOAEL, Crump ^[1984] developed a regression approach leading to the Benchmark Dose (BMD), which is the dosage that defines a predetermined response level of an adverse effect. NOAEL and LOAEL use only single points in the data and ignores the shape of the dose-response model, whereas BMD is calculated by using the curve fitted to the dose response data, so all information is considered. To use BMD, one must first specify the percentage of the response that is considered to be an acceptable adverse effect; which is known as the Benchmark Response (BMR), as shown in Figure 1 (c).

Note that the BMR depend on certain levels (percent) of control, denoted by % as BMR with particulate interest associated with 1%5%10%, for example Crump ^[1984] also introduced the Effective Dose (ED) to determine the effectiveness of a dose, combining BMD and BMR in one notation, ED . These thresholds define the tolerable region (TR) which is the region under the dose-response model that contains the tolerable dosage with acceptable adverse effects, which is determined by ED . The level of the effective dose depends on the value of , which is set according to research goals. For example, consider a 10% change in the control group to be the acceptable adverse effect, known as ED90.

Figure 1 represents ED90, which is bounded by the two dashed lines. The horizontal line corresponds to the 10% reduction in the adverse effect, while the vertical line corresponds to the associated dose, which is at 2 dosages. The rectangle area bounded by the ED90 is the TR containing the acceptable dose at a 90% adverse effect; for more, see Boone et al. ^[2015]. In our study, we use ED50, the median effective dose, which causes a 50% reduction in adverse effect on the exposed subject. Figure 1 Figure 1: The dose-response model is fitted to simulated data by the solid line. The region of the tolerable dose is bounded by the percent to control vertically at ED90 and the chemical dose horizontally at two dose, presented by dashed green lines in panel (a).

NOAEL, LOAEL and ED90 are introduced in panel (b), in red, blue and green dashed lines respectively. BMD and BMR are explained in panel (c), where the bounded rectangle area represents the tolerable region; this shows the safest dose at 10% reduction in the percent to control. There are different experiments that could be utilized for dose-response studies. In other words, toxicologists consider experimental design to determine the chemical tolerable dose and the TR, since there are undoubtedly more efficient designs of experiment. Use of the optimal design contributes to fulfilling the research goal by ensuring an accurate finding. From Chaloner and Verdinelli ^[1995] ”Experimental design involves the specification of all aspects of an experiment.

Common sense, available resources, and knowledge of the motivation for carrying out the experiment often help in selecting important features that depend on the specific problem.” Ray designs, as in Gennings et al. ^[2004], are a popular choices for toxicology experiments. They allow toxicologists to dose the experimental unit by increasing the dose incrementally in the form of a ray, where the second dose exceeds the first. Figure 2 presents a ray design in which x1 and x2 are two different chemicals. Figure 2 (a) shows a one-ray design, Figure 2 (b) shows a two-ray design. For more about ray designs, see Stork et al. ^[2007]. Figure 2 Figure 2: Ray design plots involving two different chemicals x1 and x2.

Panel (a) represents a one-ray design, panel (b) a two-ray design. Design of experiments (DOE) can be regarded as a sequential process where outcomes from an initial experiment can be used to guide them to additional experiments. Follow up experiments are often necessary to improve precision or resolve ambiguity, and several techniques exist for selecting follow-up runs. One such approach relies on alphabetic optimality criteria, which had received much attention in engineering and science applications, for example, in studies by Clyde ^[1994], M¨ uller ^[2005], Chu and Hahn ^[2008], and Guest and Curtis ^[2009]. Common criteria for optimal experimental designs are written as functions of the information matrix. For example, D-optimality minimizes the generalized variance of the parameter estimates, whereas A-optimality minimizes the average variance of the parameter estimates.

Other optimal experimental design types are mentioned in, Silvey ^[1980], Chaloner and Verdinelli ^[1995], and Atkinson et al. ^[2007]. In this article, we introduce a novel Bayesian criterion based on the variance of the area of the tolerable region to determine the optimal follow-up design points, which may then be used to better define the safest chemicals such as pesticide exposure limits. This work extends the previous work of Boone et al. ^[2015] and Farhat et al. ^[2019] who focused on determining a tolerable dose region in experiments with multiple outcomes and multiple stressors, with the former looking at additive stressor effects and the later considering where the stressors interact.

However, there is no guidance on experimental design for researchers who wish to design an experiment that considers multiple outcomes and multiple stressors. In Section 2 a review their modeling methodology is given followed by a description of the approach of Farhat et al. ^[2019] for determining the Bench Mark DoseTolerable Region in Section 3. Section 4 presents a new optimality criterion for follow up design points for these types of experiments. The remainder of the manuscript considers several illustrative examples in Section 5 and discussion and conclusion in Section 6.

Methodology

Bayesian Model To determine the relationship between chemicals and their side effects in a dose-response model, a Bayesian approach is employed model that accounts for the interaction among the chemicals. The experiment here will be conducted in a manner that experimental unit i is administered a single dose of a combination of several chemicals (stressors), then each of the outcome variables are measured. Let J be the number of outcome variables, K be the number of chemicals (stressors) that the experimental unit will be dosed with and n be the number of experimental units. The response for the ith experimental unit on the jth outcome variable is denoted yij, which is then put in vector form as a J 1 vector, Yi = (yi1 , yi2 , .

. . , yiJ)T , for i = 1, 2 , . . . , n. Let xik be the dosage of the kth chemical on the ith experimental unit, which is then formed into a vector Xi = (xi1 , xi2 , . . . , xiK ), for i =12 nwhichrepresentsthesingledosechemicalcombinationthatisadministered to the experimental unit.To relate the chemical combination dose to each out come variables there are J functions f1 f2 fJ where E[yijXXXi]=fj XXXi j and j is the vector of parameters associated with function fj. Note that for this work, all the functions fj are assumed to be monotonic and hence invertable ,i.e. f 1 exists.

In addition for this work , the j will be vectors of regression coefficients that can be organized as Research Article β = (β1 , β1 , . . . , βj) Putting all the items together the expected value of the response vector Yi for the ith subject after being exposed chemical dosage combination Xi is as follows: Here ri is a random effect to allow for the with in subject effect across the end points. For clarity, it should be noted that experimental unit i is dosed with a single dosage combination,Xi,and each of the out come variables are measured from this one dosage combination.Hence there is no index j on Xi.

The model in equation (1) allows for each outcome variable to be modeled separately for each dose combination while considering the within subject effect. Note that this formulation will allow flexibility when specifying the likelihood associated with each experi- mental outcome as it may differ from outcome to outcome. For example, one outcome may be a continuous random variable where a Gaussian likelihood is appropriate and another outcome may be a discrete random variable such as Bernoulli or Binomial, which happens to be the case in the Moser et al. ^[2005] dataset. The requirement that each of the mean functions f j be monotonic, and thus invertible, will allow for easy calculation of a tolerable dose combination given a specified level of risk.

As will be demonstrated in the Section 3 the tolerable dose combination will be a set (the preimage) not a single value and the linear component X i T β j will aid in its determination. Prior and Posterior Distribution As a Bayesian approach is being employed the likelihood must first be defined. In this case the likelihood must be defined for each of the outcome variables. Let g j be the likelihood function for outcome variable j. Specifically, the likelihood of yi j is given by g j yi j| f j, X i, β j, ri, γ j , for j = 1, 2, . . . , J.

Here γ j are parameters that may be needed to completely specify the probability distribution. For example, if yi j has a Gaussian distribu- tion then the variance parameter, σ2 is needed and hence γ j = σ2. Or in the case where yi j has a Binomial distribution then the number of trials m is needed and hence γ j = m. We are considering the Bayesian approach, therefore, we have to define the likelihood for the data, the prior distribution for parameter β and the process of sampling from the pos- terior distribution. The probability density function that represents the response measure- ments of the ith subjects on the jth outcome is given by g j(yi j|XT , β j, r j, γ j), j = 1, 2, .

. . , J. The γ j are extra parameters in the distribution. Knowing that the β′js and γ′js outcomes where the y j are conditionally independent. The likelihood is presented as follows: For consistency with previous work the prior distributions as specified the same as those by Boone et al. ^[2015] and Farhat et al. ^[2019]: β|µ, Ω ~ N(µ, Ω) µ|a, A ~ N(a, A) Ω|R, ρ ~ Wishart(R, ρ) γ j ~ p(γ j) ri|ω2 ~ N(0, ω2) ω2 ~ Gamma(1, 1) where β represents the regression coefficient vectors as defined previously, µ and Ω denote the mean vector and the accuracy matrix, respectively for β, R represents the correlation matrix with ρ being the degrees of freedom and a and A are hyper parameters for µ.

The as mentioned above γ j are additional parameters that may be needed for different outcome types the, prior distribution of γ j should be chosen accordingly. The prior distribution along with the likelihood specified produces the posterior distribution given in equation (2). Notice that the posterior distribution can not be expressed in closed form, and hence will be approximated using Markov Chain Monte Carlo (MCMC) techniques. , , , | | | exp exp exp exp exp exp exp P D r Y j j j j i i n T j j n j j i n ij x 1 1 3 2 i j # # # # # ?

b n x v n s s s b n b n v v X X X - - - - - - - - - b = = - = - - = - - = = ^ d ^ ^ ^ ^ h n h h h h ( ! ( ( & ( + 2 2 % / / / / / Different software to implement MCMC sampling may be used such as, WinBUGS, Open- BUGS, JAGS, STAN etc. For this work JAGS is employed. The analysis algorithm of the samples determines the favorable goodness of samples which convergence Gelman et al. ^[2013]. These Samples are drawn from the beginning of the process which are not accurate due to lack of stability of the Markov chain.

Trace plot analysis discard weak samples by burn-in, which is assesses the convergence of the chains where it is available for all parameters. In the literature there are different sampling tools from the posterior distribution, this sampling tool presented in Andrieu et al. ^[2004]. We are using the Bayesian experimental design for constructing more robust designs, for more information see Clyde ^[2001]. Research Article [ | ] | | ... | ... , ,..., ( ) E Yi Xi f X E yi Xi E yi Xi E yiJ Xi f X ri f X ri f X ri i n 2 1 i t i T i T j i T J 1 1 2 2 b b b b = = = + + + = J L K K K K K K K K K K K J L K K K K K K K K K K K ^ ^ ^ ^ N P O O O O O O O O O O O N P O O O O O O O O O O O h h h h 6 @ @ @ THE BMDTR EVALUATED BMDTR is the region in the chemical space, that represent the tolerable dose combinations at a specific level of interest with its corresponding tolerable exposure adverse effect.

In our research, we focus on measurements of three continuous outcomes (Exponential) as a result of exposure to the K chemicals. The Bayesian estimate of BMDTR is obtained via MCMC sampling (m = 1, 2, . . . , M) from the posterior distribution with both burn-in and thinning. In our proposed method we determine the BMDTR at chemicals dosages and the interactions of its adverse effects. For a continuous jth outcome, use the following dose response model for K = 2 chemicals, where the link function fi can be found as: Where the β0 j is the unknown intercept, βi j are the unknown slope parameters associated with each component, ri is the random effects, and (xki) are the doses of the chemicals.

The jth outcomes come from an exponential experiment. The model considered is as follows: Evaluation of the Tolerable region (TR) The BMDTR determines an acceptable level of multiple chemicals. Determining the BMDTR depends on specific measurements and thresholds that is analogous to BMD with specific level of BMRη as introduced in Boone et al. ^[2015]. To calculate the tolerable region, we fit the models of the j outcomes to the data. Then, the intersection area of all J tolerable re-gions TR = TJ j==1 TRj will define the BMDTR. BMDTR is defined across all generated MMCMC samples, therefore, for each MCMC sample of β(m) a BMDTR ,TR(m) is defined. The area A of TR is associated with the mth MCMC samples which denoted by ATR(m) and has the following formula ATR (m) = ∫ TR (m) dX, m = 1, 2, 3, .

. . , M. The tolerable area {ATR(m)}M/m=1 can be ranked and the lower bound on the area of the TR will be defined among the MCMC samples of β which corresponds to the adverse effect percentage of interest. For more explanation in Figure 3, x1 and x2 represent the two differ- ent chemicals that we are testing. We fitted the outcome models to the data as represented by the green and blue color curves, we find the intersection area of all outcomes. The small- est area with the smallest distance d from the origin across all MCMC samples will be the TR. The outcome with the smallest distance from the origin that is shown in Figure 3 by the black line, that will determine which outcomes play a role in defining the TR.

Then the area of the TR will be calculated as the red lines going across angles theta θ = θ1, θ2, . . . θi using the distance d formula d(m) (θ) = min j{d(m) j(θ) }. The area under the curve is A which defined by Where ai is the area of a triangle of each red line under the curve which has a single θ and form a triangle shape with the following equation. These equations are calculated across all samples to find the TR which has the minimum distance from the origin. In Figure 3 the curves parameters are not define here since its just for illustration of the TR calculation techniques.

Figure 3 Figure 3 : Fit of two models represented in blue and green curves. Defining the tolerable region under the fitted model by the red lines, which partition the area to calculate the ATR(m) . BMDTR Evaluation for all MCMC Samples The methods used to calculate BMDTR consider all MCMC samples (after burn-in and thinning). It follows the same step used in Boone et al. ^[2015] and Farhat ^[2014], specifi- cally in a way of finding the desired η value, determining the BMRη, and finding the link function for the appropriate model. Steps of evaluating the BMDTR using MCMC samples is as follows, S.1 Partition the chemical space by dividing the chemical space into N equally spaced angles (θ1, θ2, .

. . , θN) and increment the angle by u (for instance, π/180 or π/360 Research Article ), , ,. . . , Y N e i n 1 2 2 ij ri 2 j ijx i jx i jx ix i 0 1 2 2 3 1 2 + v = b b b b + + + + ^ ] ^ h g h exp f X ri 3 j T j j jx i jx i jx ix i 0 1 1 2 2 3 1 2 b b b b b = + + + + ^ ^ ] h h g , , . .

. , . Y i Exp e i n 3 1 4 r j jx i jx i jx ix i i 0 1 1 2 2 3 1 2 + = b b b b + + + + ] ] ^ g g h . ( ) A ai d over all area 5 i i i n = = / * * ( ) sin cos a d 2 6 i i i i = ] ] g g increments could be used), where , and on the other hand θv+1 = θv + u for all v and u. S.2 For each m=1 to M, successively with each angel: a pull a sample βm from the posterior distribution of parameters β.

b Consider the sample in a) and determine which outcomes have the minimum distance to the origin in the same time with . Now we let c Repeat step 2a and 2b for all MCMC samples. The result will be the minimum distance known as whichwillhave adistribution along the ray as shown in the figure. d Then we sort and pick the distance that is attached to the 5th quantile, which is and recalculate the value of theta associated with it. Though we considered the 5th quantile, it is possible to use a different quantile in the third step. BMDTR= with all . In this research we consider three different curves fitted with different parameters values β = (β1, β2, β3) that correspond to the curves in Figure 4.

The tolerable region is defined by the intersection of the three curves as shown in Figure 4 (a) that represent the fitted curves and Figure 4 (b) shows the curves with the shaded area as the TR.We see the black curve has with the smallest distance from the origin which define the TR, using the above Algorithm in Section 3.2. Moreover, Figure 4 show an explanation of the fitted models curve that define the TR and identify the TR in the shaded area. Figure 4 Figure 4: True Tolerable Region represented in panel (a) the fitted model where panel (b) shows the shaded area as the tolerable region. BAYESIAN OPTIMAL EXPERIMENT DESIGN In the literature, multiple optimality criteria have been developed for different goals of ex- perimental design, such as reducing the costs of an experiment, minimizing the variance of an estimator or to obtain the most useful experiment phase to perform statistical analysis.

For more see Verdinelli ^[1992], Clyde and Chaloner ^[1996] and Atkinson ^[1996]. Optimal design details and mathematical theory are introduced in Pa´zman ^[1986] and Fedorov ^[1972], with more about optimal experiment design in Pukelsheim ^[2006]. Bayesian optimal design is sensitive to the prior and posterior distribution of parameters model and in our research we are considering the posterior predictive distribution, since the posterior distribution does not count for future runs in an experiment. Posterior predictive distribution is known as the distribution of possible unobserved data conditional on the observed data, it has the following formula. , where x˜ is the new observed data, X is the existing observed data and β is the parameter β є Θ.The new observed data x˜ obtained by MCMC sampling method and that is independent of the sample data X .

For more about Bayesian optimal experimental design see Tsutakawa ^[1972], Tsutakawa ^[1980], Chaloner et al. ^[1984], Chaloner and Verdinelli ^[1995], Chevret ^[2006] and Wang et al. ^[2013]. Bayesian Follow-Up Design A follow-up experimental design is often needed to improve precision of model parameter estimation or to further explore the experimental region. Optimality criteria (such as D- optimality) are often used to select follow up runs. Here we consider finding the optimal additional control sets (doses of chemicals) to better define the TR. Recall the TR is defined as the intersection of the Benchmark dose limit BMDL as TR = BMDLi j, and the area of the tolerable region is define as ATRi = TRi dx.

In this article we propose a follow up design experiment that seeks to minimize the variance of the area of the TR. In order to define our criterion, we make use of the posterior predictive distribution, that is where Y2 is a vector of future observation, Y1 is the initial design points and e is the follow up design. Our criteria is therefore defined by Our goal here is to use this criterion in determining the Bayesian Optimal Follow-up Ex- perimental Design (BOFED) points, by finding the points that make the minimum variance of the area of the TR. To compute this criteria Equation (7); we generated data from an initial design points y1, which are the design points that have been defined from initial experiment.

After that, we fit our outcome models Equation (3) to the data. Then, we want to define the posterior predictive distribution which is a distribution that tends to be dif- ficult to obtain. Therefore, we obtain it by sampling from the posterior distribution then conditioning on these samples we then obtain new samples, repeating that a 100 times will then obtain the posterior predictive distribution. Then we predict at the future design points e which have 100 pairs of candidate points. One pair per design were examined on two rays. Then v v l i r = = / d* m i ] ] g g " , d* m i ] ] g g " , , .

. . , x x x , * * * p m 1 2 i ^ h , , , . . . , . min d d j J 1 2 ( ) *( ) ( ) *( ) m j m = = i i " " , , D * M D * M D * M d * * m i _ i , , . . . , x x x * * * p m 1 2 i ^ ] ] h g g { , , . . . , } x x x * * * p m 1 2 i ^ ] ] h g g , p x X p x X p X d y y y b b b = K ] ^ ^ g h h J R # , , , f Y Y e p Y Y e p Y d 2 1 2 1 1 ; y y b b b = b ^ ^ ^ h h h # , , , min Var A Y Y e f Y Y e dY TR Y 1 2 2 1 2 ; ; ^ h 6 @ # D e m o 7 we calculated the TR and repeated that 100 times across all MCMC samples to calculate the ATR while considering the new design point.

We calculated the variance across all 100 candidate design points and determined the minimum variance of the TR that corresponds to the optimal Bayesian follow-up experimental design e. The following algorithm explains the steps we used in our approach. Algorithm 1: Determining the optimal follow-up design points for m=1 to M do S.1 Samples are drawn to obtain posterior predictive distribution from the posterior distribution f (Y2|Y1, e). S.2 Predict future design points e for 100 candidate points along two ray design. S.3 Calculate the variance of the tolerable region across all candidate points. S.4 Determine the Bayesian optimal follow-up design points that have the minimum variance of the tolerable region using Equation (7).

end Figure 5 explains the above mathematical description. When samples are drawn from the data, a posterior predictive is generated about the new data points. These samples have variations when determining the TR, shown by blue and red curves. In Figure 5 (a) we see a variation around the fitted models which makes the TR unstable. In Figure 5 (b) the TR is more stable since we consider the minimum variance of the TR, which considers our optimal criteria. Figure 5 Figure 5: Panel (a) shows the possible variation around the curve intersect that defines the tolerable region across all MCMC samples where panel (b) represents the less variation that defines the tolerable region when using the criteria in Equation (7).

ILLUSTRATIVE EXAMPLE We need to find the minimum variance of the tolerable region with minimum dose. The ED50 process spent 42.8 hours per design of MCMC simulation. This draws a samples from the posterior distribution using the Gibbs Sampler JAGS Lunn et al. ^[2009]. The posterior distributions are used to determine the optimal fit of the study, sampling 100 times from the 5000 outputs as a pair of x1 and x2 from 1 through 10. We considered fitting the exponential dose-response models to two different cases. First Case: Design points that are generated from the initial design considering one ray design. Then, the optimal followup design points determined in a two-ray design.

Second Case: Considering design points that are generated from the initial design of two ray design. Then, the optimal follow-up design points determined in the two-ray design. One-ray designs are used in Meadows et al. ^[2002] considering the combination/mixture data of chemical. Base design and candidate design points are used to generate the study in this research as shown in Figure 2. Each ray has ten possible points and a pair of optimal points will be chosen among them. The optimal follow-up point is determined by the de- sign points that correspond to the minimum variance of the tolerable region under the fitted models. To start with an optimal base design, we replicate across the experimental points to determine the most stable area.

Among 10 replicates we found that 4 replicates has the largest drop in error, so we choose 4 replicates for our two initial designs. Replicates help ensure patterns that indicate the need for sequential design. Result Four different possible scenarios are considered in this research, as shown in Figure 6. Each scenario has different parameter values that will be introduced later, which change the tolerable region shape and the area in each scenario. Simulated data is used in this article with four replicates to study the four different scenarios. Each scenario has 10 possible dosages along two rays design. One pair of these dosages will be the BOFED point, which is determined by the criteria introduced in Section 4.1.

The different scenarios in Figure 6 are generated from the same simulated data, with a unique set of parameters for each scenario that used in the general Equation (3) which specify as follows for our study: For all scenarios the prior distributions specification are as follows: β ~ N(µ, Ω), µ ~ N(0, 100 I9), Ω ~ Wishart(I9,9), σ j ~ Gamma(1, 1), ri ~ N(0, ω2), ω2 ~ Gamma(1, 1), where I9 is a 9 × 9 identity matrix and other parameters are as defined in Equation (2). Likelihood and posterior distribution introduced in Section 2.1. Each scenario has three ex- ponential models with different parameters that contribute to the different scenarios, which change the fit of the models and changes the tolerable region for each scenario.

Panels in Figure 6 shows the three models that bound the TR, in panel (a) the TR bounded by the red line model which is the small area of that setting and the other curves affect the area but does not bound it. In panel (b) the TR is bounded by two models in red and green curves where the blue does not include in defining the TR, panel (c) shows the TR by the intersec- tion of two curves red and green and again, the blue curve affects but not including in the determination of TR, and panel (d) represent the TR as wide as possible and bounded by al- most the three curves.

Models of each scenario with the parameters settings are introduced here in Subsection 5.1.1 and 5.1.2. Figure 6 Figure 6: Scenarios, Panel (a) shows the First Scenario where the TR is the area bounded by the red curve, Panel (b) shows the Second Scenario where the TR is the area bounded by the red and green curve, panel (c) shows the Third Scenario where the TR is the area bounded by the green curve intersects with red curve and panel (d) shows the Fourth Scenario where the TR is the area bounded by the green curve. exp f X ri j T i i ix ix ix x 0 1 1 2 2 3 1 2 b b b b b = + + + + ^ ^ h h D e m o 8 Will these scenarios affect the Bayesian optimal follow-up design points?

To answer the question we investigated the two different cases for each of the four scenarios. When introducing the result for each scenario, we present a Figure with three different panels (a), (b) and (c). Panel (a) represents the scenario plot with the two ray design and two solid red points on the rays as the BOFED points and open circles that represent the candidate points for each scenario. Furthermore, visualization of the minimum variance distribution of the TR for all possible design points will be represented in a heat map in panel (b) and a heat map of the standard error distribution around the minimum variance of the TR for all possible design points is shown in panel (c).

Both heat maps have an x-axis and a y- axis ranked from 1 to 10 as it is all the possible design points. The heat maps have color gradients from dark red to light yellow, the dark color associated with lower minimum variance and light color are associated with large minimum variance. Each colored square represents points in the heat map that are read as x1 and x2 points which are the candidate design points. Our optimal point will be the point that has the minimum variance which will have the darkest color in panel (b). First Case First Case as introduced in Section 5 and shown in Figure 2 (a) which generates the follow- ing four scenarios: First Scenario, which has the following models: f1 (XT β1 ) = exp(−0.08664x1 − 0.086648x2 + 0.008x1x2 + ri) f2 (XT β2 ) = exp(−5 − 0.099x1 − 0.099x2 − 0.9x1x2 + ri) f3 (XT β3 ) = exp(−5 − 0.09x1 − 0.09x2 − 0.091x1x2 + ri) These models are defined using Equation (8).

Figure 7 (a) represents the first scenario with the BOFED points, the tolerable region bounded by the red curve and the BOFED points are indicated by the solid red points at (x1 = 9, x2 = 9). Point (x1 = 9, x2 = 9) is the point that has the minimum variance of the TR among all other possible candidate design points. We see in Figure 7 (b) the possible optimal design points will be the points with the dark red color which approximately ranges from point (4) to (10) on the x-axis and the y-axis ranging from (3) to (10). These points are determined by lower variance among other design points.

On the other hand, we see points with light yellow color gradients approximately ranging from points (1) to (3) in the x-axis with points from (1) to (10) in the y-axis, indicate high variance. Figure 7 (c) shows the heat map of the standard error around the minimum variance of the TR for each design point to check the precision of our optimal point. We see that all possible optimal points with the dark color are approximately points ranging from (4) to (10) on the x-axis to the y-axis ranging from points (1) to (10), which are around our optimal design points. That indicates (x1 = 9, x2 = 9) as an optimal follow-up design point for the first scenario since it has the minimum variance on Figure 7 (c).

Comparing the heat maps in Figure 7 (b) with (c), we see darker points in panel (c) that vary more among the heat map than in panel (b) with the same optimal design point in both heat maps. Figure 7 Figure 7: Panel (a) shows the 1st Scenario with Ray plot and the BOFED at point (x1 = 9, x2 = 9). Panel (b) represents the first scenario heat map of the possible design points and panel (c) shows the Standard error around all possible design points. Research Article (9) a. First Scenario b. Heat map of the Design points Heat map Second Scenario, has the following models: f1 (XT β1 ) = exp(−5 − 0.07x1 − 0.07x2 − 0.8x1x2 + ri) f2 (XT β2 ) = exp(−5 − 0.099x1 − 0.099x2 − 0.9x1x2 + ri) f3 (XT β3 ) = exp(−5 − 0.09x1 − 0.09x2 − 0.091x1x2 + ri) Equation (8) used to defined Equation (10) which constructs the curves in Figure 8 which defines the second scenario.

Figure 8 (a) represents the second scenario with the BOFED points (x1 = 7, x2 = 6) which has minimum variance and it indicated by the solid red points on the rays. The tolerable region bounded by the red and green curves. Figure 8 (b) represents the second scenario heat map of all possible design points with the minimum variance of TR, we see that approximately points ranging from (3) to (10) in each axis are possible optimal design points. All possible optimal design points are represented by the darker red color since it is associated with lower variance. On the other hand, we see Figure 8 (c) represents the heat map of standard error distribution around the minimum variance of the TR, which shows all possible variation of design points.

The standard error heat map in Figure 8 (c) is scattered as shown and the standard error for the BOFED points that has the minimum variance with a minimum standard error as well. Comparing the two heat map we see that Figure 8 (c) has more variation around our optimal point that indicates the precision of our point. Figure 8 Figure 8: Panel (a) represents the 2nd Scenario with Ray plot, BOFED point (x1 = 7, x2 = 6).Panel (b) shows the Second scenario heat map of possible design points and panel (c) shows standard error heat map around all possible design points. Third Scenario, has the following models: f1 (XT β1 ) = exp(−0.08664x1 − 0.08664x2 − 0.008x1x2 + ri) f2 (XT β2 ) = exp(−0.09664x1 − 0.09664x2 + 0.008x1x2 + ri) f3 (XT β3 ) = exp(−0.08664x1 − 0.08664x2 + 0.008x1x2 + ri) The third scenario has models derived from Equation (8), which constructs the plots in Figure 9.

The tolerable region is bounded by the intersection of the red and green curves and the BOFED points are indicated by the solid red points (x1 = 10, x2 = 9) as shown in Figure 9 (a). Figure 9 (b) is the third scenario heat map of all the possible optimal design points that range from (6) to (10) on both axis with lower variance. Point (x1 = 10, x2 = 9) is the BOFED point since it has the minimum variance among all other points. The points in Figure 9 (c) are the standard error distribution around the minimum variance of the TR for all possible optimal design points which has a scatter representation.

We see that our optimal point is not the same optimal point in Figure 9 (c) since it does not have the minimum variance. The standard error heat map indicates that there are other possible optimal design points. So further investigation and studies could be conducted in future research. Research Article (10) a. Second Scenario b. Design points heat map heat map (11) Figure 9 Figure 9: Panel (a) represents 3rd Scenario with ray plot, red solid points is the BOFED point (x1 = 10, x2 = 9). Panel (b) shows the third scenario heat map of possible design points and panel (c) shows the standard error around all possible design points.

Fourth Scenario, has the following models: The fourth scenario has the models that are defined in Equation (12) using Equation (8), which construct the curves in Figure 10 with the BOFED points. Figure 10 (a) shows the tolerable region for the fourth scenario bounded by the red and green curves. The BOFED points are indicated by the solid red points at point (x1 = 10, x2 = 10), which has the minimum variance. Figure 10 (b) shows the heat map of possible design points which range from (6) to (10) on the x-axis and from (1) to (10) 0n the y-axis. Figure 10 (c) represents the heat map of standard error distribution around the minimum variance of the TR around all possible design points as more scattered.

From the standard error heat map we see that there is an equivalent point to our BOFED point that could be investigated in further research. Figure 10 Figure 10 : Panel (a) represents the 4th Scenario with ray plot and the BOFED point in red solid points at (x1 = 10, x2 = 10). Panel (b) shows the fourth scenario heat map of the possible design points and panel(c) Standard error around all possible design points. In each scenario we see a variation in the BOFED points due to the considered sce- nario. The first and second scenarios have stable BOFED points, where the third and fourth scenarios show more variation around their BOFED points.

Comparing the first scenario to the third scenario, the first TR is close to the origin while the second TR is away from the origin and defined by a. Third Scenario a. Fourth Scenario b. Design points heat map b. Design points heat map heat map heat map . . . . . . ( ) . . . exp exp exp f X x x x x f X x x x x f X x x x x 0 09964 0 09964 0 008 0 09964 0 09964 0 008 12 0 08664 0 08664 0 008 T T T 1 1 1 2 1 2 2 2 1 2 1 2 3 3 1 2 1 2 b b b = - - + = - - + = - - + ^ ^ ^ ] ] ] h h h g g g the intersection of two curves causing variation we see in the heat maps and the determined BOFED.

Second Case Second Case as introduced in Section 5 and shown in Figure 2 (b) which generates the next four scenarios. First Scenario, has the same models as Equation (9). Figure 11 (a) represents the first scenario with the tolerable region bounded by the red curve and the BOFED points are indicated by the solid red points at point (x1 = 7, x2 = 9). Figure 11 (b) shows the heat map of all possible design points scattered where the darker color corresponds to minimum variance and the lighter color represents a high variance. Figure 11 (c) shows the heat map of the standard error distribution around the minimum variance of the TR around all possible design points, and we see point (x1 = 7, x2 = 9) has a dark color which has minimum standard error at this point indicating the precision of our BOFED points.

Figure 11 Figure 11: Panel (a) represents the 1st Scenario for second case with ray plot and the BOFED point at (x1 = 7, x2 = 9). Panel (b) shows the first scenario heat map of the possible design points and Panel (c) shows the Standard error around all possible design points. Second Scenario, has the same models as Equation (10). Figure 12 (a) represents the second scenario in the second case with the BOFED points in the red solid points at (x1 = 1, x2 = 10). The tolerable region is bounded by the red and green curves. Figure 12 (b) shows the second scenario heat map of all possible design points and its scattered spread and its read as the previous scenarios.

Figure 12 (c) repre- sent the heat map of standard error distribution around the minimum variance of the TR around all possible design points and we see that point (x1 = 10, x2 = 1) has the minimum standard error, and that is different from the BOFED point in panel (b) that indicate more investigation could be done to know other points. Figure 12 Figure 12: Panel (a) represent the 2nd Scenario with ray plot and the BOFED point at (x1 = 1, x2 = 10). Panel (b) shows the second scenario heat map of the possible design points and panel (c) shows the heat map of the Standard error around all possible design points.

a. First Scenario a. Second Scenario b. Design points heat map b. Design points heat map Third Scenario,has the same models as Equation (11). These models construct the curves in Figure 13 (a) represents the third scenario with the BOFED points and the solid red points at point (x1 = 9, x2 = 9) which has minimum vari- ance. The tolerable region in third scenario is bounded by the intersection of the red and green curves. Figure 13 (b) shows the heat map of all possible design points in the third scenario. Figure 13 (c) shows the heat map of the standard error distribution around the minimum variance of the TR around all possible design points.

Both heat maps are scat- tered and had more of yellow color at the design points which represent a higher variance at these points. The point with the minimum variance in Figure 13 (c) is different than the BOFED points in Figure 13 (b) which means there are equivalent points to our BOFED point and that could be investigated in future research. Figure 13 Figure 13: Panel (a) represent the 3rd Scenario with Ray plot and the BOFED point at (x1 = 9, x2 = 9). Panel (b) shows the third Scenario Heat map of the possible design points and Panel (c) shows the heat amp of the Standard error around all possible design points We see a huge difference in the heat maps comparing the first scenario heat maps to the third scenario heat maps.

That is due to the variation and the difference in the TR, we see the first scenario TR is close to the origin where the TR in the third scenario is a mixture of two curves and away from the origin and that caused the difference we see. Fourth Scenario, In the fourth scenario the models are defined by Equation (12) which construct the curves in Figure 14. Figure 14 (a) represents the fourth scenario with the BOFED points are indicated by the solid red points at (x1 = 10, x2 = 1) which correspond to the minimum variance. The tolerable region in the fourth scenario is bounded by the red and green curves.

Figure 14 (b) shows the heat map of all possible design points, it read as we mention in previous scenarios. Figure 14 (c) shows the heat map of the standard error distribution around the minimum variance of the TR around all possible design points. We see both heat maps are scattered and has different BOFED points and it indicates possible BOFED equivalent point could be conducted in future research. Figure 14 Figure 14: Panel (a) represents the 4th Scenario with Ray plot and the BOFED point (x1 = 10, x2 = 1). Panel (b) shows the fourth Scenario Heat map of the possible design points and panel (c) shows the Standard error around all possible design points a.

Third Scenario a. Third Scenario b. Design points heat map b. Design points heat map The scenarios and plots represent the location of the BOFED points. There is variation in the result at each scenario that is related to the shape of the TR that affects the BOFED points. When the TR is away from the origin the variation increases around the BOFED points as shown in the heat maps. The follow-up experimental design points are different in each scenario and have various possibilities of BOFED points. Each case represents a different BOFED point and that is due to the different data generated for each case. We could report the BOFED points to researchers to help them continue their experiment and learn about their subject by dosing at the BOFED points and have a stable tolerable region.

Discussion

AND CONCLUSION In this research, we proposed a novel method to determine the BOFED points consider- ing a toxicity experiment. Starting with a methodology to analyze datasets with multiple outcomes and multiple stressors in the dose-response setting. We estimated the tolera- ble region using MCMC techniques. Then, we find the BOFED points for future runs by minimizing the variance of the TR over the posterior predictive distribution in the fu-ture dose for a new run. This method uses the notation of area or volume to define the BMDTR for multiple chemicals. The idea discards outcomes that do not contribute to the BMDTR systematically from future consideration. We considered different possible exper- iment scenarios and determine the tolerable region for each scenario and define the results of the BOFED points in each scenario.

Researchers learn about all possible BOFED points using optimality criteria that determined the best follow-up design points, making it easier for researchers to produce more runs. That approach showed us the different impacts each scenario has on the BOFED points. This novel approach was developed using JAGS with computational effort. Computationally it takes more than 100 hours per run due to the 5 MCMC chains used, 10,000 MCMC samples taken, and 1,000 burn-in samples. Then the chains were thinned by 10 and 1,000 samples were used for inferences. Diagnostics were checked to ensure convergence and quality of MCMC samples. Thus, the optimal design criteria face some challenges in time when calculating the integral of the area, and defining the optimal follow-up design points, it spent a lot of time since it counts for all M MCMC samples.

As a future research, investigating different dose-response models by adding the Binomial distribution will extend our method to cover different distributions. One can ma- nipulate the properties of the binomial model to get an optimal fit, and find how it affects the choice of the BOFED. Applying this method to a real-life problem can help researchers to determine their next follow-up experiment to their study. This method could be applied to nonparametric models and investigate the BMDTR and the BOFED. Acknowledgement This research was funded by Princess Nourah bint Abdulrahman University and Researchers Supporting Project number (PNURSP2023R346), Princess Nourah bint Abdulrahman Uni- versity, Riyadh, Saudi Arabia.

References

1. C Andrieu, A Doucet, and CP Robert. Computational advances for and from bayesian analysis. Statistical science, pages 118–127, 2004.
2. AC Atkinson. The usefulness of optimum experimental designs. Journal of the Royal Statistical Society: Series B (Methodological), 58(1):59–76, 1996.
3. Anthony Atkinson, Alexander Donev, and Randall Tobias. Optimum experimental designs, with SAS, volume 34. Oxford University Press, 2007.
4. Edward L. Boone, J. Paul Brooks, Epiphanie Nyirabahizi, and Xi Chen. A bayesian model-based approach for determining multivariate tolerable regions. Journal of En- vironmentl Statistics, 7(6):??–??, 11 2015. ISSN 1945- 1296. URL http://jes stat.ucla.edu/v07/i06.
5. Kathryn Chaloner and Isabella Verdinelli. Bayesian experimental design: A review. Statistical Science, 10(3):273–304, 1995.
6. Kathryn Chaloner et al. Optimal bayesian experimental design for linear models. The Annals of Statistics, 12(1):283–300, 1984.
7. Chu-Chih Chen and Yin-Han Wang. Benchmark dose calculation for ordered cate- gorical responses with multiple endpoints. Stochastic Environmental Research and Risk Assessment, pages 1–9, 2018.
8. Sylvie Chevret. Statistical Methods for Dose-Finding Experiment. 06 2006. doi: 10.1002/0470861258
9. Yunfei Chu and Juergen Hahn. Integrating parameter selection with experimental design under uncertainty for nonlinear dynamic systems. AIChE Journal, 54(9):2310– 2320, 2008.
10. Merlise Clyde and Kathryn Chaloner. The equivalence of constrained and weighted designs in multiple objective design problems. Journal of the American Statistical Association, 91(435):1236–1244, 1996. Research Article
11. Merlise A Clyde. Bayesian optimal designs for approximate normality. 1994.
12. Merlise A Clyde. Experimental design: A bayesian perspective. International Ency- clopia Social and Behavioral Sciences, 2001.
13. Christopher Cox. Threshold dose-response models in toxicology. Biometrics, pages 511–523, 1987.
14. Kenny S Crump. A new method for determining allowable daily intakes. Toxicologi- cal Sciences, 4(5):854–871, 1984.
15. US EPA. Benchmark dose technical guidance document. Technical Report# EPA/100/R-12/001, 2012.
16. Naha Farhat. Design and analysis of toxicological experiments with multiple end- points and synergistic and inhibitory effects. 2014.
17. Naha J Farhat, Edward L Boone, and David J Edwards. A new method for deter- mining the benchmark dose tolerable region and endpoint probabilities for toxicology experiments. Journal of Applied Statistics, pages 1–29, 2019.
18. Valerii Fedorov. Theory of Optimal Experiments Designs. 01 1972.
19. Andrew Gelman, Hal S Stern, John B Carlin, David B Dunson, Aki Vehtari, and Donald B Rubin. Bayesian data analysis. Chapman and Hall/CRC, 2013.
20. Chris Gennings, W Hans Carter Jr, Michelle Casey, Virginia Moser, Richard Carch- man, and Jane Ellen Simmons. Analysis of functional effects of a mixture of five pesticides using a ray design. Environmental Toxicology and Pharmacology, 18(2): 115–125, 2004.
21. John P Groten, Victor J Feron, and Ju¨rgen Su¨hnel. Toxicology of simple and complex mixtures. Trends in pharmacological sciences, 22(6):316–322, 2001.
22. T Guest and A Curtis. Iteratively constructive sequential design of experiments and surveys with nonlinear parameter-data relationships. Journal of Geophysical Re- search: Solid Earth, 114(B4), 2009.
23. Anastasia Ivanova, Ken Liu, Ellen Snyder, and Duane Snavely. An adaptive design for identifying the dose with the best efficacy/tolerability profile with application to a crossover dose-finding study. Statistics in medicine, 28(24):2941–2951, 2009.
24. Marek Jakubowski and Sławomir Czerczak. A proposal for calculating the no- observed-adverse-effect level (noael) for organic compounds responsible for liver tox- icity based on their physicochemical properties. International journal of occupational medicine and environmental health, 27(4):627–640, 2014.
25. DP Lunn, N Davis-Poynter, MJBF Flaminio, DW Horohov, K Osterrieder, N Pusterla, and HGG Townsend. Equine herpesvirus-1 consensus statement. Journal of Veterinary Internal Medicine, 23(3):450–461, 2009.
26. Stephanie L Meadows, Chris Gennings, W Hans Carter Jr, and Dong-Soon Bae. Ex- perimental designs for mixtures of chemicals along fixed ratio rays. Environmental Health Perspectives, 110(Suppl 6):979, 2002.
27. VC Moser, M Casey, A Hamm, WH Carter Jr, JE Simmons, and C Gennings. Neuro- toxicological and statistical analyses of a mixture of five organophosphorus pesticides using a ray design. Toxicological Sciences, 86(1):101–115, 2005.
28. Peter Mu¨ller. Simulation based optimal design. Handbook of Statistics, 25:509–518, 2005.
29. MM Mumtaz, CT De Rosa, J Groten, VJ Feron, H Hansen, and PR Durkin. Esti- mation of toxicity of chemical mixtures through modeling of chemical interactions. Environmental health perspectives, 106(suppl 6):1353– 1360, 1998.
30. AndrejPa´zman.FoundationsofOptimumExperimental Design. 1986.
31. Friedrich Pukelsheim. Optimal design of experiments. SIAM, 2006.
32. S.D. Silvey. Optimal design. International Journal for Uncertainty Quantification, 1980.
33. LeAnna G Stork, Chris Gennings, Walter H Carter Jr, Robert E Johnson, Darcy P Mays, Jane Ellen Simmons, Elizabeth D Wagner, and Michael J Plewa. Testing for ad- ditivity in chemical mixtures using a fixed-ratio ray design and statistical equivalence testing methods. Journal of Agricultural, Biological, and Environmental Statistics, pages 514–533, 2007. Research Article
34. Robert K Tsutakawa. Design of experiment for bioassay. Journal of the American Statistical Association, 67(339):584–590, 1972.
35. Robert K Tsutakawa. Selection of dose levels for estimating a percentage point of a logistic quantal response curve. Journal of the Royal Statistical Society: Series C (Applied Statistics), 29(1):25–33, 1980.
36. Isabella Verdinelli. Advances in bayesian experimental design. Bayesian Statistics, 4:467–481, 1992.
37. Kai Wang, Feng Yang, Dale W Porter, and Nianqiang Wu. Two-stage experimental design for dose–response modeling in toxicology studies. ACS sustainable chemistry & engineering, 1(9):1119–1128, 2013. Research Article

This is a text version generated from the article. For the formatted version of record (with original tables & figures), download the PDF →