Benchmark dose profiles for bivariate exposures

1 INTRODUCTION

Benchmark dose (BMD) analysis aims to estimate a “dose” or level of exposure associated with an increase in the risk of an adverse response by a prespecified amount (called the benchmark response or BMR) (Budtz-Jørgensen et al., 2001; Crump, 1984; Erdreich & Stara, 1984). This so-called BMD is estimated by first modeling the dose–response relationship between the exposure and response and then finding the exposure level leading to the expected increase in risk. Uncertainty is addressed by the computation of a BMDL, which corresponds to a lower confidence limit on the true BMD. Government and private entities employ BMDs and BMDLs for risk assessment in both the United States and the European Union (General Accounting Office, 2001; OECD, 2006, 2014; U.S. Environmental Protection Agency, 2012). Following the introduction of the initial safe dose for methylmercury onto the US Environmental Protection Agency (EPA)'s Integrated Risk Information System (IRIS) in 1995, guidance and broad explanations regarding the application of BMD methodology emerged (European Centre for Disease Prevention and Control, 2011; European Food Safety Authority (EFSA), 2017; Klaassen, 2012; U.S. Environmental Protection Agency, 2000). Subsequently, BMD analysis became the standard method for establishing safe exposure thresholds within environmental toxicology and regulatory frameworks (National Research Council (US) Committee on the Toxicological Effects of Methylmercury, 2000). Haber et al. (2018) provide an excellent summary, including a discussion of the recent trend to use model averaging techniques to address issues of model uncertainty.

Methods for BMD analysis are well-established for settings involving a single exposure but methods for multidimensional exposures are less developed. Chen et al. (1990) developed a method to compute the BMDL for studies of the chemical mixture using a likelihood ratio technique (Chen et al., 1990), but BMDLs were obtained separately for each chemical so joint effects were not addressed. Zhu (2005) developed estimators for the BMD for quantal data, with two factors modeled in tandem, but where one of these factors was taken as time and viewed as a secondary nuisance variable (Zhu, 2005). Deutsch and Piegorsch (2012, 2013) described a methodology for conducting a full BMD analysis with two exposures for quantal and nonquantal data. All of these approaches require the specification of a parametric form for nonlinear effects and do not accommodate covariate adjustment. Several researchers have argued that parametric dose–response modeling is too restrictive for BMD estimation, but little work has been carried out for BMD analysis involving flexible dose–response models. Wheeler and Bailer (2012) use Bayesian methods and Gibbs sampling to fit monotone P-splines, while Piegorsch et al. (2012) developed nonparametric methods for BMD estimation for quantile modeling. Lin et al. (2015) developed nonparametric methods for BMD estimation with continuous responses. Wheeler et al. (2015) present a nonparametric estimator based on quantile regression. All these nonparametric and semiparametric methods are developed for assessing the risk of a single exposure.

In our motivating study, there is evidence that cognitive deficits associated with the prenatal alcohol exposure (PAE) do not simply depend on the total amount of alcohol consumed during pregnancy, but rather on both the amount of alcohol consumed on drinking occasions, and the frequency of drinking (days/week) during pregnancy (Jacobson et al., 2024). Moreover, exploratory analyses suggest the presence of nonlinear, interacting effects of these variables.

In this paper, we propose a flexible framework for BMD analysis in which the dose–response surface characterizing the effect of two continuous exposure variables on a continuous response is fitted based on a generalized additive model (GAM; Wood, 2017); confounding variables can be accommodated in the regression model directly or using propensity scores. We then show how to estimate BMD contours from the dose–response surface. An innovative parametric bootstrap method is then employed to compute corresponding BMDL contours. We use this flexible modeling approach to gain insights into the effects of dose/occasion and drinking frequency on the child's cognition and derive BMD and BMDL contours.

The remainder of the paper is organized as follows. In Section 2, we describe six US cohort studies designed to assess the relationship between maternal alcohol intake during pregnancy and cognitive function measured in the offspring. Dose–response modeling involving harmonized data from these six studies motivates our work. After briefly reviewing how BMDs are calculated for continuous outcomes assessed in epidemiological settings, Section 3 extends this standard analysis to a setting involving GAMs and describes a bootstrap approach for calculating BMDLs derived from the fitted semiparametric model. In Section 9, we describe how this can then be extended to characterize the joint effects of a bivariate exposure variable. In all our analyses, the effects of potential confounders are addressed through regression adjustment on propensity scores. In Section 10, we apply the methods to our motivating data set and draw insights, with further discussion, conclusions, and topics warranting further research presented in Section 11.

2 PRENATAL ALCOHOL AND CHILDHOOD COGNITION

