{'fraction_non_alphanumeric': 0.04410881714470463, 'fraction_numerical': 0.024663987806567825, 'mean_word_length': 4.1540109497738635, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Benchmark Dose Estimation using a Family of Link Functions

This article proposes a method of estimating benchmark dose (BMD) using a family of link functions in binomial response models dealing with model uncertainty problems. Researchers usually estimate the BMD using binomial response models with a single link function. Several forms of link function have been proposed to fit dose response models to estimate the BMD and the corresponding benchmark dose lower bound (BMDL). However, if the assumed link is not correct, then the estimated BMD and BMDL from the fitted model may not be accurate. To account for model uncertainty, model averaging (MA) methods are proposed to estimate BMD averaging over a model space containing a finite number of standard models. Usual model averaging focuses on a pre-specified list of parametric models leading to pitfalls when none of the models in the list is the correct model. Here, an alternative which augments an initial list of parametric models with an infinite number of additional models having varying links has been proposed. In addition, different methods for estimating BMDL based on the family of link functions are derived. The proposed approach is compared with MA in a simulation study and applied to a real data set. Simulation studies are also conducted to compare the four methods of estimating BMDL.


Introduction
One of the main goals in quantitative risk assessment is to estimate the risk function R(d), which is the probability of adverse events, such as death, birth defect, weight loss, cancer or mutation exhibited in a subject exposed at dose level d. Suppose n number of subjects are exposed to a dose level d and y number of adverse events are observed. Then, the response y is distributed according to a binomial distribution with parameter [n, R(d)], where R(d) is the probability of adverse events at dose level d. Methods of estimating BMD and BMDL are discussed by several researchers such as Crump (1984); Bailer et al. (2005); Morales et al. (2006); Bailer (2007, 2009) ;West et al. (2012) to name just a few. Crump (1984) introduced methods of estimating BMD and BMDL by proposing four models for discrete responses and three models for continuous responses. There are eight models Bailer, 2007, 2009;West et al., 2012) that have been identified as standard models for estimating BMD and BMDL. One of the models from the set of standard models may be chosen for fitting the data sets. However, the responses may be generated from the model other than the chosen model. Researchers Bailer, 2007, 2009;West et al., 2012) have shown that the estimation of BMD and BMDL are significantly effected if the assumed model is incorrect. So, there is a recent rise in developing methods of accounting for model uncertainty in BMD estimation.
For accounting model uncertainty in BMD and BMDL estimation, model averaging (MA) methods are proposed by Kang et al. (2000); Bailer et al. (2005); Wheeler and Bailer (2007); Shao and Small (2011); West et al. (2012); Piegorsch et al. (2013). The estimates of BMD and BMDL using model averaging methods are given by the weighted average of the estimates of BMD and BMDL using individual models belong to a set of models. Bayesian methods and Bayesian model averaging methods for estimating BMD and BMDL are also proposed by Morales et al. (2006); Shao and Small (2012); Simmons et al. (2015). The model averaging approach may solve the problems of model uncertainty, when the true model generating responses can be approximated by some members of the model space containing the assumed models. Since, the model space are always finite, there may be infinite number of other models which can not be approximated by the members of the model space. So, model averaging techniques provide a partial solution to the problem of model uncertainty.
Here, a family of link functions containing some of the standard link functions as well as infinite number of other link functions is used to fit the binomial response models. The family of link functions is parameterized by two unknown link parameters. There are infinite number of link functions can be represented by different values of link parameters. Some standard link functions correspond to some finite values of link parameters. So, we may get a better results for accounting model uncertainty in BMD and BMDL estimation using the proposed model.
The remainder of the article is organized as follows: in Section 2, the binomial response models using a family of link functions are discussed. An expression for BMD using the family of link functions is given in Section 2.2. An example with real data set is shown in Sections 3 to illustrate the proposed method of estimating BMD. Four methods of estimating BMDL are derived in Section 2.4 and a comparison study among the four methods are provided in Section 4.2. In Section 4.1, the proposed method is compared with model averaging method using simulation studies. Concluding remarks are given in Section 5.

