Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC3278194

Formats

Article sections

- Abstract
- 1 Introduction
- 2 The proposed methodology
- 3 Simulation studies
- 4 Application to hexavalent chromium data
- 5 Concluding remarks
- Supplementary Material
- References

Authors

Related links

J Stat Plan Inference. Author manuscript; available in PMC 2013 May 1.

Published in final edited form as:

J Stat Plan Inference. 2012 May 1; 142(5): 1047–1062.

doi: 10.1016/j.jspi.2011.11.003PMCID: PMC3278194

NIHMSID: NIHMS341167

See other articles in PMC that cite the published article.

Toxicologists and pharmacologists often describe toxicity of a chemical using parameters of a nonlinear regression model. Thus estimation of parameters of a nonlinear regression model is an important problem. The estimates of the parameters and their uncertainty estimates depend upon the underlying error variance structure in the model. Typically, a priori the researcher would know if the error variances are homoscedastic (i.e., constant across dose) or if they are heteroscedastic (i.e., the variance is a function of dose). Motivated by this concern, in this article we introduce an estimation procedure based on preliminary test which selects an appropriate estimation procedure accounting for the underlying error variance structure. Since outliers and influential observations are common in toxicological data, the proposed methodology uses M-estimators. The asymptotic properties of the preliminary test estimator are investigated; in particular its asymptotic covariance matrix is derived. The performance of the proposed estimator is compared with several standard estimators using simulation studies. The proposed methodology is also illustrated using a data set obtained from the National Toxicology Program.

Often toxicologists are interested in investigating the dose-response relationship when animals are exposed to varying doses of a chemical. Usually a nonlinear regression model such as a Hill model is used to describe the relationship (Gaylor and Aylward, 2004; Sand et al., 2004; Crofton et al., 2007). There may be several problems when fitting nonlinear models. Among them, one important concern is the error variance structure. Depending upon various factors, including the bioassay, dose-spacing and the endpoint of interest etc., the variability in response may not be constant across dose groups (heteroscedasticity). The standard asymptotic confidence intervals and test procedures based on the ordinary least squares (OLS) methodology may not be robust to heteroscedasticity and consequently may produce inaccurate coverage probabilities and Type I error rates. On the other hand, the standard iterated weighted least squares (IWLS) based methodology may not be efficient when the variances are approximately equal across dose groups (homoscedasticity). However, in practice one generally does not know if the data are homoscedastic or heteroscedastic.

The problem of heteroscedasticity has been extensively discussed in the literature in a wide range of contexts involving linear and nonlinear models. For instance, Hoferkamp and Peddada (2002) considered the problem of heteroscedasticity in the context of groups of experiments, as in fertilizer trials, where the error variances are ordered. Cysneiros et al. (2007) derived a joint iterative process for estimating the location and dispersion parameters in heteroscedastic linear models with symmetrical errors. Guo and Koul (2008) developed asymptotic theory for long memory time series based on heteroscedastic linear models. Recently, the problem of heteroscedasticity has also been addressed in the context of semi-parametric partially linear models (Ma et al., 2006; You et al., 2007; Lu, 2009). A Bayesian method for testing for equality of regression parameters in a heteroscedastic linear model has also been considered in the literature (Moreno et al., 2005).

Several authors modeled error variance as a function of dose in dose-response models (cf. Davidian and Carroll, 1987). Wang and Zhou (2007) even developed a nonparametric test for checking the adequacy of a given variance function. Bellio et al. (2000) proposed the use of higher order likelihood based methods for inference in heteroscedastic nonlinear models with application to dose-response models in herbicide bioassays.

Although IWLS based methods perform well under heteroscedasticity, they may lose efficiency relative to other methods when the data are homoscedastic. To illustrate this, consider the simulated data presented in Fig. 1 which is based on homoscedastic errors. The data were generated using the Hill model (Hill, 1910),
$y={\theta}_{0}+{\theta}_{1}{x}^{{\theta}_{2}}/({\theta}_{3}^{{\theta}_{2}}+{x}^{{\theta}_{2}})+\epsilon $, which is usually used to study *in vivo* concentration response relationships, where *y* is the response at dose *x*, *θ*_{0} is the intercept parameter, *θ*_{1} is the difference between the maximum effect of a drug (E_{max}) and the intercept, *θ*_{2} is the slope parameter that reflects the steepness of the effect-concentration curve, and *θ*_{3} is the sensitivity parameter, the drug concentration producing 50% of E_{max} (ED_{50}).

Example of model predictions by OLSE (solid), IWLSE_{N} (dot), IWLSE_{M} (dashes) methods for homoscedastic data.

In Fig. 1, the point estimate (with standard error in parentheses) for the ED_{50} (whose true value is 120) based on the OLS estimator (OLSE) is 149.7 (31.90). The IWLS estimator using the sample variances, denoted as IWLSE* _{N}*, is 226.4 (65.37), and the IWLS estimator using a variance model, denoted as IWLSE

Because the performance of a method relies on whether the data are homoscedastic or heteroscedastic, it is important to develop an estimation procedure which is robust to whether the error variance is homoscedastic or heteroscedastic. To make the procedure robust to the structure of the error variance, a preliminary test estimation (PTE) based methodology is developed in this paper. PTE has been well studied in the literature in a variety of contexts (Judge and Bock, 1978). For instance, Sen (1986) studied the asymptotic distributional risks for the preliminary test version of a maximum likelihood estimator. Recently, Ahmed et al. (2007) investigated the asymptotic properties of a pretest semiparametric estimator under quadratic loss and examined its performance using asymptotic analysis of quadratic risk functions in a partially linear model. Hoque et al. (2009) studied the performance of the PTE of the slope parameter of a simple linear regression model under a linex loss function and derived the risk function and the moment generating function of the PTE.

Outliers and influential observations are common in toxicological data. To make the proposed procedure robust against outliers, we use the principle of M-estimation. Thus in this paper the PTE is either the ordinary M-estimator (OME) or the weighted M-estimator (WME) depending upon the outcome of a preliminary test for heteroscedasticity.