Animal models and epidemiological studies have linked PAE to a broad range of long-term cognitive and behavioral deficits (Jacobson et al., 2024). However, there is relatively little consensus regarding the levels and patterns of PAE that are associated with clinically significant cognitive deficits. For these reasons, the US Office of the Surgeon General's (2005) recommendation (National Center on Birth Defects and Developmental Disabilities, Centers for Disease Control and Prevention (Internet), 2005) that all women abstain from alcohol during pregnancy continues to be the most appropriate guideline. However, there are compelling diagnostic reasons for studying the impact of dose/occasion and frequency of drinking on childhood cognition. This information, for example, can help clinicians differentiate between children affected by PAE and those with attention-deficit hyperactivity disorder (ADHD) and other conditions that give a similar clinical presentation (Peadon & Elliott, 2010). This distinction is critically important since the preferred treatment options are quite different (Petrenko & Alto, 2017). For example, methylphenidate is often used in treating ADHD, but it is generally ineffective for children impacted by the PAE (Mela et al., 2018). A more complete understanding of the etiology of cognitive limitations will help ensure timely and appropriate treatments are administered. This information may also provide some reassurance to women who may have inadvertently been exposed to low levels of alcohol before realizing they were pregnant. Traditionally, for diagnostic purposes, problematic exposure has been crudely characterized as a “pattern of excessive intake characterized by substantial, regular intake or heavy episodic drinking,” or the presence of social, legal, or medical problems related to drinking (Chudley et al., 2005; Cook et al., 2016; Hoyme, 2005; Stratton et al., 1996). Astley and Clarren (2000) suggested a quantitative criterion but emphasized that there is no “clear consensus on the amount of alcohol that can actually be toxic to the fetus.” More recently, risky drinking was defined as at least six drinks/week for $$>$$2 weeks, or three drinks on at least two occasions, but there is no empirical evidence of adverse effects on the offspring of women with these relatively low levels of exposure (e.g., a single drink/day over a 2-week period) or of women who did not binge drink (at least four to five drinks) on multiple occasions.

To help inform more evidence-based decision-making related to potentially risky drinking during pregnancy, the National Institutes of Health (NIH) funded a project in 2017 to combine data from six longitudinal cohort studies in order to obtain a sufficiently rich and large data set to rigorously explore the dose–response relationship between the PAE and child development. The cohorts, conducted in Detroit (Jacobson et al., 1993), Pittsburgh (Day et al., 1989; Richardson et al., 1999), Atlanta (Brown et al., 1998; Coles et al., 2006), and Seattle (Streissguth et al., 1981) had all been funded by the NIH between 1975 and 1993 and each recruited expectant women who were interviewed regarding their drinking behavior during pregnancy. Their offspring were then followed longitudinally from infancy through childhood and assessed periodically using a battery of neuropsychological tests to assess IQ and four domains of cognitive function: learning and memory, executive function, and academic achievement in reading and mathematics. The number of mother/infant pairs recruited for these cohorts ranged from 138 to 720, with a total combined sample size of 2226. Over the past several years, our team has published several articles focused on different aspects of the project. There have been several statistical methods publications addressing some of the analysis challenges, which have required nonstandard solutions. Topics addressed include how to apply meta-analysis methods to combine multiple outcome data from multiple studies (Akkaya Hocagil et al., 2022), how to use propensity scores in meta-analysis (Akkaya Hocagil et al., 2024) as well as several others (Dang, Ryan, Akkaya Hocagil et al., 2023; Dang, Ryan, Cook et al., 2023; Li et al., 2023). Akkaya Hocagil et al. (2021) extended the computation of propensity scores to account for the possibility that abstainers might differ from even very low-level drinkers and also allowed for long-tailed exposure distributions. All these various methodologies have been used to report on the substantive findings in the alcohol epidemiology literature. Jacobson et al. (2021) used the meta-analysis framework to establish that alcohol exposure had a similar magnitude of effect across a range of outcomes measured within all four cognitive domains. To address the dose–response question, Jacobson et al. (2024) describe the use of second-order confirmatory factor analysis, conducted separately within each cohort, to obtain an overall cognitive function score for each child, which was then used as the outcome measure in a dose–response analysis approach. This paper found significant evidence that cognitive deficits associated with the PAE depend not simply on the total amount of alcohol consumed during pregnancy but on both the amount of alcohol consumed on each drinking occasion and the frequency of drinking (days/week) during pregnancy. Furthermore, the paper found evidence of nonlinearities in the dose–response pattern. This present paper builds on that previous work to derive a BMD that can be used to guide clinical guidelines.

Here, we summarize some of the key features of the six cohorts and refer readers to two earlier articles from our group (Jacobson et al., 2021, 2024) for further remarks on important epidemiological issues and limitations of these studies. A critically important issue is the reliability of the methods used to assess maternal alcohol usage during pregnancy. All interviews took place during pregnancy, with the exception of one study in which mothers were interviewed shortly after delivery. The number of maternal interviews varied by study. In these interviews, detailed information was obtained regarding their alcohol consumption, as well as various other baseline characteristics, including family socioeconomic status, maternal education (years), and cognitive ability. Because the covariates were assessed and reported slightly differently among the cohorts, we computed propensity scores for each study by fitting separate models for alcohol exposure (Akkaya Hocagil et al., 2024). In computing these propensity scores, we followed the recommended practice (Austin, 2011) of including all relevant, available covariates. To avoid potential overadjustment bias (Schisterman et al., 2009), only covariates evaluated at the time of study enrollment, during pregnancy, or around the time of the child's birth were included (See Table S1 in Supplementary Materials). Propensity scores were then used as a covariate in a global response model characterizing the effects of alcohol on cognitive function. To reflect the fact that propensity scores were based on slightly different sets of covariates for each study, we also allowed the effect of the propensity score on the outcome to vary by the study by including a study by propensity score interaction (Akkaya Hocagil et al., 2024). Additional modeling details are provided in Section 10, as well as in Jacobson et al. (2024) who discuss the goodness-of-fit considerations and our rationale for using nonlinear models via GAMs.

All six cohorts reported continuous measures of maternal alcohol consumption during pregnancy. In particular, two exposure variables were reported: (1) alcohol dose per occasion, measured as the average number of standard drinks consumed per day on which drinking occurred, and (2) frequency of drinking, measured as the proportion of days during pregnancy when alcohol was consumed. Average absolute alcohol measured in ounces (oz) and consumed per day during pregnancy is then computed from the information collected directly in the interim for the other two exposure variables. For the purpose of our discussion, we consider all three of these exposure variables as measures of “dose” and retain the terminology benchmark dose.

In this paper, we present a new procedure for BMD estimation that leads to a benchmark dose contour quantifying levels and patterns of exposure that are associated with the development of clinically important cognitive deficits in children born to mothers who consume alcohol during pregnancy.

3 BMD ANALYSIS, SCALAR EXPOSURE VARIABLE

We begin with a review of a traditional BMD analysis involving a continuous response, a continuous exposure variable of interest, and several confounding variables. We formulate a linear model for the response as a function of the exposure while adjusting for confounders through the inclusion of a propensity score as an additional term in the regression. We then consider BMD analysis based on a GAM (Wood, 2017).

3.1 BMD estimation based on a linear model

3.1.1 Point estimation of the BMD

Let $$Y$$ denote the response of interest, $$X$$ a scalar continuous exposure variable, and $$\mbox{\boldmath $Z$}$$ a vector of potential confounding variables. To mitigate the effect of the potential confounders, we estimate propensity scores for each individual by modeling the exposure $$X$$ as a function of the potential confounders $$\mbox{\boldmath $Z$}$$ through a linear regression model:

$$$$\begin{equation} X = \alpha _{0} + \mbox{\boldmath $\alpha $}_{1} \mbox{\boldmath $Z$}+ R_x, \end{equation}$$$$(1)

with $$\alpha _{0}$$ and $$\mbox{\boldmath $\alpha $}_{1}$$ parameters to be estimated and $$R_x$$ a mean zero error term, independent of $$Z$$ and with constant variance. While propensity score methodology was originally developed for the context of a binary exposure (e.g., treated vs. nontreated), a number of relatively recent developments have extended the methods to the context of continuous exposures (Hirano & Imbens, 2004; Kosuke & van Dyk, 2004). Upon fitting (1), let $$S = \widehat{E}(X|\mbox{\boldmath $Z$}) = \hat{\alpha }_{0} + \hat{\mbox{\boldmath $\alpha $}}_{1} \mbox{\boldmath $Z$}$$ denote the propensity score, corresponding to the predicted value of the exposure $$X$$, given the measured covariates $$\mbox{\boldmath $Z$}$$. While there are a number of ways propensity scores can be used, we regress upon the propensity score in the linear outcome model.

$$$$\begin{equation} Y= \beta _0 + \beta _1 X + \beta _2 S + R_y, \end{equation}$$$$(2)

where $$R_y$$ is an error term, independent of both $$X$$ and $$S$$ and variance $$ {\rm var}(Y|X, S) = {\rm var}(R_y)=\sigma ^2$$. Provided there are no unmeasured confounders and that the propensity score model is correct, fitting model (2) yields a consistent estimate for the causal effect of $$X$$, represented by $$\beta _1$$. See Kosuke and van Dyk (2004) for further discussion about this approach. In our motivating study, $$Y$$ is a measure of cognitive function. High values of $$Y$$ are therefore desirable and hence if X represents some measure of alcohol intake we anticipate $$\beta _1 < 0$$. For convenience in what follows, we let $$\mu (X, S) = \beta _0 + \beta _1 X + \beta _2 S$$ denote the predicted outcome given exposure $$X$$ and propensity score $$S$$.

Crump (1984) originally proposed the estimation of a BMD in the context of a binary outcome but generalized the method to accommodate continuous outcomes. Specifically, he assumed that an individual has had an adverse effect if their outcome falls in the extreme tail (in our application this will be the lower tail) of the response distribution (Crump, 1995). In practice, one can specify the tail area directly (e.g., IQ scores below 75 are considered adverse) or indirectly (e.g., the lower 5% of the response distribution of unexposed individuals are considered to have an adverse response). There are computational and analytical advantages of the indirect approach (Budtz-Jørgensen et al., 2001), especially in settings where the definition of an abnormal response depends on the confounders and we will use this approach. We let $$p_\circ$$ denote the proportion of unexposed individuals considered to be in the adverse range. The threshold $$\tau _\circ (S)$$ is the $$100 p_{\circ }$$th percentile of the response distribution for a subject with propensity score $$S$$; thus the probability a response falls below $$\tau _\circ (S)$$ for unexposed individuals ($$X=0$$) with propensity score $$S$$ is $$p_\circ$$. In Figure 1, we display quantities $$p_\circ$$ and $$\tau _\circ (S)$$ on the left side with the fitted dose–response line depicting a decreasing mean response with increasing levels of exposure, consistent with what would be expected when $$Y$$ is cognitive function and $$X$$ is PAE.

Under the response model (2) and under the additional assumption that the error terms $$(R_{y})$$ are normally distributed,

$$$$\begin{equation*} P(Y \!<\! \tau _\circ (S) | X = 0, S) = P{\left(\frac{Y - \mu (0, S)}{\sigma } < \frac{\tau _\circ (S) - \mu (0, S)}{\sigma } \right)}, \end{equation*}$$$$

which if $$\Phi (\cdot)$$ is the cumulative distribution function for a standard normal random variable, gives the following relation between $$p_{\circ }$$ and $$\tau _{\circ }(S)$$:

$$$$\begin{equation*} p_\circ = \Phi ((\tau _\circ (S) - \mu (0, S))/\sigma). \end{equation*}$$$$

Solving for $$\tau _\circ (S)$$ gives

$$$$\begin{equation} \tau _\circ (S) = \mu (0, S) + \sigma \Phi ^{-1}(p_\circ) \,. \end{equation}$$$$(3)

Computing a BMD requires specification of the BMR, or the specified increase in the proportion of individuals with responses in the abnormal range (i.e., with $$Y< \tau _{\circ }(S)$$). Choosing a suitable value for the BMR generally requires judgment and consideration regarding the nature of the data set and the context in which the BMD will be used. A BMR will usually be chosen so that the corresponding BMD estimate can be reliably estimated, but then additional uncertainty factors might be applied to linearly extrapolate to the risk levels of specific regulatory interest (U.S. Environmental Protection Agency, 2012). When analysis is based on epidemiological data, a value of $$1\%$$ or $$5\%$$ is typically used for the BMR. The BMD, denoted here by $$X_B$$, is the dose of exposure that results in an increase in the proportion of responses in the abnormal range in the population corresponding to the BMR. That is, the BMD satisfies:

$$$$\begin{equation} P(Y < \tau _\circ (S) | X = X_B, S) = p_\circ + \rm {BMR}, \end{equation}$$$$(4)

where the right-hand side of (4) is the proportion of the unexposed population ($$X=0$$) with abnormally low responses plus the additional proportion in this range (BMR) due to the BMD of exposure. The normal distribution displayed on the right-hand side of Figure 1 is the response distribution at $$X_B$$, the BMD. The increase in risk of falling below the threshold $$\tau _{\circ }(S)$$ is evident from the larger shaded area. The decreasing mean response with increasing exposure $$X$$ is depicted on the right is determined by the linear dashed plot in Figure 3.

Again based on the linear model, if we further assume Gaussian errors we have

$$$$\begin{equation*} P(Y < \tau _\circ (S) | X = X_B, S) = \Phi {\left(\frac{\tau _\circ (S) - \mu (X_B, S)}{\sigma } \right)} \,, \end{equation*}$$$$

giving the relation

$$$$\begin{equation} \Phi ^{-1}(p_\circ + {\rm BMR}) = (\tau _\circ (S) - \mu (X_B, S))/\sigma \,. \end{equation}$$$$(5)

Replacing $$\tau _{\circ }(S)$$ in (5) with the expression in (3), we obtain

$$$$\begin{equation} \Phi ^{-1}(p_\circ + {\rm BMR}) - \Phi ^{-1}(p_\circ) = \frac{\mu (0, S) - \mu (X_B, S)}{\sigma }. \; \end{equation}$$$$(6)

For the additive linear model in (2), the right-hand side reduces to $$-\beta _1X_B/\sigma$$ which does not depend on $$S$$, and hence the confounding variables. It is this simplification which gives the appeal of the indirect approach to computing $$\tau _{\circ }(S)$$ (Budtz-Jørgensen et al., 2001). Letting

$$$$\begin{equation} A = A(p_\circ, \rm BMR) = \Phi ^{-1}(p_\circ + {\rm BMR}) - \Phi ^{-1}(p_\circ), \end{equation}$$$$(7)

we can explicitly write

$$$$\begin{equation} X_B = - A {\sigma }/{\beta _1}. \end{equation}$$$$(8)

The quantity in (8) is estimated by replacing the unknown parameters, $$\beta _1$$ and $$\sigma$$, by their respective estimated values to obtain $$\widehat{X}_B = - A \widehat{\sigma }/\widehat{\beta }_1$$.

3.1.2 Estimation of the BMDL in the linear model

The BMDL is a lower confidence limit for the BMD and reflects the statistical uncertainty associated with the use of estimated parameters instead of their true values. In practice, a lower one-sided $$95\%$$ or $$97.5\%$$ confidence interval (CI) is often used. This can be obtained in a number of ways. Crump (1984, 1995) proposed using likelihood-based limits to calculate the BMDL. Budtz-Jørgensen et al. (2001) derive the exact sampling distribution of the estimated BMD and from this are able to find an exact formula for a BMDL in the context of a linear model. For more complex dose–response functions, there will generally be no simple closed-form expression for the BMD. We therefore describe and implement an alternative approach that exploits the duality between hypothesis testing and CI estimation. CIs for a parameter of interest can be constructed by the set of values that are not rejected by a corresponding hypothesis test. All parameter values not rejected by profile likelihood ratio tests at the 5% significance level, for example, form 95% likelihood ratio-based CIs (Sprott, 2008); profile likelihood ratio intervals can be likewise constructed in multiparameter problems. An exact CI for a binomial probability is another example, and in failure time analysis, a 95% CI for the median survival time is defined implicitly by the survival times at which a test that the survival probability is 0.5 is not rejected. More generally, Hu and Kalbfleisch (2000) discuss the use of estimating functions to test hypotheses, and use resampling methods to carry out such tests and define corresponding CIs.

The BMD for the linear model is the solution to $$U(X, \mbox{\boldmath $\theta $})=0$$ where

$$$$\begin{equation} U(X; \mbox{\boldmath $\theta $}) = \displaystyle \frac{(- X) \, \beta _1}{\sigma } - A \end{equation}$$$$(9)

is a functional, with $$\mbox{\boldmath $\theta $}$$ a vector containing the unknown parameters $$\beta _1$$ and $$\sigma$$. The BMD is the value of $$X$$ such that $$X_B \beta _1 = -\sigma A$$ and the estimated BMD, $$\widehat{X}_B$$, is the solution to $$U(X; \widehat{\mbox{\boldmath $\theta $}}) = 0$$, where $$\theta$$ has been replaced by $$\hat{\theta }$$. It is the use of these estimated parameters in $$U(X; \mbox{\boldmath $\theta $})$$ that leads to uncertainty regarding the BMD and, hence, the need to compute a lower confidence limit defining the BMDL. The BMD defined in terms of (9) is a function of the maximum likelihood estimates and so is invariant to reparameterization.

Note that the null hypothesis $$H_0: X_B = X^\circ _B$$ is equivalent to the null hypothesis $$H_0: U(X^\circ _B;\mbox{\boldmath $\theta $})=0$$. It follows that a 95% CI for $$X_B$$ can be defined as the set of $$X^\circ _B$$ values such that this latter null hypothesis is not rejected based on the test statistic $$U(X^\circ _B; \widehat{\mbox{\boldmath $\theta $}})$$. There are several ways to do this. If $$U(X^\circ _B; \widehat{\mbox{\boldmath $\theta $}})$$ is assumed to be asymptotically normal and a variance estimate of $$U(X^\circ _B; \widehat{\mbox{\boldmath $\theta $}})$$ is known, then the null hypothesis $$H_0: U(X^\circ _B;\mbox{\boldmath $\theta $})=0$$ would not be rejected if

$$$$\begin{equation} -1.96 \sqrt {{\rm var} \lbrace U(X^\circ _B; \widehat{\mbox{\boldmath $\theta $}}) \rbrace }\le U(X^\circ _B; \widehat{\mbox{\boldmath $\theta $}}) \le 1.96 \sqrt {{\rm var} \lbrace U(X^\circ _B; \widehat{\mbox{\boldmath $\theta $}}) \rbrace }. \end{equation}$$$$(10)

That is a two-sided 95% CI for the BMD corresponds to the set of

$$X^\circ _B$$ values such that (10) is satisfied and the lower BMDL (lower 2.5% confidence limit) corresponds to the smallest value within this set. It is straightforward to establish that (10) holds if the following interval includes 0:

$$$$\begin{eqnarray} && \left(U(X^\circ _B; \widehat{\mbox{\boldmath $\theta $}}) -1.96 \sqrt {{\rm var} \lbrace U(X^\circ _B; \widehat{\mbox{\boldmath $\theta $}}) \rbrace } \,, \,U(X^\circ _B; \widehat{\mbox{\boldmath $\theta $}})\right.\nonumber\\ &&\quad +\;\left.1.96 \sqrt {{\rm var} \lbrace U(X^\circ _B; \widehat{\mbox{\boldmath $\theta $}}) \rbrace }\right). \end{eqnarray}$$$$(11)

Because the functional $$U(X; \mbox{\boldmath $\theta $})$$ equals $$-A$$ by definition when $$X$$ equals zero, the right-hand side of interval (11) will first cross 0. Hence, the BMDL will correspond to the smallest value of $$X_B^\circ$$ such that $$U(X^\circ _B; \widehat{\mbox{\boldmath $\theta $}}) +1.96 \sqrt {{\rm var} \lbrace U(X^\circ _B; \widehat{\mbox{\boldmath $\theta $}}) \rbrace }=0$$.

In complex semiparametric models, there are considerable computational advantages to deriving CIs for the BMD indirectly through the functional $$U(X^\circ _B; \widehat{\mbox{\boldmath $\theta $}})$$ since this quantity may generally be much easier to calculate than the BMD itself. As we show shortly, however, it may be difficult to derive an analytic expression for $${\rm var}\lbrace U(X^\circ _B; \widehat{\mbox{\boldmath $\theta $}})\rbrace$$. We circumvent this problem through use of a parametric bootstrap procedure (Efron, 2000). For the linear model, the first step is to fit (2) and compute the predicted response $$\hat{\mu }(X_{i}, S_{i})=\hat{\beta }_{0}+\hat{\beta }_{1}X_{i}+\hat{\beta }_{2}S_{i}$$ for each individual $$i$$, where $$X_i$$ and $$S_i$$ are the observed values of exposure and propensity score for that individual. Next, using a parametric bootstrap, we generate $$M$$ samples of size $$n$$ where for the $$m$$th sample we:

• 1. For each individual $$i$$, we generate a new residual $$R_{y_i}^{(m)}$$ from a Gaussian distribution with mean zero and variance $$\hat{\sigma }^2$$ where $$\hat{\sigma }^2$$ is the estimated error variance based on the model fitted to the original data;
• 2. Form a new outcome for each individual by adding the predicted value and residual from Step 1 to get $$Y^{(m)}_{i}= \hat{\mu }(X_{i}, S_{i})+R^{(m)}_{y_{i}}$$;
• 3. Given the $$m$$th bootstrap sample we fit the model (2) and record $$\widehat{\mbox{\boldmath $\theta $}}^{(m)}=(\hat{\beta }^{(m)}, \hat{\sigma }^{(m)})$$ and compute the value of the functional $$U(X^\circ _B;\widehat{\mbox{\boldmath $\theta $}}^{(m)})$$ at a candidate BMD value $$X^\circ _B$$ on a grid of points.

Once the $$M$$ boostrap samples have been generated a bootstrap variance estimate $$\widehat{\rm var} \lbrace U(X^\circ _B; \widehat{\mbox{\boldmath $\theta $}})\rbrace$$ can be obtained and (10) or (11) can be used.

An alternative approach to interval estimation is to base it on the quantiles of the bootstrap distribution of the functional—this has the advantage of relaxing the normality assumption required for CI construction based on $$\widehat{\rm var} \lbrace U(X^\circ _B; \widehat{\mbox{\boldmath $\theta $}})\rbrace$$ and therefore handles a wider range of problems. The procedure is as follows. Let $$Q_{p}(X_B^\circ; \widehat{\mbox{\boldmath $\theta $}})$$ be the $$100 p$$th percentile of the bootstrap distribution of $$U(X_B^\circ; \widehat{\mbox{\boldmath $\theta $}})$$ and consider an approximate 95% CI for the functional $$U(X^\circ _B; \mbox{\boldmath $\theta $})$$ given by

$$$$\begin{equation*} \lbrace Q_{2.5}(X^\circ _B) < U(X^\circ _B; \hat{\mbox{\boldmath $\theta $}}) < Q_{97.5}(X^\circ _B) \rbrace \,. \end{equation*}$$$$

Because the quantity $$U(X; \mbox{\boldmath $\theta $})$$ equals $$-A$$ when $$X=0$$, the BMDL corresponds to the smallest value of $$X^{\circ }_{B}$$ for which $$Q_{97.5}(X^\circ _{B})=0$$. This value can be obtained by a simple grid search. A large number of bootstrap samples is recommended when estimating extreme quantiles to ensure accurate coverage probabilities.

In Section 10, we will see that in the linear model setting, the parametric bootstrap approach yields almost identical BMDL estimates as the analytic approach of Budtz-Jørgensen et al. (2001), giving confidence that the bootstrap approach will also work well in settings where the simpler solution is not available. In particular, the next section describes how the methodology described above can be used to estimate a BMD and BMDL based on a flexible nonlinear model where no simple analytic solution is available. Section  7 will report on a small simulation study designed to assess the performance of the approach.

4 BMD VIA GAMS WITH A BIVARIATE EXPOSURE

5 APPLICATION TO THE PAE STUDY