Method
In this section, we discuss the binomial response models with a family of link functions and provide methods of estimating BMD and BMDL using the models.

Binomial Response Models
The Binomial Response Models are members of the Generalized Linear Models (GLMs) described by three components given below.
1. Distributional components: let y 1 , y 2 , . . . , y n be n random samples of adverse events at dose levels d 1 , d 2 , . . . , d n , where for each i ∈ {1, 2, . . . , n}, y i has binomial distribution with parameter (n i , r i ), r i ∈ [0, 1], andȳ i = y i n i has scaled binomial distribution belongs to the exponential family having the form of probability mass function (pmf) given by (Fahrmeir and Tutz, 2001) , where g is called parametric link function and h is the inverse of g. We usually assume that the inverse of g exists.
For dose-response studies, the linear predictor is usually assumed as η(d) = β 0 + β 1 d, or η(d) = β 0 + β 1 d + β 2 d 2 , and a single link function such as logistic, probit, log-log, complementary log-log or some other link functions are assumed to fit the models. Here, instead of a single link function, we are using a family of link functions (parametric link function) parameterized by a link parameter vector α to fit the models. So, we are denoting the link function as g(α, ·) in place of g(·).
For estimating the risk function using the above models, we need to estimate the unknown parameters using a available data sets. Let us denote δ = [β ′ , α ′ ] ′ for the combined parameter vectors including the unknown regression parameter vector β and the link parameter vector α. The unknown parameter vector δ can be estimated using the Maximum Likelihood Estimation (MLE) methods given in Stukel (1988). Due to estimation of the link parameters along with regression parameters, the variances of the estimated regression parameters are increased (Taylor, 1988). The variance inflations of the regression parameters are asymptotically zero if the link parameters are orthogonal to the regression parameters (Cox and Reid, 1987). Czado (1997) proposed some conditions on the family of link functions providing local orthogonality between link and regression parameter vectors. A family of link functions Λ = {h(α, ·) : α ∈ Ω} provides local orthogonality between link and regression parameter vectors around a point η 0 asymptotically, if the following conditions are satisfied.
1. There exists η 0 and r 0 such that and 2. There exists s 0 such that ∂h(α, η) ∂η where Ω is denoted for the parameter space of α. Such a family Λ = {h(α, ·) : α ∈ Ω} satisfying conditions (2.2) and (2.3) is called (r 0 , s 0 ) − standardized at η 0 (Czado, 1997). Now, for estimating risk function R(d) using (r 0 , s 0 )−standardized family at η 0 , we need to estimate extra three parameters r 0 , s 0 , and η 0 . For avoiding estimating extra three parameters, Czado (1997) proposed to choose r 0 = β 0 , s 0 = 1, and η 0 = β 0 . By choosing the values such a way, the variance inflations of β are reduced as the values of η vary around the point η 0 = β 0 , when centered covariates (Czado, 1997). For this, if dose levels are not centered, we need to transfer the available dose levels as x i = d i −d, and after estimating BMD/BMDL from the model, we make the inverse transformation to get the estimates of BMD/BMDL in the true range of dose levels. For constructing (r 0 = β 0 , s 0 = 1) − standardized family at η 0 = β 0 , we adopt the methodologies given by Czado (1997). Here, we use Stukel (1988)'s generating family to construct the family of link functions, and the (r 0 = β 0 , s 0 = 1) − standardized at η 0 = β 0 generating family is given by: using the binomial response model with (r 0 = β 0 , s 0 = 1) − standardized at η 0 = β 0 generating family is given by where η ≡ η(d), and G c (α, ·) is given by equations (2.4) & (2.5). In the next section, we provide a expression for the Benchmark dose (BMD) using the above model.