The PTE methodology based on OME and WME is proposed in Section 2. Results of a sample of simulation studies are provided in Section 3 and the proposed methodology is illustrated using a data set from the National Toxicology Program (NTP) in Section 4. Proofs of the main results are provided in Appendix while the theorems needed for proving these main results are provided in the online supplementary material.

Let *y _{i}* denote an

$${y}_{i}=f({x}_{i},\theta )+{\sigma}_{i}{\epsilon}_{i},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}i=1,\dots ,k,$$

denote the nonlinear regression model, where *f*(*x _{i}*,

It is assumed that *σ _{i}*

The definition of an M-estimator depends upon the Huber score function which is defined as follows. For a pre-specified positive constant *k*_{0}, the Huber score function *h*(*u*) is given by:

$$h(u)=\{\begin{array}{ll}u/\sqrt{2},\hfill & \text{if}u\phantom{\rule{0.16667em}{0ex}}{k}_{0}{\{{k}_{0}(u-{k}_{0}/2)\}1/2,}^{}\text{otherwise}.\hfill \hfill \hfill \end{array}$$

Throughout this paper we took *k*_{0} to be 1.5. Then the ordinary M-estimator for *θ* is obtained by solving the following minimization problem:

$${S}_{o}(\theta )=\sum _{i,j}{h}^{2}({y}_{ij}-f({x}_{i},\theta )).$$

If we take *h*(*u*) = *u* then we obtain the classical OLSE. Note that the OME does not account for heteroscedasticity in the data. The estimating equations for solving the above optimization problem are given by (Sanhueza and Sen, 2001):

$$\sum _{i,j}{\lambda}_{o}({x}_{i},{y}_{ij},{\stackrel{~}{\theta}}_{n})=0,$$

where

$${\lambda}_{o}({x}_{i},{y}_{ij},\theta )=\psi ({y}_{ij}-f({x}_{i},\theta )){f}_{\theta}({x}_{i},\theta ),$$

*f _{θ}*(

To deal with heteroscedasticity, one may define a weighted M-estimation procedure similar in spirit to the popular IWLSE methodology using a variance function. Thus the WME is obtained by solving the following minimization problem:

$$\left(\begin{array}{l}{\widehat{\theta}}_{n}\hfill \\ {\widehat{\tau}}_{n}\hfill \end{array}\right)=\mathit{Argmin}[\sum _{i,j}\left\{{h}^{2}\left(\frac{{y}_{ij}-f({x}_{i},\theta )}{\sigma ({z}_{i},\tau )}\right)+log\sigma ({z}_{i},\tau )\right\}:\theta {\mathfrak{R}}^{p},\tau {\mathfrak{R}}^{q}],$$

where log *σ*(*z _{i}*,

$$\sum _{i,j}\lambda ({x}_{i},{y}_{ij},{\widehat{\theta}}_{n},{\widehat{\tau}}_{n})=0,$$

(1)

where

$$\lambda ({x}_{i},{y}_{ij},\theta ,\tau )=\left(\begin{array}{c}k({z}_{i},\tau )\psi ({\epsilon}_{ij}){f}_{\theta}({x}_{i},\theta )\\ k({z}_{i},\tau )\phantom{\rule{0.16667em}{0ex}}\{\psi \phantom{\rule{0.16667em}{0ex}}({\epsilon}_{ij})\phantom{\rule{0.16667em}{0ex}}{\epsilon}_{ij}-1\}{\sigma}_{\tau}({z}_{i},\tau )\end{array}\right),$$

(2)

*σ _{τ}* (

We now derive the asymptotic normality of the WME. It is important to note that all the asymptotic results discussed in this paper are valid as long as *n*, the total sample size, goes to infinity. Two types of asymptotics often discussed in the literature are (a) the *n _{i}* is finite (e.g.,

We require the following sets of regularity conditions concerning (A) the score function *ψ*, (B) the function *f* and (C) the function *σ*.

- [A1]Let
*ε*= {*y*−*f*(*x*,*θ*)}/*σ*(*z*,*τ*). Then*γ*_{1}=*E*{*ψ*(*ε*)*ε*}(≠ 0);*Eψ*′(*ε*) =*γ*_{2}(≠ 0);*E*{*ψ*′(*ε*)*ε*^{2}} =*γ*_{3}(≠ 0); $E{\psi}^{2}(\epsilon )={\sigma}_{\psi 1}^{2}<\infty ;var\{\psi (\epsilon )\epsilon \}={\sigma}_{\psi 2}^{2}<\infty $. - [B1]
- lim
_{n}_{→∞}*n*^{−1}Γ_{1}(_{n}*θ*,*τ*) = Γ_{1}(*θ*,*τ*), where$${\mathrm{\Gamma}}_{1n}(\theta ,\tau )={\gamma}_{2}\sum _{i=1}^{n}{k}^{2}({z}_{i},\tau ){f}_{\theta}({x}_{i},\theta ){f}_{\theta}^{\text{T}}({x}_{i},\theta ).$$ - lim
_{n}_{→∞}*n*^{−1}Γ_{31}(_{n}*θ*,*τ*) = Γ_{31}(*θ*,*τ*), where$${\mathrm{\Gamma}}_{31n}(\theta ,\tau )={\sigma}_{\psi 1}^{2}\sum _{i=1}^{n}{k}^{2}({z}_{i},\tau ){f}_{\theta}({x}_{i},\theta ){f}_{\theta}^{\text{T}}({x}_{i},\theta ),$$and Γ_{31}(*θ*,*τ*) is a positive definite matrix. - ${max}_{i}\{{k}^{2}({z}_{i},\tau ){f}_{\theta}^{\text{T}}({x}_{i},\theta ){\mathrm{\Gamma}}_{31n}^{-1}(\theta ,\tau ){f}_{\theta}({x}_{i},\theta )\}\to 0$, as
*n*→ ∞

- [C]
- lim
_{n}_{→∞}*n*^{−1}Γ_{2}(_{n}*θ*,*τ*) = Γ_{2}(*θ*,*τ*), where$${\mathrm{\Gamma}}_{2n}(\theta ,\tau )=\sum _{i=1}^{n}\left\{\frac{2{\gamma}_{1}+{\gamma}_{3}-1}{{\sigma}^{2}({z}_{i},\tau )}{\sigma}_{\tau}({z}_{i},\tau ){\sigma}_{\tau}^{\text{T}}({z}_{i},\tau )+\frac{1-{\gamma}_{1}}{\sigma ({z}_{i},\tau )}{\mathrm{\sum}}_{\tau}({z}_{i},\tau )\right\},$$and Σ(_{τ}*z*,_{i}*τ*) = (^{2}/*ττ*^{T})*σ*(*z*,_{i}*τ*). - lim
_{n}_{→∞}*n*^{−1}Γ_{32}(_{n}*θ*,*τ*) = Γ_{32}(*θ*,*τ*), where$${\mathrm{\Gamma}}_{32n}(\theta ,\tau )={\sigma}_{\psi 2}^{2}\sum _{i=1}^{n}{k}^{2}({z}_{i},\tau ){\sigma}_{\tau}({z}_{i},\tau ){\sigma}_{\tau}^{\text{T}}({z}_{i},\tau ),$$and Γ_{32}(*θ*,*τ*) is a positive definite matrix. - ${max}_{i}\{{k}^{2}({z}_{i},\tau ){\sigma}_{\tau}^{\text{T}}({z}_{i},\tau ){\mathrm{\Gamma}}_{32n}^{-1}(\theta ,\tau ){\sigma}_{\tau}({z}_{i},\tau )\}\to 0$, as
*n*→ ∞

The asymptotic normality of the WME is established in Theorem 1 under the above regularity conditions.

Under the conditions [A1], [B1] and [C]; [S1] – [S9] in the supplementary material,

$${\widehat{\mathrm{\Gamma}}}^{-{\scriptstyle \frac{1}{2}}}\surd n\left(\begin{array}{c}{\widehat{\theta}}_{n}-\theta \\ {\widehat{\tau}}_{n}-\tau -{\nu}_{n}(\theta ,\tau )\end{array}\right)\to {N}_{p+q}(0,{I}_{p+q}),$$

(3)

where

$$\begin{array}{c}{\nu}_{n}(\theta ,\tau )={\left(\frac{1}{n}{\mathrm{\Gamma}}_{2n}(\theta ,\tau )\right)}^{-1}\frac{{\gamma}_{1}-1}{n}\sum _{i=1}^{n}k({z}_{i},\tau ){\sigma}_{\tau}({z}_{i},\tau ),\\ \widehat{\mathrm{\Gamma}}={\left(\frac{1}{n}{\mathrm{\Gamma}}_{5n}({\widehat{\theta}}_{n},{\widehat{\tau}}_{n})\right)}^{-1}\left(\frac{1}{n}{\mathrm{\Gamma}}_{3n}({\widehat{\theta}}_{n},{\widehat{\tau}}_{n})\right)\phantom{\rule{0.16667em}{0ex}}{\left(\frac{1}{n}{\mathrm{\Gamma}}_{5n}({\widehat{\theta}}_{n},{\widehat{\tau}}_{n})\right)}^{-1},\\ {\mathrm{\Gamma}}_{3n}(\theta ,\tau )=\left(\begin{array}{cc}{\mathrm{\Gamma}}_{31n}(\theta ,\tau )& 0\\ 0& {\mathrm{\Gamma}}_{32n}(\theta ,\tau )\end{array}\right),\end{array}$$

and

$${\mathrm{\Gamma}}_{5n}(\theta ,\tau )=\left(\begin{array}{cc}{\mathrm{\Gamma}}_{1n}(\theta ,\tau )& 0\\ 0& {\mathrm{\Gamma}}_{2n}(\theta ,\tau )\end{array}\right).$$

Note that the above theorem also provides the asymptotic covariance matrix of the WME.

We now describe the PTE procedure in the context of dose-response studies. Based on our experience with dose-response studies in toxicology, it is reasonable to assume that log-variance is a linear function of dose. Thus we assume log *σ _{i}* =

Then, the PTE is defined as

$${\widehat{\theta}}_{n}^{\text{PT}}=\{\begin{array}{ll}{\stackrel{~}{\theta}}_{n}\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}{T}_{n}\le {t}_{\alpha ,n-2}\hfill \\ {\stackrel{~}{\theta}}_{n}\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}{T}_{n}>{t}_{\alpha ,n-2},\hfill \end{array}$$

(4)

where *t _{α}*

In order to derive asymptotic results regarding PTE, we require the following sets of regularity conditions.

- [A2]Let
*ε*= {*y*−*f*(*x*,*θ*)}/*σ*(*z*,*τ*). Then,*Eψ*′(*σ*(*z*,*τ*)*ε*) =*γ*_{4}(≠ 0), $E{\psi}^{2}(\sigma (z,\tau )\epsilon )={\sigma}_{\psi 3}^{2}{w}_{1}(x)<\infty $ and $E\{\psi (\epsilon )\psi (\sigma (z,\tau )\epsilon )\}={\sigma}_{\psi 4}^{2}{w}_{2}(x)<\infty $. - [B2]
- lim
_{n}_{→∞}*n*^{−1}Γ_{4}(_{n}*θ*) = Γ_{4}(*θ*), where$${\mathrm{\Gamma}}_{4n}(\theta )={\gamma}_{4}\sum _{i=1}^{k}{n}_{i}{f}_{\theta}({x}_{i},\theta ){f}_{\theta}^{\text{T}}({x}_{i},\theta ).$$ - lim
_{n}_{→∞}*n*^{−1}Γ_{33}(_{n}*θ*) = Γ_{33}(*θ*), where$${\mathrm{\Gamma}}_{33n}(\theta )={\sigma}_{\psi 3}^{2}\sum _{i=1}^{k}{n}_{i}{w}_{1}({x}_{i}){f}_{\theta}({x}_{i},\theta ){f}_{\theta}^{\text{T}}({x}_{i},\theta ),$$and Γ_{33}(*θ*) is a positive definite matrix. - lim
_{n}_{→∞}*n*^{−1}Γ_{34}(_{n}*θ*,*τ*) = Γ_{34}(*θ*,*τ*), where$${\mathrm{\Gamma}}_{34n}(\theta ,\tau )={\sigma}_{\psi 4}^{2}\sum _{i=1}^{k}\frac{{n}_{i}{w}_{2}({x}_{i})}{{\sigma}_{i}}{f}_{\theta}({x}_{i},\theta ){f}_{\theta}^{\text{T}}({x}_{i},\theta ),$$ - lim
_{n}_{→∞}*n*^{−1}*G*_{2}(_{n}*θ*,*τ*) =*G*_{2}(*θ*,*τ*), where$${G}_{2n}(\theta ,\tau )=\left(\begin{array}{ccc}{\mathrm{\Gamma}}_{31n}(\theta ,\tau )& {\mathrm{\Gamma}}_{34n}(\theta ,\tau )& 0\\ {\mathrm{\Gamma}}_{34n}(\theta ,\tau )& {\mathrm{\Gamma}}_{33n}(\theta )& 0\\ 0& 0& 2{n}^{2}\sum _{i=1}^{k}{n}_{i}{w}_{i12}^{2}\end{array}\right),$$(5)*w*_{i}_{12}is the second element of*w*_{i}_{1}= (*Z*^{T}*Z*)^{−1}*z*_{i}_{1}, and*G*_{2}(*θ*,*τ*) is a positive definite matrix. - max
_{i}ch_{1}{*G*_{1}(*x*,_{i}*θ*,*τ*)*G*_{2}(_{n}*θ*,*τ*))^{−1}} → 0, as*n*→ ∞, where$${G}_{1}({x}_{i},\theta ,\tau )=\left(\begin{array}{ccc}{\sigma}_{\psi 1}^{2}{\sigma}_{i}^{-2}{H}_{i}& {\sigma}_{\psi 4}^{2}{w}_{2}({x}_{i}){\sigma}_{i}^{-1}{H}_{i}& 0\\ {\sigma}_{\psi 4}^{2}{w}_{2}({x}_{i}){\sigma}_{i}^{-1}{H}_{i}& {\sigma}_{\psi 3}^{2}{w}_{1}({x}_{i}){H}_{i}& 0\\ 0& 0& 2{n}^{2}{w}_{i12}^{2}\end{array}\right),$$and ${H}_{i}={f}_{\theta}({x}_{i},\theta ){f}_{\theta}^{\text{T}}({x}_{i},\theta )$.

We begin by proving the asymptotic joint normality of OME and WME which is then used for deriving the asymptotic covariance matrix of the PTE. The proof, provided in the appendix, follows arguments similar to that of Theorem 1.

Let *β* = (*θ*^{T}, *θ*^{T}, *τ*_{1})^{T} and
${\widehat{\beta}}_{n}={\left({\widehat{\theta}}_{n}^{\text{T}},{\stackrel{~}{\theta}}_{n}^{\text{T}},{\widehat{\tau}}_{1n}\right)}^{\text{T}}$. Then, under the conditions [A1], [A2], [B1], [B2] and [C]; [S1] – [S9] in the supplementary material,

$$\surd n\left({\widehat{\beta}}_{n}-\beta \right)\to {N}_{2p+1}(0,G(\theta ,\tau ))\phantom{\rule{0.38889em}{0ex}}as\phantom{\rule{0.38889em}{0ex}}n\to \infty ,$$

(6)

where

$$\begin{array}{l}G(\theta ,\tau )={G}_{3}^{-1}(\theta ,\tau ){G}_{2}(\theta ,\tau ){G}_{3}^{-1}(\theta ,\tau ),\\ {G}_{3}(\theta ,\tau )=\left(\begin{array}{ccc}{\mathrm{\Gamma}}_{1}(\theta ,\tau )& 0& 0\\ 0& {\mathrm{\Gamma}}_{4}(\theta )& 0\\ 0& 0& 2\end{array}\right).\end{array}$$

From the above theorem we deduce the asymptotic covariance matrix of PTE in the following theorem.

Under the conditions [A1], [A2], [B1], [B2] and [C]; [S1] – [S9] in the supplementary material,

$$\begin{array}{l}E\left[n\left({\widehat{\theta}}_{n}^{\text{PT}}-\theta \right){\left({\widehat{\theta}}_{n}^{\text{PT}}-\theta \right)}^{\text{T}}\right]\\ ={F}_{t}\left({t}_{\alpha ,n-2}-\frac{{\tau}_{1}}{\surd var({\widehat{\tau}}_{1n})}\right)E\left[n\left({\stackrel{~}{\theta}}_{n}-\theta \right){\left({\stackrel{~}{\theta}}_{n}-\theta \right)}^{\text{T}}\right]\\ +\left\{1-{F}_{t}\left({t}_{\alpha ,n-2}-\frac{{\tau}_{1}}{\surd var({\widehat{\tau}}_{1n})}\right)\right\}E\left[n\left({\widehat{\theta}}_{n}-\theta \right){\left({\widehat{\theta}}_{n}-\theta \right)}^{\text{T}}\right]\\ ={F}_{t}\left({t}_{\alpha ,n-2}-\frac{{\tau}_{1}}{\surd var({\widehat{\tau}}_{1n})}\right){\left(\frac{1}{n}{\mathrm{\Gamma}}_{4n}(\theta )\right)}^{1}\left(\frac{1}{n}{\mathrm{\Gamma}}_{33n}(\theta )\right){\left(\frac{1}{n}{\mathrm{\Gamma}}_{4n}(\theta )\right)}^{-1}\\ +\left\{1-{F}_{t}\left({t}_{\alpha ,n-2}-\frac{{\tau}_{1}}{\surd var({\widehat{\tau}}_{1n})}\right)\right\}{\left(\frac{1}{n}{\mathrm{\Gamma}}_{1n}(\theta )\right)}^{-1}\left(\frac{1}{n}{\mathrm{\Gamma}}_{31n}(\theta ,\tau )\right){\left(\frac{1}{n}{\mathrm{\Gamma}}_{1n}(\theta )\right)}^{-1},\end{array}$$

(7)

where *F _{t}* is the cdf of the t-distribution with

As we see from the above theorem, the asymptotic covariance matrix of PTE is a weighted average of those of OME and WME. The weights are directly related to *τ*_{1}. In particular, the weight corresponding to the covariance of OME is monotonically decreasing in *τ*_{1}. It equals 1 − *α* when *τ*_{1} = 0.

We simulated data using the following Hill model under three different error variance structures. In each case the errors are normally distributed with mean 0.

$${y}_{ij}=f({x}_{i},\theta )+{\epsilon}_{ij}={\theta}_{0}+\frac{{\theta}_{1}{x}_{i}^{{\theta}_{2}}}{{x}_{3}^{{\theta}_{2}}+{x}_{i}^{{\theta}_{2}}}+{\epsilon}_{ij},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}i=1,\dots ,8,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}j=1,\dots ,5.$$

(8)

The values of *x _{i}* were set to be 0, 1, 3, 10, 30, 100, 400, 600, and (

There are two parts to the simulation study. Firstly, for illustration purposes, we generated one data set of each type (i.e., Data 1, Data 2 and Data 3) and the parameters were estimated using OLSE, IWLSE* _{N}*, IWLSE

Secondly, using 1,000 simulation runs, we compared the performance of the estimators in terms of three standard criteria: (i) mean squared error (MSE) of individual parameters as well as total MSE, (ii) the coverage probabilities of the 95% confidence intervals (CIs) of individual parameters as well as the simultaneous confidence ellipsoid defined below, and (iii) the length of the 95% CI of individual parameters as well as the volume of 95% confidence ellipsoid for each estimator. The 100(1 − *α*)% confidence ellipsoid centered at an estimator is defined as

$${(\widehat{\theta}-\theta )}^{\text{T}}{[\widehat{var}(\widehat{\theta})]}^{-1}(\widehat{\theta}-\theta )\le {pF}_{p,n-p}(\alpha ),$$

where *p* is the number of parameters and
$\widehat{var}(\widehat{\theta})$ is the appropriate variance estimator. Note that for simplicity we are using the critical values based on WME for the critical values for PTE because the exact critical values for PTE are not available. This is because the asymptotic distribution of PTE is not normal but a mixture of two normals.

The results of various estimates (and their standard errors) for the three simulated data sets without outliers are summarized in Table 1. As described in the introduction, for homoscedastic data (Data 1), although the fitted curves using IWLSE* _{N}*, IWLSE

Model predictions by OLSE (solid), IWLSE_{N} (dot), IWLSE_{M} (dashes), OME (long dashes), WME (dot-dash) methods when there are no outliers for: (a) homoscedastic data (Data 1), (b) heteroscedastic data (Data 2), and (c) mismodeled heteroscedastic data (Data **...**

Estimate and Standard Error for parameters of the models for Data 1, 2 and 3 without outliers using OLSE, IWLSE_{N}, IWLSE_{M}, OME, WME and PTE methods.

Note that if the data are homoscedastic (Data 1), then the “correct” choice of estimator is OME (OLSE when there are no outliers), whereas for heteroscedastic data (Data 2), the “correct” choice is WME (IWLSE* _{M}* when there are no outliers). However, in a practical setting, for a given data set one does not know

The results of various estimates (and their standard errors) for the three data sets with outliers are summarized in Table 2. Because there were outliers in the data, *θ*_{1} (E_{max}) and *θ*_{3} (ED_{50}) were severely underestimated using the least square estimators (true values of *θ*_{1} and *θ*_{3} are 4 and 120, respectively), while OME and WME were closer to the true values. Figure 3 also reflects this result. It is noted that because of the outliers the preliminary test rejected the null hypothesis of homoscedasticity and PTE selected WME, which was closer to the true value than OME, even though the data were homogeneous.

Model predictions by OLSE (solid), IWLSE_{N} (dot), IWLSE_{M} (dashes), OME (long dashes), WME (dot-dash) methods when there are outliers for: (a) homoscedastic data (Data 1), (b) heteroscedastic data (Data 2), and (c) mismodeled heteroscedastic data (Data **...**

Estimate and Standard Error for parameters of the models for Data 1, 2 and 3 with outliers using OLSE, IWLSE_{N}, IWLSE_{M}, OME, WME and PTE methods.

Table 3 shows the results of the 1,000 simulations using data without outliers. When data were generated from homoscedastic model (Data 1), as expected, the estimated mean squared errors of OLSE and OME were smaller than those of the other estimators (except *θ*_{0}, the intercept) and the estimated MSEs of PTE were slightly larger than those of OME and much smaller than those of IWLSE* _{N}* and IWLSE

Simulation results based on 1,000 replications generated from Data 1, 2 and 3 without outliers for OLSE, IWLSE_{N}, IWLSE_{M}, OME, WME and PTE.

In the case of heteroscedastic data (Table 3) we observe that OLSE and OME could potentially perform extremely poorly. This is because when there is a large variation in the higher dose groups the observed data may fail to “plateau” at higher dose groups. Consequently, estimators such as the OLSE and OME would tend to overestimate *θ*_{1} (E_{max}) and *θ*_{3} (ED_{50}). For this reason, for some random samples, the estimates of ED_{50} became arbitrarily large. As a consequence the estimated MSEs of OLSE and OME became extremely large! However, as one would expect, estimators such as IWLSE and WME performed well for such data with smaller MSE, and approximately correct coverage probability. The PTE performed as well as IWLSE and WME in terms of MSE and coverage probability. It competed very well in terms of efficiency relative to the weighted estimators (in terms of total MSE). For example, in the case of Data 2, it was 23% more efficient than IWLSE* _{N}*, about 4% less efficient than IWLSE

Table 4 shows the results of the 1,000 simulations using data with outliers. As expected, the least squares based methods performed very poorly in terms of MSE and coverage probability in the presence of outliers, while the M-estimation based methods performed much better. In some cases the coverage probability of CIs centered at the least squares estimators were substantially smaller than the nominal level. For example, in the case of Data 2, the coverage probability of CIs centered at IWLSE* _{M}* for parameter

Simulation results based on 1,000 replications generated from Data 1, 2 and 3 with outliers for OSLE, IWLSE_{N}, IWLSE_{M}, OME, WME and PTE.

In all the cases, the gains in efficiency (in terms of total MSE) of WME and PTE relative to IWLSE* _{N}* and IWLSE

Our simulation studies made a strong case for the use of the proposed methodology. The gains in terms of MSE were generally substantial.

As commonly understood, the performance of an estimator may be affected much by dose placement. To illustrate this point, we generated a data set from the Hill model with homoscedastic error. The true values of the parameters are (*θ*_{0}, *θ*_{1}, *θ*_{2}, *θ*_{3}) = (2, 4, 2, 30) and the values of the dose are *x* = 0, 1, 3, 10, 30, 100, 300. Because all estimators considered in this paper are affected, the OLSE was used to fit the curve for simplicity. The fitted curve and the estimation result are presented in Fig. 4(a) and Table 5, respectively. Even though the fitted curve is visually reasonable based on the data, the estimate of the parameter *θ*_{2} (slope), its standard error and the standard error of the estimate of *θ*_{3} (ED_{50}) are extremely large. We then chose additional dose arbitrarily, that is, *x* = 70, to generate observations at the dose and estimate the parameters and their standard errors based on the new data set. The fitted curve and the estimation result are presented in Fig. 4(b) and Table 5, respectively. Now the estimate of *θ*_{2}, its standard error and the standard error of the estimate of *θ*_{3} are all reasonably small. This illustration suggests the same argument presented in Lim et al. (2011) that “dose-spacing plays a major role when estimating parameters of nonlinear models, especially the ED_{50} and the slope parameters of a Hill model”.

Data generated from the Hill model and the fitted curve (a) without *x* = 70 and (b) with *x* = 70 using OLSE.

In this section the proposed PTE methodology is applied to a data set from a toxicological study that was designed to examine the relationship between concentrations of Hexavalent Chromium (CrVI), as sodium dichromate dihydrate, in drinking water and accumulation of total chromium in tissue for three species (rats, mice, and guinea pigs) (NTP, 2007).

As commonly done in toxicology, we use the Hill model (8) to describe the dose-response relationship. In our model, *x* denotes the dose (in mg/L), ranging from 0 to 300, and *y* denotes the total chromium concentration (in mg/L) in the tissues. The proposed methodology is illustrated using rat kidney and blood data sets where chromium concentration (*y*) is modeled. There were 7 dose levels and 4 observations at each dose level except *x* = 0 (3 observation). Thus, total sample size is 27. Based on each of the data sets the parameters (and their standard errors) are estimated using OLSE, IWLSE* _{N}*, IWLSE

Chromium concentration in (a) rat kidney and (b) rat blood using OLSE (solid), IWLSE_{N} (dot-dash), IWLSE_{M} (long dashes), OME (dot), WME (dashes) methods.

Estimate and Standard Error for parameters of the models for chromium rat kidney and blood data using OLSE, IWLSE_{N}, IWLSE_{M}, OME, WME and PTE methods.

Visually the scatter plot in Fig. 5(a) seems to suggest that there is some amount of heteroscedasticity in the data. The sample variances of kidney chromium for the seven dose groups were 0.0012, 0.0046, 0.0006, 0.054, 0.720, 0.213, and 6.507, respectively. Thus they ranged from about 0.0006 to 6.507, which indicates potential heteroscedasticity in the data. When the log-linear model for the absolute residual based on the OME was fitted against doses, the estimated value of *τ*_{1} was 0.005 with a standard error of 0.002. The slope was significant at the 5% level (*p* = 0. 009). Thus the data appears to be heteroscedastic.

The data in Fig. 5(a) appear to be heteroscedastic and hence it is not surprising that the point estimates and their standard errors are quite different between ordinary estimates (OLSE and OME) and weighted estimates (IWLSE* _{N}*, IWLSE

On the other hand, Fig. 5(b) shows that the data might be homoscedastic. The sample variances of blood chromium were 0.0006, 0.0002, 0.0003, 0.0001, 0.0007, 0.0016, and 0.0013, respectively. Thus they ranged from about 0.0001 to 0.0016. The result of the preliminary test using the absolute residuals based on the OME revealed that the slope (*τ*_{1}) was not significantly greater than zero at the 5% level (*p* = 0. 531). Thus the data appear homoscedastic.

Visually the fitted curves are almost identical except IWLSE* _{N}*. However, Table 6 shows that point estimates and their standard errors are not similar, although OLSE and OME are exactly the same. Again, point estimates of

In this paper, PTE based methodology has been developed for analyzing nonlinear models that are possibly subject to heteroscedastic variance structure. The methodology proposed here allows researchers to use estimation procedures that are robust to the error variance structure in nonlinear models. We demonstrated its utility using simulation studies and a real data set obtained by the NTP on chromium VI.

The proposed methodology depends on the model used for describing the heteroscedasticity. In our experience, a log-linear model for variance is plausible for data observed in these toxicological experiments because of the underlying sigmoidal shaped dose-response curve and variance being monotonic with mean response. Also, the log-linear model offers a simple interpretation of variability with fewer parameters. If, however, an experimenter is interested in using a different parametric model to describe the variance, then the proposed methodology can be modified easily. For example, Lim et al. (2011a) have studied WME using a different variance model, where the error standard deviation was assumed to be a nonlinear function of three unknown parameters.

According to their document “Benchmark Dose Technical Guidance Document” (USEPA, 2000), for independent continuous outcome variables, the EPA determines BMD using nonlinear least squares methodology. They typically estimate various parameters of the model, including BMD, using the ordinary least squares estimator (OLSE) for homoscedastic data and the weighted least squares estimator (WLSE) for heteroscedastic data. According to the above document, they determine whether the data are homoscedastic or heteroscedastic by visual judgment using a scatter plot of the data. After making decisions regarding the error variance on the basis of the observed data, the EPA either uses the OLSE or the WLSE. However, since the choice of WLSE and OLSE depends upon the observed data, the standard errors of the estimators of various parameters (including BMD) should account for the uncertainty in decision made regarding the error variance. This is not done by the EPA’s methodology, which is the point of our paper where we account for the uncertainty induced by the preliminary evaluation of data regarding the error variance. Hence we believe that the methodology developed in this paper can be extended to estimate the BMD.

This research was supported, in part, by the Intramural Research Program of the NIH, National Institute of Environmental Health Sciences [Z01 ES101744-04]. We thank Drs Gregg Dinse and Paramita Saha as well as the editors and the referees for many important comments which helped improve the presentation of the manuscript.

The proof relies on the asymptotic linearity of the WME (Lim et al., 2011b) and the existence of a solution to (1) which yields a √*n*-consistent estimator of (*θ*^{T}, *τ*^{T})^{T} (Theorem 4 in the supplementary material).

From Theorem 5 in the supplementary material we have that:

$$\surd n\left(\begin{array}{c}{\widehat{\theta}}_{n}-\theta \\ {\widehat{\tau}}_{n}-\tau \end{array}\right)={\left(\frac{1}{n}{\mathrm{\Gamma}}_{5n}(\theta ,\tau )\right)}^{-1}\frac{1}{\surd n}\sum _{i=1}^{n}\lambda ({x}_{i},{y}_{i},\theta ,\tau )+{o}_{p}(1).$$

Then from Theorem 6 in the supplementary material and the Slutsky Theorem we have the expression in (3).

From the estimating equation (1) for the WME, we let

$${\lambda}_{w}({x}_{i},{y}_{ij},\theta )=\frac{1}{{\sigma}_{i}}\psi \left(\frac{{y}_{ij}-f({x}_{i},\theta )}{{\sigma}_{i}}\right){f}_{\theta}({x}_{i},\theta ).$$

Then, from (13) in Theorem 5 in the supplementary material, the following asymptotic representation of * _{n}* is obtained:

$$\surd n({\widehat{\theta}}_{n}-\theta )={\left(\frac{1}{n}{\mathrm{\Gamma}}_{1n}(\theta ,\tau )\right)}^{-1}\frac{1}{\surd n}\sum _{i,j}{\lambda}_{w}({x}_{i},{y}_{ij},\theta )+{o}_{p}(1).$$

(9)

Now, similarly as in Theorem 4 and 5 in the supplementary material, the uniform asymptotic linearity on the OME can be shown as:

$$\underset{t\le C\Vert \frac{1}{\surd n}\sum _{i,j}\{{\lambda}_{o}({x}_{i},{y}_{ij},\theta +{n}^{-{\scriptstyle \frac{1}{2}}}t)-{\lambda}_{o}({x}_{i},{y}_{ij},\theta )\}+\frac{1}{n}{\mathrm{\Gamma}}_{4n}(\theta )t\Vert ={o}_{p}(1),}{sup}$$

and hence the following asymptotic representation of * _{n}* is obtained:

$$\surd n({\widehat{\theta}}_{n}-\theta )={\left(\frac{1}{n}{\mathrm{\Gamma}}_{4n}(\theta )\right)}^{-1}\frac{1}{\surd n}\sum _{i,j}{\lambda}_{o}({x}_{i},{y}_{ij},\theta ,)+{o}_{p}(1).$$

(10)

Then, from (9), (10) and (16) in the supplementary material,

$$\surd n\left({\widehat{\beta}}_{n}-\beta \right)={G}_{3n}^{-1}(\theta ,\tau )\frac{1}{\surd n}\sum _{i,j}{\lambda}^{},$$

where

$${G}_{3n}(\theta ,\tau )=\left(\begin{array}{ccc}{n}^{-1}{\mathrm{\Gamma}}_{1n}(\theta ,\tau )& 0& 0\\ 0& {n}^{-1}{\mathrm{\Gamma}}_{4n}(\theta )& 0\\ 0& 0& 2\end{array}\right).$$

Then from Theorem 7 in the supplementary material and the Slutsky Theorem the expression in (5) is shown.

From (4), for arbitrary *x* * ^{p}*,

$$\begin{array}{l}P\left\{\surd n\left({\widehat{\theta}}_{n}^{\text{PT}}-\theta \right)\le x\right\}=P\left\{\surd n\left({\stackrel{~}{\theta}}_{n}-\theta \right)\le x,{T}_{n}\le {t}_{\alpha ,n-2}\right\}+P\left\{\surd n\left({\widehat{\theta}}_{n}-\theta \right)\le x,{T}_{n}>{t}_{\alpha ,n-2}\right\}\\ =P\left\{\surd n\left({\stackrel{~}{\theta}}_{n}-\theta \right)\le x\right\}P\{{T}_{n}\le {t}_{\alpha ,n-2}\}+P\left\{\surd n\left({\widehat{\theta}}_{n}-\theta \right)\le x\right\}P\{{T}_{n}>{t}_{\alpha ,n-2}\}\\ ={F}_{t}\left({t}_{\alpha ,n-2}-\frac{{\tau}_{1}}{\surd var({\widehat{\tau}}_{1n})}\right)P\left\{\surd n\left({\stackrel{~}{\theta}}_{n}-\theta \right)\le x\right\}+\left\{1-{F}_{t}\left({t}_{\alpha ,n-2}-\frac{{\tau}_{1}}{\surd var({\widehat{\tau}}_{1n})}\right)\right\}P\left\{\surd n\left({\widehat{\theta}}_{n}-\theta \right)\le x\right\}.\end{array}$$

The second equality above holds because of the asymptotic independence of * _{n}* and

$$\begin{array}{l}E\left[n\left({\widehat{\theta}}_{n}^{\text{PT}}-\theta \right){\left({\widehat{\theta}}_{n}^{\text{PT}}-\widehat{\theta}\right)}^{\text{T}}\right]\\ ={F}_{t}\left({t}_{\alpha ,n-2}-\frac{{\tau}_{1}}{\surd var({\widehat{\tau}}_{1n})}\right)E\left[n\left({\stackrel{~}{\theta}}_{n}-\theta \right){\left({\stackrel{~}{\theta}}_{n}-\theta \right)}^{\text{T}}\right]\\ +\left\{1-{F}_{t}\left({t}_{\alpha ,n-2}-\frac{{\tau}_{1}}{\surd var({\widehat{\tau}}_{1n})}\right)\right\}E\left[n\left({\widehat{\theta}}_{n}-\theta \right){\left({\widehat{\theta}}_{n}-\theta \right)}^{\text{T}}\right],\end{array}$$

and since from (6),

$$E\left[n\left({\stackrel{~}{\theta}}_{n}-\theta \right){\left({\stackrel{~}{\theta}}_{n}-\theta \right)}^{\text{T}}\right]={\left(\frac{1}{n}{\mathrm{\Gamma}}_{4n}(\theta )\right)}^{-1}\left(\frac{1}{n}{\mathrm{\Gamma}}_{33n}(\theta )\right){\left(\frac{1}{n}{\mathrm{\Gamma}}_{4n}(\theta )\right)}^{-1}$$

and

$$E\left[n\left({\widehat{\theta}}_{n}-\theta \right){\left({\widehat{\theta}}_{n}-\theta \right)}^{\text{T}}\right]={\left(\frac{1}{n}{\mathrm{\Gamma}}_{1n}(\theta ,\tau )\right)}^{-1}\left(\frac{1}{n}{\mathrm{\Gamma}}_{31n}(\theta ,\tau )\right){\left(\frac{1}{n}{\mathrm{\Gamma}}_{1n}(\theta ,\tau )\right)}^{-1},$$

the expression in (7) is proved.

Supplementary material associated with this article can be found in the online version.

**Publisher's Disclaimer: **This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

- Ahmed SE, Doksum KA, Hossain S, You J. Shrinkage, pretest and absolute penalty estimators in partially linear models. Aust N Z J Stat. 2007;49:435–454.
- Bellio R, Jensen JE, Seiden P. Applications of likelihood asymptotics for nonlinear regression in herbicide bioassays. Biometrics. 2000;56:1204–1212. [PubMed]
- Carroll RJ, Ruppert D. Transformation and Weighting in Regression. Chapman and Hall; New York: 1988.
- Crofton KM, Paul KB, DeVito MJ, Hedge JM. Short-term in vivo exposure to the water contaminant triclosan: evidence for disruption of thyroxine. Environ Toxicol Pharmacol. 2007;24:194–197. [PubMed]
- Cysneiros FJA, Paula GA, Galea M. Heteroscedastic symmetrical linear models. Statist Prob Lett. 2007;77:1084–1090.
- Davidian M, Carroll RJ. Variance function estimation. J Amer Statist Assoc. 1987;82:1079–1091.
- Gaylor DW, Aylward LL. An evaluation of benchmark dose methodology for non-cancer continuous-data health effects in animals due to exposures to dioxin (TCDD) Regul Toxicol Pharm. 2004;40:9–17. [PubMed]
- Guo H, Koul HL. Asymptotic inference in some heteroscedastic regression models with long memory design and errors. Ann Statist. 2008;36:458–487.
- Hill AV. The possible effects of the aggregation of the molecules of haemoglobin on its dissociation curves. J Physiol. 1910;40(Suppl):iv–vii.
- Hoferkamp C, Peddada SD. Parameter estimation in linear models with heteroscedastic variances subject to order restrictions. J Multivariate Anal. 2002;82:65–87.
- Hoque Z, Khan S, Wesolowski J. Performance of preliminary test estimator under linex loss function. Commun Stat–Theor M. 2009;38:252–261.
- Judge GG, Bock ME. The Statistical Implications of Pre-test and Stein-rule Estimators in Econometrics. North-Holland, Amsterdam: 1978.
- Lim C, Sen PK, Peddada SD. Statistical inference in nonlinear regression under heteroscedasticity. Sankhya B. 2011a;72:202–218.
- Lim C, Sen PK, Peddada SD. Asymptotic linearity of an M-estimator in heteroscedastic nonlinear regression models. 2011b submitted.
- Lu X. Empirical likelihood for heteroscedastic partially linear models. J Multivariate Anal. 2009;100:387–396.
- Ma Y, Chiou JM, Wang N. Efficient semiparametric estimator for heteroscedastic partially linear models. Biometrika. 2006;93:75–84.
- Moreno E, Torres F, Casella G. Testing equality of regression coefficients in heteroscedastic normal regression models. J Statist Plann Inference. 2005;131:117–134.
- National Toxicology Program. Toxicity Report Series 72. U.S. Department of Health and Human Services, Public Health Service, National Institutes of Health, RTP; North Carolina, U.S.A: 2007. NTP Toxicity Studies of Sodium Dichromate Dihydrate (CAS No. 7789-12-0) Administered in Drinking Water to Male and Female F344/N Rats and B6C3F1 Mice and Male BALB/c and am3-C57BL/6 Mice; pp. 1–G4. [PubMed]
- Sand S, von Rosen D, Eriksson P, Fredriksson A, Viberg H, Victorin K, Filpsson AF. Dose-response modeling and benchmark calculations from spontaneous behavior data on mice neonatally exposed to 2,2′,4,4′,5-pentabromodiphenyl ether. Toxicol Sci. 2004;81:491–501. [PubMed]
- Sanhueza AI, Sen PK. M-methods in generalized nonlinear models. In: Charalambides ChA, Koutras MV, Balakrishnan N., editors. Probability and Statistical Models with Applications. Chapman & Hall/CRC; New York: 2001. pp. 359–375.
- Sen PK. On the asymptotic distributional risks of shrinkage and preliminary test versions of maximum likelihood estimators. Sankhya Ser A. 1986;48:354–371.
- Sen PK, Singer JM, Pedroso de Lima AC. From Finite Samples to Asymptotic Methods in Statistics. Cambridge Univ. Press; New York: 2009.
- Peddada SD, Smith T. Consistency of a class of variance estimators in linear models under heteroscedasticity. Sankhya Ser B. 1997;59:1–10.
- US Environmetal Protection Agency. Benchmark Dose Technical Guidance Document. 2000. EPA/630/R-00/001.
- Wang L, Zhou XH. Assessing the adequacy of variance function in heteroscedastic regression models. Biometrics. 2007;63:1218–1225. [PubMed]
- Wu CFJ. Jackknife, bootstrap and other resampling methods in regression analysis. Ann Stat. 1986;14:1261–1295.
- You J, Chen G, Zhou Y. Statistical inference of partially linear regression models with heteroscedastic errors. J Multivariate Anal. 2007;98:1539–1557.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |