https://orcid.org/0000-0002-4133-6884
Abstract
Benchmark dose analysis aims to estimate the level of exposure to a toxin associated with a clinically significant adverse outcome and quantifies uncertainty using the lower limit of a confidence interval for this level. We develop a novel framework for benchmark dose analysis based on monotone additive dose-response models. We first introduce a flexible approach for fitting monotone additive models via penalized B-splines and Laplace-approximate marginal likelihood. A reflective Newton method is then developed that employs de Boor’s algorithm for computing splines and their derivatives for efficient estimation of the benchmark dose. Finally, we develop a novel approach for calculating benchmark dose lower limits based on an approximate pivot for the nonlinear equation solved by the estimated benchmark dose. The favorable properties of this approach compared to the Delta method and a parameteric bootstrap are discussed. We apply the new methods to make inferences about the level of prenatal alcohol exposure associated with clinically significant cognitive defects in children using data from six NIH-funded longitudinal cohort studies. Software to reproduce the results in this paper is available online and makes use of the novel semibmd R package, which implements the methods in this paper.
1 INTRODUCTION
1.1 Benchmark dose analysis
Benchmark dose analysis is a statistical technique used by environmental toxicologists to quantify the level of exposure to a harmful substance that is associated with an adverse response. The US Environmental Protection Agency (EPA, 2012) and the European Food Safety Authority (EFSA et al., 2022) use the lower limit of a |$95\%$| confidence interval of the benchmark dose—called the benchmark dose, lower (BMDL)—to set limits on acceptable levels of exposure to a wide variety of toxins. In this paper, we present an innovative general framework for benchmark dose analysis and inference based on monotone additive models. We apply the new methodology to the problem of estimating levels of prenatal alcohol exposure that are associated with clinically significant cognitive defects in children.
Crump (1984; 1995) introduced the concept of a benchmark dose (BMD) and BMDL. Methods were developed first for parametric dose-response modeling of data obtained from designed experiments. Budtz-Jorgensen et al. (2001) proposed use of a dose-response model with a continuous exposure variable arising in epidemiological studies. Such a model can adjust for confounders through regression on covariates or propensity scores. They derive an exact confidence interval for the BMD estimator based on a linear dose-response model with Gaussian errors. In recent work on parametric dose-response models, Aerts et al. (2020) defined a family of parametric dose-response functions and recommended model-averaged estimates based on a computationally intensive bootstrap procedure.
Parametric models are needed for data with few unique exposure values, such as those often arising from experimental studies. Epidemiological studies often involve a broad range of exposure with many unique values. This type of data supports more flexible non- and semi-parametric modeling to characterize the dose-response function. In epidemiological setups, such models should also enable one to mitigate biases due to the effect of confounding variables. Accordingly, various semi- and non-parametric approaches to benchmark dose estimation have been developed for epidemiological data with continuous and binary outcomes, although they may also be useful for experimental data with sufficiently rich design. Wheeler and Bailer (2012) fit models based on monotone P-splines in a Bayesian framework using MCMC. For binary outcomes, non-parametric methods for dose-response modeling were developed by Piegorsch et al. (2012; 2014) using isotonic regression, again relying on intensive bootstrap calculations for BMDL estimation; methods to deal with continuous responses were developed by Lin et al. (2015). Wheeler et al. (2015) present a non-parametric quantile regression-based method. These methods are all based on non-parametric estimation of the dose-response relation and associated computationally-intensive inference procedures, which typically require the dose-response model to be fit many times.
1.2 Contributions
We take a semi-parametric approach to benchmark dose analysis, combining the efficiency of the parameteric approaches with the flexibility of the non-parameteric methods. We base our approach on monotone generalized additive models to estimate the benchmark dose and calculate the BMDL; restricting attention to monotone dose-response curves means that in this framework, an increased exposure cannot be associated with a preferred response (Wheeler and Bailer, 2012). We describe a novel procedure for fitting monotone generalized additive models based on B-splines with coefficients parameterized in such a way that they yield estimates of a monotone curve. Laplace-approximate marginal likelihood (Wood, 2011) is applied for selection of the smoothing parameter (Wood et al., 2016). We then develop a fast and stable method for solving the non-linear equation defining the BMD estimate based on de Boor’s algorithm for B-splines and their derivatives (de Boor, 2001). This is made feasible by deriving bounds on the estimated BMD and then applying a reflective Newton line search (Coleman and Li, 1994). Finally, we introduce a novel method for BMDL calculation based on the inversion of a hypothesis test based on an approximate pivot derived from the nonlinear equation solved by the estimated BMD. This BMDL is shown to have superior empirical and computational properties to those of the usual Delta method and bootstrap-based BMDLs in the settings considered. Significant computational improvements for the bootstrap are achieved by sampling from an approximate posterior distribution for the parameters, which enables full bootstrap-based BMDL calculation using only a single fit of the dose-response model. The speed of the novel BMD estimation algorithm makes this Bayesian parametric bootstrap approach computationally attractive in typical applications.
1.3 Motivating application
Prenatal alcohol exposure (PAE) has been linked to a broad range of long-term cognitive and behavioral deficits (Jacobson et al., 2024). However, there is relatively little information about the level of PAE that is associated with clinically significant cognitive deficits. Attempts to define what constitutes excessive drinking have been both qualitative (Stratton et al., 1996; Hoyme, 2005; Chudley et al., 2005; Cook et al., 2016) and quantitative (Astley and Clarren, 2000; Hoyme et al., 2016), although Astley and Clarren (2000) emphasize that there is no “clear consensus on the amount of alcohol that can actually be toxic to the fetus.” Addressing this question has important clinical implications in terms of diagnosing children who have been adversely affected by prenatal alcohol exposure. We apply our method for semi-parameteric benchmark dose analysis to this problem using data from six United States National Institutes of Health-funded longitudinal studies in which expectant women were interviewed regarding their drinking behavior during pregnancy and their children followed through young adulthood and assessed using a variety of cognitive tests.
2 BENCHMARK DOSE ANALYSIS AND MONOTONE SPLINES
2.1 Benchmark dose analysis
Consider a response |$Y_i\in \mathbb {R}$| and let |$x_i\in \mathbb {R}$| represent the exposure for individual |$i=1,\ldots ,n$|. We consider the following dose-response model:
where |$\alpha$| is an intercept, |$f(x)$| is a strictly monotone decreasing function of the exposure, and |$g_j(z_j)$| are |$m$| functions of other covariates, |$z_j$|, for |$j=1,\ldots ,m$|. This specification assumes that an increase in |$x$| reduces the mean response, appropriate for settings where smaller values of the response are undesirable, as is the case with cognition scores; see Sections 1.3 and 5. Adaptations to settings where an increase in exposure leads to a greater response (and hence |$f(x)$| would be monotone increasing) are straightforward.
Let |$p_0\in (0,1)$|, let |$p_{+}\in (0,1-p_0)$|, and let |$\tau _0(x_{0},z)$| satisfy |$\mathbb {P}(Y\lt \tau _0(x_{0},z) ; x_{0},z) = p_0$| where |$\mathbb {P}(\cdot ;x,z)$| is the probability distribution for |$Y$| following the model (1) with covariates |$(x,z) \equiv (x,z_1,\ldots ,z_m)$|. Here |$x_{0}\in \mathbb {R}$| is a baseline exposure representing an unexposed subject, and is almost always set at |$x_{0}= 0$| in practice if the data contain unexposed subjects. The benchmark dose (BMD, denoted as |$x_{\tt {b}}$|) is defined as the value of |$x$| satisfying
The value |$p_{+}$| is called the benchmark response (BMR); see Crump (1995) and Akkaya Hocagil et al. (2024) for further details.
Under the additive model (1), |$x_{\tt {b}}$| is defined by the nonlinear equation |$U(x_{\tt {b}}) = 0$|, where
and |$c(p_0,p_{+})= \Phi ^{-1}(p_0+p_{+}) - \Phi ^{-1}(p_0)$| where |$\Phi (\cdot )$| is the standard normal cumulative distribution function. It is clear from (3) that defined this way, |$x_{\tt {b}}$| does not depend on |$\tau _0(x_{0},z)$| or |$z$|. An estimate, |$\widehat{x}_{\tt {b}}$|, of |$x_{\tt {b}}$| solves the nonlinear equation |$U_n(\widehat{x}_{\tt {b}}) = 0$|, where
and |$\widehat{f}(x)$| and |$\widehat{\sigma }$| are estimates of |$f(x)$| and |$\sigma$|.
2.2 Quantifying uncertainty in the benchmark dose
To address the uncertainty in estimation of the benchmark dose, it is common to base guidelines on the lower endpoint of an approximate |$95\%$| confidence interval for |$x_{\tt {b}}$|; this is referred to as the benchmark dose-lower, denoted by BMDL. Budtz-Jorgensen et al. (2001) derive the exact sampling distribution of |$\widehat{x}_{\tt {b}}$| in linear dose-response models and use it to give a formula for a BMDL. More generally, an explicit sampling distribution of |$\widehat{x}_{\tt {b}}$| is not available, so a common approach is to use the Delta method for variance estimation, and rely on asymptotic normality of |$\widehat{x}_{\tt {b}}$|. Computationally-intensive bootstrap-based BMDL estimation is also possible by refitting the dose-response model many times; see Moerbeek et al. (2004) for a review of existing methods. In Section 3.4, we derive a novel method for uncertainty quantification of the BMD based on semi-parametric dose-response models and asymptotic normality of the approximate pivot in Equation 4.
2.3 Monotone B-splines
The dose-response function is represented as a monotone spline function. This proves useful for both fitting the dose-response model and rapid computation of the BMD and BMDL, and ensures that a BMD exists for some choice of the BMR; see Section 3. The function |$f(x)$| can be written as
where |$b_{l,p}(x)$| is the |$l^{th}$| B-spline function of order |$p$| on the interval |$[a,b]\subset \mathbb {R}$| with knots |$\boldsymbol{t}= (t_1,\ldots ,t_{L+p})$| such that |$a \le t_1\le \cdots \le t_{L+p}\le b$|. Using cubic B-splines having |$p=4$|, we follow the standard practice of choosing a large value of |$L$| and controlling over-fitting through penalized estimation; see Wood et al. (2016) and Wood (2017) for details.
B-splines are defined recursively with |$b_{l,1}(x) = I(\left[t_l,t_{l+1}\right))$| and |$b_{l,p}(x) = \omega _{l,p-1}b_{l,p-1}(x) + (1-\omega _{l+1,p-1})b_{l+1,p-1}(x)$| where |$\omega _{l,p}(x) = (x-t_l)/(t_{l+p}-t_l) I(t_{l+p}\ne t_l)$|. These definitions are used during model fitting (Section 3.2) to compute the design matrix where the |$x$| at which the splines are to be evaluated are specified in advance.
The numerical methods of Section 3 (see also Algorithms 2 and 3 in Supplementary Materials Section A) require computation of the spline function (5) within an iterative root-finding procedure at many values of |$x$| that cannot be known in advance. In this setting, computing |$b_{1,p}(x),\ldots ,b_{L,p}(x)$| naively and computing |$f(x)$| using (5) is too computationally intensive. For each fixed value of |$x$|, the sum defining the value of the spline function |$f(x)$| has only |$p+1$| nonzero terms. de Boor’s algorithm (de Boor, 2001) computes the spline function |$f(x)$| efficiently by ignoring computations known to add zero to the overall sum because the associated basis functions evaluate to zero; see Algorithm 1 in Supplementary Materials Section A. This greatly increases the speed of the computations involved in estimating the BMD and hence makes the parametric bootstrap of Section 3.4 practically feasible.
The derivative of a B-spline basis function is given by
and the derivative of the spline function (5) is:
From this, it can be seen that |$\beta _{1} \lt \cdots \lt \beta _{L}$| is a sufficient condition for |$f(x)$| to be strictly monotone decreasing, and this is straightforward to ensure during model fitting via re-parameterization as introduced by Pya and Wood (2015); see Section 3.2. Further, because the derivatives of a spline function are also spline functions (linear combinations of spline basis functions), de Boor’s algorithm is also applied to obtain the derivatives required to implement Newton’s method, leading to further computational benefits.
3 INFERENCES ABOUT THE BENCHMARK DOSE
3.1 Existence of the benchmark dose
The benchmark dose |$x_{\tt {b}}$|, defined as the solution to |$U(x)=0$| where |$U(x)$| is given by (2), is not guaranteed to exist for every choice of |$x_{0},p_0,p_{+}$|. The specification of |$x_{0},p_0,p_{+}$| should be context-dependent; see Haber et al. (2018). If the data involve a large number of unexposed subjects, as is the case in our motivating example on prenatal alcohol exposure, then setting |$x_{0}=0$| is reasonable. Otherwise, this is an important modeling choice; see Whitney and Ryan (2013) for a detailed discussion.
Assumption 1 is sufficient to ensure that for some |$p_0,p_{+}$|, there exists a value of |$x_{\tt {b}}$| satisfying |$U(x_{\tt {b}})=0$|:
(a) The unknown function |$f(x)$| is continuously differentiable with
|$\sup _{x\in \mathbb {R}}|f^\prime (x)|\lt \infty$| and |$f^\prime (x)\lt 0$| for every |$x\in \mathbb {R}$|; and (b) with probability 1 the estimated function |$\widehat{f}(x)$| is continuously differentiable with |$\sup _{x\in \mathbb {R}}|\widehat{f}^\prime (x)|\lt \infty$| and |$\widehat{f}^\prime (x)\lt 0$| for every |$x\in \mathbb {R}$|.
Assumption 1 (a) states that higher exposure cannot be associated with an equal or better expected response. Assumption 1 (b) can be satisfied by applying a monontonic smoother, as we do in Sections 2.3 and 3.2. Since |$U(x_{0}) \lt 0$| from Assumption 1, by the intermediate value theorem |$U(x_{\tt {max}})\gt 0$| is a sufficient condition for the existence of |$x_{\tt {b}}\in (x_{0},x_{\tt {max}})$| satisfying |$U(x_{\tt {b}})=0$|. Further, a sufficient condition for the estimate |$\widehat{x}_{\tt {b}}$| to exist is |$U_n(x_{\tt {max}})\gt 0$|. In cases where the underlying dose-response curve is not monotonic or when |$U_n(x_{\tt {max}}) \lt 0$|, a benchmark dose may not exist. This is important to consider in applications.
3.2 Fitting the dose-response model
Under the B-spline representation (5), the dose-response model (1) is:
Monotonicity can be enforced by reparameterization; see Pya and Wood (2015) for transformations that yield monotone increasing and decreasing functions as well as other shape constraints. We set |$\beta ^\tt {c}_1 = \beta _1$| and |$\beta ^\tt {c}_l = \beta ^\tt {c}_{l-1} - \exp \left(\beta _l\right),l=2,\ldots ,L$|, which guarantees that |$\beta ^\tt {c}_{1} \lt \cdots \lt \beta ^\tt {c}_{L}$| and hence that |$f(x)$| and its estimate, |$\widehat{f}(x) = b_{1,4}(x)\widehat{\beta ^\tt {c}}_1 + \cdots + b_{L,4}(x)\widehat{\beta ^\tt {c}}_L$|, are monotone decreasing. The unknown parameters to be estimated are |$\boldsymbol{\Psi }= (\alpha ,\boldsymbol{\beta },\boldsymbol{\gamma }, \sigma )\in \mathbb {R}^D$| where |$\alpha \in \mathbb {R},\boldsymbol{\beta }= (\beta _1,\ldots ,\beta _L)\in \mathbb {R}^L,$| |$\boldsymbol{\gamma }= (\gamma _{11},\ldots ,\gamma _{mL})\in \mathbb {R}^{mL}$| and |$\sigma \gt 0$|, so that |$D= 2 + L(1 + m)$|.
Define the design matrices |$\boldsymbol{B}= \left(B_{il}\right)$| with |$B_{il} = b_{l,4}(x_i)$| and |$\boldsymbol{Z}= \left[\boldsymbol{Z}_1 : \cdots : \boldsymbol{Z}_m\right]$| with |$\boldsymbol{Z}_{j} = \left(Z_{jil}\right)$| and |$Z_{jil} = b_{l,4}(z_{ij})$|. Further define |$\boldsymbol{\phi }= (\tau ,\boldsymbol{\lambda }) \in \mathbb {R}^{s}$| where |$s= m+2,\tau = -2\log \sigma$|, and |$\boldsymbol{\lambda }= (\lambda _1,\ldots ,\lambda _{m+1})$| where |$\lambda _j\in \mathbb {R}$| are (log) smoothing penalty parameters to be estimated. Finally, let |$\boldsymbol{S}_1,\boldsymbol{S}_2,\ldots ,\boldsymbol{S}_{m+1}$| be matrices of integrated squared B-spline second derivatives, computed according to the algorithm of Wood (2017), and let |$\boldsymbol{S}_{\boldsymbol{\lambda }} = \text{diag}\left(e^{\lambda _2}\boldsymbol{S}_2,\ldots ,e^{\lambda _{m+1}}\boldsymbol{S}_{m+1}\right)$|. A penalized joint (negative) log-likelihood for the unknown parameters, |$(\boldsymbol{\Psi },\boldsymbol{\phi })$|, is
and a marginal likelihood for |$\boldsymbol{\phi }$| is
We define the profile maximum likelihood estimator as |$\widehat{\boldsymbol{\Psi }}(\boldsymbol{\phi }) = \text{argmin}_{\boldsymbol{\Psi }}\ell (\boldsymbol{\Psi },\boldsymbol{\phi })$|, and the Hessian matrix |$\boldsymbol{H}(\boldsymbol{\phi }) = -\partial ^2_{\boldsymbol{\Psi }}\ell (\widehat{\boldsymbol{\Psi }}(\boldsymbol{\phi }),\boldsymbol{\phi })$|. Inferences about |$\boldsymbol{\phi }$| are to be based on the Laplace-approximate maximum marginal likelihood estimator, |$\widehat{\boldsymbol{\phi }}= \text{argmax}_{\boldsymbol{\phi }}\mathcal {L}_{\tt {LA}}(\boldsymbol{\phi })$|, where
is a Laplace approximation to |$\mathcal {L}(\boldsymbol{\phi })$|. Point estimates of |$\boldsymbol{\Psi }$| are then given as |$\widehat{\boldsymbol{\Psi }}(\widehat{\boldsymbol{\phi }})$|, and a point estimate of the unknown dose-response function at any |$x\in \mathbb {R}$| is |$\widehat{f}(x) = b(x)^{\textsf {T}}\widehat{\boldsymbol{\beta }}_\tt {c}$| with |$b(x) = (b_{1,4}(x),\ldots ,b_{L,4}(x))^{\textsf {T}}$|. Estimates |$\widehat{g}_j$| for each |$g_j,j=1,\ldots ,m$| are analogously obtained. Uncertainty quantification is discussed in Section 3.4.
Spline calculations are carried out using the Rcpp framework (Eddelbuettel and Francois, 2011). The Template Model Builder (TMB) framework (Kristensen et al., 2016) is used for automatic differentiation and Laplace approximation, yielding efficient computation of |$\log \mathcal {L}_{\tt {LA}}(\boldsymbol{\phi })$| and |$\partial _{\boldsymbol{\phi }}\log \mathcal {L}_{\tt {LA}}(\boldsymbol{\phi })$|. We compute |$\widehat{\boldsymbol{\phi }}$| via quasi-Newton optimization based on these quantities. The linear constraints |$\sum _{i=1}^n\widehat{f}(x_i) = 0$| and |$\sum _{i=1}^n\widehat{g}(z_{ij}) = 0,j=1,\ldots ,m$| are imposed directly, without “constraint-absorbing” reparameterizations (Stringer, 2024). Note that because |$U_n(\cdot )$| is invariant to the addition of constant terms to |$\widehat{f}$|, estimation of |$x_{\tt {b}}$| is invariant to the choice of linear constraints. However, such constraints are required to ensure the identifiability of |$\alpha$| as well as more than one |$f$| and |$g$| together, so they still must be used.
3.3 Computational method for benchmark dose estimation
Given a fitted dose-response model, we obtain the estimate, |$\widehat{x}_{\tt {b}}$|, defined as the solution to a non-linear equation, |$U_n(\widehat{x}_{\tt {b}})=0$|, via a reflective Newton line search. The efficiency and stability of the method are critically important in the parametric bootstrap of Section 3.4, where we repeat the procedure thousands of times to calculate a benchmark dose lower limit. Efficient computation of the required B-spline functions and derivatives to implement the Newton iteration is facilitated by de Boor’s algorithm (de Boor, 2001), and a stable iteration is achieved by bounding the solution to a known interval.
Bounds on |$\widehat{x}_{\tt {b}}$| are obtained by noting that Assumption 1 (b) guarantees that |$U_n(\cdot )$| is continuous and strictly monotonic, and further that |$U_n(x_{0}) = -c(p_0,p_{+})\lt 0$|. We therefore check the condition that |$U_n(x_{\tt {max}}) \gt 0$| in applications. If this constraint is satisfied, then by definition of |$\widehat{x}_{\tt {b}}$| and the intermediate value theorem, |$x_{0}\lt \widehat{x}_{\tt {b}}\lt x_{\tt {max}}$|. We then augment the typical Newton iteration with the reflective transformation given by Coleman and Li (1994), using |$(x_{0},x_{\tt {max}})$| as the required interval within which the solution is known to lie. The full algorithm is given in Algorithm 2 in Supplementary Materials Section A, which depends further on de Boor’s algorithm (Algorithm 1 in Supplementary Materials Section A).
3.4 Computation of benchmark dose lower limits
A method for computing a |$95\%$| confidence interval for |$x_{\tt {b}}$| is required since the lower endpoint of a |$95\%$| confidence interval for |$x_{\tt {b}}$| defines the benchmark dose lower limit (BMDL); see Section 2.2. We introduce a novel BMDL based on an approximate pivot obtained from |$U_n$|, and give a fast and stable method for computing it. This computational method draws on and extends the efficient framework for computation of |$\widehat{x}_{\tt {b}}$| (Section 3.3). Accordingly, we also provide a very efficient Bayesian parametric bootstrap that uses this framework to drastically speed up bootstrap BMDL calculation by avoiding the re-fitting of the dose-response model. We further compare to the more usual but undesirable Delta method based on asymptotic normality of |$\widehat{x}_{\tt {b}}$| directly. Further details for the bootstrap and Delta method are given in Supplementary Materials Section C.
Inspection of Equation 4 yields the approximate pivot |$U_n(x_{\tt {b}})^2/V_n(\widehat{x}_{\tt {b}})\overset{\cdot }{\sim }\chi ^2_1$|, where
and |$\Sigma (\widehat{\boldsymbol{\beta }}_\tt {c})= \text{Cov}(\widehat{\boldsymbol{\beta }})$|. Denote the entire vector of spline weights and variance/smoothing parameters by |$\boldsymbol{\theta }= (\boldsymbol{\Psi },\boldsymbol{\phi })$|, with estimate |$\widehat{\boldsymbol{\theta }}= (\widehat{\boldsymbol{\Psi }},\widehat{\boldsymbol{\phi }})$|. We follow a similar method to Wood et al. (2016) and draw samples, |$\lbrace \boldsymbol{\theta }_j\rbrace _{j=1}^{M}$| where |$M\in \mathbb {N}$|, from the approximate posterior |$\boldsymbol{\theta }| \boldsymbol{Y} \overset{\cdot }{\sim }\text{N}(\widehat{\boldsymbol{\theta }},\boldsymbol{H}^{-1}(\widehat{\boldsymbol{\theta }}))$| using the method of Rue (2001). We use the sample covariance matrix of |$\boldsymbol{\beta }_\tt {c}^1,\ldots ,\boldsymbol{\beta }_\tt {c}^M$|, |$\widehat{\Sigma }(\widehat{\boldsymbol{\beta }}_\tt {c})= \text{Cov}(\boldsymbol{\beta }_\tt {c}^1,\ldots ,\boldsymbol{\beta }_\tt {c}^M)$|, to estimate |$\Sigma (\widehat{\boldsymbol{\beta }}_\tt {c})$|. With these quantities available, we define the pivot-based BMDL, |$\widehat{x}_{\tt {l}}^\tt {piv}$|, as |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}\in \inf \widehat{\mathcal {C}}_n(\alpha )$|, where
Computation of |$\widehat{x}_{\tt {l}}^\tt {piv}$| is more involved than that of |$\widehat{x}_{\tt {b}}$|. We utilize another application of the reflective Newton line search coupled with de Boor’s algorithm as follows: Define the function |$\kappa _n(x) = U_n(x)^2 - V_n(x)\chi ^{2}_{1,1-\alpha /2}$| and note that |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}\in \inf \left\lbrace x:\kappa _n(x)=0\right\rbrace$|. Accordingly, we again apply a reflective Newton line search to find an appropriate zero of |$\kappa _n$|. The required derivative is |$\kappa _n^\prime (x) = 2U_n(x)U_n^\prime (x) - V_n^\prime (x)\chi ^{2}_{1,\alpha }$| where |$V_n^\prime (x) = -2b^\prime (x)^{\textsf {T}}\Sigma (\widehat{\boldsymbol{\beta }}_\tt {c})\left(b(x_{0})-b(x)\right)$|. Both the B-spline basis function vectors, |$b(x)$|, and their vectors of derivatives, |$b^\prime (x)$|, are computed directly using the recursions described in Section 2.3. Bounds are obtained by observing that |$\kappa _n(x_{0}) = c(p_0,p_{+})^2\gt 0$| and |$\kappa _n(\widehat{x}_{\tt {b}}) = - V_n(x)\chi ^{2}_{1,\alpha }\lt 0$|, and applying the intermediate value theorem to conclude that there exists |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}\in (x_{0},\widehat{x}_{\tt {b}})$|. Although such a root may or may not not be unique, we have not observed problems with convergence or stability in our experiments (Section 4) or data analysis (Section 5). The full algorithm is given in Algorithm 3 in Supplementary Materials Section A.
With the sample |$\lbrace \boldsymbol{\theta }_j\rbrace _{j=1}^{M}$| already efficiently drawn from the approximate parameter posterior |$\boldsymbol{\theta }| \boldsymbol{Y} \sim \text{N}(\widehat{\boldsymbol{\theta }},\boldsymbol{H}^{-1}(\widehat{\boldsymbol{\phi }}))$|, an efficient an parametric bootstrap BMDL, |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$|, is obtained immediately by repeating the procedure introduced in Section 3.3 for each sampled |$\boldsymbol{\theta }_j$|. This yields a bootstrap sample |$\lbrace (\widehat{x}_{\tt {b}})_j\rbrace _{j=1}^{M}$|, and |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| is defined as the lower |$2.5\%$| quantile of this sample. A BMDL estimated by the Delta method is given by
Further details about |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| and |$\widehat{x}_{\tt {l}}^\Delta$| are given in Supplementary Materials Section C.
4 EMPIRICAL PERFORMANCE AND COMPUTATIONAL CONSIDERATIONS
We compare via simulations the empirical performance of our proposed BMDL |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| with that of the parametric bootstrap |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| and the Delta method |$\widehat{x}_{\tt {l}}^\Delta$|; see Section 3.4. All computations are performed using the semibmd R package, which implements the methods described in Section 3. We simulate 100 000 replicates from the true dose-response function |$f(x) = \exp (-sx)$| for varying |$s\gt 0$|. This gives a monotone curve of varying steepness, where smaller |$s$| yields a flatter curve and hence a more difficult estimation problem. We report here selected results on empirical coverage and computation times; extended results are available in Supplementary Materials Section B.
We report empirical coverage probabilities and relative computation times for each simulation. Empirical coverage probabilities are defined as the proportion of simulated BMDLs that were less than the true BMD. The expected such proportion is |$97.5\%$|, reflecting the definition of the BMDL as the lower limit of an equal-tailed |$95\%$| confidence interval for |$x_{\tt {b}}$|. Monte Carlo confidence intervals for the empirical coverages are also shown, and the number of simulations was chosen to be high enough for these to be very narrow.
Relative computation times are defined as the number of seconds of total computation time required to compute each of |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| and |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| divided by the number of seconds of total computation time required to compute |$\widehat{x}_{\tt {l}}^\Delta$| for each simulated set of data. Such times are therefore unitless. Monte Carlo confidence intervals are based on the central limit theorem and hence we simply report means and standard errors.
Figure 1 shows the empirical coverage rates of |$\widehat{x}_{\tt {l}}^\Delta ,\widehat{x}_{\tt {l}}^\tt {\,\,piv}$|, and |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$|. The latter based on 1000 parametric bootstrap replications for each of the 100 000 simulated replicates; this massive amount of computation was made feasible by the speed and stability of the |$\widehat{x}_{\tt {b}}$| computations described in Sections 3.3 and 3.4. The performance of |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| is satisfactory: coverage is at or just above the nominal level in most cases. In the more difficult cases of flat dose-response function (low |$s$|, top 2 rows) or higher residual standard deviation |$\sigma$| (rightmost column), undercoverage is observed at the lower sample sizes |$n=200$| and |$n=500$|, but this corrects in all cases for |$n=1000$|. The bootstrap |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| is sometimes more conservative than |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| but otherwise generally agrees with it, at a cost of roughly |$30\times$| more computation time (see Table 1). The Delta method |$\widehat{x}_{\tt {l}}^\Delta$| is highly conservative, and often gives coverage of |$100\%$| or very near. All methods were too conservative in the case of a flat dose-response curve with higher |$\sigma$| (top right corner).
Empirical coverage proportions, |$\widehat{x}_\tt {l}\le x_{\tt {b}}$|, from 100 000 simulations described in Section 4.
Average and empirical standard error of computation time for |$\widehat{x}_{\tt {l}}^\tt {\,\,piv},\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| relative to |$\widehat{x}_{\tt {l}}^\Delta$|. Numbers shown are average and (standard error) of the computational time in seconds of |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| or |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| divided by that of |$\widehat{x}_{\tt {l}}^\Delta$|, and are hence unitless relative computation times.
| |$n = 200$| | |$n = 500$| | |$n = 1000$| | |||||
|---|---|---|---|---|---|---|---|
| |$s$| | |$\sigma$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| |
| .1 | .1 | 1.02(0.00) | 33.80(0.05) | 1.01(0.00) | 33.47(0.05) | 1.01(0.00) | 33.34(0.05) |
| .2 | 1.02(0.00) | 34.60(0.06) | 1.02(0.00) | 33.98(0.06) | 1.02(0.00) | 33.71(0.05) | |
| .5 | 1.02(0.00) | 35.50(0.07) | 1.02(0.00) | 35.58(0.06) | 1.02(0.00) | 35.78(0.06) | |
| .5 | .1 | 1.01(0.00) | 33.43(0.05) | 1.01(0.00) | 33.29(0.05) | 1.01(0.00) | 33.43(0.05) |
| .2 | 1.01(0.00) | 33.41(0.05) | 1.01(0.00) | 33.31(0.05) | 1.01(0.00) | 33.35(0.05) | |
| .5 | 1.02(0.00) | 33.90(0.05) | 1.01(0.00) | 33.55(0.05) | 1.01(0.00) | 33.44(0.05) | |
| 1 | .1 | 1.01(0.00) | 33.54(0.05) | 1.01(0.00) | 33.38(0.05) | 1.01(0.00) | 33.39(0.05) |
| .2 | 1.01(0.00) | 33.44(0.05) | 1.01(0.00) | 33.36(0.05) | 1.01(0.00) | 33.22(0.05) | |
| .5 | 1.01(0.00) | 33.52(0.04) | 1.01(0.00) | 33.46(0.05) | 1.01(0.00) | 33.48(0.05) | |
| 2 | .1 | 1.01(0.00) | 33.71(0.05) | 1.01(0.00) | 33.53(0.05) | 1.01(0.00) | 33.42(0.05) |
| .2 | 1.01(0.00) | 33.76(0.05) | 1.01(0.00) | 33.46(0.05) | 1.01(0.00) | 33.70(0.05) | |
| .5 | 1.01(0.00) | 33.66(0.05) | 1.01(0.00) | 33.69(0.05) | 1.01(0.00) | 33.51(0.05) | |
| 5 | .1 | 1.01(0.00) | 33.91(0.05) | 1.01(0.00) | 33.79(0.05) | 1.01(0.00) | 34.01(0.05) |
| .2 | 1.01(0.00) | 34.14(0.05) | 1.01(0.00) | 33.77(0.05) | 1.02(0.01) | 33.88(0.05) | |
| .5 | 1.01(0.00) | 34.07(0.05) | 1.01(0.00) | 33.97(0.05) | 1.01(0.00) | 34.03(0.05) | |
| |$n = 200$| | |$n = 500$| | |$n = 1000$| | |||||
|---|---|---|---|---|---|---|---|
| |$s$| | |$\sigma$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| |
| .1 | .1 | 1.02(0.00) | 33.80(0.05) | 1.01(0.00) | 33.47(0.05) | 1.01(0.00) | 33.34(0.05) |
| .2 | 1.02(0.00) | 34.60(0.06) | 1.02(0.00) | 33.98(0.06) | 1.02(0.00) | 33.71(0.05) | |
| .5 | 1.02(0.00) | 35.50(0.07) | 1.02(0.00) | 35.58(0.06) | 1.02(0.00) | 35.78(0.06) | |
| .5 | .1 | 1.01(0.00) | 33.43(0.05) | 1.01(0.00) | 33.29(0.05) | 1.01(0.00) | 33.43(0.05) |
| .2 | 1.01(0.00) | 33.41(0.05) | 1.01(0.00) | 33.31(0.05) | 1.01(0.00) | 33.35(0.05) | |
| .5 | 1.02(0.00) | 33.90(0.05) | 1.01(0.00) | 33.55(0.05) | 1.01(0.00) | 33.44(0.05) | |
| 1 | .1 | 1.01(0.00) | 33.54(0.05) | 1.01(0.00) | 33.38(0.05) | 1.01(0.00) | 33.39(0.05) |
| .2 | 1.01(0.00) | 33.44(0.05) | 1.01(0.00) | 33.36(0.05) | 1.01(0.00) | 33.22(0.05) | |
| .5 | 1.01(0.00) | 33.52(0.04) | 1.01(0.00) | 33.46(0.05) | 1.01(0.00) | 33.48(0.05) | |
| 2 | .1 | 1.01(0.00) | 33.71(0.05) | 1.01(0.00) | 33.53(0.05) | 1.01(0.00) | 33.42(0.05) |
| .2 | 1.01(0.00) | 33.76(0.05) | 1.01(0.00) | 33.46(0.05) | 1.01(0.00) | 33.70(0.05) | |
| .5 | 1.01(0.00) | 33.66(0.05) | 1.01(0.00) | 33.69(0.05) | 1.01(0.00) | 33.51(0.05) | |
| 5 | .1 | 1.01(0.00) | 33.91(0.05) | 1.01(0.00) | 33.79(0.05) | 1.01(0.00) | 34.01(0.05) |
| .2 | 1.01(0.00) | 34.14(0.05) | 1.01(0.00) | 33.77(0.05) | 1.02(0.01) | 33.88(0.05) | |
| .5 | 1.01(0.00) | 34.07(0.05) | 1.01(0.00) | 33.97(0.05) | 1.01(0.00) | 34.03(0.05) | |
Average and empirical standard error of computation time for |$\widehat{x}_{\tt {l}}^\tt {\,\,piv},\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| relative to |$\widehat{x}_{\tt {l}}^\Delta$|. Numbers shown are average and (standard error) of the computational time in seconds of |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| or |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| divided by that of |$\widehat{x}_{\tt {l}}^\Delta$|, and are hence unitless relative computation times.
| |$n = 200$| | |$n = 500$| | |$n = 1000$| | |||||
|---|---|---|---|---|---|---|---|
| |$s$| | |$\sigma$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| |
| .1 | .1 | 1.02(0.00) | 33.80(0.05) | 1.01(0.00) | 33.47(0.05) | 1.01(0.00) | 33.34(0.05) |
| .2 | 1.02(0.00) | 34.60(0.06) | 1.02(0.00) | 33.98(0.06) | 1.02(0.00) | 33.71(0.05) | |
| .5 | 1.02(0.00) | 35.50(0.07) | 1.02(0.00) | 35.58(0.06) | 1.02(0.00) | 35.78(0.06) | |
| .5 | .1 | 1.01(0.00) | 33.43(0.05) | 1.01(0.00) | 33.29(0.05) | 1.01(0.00) | 33.43(0.05) |
| .2 | 1.01(0.00) | 33.41(0.05) | 1.01(0.00) | 33.31(0.05) | 1.01(0.00) | 33.35(0.05) | |
| .5 | 1.02(0.00) | 33.90(0.05) | 1.01(0.00) | 33.55(0.05) | 1.01(0.00) | 33.44(0.05) | |
| 1 | .1 | 1.01(0.00) | 33.54(0.05) | 1.01(0.00) | 33.38(0.05) | 1.01(0.00) | 33.39(0.05) |
| .2 | 1.01(0.00) | 33.44(0.05) | 1.01(0.00) | 33.36(0.05) | 1.01(0.00) | 33.22(0.05) | |
| .5 | 1.01(0.00) | 33.52(0.04) | 1.01(0.00) | 33.46(0.05) | 1.01(0.00) | 33.48(0.05) | |
| 2 | .1 | 1.01(0.00) | 33.71(0.05) | 1.01(0.00) | 33.53(0.05) | 1.01(0.00) | 33.42(0.05) |
| .2 | 1.01(0.00) | 33.76(0.05) | 1.01(0.00) | 33.46(0.05) | 1.01(0.00) | 33.70(0.05) | |
| .5 | 1.01(0.00) | 33.66(0.05) | 1.01(0.00) | 33.69(0.05) | 1.01(0.00) | 33.51(0.05) | |
| 5 | .1 | 1.01(0.00) | 33.91(0.05) | 1.01(0.00) | 33.79(0.05) | 1.01(0.00) | 34.01(0.05) |
| .2 | 1.01(0.00) | 34.14(0.05) | 1.01(0.00) | 33.77(0.05) | 1.02(0.01) | 33.88(0.05) | |
| .5 | 1.01(0.00) | 34.07(0.05) | 1.01(0.00) | 33.97(0.05) | 1.01(0.00) | 34.03(0.05) | |
| |$n = 200$| | |$n = 500$| | |$n = 1000$| | |||||
|---|---|---|---|---|---|---|---|
| |$s$| | |$\sigma$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| | |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| |
| .1 | .1 | 1.02(0.00) | 33.80(0.05) | 1.01(0.00) | 33.47(0.05) | 1.01(0.00) | 33.34(0.05) |
| .2 | 1.02(0.00) | 34.60(0.06) | 1.02(0.00) | 33.98(0.06) | 1.02(0.00) | 33.71(0.05) | |
| .5 | 1.02(0.00) | 35.50(0.07) | 1.02(0.00) | 35.58(0.06) | 1.02(0.00) | 35.78(0.06) | |
| .5 | .1 | 1.01(0.00) | 33.43(0.05) | 1.01(0.00) | 33.29(0.05) | 1.01(0.00) | 33.43(0.05) |
| .2 | 1.01(0.00) | 33.41(0.05) | 1.01(0.00) | 33.31(0.05) | 1.01(0.00) | 33.35(0.05) | |
| .5 | 1.02(0.00) | 33.90(0.05) | 1.01(0.00) | 33.55(0.05) | 1.01(0.00) | 33.44(0.05) | |
| 1 | .1 | 1.01(0.00) | 33.54(0.05) | 1.01(0.00) | 33.38(0.05) | 1.01(0.00) | 33.39(0.05) |
| .2 | 1.01(0.00) | 33.44(0.05) | 1.01(0.00) | 33.36(0.05) | 1.01(0.00) | 33.22(0.05) | |
| .5 | 1.01(0.00) | 33.52(0.04) | 1.01(0.00) | 33.46(0.05) | 1.01(0.00) | 33.48(0.05) | |
| 2 | .1 | 1.01(0.00) | 33.71(0.05) | 1.01(0.00) | 33.53(0.05) | 1.01(0.00) | 33.42(0.05) |
| .2 | 1.01(0.00) | 33.76(0.05) | 1.01(0.00) | 33.46(0.05) | 1.01(0.00) | 33.70(0.05) | |
| .5 | 1.01(0.00) | 33.66(0.05) | 1.01(0.00) | 33.69(0.05) | 1.01(0.00) | 33.51(0.05) | |
| 5 | .1 | 1.01(0.00) | 33.91(0.05) | 1.01(0.00) | 33.79(0.05) | 1.01(0.00) | 34.01(0.05) |
| .2 | 1.01(0.00) | 34.14(0.05) | 1.01(0.00) | 33.77(0.05) | 1.02(0.01) | 33.88(0.05) | |
| .5 | 1.01(0.00) | 34.07(0.05) | 1.01(0.00) | 33.97(0.05) | 1.01(0.00) | 34.03(0.05) | |
Table 1 summarizes computational aspects of the procedure. Shown are the (unitless) relative computation times for |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| and |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| relative to |$\widehat{x}_{\tt {l}}^\Delta$|, which should be the fastest. The proposed pivot-based BMDL |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| is essentially just as fast as the Delta method with relative computation time only |$1\%$| higher; the raw computation times for both are on the order of microseconds on the hardware used. Because the Bayesian parametric bootstrap given in Section 3.4 does not require re-fitting the dose-response model, it is also fast, and with 1000 samples, it takes only around 30 times more computing time than |$\widehat{x}_{\tt {l}}^\Delta$| and |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$|. However, its performance is broadly comparable to that of |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$|. The conclusion is that |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| achieves better results than |$\widehat{x}_{\tt {l}}^\Delta$| in the same computational time, and equally satisfactory results to |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| in about |$30\times$| less computational time.
5 BENCHMARK DOSE ANALYSIS OF PRENATAL ALCOHOL EXPOSURE
We consider data from six longitudinal cohort studies of alcohol consumption by pregnant women conducted in Detroit (Jacobson et al., 1993), Pittsburgh (Day et al., 1991; Richardson et al., 1999), Atlanta (Coles et al., 2006; Brown et al., 1998), and Seattle (Streissguth et al., 1981). The research question of interest is to quantify levels of prenatal alcohol exposure (PAE) associated with the development of clinically important cognitive deficits in children. In these six cohort studies, children were followed from infancy through young adulthood and investigators administered a range of neuropsychological tests to assess IQ and four domains of cognitive function: learning and memory, executive function, and academic achievement in reading and mathematics. To obtain an overall cognitive function score for each child, a structural equation model was fitted for each cohort. We then used the estimated cognitive function score as the outcome measure in the below analysis. For details on this approach, see Jacobson et al. (2024) and Akkaya Hocagil et al. (2024).
Data on maternal alcohol consumption were summarized in terms of average alcohol intake per day (ounces of absolute alcohol, (AA)/day) during pregnancy, and data on a broad range of potential confounders were collected in these cohorts. Since each cohort provided a somewhat different set of confounding variables, for these data Jacobson et al. (2024) modeled the exposure variable, average alcohol intake per day, as a function of the potential confounders via a linear regression model and use this to estimate a propensity score (Rosenbaum and Rubin, 1983). We used these propensity scores as covariates (Rosenbaum and Rubin, 1983; Imbens and Hirano, 2004) in the following dose-response model.
Let |$Y_i$| represent the cognitive function score, |$x_i = \log (a_i+1)$| where |$a_i$| is the average alcohol intake per day during pregnancy, |$j_i\in \left\lbrace 1,\ldots ,6\right\rbrace$| index the cohort to which subject |$i$| belongs, and |$z_{i} = s_{j_i}$| denote the computed propensity scores for subject |$i=1,\ldots ,2226$| in cohort |$j_i$|. A log-transformation was applied to |$a_i$| to reduce the influence of very high exposure values on the fitted model. The dose-response model is then
where |$f$| is an unknown, smooth monotone dose-response function, and |$g_j,j=1,\ldots ,6$| are unknown smooth functions. We apply the basis expansions of Section 2.3 and fit the model using the methods described in Section 3.2. We compute the estimated BMD, |$\widehat{x}_{\tt {b}}$|, using the method of Section 3.3, and the BMDL, |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$|, using the methods of Section 3.4. We also compute |$\widehat{x}_{\tt {l}}^\Delta$| and |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| for comparison. We used |$p_{+}=.01$| and |$p_0=.025$|; see Akkaya Hocagil et al. (2024) for the rationale for these choices.
Figure 2 represents the estimated dose-response curve obtained from the fitted monotone additive model along with the estimated BMD and corresponding BMDLs; see Figure 1 in Supplementary Materials Section D for model diagnostics. The estimated BMD was |$\widehat{x}_{\tt {b}}=0.764$| and the three BMDLs were |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}= 0.356$|, |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}= 0.310$|, and |$\widehat{x}_{\tt {l}}^\Delta = 0.175$|. The Delta method estimate, |$\widehat{x}_{\tt {l}}^\Delta$|, is much closer to zero than the other two, and hence communicates a much more conservative clinical recommendation; combined with its empirical coverage being too high in simulations, we remark that this can be regarded as a disadvantage of the Delta method over the two proposed BMDLs.
(a) Estimated dose-response curve (—) and uncertainty bands computed as pointwise |$95\%$| credible intervals from 100 000 samples from the posterior of the fitted curve (- - -). Vertical lines are the estimated BMD, |$\widehat{x}_{\tt {b}}$| (—) and the BMDLs: |$\widehat{x}_{\tt {l}}^\Delta$| (|$\cdots$|), |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| (- - -), and |$\widehat{x}_{\tt {l}}^\tt {\,\,boot}$| (|$-\cdot -$|). (b) 100 000 samples from the approximate posterior distribution of |$\widehat{x}_{\tt {b}}$| (|$\blacksquare$|) and normal approximation to this distribution (—) upon which |$\widehat{x}_{\tt {l}}^\Delta$| is based.
The total amount of time that the procedure took to fit the model, estimate the BMD, and compute the three BMDLs based on 100 000 bootstrap iterations was 396 seconds on a 2021 M1 Macbook Pro with 64 GB of RAM. Model fitting and estimating the BMD took 3.93 seconds and 30 microseconds, respectively. Estimating the covariance matrix via posterior sampling and computing the pivot BMDL took 1.18 and 0.4 seconds, respectively. The computation time for 100 000 iterations of the parametric bootstrap for |$\widehat{x}_{\tt {b}}$| was the longest at 390 seconds, and produced a result close to |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| which was computed nearly |$1000\times$| faster. Fitting the model, estimating the BMD, and computing the pivot-based BMDL |$\widehat{x}_{\tt {l}}^\tt {\,\,piv}$| took a total of 5.53 seconds.
6 DISCUSSION
The method described in Section 3.2 for fitting the dose-response model can be used to fit monotone generalized additive models in which the response is binary, count, categorical, and generally has any smooth density; see Pya and Wood (2015). The quadratic term in Equation 8 would simply be replaced by the negative log-likelihood, and then the Laplace approximation described in Section 3.2 applies unchanged; see Wood et al. (2016), Kristensen et al. (2016), and Brooks et al. (2017). For continuous data, Wood et al. (2016) show how this approach can also be used to relax the assumption of constant variance in additive models with Gaussian responses. Our proposed approach for fitting dose-response models is easily adapted to other settings.
Extension of the present methods for BMD(L) determination for continuous outcomes to binary, count, and categorical outcomes is complicated only by the BMD itself having different definitions in each of these types of data. For example, for a binary outcome |$Y_i$|, Crump (1995, Equation (1)) defines the BMD as the excess risk; if |$P(Y_i = 1;x) = \mu (x)$| then the BMD |$x_{\tt {b}}$| is defined by |$ \mu (x_{\tt {b}}) - p_0= p_{+}.$| If our approach is applied to fit the logistic dose-response model |$\log [\mu (x)/\lbrace 1-\mu (x)\rbrace ] = f(x) + g_1(z_1) + \cdots + g_d(z_d)$|, then the reflective Newton approach of Section 3.3 could be applied to solve |$\widehat{\mu }(x_{\tt {b}}) - p_0= p_{+}$|. A pivot-based BMDL may also be possible to derive. These questions are interesting for future work; one challenge to overcome is that the BMD will depend on the values of any other covariates in the model, complicating analysis and reporting.
We recommend avoiding the Delta method (|$\widehat{x}_{\tt {l}}^\Delta$|; Section 3.4) for forming confidence intervals for |$\widehat{x}_{\tt {b}}$| in general. The estimate |$\widehat{x}_{\tt {b}}$| is a highly non-linear function of the MLE, and should not be expected to have a sampling distribution that is close to Gaussian. An example of this is shown in Figure 2(b) for the PAE data analysis, in which the approximate posterior distribution of |$\widehat{x}_{\tt {b}}$| is extremely non-normal, leading to a substantial overestimation of uncertainty and a low value of |$\widehat{x}_{\tt {l}}^\Delta$|. This recommendation echoes that of Moerbeek et al. (2004).
A substantial motivation for using semi-parametric approaches is to avoid specifying the form of the dose-response function. In designed experiments, there may be scientific reasons to favor a specific parametric model, but this is less often the case in epidemiological studies. Regarding the motivating study of the effects of PAE on child cognition, Jacobson et al. (2024) argue that the dose-response curve is plausibly non-linear, and its shape should be flat at low levels of exposure and then steeper beyond a certain level. Such situations in which the shape of the dose-response curve is partially specified are particularly amenable to analysis using the semi-parametric methods we develop in this paper.
The estimated dose-response curve shown in Figure 2(a) for the PAE data analysis is very flat at low exposure values. This appears to be due to a large amount of noise relative to signal in the PAE data themselves, suggesting that a single exposure measure may not be sufficiently informative (Jacobson et al., 2024). Akkaya Hocagil et al. (2024) perform a benchmark dose analysis for PAE using a bivariate exposure that incorporates both frequency and severity of drinking. In future work, we plan to develop an efficient computational method for the semi-parametric bivariate exposure case, making use of the fast computational methods developed in the present manuscript to improve the practical application of benchmark dose estimation with bivariate exposure.
Although the present paper presents methodology that is specifically tied to computations related to benchmark dose estimation and its associated lower confidence limit, the methods could be applied to more general inverse estimation problems. In such problems, the object of interest is a point on the |$x$|-axis in a semi-parametric regression model, and is hence obtained by solving a nonlinear equation. Inverse problems arise in a wide range of applied settings, for example, chemical calibration or medical imaging, and it would be interesting to apply similar approaches to those used in the present paper in other such problems.
ACKNOWLEDGEMENT
The authors wish to recognize the principal investigators, Claire D. Coles (Atlanta Cohort Study 1 and Atlanta Cohort Study 2), Gale A. Richardson (Pittsburgh Cohort Study 2), the late Ann P. Streissguth, Heather Carmichael Olson (Seattle Cohort Study), Nancy L. Day (Pittsburgh Cohort Study 1), and their study team members for their generous contribution of the PAE data used in Section 5.
FUNDING
This research was supported by Grant No. R01 AA025905 from the National Institute on Alcohol Abuse and Alcoholism (NIAAA)/National Institutes of Health (NIH) and the Lycaki-Young Fund from the State of Michigan. Data collection for the Detroit cohort was supported by NIAAA/NIH grants R01 AA06966, R01 AA09524, and P50 AA07606 and National Institute on Drug Abuse (NIDA)/NIH grant R21 DA021034; for the Pittsburgh cohorts, by NIAAA/NIH grants R01 AA06390, R01 AA06666, R01 AA14215, R01 AA18116, National Institute on Child Health and Human Development/NIH grant HD036890, and NIDA/NIH grants R01 DA00090, R01 DA03209, R01 DA03874, and R01 DA17786, R01 DA05460, R01 DA06839, R01-DA08916, and R01 DA12401; for the Atlanta cohorts, by NIAAA/NIH grants R01 AA08105, R01 AA10108, and R01 AA13272, NIDA/NIH grants R01 DA0737362 and R01 DA0737362, and Georgia State Contract DHR 427 93 05050762; and for the Seattle cohort, by NIAAA/NIH grant R01 AA001455. This work also supported by grants from the Australian Research Council Centre of Excellence for Mathematical and Statistical Frontiers to Louise M. Ryan (CE140100049), the Natural Sciences and Engineering Research Council of Canada to Richard J. Cook (RGPIN2017-04207), the Canadian Institutes of Health Research to Richard J. Cook and Louise M.Ryan (FRN PJT-180551). Alex Stringer was supported by grant RGPIN-03331-2023 from the Natural Sciences and Engineering Research Council of Canada..
CONFLICT OF INTEREST
None declared.
DATA AVAILABILITY
Data analyzed in Section 5 is available from the Harvard DataVerse (Akkaya Hocagil, 2024) at https://doi.org/10.7910/DVN/SCLH9W.
References
Open Data
Open Materials