Benchmark Dose Estimation
The Benchmark dose (BMD) is defined by the dose level having extra risk R E (BMD) = BMR, where BMR is the Benchmark risk usually pre-specified as 0.01, 0.05, and 0.1. So, BMD is the solution of the equation where p 0 is the Background risk, i.e, the risk at the minimum dose level d 1 . We denote δ = [β ′ , α ′ ] ′ for the joint parameter vector including the regression parameter vector β, and the link parameter vector α. For a fixed value of BMR ∈ [0, 1], the BMD can be expressed as a function of δ, S(δ), say. From equations (2.6) and (2.7), we get a expression for BMD as where LBMR = log( BM RE 1−BM RE ), and I {LBM R≥β 0 } is the indicator function taking value 1 if LBMR ≥ β 0 , and 0 otherwise. The functions S 1 (δ) and S 2 (δ) are given by and, Note that we require centered dose levels (i.e,d = 1 n n i=1 d i = 0) for using the model (2.6). If the dose levels are not centered, we make the transformation x i = d i −d to have the centered dose levels.
After estimating BMD from the model we make the inverse transformation to get the estimated value of BMD within the true range of dose levels.

Asymptotic Results
The asymptotic distributions of unknown parameters for q dimensional multinomial response models with a family of link functions are discussed in Das and Mukhopadhyay (2014). For q = 1, we get the binomial response models using a family of link functions. So, the similar results can be applicable for binomial response models using a family of link functions. However, for making this article self contained, we provide the required asymptotic results here. We denoteδ for the log-likelihood function, and ∂l ∂δ for the score function for the observed responses. Also, J n is denoted for the Fisher's information matrix. The asymptotic results are given by the following Lemmas.
Lemma 1: The score function ∂l ∂δ has an asymptotic multivariate normal distribution with mean 0 and variance J n .
Proof: From Section 2, the risk function is given by, ( 2.9) where η(d) = f (d)β, β is an unknown regression parameter vector and α is a vector of unknown link parameters. Also, where g is the inverse of h. Now, from Section 2, the log-likelihood function for the sample y 1 , . . . , y n is given by (2.11) Thus, the score function is (Fahrmeir and Tutz, 2001, p 436), .12) and (Fahrmeir and Tutz, 2001, p From equation (2.13), we get the Fisher information matrix is (2.14) From equation (2.12), using the central limit theorem we have ∂l(δ) ∂δ has asymptotic normal distribution with mean 0 and variance J n .
Lemma 2: The MLE of δ,δ has an asymptotic multivariate normal distribution with mean δ and variance J −1 n .
Proof: By Taylor series expansion and approximating up to first order term, we have which gives (Fahrmeir and Tutz, 2001, p 439), Thus, the MLE of δ,δ has an asymptotic normal distribution with mean δ and variance J −1 n .
Lemma 3: The estimate BMD = S(δ) is a consistent estimator for BMD.
Proof: The proof is trivial from the result that the MLE of δ,δ is a consistent estimator of δ, and S(δ) is a continuous function of δ. Hence, S(δ) is a consistent estimator for S(δ), i.e., BMD is a consistent estimator for BMD.
In the next section, we provide confidence intervals for BMD to find BMDL from the proposed model.

Confidence Intervals
Here, we provide four methods of constructing confidence intervals for BMD for a particular value of BMR = BMR 0 . The methods are discussed as follows:

Confidence interval using ML estimates
Here, we use the asymptotic result of the distribution ofδ for constructing the confidence interval for BMD. From Lemma 2, we haveδ has an asymptotic multivariate normal distribution with mean δ and variance Σ = J −1 n . Hence, (δ − δ) ′ Σ −1 (δ − δ) has an asymptotic χ 2 -distribution with p degrees of freedom, where p is the order of the vector δ. Hence, the 100(1 − τ )% confidence region for δ is given by where χ 2 p,(1−τ ) , is the (1 − τ )th quantile of the χ 2 distribution with p degrees of freedom. For BMR = BMR 0 , we compute BMD = S(δ), for δ ∈ C using equation (2.8). Let us denote S L = Min{S(δ) : δ ∈ C}, and Hence, the 100(1 − τ )% conservative confidence interval for BMD is given by [S L , S U ].

Confidence interval using LR test
Here, we test the null hypothesis where R E (d) = R(d)−p 0 1−p 0 , with p 0 is the background risk. Let D(d) be the deviance (Fahrmeir and Tutz, 2001, p 108) under null hypothesis and D(d) be the deviance of the fitted model. Then, L(d) = D(d) − D(d) has an asymptotic χ 2 -distribution with 1 degree of freedom. Let us denote Then, the 100(1 − τ )% confidence interval for BMD is given by [L min , L max ].

Confidence interval using score test
Let us denote u 0 = ∂l ∂β 0 δ 0 , whereδ 0 is the MLE of δ under H 0 given in equation (2.17). Letσ 2 0 be the estimated variance of u 0 at δ =δ 0 . Then, T (d) = u 2 0 /σ 2 0 has an asymptotic χ 2 distribution with 1 degree of freedom (Fahrmeir and Tutz, 2001, p 48) Note that the above confidence intervals are by nature two sided. To get one sided confidence interval (BMDL), some researchers (Buckley et al., 2009;Nitcheva et al., 2005) proposed an adjustment by doubling the significance level of the test and then ignoring the upper limit. So, we construct 100(1 − 2τ )% two sided confidence intervals for BMD using the above methods and then the lower limits of that intervals are taken as the one sided 100(1 − τ )% lower confidence bound for BMD.
Since, all of the above confidence intervals are constructed using asymptotic results, here we provide a confidence interval using bootstrap technique.

A bootstrap lower confidence bound
For constructing bootstrap lower confidence bound, we generate l responses y k = [y 1k , y 2k , y 3k , y 4k ] ′ using the fitted modelr where y ik has binomial distribution with parameter (n,r i ) for i = 1, . . . , 4 and k = 1, . . . , l. From the generated l data sets, we estimate BMDs using the proposed method to have samples {BMD 1 , BMD 2 , . . . , BMD l } for BMD. The 100(1 − τ )% bootstrap lower confidence limit for BMD is given by the τ th quantile of the sample {BMD 1 , BMD 2 , . . . , BMD l }.
Let us denote ML, LR, ST, and BT for the methods of estimating BMDL using ML estimates, LR 3 Example: Experiment on Rats Exposed to 1-Bromopropane For illustrating the proposed method, we present an example of estimating BMD using a real data set on lung cancer incidence of rats exposed to 1-Bromopropane given in the NTP Technical Report TR-569 (Program et al., 2011). In this study, four groups of rats with each group contains 50 rats are exposed to four dose levels of 1-Bromopropane. After two years of studies, the observed lung cancer incidence of rats at four dose levels 0 ppm , 62.5 ppm, 125 ppm, and 250 ppm are recorded as 1/50, 9/50, 8/50, and 14/50 respectively as given in Table 1.
Wheeler and Bailer (2012) also analyzed the same data set noting that "the data, given in Table 1 Quantal-Linear for BMR=0.01, and 0.1 from Wheeler and Bailer (2012) in Table 2.
From  of BMD using MA diverge and smaller than those using other methods. We observe that the estimated values of BMDLs using FL are higher than those by the methods given in Wheeler and Bailer (2012).
So, the proposed methods may provide a better estimates of BMDL if the estimated confidence intervals have the expected coverage probabilities for small samples. So, we need to do simulation studies for verifying the coverage probabilities of the proposed confidence intervals for small samples.

Simulation Studies
In this section, we conduct simulation studies for testing the performance of the proposed methods of estimating BMD and BMDL considering all types of possible cases of generating data sets. Let us denote the proposed method of estimating BMD using the family of link functions as FL. The model averaging method is usually denoted as MA. For testing the performance of FL compare to MA, we provide a simulation study by estimating BMD using FL and MA considering different simulation scenarios with varying sample sizes. Simulation studies are also conducted for testing the performance of the proposed methods of estimating BMDL with respect to their coverage probabilities for small samples.

Comparison between FL and MA
We compare the proposed method FL with MA considering the following simulation set up with experimental design consists of four dose levels as d 1 = 0, d 2 = 0.25, d 3 = 0.5, and d 4 = 1.00 mimicking the design considered by West et al. (2012). Six scenarios for dose response relationships have been considered to represent all types of possibilities of having probability of adverse events at dose levels varying from shallow to steep curves. The scenarios with true parameter values and the probabilities of adverse events at dose levels are given in Table 3.  Table 4 are considered.
The expression for model averaging estimates of BMD, denoted as BMD M A is given by whereB k is the estimate of BMD using model k, and w k = exp(−0.5A k ) 8 k=1 exp(−0.5A k ) with A k is the Akaike Information Criteria (Akaike, 1973) given by A k = −2L k +2p k , whereL k is the maximized log-likelihood value and p k is the number of parameters in model k.
We generate l = 2000 data sets for each simulation set up. For some cases the simulated responses Log-logistic produce virtually flat dose repones curve (Wheeler and Bailer, 2009) which does not give any finite estimate of BMD. So, we mimic the methodology given in Wheeler and Bailer (2009) of screening the data sets using Kendall correlation test (Kendall, 1955). We regenerate responses until the responses exhibit the Kendall p-value less than or equal to 0.15.  Table 5.
From in Scenario 1 has very slowly increasing probability of adverse events at dose levels with R(d 1 ) = 0,

Comparison among Four BMDL Estimation Methods
Here, we conduct simulation studies to compare four methods of estimating BMDL using ML estimates  Table 3 to generate data sets. We also consider two values of BMR = 0.01, & 0.1 and three sample sizes n = 25, 50, & 100 for each scenario. For each simulation set up, we generate l = 1000 data sets which are also screened by Kendall correlation test (Kendall, 1955) as discussed in Section 4.1.
The simulated data sets are used to estimate 95% BMDL using the four methods ML, LR, ST, and BT. After estimating BMDL using a method, we find an approximate value of coverage probability given by N l l , where N l is the number of times the estimated values of BMDL are less than or equal to BMD out of l data sets generated. The coverage probabilities by four methods ML, LR, ST, and BT for each simulation set-up are reported in Table 6.
From Table 6, we see that the observed coverage probabilities by LR and ST are greater than 0.95 for all scenarios and BMR values with all sample sizes. The method BT fails to provide the expected coverage probabilities for Scenario 2 with n = 25, 50 and Scenarios 5 & 6 for all sample sizes when BMR=0.01. The observed coverage probabilities by ML also exceed the expected probability 0.95 for all the cases except for Scenario 5 with n = 50, 100, when BMR=0.1. We also studied the observed average length of one sided confidence interval (average of [BMD-BMDL]) by four methods. We observe that BT provides smallest values and ML & ST provide largest values for the average of (BMD-BMDL) for all the cases considered. Hence, we conclude that LR is best among all the methods of estimating BMDL with respect to coverage probabilities and lengths of the confidence intervals. The proposed method is illustrated by an example with a real data set observing that FL is consistent with the existing results in literature. By comparing FL with MA using simulation studies considering different simulation scenarios, we see that FL outperforms MA for most of the scenarios. Simulation studies are also conducted to compare the four methods of estimating BMDL and we see that LR is best among the four methods of estimating BMDL considering both the coverage probability as well as the length of the confidence intervals.

Conclusions
There are other methods exist in literature using Bayesian and non parametric approach for accounting model uncertainty in BMD estimation. The frequentist methods are usually easy to implement and require less time for computations than other non frequentist approach. We compared FL with MA as both the methods are based on frequentist approach to deal with the model uncertainty problems. In future, the proposed method may be compared with other non frequentist approach to estimate BMD to test the performance of FL.