In pharmacological and toxicological dose-response experiments, a primary and important goal is to identify a threshold value for a specific purpose. For example, identifying the minimum effective dose (MED) of a medicine in the beginning stage of a dose-response experiment is important, as any dose level that is less than the identified level does not need to be considered in further studies. As another related example, the identification of the noobserved-adverse-effect-level (NOAEL) is equally important. It is because in general, as the dosage increases, the possibility of having adverse effects also increases. Therefore, once NOAEL is identified, one can minimize the adverse side effects from the medicine and also the cost of the experiments (Bretz et al. 2008;Guo and Li 2014;Kim et al. 2016;Schwartz et al. 2001;West and Kodell 2005;Coffey and Gennings 2007). In this context, the detection of MED or NOAEL can be considered as a threshold estimation problem (e.g., Mallik et al. 2011). In addition to dose-response studies, threshold estimation problems can also be found in environmental studies such as detecting changes in air pollutant levels (Holst et al. 1996).
Existing statistical methods for threshold estimation in dose-response experiments can be broadly grouped into two different categories (Bretz et al. 2008): multiple comparison methods (Williams 1972;Ruberg 1989;Hsu and Berger 1999;Tamhane and Logan 2002;Guo and Li 2014;Otava et al. 2017) and model-based methods (Cox 1987;Muggeo 2003;Schwartz et al. 2001;West and Kodell 2005;Coffey and Gennings 2007;Lutz and Lutz 2009;Mallik et al. 2011), while multiple comparison methods can be further divided into Bayesian (e.g., Guo and Li 2014;Otava et al. 2017) and frequentist (e.g., Williams 1972;Ruberg 1989;Hsu and Berger 1999;Tamhane and Logan 2002) approaches.
Multiple comparison methods require fewer assumptions, but they restrict the estimated MED to be one of the investigated levels, i.e., they will not produce any estimated MED value that is different from those dose levels that were experimented with-not even with interpolation between doses levels or similar techniques. On the other hand, model-based methods are capable of interpolating and extrapolating dose levels as they impose a stronger assumption that there exists a model that adequately captures the relationship between the dose levels and the response. Most drug development experiments assume that such a relationship is non-decreasing, which means the higher the dosage the stronger the therapeutical effects (Bretz et al. 2008;Lutz and Lutz 2009;Mallik et al. 2011). Of course, if the chosen model does not adequately capture the relationship, the estimated threshold value is unlikely to be close to the true value. This is a potential drawback of the model-based methods.
Both Bayesian (e.g., Guo and Li 2014;Otava et al. 2017) and frequentist (e.g., Williams 1972;Hsu and Berger 1999) multiple comparison methods typically define the hypotheses of interest according to the investigated dose levels. For example, if there are two dose levels, four hypotheses will be investigated: (i) none of the doses have any effect, (ii) only one dose has effect, (iii) both doses have the same effect, and (iv) the two doses have different effects. As the number of investigated dose levels increases, the number of hypotheses increases rapidly. Therefore, such methods (e.g., Williams 1972;Hsu and Berger 1999;Guo and Li 2014;Otava et al. 2017) only considered the case when the number of dose levels is small (4 or 5). The Bayesian methods typically assign prior probabilities to the hypotheses and make final conclusions using the posterior probabilities.
For model-based methods, if one believes the relationship between the dosage and response is monotonic, one can model this problem with a regression function that is constant to the left of a certain covariate value t 0 , say, and is increasing to the right of t 0 (Bretz et al. 2008;Lutz and Lutz 2009;Mallik et al. 2011). Under this regression model, the main purpose is to make an inference about the value of t 0 , the threshold value. In dose-response experiments, it can be a MED, which starts showing a discernible effect compared to a control. A notable example is the break point estimation method developed by Muggeo (2003).
With a parametric segmented regression model commonly used in the threshold estimation problem, this method first detects if a threshold t 0 exists, and if yes, it then provides the maximum likelihood estimate (MLE) as well as a confidence interval for t 0 . Consequently, this method does not always produce an estimate for t 0 . More recently, in Mallik et al. (2011) a method based on a novel p value framework was developed. The idea is to first conduct a sequence of tests that, for each covariate value, assesses whether the regression function is at its baseline, then calculate the corresponding p values, and lastly estimate t 0 by fitting a piecewise function to those calculated p values. This method, however, does not provide confidence intervals.
The main contribution of this paper is a new method for performing statistical inference for t 0 . This method is model-based and has the following properties:
• It is neither Bayesian nor frequentist; rather, it follows the generalized fiducial inference methodology of Hannig et al. (2016).
• Model averaging is employed to avoid the above-mentioned danger of selecting a wrong model.
• Unlike most existing methods, it provides point estimates as well as confidence intervals for t 0 .
To be more specific, four different parametric forms are used for model averaging: the first one is a step function with its jump point at t 0 , the second one begins with a constant function until it hits t 0 and then has a linear trend afterward, while the third and fourth are similar, except they have, respectively, a quadratic and a cubic trend after t 0 ; some of these parametric forms were used by previous authors (more below). The generalized fiducial inference methodology is used to compute a probability density function for these four parametric forms. With such a density function, the proposed method pools multiple threshold estimates from each of the four parametric forms to form a single point estimate for t 0 , as well as constructs a confidence interval for its value. The rest of this paper is organized as follows. Section 2 provides a brief introduction of generalized fiducial inference, while Sect. 3 states the problem formulation. In Sect. 4, the empirical performance of the proposed method is illustrated under various experimental settings, while in Sect. 5 the proposed method is applied to two real datasets. Lastly, concluding remarks are offered Sect. 6.
Bayesian inference is an important statistical methodology that offers both point and uncertainty estimates for the unknown parameters via the calculations of posterior distributions. This provides more information compared to frequentist inference (Xie and Singh 2013). However, when no prior information is available, any naive application of the Bayesian methodology could cause concerns; see Efron (2013) for an example in which different priors lead to very different results. To address this issue, Fisher (1930) introduced fiducial inference, where the main idea is to transfer the randomness from the observed data into the uncertainties of the parameter estimates via the use of a switching principle. This switching principle is also used in maximum likelihood estimation and will be described further below. The uncertainties of the parameter estimates are summarized in the form of a distribution function, termed generalized fiducial density, which plays a similar role as the posterior distribution in the Bayesian setting.
Soon it was realized by other researchers that Fisher's fiducial idea only works for univariate problems, and hence, it did not receive much attention from the statistics community for a few decades. However, there has been a recent resurgence of interest in different variants of Fisher's idea, including Dempster-Shafer theory (Dempster 2008), inferential models (Martin and Liu 2015), confidence distribution (Xie and Singh 2013), and more generally fusion learning Cheng et al. (2014). All these variants use similar but different ideas to extend Fisher's idea to multivariate cases. The particular variant that this paper considers is the so-called generalized fiducial inference (Hannig et al. 2016).
Generalized fiducial inference (GFI) repairs many drawbacks of Fisher's original proposal and has been applied to various inference problems with great success, including wavelet regression, ultrahigh-dimensional models, nonparametric survival function estimation; e.g., see Hannig and Lee (2009), Lai et al. (2015), Hannig et al. (2016), Cui and Hannig (2019) and references therein. GFI assumes that the data are generated from a data generating equation G and expresses the relationship between the data Y and the parameters θ as
where U is the random component whose distribution is assumed to be completely known. For many practical problems, it is the distribution of the (normalized) residuals. For the current problem, it is the distribution of u in (4) below, which is assumed to be N (0, 1). Similar to Fisher's maximum likelihood method, an important idea behind GFI is the so-called switching principle. That is, the roles of the random data Y and the deterministic parameters θ are switched after Y is observed: the observed data, denoted as y, are now treated as deterministic while θ is treated as random. By using (1) and this switching principle, we can derive a probability distribution of θ as follows.
Suppose for now, the data generating equation G is invertible for any u ∈ U and any observed y. That is, the inverse G -1 always exists:
(2)
Recall that the distribution of U is completely known, thus one can generate a random sample {u 1 , u 2 , . . . , u n } of U and plug them into (2) and obtain
We shall call { θ1 , . . . , θ n } a fiducial sample of θ , which can be used to form point estimates and confidence intervals for θ, as similar to a posterior sample in the Bayesian context. One can see that through G -1 , the randomness in the data y is transformed into the randomness in the parameter θ .
Of course in practice, the inverse function G -1 in (2) does not always exist, and Hannig (2013) proposed a solution to deal with this issue.
From the above discussion, one can see that there is a density function behind the fiducial sample { θ1 , . . . , θ n }. We shall call this density r (θ|y) the generalized fiducial density. When θ ∈ ⊂ R p and y ∈ R n , and under some regularity conditions, it can be shown that r (θ|y) admits the following expression (Hannig 2013;Hannig et al. 2016):
where f (y, θ ) is the likelihood, J (y, θ ) is a Jacobian-like quantity that admits the expression:
with
,
where A ∈ R n× p . The derivation of (3) is not straightforward, but it was obtained by applying the implicit function theorem and Jacobian transformation. Note also its similarity with the Jeffreys prior in the Bayesian methodology.
Depending on the problem settings, it may not be feasible to obtain closed-form expressions for J (y, θ ) or f (y, θ )J (y, θ )dθ . Especially, with the large sample size, the number of combinations of tuples is too large to consider all possible combinations. This problem can sometimes be solved by sampling the combinations of tuples and taking the mean of them (e.g., Hannig et al. 2014). Some other times, one could use Markov Chain Monte Carlo (MCMC) methods to generate fiducial samples without the need to evaluate the normalizing constant f (y, θ )J (y, θ )dθ . In fact, as described below, for the present problem we need to use MCMC to generate a fiducial sample.
As mentioned before, many dose-response studies assume a monotonic relationship between the dose levels x and the pharmacological effects y of the medicine (e.g., Bretz et al. 2008Bretz et al. , 2010;;Pinheiro et al. 2006;Mallik et al. 2011). Hence, one can consider the following regression model:
where t 0 is the threshold, b 0 ∈ R is the baseline, u is the error term (usually standard normal), and σ is the noise standard deviation. Also, we use the following to denote the observed data from a study:
where y i j is the measured pharmacological effect of i-th dosage from subject j, x i is the i-th dosage (x 0 represents control), n is the number of dosages, m i is the number of subjects with i-th dosage, and u i j iid ∼ N (0, 1) The goal is to conduct statistical inference on t 0 , which is the MED in this dose-response setting. As mentioned in the introduction, we assume a constant value (β 0 ) for the first portion of μ(x) (i.e., μ(x) for x ≤ t 0 ) while we use four different non-decreasing functions to model μ(x) for x > t 0 :
where I E is the indicator function for the event E.
At least for the cases of p = 1 and 2, these functions have been used by previous authors (Muggeo 2003; Mallik et al. 2011) to model μ(x). We incorporate p = 3 and 4 with the expectation that ( 5) is rich enough to capture the behaviors of the responses in most practical situations.
To proceed with inference on t 0 , we need to calculate the generalized fiducial density r (θ|y) for the above setup, where θ = {t 0 , p, σ, β 0 , β 1 } and y = (y i j ; ∀i, j).
Due to the indicator function, the four regression functions above do not satisfy a smoothness assumption required by Hannig (2013), and therefore (3) cannot be applied. To overcome this issue, the indicator function is replaced by a sigmoid function:
where C is a large constant that we set it to be C = 10 3 . In the sigmoid function, λ controls the slope near t 0 . See Figure S.1 for some examples of (5) with the indicator function replaced by the sigmoid function.
The data generating equation G for the current problem is:
for p = 1, 2, 3, 4. Now with this data generating equation ( 7), one could attempt to calculate the generalized fiducial density r (θ|y) using (3). However, the analytic evaluation of the normalizing constant f (y, θ )J (y, θ )dθ in (3) is infeasible. Consequently, we need to resort to MCMC methods for generating fiducial samples for θ , as to be described below.
This subsection presents an MCMC algorithm for generating a fiducial sample for θ = {t 0 , p, σ, β 0 , β 1 }. We first provide a summary of the algorithm, then followed by the details. The steps of the algorithm are:
1. generate a t0 using (8); 2. given t0 obtained in the previous step, generate a p from (9); 3. given t0 and p obtained in the previous steps, generate a σ from (10); 4. given t0 , p and σ obtained in the previous steps, generate β0 and β1 using (11); 12) is larger than a realization from Uniform(0, 1); otherwise, set θ new = θ old . Now, we present the detailed steps. Generating t 0 : To generate a value for t 0 , we use the truncated normal distribution with x 0 and x n as, respectively, the lower and upper bounds:
where t 0 old is the t 0 value from the previous iteration, and σ 2 t is used to perturb its value. In our practical implementation, we set σ t = (x nx 0 )/10 and chose the initial value of t 0 at random uniformly from [x 0 , x n ].
Generating p given t 0 : Suppose now we have a value for t 0 . We can calculate the conditional generalized fiducial density for each of the four models μ p (x) in (5), indexed by p = 1, 2, 3, 4. Let ξ = {β 0 , β 1 , σ } ∈ and denote the conditional generalized fiducial density as r ( p|t 0 , y):
Direct calculations give
where X p is the design matrix for μ p (x), RSS p is the corresponding residual sum of squares, and N is total number of observations (N = i m i ). Therefore,
Thus, given a value for t 0 , using (9) one can calculate the fiducial probabilities for the four models in (5), and from these probabilities, one can generate a model out of the four, or equivalently, generate a value for p. Generating σ 2 given (t 0 , p): Next, given t 0 and p, it can be shown that the conditional generalized fiducial distribution of σ 2 is inverse χ 2 :
where
Therefore, σ 2 can be sampled from (10).
Generating β p = (β 0 , β 1 ) given (t 0 , p, σ ): It is straightforward to show that the generalized fiducial distribution of β p condition on t 0 , p, σ 2 is
from which a value for β p can be sampled from.
Accept/Reject: Executing the above steps will produce a value for θ = (t 0 , p, σ, β 0 , β 1 ) which we need to decide if it should be accepted or rejected. To do so, we need the generalized fiducial density r (θ|y):
where
.
To speed up the computations, we use the sampling technique proposed in Hannig et al. (2014) to approximate the summation i | det(A p ) i |.
Denote those new parameter values generated using the above steps as θ new , and those parameter values in the previous MCMC iteration as θ old . Let T (θ new |θ old ) be the proposal distribution which is the product of (10), ( 11) and ( 8). If the ratio
is larger than a realization from Uniform(0, 1), we accept θ new . Otherwise, we use θ old again for the next iteration.
Repeating the above procedure, a fiducial sample for θ can be generated, and the methods implemented by Plummer et al. (2006) can be applied to determine if the sample size is large enough. Once a sufficiently large fiducial sample is available, statistical inference can then be conducted. For example, one could use the average of all t0 's as an estimate of t 0 . Also, for a 95% confidence interval, one could take the quantiles 2.5% and 97.5 % of the t0 's.
Lastly, we stress that in our approach more than one model (see ( 5)) are used to generate the above fiducial sample t0 's, and therefore reducing the risk of incorrect model selection.
This section studies the empirical performance of the proposed method through a sequence of numerical experiments. In particular, the proposed method is compared to two existing methods in the literature: the MLE break point estimation method of Muggeo (2003) and the p value method of Mallik et al. (2011). Further details of these two methods can be found in the introduction. In sequel, we shall refer the proposed method as GFI, the method of Mallik et al. (2011) as p value, and the method of Muggeo (2003) as MLE. The code for GFI can be downloaded from https://github.com/vic-dragon/GFI_threshold Altogether five different test functions were used for μ(x):
• Test Function 1: μ(x) = 0.5 × I {x≥0.5} ,
• Test Function 2: μ(x) = 2(x -0.5) 2 × I {x≥0.5} ,
• Test Function 4: μ(x) = 0.05 sin(50x)I {x<0.5} + 2(x -0.5) 2 × I {x≥0.5} ,
• Test Function 5: μ(x) = 0.1 -0.1(x -0.5)I {x<0.5} + 0.5(x -0.5) × I {x≥0.5} .
These functions are displayed in Fig. 1. Notice that only the first two functions satisfy the parametric assumption (5). Also notice that all five functions share the same threshold value t 0 = 0.5, although their shapes are quite different. However, for Test Function 5, one should treat t 0 = 0.6 as the threshold value (instead of 0.5), as it is the point that starts showing a discernible effect compared to the baseline (Agathokleous et al. 2019). Some of these test functions have been used by other authors (e.g., Muggeo 2003; Mallik et al. 2011). The data were generated as follows. First, we set x i = i/(n + 1) for i = 1, . . . , n. Then, we calculated y i j = μ(x i ) + i j , where μ(x) is one of the five test functions, i j iid ∼ N 0, σ 2 , for j = 1, . . . , m and i = 1, . . . , n. Two values were used for σ = {0.1, 0.3}, while five different combinations of values were used for (m, n) = {(6, 5), (10, 10), (10, 20), (20, 50), (50, 100)}. We also investigated the performances of the methods under heteroscedastic noise: we generated data with the gradually increasing noise variance (σ + 0.2(x -0.5)I {x≥0.5} ).
For each combination of experimental parameters, 200 data sets were generated. For each of these data sets, the three methods mentioned above were applied to estimate the threshold value t 0 = 0.5 for Test Functions 1 to 4 and t 0 = 0.6 for Test Function 5. The averaged biases and mean-squared-errors (MSEs) were calculated. Also, for the GFI and MLE methods, we constructed 95% confidence intervals for the threshold value. The empirical coverage rates (ECRs) and their averaged widths of these confidence intervals were calculated. The results can be found in Tables 1,2, S.1, and S.2. Note that the p value method does not produce confidence intervals and hence, is omitted for this comparison.
With the largest sample size that we tested (n = 50, m = 100), the GFI method typically took less than 2 min to process one data set with a Dell PowerEdge R360 server, 2015 model. The acceptance rates of the MCMC algorithm depend on the experimental configurations. They ranged from 0.022 to 0.572 for homoscedastic noise, and from 0.024 to 0.559 for heteroscedastic noise. More detailed information can be found in Tables S.3 and S.4.
Recall that the MLE method does not always produce an estimate t0 for t 0 (see Sect. 1). Thus, all the results in the previous tables are based on those cases where t0 were available. The percentages of time that MLE produced estimates are given in Tables S.5 and S.6. These percentages ranged from 42% to 100%, with most cases larger than 90%.
From the experimental results, the following empirical conclusions can be drawn. First, in terms of point estimation, no method is uniformly the best. The relative performances of the three methods depend largely on the test function. Nevertheless, it seems GFI provided the best results for a majority of the experimental configurations (72% of the experimental configurations for homoscedastic noise and in 66% of the experimental configurations for heteroscedastic noise), and that p value suffered from bias issues. For confidence intervals, one cannot conclude if GFI or MLE is uniformly better, although GFI gave better results in 88% of the experimental configurations for homoscedastic noise and in 96% of the experimental configurations for heteroscedastic noise. Therefore, in practice, our recommendation is that, if nothing is known about the shape of the underlying function, GFI would be the first method to apply.
Recall that the proposed GFI method uses four different parametric functions to model the underlying true function (see ( 5)), and that the fiducial samples for p were sampled from the conditional distribution (9). Figures S.2 to S.5 report the relative percentages of the sampled values of p. As Test Functions 1 and 2 are special cases of (5), one can see that when m and n are large, GFI has a high tendency of selecting the true model.
In this section, two real data sets are analyzed with the proposed GFI method. We note that for one data set the monotonicity is in the opposite direction to model (5), while for the Reported are the averages (×10 2 ) and standard errors (×10 2 , in parentheses) of the biases and the MSEs for the threshold estimates obtained by the GFI, MLE, and p value methods (smallest values are bolded). Also reported are the empirical coverage rates (ECRs) and the average widths (in parentheses) of the 95% confidence intervals produced by the GFI and MLE methods (ECRs that are closest to the nominal target rate of 95% are bolded). Recall that n is the number of dosages and m is the number of subjects for each dosage The first data set is related to a dose response experiment originated from a biological study of the effect of the 1-methyl-3-butylimidazolium tetrafluoroborate treatment on the IPC-81 rat cell lines; see Ranke et al. (2004). The purpose of this study is to assess environmental hazards through measuring toxicity of ionic liquids on the mammalian cell lines. It is of particular interest to estimate the lethal dose level which makes cell lines stop responding. The cell responses were measured as percentages of cell viability compared to control, taken at different dose levels (measured in μM) of the treatment. The sample size is N = n × m = 9 × 9 = 81. The GFI method provides a threshold value estimate of t0 = 5.22 μM, with a 95% confidence interval as (4.92 μM, 5.34 μM). The MLE method provides a threshold value estimate of t0 = 5.10 μM, with a 95% confidence interval as (4.58 μM, 5.62 μM), while the p value method provides a threshold value estimate of t0 = 5.52 μM. The observed data, as well as these results, are displayed in the left panel of Fig. 2. One can see that, after the GFI estimated threshold 5.22 μM, the cell culture dose not show any significant response. Also, the GFI point estimate for t 0 is not in the middle of the confidence interval. Instead, it is closer to the right side, which is very reasonable when inspecting the left panel of Fig. 2. Lastly, as 100% of the fiducial sampled values from p = 2, it suggests that the underlying relationship can be well modeled with a piecewise linear function.
The second data set concerns a light detection and ranging (LiDAR) experiment. LiDAR is a technique for measuring the concentrations of various atmospheric particulates. Very often, it is used to monitor the air pollution status in an area. This second data set contains atmospheric atomic mercury measurements taken at the geothermal power plan Bella Vista in Italy (e.g., Holst et al. 1996;Mallik et al. 2011). There are 221 observations. The predictor is the distance from the measuring equipment to the source that reflected the light, with values ranging from 390 m to 720 m. The response is the logarithm of the received light ratio from two different wavelengths; see the right panel of Fig. 2. The decreasing trend starting around range=540 m indicates the existence of mercury. And it is of interest in knowing exactly when this decreasing trend started. The GFI method gives a threshold value estimate of t0 = 543.55 m, with a 95% confidence interval as (541.47 m, 548.92 m). The MLE method provides a threshold value estimate of t0 = 544.66 m, with a 95% confidence interval as (538.04 m, 551.29 m), while the p value method provides a threshold value estimate of t0 = 541 m. These results are also displayed in the right panel of Fig. 2. As similar to above, the point estimate from GFI is not in the middle of the confidence interval. In addition, 24% and 76% of the fiducial sampled values for p are 1 and 2, respectively. This suggests that the underlying relationship may not be simply modeled by either a piecewise step function or linear function.
From Fig. 2, for both real data examples, one can see that the GFI and MLE methods gave very similar point estimates for t 0 , while the GFI confidence intervals are narrower than those from MLE. The point estimates from the p value method are outside the GFI confidence intervals: it seems to suggest that the p value estimates are biased, which are in agreement with the simulation results.
In this paper, we developed a method for estimating threshold under the a general setting used in many biomedical studies such as toxicology and dose-response studies. The proposed method is based on the generalized fiducial inference framework. Two noteworthy features of the proposed method are (i), in addition to point estimates, it also provides confidence intervals for the parameters of interest, and (ii) it allows the possibility of having multiple parametric models to capture the underlying pattern of the data so that the risk of model mis-specification is reduced. Results from numerical experiments and applications to real data suggest that the proposed method possesses excellent empirical properties. Some extensions of the current work are worth considering. One possibility would be to adopt a richer set of non-decreasing functions to model μ(x), as opposed to (5). This could be done quite straightforwardly. Another possibility would be to allow the occurrence of multiple change points, where the number of the change points is unknown a priori. However, this extension could cause some complications, as for example, one potentially would need to handle the issue that different fiducial samples could give different number of change points. At the moment, if there is no change point, our numerical experience suggests that our method will often return x 1 as the estimated threshold, while its performance is less predictable if there are multiple change points.
One could also investigate the theoretical properties of the proposed method. For example, if one is willing to make the unrealistic assumption that t 0 is given, one can show that as the sample size increases, the probability of selecting the correct model from (5) goes to 1, which leads to the following proposition. Gao et al. (2020) and that t 0 is given. Denote the true value of p in (5) as p 0 . Then, the fiducial probability (9) follows r ( p 0 |t 0 , y) → 1 as n → ∞. This proposition can be proved following the arguments of Gao et al. (2020). Essentially, it involves showing that RSS p /RSS p 0 → 0 for any p = p 0 , which leads to r ( p 0 |t 0 , y) = → 1 as n → ∞.
With Proposition 1, the results of Hannig et al. (2016) also guarantee that inferences made by GFI enjoy asymptotically correct frequentist properties, i.e., confidence intervals constructed by r (θ|y) will achieve the nominated coverage levels when the sample size is large enough. However, the case of unknown t 0 is substantially more complicated and we plan to report our results in a future paper.