Here, we apply the proposed methodology in Section 2 to examine the effect of PAE on child cognition. The data are from six US longitudinal cohort studies aiming to quantify levels and patterns of exposure that are associated with the development of clinically important cognitive deficits in children born to mothers who consume alcohol during pregnancy. In all six cohorts, data on maternal alcohol consumption were quantified in terms of daily alcohol intake (standard drinks or oz of absolute alcohol per day averaged across pregnancy as well as the frequency of drinking days and the amount of alcohol consumed on each drinking day during pregnancy. The distributions of daily alcohol intake are summarized by cohort in the Supplementary Materials. The daily alcohol intake (standard drinks or oz of absolute alcohol per day averaged across pregnancy and the amount of alcohol consumed on each drinking day during pregnancy were highly skewed (skewness $$>$$ 3.0). We adopted a log($$X + 1$$) transformation to improve model stability. Further discussion related to our decision to use this logarithmic transformation may be found in Supplementary Materials. We also provide there some more detailed discussion regarding the distributions of daily alcohol intake averaged across pregnancy by cohort, showing histograms of these variables in the original scale and in the log scale (see Figures S1 and S2). In our recent article, Jacobson et al. (2024) discuss the results of exploratory analyses on the harmonized data set and highlight the need for models accommodating nonlinear interaction effects when considering bivariate exposure.

We consider a BMD analysis for a continuous response, cognitive function score, based on a linear model involving a single continuous exposure, daily alcohol intake averaged across pregnancy (log scale), and the estimated propensity score. We set $$BMR = 0.01$$ and $$p_{0} = 0.025$$ and obtain a point estimate of the BMD and corresponding $$97.5\%$$ lower confidence limit BMDL. We estimated the BMDL using the two approaches described in Section 4. In the first approach, we approximated the $$97.5\%$$ lower confidence limit for the BMD using the analytical formula derived by Budtz-Jørgensen et al. (2001). For the parametric bootstrap approach, a lower confidence limit for the BMD is approximated by the 2.5th percentile from the bootstrapped distribution of the estimated BMD, based on simulating 10,000 new parametric bootstrap samples from the fitted model. Since we use a parametric bootstrap approach, cohort sample sizes for all six cohort studies and exposure distributions are unchanged across the bootstrapped data sets. The first row of Table 1 shows the estimated BMD and two BMDL estimates which are seen to be very similar.

To relax the assumption of a linear effect of the average alcohol consumption per day across pregnancy and the propensity score, we then fitted a GAM using the gam function from the mgcv package in R (Wood, 2017) and a BMDL was estimated using the parametric bootstrap approach described in Section 6, setting $$M=10,000$$. We considered a grid search consisting of 100 candidate values for $$X_B^\circ$$, evenly distributed between the minimum (0.00) and maximum (3.54) amount of alcohol consumed on each drinking day during pregnancy recorded in the data set and generated a bootstrap distribution of the functional $$U(X^{\circ };\widehat{\mbox{\boldmath $\theta $}})$$ at each grid point. As discussed previously, the BMDL was chosen as smallest exposure level where the 97.5th percentile of the bootstrap distribution equaled 0. The resulting BMD and BMDL estimates appear in the second row of Table 1. The BMD estimate from the GAM is higher than what those of the linear model but the corresponding BMDL estimates are similar.

To provide a visualization of how the bootstrap approach works, Figure 3 shows smoothed histograms of the bootstrap distribution of $$U(X^{\circ }; \widehat{\mbox{\boldmath $\theta $}})$$ at a discrete set of values of $$X^{\circ }$$. The vertical line represents the point 0 and for each histogram, tail areas have been shaded below the 2.5th percentile, and above the 97.5th percentile. The BMDL corresponds to the lowest value of $$X_B^\circ$$ for which 97.5th percentile of the bootstrap distribution equals 0. In the figure, we can see that this indeed happens right around $$X^{\circ }_B=0.32$$ and this is consistent with our analysis which yields the value of 0.31 for the BMDL.

We next consider estimation of a BMD curve and associated BMDL in terms of a bivariate exposure. To demonstrate our proposed method, we investigate the role of the pattern of exposure where we distinguish between consistent low levels of exposure and episodic periods of higher exposure (binge drinking) that may result in a similar total alcohol intake over a given period of time. We fitted the model described in Section 9 to predict cognition, with $$X_1$$ and $$X_2$$ representing alcohol dose per occasion (log scale) and frequency of drinking, respectively. Using the approach described in Section 9, we used a bivariate GAM in order to construct a flexible three-dimensional surface predicting cognitive function in relation to alcohol dose per occasion (log scale) and drinking frequency, while adjusting for confounders via propensity scores. The bivariate GAM approach was justified on the basis of previous exploratory data analysis that suggested the presence of both nonlinearities and interactions between these two different dimensions of alcohol exposure (Jacobson et al., 2024). Further discussion may be found in Supplementary Materials.

Figure 4 presents the estimated dose–response surface obtained from the fitted GAM. The curve depicted at the front edge of the surface plot shows a nonlinear effect of dose per occasion among infrequent drinkers. At low levels of dose per occasion, the curve is quite flat, suggesting little adverse effect at these drinking levels, for infrequent drinkers. As the dose per occasion increases beyond about 1.7 oz, the curve decreases as the frequency increases. At the highest levels of drinking frequency (i.e., daily or 30 days per month), the curve slopes sharply downward in a relatively linear manner, suggesting that for frequent drinkers, any increase in dose/occasion is associated with an appreciable adverse effect on cognitive function.

Figure 5 shows the estimated dose–response surface of the exposure effects with CI surfaces at $$\pm$$2 standard errors. These wider pointwise confidence surfaces reflect greater uncertainty in the model at higher levels of drinking, where there were fewer women in the cohorts.

The solid line in Figure 6 shows the estimated BMD curve overlayed on a scatter plot of alcohol per occasion (y-axis) versus drinking frequency in the proportion of days per month (x-axis). This BMD curve shows the combinations of dose per occasion and drinking frequency associated with the effects defining the BMD. We obtained a point estimate for the BMD curve for the combination of values of alcohol per occasion and drinking frequency by a grid search in the $$X$$-space. We estimated the BMDL contour using the parametric bootstrap approach described in Section 9. We searched over a grid containing 100 different candidate values of $$X_{1B}$$ and $$X_{2B}$$, yielding a total of 10,000 grid points to search. We used $$M=10,000$$ bootstrap samples for the estimated functional $$U(\mbox{\boldmath $X$}_B^{\circ };\widehat{\theta })$$ at each candidate point in the grid, $$\mbox{\boldmath $X$}_B^{\circ }=(X_{1B},X_{2B})$$. Next, for each bootstrapped sample, we computed the set of indicators, $$\Delta ^m(\mbox{\boldmath $X$}^{\circ }_{B})=I(U(\mbox{\boldmath $X$}^{\circ }_{B},\widehat{\mbox{\boldmath $\theta $}}^{(m)})>0)$$ of whether the $$m$$th bootstrapped estimated estimating function at each grid point exceeded zero. Finally, we fitted a generalized additive logistic regression model with a tensor product of the logarithm of absolute alcohol per occasion and drinking frequency to find the combinations of the logarithm of absolute alcohol per occasion and drinking frequency where $$97.5\%$$ of the bootstrap distribution exceeds the threshold. Figure 6 shows the estimated BMDL (dashed line) that is obtained from the parametric bootstrap approach. The results suggest an estimated BMDLs range from around 1.5 oz for infrequent drinkers to just over 0.6 oz for frequent or daily drinkers.

All analyses were conducted in R (R Core Team, 2021). The R Code for our bivariate modeling can be found in the Appendix and the data used in the analysis have been uploaded to the Harvard Dataverse (https://doi.org/10.7910/DVN/SCLH9W).

6 DISCUSSION AND CONCLUSIONS

We have developed a flexible joint dose–response model for estimating a BMD contour for bivariate exposures in quantitative risk assessment. Specifically, we estimated joint BMD points via GAMs, while adjusting for confounding variables through regression adjustment based on propensity scores, and derived the corresponding lower confidence limit using an innovative parametric bootstrap approach based on generating a bootstrap distribution of a functional, then inverting this to obtain the lower limit of interest, which corresponds to the BMDL. Our methodology has provided important insights into the effects of PAE on child cognition and in particular regarding the interacting effects of dose per occasion and frequency of drinking. It also can be used in settings directed at more nuanced questions regarding the timing and intensity of exposure (Kriebel et al., 2007) when such detailed information is available.

An important caveat is that our semiparametric analysis based on generalized additive modeling or GAMs may encounter problems for data sets that have only a small number of control individuals, or those with zero exposure. This is because the method requires having a reliable estimate of the predicted cognitive score for unexposed individuals (see Equation 16). In our application, at least from a technical perspective (see discussion below), this is not an issue since our data include a large number of nondrinkers or abstainers. In practice, however, it will be a problem for settings involving ubiquitous exposures such as methylmercury (Whitney & Ryan, 2013) where the majority of subjects will have at least low-level exposure levels. In such settings, it might be useful to consider modified approaches such as isotonic regression or monotone splines that impose some mild conditions on the shape of the dose–response for low levels of exposure (Piegorsch et al., 2012).

We demonstrated our approach using data from six large longitudinal cohort studies that aimed to assess the effects of PAE on child cognitive function. In order to assess the performance of the bootstrap methodology, we started with an analysis based on a single exposure variable, the logarithm of average daily alcohol consumption across pregnancy, estimating a BMD based on a linear model. This allowed us to compare BMDL estimates based on our bootstrap approach with the analytic solution described by Budtz-Jørgensen et al. (2001). We found the two BMDL estimates to be very close, suggesting that our bootstrap approach worked well. We then applied the extended methodology to the bivariate exposure setting where alcohol exposure was characterized by both dose per occasion and frequency of drinking. The estimated dose–response surface suggests a nonlinear interaction. Specifically, the dose–response curve appears to be quite flat at low levels of dose per occasion (less than 1.7 oz), for infrequent drinkers. But at higher levels of dose per occasion, the estimated curve starts sloping downward. By contrast, the dose–response curve at the highest levels of drinking frequency (i.e., daily or 30 days per month), the curve slopes sharply downward in a relatively linear manner, suggesting that for frequent drinkers, any increase in dose/occasion is associated with an increase in the adverse effect on cognitive function. This visual interpretation is supported by the BMD analysis. For infrequent drinkers, the estimated BMD is quite high (around 5.7 oz per occasion), falling to lower than 2.0 oz for frequent or daily drinkers. The BMDLs range from around 1.5 oz for infrequent drinkers to just over 0.6 oz for frequent or daily drinkers. These findings suggest that whether there are clinically significant adverse effects on a child's cognitive development depends in many cases on the pattern of the mother's drinking during pregnancy such that there appears to be little effect if the mother drinks $$<$$3 drinks/occasion infrequently (i.e., $$<$$2 days/month). In this context, it is important to emphasize that there are several critical caveats that need to be kept in mind. First of all, there are always uncertainties associated with epidemiological analysis. Also, a BMD analysis is not intended to identify levels where there is no effect, but rather where risk levels are not excessively raised relative to background. For these reasons, the US Office of the Surgeon General's (2005) recommendation that all women abstain from alcohol during pregnancy continues to be the most appropriate guideline. It also important to emphasize that the present paper has not set out to present a definitive risk analysis for the PAE, but rather to present a methodology for BMD and BMDL computation in the context of a bivariate exposure. A more detailed applied analysis focusing on the interpretation and clinical implications of our analysis is underway. This follow-up paper will provide an opportunity to address in more depth some of the nuances and challenges associated with epidemiological explorations of the health effects of alcohol. For example, some argue (Mugavin et al., 2020) that abstainers can differ from low-level drinkers for a variety of reasons that are not always well captured by standard covariate adjustment.

In our motivating example, we addressed the challenge of modeling an outcome of interest as a function of exposure that can be measured in two dimensions. The methodology can also potentially apply in settings involving mixtures of different exposures. The methodology extends naturally beyond the two-dimensional setting to higher-order exposures, though there may be practical limitations in such cases. For example, tensor product terms with dimensions higher than two may be unstable due to sparse data. There are opportunities for further methodological developments in this area.